meshgrid6.m

function [node, elem] = meshgrid6(varargin)
%
% [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 (q.fang at neu.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)
%

% dimension of the lattice
n = length(varargin);

% create a single n-d hypercube
% list of node of the cube itself
vhc = ('1' == dec2bin(0:(2^n - 1)));
% permutations of the integers 1:n
p = perms(1:n);
nt = factorial(n);
thc = zeros(nt, n + 1);
for i = 1:nt
    thc(i, :) = find(all(diff(vhc(:, p(i, :)), [], 2) >= 0, 2))';
end

% build the complete lattice
nodecount = cellfun('length', varargin);
if any(nodecount < 2)
    error 'Each dimension must be of size 2 or more.';
end
node = lattice(varargin{:});

% unrolled index into each hyper-rectangle in the lattice
ind = cell(1, n);
for i = 1:n
    ind{i} = 0:(nodecount(i) - 2);
end
ind = lattice(ind{:});
k = cumprod([1, nodecount(1:(end - 1))]);
ind = 1 + ind * k';
nind = length(ind);

offset = vhc * k';
elem = zeros(nt * nind, n + 1);
L = (1:nind)';
for i = 1:nt
    elem(L, :) = repmat(ind, 1, n + 1) + repmat(offset(thc(i, :))', nind, 1);
    L = L + nind;
end

if (n == 2 || n == 3)
    elem = meshreorient(node, elem);
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