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)