meshgrid5

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

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

 author: Qianqian Fang, <q.fang at neu.edu>
 inspired by John D'Errico
 URL: http://www.mathworks.com/matlabcentral/newsreader/view_thread/107191

 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]=meshgrid5(0:5,0:6,0:4);
     plotmesh(node,elem);

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

0001 function [node,elem]=meshgrid5(varargin)
0002 %
0003 % [node,elem]=meshgrid5(v1,v2,v3,...)
0004 %
0005 % mesh an ND rectangular lattice by splitting
0006 % each hypercube into 5 tetrahedra
0007 %
0008 % author: Qianqian Fang, <q.fang at neu.edu>
0009 % inspired by John D'Errico
0010 % URL: http://www.mathworks.com/matlabcentral/newsreader/view_thread/107191
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]=meshgrid5(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 if(n~=3) error('only works for 3D case!'); end
0033 
0034 for i=1:n
0035     v=varargin{i};
0036     if(mod(length(v),2)==0)
0037         varargin{i}=linspace(v(1),v(end),length(v)+1);
0038     end
0039 end
0040 
0041 % create a single n-d hypercube
0042 % list of node of the cube itself
0043 
0044 cube8=...
0045 [1 4 5 13;1 2 5 11;1 10 11 13;11 13 14 5;11 13 1 5;...
0046  2 3 5 11;3 5 6 15;15 11 12 3;15 11 14 5;11 15 3 5;...
0047  4 5 7 13;5 7 8 17;16 17 13 7;13 17 14 5;5 7 17 13;...
0048  5 6 9 15;5 8 9 17;17 18 15 9;17 15 14 5;17 15 5 9;...
0049  10 13 11 19;13 11 14 23;22 19 23 13;19 23 20 11;13 11 19 23;...
0050  11 12 15 21;11 15 14 23;23 21 20 11;23 24 21 15;23 21 11 15;...
0051  16 13 17 25;13 17 14 23;25 26 23 17;25 22 23 13;13 17 25 23;...
0052  17 18 15 27;17 15 14 23;26 27 23 17;27 23 24 15;23 27 17 15]';
0053 
0054 % build the complete lattice
0055 nodecount = cellfun('length',varargin);
0056 if any(nodecount<2)
0057     error 'Each dimension must be of size 2 or more.'
0058 end
0059 node = lattice(varargin{:});
0060 
0061 [ix,iy,iz]=meshgrid(1:2:nodecount(1)-2,1:2:nodecount(2)-2,1:2:nodecount(3)-2);
0062 ind=sub2ind(nodecount,ix(:),iy(:),iz(:));
0063 
0064 nodeshift=[0 1 2 nodecount(1) nodecount(1)+1 nodecount(1)+2 ...
0065 2*nodecount(1) 2*nodecount(1)+1 2*nodecount(1)+2];
0066 nodeshift=[nodeshift,nodeshift+nodecount(1)*nodecount(2),nodeshift+2*nodecount(1)*nodecount(2)];
0067 
0068 nc=length(ind);
0069 elem=zeros(nc*40,4);
0070 for i=1:nc
0071     elem((1:40)+(i-1)*40,:)=reshape(nodeshift(cube8(:)),4,40)'+ind(i);
0072 end
0073 
0074 % ======== subfunction ========
0075 function g = lattice(varargin)
0076 % generate a factorial lattice in n variables
0077 n=nargin;
0078 sizes = cellfun('length',varargin);
0079 c=cell(1,n);
0080 [c{1:n}]=ndgrid(varargin{:});
0081 g=zeros(prod(sizes),n);
0082 for i=1:n
0083     g(:,i)=c{i}(:);
0084 end