Analyze triangulation

This function returns information about a triangulation, which may be useful when setting up a partial differential equation problem.

Function analyzeTriangulation() was written to diagnose problems that can occur when using built-in functions alphaShape() and alphaTriangulation(). analyzeTriangulation() works on 2D and 3D (tetrahedral) triangulations.

function [mtri,N,vol,totV,aspRatio,sideL,minL,maxL,minV,maxV,minAR,maxAR,numZ,numNZ]...
    = analyzeTriangulation(tri,P)

Two scripts are provided that demonstrate the use of analyzeTriangulation() with triangles and tetrahedra respectyively: analyzeTriangulationTest2d.m and analyzeTriangulationTest3d.m.

The inputs are

%   tri=matrix of vertex numbers for each tri or tetra (mtri-by-nv) where
%       mtri=number of tris or tetras; nv=3 or 4 for 2D or 3D triangulations
%   P  =matrix of vertex cordinates (N-by-dim) where 
%       N=number of vertices; dim=2 for triangles, 3 for tetrahedra

Inputs tri,P may be obtained from [tri,P]=alphaTriangulation(...) or tri=delaunay(P).

The outputs are

%   mtri=number of triangles or tetrahedra
%   N=number of nodes
%   vol=area or volume of each triangle or tetra (mtri x 1)
%   totV=total area or volume of all tris or tetras
%   aspRatio=ratio of largest to smallest side length for each tri or tetra (mtri x 1)
%   sideL=lengths of sides of each tri or tetra (mtri x 3 or mtri x 6)
%   minL=shortest side length among all
%   maxL=longest side length among all
%   minV=area or volume of smallest tri or tetra
%   maxV=area or volume of largest tri or tetra
%   minAR=smallest aspect ratio of all tris or tetras
%   maxAR=largest aspect ratio
%   numZ = number of tris or tetras with area or volume=0 (<=1e-9)
%   numNZ= number of tris or tetras with area or volume>0 (>1e-9)

Here are three examples of problems that can occur with a triangulation, which can be identified and fixed with analyzeTriangulation():

1. triangles or tetrahedra with zero area or volume
2. incorrect filling in of internal voids in a shape
3. ill-conditioned (i.e. long and skinny) triangles or tetrahedra

Discussion of each example:

1. Triangles or tetrahedra with zero area or volume

alphaTriangulation() can produce a triangulation that include triangles or tetrahedra with area=0 or volume=0 (as seen in the example scripts). Such triangles can cause problems in a partial differential equation model. Use the "vol" output from analyzeTriangulation to quickly eliminate such triangles or tetrahedra from the triangulation. Suppose you have a triangulation defined by matrices tri, P (see above), and you wish to eliminate any triangles or tetrahedra with area or volume=0. Do this:

[mtri,N,vol] = analyzeTriangulation(tri,P);   % vol is a vector
tri2=tri(vol>1e-9,:)

The triangulation triA2, P has no triangles or tetrahedra with area or volume less than 1e-9.

2. Incorrect filling in of internal voids in a shape

When using alphaShape() and alphaTriangulation() on a shape that includes concave regions or internal voids, the resulting triangulation could fill in a region that should be a void. Or, with any shape, if the nodes are not close enough together, the triangulaiton could include a void where there should not be one. (This depends on the values of parameters alpha radius and HoleThreshold - see the help for alphaShape.) It can be difficult to tell if this kind of error has occurred, especially for 3D shapes with internal voids, because it is hard to see the inside of the triangulation. Therefore use the totV output to compare the total volume of the tetrahedra to the expected total volume. For example, if you have a cube with side length L that contains an internal spherical cavity of radius R, you can check that the total volume of tetrahedra has the expected value (which is ):

[mtri,N,vol,totV] = analyzeTriangulation(tri,P)
if abs(totV-(L^3-pi*4*R^3/3))>1e-6
    disp('Total volume of tetrahedra not equal to expected.')
end

You can fix the problem by adjusting the value of alpha radius or HoleThreshold passed to alphaShape(). The returned values maxL and maxV can be checked to make sure there are no unexpectedly large triangles - which could indicate that a concave region has been filled in when it should not have been.

3. Ill conditioned triangles or tetrahedra

Ill-conditioned (long and skinny) triangles or tetrahedra lead to numerical instability and unreliable results when solving PDEs, especially in finite element analysis. One can identify ill conditioned triangles or tetrahedra with the output variable aspRatio (aspect ratio): a vector containing the ratio of the longest to the shortest edge length for every triangle or tetrahedron. This is similar to condition number for a matrix, defined as the ratio of largest to smallest eigenvalue.

[mtri,N,vol,totV,aspRatio] = analyzeTriangulation(tri,P)
badTri=find(aspRatio>5);
if length(badTri)>0
    fprintf('The following rows of tri contain ill-conditioned triangles:\n')
    disp(badTri)
end

You can then determine the region(s) with ill-conditioned triangles, by checking P to find the coordinates of the nodes of the bad triangles, and you can fix the problem by adding more nodes in that region.