cdfm.m

C++编译指导

function u =CDFM(mesh, eps,dta,a1, a2, b, f, g_D)
% 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); 

%--------------------------------------------------------------------------
% 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;

%--------------------------------------------------------------------------
% 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))...
           +(-feval(a1,center).*ve(:,2,i)+feval(a2,center).*ve(:,1,i))/6...
                                              +(area.*feval(b,center))/9...
+(feval(dta,center).*(-feval(a1,center).*ve(:,2,i)+feval(a2,center)...
.*ve(:,1,j)).*(-feval(a1,center).*ve(:,2,j)+feval(a2,center).*ve(:,1,j)))...
+(feval(dta,center).*feval(b,center).*(-feval(a1,center).*ve(:,2,i)...
                                       +feval(a2,center).*ve(:,1,i)))/6...
+(feval(dta,center).*feval(b,center).*(-feval(a1,center).*ve(:,2,j)...
                                       +feval(a2,center).*ve(:,1,j)))/6...
           +(area.*feval(b,center).*feval(b,center).*feval(dta,center))/9;   
   A = A + sparse(mesh.elem(:,i),mesh.elem(:,j),Aij,N,N);
    end
end
%--------------------------------------------------------------------------
% Assemble right-hand-side
%--------------------------------------------------------------------------
for i=1:3
    bi=(area.*feval(f,center))/3 + (feval(dta,center).*feval(f,center)...
        .*(-feval(a1,center).*ve(:,2,i)+feval(a2,center).*ve(:,1,i)))/2...
        +(area.*feval(dta,center).*feval(b,center).*feval(f,center))/3;
    b1=b1+sparse(mesh.elem(:,i),bi,N) ; 
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);
if size(freeNode)>0
   u(freeNode) = A(freeNode, freeNode) \ b1(freeNode);
end

% 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
%--------------------------------------------------------------------------