Search code examples
3dwolfram-mathematicatetgen

How to use TetGen for this simple 3D geometry


Here are the points defining my simple 3D geometry.

datN = {{{-0.47150460764747554`, 0.29559274991660417`, 
 0.010131794240974218`}, {-0.4762714873728534`, 
 0.2955927499166042`, 
 0.010567957416020535`}, {-0.4835042628911566`, 
 0.29559274991660417`, 
 0.01066658601048008`}, {-0.49133736140975415`, 
 0.29559274991660417`, 
 0.01010572204377315`}, {-0.4974365622729896`, 
 0.29559274991660417`, 
 0.009602808597554033`}, {-0.4999590574180981`, 
 0.2955927499166041`, 
 0.010150149141898643`}, {-0.497870343592127`, 
 0.2955927499166042`, 
 0.011728012221066566`}, {-0.491634397829927`, 
 0.2955927499166041`, 
 0.013089897457762985`}, {-0.4834169387190052`, 
 0.2955927499166042`, 
 0.013009607974106477`}, {-0.47609963350102275`, 
 0.2955927499166043`, 
 0.011622413291940486`}, {-0.471504606936728`, 
 0.2955927499166041`, 
 0.010131794240974216`}}, {{-0.5619323339485054`, 
 0.13709660728856332`, 
 0.010131794240974218`}, {-0.5878076066290028`, 
 0.13709660728856335`, 
 0.01249934738636439`}, {-0.6270680976744502`, 
 0.13709660728856332`, 
 0.0130347168361427`}, {-0.6695872237650179`, 
 0.13709660728856332`, 
 0.00999027080199048`}, {-0.7026945171227986`, 
 0.13709660728856332`, 
 0.007260388089336815`}, {-0.7163869644835803`, 
 0.13709660728856332`, 
 0.010231427144215837`}, {-0.705049141229765`, 
 0.13709660728856338`, 
 0.018796282936276536`}, {-0.6711995779276564`, 
 0.13709660728856332`, 
 0.02618878157043711`}, {-0.6265940901692914`, 
 0.13709660728856332`, 
 0.02575295931296998`}, {-0.5868747603960375`, 
 0.13709660728856335`, 
 0.018223077560156144`}, {-0.5619323300904714`, 
 0.1370966072885633`, 0.010131794240974216`}}};

Now we prepare the facets and vertices

pt = Flatten[{datN[[1]], datN[[2]]}, 1];
facets = Join[{{Flatten@Map[Position[pt, #] &, datN[[1]]]}}, 
Table[{Flatten@
  Map[Position[pt, #] &, {datN[[1]][[i]], datN[[2]][[i]], 
    datN[[2]][[i + 1]], datN[[1]][[i + 1]]}]}, {i, 1, 
 10}], {{Flatten@Map[Position[pt, #] &, datN[[2]]]}}];

Then we use TetGen in the same line as it is described in the documentation for the simplest example with a box.

Needs["TetGenLink`"]
inInst = TetGenCreate[];
TetGenSetPoints[inInst, pt];
TetGenSetFacets[inInst, facets];
outInst = TetGenTetrahedralize[inInst, "pq1.414a0.01"];
coords = TetGenGetPoints[outInst];
surface = TetGenGetFaces[outInst];

We can see that no mesh is generated and also TetGenGetPoints failed to retun the vertices. Result very disappointing.

GraphicsGrid@{{Graphics3D[GraphicsComplex[pt, Map[Polygon, facets]], 
Boxed -> False], 
Graphics3D[GraphicsComplex[coords, Polygon[surface]]]}}

enter image description here

Why is this strange thing happening. TetGen documentation is not satisfactory either.


Solution

  • While in datN the begin and end points of the two sublists are effectively the same, they count as different points as far as Mathematica is concerned. This means that facets doesn't actually represent a closed polyhedron (there is a tiny gap between the edges {datN[[1,1]], datN[[2,1]]} and {datN[[1,-1]], datN[[2,-1]]}).

    To solve this, you could drop the end points from datN[[1]] and dat[[2]] and replace any occurrence of datN[[i]][[11]] with datN[[i]][[1]] in the definition of facets, e.g.

    datN2 = Drop[#, -1]& /@ datN;
    pt = Flatten[datN2, 1];
    facets = Join[
        {{Flatten@Map[Position[pt, #] &, datN2[[1]]]}}, 
        Table[{Flatten@
           Map[Position[pt, #] &, {datN2[[1]][[i]], datN2[[2]][[i]], 
             datN2[[2]][[Mod[i, 10] + 1]], 
             datN2[[1]][[Mod[i, 10] + 1]]}]}, {i, 1, 10}], 
        {{Flatten@Map[Position[pt, #] &, datN2[[2]]]}}];
    

    The rest of the code stays the same, i.e.

    Needs["TetGenLink`"]
    inInst = TetGenCreate[];
    TetGenSetPoints[inInst, pt];
    TetGenSetFacets[inInst, facets];
    outInst = TetGenTetrahedralize[inInst, "pq1.414a0.01"];
    coords = TetGenGetPoints[outInst];
    surface = TetGenGetFaces[outInst];
    

    Then plotting the surface gives the following result:

    GraphicsGrid@{{Graphics3D[GraphicsComplex[pt, Map[Polygon, facets]], 
      Boxed -> False], 
    Graphics3D[GraphicsComplex[coords, Polygon[surface]]]}}
    

    tetrahedralization