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


function h=getElementSizeFactor(node,conn,elemType)

% function h=getElementSizeFactor(node,conn,elemType)
% 
% Returns the element characteristic size factor.

numNode=size(node,1);
numElem=size(conn,1);

if ( size(node,2) == 1 )
  node=[node,zeros(numNode,2)];
elseif ( size(node,2) == 2 )
  node=[node,zeros(numNode,1)];
end

X=node(:,1);
Xconn=X(conn);
Y=node(:,2);
Yconn=Y(conn);
Z=node(:,3);
Zconn=Z(conn);

switch elemType
case 'L2'
  topo='line';
  factor=1;
case 'L3'
  topo='line';
  factor=0.5;
case 'T3'
  topo='tri';
  factor=1;
case 'T6'
  topo='tri';
  factor=0.5;
case 'T4'
  topo='tri';
  factor=1;
case 'Q4'
  topo='quad';
  factor=1;
case 'Q9'
  topo='quad';
  factor=0.5;
% case 'H4'
%   topo='tet';
%   factor=1.0;
% case 'H10'
%   topo='tet';
%   factor=0.5;
% case 'B8'
%   topo='brick';
%   factor=1.0;
% case 'B27'
%   topo='brick';
%   factor=0.5;
otherwise
  disp('ERROR: Element type not found in getElementSizeFactor')
  h=0;
  return
end

% Get characteristic edge lengths
switch topo

case 'line'
  h=factor*sqrt((Xconn(:,2)-Xconn(:,1)).^2+(Yconn(:,2)-Yconn(:,1)).^2+...
                                       (Zconn(:,2)-Zconn(:,1)).^2);
case 'tri'
  cls=zeros(numElem,3);
  
  mdPtX=0.5*(Xconn(:,2)+Xconn(:,3));
  mdPtY=0.5*(Yconn(:,2)+Yconn(:,3));
  mdPtZ=0.5*(Zconn(:,2)+Zconn(:,3));
  cls(:,1)=factor*sqrt((Xconn(:,1)-mdPtX).^2+(Yconn(:,1)-mdPtY).^2+...
                                               (Zconn(:,1)-mdPtZ).^2);
  
  mdPtX=0.5*(Xconn(:,3)+Xconn(:,1));
  mdPtY=0.5*(Yconn(:,3)+Yconn(:,1));
  mdPtZ=0.5*(Zconn(:,3)+Zconn(:,1));
  cls(:,2)=factor*sqrt((Xconn(:,2)-mdPtX).^2+(Yconn(:,2)-mdPtY).^2+...
                                               (Zconn(:,2)-mdPtZ).^2);
  
  mdPtX=0.5*(Xconn(:,2)+Xconn(:,1));
  mdPtY=0.5*(Yconn(:,2)+Yconn(:,1));
  mdPtZ=0.5*(Zconn(:,2)+Zconn(:,1));
  cls(:,3)=factor*sqrt((Xconn(:,3)-mdPtX).^2+(Yconn(:,3)-mdPtY).^2+...
                                                 (Zconn(:,3)-mdPtZ).^2);
  
  h=min(cls')';   
  
case 'quad'
  
  mdPt1X=0.5*(Xconn(:,1)+Xconn(:,2));
  mdPt1Y=0.5*(Yconn(:,1)+Yconn(:,2));
  mdPt1Z=0.5*(Zconn(:,1)+Zconn(:,2));
  
  mdPt2X=0.5*(Xconn(:,3)+Xconn(:,2));
  mdPt2Y=0.5*(Yconn(:,3)+Yconn(:,2));
  mdPt2Z=0.5*(Zconn(:,3)+Zconn(:,2));
  
  mdPt3X=0.5*(Xconn(:,3)+Xconn(:,4));
  mdPt3Y=0.5*(Yconn(:,3)+Yconn(:,4));
  mdPt3Z=0.5*(Zconn(:,3)+Zconn(:,4));
  
  mdPt4X=0.5*(Xconn(:,1)+Xconn(:,4));
  mdPt4Y=0.5*(Yconn(:,1)+Yconn(:,4));
  mdPt4Z=0.5*(Zconn(:,1)+Zconn(:,4));
  
  l1=factor*sqrt((mdPt1X-mdPt3X).^2+(mdPt1Y-mdPt3Y).^2+...
                                                 (mdPt1Z-mdPt3Z).^2);
  l2=factor*sqrt((mdPt2X-mdPt4X).^2+(mdPt4Y-mdPt2Y).^2+...
                                                 (mdPt2Z-mdPt4Z).^2);
                                               
  h=min(l1,l2);  
  
% case 'tet'  
%   ls=zeros(numElem,4);
%   
%   for e=1:4
%   
%     a=zeros(3,3,numElem);
%     b=zeros(3,3,numElem);
%     c=zeros(3,3,numElem);
%     
%     ls(:,e)= abs( a*Xconn(e)+b*Yconn(e)+c*Zconn(e) +d )./ ...
%                                     sqrt( a.^2 + b.^2 + c.^2);
%   end
%                            
%   h=min(cls')';  
%   
% case 'brick' 
  
end