poisson.m

C++编译指导

function u = Poisson(node,elem,Dirichlet,Neumann,f,u_D,g)
% POISSON solve the 2-D Poisson equation 
%     -\Delta u = f, 
% in a domain described by node and elem,
% with boundary edges Dirichlet, Neumann.
%     u = u_D on the Dirichelet boundary  
%    du/dn = g on the Neumann boundary 
%
% Input: 
%   node, elem: standard mesh data;
%   Dirichlet, Neumann: boundary edges;
%   f:right side or data
%   u_D: for Dirichelet condition
%   g: for Neumann condition
%
% Output:
%   u: solution on the current mesh
%
% L. Chen 03-26-2006
%
% The code is based on the following reference.
% NOTE: it is optimized using vectorizing lauange 
%       to avoid for loops
%
% Reference:
%
%    Jochen Alberty, Carsten Carstensen, Stefan Funken,
%    Remarks Around 50 Lines of MATLAB:
%    Short Finite Element Implementation,
%    Numerical Algorithms,
%    Volume 20, pages 117-137, 1999.

N = size(node,1); NT = size(elem,1);
A = sparse(N,N); u = sparse(N,1); b = sparse(N,1);

%  Assembing stiffness matrix
loce(:,1,:) = node(elem(:,3),:)-node(elem(:,2),:);
loce(:,2,:) = node(elem(:,1),:)-node(elem(:,3),:);
loce(:,3,:) = node(elem(:,2),:)-node(elem(:,1),:);
area = 0.5*abs(-loce(:,3,1).*loce(:,2,2)+loce(:,3,2).*loce(:,2,1));
for i = 1:3
    for j = 1:3
        Aij = (loce(:,i,1).*loce(:,j,1)+loce(:,i,2).*loce(:,j,2))./(4*area);
        A = A + sparse(elem(:,i),elem(:,j),Aij,N,N);
    end
end

% assembing right side
center = (node(elem(:,1),:)+node(elem(:,2),:)+node(elem(:,3),:))/3;
b = sparse(elem,ones(NT,3),feval(f,center)'.*area/3*[1,1,1],N,1);

% Neumann boundary conditions
if(~isempty(Neumann))
    Nedge = node(Neumann(:,1),:) - node(Neumann(:,2),:);
    d = sqrt(sum(Nedge.^2,2)); 
    mid = (node(Neumann(:,1),:) + node(Neumann(:,2),:))/2;
    b = b + sparse(Neumann,ones(size(Neumann,1),2),d.*g(mid)/2*[1,1],N,1); 
end

%  Dirichlet boundary conditions
%  Assign the corresponding entries of U,and adjust the right hand side. 
% u = sparse(N,1);
BoundNodes = unique(Dirichlet);
u(Dirichlet) = feval(u_D,node(Dirichlet,:));
b = b - A * u;

% Compute the solution by solving A * U = B for the remaining unknown values of U.
FreeNodes = setdiff(1:N,BoundNodes);
if size(FreeNodes)>0
	u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);
end

%  Graphic representation.
hold off
trisurf(elem,node(:,1),node(:,2),u','facecolor','interp');
view(67.5,30);
title('Solution to the Problem')
%view(2), axis equal