latticegrid.m

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