Documentation

This is machine translation

Translated by Microsoft
Mouse over text to see original. Click the button below to return to the English verison of the page.

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.

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

g = @lshapeg;
c = 1;
a = 0;
f = 1;

% Create a pde entity for a PDE with a single dependent variable
numberOfPDE = 1;
pdem = createpde(numberOfPDE);

Boundary Conditions

Plot the geometry and display the edge labels for use in the boundary condition definition.

figure;
pdegplot(g,'EdgeLabels','on');
axis equal
title 'Geometry With Edge Labels Displayed'

% Create a geometry entity
pg = geometryFromEdges(pdem,g);
% Solution is zero at all outer edges
applyBoundaryCondition(pdem,'dirichlet','Edge',(1:10),'u',0);

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(pdem,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(pdem,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(pdem,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(pdem,p,e,t,1,0,1);
u1 = K\F;
fprintf('Difference between solution vectors = %g\n', norm(u-u1,'inf'));
Difference between solution vectors = 0.000231536

Was this topic helpful?