meshgrid6

 [node,elem]=meshgrid6(v1,v2,v3,...)

 mesh an ND rectangular lattice by splitting 
 each hypercube into 6 tetrahedra

 author: John D'Errico
 URL: http://www.mathworks.com/matlabcentral/newsreader/view_thread/107191
 modified by Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)

 input:
    v1,v2,v3,... - numeric vectors defining the lattice in
                   each dimension.
                   Each vector must be of length >= 1

 output:
    node - factorial lattice created from (v1,v2,v3,...)
           Each row of this array is one node in the lattice
    elem - integer array defining simplexes as references to
           rows of "node".

 example:
     [node,elem]=meshgrid6(0:5,0:6,0:4);
     plotmesh(node,elem);

 -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)

0001 function [node,elem]=meshgrid6(varargin)
0002 %
0003 % [node,elem]=meshgrid6(v1,v2,v3,...)
0004 %
0005 % mesh an ND rectangular lattice by splitting
0006 % each hypercube into 6 tetrahedra
0007 %
0008 % author: John D'Errico
0009 % URL: http://www.mathworks.com/matlabcentral/newsreader/view_thread/107191
0010 % modified by Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
0011 %
0012 % input:
0013 %    v1,v2,v3,... - numeric vectors defining the lattice in
0014 %                   each dimension.
0015 %                   Each vector must be of length >= 1
0016 %
0017 % output:
0018 %    node - factorial lattice created from (v1,v2,v3,...)
0019 %           Each row of this array is one node in the lattice
0020 %    elem - integer array defining simplexes as references to
0021 %           rows of "node".
0022 %
0023 % example:
0024 %     [node,elem]=meshgrid6(0:5,0:6,0:4);
0025 %     plotmesh(node,elem);
0026 %
0027 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0028 %
0029 
0030 % dimension of the lattice
0031 n = length(varargin);
0032 
0033 % create a single n-d hypercube
0034 % list of node of the cube itself
0035 vhc=('1'==dec2bin(0:(2^n-1)));
0036 % permutations of the integers 1:n
0037 p=perms(1:n);
0038 nt=factorial(n);
0039 thc=zeros(nt,n+1);
0040 for i=1:nt
0041   thc(i,:)=find(all(diff(vhc(:,p(i,:)),[],2)>=0,2))';
0042 end
0043 
0044 % build the complete lattice
0045 nodecount = cellfun('length',varargin);
0046 if any(nodecount<2)
0047   error 'Each dimension must be of size 2 or more.'
0048 end
0049 node = lattice(varargin{:});
0050 
0051 % unrolled index into each hyper-rectangle in the lattice
0052 ind = cell(1,n);
0053 for i=1:n
0054 ind{i} = 0:(nodecount(i)-2);
0055 end
0056 ind = lattice(ind{:});
0057 k = cumprod([1,nodecount(1:(end-1))]);
0058 ind = 1+ind*k';
0059 nind = length(ind);
0060 
0061 offset=vhc*k';
0062 elem=zeros(nt*nind,n+1);
0063 L=(1:nind)';
0064 for i=1:nt
0065   elem(L,:)=repmat(ind,1,n+1)+repmat(offset(thc(i,:))',nind,1);
0066 L=L+nind;
0067 end
0068 
0069 % ======== subfunction ========
0070 function g = lattice(varargin)
0071 % generate a factorial lattice in n variables
0072 n=nargin;
0073 sizes = cellfun('length',varargin);
0074 c=cell(1,n);
0075 [c{1:n}]=ndgrid(varargin{:});
0076 g=zeros(prod(sizes),n);
0077 for i=1:n
0078 g(:,i)=c{i}(:);
0079 end