MATLAB Examples

Domain Decomposition Problem

This example shows how to perform one-level domain decomposition for complicated geometries, where you can decompose this geometry into the union of more subdomains of simpler structure. Such structures are often introduced by the PDE Modeler app.

Assume now that $\Omega$ is the disjoint union of some subdomains $\Omega_1, \Omega_2, \ldots, \Omega_n$. Then you could renumber the nodes of a mesh on $\Omega$ such that the indices of the nodes of each subdomain are grouped together, while all the indices of nodes common to two or more subdomains come last. Since K has nonzero entries only at the lines and columns that are indices of neighboring nodes, the stiffness matrix is partitioned as follows:

$$ K =
\left(
\begin{array}{ccccc}
K_1 & 0 & \cdots & 0 & B_1^T \\
0 & K_2 & \cdots & 0 & B_2^T \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & \cdots & 0 & B_n^T \\
B_1 & B_2 & \cdots & B_n & C \\
\end{array}
\right) $$

while the right side is

$$ F =
\left(
\begin{array}{c}
f_1 \\
f_2 \\
\vdots \\
f_n \\
f_c \\
\end{array}
\right) $$

The Partial Differential Equation Toolbox™ function assempde can assemble the matrices $K_j$, $B_j$, and $C$ separately. You have full control over the storage and further processing of these matrices.

Furthermore, the structure of the linear system

$$Ku = F$$

is simplified by decomposing $K$ into the partitioned matrix.

Now consider the geometry of the L-shaped membrane. You can plot the geometry of the membrane by typing

pdegplot('lshapeg')

Notice the borders between the subdomains. There are three subdomains. Thus the matrix formulas with $n = 3$ can be used. Now generate a mesh for the geometry:

[p,e,t] = initmesh('lshapeg');
[p,e,t] = refinemesh('lshapeg',p,e,t);
[p,e,t] = refinemesh('lshapeg',p,e,t);

So for this case, with $n = 3$, you have

$$ \left(
\begin{array}{cccc}
K_1 & 0 & 0 & B_1^T \\
0 & K_2 & 0 & B_2^T \\
0 & 0 & K_3 & B_3^T \\
B_1 & B_2 & B_3 & C \\
\end{array}
\right)
\left(
\begin{array}{c}
u_1 \\
u_2 \\
u_3 \\
u_c \\
\end{array}
\right) =
\left(
\begin{array}{c}
f_1 \\
f_2 \\
f_3 \\
f_c \\
\end{array}
\right) $$

and the solution is given by block elimination:

$$ \left( C - B_1 K_1^{-1}B_1^T - B_2 K_2^{-1} B_2^T - B_3 K_3^{-1} B_3^T
\right) = f_c - B_1 K_1^{-1} f_1 - B_2 K_2^{-1} f_2 - B_3 K_3^{-1} f_3 $$

$$ u_1 = K_1^{-1} \left( f_1 - B_1^T u_c \right) $$

$$ \cdots $$

In the following MATLAB™ solution, a more efficient algorithm using Cholesky factorization is used:

time = [];
np = size(p,2);
% Find common points
c = pdesdp(p,e,t);

nc = length(c);
C = zeros(nc,nc);
FC = zeros(nc,1);

[i1,c1] = pdesdp(p,e,t,1);
ic1 = pdesubix(c,c1);
[K,F] = assempde('lshapeb',p,e,t,1,0,1,time,1);
K1 = K(i1,i1);
d = symamd(K1);
i1 = i1(d);
K1 = chol(K1(d,d));
B1 = K(c1,i1);
a1 = B1/K1;
C(ic1,ic1) = C(ic1,ic1)+K(c1,c1)-a1*a1';
f1 = F(i1);
e1 = K1'\f1;
FC(ic1) = FC(ic1)+F(c1)-a1*e1;

[i2,c2] = pdesdp(p,e,t,2);
ic2 = pdesubix(c,c2);
[K,F] = assempde('lshapeb',p,e,t,1,0,1,time,2);
K2 = K(i2,i2);
d = symamd(K2);
i2 = i2(d);
K2 = chol(K2(d,d));
B2 = K(c2,i2);
a2 = B2/K2;
C(ic2,ic2) = C(ic2,ic2)+K(c2,c2)-a2*a2';
f2 = F(i2);
e2 = K2'\f2;
FC(ic2) = FC(ic2)+F(c2)-a2*e2;

[i3,c3] = pdesdp(p,e,t,3);
ic3 = pdesubix(c,c3);
[K,F] = assempde('lshapeb',p,e,t,1,0,1,time,3);
K3 = K(i3,i3);
d = symamd(K3);
i3 = i3(d);
K3 = chol(K3(d,d));
B3 = K(c3,i3);
a3 = B3/K3;
C(ic3,ic3) = C(ic3,ic3)+K(c3,c3)-a3*a3';
f3 = F(i3);
e3 = K3'\f3;
FC(ic3) = FC(ic3)+F(c3)-a3*e3;

% Solve
u = zeros(np,1);
u(c) = C\ FC;
u(i1) = K1\(e1-a1'*u(c1));
u(i2) = K2\(e2-a2'*u(c2));
u(i3) = K3\(e3-a3'*u(c3));

The problem can also be solved by typing

% Compare with solution not using subdomains
[K,F] = assempde('lshapeb',p,e,t,1,0,1);
u1 = K\F;
norm(u-u1,'inf')
pdesurf(p,t,u)
ans =

   1.7397e-04

For the entire example, see docid:pde_examples.example-ex59440738.