Quantcast

Partial Differential Equation Toolbox

Poisson's Equation Using Domain Decomposition

This example shows how to numerically solve a Poisson's equation using the assempde function in the Partial Differential Equation Toolbox™ in conjunction with domain decomposition.

The Poisson's equation we are solving is

$-\Delta u = 1$

on the L-shaped membrane with zero-Dirichlet boundary conditions.

The Partial Differential Equation Toolbox™ is designed to deal with one-level domain decomposition. If the domain has a complicated geometry, it is often useful to decompose it into the union of two or more subdomains of simpler structure. In this example, an L-shaped domain is decomposed into three subdomains. The FEM solution is found on each subdomain by using the Schur complement method.

Problem Definition

The following variables will define our problem:

  • g: A specification function that is used by initmesh and refinemesh. For more information, please see the documentation page for lshapeg and pdegeom.

  • b: A boundary file used by assempde. For more information, please see the documentation pages for lshapeb and pdebound.

  • c, a, f: The coefficients and inhomogeneous term.

g='lshapeg';
b='lshapeb';
c=1;
a=0;
f=1;

Generate Initial Mesh

[p,e,t]=initmesh(g);
[p,e,t]=refinemesh(g,p,e,t);
[p,e,t]=refinemesh(g,p,e,t);
figure;
pdemesh(p,e,t); axis equal

Find Common Points

np=size(p,2);
cp=pdesdp(p,e,t);

Allocate Space

Matrix C will hold a Schur complement.

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

Assemble First Domain and Update Complement

[i1,c1]=pdesdp(p,e,t,1);ic1=pdesubix(cp,c1);
[K,F]=assempde(b,p,e,t,c,a,f,[],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;

Assemble Second Domain and Update Complement

[i2,c2]=pdesdp(p,e,t,2);ic2=pdesubix(cp,c2);
[K,F]=assempde(b,p,e,t,c,a,f,[],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;

Assemble Third Domain and Update Complement

[i3,c3]=pdesdp(p,e,t,3);ic3=pdesubix(cp,c3);
[K,F]=assempde(b,p,e,t,c,a,f,[],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 For Solution on Each Subdomain.

u=zeros(np,1);
u(cp)=C\FC; % Common points
u(i1)=K1\(e1-a1'*u(c1)); % Points in SD 1
u(i2)=K2\(e2-a2'*u(c2)); % Points in SD 2
u(i3)=K3\(e3-a3'*u(c3)); % Points in SD 3

Plot FEM Solution

figure
pdesurf(p,t,u)

Compare with Solution Found without Domain Decomposition

[K,F]=assempde(b,p,e,t,1,0,1);u1=K\F;
norm(u-u1,'inf')
ans =

   1.7397e-04