This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English version of the page.

Note: This page has been translated by MathWorks. Click here to see
To view all translated materials including this page, select Country from the country navigator on the bottom of this 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

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;
model = createpde(numberOfPDE);

Boundary Conditions

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

axis equal
title 'Geometry With Edge Labels Displayed'

Create a geometry entity.

pg = geometryFromEdges(model,g);

Solution is zero at all outer edges.


Generate Initial Mesh

[p,e,t] = initmesh(g);
[p,e,t] = refinemesh(g,p,e,t);
[p,e,t] = refinemesh(g,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(model,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(model,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(model,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


Compare with Solution Found without Domain Decomposition

[K,F] = assempde(model,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?