function [vertices,tess]=genT5mesh(varargin)
% tess_lat: simplicial tessellation of a rectangular lattice
% usage: [tessellation,vertices]=genT5mesh(v1,v2,v3,...)
%
% arguments: input
% v1,v2,v3,... - numeric vectors defining the lattice in
% each dimension.
% Each vector must be of length >= 1
%
% arguments: (output)
% vertices - factorial lattice created from (v1,v2,v3,...)
% 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);
if(n~=3) error('only works for 3D case!'); end

for i=1:n
	v=varargin{i};
	if(mod(length(v),2)==0)
		varargin{i}=linspace(v(1),v(end),length(v)+1);
	end
end

% create a single n-d hypercube
% list of vertices of the cube itself

cube8=...
[1 4 5 13;1 2 5 11;1 10 11 13;11 13 14 5;11 13 1 5;...
 2 3 5 11;3 5 6 15;15 11 12 3;15 11 14 5;11 15 3 5;...
 4 5 7 13;5 7 8 17;16 17 13 7;13 17 14 5;5 7 17 13;...
 5 6 9 15;5 8 9 17;17 18 15 9;17 15 14 5;17 15 5 9;...
 10 13 11 19;13 11 14 23;22 19 23 13;19 23 20 11;13 11 19 23;...
 11 12 15 21;11 15 14 23;23 21 20 11;23 24 21 15;23 21 11 15;...
 16 13 17 25;13 17 14 23;25 26 23 17;25 22 23 13;13 17 25 23;...
 17 18 15 27;17 15 14 23;26 27 23 17;27 23 24 15;23 27 17 15]';

% 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{:});

[ix,iy,iz]=meshgrid(1:2:nodecount(1)-2,1:2:nodecount(2)-2,1:2:nodecount(3)-2);
ind=sub2ind(nodecount,ix(:),iy(:),iz(:));

nodeshift=[0 1 2 nodecount(1) nodecount(1)+1 nodecount(1)+2 ...
2*nodecount(1) 2*nodecount(1)+1 2*nodecount(1)+2];
nodeshift=[nodeshift,nodeshift+nodecount(1)*nodecount(2),nodeshift+2*nodecount(1)*nodecount(2)];

nc=length(ind);
tess=zeros(nc*40,4);
for i=1:nc
	tess((1:40)+(i-1)*40,:)=reshape(nodeshift(cube8(:)),4,40)'+ind(i);
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
