Jonathan's research

The importance of mesh generation. One of the central tools of scientific computing is the forty-year old finite element method--a numerical method for approximating solutions to partial differential equations. The finite element method is used to develop computer simulations of a wide variety of physical phenomena, including fluid flow, heat transfer, mechanical deformation, and electromagnetic wave propagation. It is used heavily in industry and science for marvelously diverse purposes--for instance, determining optimal pumping strategies for petroleum extraction, anticipating the electromagnetic behavior of integrated circuits, and studying phenomena from quantum mechanics to earthquakes to black holes. The major automakers invest heavily in finite element software and expertise to help design engines and study the aerodynamics of car bodies. The aerodynamic response of the Boeing 777 airplane was simulated by the finite element method before even the first reduced-scale physical model was constructed for wind tunnel tests.

Joe F. Thompson, a celebrated engineer who commanded a multi-institutional mesh generation project known as the National Grid Project wrote in 1992 that

An essential element of the numerical solution of partial differential equations (PDEs) on general regions is the construction of a grid (mesh) on which to represent the equations in finite form. The grid must be generated for the regions of interest, and this is far from being a trivial problem. In fact, at present it can take orders of magnitude more man-hours to construct the grid than it does to perform and analyze the PDE solution on the grid. This is especially true now that PDE codes of wide applicability are becoming available, and grid generation has been cited repeatedly as being a major pacing item. The PDE codes now available typically require much less esoteric expertise of the knowledgeable user than do the grid generation codes.

Twelve years later, mesh generation is still usually the bottleneck in solving PDEs on objects or domains with complicated geometry. The automatic mesh generation problem is to program a computer to divide a complicated geometry--say, a car engine, the Golden Gate bridge, or the air around an airplane--into small, simple pieces called elements, such as triangles or rectangles (for two-dimensional geometries), or tetrahedra or rectangular prisms (for three-dimensional geometries). The division may yield millions or billions of such elements, in which case it is infeasible for a human to interactively help the computer with its decisions; the solution must be found entirely automatically.

To devise such a computer program is a difficult mathematical challenge. Mesh generation is something like a reverse jigsaw puzzle. A mesh must satisfy nearly contradictory requirements: it must conform to the shape of the object or domain; it may have to grade from small to large elements over a relatively short distance; and it must be composed of elements that are of the right sizes and shapes. ``The right shapes'' typically include elements that are nearly equilateral, and exclude elements that are long and thin (e.g. shaped like a needle or a kite). Some versions of the mesh generation problem--for example, dividing an arbitrary three-dimensional domain into ``nicely shaped'' hexahedra (elements shaped like rectangular prisms, though the sides need not be parallel)--are widely believed to have no general solution.

Since Thompson wrote his article, the resources committed to mesh generation have expanded rapidly. The field has seen the birth of three dedicated conferences: the annual International Meshing Roundtable, the biennial Symposium on Trends in Unstructured Mesh Generation, and the biennial International Conference on Numerical Grid Generation in Computational Field Simulations. Dozens of companies market mesh generation software. Many of the national research laboratories and most of the large automakers have in-house efforts.

Mesh generation is truly an interdisciplinary topic. Most of the early work was done by researchers from several branches of engineering, especially mechanics and fluid dynamics. Unfortunately, mesh generation is such a difficult problem that all the algorithms developed during this period are fragile, and often fail to function when confronted by small domain features, faces meeting at small angles, or complicated topologies that were not anticipated by the programmer.

Around the time of Thompson's article, these problems attracted the interest of researchers in computational geometry. Whereas engineers were once satisfied with mesh generators that (usually) work for their particular domains, those of us in computational geometry set a loftier (and more difficult) goal: provably good mesh generation, the design of algorithms that are mathematically guaranteed to produce a satisfying mesh of any geometry--even one unimagined by the algorithm designer.

Meanwhile, the applications of mesh generation software have expanded far beyond their provenance. As computers have grown more powerful, triangulated surface models have become widely used in computer graphics, and graphics conferences have seen an explosion of work in mesh generation and mesh processing. In economic terms, the computer graphics and animation industry probably now exceeds the finite element industry as a user of meshes.

More generally, meshes are used for hundreds of applications of multivariate interpolation. For example, if a cartographer knows the altitudes of the ground at a finite number of locations, she can generate a mesh of triangles whose corners lie at those locations, estimate the altitude at any other point (where an altitude measurement was not taken) by interpolation, and generate a contour map of the terrain. Meshes are used heavily in geographic information systems (GIS), image processing, radio propagation analysis, and countless other applications.

My work in meshing. I developed one of the earliest provably good three-dimensional mesh generation algorithms (the first to be as practical as the non-guaranteed algorithms); implemented a popular two-dimensional mesh generator; and expanded the range of geometries that mesh generation algorithms reliably work on (thereby converting several algorithms from hit-and-miss to truly practical). In collaboration with my student François Labelle, we have just developed the first provably good algorithm for anisotropic mesh generation. My work extends from fundamental combinatorial geometry to implementation.

Two-dimensional triangular mesh generation. The 1990s brought about two-dimensional Delaunay refinement algorithms that not only work well in practice, but have provable bounds on element quality and mesh grading. My mesh generator Triangle is an industrial-strength implementation of two-dimensional Delaunay refinement, based on a combination of algorithms by Jim Ruppert, Paul Chew, and myself. Triangle has been freely available to the public since 1995, and has undergone several revisions, most recently in 2002.

Part of my recent work has been to remove the last barrier to the practicality of these algorithms--their tendency to fail in the presence of small input angles. Unfortunately, the algorithms of Ruppert and Chew rely on an unlikely assumption: that no angle in the domain is less than 60^o. In the real world, objects have sharp corners. The fault is not just theoretical; in practice, these algorithms fail to terminate on many input domains that have small interior angles. In fact, there is a negative result that shows that some domains will foil any algorithm that isn't willing to compromise on mesh quality--any triangulation of such a domain must have one or more elements with a new small angle (besides the small angles in the domain geometry itself, which cannot be removed.) This leaves the algorithm designer with a difficult dilemma: a mesh generator that too aggressively eradicates bad elements may fail to produce a mesh at all.

I have shown how to amend several mesh generation algorithms (including Ruppert's algorithm, Chew's algorithm, and my three-dimensional algorithms discussed later) to use a ``look-ahead'' step to decide when the algorithms must locally compromise the quality of the mesh. The improved algorithms cope well with small angles, they never fail to produce a mesh, and the quality of the mesh degrades gracefully with the intractability of the input domain--and all of these qualities are mathematically guaranteed. These improvements are incorporated into the latest version of Triangle, and work even better in practice than in theory. In conjunction with my work on numerical robustness (discussed below), they make Triangle the most robust mesh generation software available. The figures at right show an example of a difficult domain (whose nine small angles confound most mesh generators) and a mesh thereof produced by Triangle, wherein most angles are greater than 20^o. There are fourteen exceptions (angles less than 20^o), and these are essentially unavoidable.

Triangle is downloaded from the Netlib repository of mathematical software an average of 31 times a day, and has been downloaded over 75,000 times since it was first placed on Netlib in 1995. The software has thousands of users (including many hundreds I have heard from personally), and has been licensed for inclusion in thirteen commercial software products. Commercial licensees use Triangle for tasks like digital ink and paint and compositing for cartoon animation, cloth animation (in Alias|Wavefront's well-known Maya computer animation package), surveying and geographical data processing, visualizing three-dimensional ocean floor models, mapping and contouring of subsurface geology and geophysics, and finite element simulation of groundwater flow, contamination transport, and thermal modeling. Non-commercial users use Triangle (free of charge) for discontinuity meshing for global illumination; terrain databases for real-time simulations; stereographic vision; interpolation of speech signals; and finite element methods for electromagnetic modeling, integrated circuit design, propagation of electric currents through myocardial tissue, transport processes in estuaries and aquifers, and earthquake rupture dynamics; plus hundreds of other applications. In 2003, Triangle received the James Hardy Wilkinson Prize in Numerical Software.

Three-dimensional tetrahedral mesh generation. I have developed provably good three-dimensional mesh generation algorithms based on Delaunay refinement. The jump from two dimensions to three is surprisingly difficult. Nonetheless, one of my algorithms provably generates a nicely graded mesh whose tetrahedral elements have circumradius-to-shortest edge ratios bounded below 1.63. Unlike previous provably good tetrahedral mesh generation algorithms, my algorithm can mesh general straight-line nonmanifold domains with holes, interior boundaries, and dangling edges and facets--and it works well in practice! (Previous provably good algorithms either generated far too many tetrahedra, or could only mesh convex objects.)

My theoretical results ensure that most types of bad tetrahedra cannot appear in the mesh, but do not rule out the possibility of a bad tetrahedron known as a sliver. Fortunately, Delaunay refinement algorithms outperform their worst-case bounds in practice. Pyramid, my implementation of three-dimensional Delaunay refinement, can usually generate meshes of tetrahedra whose dihedral angles are bounded between 20^o and 150^o (including the tire incinerator mesh above).

When my three-dimensional Delaunay refinement algorithm is combined with my extensions for handling small input angles (which rely upon constrained Delaunay tetrahedralizations, discussed below), the algorithm becomes truly practical, guaranteed to handle any piecewise linear domain.

Anisotropic mesh generation. This work combines a formerly unsolved problem of real importance to mesh generation; a novel theoretical approach; and an algorithm that is practical and the first to be provably good.

Most mesh generation algorithms try to create triangles or tetrahedra that are close to equilateral. But for some applications of interpolation and numerical modeling, the mesh must be anisotropic: having long, skinny triangles with orientations and aspect ratios dictated by the function being approximated. The importance of anisotropic meshes is stressed more loudly every year, especially by engineers using finite elements for applied problems. Functions with strongly anisotropic curvature are best interpolated with anisotropic meshes, and partial differential equations that are inherently anisotropic yield the best-conditioned numerical methods if anisotropic elements are used. For example, laminar air flow over an airplane wing is best modeled with extremely thin slab-shaped elements aligned with the surface of the wing. The ideal orientations and aspect ratios of the elements vary (smoothly) from one point in space to another, making anisotropic meshing more difficult than the traditional mesh generation problem.

Delaunay triangulations have mathematical properties that make them excellent for mesh generation, especially provably good mesh generation. However, it appears to be difficult or impossible to adapt Delaunay triangulations to anisotropic metric spaces without losing some of the properties necessary to create an algorithm that is both practical and provably good. Hence, François Labelle and I have developed a novel approach: we eschew the usual idea of constructing and improving a triangulation incrementally. Instead, we have invented a geometric construction that we call an anisotropic Voronoi diagram (a generalization of the well-known multiplicatively weighted Voronoi diagram). See above for an example. Whereas the geometric dual of an ordinary Voronoi diagram is a Delaunay triangulation, the geometric dual of an anisotropic Voronoi diagram is not always a triangulation at all; anisotropic Voronoi diagrams can be quite badly behaved. We can tame an anisotropic Voronoi diagram by inserting new sites at carefully chosen positions, until the geometric dual of the diagram is a guaranteed-quality anisotropic mesh, as illustrated above.

Our algorithm produces a mesh in which no triangle has an angle smaller than 20^o, as measured from the skewed perspective (i.e. metric tensor) of any point in the triangle. This is the first provable guarantee of any kind in anisotropic meshing. François has implemented the algorithm and verified that it produces excellent meshes in practice. The figures above show an anisotropic mesh and the anisotropic Voronoi diagram used to generate it.

• My papers on mesh generation by Delaunay refinement.
• The Triangle page.

My predicates (or variants thereof) are incorporated in several software programs, including the Cart3D inviscid analysis package by Michael Aftosmis, Marsha Berger, and John Melton; the Power Crust surface reconstruction package by Nina Amenta, Sunghee Choi, and Ravi Kolluri; the GNU Triangulated Surface Library (GTS); and the Manifold Code (MC) adaptive multilevel finite element kernel by Michael Holst.

Cart3D is a notable example, used for computational fluid dynamics modeling of aerodynamic behavior in aviation and launch vehicle projects at NASA. Cart3D uses my geometric predicates to robustly prepare surface meshes for simulations that integrate many different physical components. Cart3D played a central role in debris transport analysis of the Columbia space shuttle disaster, as illustrated at right. Cart3D's fluid dynamics simulation was a primary means of determining the size, weight, and impact velocity of the foam that damaged the Columbia's wing.

Triangle and Pyramid both use the robust predicates, which are indispensable in allowing them to generate meshes that had previously been unattainable because roundoff errors had caused the mesh generators to fail.

Happily, Nanevski et al. designed and implemented an expression compiler that generates predicates automatically (in the style of my hand-coded predicates).

• My papers on adaptive precision floating-point arithmetic and robust geometric predicates.
• The Robust Predicates page.

A Delaunay triangulation is the most ``natural'' triangulation of a point set. Its enduring popularity arises from many favorable properties such as its tendency to consist of ``nicely rounded'' triangles or tetrahedra, the fact that it is optimal by several useful criteria, the existence of many well-studied algorithms for constructing it quickly, and the fact that modifications are usually local in effect (and thus fast). Delaunay triangulations can be defined in any dimensionality, though the two- and three-dimensional versions are used most heavily in practice.

However, many applications (particularly in graphics, scientific computing, and interpolation) need triangulations of nonconvex objects, or objects that have internal boundaries and discontinuities. For this reason, a two-dimensional geometric construction called the constrained Delaunay triangulation (CDT) was developed in 1986. A CDT is not merely a triangulation of a point set; it is also guaranteed to respect specified segments. Therefore, you can produce a CDT of an object with a complicated nonconvex shape, holes, and interior boundaries (such as the interface between two different materials in a heat conduction simulation). Yet CDTs preserve most of the optimality properties that make Delaunay triangulations so popular.

In 1997, most geometers believed that the CDT could not generalize to three dimensions, however desirable it might be. (CDTs are, in several senses, optimal for piecewise linear interpolation.) A major hurdle is the fact that not all polyhedra can be tetrahedralized without additional vertices. Furthermore, it is NP-hard to determine whether or not a polyhedron can be tetrahedralized without extra vertices, or how many vertices must be added to make it possible.

In 1998, I published a proof that a three-dimensional domain whose boundaries are flat facets (that are not necessarily convex) has a CDT if a certain mathematical condition holds. The figure at right shows a domain and its CDT. This condition is easily tested by a simple algorithm, and the existence result extends to higher dimensions. When the condition does not hold, it can be enforced by a simple algorithm that inserts extra, carefully chosen vertices in a way that guarantees that unnecessarily short edges are not created--thereby making CDTs an ideal foundation for three-dimensional mesh generation. Recall that two-dimensional domains with small angles are difficult to mesh well. Small angles present even greater problems in three dimensions. CDTs provide a way to extend provably good Delaunay mesh generation to domains with small angles (both planar and dihedral angles). Because real-world domains have small angles, this makes all the difference between a merely theoretical algorithm and a practical one.

The next problem was to find a good algorithm for constructing CDTs--one that is both fast and easy to implement. In 2000, I published a sweep algorithm, the first algorithm to be fast enough for practical use. In 2003, I published an algorithm based on bistellar flips; it is simpler and even faster than the sweep algorithm. I am currently incorporating the flip algorithm into Pyramid. I believe the inclusion of CDTs will make the software robust enough for public distribution.

I have also discovered that my existence proofs and construction algorithms can be adapted to a more general class of triangulations known as constrained regular triangulations, or weighted CDTs. Other innovations include simple algorithms for inserting a vertex or facet into a CDT, deleting a vertex or facet from a CDT, and finding a single simplex of a CDT through linear programming. A surprising result, discovered in collaboration with visiting student Nicolas Grislain (of the École Normale Supérieure de Lyon), is that it is NP-hard to determine whether a domain has a CDT, even though it can be determined in polynomial time under the ``general position'' assumption that no five vertices lie on a common sphere!

• My papers on higher-dimensional constrained Delaunay triangulations.

When I began to learn about scientific computing and mesh generation, I wanted to find an article that would explain to me the connection between the geometry of a mesh and the results obtained by an application that uses the mesh. To my surprise, I could not find such an article. After eight years of searching, I decided to write one myself.

When I began the article, I was further surprised to discover that many of the results I wanted to survey do not seem to be in the literature. Hence, my ``survey'' has entailed much original research too, primarily on the error analysis of piecewise linear interpolation applied to isotropic and anisotropic functions over triangles and tetrahedra.

Some of my new results are on quality measures for meshes. A quality measure is a formula used to determine the quality of a triangle or tetrahedron; it is used as an objective function in algorithms for improving meshes by numerical or combinatorial optimization. Strange as it may seem, most publications on quality measures in the mesh generation literature establish no direct mathematical connection between the quality measures they propose and the goals of the application! In the 2002 International Meshing Roundtable, I published an excerpt that derives quality measures directly from the worst-case interpolation error that an element may suffer (and from other considerations, like stiffness matrix condition number). My error bounds for gradients interpolated over triangles and tetrahedra are tighter than any that have previously appeared in the published literature.

The full-length paper, still in progress, is a survey of the relationship between mesh geometry, interpolation errors, stiffness matrix conditioning, and discretization errors, as well as quality measures. It also treats the circumstances where anisotropic meshes are needed, either because a solution function has anisotropic curvature, or because a partial differential equation is inherently anisotropic. These results played a crucial role in defining the problem solved by our anisotropic meshing algorithm.

• My papers on element quality.

• My papers on the Quake Project.
• All papers for the Quake Project.
• The Quake Project page.
• The Archimedes page.

Jonathan Shewchuk