wmcdfm.m

C++编译指导

function mesh=WMCDFM(mesh, eps, a1, a2, b, f, g_D,g_N)
% CDFM solve the 2-D convection-dominated diffusion equation 
%         -\epslon\delta u +a*grad(u)+bu=0, 
% in the current mesh with boundary conditions
%           u = 0  on the Dirichelet boundary  
%    
%
% USAGE
%            [mesh] =CDFM(mesh, f, g_D, )
%
% INPUT 
%    \epslon
%    \delta
%    mesh:  current mesh
%    f:     right side or data
%    g_D:   Dirichelet condition
%    g_N:   Neumann condition
%
% OUTPUT
%    u:     solution on the current mesh
%    
%
% 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 10-08-2006

%--------------------------------------------------------------------------
% Initialize the data 
%--------------------------------------------------------------------------
N = size(mesh.node,1); NT = size(mesh.elem,1);
A = sparse(N,N); b1 = sparse(N,1); u = sparse(N,1); 
A1=sparse(NT,2);vb=sparse(NT,1);tlength=sparse(NT,1);area=sparse(NT,1);

%--------------------------------------------------------------------------
% Compute vedge: edge as a vector and area of each element
%--------------------------------------------------------------------------
center = ( mesh.node(mesh.elem(:,1),:) + mesh.node(mesh.elem(:,2),:)...
                                  + mesh.node(mesh.elem(:,3),:) )/3;
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));
tlength=(ve(:,1,1).^2 + ve(:,2,1).^2 + ve(:,1,2).^2 + ve(:,2,2).^2 ...
                                     +ve(:,1,3).^2 + ve(:,2,3).^2)/3;
A1(:,1)=feval(a1,center);A1(:,2)=feval(a2,center);
vb=(feval(b,center))';



%--------------------------------------------------------------------------
% Assemble Stiffness matrix
%--------------------------------------------------------------------------
Aij = zeros(NT,1); % 
for i = 1:3
    for j = 1:3
        Aij = eps*((ve(:,1,i).*ve(:,1,j)+ve(:,2,i).*ve(:,2,j))./(4*area)) ...
             +(-A1(:,1).*ve(:,2,i)+A1(:,2).*ve(:,1,i))/6  +  (area.*vb)/9; 
                  
        A = A + sparse(mesh.elem(:,j),mesh.elem(:,i),Aij,N,N);
    end
end
%--------------------------------------------------------------------------
% Assemble right-hand-side
%--------------------------------------------------------------------------
for i=1:3
     bi=(area.*(feval(f,center))')/3;
     b1=b1+sparse(mesh.elem(:,i),1,bi,N,1); 
end  



%assemble neumann condition
if(size(mesh.Neumann,1)~=0)
   for i=1:2
       bN=(mesh.node(mesh.Neumann(:,1),1)-mesh.node(mesh.Neumann(:,2),1)).^2+(mesh.node(mesh.Neumann(:,1),2)-mesh.node(mesh.Neumann(:,2),2)).^2;
       bN=sqrt(bN).*feval(g_N,(mesh.node(mesh.Neumann(:,1),:)+mesh.node(mesh.Neumann(:,2),:))/2)'/2;
       b1=b1+sparse(mesh.Neumann(:,i),1,bN,N,1); 
   end
end

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

%--------------------------------------------------------------------------
% Solve the linear system A * U = b1 for the free nodes
%--------------------------------------------------------------------------
freeNode = find(mesh.type>0); 
freeNode = setdiff(freeNode, bdNode);B=A(freeNode, freeNode);
if size(freeNode)>0
   u(freeNode) =bicg(A(freeNode, freeNode),b1(freeNode),0.0000001,200000);
end
mesh.solu=u(:);
% Graphic representation
%--------------------------------------------------------------------------
subplot(1,2,1)
hold off
trisurf(mesh.elem, mesh.node(:,1), mesh.node(:,2), u', ...
        'FaceColor', 'interp', 'EdgeColor', 'interp');
view(40,15);
title('Solution to the Poisson equation', 'FontSize', 14)

%--------------------------------------------------------------------------
% end of function CDFM
%--------------------------------------------------------------------------