poisson.m

五点差分型多重网格方法：各种插值算子的比较）

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 = \int_T grad u_i * grad u_j
        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
% More readable code is
%   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);
%--------------------------------------------------------------------------
% More readable code is
%   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
%--------------------------------------------------------------------------