distmeshnd

DISTMESHND N-D Mesh Generator using Distance Functions.
   [P,T]=DISTMESHND(FDIST,FH,H,BOX,FIX,FDISTPARAMS)

      P:           Node positions (NxNDIM)
      T:           Triangle indices (NTx(NDIM+1))
      FDIST:       Distance function
      FH:          Edge length function
      H:           Smallest edge length
      BOX:         Bounding box [xmin,xmax;ymin,ymax; ...] (NDIMx2)
      FIX:         Fixed node positions (NFIXxNDIM)
      FDISTPARAMS: Additional parameters passed to FDIST

   Example: Unit ball
      dim=3;
      d=inline('sqrt(sum(p.^2,2))-1','p');
      [p,t]=distmeshnd(d,@huniform,0.2,[-ones(1,dim);ones(1,dim)],[]);

   See also: DISTMESH2D, DELAUNAYN, TRIMESH, MESHDEMOND.

0001 function [p,t]=distmeshnd(fdist,fh,h,box,fix,varargin)
0002 %DISTMESHND N-D Mesh Generator using Distance Functions.
0003 %   [P,T]=DISTMESHND(FDIST,FH,H,BOX,FIX,FDISTPARAMS)
0004 %
0005 %      P:           Node positions (NxNDIM)
0006 %      T:           Triangle indices (NTx(NDIM+1))
0007 %      FDIST:       Distance function
0008 %      FH:          Edge length function
0009 %      H:           Smallest edge length
0010 %      BOX:         Bounding box [xmin,xmax;ymin,ymax; ...] (NDIMx2)
0011 %      FIX:         Fixed node positions (NFIXxNDIM)
0012 %      FDISTPARAMS: Additional parameters passed to FDIST
0013 %
0014 %   Example: Unit ball
0015 %      dim=3;
0016 %      d=inline('sqrt(sum(p.^2,2))-1','p');
0017 %      [p,t]=distmeshnd(d,@huniform,0.2,[-ones(1,dim);ones(1,dim)],[]);
0018 %
0019 %   See also: DISTMESH2D, DELAUNAYN, TRIMESH, MESHDEMOND.
0020 
0021 %   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.
0022 % This file is taken from the 3.2 Release by Per-Olof Persson,
0023 % available under the GPL version 2 or any later version.
0024 % Source is from http://www-math.mit.edu/~persson/mesh/
0025 
0026 global distmesh_do_graphics; % flag do decide if we do graphics
0027 maxiter = 500;
0028 
0029 dim=size(box,2);
0030 ptol=.01;
0031 ttol=.1;
0032 L0mult=1+.4/2^(dim-1);
0033 deltat=.10; % Stimulation Time step - Speed or slow simulation
0034 geps=1e-1*h; deps=sqrt(eps)*h;
0035 
0036 % 1. Create initial distribution in bounding box
0037 if dim==1
0038   p=(box(1):h:box(2))';
0039 else
0040   cbox=cell(1,dim);
0041   for ii=1:dim
0042     cbox{ii}=box(1,ii):h:box(2,ii);
0043   end
0044   pp=cell(1,dim);
0045   [pp{:}]=ndgrid(cbox{:});
0046   for ii=1:dim; pp{ii} = pp{ii} + h/2*(rand(size(pp{ii}))-0.5); end
0047   p=zeros(prod(size(pp{1})),dim);
0048   for ii=1:dim
0049     p(:,ii)=pp{ii}(:);
0050   end
0051 end
0052 
0053 % 2. Remove points outside the region, apply the rejection method
0054 p=p(feval(fdist,p,varargin{:})<geps,:);
0055 r0=feval(fh,p,varargin{:});
0056 p=[fix; p(rand(size(p,1),1)<min(r0)^dim./r0.^dim,:)];
0057 N=size(p,1);
0058 
0059 count=0;
0060 p0=inf;
0061 while 1
0062   % 3. Retriangulation by Delaunay
0063   if max(sqrt(sum((p-p0).^2,2)))>ttol*h
0064     p0=p;
0065     t=delaunayn(p);
0066     pmid=zeros(size(t,1),dim);
0067     for ii=1:dim+1
0068       pmid=pmid+p(t(:,ii),:)/(dim+1);
0069     end
0070     t=t(feval(fdist,pmid,varargin{:})<geps,:); %BUGFIX -geps to geps
0071     % 4. Describe each edge by a unique pair of nodes
0072     pair=zeros(0,2);
0073     localpairs=nchoosek(1:dim+1,2);
0074     for ii=1:size(localpairs,1)
0075       pair=[pair;t(:,localpairs(ii,:))];
0076     end
0077     pair=unique(sort(pair,2),'rows');
0078     % 5. Graphical output of the current mesh
0079     if distmesh_do_graphics
0080        if dim==2
0081          hh=trimesh(t,p(:,1),p(:,2),zeros(N,1),'edgecolor','black');
0082          view(2),axis equal,axis off,drawnow
0083        elseif dim==3
0084          if mod(count,5)==0
0085            simpplot(p,t,'p(:,2)>0');
0086            title(['Retriangulation #',int2str(count)])
0087            drawnow
0088          end
0089       end
0090     else
0091       if count ==0
0092          opt.final_msg = '';
0093          progress_msg('distmeshng: retriangulation #:',0,0,opt);
0094       else
0095          progress_msg(count,0);
0096       end
0097     end
0098     count=count+1;
0099   end
0100 
0101   
0102   % 6. Move mesh points based on edge lengths L and forces F
0103   bars=p(pair(:,1),:)-p(pair(:,2),:);
0104   L=sqrt(sum(bars.^2,2));
0105   L0=feval(fh,(p(pair(:,1),:)+p(pair(:,2),:))/2,varargin{:});
0106   L0=L0*L0mult*(sum(L.^dim)/sum(L0.^dim))^(1/dim);
0107   F=max(L0-L,0);
0108   Fbar=[bars,-bars].*repmat(F./L,1,2*dim);
0109   dp=full(sparse(pair(:,[ones(1,dim),2*ones(1,dim)]), ...
0110                  ones(size(pair,1),1)*[1:dim,1:dim], ...
0111                  Fbar,N,dim));
0112   dp(1:size(fix,1),:)=0;
0113   p=p+deltat*dp;
0114 
0115   % 7. Bring outside points back to the boundary
0116   d=feval(fdist,p,varargin{:}); ix=d>0;
0117   ix(1:size(fix,1))= 0; % EDIT - Don't move fixed points
0118   gradd=zeros(sum(ix),dim);
0119   for ii=1:dim
0120     a=zeros(1,dim);
0121     a(ii)=deps;
0122     d1x=feval(fdist,p(ix,:)+ones(sum(ix),1)*a,varargin{:});
0123     gradd(:,ii)=(d1x-d(ix))/deps;
0124   end
0125   p(ix,:)=p(ix,:)-d(ix)*ones(1,dim).*gradd;
0126 
0127   % 8. Termination criterion
0128   maxdp=max(deltat*sqrt(sum(dp(d<-geps,:).^2,2)));
0129   if maxdp<ptol*h; break; end
0130   % maximum iterations
0131   if count>maxiter; break; end
0132 end
0133 if ~distmesh_do_graphics
0134     progress_msg(Inf);
0135 end
0136 
0137 % final delaunayn just to make sure triangularization returned is good
0138 t=delaunayn(p); % New code - AA
0139 pmid=zeros(size(t,1),dim);
0140 for ii=1:dim+1
0141   pmid=pmid+p(t(:,ii),:)/(dim+1);
0142 end
0143 t=t(feval(fdist,pmid,varargin{:})<geps,:); %BUGFIX -geps to geps