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


function F=projectInterfaceSpeedT6(node,conn,phi,zlsNode,speed,betaFact)

% function F=projectInterfaceSpeedT6(node,conn,phi,zlsNode,speed,penalty)
%
% 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);

if ( nargin == 5 )
  betaFact=10e3;
end

% GET SYSTEM MATRIX A
A=sparse(numnode,numnode);
for e=1:numelem
  
  sctr=conn(e,:);
  [W,Q]=quadrature( 5, 'TRIANGULAR', 2 );
  
  for q=1:length(W)
    
    pt=Q(q,:);
    wt=W(q);
    [N,dNdxi]=lagrange_basis('T6',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);
xis=xis(:,2:3);

% loop over zleNodes
for n=1:nn
  
  % find edge cut by zls
  e=es(n);
  sctr=conn(e,:);
  phis=phi(sctr);
  pt=xis(n,:);
  N=lagrange_basis('T6',pt);
  
  gn=N';
  
  G(n,sctr)=G(n,sctr)+gn; 
  l(n)=speed(n);
  
end % loop over zlsNodes

% solve system
beta=betaFact*trace(A)/numnode;
% A=[A+beta*G'*G             G';  % Augmented Lagrangian formulation
%              G  sparse(nn,nn)];
% l=[beta*G'*l; l];

A=A+beta*G'*G;  % penalty method formulation
l=beta*G'*l;

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