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.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