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

function F=projectInterfaceSpeed(node,conn,phi,zlsNode,speed)

% function F=projectInterfaceSpeed(node,conn,phi,zlsNode,speed)
%
% This funciton projects a interfacial speed defined on an iterface to
% the rest of the nodes in the computational domain.

numnode=size(node,1);
numelem=size(conn,1);

% GET SYSTEM MATRIX A
A=sparse(numnode,numnode);
for e=1:numelem
  
  sctr=conn(e,:);
  [W,Q]=quadrature( 1, 'TRIANGULAR', 2 );
  
  for q=1:length(W)
    
    pt=Q(q,:);
    wt=W(q);
    [N,dNdxi]=lagrange_basis('T3',pt);
    J=node(sctr,:)'*dNdxi;
    detJ=det(J); 
    dNdx=dNdxi*inv(J);
    
    phiPt=N'*phi(sctr);
    gradPhi=dNdx'*phi(sctr);
    tau=detJ/(gradPhi'*gradPhi);
    
    % add gradPhi . gradF operator
    A(sctr,sctr)=A(sctr,sctr)+N*(gradPhi'*dNdx')*wt*detJ;
    % add SUPG operator
    A(sctr,sctr)=A(sctr,sctr)+(dNdx*gradPhi)*tau*(gradPhi'*dNdx')*wt*detJ;
    
  end % of quadrature loop
  
end % of element loop

% GET CONSTRIANT MATRIX G
nn=size(zlsNode,1);
G=sparse(nn,numnode);
l=zeros(nn,1);

% get elements zlsNodes are in
[es,xis] = tsearchn(node,conn(:,1:3),zlsNode);
edgeLookUp=[2 3;3 1;1 2];

% loop over zleNodes
for n=1:nn
  
  % find edge cut by zls
  e=es(n);
  xi=xis(n,:);
  [val,ed]=min(xi);
  [ns]=edgeLookUp(ed,:);
  n1=conn(e,ns(1));
  n2=conn(e,ns(2));
  
  phi1=phi(n1);
  phi2=phi(n2);
  
  G(n,[n1 n2])=G(n,[n1 n2])+[phi2 -phi1]/(phi2-phi1);
  l(n)=speed(n);
  
end % loop over zlsNodes

% solve system
beta=(10e3)*trace(A)/numnode;
A=[A+beta*G'*G             G';
             G  sparse(nn,nn)];
l=[beta*G'*l; l];

F=A\l;
F=F(1:numnode);