distmesh2d

PURPOSE ^

DISTMESH2D 2-D Mesh Generator using Distance Functions.

SYNOPSIS ^

function [p,t]=distmesh2d(fd,fh,h0,bbox,pfix,varargin)

DESCRIPTION ^

DISTMESH2D 2-D Mesh Generator using Distance Functions.
   [P,T]=DISTMESH2D(FD,FH,H0,BBOX,PFIX,FPARAMS)

      P:         Node positions (Nx2)
      T:         Triangle indices (NTx3)
      FD:        Distance function d(x,y)
      FH:        Scaled edge length function h(x,y)
      H0:        Initial edge length
      BBOX:      Bounding box [xmin,ymin; xmax,ymax]
      PFIX:      Fixed node positions (NFIXx2)
      FPARAMS:   Additional parameters passed to FD and FH

   Example: (Uniform Mesh on Unit Circle)
      fd=inline('sqrt(sum(p.^2,2))-1','p');
      [p,t]=distmesh2d(fd,@huniform,0.2,[-1,-1;1,1],[]);

   Example: (Rectangle with circular hole, refined at circle boundary)
      fd=inline('ddiff(drectangle(p,-1,1,-1,1),dcircle(p,0,0,0.5))','p');
      fh=inline('min(4*sqrt(sum(p.^2,2))-1,2)','p');
      [p,t]=distmesh2d(fd,fh,0.05,[-1,-1;1,1],[-1,-1;-1,1;1,-1;1,1]);

   See also: MESHDEMO2D, DISTMESHND, DELAUNAYN, TRIMESH.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [p,t]=distmesh2d(fd,fh,h0,bbox,pfix,varargin)
0002 %DISTMESH2D 2-D Mesh Generator using Distance Functions.
0003 %   [P,T]=DISTMESH2D(FD,FH,H0,BBOX,PFIX,FPARAMS)
0004 %
0005 %      P:         Node positions (Nx2)
0006 %      T:         Triangle indices (NTx3)
0007 %      FD:        Distance function d(x,y)
0008 %      FH:        Scaled edge length function h(x,y)
0009 %      H0:        Initial edge length
0010 %      BBOX:      Bounding box [xmin,ymin; xmax,ymax]
0011 %      PFIX:      Fixed node positions (NFIXx2)
0012 %      FPARAMS:   Additional parameters passed to FD and FH
0013 %
0014 %   Example: (Uniform Mesh on Unit Circle)
0015 %      fd=inline('sqrt(sum(p.^2,2))-1','p');
0016 %      [p,t]=distmesh2d(fd,@huniform,0.2,[-1,-1;1,1],[]);
0017 %
0018 %   Example: (Rectangle with circular hole, refined at circle boundary)
0019 %      fd=inline('ddiff(drectangle(p,-1,1,-1,1),dcircle(p,0,0,0.5))','p');
0020 %      fh=inline('min(4*sqrt(sum(p.^2,2))-1,2)','p');
0021 %      [p,t]=distmesh2d(fd,fh,0.05,[-1,-1;1,1],[-1,-1;-1,1;1,-1;1,1]);
0022 %
0023 %   See also: MESHDEMO2D, DISTMESHND, DELAUNAYN, TRIMESH.
0024 
0025 %   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.
0026 % This file is taken from the 3.2 Release by Per-Olof Persson,
0027 % available under the GPL version 2 or any later version.
0028 % Source is from http://www-math.mit.edu/~persson/mesh/
0029 
0030 dptol=.01; ttol=.1; Fscale=1.2; deltat=.2; geps=.001*h0; deps=sqrt(eps)*h0;
0031 
0032 % 1. Create initial distribution in bounding box (equilateral triangles)
0033 [x,y]=meshgrid(bbox(1,1):h0:bbox(2,1),bbox(1,2):h0*sqrt(3)/2:bbox(2,2));
0034 x(2:2:end,:)=x(2:2:end,:)+h0/2;                      % Shift even rows
0035 p=[x(:),y(:)];                                       % List of node coordinates
0036 
0037 % 2. Remove points outside the region, apply the rejection method
0038 p=p(feval(fd,p,varargin{:})<geps,:);                 % Keep only d<0 points
0039 r0=1./feval(fh,p,varargin{:}).^2;                    % Probability to keep point
0040 p=[pfix; p(rand(size(p,1),1)<r0./max(r0),:)];        % Rejection method
0041 
0042 pold=inf;                                            % For first iteration
0043 iteration=1;
0044 while 1
0045   N=size(p,1);                                       % Number of points N
0046   % 3. Retriangulation by the Delaunay algorithm
0047   if max(sqrt(sum((p-pold).^2,2))/h0)>ttol           % Any large movement?
0048     pold=p;                                          % Save current positions
0049     t=delaunayn(p);                                  % List of triangles
0050     pmid=(p(t(:,1),:)+p(t(:,2),:)+p(t(:,3),:))/3;    % Compute centroids
0051     t=t(feval(fd,pmid,varargin{:})<-geps,:);         % Keep interior triangles
0052     % 4. Describe each bar by a unique pair of nodes
0053     bars=[t(:,[1,2]);t(:,[1,3]);t(:,[2,3])];         % Interior bars duplicated
0054     [bars,i_bars,j_bars]=unique(bars,'rows');% Bars as node pairs
0055 %   [bars,i_bars,j_bars]=unique(sort(bars,2),'rows');% Bars as node pairs
0056     % 5. Graphical output of the current mesh
0057     trimesh(t,p(:,1),p(:,2),zeros(N,1),'edgecolor','black')
0058     view(2),axis equal,axis off,drawnow
0059   end
0060 
0061   % 6. Move mesh points based on bar lengths L and forces F
0062   barvec=p(bars(:,1),:)-p(bars(:,2),:);              % List of bar vectors
0063   L=sqrt(sum(barvec.^2,2));                          % L = Bar lengths
0064 
0065   if 1
0066 % Get length function at average
0067      hbars=feval(fh,(p(bars(:,1),:)+p(bars(:,2),:))/2,varargin{:});
0068   else
0069 % Get minimum length function at both ends
0070      hbars1=feval(fh,p(bars(:,1),:),varargin{:});
0071      hbars2=feval(fh,p(bars(:,2),:),varargin{:});
0072      hbars= min([hbars1,hbars2],[],2);
0073   end
0074 
0075   L0=hbars*Fscale*sqrt(sum(L.^2)/sum(hbars.^2));     % L0 = Desired lengths
0076 
0077   add_p= zeros(0,2);
0078 if 0
0079   % 6.3 Desired length cannot be larger than tshape * connecting ones
0080   nt= size(t,1); % number of triangles
0081   min_L= zeros(nt,1);
0082   for i= 1:nt
0083      Lt = L(j_bars([i,i+nt,i+nt*2]));
0084      min_L(i) = min(Lt);
0085   
0086      if min(Lt)*5 < max(Lt); add_p= [add_p;pmid(i,:)];end
0087   end
0088   min_L = 2*[min_L;min_L;min_L];
0089   L0i= L0;
0090   L0 = min([L0,min_L(i_bars)],[],2);
0091 end
0092 
0093   F=max(L0-L,0);                                     % Bar forces (scalars)
0094   Fvec=F./L*[1,1].*barvec;                           % Bar forces (x,y components)
0095 
0096   % 6.5. Points get sucked to the centroid to make them equilateral
0097 if 0
0098   for i= 1:size(t,1)
0099      idx= t(i,:);
0100      dist_vtx = p(idx,:) - ones(3,1)*pmid(i,:);
0101      dist_mn  = sqrt(sum(dist_vtx.^2,2));
0102      dist_mn  = dist_mn/mean(dist_mn) - 1;
0103      dist_mn  = dist_mn.*(abs(dist_mn)>.4); %zero if too small
0104      F_adj    = (dist_mn*[1,1]) .* dist_vtx; 
0105      Fvec(idx,:) = Fvec(idx,:) + .3*F_adj;
0106   end
0107 end
0108 
0109   Ftot=full(sparse(bars(:,[1,1,2,2]),ones(size(F))*[1,2,1,2],[Fvec,-Fvec],N,2));
0110   Ftot(1:size(pfix,1),:)=0;                          % Force = 0 at fixed points
0111   p=p+deltat*Ftot;                                   % Update node positions
0112 
0113 
0114   % 7. Bring outside points back to the boundary
0115   d=feval(fd,p,varargin{:}); ix=d>0;                 % Find points outside (d>0)
0116   dgradx=(feval(fd,[p(ix,1)+deps,p(ix,2)],varargin{:})-d(ix))/deps; % Numerical
0117   dgrady=(feval(fd,[p(ix,1),p(ix,2)+deps],varargin{:})-d(ix))/deps; %    gradient
0118   p(ix,:)=p(ix,:)-[d(ix).*dgradx,d(ix).*dgrady];     % Project back to boundary
0119 
0120   % 8. Termination criterion: All interior nodes move less than dptol (scaled)
0121   if max(sqrt(sum(deltat*Ftot(d<-geps,:).^2,2))/h0)<dptol, break; end
0122 
0123   % 9. Add new points
0124   if size(add_p,1)>0 && rem(iteration,10)==1
0125      p= [p;add_p];
0126      pold=inf; 
0127   end
0128   iteration=iteration+1;
0129 end

Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005