www.pudn.com > AFEM@matlab.rar > Poisson.m, change:2006-11-16,size:5147b

function [u, energy] = Poisson(mesh, f, g_D, g_N)
% POISSON solve the 2-D Poisson equation
%          -\Delta u = f,
% in the current mesh with boundary conditions
%           u = g_D  on the Dirichelet boundary
%       du/dn = g_N  on the Neumann boundary
%
% USAGE
%            [u] = Poisson(mesh, f, g_D, g_N)
%    [u, energy] = Poisson(mesh, f, g_D, g_N)
%
% INPUT
%    mesh:  current mesh
%    f:     right side or data
%    g_D:   Dirichelet condition
%    g_N:   Neumann condition
%
% OUTPUT
%    u:     solution on the current mesh
%    energy:energy of the discrete solution u
%
% 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.
%
% NOTE
%    It is optimized using vectorizing lauange to avoid for loops
%

% L. Chen & C. Zhang 11-15-2006

%--------------------------------------------------------------------------
% Initialize the data
%--------------------------------------------------------------------------
N = size(mesh.node,1);
A = sparse(N,N);
u = zeros(N,1);

%--------------------------------------------------------------------------
% Compute vedge: edge as a vector and area of each element
%--------------------------------------------------------------------------
ve(:,:,1) = mesh.node(mesh.elem(:,3),:)-mesh.node(mesh.elem(:,2),:);
ve(:,:,2) = mesh.node(mesh.elem(:,1),:)-mesh.node(mesh.elem(:,3),:);
ve(:,:,3) = mesh.node(mesh.elem(:,2),:)-mesh.node(mesh.elem(:,1),:);
area = 0.5*abs(-ve(:,1,3).*ve(:,2,2)+ve(:,2,3).*ve(:,1,2));

%--------------------------------------------------------------------------
% Assemble Stiffness matrix
%--------------------------------------------------------------------------
for i = 1:3
for j = 1:3
Aij = (ve(:,1,i).*ve(:,1,j)+ve(:,2,i).*ve(:,2,j))./(4*area);
A = A + sparse(mesh.elem(:,i),mesh.elem(:,j),Aij,N,N);
end
end
%   for j = 1 : NT
%       A(elem(j,:),elem(j,:)) = A(elem(j,:),elem(j,:)) ...
%    	                           + localstiffness(node(elem(j,:),:));
%   end
% with localstiffness is a function to compute lcoal stiffness matrix.
% To avoide loop for large NT, we use 'sparse' command here.

%--------------------------------------------------------------------------
% Assemble Mass matrix
%--------------------------------------------------------------------------
M = sparse(mesh.elem(:,[1,1,1,2,2,2,3,3,3]), ...
mesh.elem(:,[1,2,3,1,2,3,1,2,3]), ...
area*[2,1,1,1,2,1,1,1,2]/12, N, N);
%--------------------------------------------------------------------------
%   for j = 1 : NT
%       M(elem(j,:),elem(j,:)) = area(j)/12*[2,1,1;
%                                            1,2,1;
%                                            1,1,2];
%   end
% To avoide loop for large NT, we use 'sparse' command here.
% It is also very neat, by the way.
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Assemble right-hand-side
%--------------------------------------------------------------------------
b = M*feval(f,mesh.node);
% feval(f,mesh.node) is the nodal interpolation of f_I and M*f_I gives the
% right side. It is equivalent to the middle point rule

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

%--------------------------------------------------------------------------
% Dirichlet boundary conditions
%--------------------------------------------------------------------------
bdNode = unique(mesh.Dirichlet);
u(bdNode) = feval(g_D,mesh.node(bdNode,:));
b = b - A * u; % adjust the right hand side

%--------------------------------------------------------------------------
% Solve the linear system A * U = B for the free nodes
%--------------------------------------------------------------------------
freeNode = find(mesh.type>0);
freeNode = setdiff(freeNode, bdNode);
if size(freeNode)>0
u(freeNode) = A(freeNode, freeNode) \ b(freeNode);
end

%--------------------------------------------------------------------------
% Compute the energy of u
%--------------------------------------------------------------------------
if (nargout > 1)
energy = full(0.5*u'*(A*u)- u'*M*feval(f,mesh.node));
end % otherwise do not need to compute the discrete energy

%--------------------------------------------------------------------------
% end of function POISSON
%--------------------------------------------------------------------------