function [node,face,centroids]=latticegrid(varargin)
%
% [node,face,centroids]=latticegrid(xrange,yrange,zrange,...)
%
% generate a 3D lattice
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
% xrange, yrange, zrange ...: 1D vectors specifying the range of each
% dimension of the lattice
%
% output:
% node: the vertices of the 3D lattice
% face: the list of cell faces of the lattice, including both internal
% and external facets. By default, face is in the form of a cell
% array, with each row representing a face. One can use
% cell2mat(face) to convert it to an array
% centroids: the centroids of each lattice cell
%
% example:
% % generate a 3D lattice
% [node,face,c0]=latticegrid([1 2 4],1:3,1:4);
% plotmesh(node,face)
%
% % mesh the 3D lattice based on the face info
% [no,el]=surf2mesh(node,face,[],[],1,0.01,c0);
% figure; plotmesh(no,el)
%
% % mesh a 2-layer structure using a simple lattice
% [node,face,c0]=latticegrid([0 10],[0 5],[0 3.5 4]);
% c0(:,4)=[0.01;0.001];
% [no,el]=surf2mesh(node,face,[],[],1,[],c0);
% figure; plotmesh(no,el)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%
n=length(varargin);
p=cell(n,1);
[p{:}]=ndgrid(varargin{:});
node=zeros(length(p{1}(:)),n);
for i=1:n
node(:,i)=p{i}(:);
end
if(nargout==1)
return;
end
dim=size(p{1});
dd=[dim(1) dim(1)*dim(2)];
onecube=[0 dd(1) dd(1)+1 1; ...
0 1 dd(2)+1 dd(2); ...
0 dd(2) dd(2)+dd(1) dd(1)];
onecube=[onecube;onecube+repmat([dd(2);dd(1);1],1,4)];
len=prod(dim(1:3)-1);
face=repmat(onecube,len,1);
[xx,yy,zz]=ndgrid(1:dim(1)-1,1:dim(2)-1,1:dim(3)-1);
idx=sub2ind(dim,xx(:),yy(:),zz(:))';
orig=repmat(idx,size(onecube,1),1);
for i=1:size(onecube,2)
face(:,i)=face(:,i)+orig(:);
end
face=unique(face,'rows');
face=mat2cell(face,ones(size(face,1),1));
if(nargout>=3)
diffvec=cellfun(@diff,varargin,'UniformOutput',false);
[xx,yy,zz]=ndgrid(diffvec{:});
centroids=node(idx,:)+[xx(:) yy(:) zz(:)]*0.5;
end