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


function a=mapmesh2D(e1,e2,e3,e4)

% MAPMESH(e1,e2,e3,e4) Creates a mapped mesh along 2D surface
%
%  

if ( size(e1,1) ~= size(e3,1) |  size(e2,1) ~= size(e4,1) )
  disp('ERROR: opposite edges must be of same length')
end

nx=size(e1,1);
ny=size(e2,1);

if ( nargin == 3 )   % TRIANGULAR SURFACE

  

else                 % RULED SURFACE  
   
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % interpolate cubic serindipity mapping nodes from boundary data
  pt1=[e1(1,:)+e4(1,:)]/2;
  pt2=[e1(nx,:)+e2(1,:)]/2;
  pt3=[e2(ny,:)+e3(nx,:)]/2;
  pt4=[e3(1,:)+e4(ny,:)]/2;

  % on lower edge
  lseg=zeros(nx,1);
  for i=2:nx
    lseg(i)=lseg(i-1)+norm(e1(i,:)-e1(i-1,:));
  end

  L1=lseg(nx)/3; L2=2*lseg(nx)/3;
  temp=[interp1(lseg,e1(:,1),[L1,L2],'cubic');
        interp1(lseg,e1(:,2),[L1,L2],'cubic')]';   

  pt5=temp(1,:);  pt6=temp(2,:);
  E1=lseg*2/lseg(nx)-1;   % boundary nodes in parent domain

  % on right edge
  lseg=zeros(ny,1);
  for i=2:ny
    lseg(i)=lseg(i-1)+norm(e2(i,:)-e2(i-1,:));
  end

  L1=lseg(ny)/3; L2=2*lseg(ny)/3;
  temp=[interp1(lseg,e2(:,1),[L1,L2],'cubic');
        interp1(lseg,e2(:,2),[L1,L2],'cubic')]';   

  pt7=temp(1,:);  pt8=temp(2,:);
  E2=lseg*2/lseg(ny)-1;   % boundary nodes in parent domain

  % on top edge
  lseg=zeros(nx,1);
  for i=2:nx
    lseg(i)=lseg(i-1)+norm(e3(i,:)-e3(i-1,:));
  end

  L1=lseg(nx)/3; L2=2*lseg(nx)/3;
  temp=[interp1(lseg,e3(:,1),[L1,L2],'cubic');
        interp1(lseg,e3(:,2),[L1,L2],'cubic')]';   

  pt10=temp(1,:);  pt9=temp(2,:);
  E3=lseg*2/lseg(nx)-1;   % boundary nodes in parent domain    

  % on left edge
  lseg=zeros(ny,1);
  for i=2:ny
    lseg(i)=lseg(i-1)+norm(e4(i,:)-e4(i-1,:));
  end

  L1=lseg(ny)/3; L2=2*lseg(ny)/3;
  temp=[interp1(lseg,e4(:,1),[L1,L2],'cubic');
        interp1(lseg,e4(:,2),[L1,L2],'cubic')]';   

  pt12=temp(1,:);  pt11=temp(2,:);
  E4=lseg*2/lseg(ny)-1;   % boundary nodes in parent domain 

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % setup interior node grid on parent domain
  x=zeros(nx,ny); y=zeros(nx,ny);

  X=[pt1(1); pt2(1); pt3(1); pt4(1); pt5(1); pt6(1);
      pt7(1); pt8(1); pt9(1); pt10(1); pt11(1); pt12(1)]; 
 
  Y=[pt1(2); pt2(2); pt3(2); pt4(2); pt5(2); pt6(2);
      pt7(2); pt8(2); pt9(2); pt10(2); pt11(2); pt12(2)]; 


  for i=2:nx-1
    for j=2:ny-1

      A=[E4(j)-E2(j) 2;2 E1(i)-E3(i)];
      b=[E2(j)+E4(j);E1(i)+E3(i)];
      inter=(A\b);
      xi=inter(1); eta=inter(2);
   
      N(1)=(1-xi)*(1-eta)*(9*(xi^2+eta^2)-10)/32;     % mapping functions
      N(2)=(1+xi)*(1-eta)*(9*(xi^2+eta^2)-10)/32; 
      N(3)=(1+xi)*(1+eta)*(9*(xi^2+eta^2)-10)/32;  
      N(4)=(1-xi)*(1+eta)*(9*(xi^2+eta^2)-10)/32;

      N(5)=9/32*(1-eta)*(1-xi^2)*(1-3*xi);
      N(6)=9/32*(1-eta)*(1-xi^2)*(1+3*xi);
       
      N(7)=9/32*(1+xi)*(1-eta^2)*(1-3*eta);
      N(8)=9/32*(1+xi)*(1-eta^2)*(1+3*eta);
  
      N(9)= 9/32*(1+eta)*(1-xi^2)*(1+3*xi);
      N(10)=9/32*(1+eta)*(1-xi^2)*(1-3*xi); 
       
      N(11)=9/32*(1-xi)*(1-eta^2)*(1+3*eta);
      N(12)=9/32*(1-xi)*(1-eta^2)*(1-3*eta);        

      x(i,j)=N*X; y(i,j)=N*Y;
      
    end
  end

  x(:,1)  = e1(:,1);
  x(nx,:) = e2(:,1)';
  x(:,ny) = e3(:,1);
  x(1,:)  = e4(:,1)';

  y(:,1)  = e1(:,2);
  y(nx,:) = e2(:,2)';
  y(:,ny) = e3(:,2);
  y(1,:)  = e4(:,2)';
   
  a=zeros(nx*ny,2);
  for j=1:ny
    for i=1:nx
      a((j-1)*nx+i,1)=x(i,j); 
      a((j-1)*nx+i,2)=y(i,j);
    end
  end 

end