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


function [zlsNode,zlsConn]=lsmeshT6(node,conn,phi)

% function [zlsNode,zlsConn]=lsmeshT6(node,conn,phi)
%
% This function finds the zero level set mesh for T6 elements

numNode=size(node,1);   
numElem=size(conn,1);
zlsTol=.0005;
maxIter=10;

%% find the elements cut by the level set and near nodes
cutElements=getCutElements(conn,phi);
numCutElem=length(cutElements);

cutConn=conn(cutElements,:);
xConn=node(:,1);
xConn=xConn(cutConn);
yConn=node(:,2);
yConn=yConn(cutConn);
cutPhi=phi(cutConn);

% clear('cutConn');
% clear('cutElements');

% find xi for intersection of level set with edge
a1=1/2*(cutPhi(:,2)+cutPhi(:,3))-cutPhi(:,5);
a2=1/2*(cutPhi(:,3)+cutPhi(:,1))-cutPhi(:,6);
a3=1/2*(cutPhi(:,1)+cutPhi(:,2))-cutPhi(:,4);
b1=1/2*(-cutPhi(:,2)+cutPhi(:,3));
b2=1/2*(-cutPhi(:,3)+cutPhi(:,1));
b3=1/2*(-cutPhi(:,1)+cutPhi(:,2));
c1=cutPhi(:,5);
c2=cutPhi(:,6);
c3=cutPhi(:,4);

edge1Xi1=-c1./b1;
nzs=find(a1~=0);
edge1Xi1(nzs)=(-b1(nzs)+sqrt(b1(nzs).^2-4*a1(nzs).*c1(nzs)))./(2*a1(nzs));

edge1Xi2=-c1./b1;
edge1Xi2(nzs)=(-b1(nzs)-sqrt(b1(nzs).^2-4*a1(nzs).*c1(nzs)))./(2*a1(nzs));

edge2Xi1=-c2./b2;
nzs=find(a2~=0);
edge2Xi1(nzs)=(-b2(nzs)+sqrt(b2(nzs).^2-4*a2(nzs).*c2(nzs)))./(2*a2(nzs));

edge2Xi2=-c2./b2;
edge2Xi2(nzs)=(-b2(nzs)-sqrt(b2(nzs).^2-4*a2(nzs).*c2(nzs)))./(2*a2(nzs));

edge3Xi1=-c3./b3;
nzs=find(a3~=0);
edge3Xi1(nzs)=(-b3(nzs)+sqrt(b3(nzs).^2-4*a3(nzs).*c3(nzs)))./(2*a3(nzs));

edge3Xi2=-c3./b3;
edge3Xi2(nzs)=(-b3(nzs)-sqrt(b3(nzs).^2-4*a3(nzs).*c3(nzs)))./(2*a3(nzs));

% pick right quadratic
edge1Xi1( find(imag(edge1Xi1)) ) = 1000;
edge1Xi2( find(imag(edge1Xi2)) ) = 1000;
goodQuadElems=intersect(find(edge1Xi1 >= -1 ),find(edge1Xi1 < 1 ));
edge1Xi=edge1Xi2;
edge1Xi(goodQuadElems)=edge1Xi1(goodQuadElems);

edge2Xi1( find(imag(edge2Xi1)) ) = 1000;
edge2Xi2( find(imag(edge2Xi2)) ) = 1000;
goodQuadElems=intersect(find(edge2Xi1 >= -1 ),find(edge2Xi1 < 1 ));
edge2Xi=edge2Xi2;
edge2Xi(goodQuadElems)=edge2Xi1(goodQuadElems);

edge3Xi1( find(imag(edge3Xi1)) ) = 1000;
edge3Xi2( find(imag(edge3Xi2)) ) = 1000;
goodQuadElems=intersect(find(edge3Xi1 >= -1 ),find(edge3Xi1 < 1 ));
edge3Xi=edge3Xi2;
edge3Xi(goodQuadElems)=edge3Xi1(goodQuadElems);

% get x and y coords of edge intersection
edge1X=-xConn(:,2).*(1-edge1Xi).*edge1Xi/2+...
     xConn(:,3).*(1+edge1Xi).*edge1Xi/2-xConn(:,5).*(edge1Xi+1).*(edge1Xi-1);
edge1Y=-yConn(:,2).*(1-edge1Xi).*edge1Xi/2+...
     yConn(:,3).*(1+edge1Xi).*edge1Xi/2-yConn(:,5).*(edge1Xi+1).*(edge1Xi-1);
   
edge2X=-xConn(:,3).*(1-edge2Xi).*edge2Xi/2+...
     xConn(:,1).*(1+edge2Xi).*edge2Xi/2-xConn(:,6).*(edge2Xi+1).*(edge2Xi-1);
edge2Y=-yConn(:,3).*(1-edge2Xi).*edge2Xi/2+...
     yConn(:,1).*(1+edge2Xi).*edge2Xi/2-yConn(:,6).*(edge2Xi+1).*(edge2Xi-1);

edge3X=-xConn(:,1).*(1-edge3Xi).*edge3Xi/2+...
     xConn(:,2).*(1+edge3Xi).*edge3Xi/2-xConn(:,4).*(edge3Xi+1).*(edge3Xi-1);
edge3Y=-yConn(:,1).*(1-edge3Xi).*edge3Xi/2+...
     yConn(:,2).*(1+edge3Xi).*edge3Xi/2-yConn(:,4).*(edge3Xi+1).*(edge3Xi-1);

plot_mesh(node,conn(cutElements,:),'T6','g')
hold on
plot(edge1X,edge1Y,'r.')
hold on
plot(edge2X,edge2Y,'b.')
hold on
plot(edge3X,edge3Y,'m.')
   
% filter out xis such that -1 <= xi < 1 ( edge is not intersected )
goodEdge1Elems=intersect(find(edge1Xi >= -1 ),find(edge1Xi < 1 ));
goodEdge2Elems=intersect(find(edge2Xi >= -1 ),find(edge2Xi < 1 ));
goodEdge3Elems=intersect(find(edge3Xi >= -1 ),find(edge3Xi < 1 ));

% clear('edge1Xi');
% clear('edge2Xi');
% clear('edge3Xi');
% clear('edge1Xi1');
% clear('edge2Xi1');
% clear('edge3Xi1');
% clear('edge1Xi2');
% clear('edge2Xi2');
% clear('edge3Xi2');

% initial zlsNode
zlsNode=(-1e10)*ones(3*numCutElem,2);
zlsNode(3*goodEdge1Elems-2,1)=edge1X(goodEdge1Elems);
zlsNode(3*goodEdge1Elems-2,2)=edge1Y(goodEdge1Elems);
zlsNode(3*goodEdge2Elems-1,1)=edge2X(goodEdge2Elems);
zlsNode(3*goodEdge2Elems-1,2)=edge2Y(goodEdge2Elems);
zlsNode(3*goodEdge3Elems,1)=edge3X(goodEdge3Elems);
zlsNode(3*goodEdge3Elems,2)=edge3Y(goodEdge3Elems);

goodNodes=find(zlsNode(:,1)~=(-1e10));
zlsNode=zlsNode(goodNodes,:);

% initial zls connectivity
zlsConn=[1:2:2*numCutElem-1;2:2:2*numCutElem]';

% eliminate redundant nodes
[zlsNode,inverseMap,newNumberMap]=unique(zlsNode,'rows');
zlsConn=newNumberMap(zlsConn);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find mid-point nodes on zls
% get initial guess
xNode=zlsNode(:,1);
yNode=zlsNode(:,2);
xConn=xNode(zlsConn);
yConn=yNode(zlsConn);
p=1/2*[ xConn(:,1)+xConn(:,2), yConn(:,1)+yConn(:,2) ];
iterAt=(1:size(p,1))';

err=1000;
errors=err*ones(size(p,2),1);
cnt=0;
while ( err >= zlsTol & cnt <= maxIter ) % iteration loop
  
%   clf
%   plot_mesh(node,cutConn,'T6','b')
%   hold on
%   plot_mesh(zlsNode,zlsConn,'L2','k-')
%   hold on
%   plot(p(iterAt,1),p(iterAt,2),'r.')
  
  [es,xis]=tsearchn( node, conn(:,1:3), p(iterAt,:) );
  if ( length(iterAt)==1 )
    xis=xis';
  end
  
  for n=1:length(iterAt)
    e=es(n);
    sctr=conn(e,:);
    phis=phi(sctr);
    
    pt=xis(n,2:3);
    [N,dNdxi]=lagrange_basis('T6',pt);
    J=node(sctr,:)'*dNdxi;
    dNdx=dNdxi*inv(J);
    
    phiPt=N'*phis;
    gradPhi=phis'*dNdx;
    
    delta=-phiPt*gradPhi/norm(gradPhi);
    p(iterAt(n),:)=p(iterAt(n),:)+delta;
    errors(iterAt(n))=phiPt;
    
    if ( abs(phiPt) <= zlsTol )
      iterAt(n)=0;
    end
    
  end
  
  iterAt=setdiff(iterAt,0);
  
  %length(iterAt)
  err=max(abs(errors));
  cnt=cnt+1;
end

nn=size(zlsNode,1);
ee=size(zlsConn,1);
zlsNode=[zlsNode;p];
zlsConn=[zlsConn,(1:ee)'+nn];