This post is about generating meshes for finite element simulations. I'll be covering other aspects of FEM based simulation in a later post, until then I recommend checking out Matthias Müller's very good introduction in the SIGGRAPH 2008 Real Time Physics course .
After spending the last few weeks reading, implementing and debugging meshing algorithms I have a new-found respect for people in this field. It is amazing how many ways meshes can "go wrong", even the experts have it tough:
“I hate meshes. I cannot believe how hard this is. Geometry is hard.” — David Baraff, Senior Research Scientist, Pixar Animation Studios
Meshing algorithms are hard, but unless you are satisfied simulating cantilever beams and simple geometric shapes you will eventually need to deal with them.
My goal was to find an algorithm that would take an image as input, and produce as output a good quality triangle mesh that conformed to the boundary of any non-zero regions in the image.
My first attempt was to perform a coarse grained edge detect and generate a Delaunay triangulation of the resulting point set. The input image and the result of a low-res edge detect:
This point set can be converted to a mesh by any Delaunay triangulation method, the Bowyer-Watson algorithm is probably the simplest. It works by inserting one point at a time, removing any triangles whose circumcircle is encroached by the new point and re-tessellating the surrounding edges. A nice feature is that the algorithm has a direct analogue for tetrahedral meshes, triangles become tetrahedra, edges become faces and circumcircles become circumspheres.
Here's an illustration of how Bowyer/Watson proceeds to insert the point in red into the mesh:
And here is the Delaunay triangulation of the Armadillo point set:
As you can see, Delaunay triangulation algorithms generate the convex hull of the input points. But we want a mesh that conforms to the shape boundary - one way to fix this is to sample the image at each triangle's centroid, if the sample lies outside the shape then simply throw away the triangle. This produces:
Much better! Now we have a reasonably good approximation of the input shape. Unfortunately, FEM simulations don't work well with long thin "sliver" triangles. This is due to interpolation error and because a small movement in one of the triangle's vertices leads to large forces, which leads to inaccuracy and small time steps .
Before we look at ways to improve triangle quality it's worth talking about how to measure it. One measure that works well in 2D is the ratio of the triangle's circumradius to it's shortest edge. A smaller ratio indicates a higher quality triangle, which intuitively seems reasonable, long skinny triangles have a large circumradius but one very short edge:
The triangle on the left, which is equilateral, has a ratio ~0.5 and is the best possible triangle by this measure. The triangle on the right has a ratio of ~8.7, note the circumcenter of sliver triangles tend to fall outside of the triangle itself.
Methods such as Chew's algorithm and Ruppert's algorithm are probably the most well known refinement algorithms. They attempt to improve mesh quality while maintaining the Delaunay property (no vertex encroaching a triangle's circumcircle). This is typically done by inserting the circumcenter of low-quality triangles and subdividing edges.
Unfortunately these algorithms require an accurate polygonal boundary as input as the output is sensitive to the input segment lengths. They are also famously difficult to implement robustly and efficiently, I spent most of my time implementing Ruppert's algorithm only to find the next methods produced better results with much simpler code.
Variational (energy based) algorithms improve the mesh through a series of optimization steps that attempt to minimize a global energy function. I adapted the approach in Variational Tetrahedral Meshing  to 2D and found it produced great results, this is the method I settled on so I'll go into some detail.
The algorithm proceeds as follows:
- Generate a set of uniformly distributed points interior to the shape P
- Generate a set of points on the boundary of the shape B
- Generate a Delaunay triangulation of P
- Optimize boundary points by moving them them to the average of their neighbours in B
- Optimize interior points by moving them to the centroid of their Voronoi cell (area weighted average of connected triangle circumcenters)
- Unless stopping criteria met, go to 3.
- Remove boundary sliver triangles
The core idea is that of repeated triangulation (3) and relaxation (4,5), it's a somewhat similar process to Lloyd's clustering, conincidentally the same algorithm I had used to generate surfel hierarchies for global illumination sims in the past.
Here's an animation of 7 iterations on the Armadillo, note the number of points stays the same throughout (another nice property):
It's interesting to see how much the quality improves after the very first step. Although Alliez et al.  don't provide any guarantees on the resulting mesh quality I found the algorithm works very well on a variety of input images with a fixed number of iterations.
This is the algorithm I ended up using but I'll quickly cover a couple of alternatives for completeness.
These algorithms typically start by tiling interior space using a BCC (body centered cubic) lattice which is simply two interleaved grids. They then generate a Delaunay triangulation and throw away elements lying completely outside the region of interest.
As usual, handling boundaries is where the real challenge lies, Molino et al.  use a force based simulation to push grid points towards the boundary. Isosurface Stuffing  refines the boundary by directly moving vertices to the zero-contour of a signed distance field or inserts new vertices if moving the existing lattice nodes would generate a poor quality triangle.
Lattice based methods are typically very fast and don't suffer from the numerical robustness issues of algorithms that rely on triangulation. However if you plan on fracturing the mesh along element boundaries then this regular nature is exposed and looks quite unconvincing.
Another approach is to start with a very fine-grained mesh and progressively simplify it in the style of Progressive Meshes . Barbara Cutler's thesis and associated paper discusses the details and very helpfully provides the resulting tetrahedral meshes, but the implementation appears to be considerably more complex than variational methods and relies on quite a few heuristics to get good results.
Now the mesh is ready it's time for the fun part (apologies if you really love meshing). This simple simulation is using co-rotational linear FEM with a semi-implicit time-stepping scheme:
(Armadillo and Bunny images courtesy of the Stanford Scanning Respository)
Pre-built binaries for OSX/Win32 here: http://mmacklin.com/fem.zip
Source code is available on Github: https://github.com/mmacklin/sandbox/tree/master/projects/fem.
 Matthias Müller, Jos Stam, Doug James, and Nils Thürey. Real time physics: class notes. In ACM SIGGRAPH 2008 classes http://www.matthiasmueller.info/realtimephysics/index.html
 Jonathan Richard Shewchuk. 2002. What Is a Good Linear Finite Element? Interpolation, Conditioning, Anisotropy, and Quality Measures, unpublished preprint. http://www.cs.berkeley.edu/~jrs/papers/elemj.pdf
 Pierre Alliez, David Cohen-Steiner, Mariette Yvinec, and Mathieu Desbrun. 2005. Variational tetrahedral meshing. ftp://ftp‑sop.inria.fr/prisme/alliez/vtm.pdf
 Molino, Bridson, et al. - 2003. A Crystalline, Red Green Strategy for Meshing Highly Deformable Objects with Tetrahedra http://www.math.ucla.edu/~jteran/papers/MBTF03.pdf
 François Labelle and Jonathan Richard Shewchuk. 2007. Isosurface stuffing: fast tetrahedral meshes with good dihedral angles. In ACM SIGGRAPH 2007 papers http://www.cs.berkeley.edu/~jrs/papers/stuffing.pdf
 Hugues Hoppe. 1996. Progressive meshes. http://research.microsoft.com/en-us/um/people/hoppe/pm.pdf