Which algorithms exist to create a tetrahedral volume mesh from an STL file?

Question

Basically just that. I am trying to write a C++ program that reads in an STL and should compute an unstructured tetrahedral volume mesh based on the surface triangulation given by the STL file. I would like to implement the mesh generation step from scratch, so I am looking for algorithms or literature recommendations on algorithms, not existing library implementations or software like gmsh.

Answer 1

I don't agree with Wolfgang that you should never write things yourself, but writing a tetrahedral mesh generator is a big task and his estimate of a year's worth of work sounds about right. My favorite reference on this is the book by Cheng, Dey, and Shewchuk.

If you're not already familiar with it, most mesh generation algorithms rely on some variant of the Delaunay triangulation of the input points. In 2D, the (constrained) Delaunay triangulation of a point set maximizes the minimum angle among all possible triangulations. Proving that the approximate solutions from the finite element method converge to the true solution depends on the angles of the mesh cells being bounded above, so this optimality property is important. There are several algorithms to compute the Delaunay triangulation; the Bowyer-Watson algorithm is a common one but there are others based on local edge flips, others based on divide-and-conquer, etc.

Now, the constrained Delaunay triangulation is optimal for a given set of boundary and interior points, but you're starting with an STL file or some other description of just the boundary of the domain. The next challenging part of the procedure is to generate the interior points. Just because the Delaunay triangulation maximizes the minimum angle doesn't mean that the minimum angle you obtain is good enough. Ruppert's algorithm is an approach to generating the interior points such that the minimum angle stays above about 28.6${}^\circ$. Very often the minimum is much greater than that and improving on this lower bound is an active area of research.

All of what I've said applies just to 2D triangulations. In 3D, things become so, so much harder because of slivers -- tetrahedra that don't have large or small dihedral angles but that nonetheless have vanishingly small volume. In other words, the guarantees afforded by the Delaunay triangulation no longer get you what you want in 3D. Mesh generation algorithms in 3D then have to include a special step called sliver exudation.

Leaving aside the algorithms, I'll add that it is not at all obvious what a good data structure is to describe triangulations or more general unstructured meshes in 2D and 3D. Enforcing the mathematical invariants expected of a triangulation while also enabling the kinds of transformations that are needed for mesh generation is hard. I've tried to work on this before and it was more difficult than anything else I've ever done, by far. This is very different from, say, sparse matrices, where there's just about everyone would probably agree that the compressed sparse row / Harwell-Boeing format is a reasonable choice for most use cases.

When tackling big problems like this, I always try to think of what are some smaller problems that you can break it down into. For 2D Delaunay mesh generation, those steps might be:

1. unconstrained Delaunay triangulation of an input point set, including both the boundary and the interior points
2. constrained Delaunay triangulation, again of the boundary polygon and interior points
3. generating the interior mesh points through Delaunay refinement

Doing one or more of these is a reasonable project. You could then decide whether you want to roll your own 3D mesh generator or just use gmsh.

Answer 2

Creating meshes is a non-trivial enterprise and you will probably spend a year or two implementing the relevant algorithms yourself. Instead, just use what others have used -- specifically, use the gmsh program that already does everything you want.

As for the specific question: STL surfaces are really just one specific way of describing a surface. The algorithms to create meshes will simply have to query properties of the boundary, which can be implemented in a generic way so that STL is one way and other descriptions of the boundary are other ways to implement the generic interface. The mesh generation algorithm itself shouldn't care which one you use.