function [ p, t ] = distmesh_nd ( fdist, fh, h, box, fix, varargin ) %*****************************************************************************80 % %% DISTMESH_ND N-D Mesh Generator using Distance Functions. % % Example: % % For the unit ball % % dim = 3; % d = inline ( 'sqrt(sum(p.^2,2))-1', 'p' ); % [ p, t ] = distmesh_nd ( d, @huniform, 0.2, [-ones(1,dim);ones(1,dim)], [] ); % % Licensing: % % (C) 2004 Per-Olof Persson. % See COPYRIGHT.TXT for details. % % Reference: % % Per-Olof Persson and Gilbert Strang, % A Simple Mesh Generator in MATLAB, % SIAM Review, % Volume 46, Number 2, June 2004, pages 329-345. % % Parameters: % % Input, FDIST: Distance function % % Input, FH: Edge length function % % Input, real H, the smallest edge length. % % Input, real BOX(NDIM,2), a bounding box for the region. % % Input, real FIX(NFIX,NDIM), the coordinates of nodes that are % required to be included in the mesh. % % Input, VARARGIN: Additional parameters passed to FDIST % % Output, P: Node positions (NxNDIM) % % Output, T: Triangle indices (NTx(NDIM+1)) % dim = size ( box, 2 ); ptol = 0.001; ttol = 0.1; L0mult = 1 + 0.4 / 2^(dim-1); deltat = 0.1; geps = 0.1 * h; deps = sqrt ( eps ) * h; % % 1. Create the initial distribution in the bounding box. % if ( dim == 1 ) p = ( box(1) : h : box(2) )'; else cbox = cell(1,dim); for ii = 1 : dim cbox{ii} = box(1,ii):h:box(2,ii); end pp = cell(1,dim); [pp{:}] = ndgrid(cbox{:}); p = zeros(prod(size(pp{1})),dim); for ii = 1 : dim p(:,ii) = pp{ii}(:); end end % % 2. Remove points outside the region, apply the rejection method % p = p ( feval(fdist,p,varargin{:})ttol*h p0=p; t = delaunayn ( p ); pmid=zeros(size(t,1),dim); for ii=1:dim+1 pmid=pmid+p(t(:,ii),:)/(dim+1); end t=t(feval(fdist,pmid,varargin{:})<-geps,:); % % 4. Describe each edge by a unique pair of nodes % pair=zeros(0,2); localpairs=nchoosek(1:dim+1,2); for ii=1:size(localpairs,1) pair=[pair;t(:,localpairs(ii,:))]; end pair=unique(sort(pair,2),'rows'); % % 5. Graphical output of the current mesh % if ( dim == 2 ) trimesh(t,p(:,1),p(:,2),zeros(N,1)) view(2),axis equal,axis off,drawnow elseif ( dim == 3 ) if mod(count,5)==0 simpplot(p,t,'p(:,2)>0'); title(['Retriangulation #',int2str(count)]) drawnow end else disp(sprintf('Retriangulation #%d',count)) end count=count+1; end % % 6. Move mesh points based on edge lengths L and forces F % bars=p(pair(:,1),:)-p(pair(:,2),:); L=sqrt(sum(bars.^2,2)); L0=feval(fh,(p(pair(:,1),:)+p(pair(:,2),:))/2); L0=L0*L0mult*(sum(L.^dim)/sum(L0.^dim))^(1/dim); F=max(L0-L,0); Fbar=[bars,-bars].*repmat(F./L,1,2*dim); dp=full(sparse(pair(:,[ones(1,dim),2*ones(1,dim)]), ... ones(size(pair,1),1)*[1:dim,1:dim], ... Fbar,N,dim)); dp(1:size(fix,1),:)=0; p=p+deltat*dp; % % 7. Bring outside points back to the boundary % d=feval(fdist,p,varargin{:}); ix=d>0; gradd=zeros(sum(ix),dim); for ii=1:dim a=zeros(1,dim); a(ii)=deps; d1x=feval(fdist,p(ix,:)+ones(sum(ix),1)*a,varargin{:}); gradd(:,ii)=(d1x-d(ix))/deps; end p(ix,:)=p(ix,:)-d(ix)*ones(1,dim).*gradd; % % 8. Termination criterion % maxdp=max(deltat*sqrt(sum(dp(d<-geps,:).^2,2))); if maxdp < ptol * h, break; end end return end