From: James Ramm on
"John D'Errico" <woodchips(a)rochester.rr.com> wrote in message <ef26583.0(a)webx.raydaftYaTP>...
> Mayur wrote:
> >
> >
> > Hello,
> >
> > I am generating a Tetrahedral mesh and I am using delaunay3
> > function
> > for this. I do it as follows
> >
> > x=[(0):1/(numx):(1)]; % numx=numy=numz = 1 (typical cube)
> > y=[(0):1/(numy):(1)];
> > z=[(0):1/(numz):(1)];
> >
> > [X,Y,Z] = meshgrid(x,y,z);
> > X = [X(:)];
> > Y = [Y(:)];
> > Z = [Z(:)];
> > Tes = delaunay3(X,Y,Z);
> > L = [X(:) Y(:) Z(:)];
> > tetramesh1(Tes,L);camorbit(20,0)
> >
> > This gives me a tetrahedral mesh with 'Tes' having the nodal
> > connectivity matrix for a typical case (a cube with 6-tetrahedra, 8
> > nodes)it looks like
> > Tes =
> > 4 7 8 6
> > 1 5 7 6
> > 2 4 7 6
> > 2 1 7 6
> > 2 4 7 3
> > 2 1 7 3
> >
> > Now my problem is: I want to know how is the above numbering
> > associated with the individual tetrahedra. As when numx=numy=numz=
> > 2
> > it results in a 48 tetrahedra mesh and I want to know the boundary
> > nodes. I know that while using delauny function for creating 2-D
> > triangular mesh one needs to reorder the vertices given as output
> > by
> > the function. Does the same thing has to be done here also and if
> > yes
> > does any body have any clue how to do that.
>
> Don't use delaunay to build that tessellation.
> This code will do it in n dimensions, without
> any of the irritating problems that use of
> delaunay will entail. (Yet another piece of
> code that belongs on the FEX.)
>
> I'm not sure what you mean about re-ordering
> vertices though. Are you looking for a surface
> triangulation of the tetrahedral mesh?
>
> HTH,
> John D'Errico
>
> function [vertices,tess]=tess_lat(varargin)
> % tess_lat: simplicial tessellation of a rectangular lattice
> % usage: [tessellation,vertices]=tess_lat(p1,p2,p3,...)
> %
> % arguments: input
> % p1,p2,p3,... - numeric vectors defining the lattice in
> % each dimension.
> % Each vector must be of length >= 1
> %
> % arguments: (output)
> % vertices - factorial lattice created from (p1,p2,p3,...)
> % Each row of this array is one node in the lattice
> % tess - integer array defining simplexes as references to
> % rows of "vertices".
>
> % dimension of the lattice
> n = length(varargin);
>
> % create a single n-d hypercube
> % list of vertices of the cube itself
> vhc=('1'==dec2bin(0:(2^n-1)));
> % permutations of the integers 1:n
> p=perms(1:n);
> nt=factorial(n);
> thc=zeros(nt,n+1);
> for i=1:nt
> thc(i,:)=find(all(diff(vhc(:,p(i,:)),[],2)>=0,2))';
> end
>
> % build the complete lattice
> nodecount = cellfun('length',varargin);
> if any(nodecount<2)
> error 'Each dimension must be of size 2 or more.'
> end
> vertices = lattice(varargin{:});
>
> % unrolled index into each hyper-rectangle in the lattice
> ind = cell(1,n);
> for i=1:n
> ind{i} = 0:(nodecount(i)-2);
> end
> ind = lattice(ind{:});
> k = cumprod([1,nodecount(1:(end-1))]);
> ind = 1+ind*k';
> nind = length(ind);
>
> offset=vhc*k';
> tess=zeros(nt*nind,n+1);
> L=(1:nind)';
> for i=1:nt
> tess(L,:)=repmat(ind,1,n+1)+repmat(offset(thc(i,:))',nind,1);
> L=L+nind;
> end
>
> % ======== subfunction ========
> function g = lattice(varargin)
> % generate a factorial lattice in n variables
> n=nargin;
> sizes = cellfun('length',varargin);
> c=cell(1,n);
> [c{1:n}]=ndgrid(varargin{:});
> g=zeros(prod(sizes),n);
> for i=1:n
> g(:,i)=c{i}(:);
> end



Hi ...sorry for the many years late reply, but I just came across this great piece of code, and wondered if you had:
a) made any file submissions based on this?
b) given any thought to/solved the problem that the use of ndgrid means only relatively small xyz vectors can be used (I tried using it with 20000 points to no avail)..

Thanks!!
James