www.pudn.com > 87361026FEA.rar > lsupdateT6.m, change:2002-10-01,size:1596b


function [phi,M]=lsupdateT6(node,conn,phi,F,delt,M)

% function phi=lsupdateT6(node,conn,phi,F,delt)
% function [phi,M]=lsupdateT6(node,conn,phi,F,delt)
% function phi=lsupdateT6(node,conn,phi,F,delt,M)
%
% This function performs an update of a level set field phi on a finite
% element descretization using a Taylor-Galerkin stabilized finite element
% scheme.  Node is the nodal coordinate matrix, conn is the element 
% connectivity matrix, F is the nodal speed function, delt is the 
% time step and elemType it the element type.

% INITIALIZE DATA STRUCTURES
numNode=size(node,1);
numElem=size(conn,1);

% COMPUTE M and fint
fint=zeros(numNode,1);

if ( nargin==5 )
  M=sparse(numNode,numNode);
  [W,Q]=quadrature(4,'TRIANGULAR',2);
  flag=1;
else
  [W,Q]=quadrature(4,'TRIANGULAR',2);
  flag=0;
end

if ( size(F,2) == 1 )
  speedFlag=1;
else
  speedFlag=0;
end

for e=1:numElem

  sctr=conn(e,:); 
   
  for q=1:size(W)
    
    [N,dNdxi]=lagrange_basis('T6',Q(q,:));
    J=node(sctr,:)'*dNdxi;
    detJ=det(J);
    dNdx=dNdxi*inv(J);
    
    gradPhi=phi(sctr)'*dNdx;
    normal=gradPhi/norm(gradPhi+.00001);
    dNdn=dNdx*normal';
    dPhidn=phi(sctr)'*dNdn;
    
    Fpt=N'*F(sctr);
    
    % advection part
    fint(sctr) = fint(sctr) - N*Fpt*norm(gradPhi)*W(q)*detJ;  
    
    % stabilization part
    fint(sctr) = fint(sctr) - delt/2*dNdn*Fpt^2*dPhidn'*W(q)*detJ;
    
    % lumped mass
    if ( flag )
      M(sctr,sctr) = M(sctr,sctr) + N*N'*W(q)*detJ*flag;
    end
    
  end % of quadrature loop
  
end % of element loop

% UPDATE LEVEL SET
dphi=delt*(M\fint);
phi = phi + dphi;