www.gusucode.com > pde 案例源码 matlab代码程序 > pde/pdedemo4.m

    %% 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(TM) 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(TM) 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.

%       Copyright 1994-2014 The MathWorks, Inc.

%% 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'));


displayEndOfDemoMessage(mfilename)