Mesh Generation Using MATLAB Delaunay Function

If you are using Matlab Delaunay function to generate a unstructured 2D or 3D mesh you will land up in some sort of trouble. For very simple reasons. Have a look at the following piece of MATLAB code.

This code is to generate 3D tetrahedral mesh using the Matlab function Delaunay

x=[(0):1/(numx):(1)]; % numx=numy=numz = 1 (typical cube)

y=[(0):1/(numy):(1)];

z=[(0):1/(numz):(1)];
[X,Y,Z] = meshgrid(x,y,z);

X = [X(:)];

Y= [Y(:)];

Z = [Z(:)];

Tes = delaunay3(X,Y,Z);

L = [X(:) Y(:) Z(:)];

tetramesh1(Tes,L);camorbit(20,0) ;

This gives me a tetrahedral mesh with 'Tes' having the nodal connectivity matrix for a typical case (a cube with 6-tetrahedra, 8 nodes)it looks like

Tes = 4 7 8 6

1 5 7 6

2 4 7 6

2 1 7 6

2 4 7 3

2 1 7 3

Now if you are not care full you might land up in trouble. The reason being that output of delaunay function is in a random order some of the tessellation are arranged in clock wise direction and some in anticlockwise direction. So you need to reorder the tessellation in a particluar oder keeping one/summit node fixed and arranging others in clockwise or anticlcokwise direction.

Same thing happens in 2D. Now you can either do that if you are fond of MATLAB function Delaunay or you can use fowllowing piece of code given to me by some one who is master in tetrahedral meshes.

function [vertices,tess]=tess_lat(varargin)
% tess_lat: simplicial tessellation of a rectangular lattice
% usage: [tessellation,vertices]=tess_lat(p1,p2,p3,...)
%
% arguments: input
% p1,p2,p3,... - numeric vectors defining the lattice in
% each dimension.
% Each vector must be of length >= 1
%
% arguments: (output)
% vertices - factorial lattice created from (p1,p2,p3,...)
% 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);


% create a single n-d hypercube
% list of vertices 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
vertices = 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';
tess=zeros(nt*nind,n+1);
L=(1:nind)';
for i=1:nt
tess(L,:)=repmat(ind,1,n+1)+repmat(offset(thc(i,:))',nind,1);
L=L+nind;
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

In the next post I will discuss how you can get boundary nodes in a much easier fashion for any kind of meshes by using Faces in 3D.