Broadly speaking, there are two kinds of ways to describe geometric regions: explicitly and implicitly. Meshes are explicit representations consisting of lists of vertex coordinates and connectivity information (e.g. edges, triangles, ... ). Implicit representations define regions through a function, where the sign of the output classifies points as being inside, outside or on the boundary. Some tasks are easy with explicit representations but hard for implicit, like enumerating points on the boundary or calculating the center of mass. Other tasks are trivial for implicit representations, but tricky with explicit, like determining if a point belongs to the region or not.

Because of the inherent benefits and drawbacks of each representation, it is frequently useful to be able to convert between them. This document goes over an algorithm for creating a simplex mesh (an explicit representation made up of triangles / tetrahedra) of an implicitly-defined region, as originally described in this paper. Along the way, we'll note some important modifications to that algorithm that make it more effective and performant.

Defining Regions from Level Set Functions

The main idea behind this algorithm is that if you have a function, $f(x,y)$$f(x, y)$, you can define a region from it. We'll use the following arbitrarily-chosen function pictured below for the examples in this document:

If we also draw a plane at $z=0$$z=0$, then the 2D region defined by $f(x,y)\le 0$$f(x, y) \leq 0$ is the surface of the "lake" in the picture

If we just look at it from above, that region looks like

Great, so we can generate a squiggle-shaped region from this function $f(x,y)$$f(x,y)$, but how do we create a mesh of that region?

Sampling and Background Mesh

Our algorithm starts by creating a background mesh of simplex elements around the region and evaluating $f$$f$ at each vertex. From the sign of the vertex values, we classify points as inside, outside and on the boundary (blue, red and purple, respectively).

Clearly, the triangles made up of only boundary and exterior points don't belong in our mesh, so we discard them

Now our mesh has started to capture the overall shape of the region, but it still needs work. If we only cared about discretizing the boundary of the region, we could apply a generalization of the marching cubes algorithm on this data. However, I'm interested in generating meshes suitable for finite element analysis, and that means we need quality elements (i.e. as close to equilateral triangles as possible) of the whole region.

Boundary Fitting

The background mesh has good quality elements, but it doesn't conform to the boundary of the region. The original algorithm addresses this by moving all of the exterior vertices inward until they lie on the boundary.

Now the mesh matches the actual shape of the region closely, but the element shapes near the boundary are really distorted (or even inverted). That means this mesh is unsuitable for finite element analysis in its current state. Later on, we'll see how the original algorithm in the paper mitigates this problem, but I recommend a different approach for this step.

The really misshapen elements seem to show up in places where there was an interior vertex close to the boundary. Instead, if we start by moving the "barely-inside" vertices outward onto the boundary, then we won't have to move vertices as far, so the element qualities stay high. This two pass approach is shown below

In addition to producing higher quality elements, this approach also lets us discard additional triangles after the first pass (where we moved the "barely inside" vertices onto the boundary), so the resulting mesh needs fewer vertices and elements.

Vertex Relaxation

The last step of the algorithm is an iterative procedure to improve element qualities. At a high level, this involves moving the mesh vertices so that the minimum element quality goes up. The procedure described in the paper for improving the vertex positions is more complicated and expensive than it needs to be.

Instead, the same effect can be accomplished in a fraction of the time and effort with simple gradient ascent. If we define the "quality" of each vertex position to be the minimum of the qualities of elements it belongs to, then we can just push the vertex along the gradient of that function to improve the mesh quality

In practice, even severely misshapen or inverted elements are fixed in just a few iterations. More details about this step, including all the relevant mathematical expressions and C++ code, are available in a related document here.

Summary

This algorithm is remarkably simple, and tends to work well in practice, but there are a few important limitations:

1. It only works with simplex elements. Many finite element codes prefer quadrilateral or hexahedral meshes, but the boundary-fitting step doesn't work robustly with those kinds of elements.

2. It only works for smooth regions.

3. Small geometric features require refining the background mesh to resolve.

For 3D, the algorithm plays out almost exactly the same way, except that it's a little less clear what the background mesh should be since regular tetrahedra do not tile space. In my implementation, I use the A15 lattice for the background grid, since it is easy to implement and the tetrahedrons in that arrangement are high quality.