The QMG Mesh Generator

QMG 2.0 provides a quadtree/octree mesh generator. (The Delaunay mesh generator available in QMG 1.1 has been discontinued.) The mesh generator takes as input a brep, that is, a boundary representation of a two- or three-dimensional geometric object and produces as output a triangulation of that brep. The mesh generator operates only on full-dimensional breps (that is, breps whose intrinsic and embedded dimension are equal). The triangulation is stored as a simplicial complex. The mesh generator introduces Steiner points, that is, triangulation vertices that are not necessarily vertices of the original brep. (Indeed, in three dimensions, it is not possible in general to triangulate a domain without introducing Steiner points.)

The mesh generator is invoked as follows in Matlab:

s = gmmeshgen(brep, opt₁, val₁, ..., ,opt_k, val_k);

or in Tcl/Tk:

gmset s [gmmeshgen brep opt₁ val₁ ... opt_k val_k]

The return variable s is a simplicial complex. The first argument to the routine is the brep to triangulate. The remaining arguments are keyword-value pairs. Each opt is a string that specifies a keyword from the list below, and each val is its corresponding value. The keywords may be abbreviated.

Mesh grading control

Two keyword arguments to the mesh generator are used for mesh grading control: userdata and sizecontrol. The mesh grading control function determines an upper bound on how big the mesh elements should be at every point. As the mesh generator executes, it repeatedly calls the function named gm_sizecontrol and gm_sizecontrol_interior. (These functions are in m-files with the same names in Matlab and in file qmg_sizecontrol.tcl for Tcl/Tk.) The arguments to gm_sizecontrol are the real coordinates of the point in the domain, the userdata value (currently unused by the routine) passed into the mesh generator, the sizecontrol value passed into the mesh generator (either using the keyword argument or as a property-value pair on a face), the face containing the point, the geometric entity index within that face and the parametric coordinates within that geometric entity. This function is invoked for points on the boundary. In the case of points in the interior, there is a similar function gm_sizecontrol_interior which has four arguments: real coordinates, the sizecontrol value, the userdata value (currently unused by the routine) passed into the mesh generator, and the topological entity (a region) containing the point.

The two functions gm_sizecontrol and gm_sizecontrol_interior use this information to compute a positive number, which is the maximum allowable size of elements at the input point. The mesh generator uses this value when deciding whether to refine the quadtree/octree. Thus, the sizecontrol keyword option is interpreted by these two functions and controls the mesh upper bound size. The option's default value is (const 1e307), a huge upper bound that essentially means no upper limit is imposed on the mesh size. The quadtree/octree mesh generator uses the number as an upper bound on its largest box size. Therefore, this size control is somewhat imprecise since box sizes occur only in powers of two. See the page on user-specified functions for the details on how the mesh grading control works, along with examples.

If every size control value is constant, then gm_sizecontrol in Matlab or Tcl is not invoked. Constant size control functions are handled entirely within the C++ code for efficiency. A nonconstant size-control function should vary smoothly over the domain, because sharp variations may be missed by the mesh generator. Furthermore, sharp variations may lead to elements with bad aspect ratio. A good way to obtain very small elements in a specific local area in the domain is to put a boundary vertex, edge or surface (internal or external) in the location where you want the small elements and tag it with a sizecontrol property.

As an example of a nonconstant interesting size-control function, take a look at test1 in the test set, which is a size function based a rule in Claes Johnson's book Numerical solutions of partial differential equations by the finite element method, Cambridge University Press, 1987.

Curvature control

The curvecontrol value controls how much the boundary is allowed to deviate from being flat over each simplex in the final mesh. A small value of curvature control forces the mesh generator to refine the mesh at the boundary.

The interpretation of the curvature control argument is as follows. Suppose the value is p. Let x₁, x₂ be two distinct points lying on the same a boundary face F in an octree box. Let n₁, n₂ be the normals at x₁, x₂ if F is 2-dimensional, else let n₁, n₂ be the tangents at x₁, x₂ if F is 1-dimensional. Then if ||n₁−n₂||>p, the box must be split. The norm used in this comparison is like the l² norm but is based on a 26-sided polyhedron rather than a sphere (for efficiency).

Thus, if p is sqrt(2), this allows the boundary to bend by about 90 degrees over one simplex (because two unit vectors at right angles are distance sqrt(2) from each other in the usual 2-norm.) The default value of the curvature setting is 0.5. Values larger than 1.0 are not recommended.

The format for this setting is a string that obeys similar rules to the sizecontrol setting described above. For instance, you can declare (const X) or (formula X) as values for curvature control. You can also modify the gm_curvecontrol function (which is in gm_curvecontrol.m under Matlab and in qmg_sizecontrol.tcl under Tcl/Tk) which is responsible for interpreting these strings. You can also tag individual topological entities with their own curvecontrol property-value pairs. In a situation where multiple curvature control strings apply, the most stringent one determines the curvature control.

For a brep containing a face that is not G¹ (i.e., the normal does not vary continuously over the patches of that face, for instance, because the face is made up of noncoplanar linear patches), there is a minimum positive allowable value for the curvature control. Below this value, the mesh generator will fail. In particular, it will terminate with an error message that the quadtree has been split too finely because a curvature setting could not be met. To help determine this lower bound for your brep, use the gmchecknormals function.

Tolerance

The value of the tol keyword controls the tolerance used by the mesh generator. This tolerance is an indication of how accurate the input brep's geometric data is.

The tolerance is a number between 0 and 1 and is used by QMG to determine when two items are the same and when they are different. The default value of the tol argument is the value of the global variable gm_default_tol in Tcl or GM_DEFAULT_TOL in Matlab. This global variable is initialized to 1e-14 by the Matlab or Tcl/Tk startup routines.

QMG uses double precision arithmetic, so the smallest reasonable setting for this parameter is about 1e-16. But we found that roundoff error quickly builds, so 1e-14 seems reasonable. An error message from the mesh generator like “Ray missed starting patch” indicates tolerance is set too low.

Note that QMG routines automatically take into account the inherent size of the object you are working with and use this number to scale tolerance when necessary. For example, given a cube of side length 1e10 and a tolerance of 1e-11, QMG will assume items in this brep are the same if they are within a distance 1e-1.

Log file

Because the quadtree/octree mesh generator is still in an experimental stage, we have found it very handy to have log of its computations. The value associated with keyword logfile specifies the name of the file in which the log is written. Under Tcl/Tk, the default is mglogNNN, where NNN is the PID of the process running the mesh generator. Under matlab, the default is mglog. Set this value to /dev/null in unix or nul in Windows to prevent generation of the log.

Log file verbosity

The value for verbosity determines how much information to log. This is a nonnegative integer; larger integers mean that more information is logged. For larger 3D problems, mesh generation with the verbosity level of 4 or higher could generate a log that fills up your disk. The default value is 1. The logfile begins with a preamble that confirms the calling arguments to the routine. For verbosity 2 or higher, it traces the creation and deletion of every single box.

If you discover a problem with the mesh generator, please save the log file when you report the problem.

Displaying the gui

The keyword showgui takes either 0 or 1 as a value. Setting showgui to 1 (the default) causes the quadtree/octree mesh generator to put up an informational display as it runs to update you on its progress. This GUI is described below.

Expansion factors

You can specify a vector of expansion factors with keyword expansion. The expansion factors are a decreasing sequence of d+1 numbers, where d is the dimension of the brep. All must be positive and no more than 0.5. In the quadtree/octree mesh generator, boxes are expanded when the algorithm looks for crowding, and these boxes control the degree of expansion. They also serve as upper bounds on the warping distances. Generally these number should be left at their defaults. Increasing these numbers in some cases may improve aspect ratio.

Orientation

You can specify values of either 1 (default) or 2 for keyword orientation. A value of 1 tells the mesh generator to generate all right-handed elements. A value of 2 tells the mesh generator to generate all left-handed elements.

Aspect ratio degrade factors for warping

A tuning parameter aspdegrade should probably be left at default values. This parameter is an array of numbers in the interval (0,1) (three numbers for 2D, four numbers for 3D).

Checking the output

The function gmchecktri checks the output of the mesh generator. The calling sequence is:

Matlab: gmchecktri(b, s {, checko {, tol}}) Tcl/Tk: gmchecktri $b $s {$checko {$tol}}

where b is the brep and s is the simplicial complex. It returns the worst aspect ratio among simplicies in the triangulation, or a negative number if there is an error in a mesh. Because the mesh generator is still experimental, it is recommended that you run this utility after every single run of the mesh generator. The optional argument checko indicates whether you want to check simplex orientation; its value is 0, 1 or 2. If checko=0, then correct orientation of simplices is not checked. If checko=1 (default), then the program checks that each simplex has right-handed orientation. If checko=2, then the program checks that each simplex has left-handed orientation.

This utility performs a series of correctness tests on a mesh generated for a brep. Here is a list of the tests conducted by the gmchecktri program.

Each vertex in the simplicial complex must appear in at least one simplex.
Each boundary vertex must lie in the face of the brep with which it is labeled at the correct parametric coordinates.
Every vertex of the brep must have exactly one vertex of the simplicial complex assigned to it.
Each simplex must have the correct orientation (if checko>0).
Each simplex face must lie on the boundary as specified by the facespec of that simplex face in the mesh. It is not checked whether the face actually lies on the boundary, but it is checked that the vertices of the face all lie on the boundary.
A simplex facet that lies on an exterior boundary must occur in exactly one simplex.
A simplex facet that lies in a region or on an internal boundary must be adjacent to exactly two simplices in the complex, and furthermore, these two simplices must lie on opposite sides of the plane through the facet.

These last two conditions are verified using a hashtable on facets. The gmchecktri program also computes the worst aspect ratio in the simplicial complex. It also computes the globally longest edge length and globally minimum altitude.

How the quadtree/octree mesh generator works

The mesh generator is based on an algorithm due to Mitchell and Vavasis. An early version of this algorithm is described in:

S. A. Mitchell and S. A. Vavasis, “Quality mesh generation in three dimensions,” Proceedings of the ACM Computational Geometry Conference, 1992, pp. 212–221. Also appeared as Cornell C.S. TR 92-1267.

A more recent version of the algorithm is described by the following two papers by Mitchell and Vavasis:

“Quality mesh generation in higher dimensions,” SIAM J. Computing, to appear.
“An aspect ratio bound for triangulating a mesh cut by an affine set,” Proceedings of the ACM Computational Geometry Conference, 1996. See the preprint online in compressed postscript format.

The algorithm is a quadtree/octree algorithm, meaning that it covers the brep with a large box and then subdivides the box into 2^d subboxes, where d is the dimension of the brep (either 2 or 3). The subdivision into boxes continues until no box is crowded. A box is said to be crowded in the following situation. A neighborhood of the box b called ex(b) is constructed; this is a neighborhood in the l¹ norm. The size of the neighborhood is given by the expansion factor multiplied by the size of the box. Let f be the lowest-dimensional face of the brep in ex(b). If there is more than one such f, or if there is any other brep face g in ex(b) such that g is not a superface of f, then b is crowded, else it is not crowded.

The splitting of boxes takes place in d+1 phases. In phase 0, boxes are split if they are crowded and if ex(b) contains a 0-dimensional brep face (a vertex). But if ex(b) does not contain a vertex, then b is made idle until the next phase. The expansion factors vary from phase to phase.

During splitting for separation, the size control and curvature control conditions are also checked.

The operation in the preceding paragraphs is called splitting for separation and there are d+1 phases of this. There are also d+1 phases of alignment which are interleaved with the separation phases. For each phase k, the separation precedes the alignment. The alignment part of phase k operates only on active boxes that have a k-dimensional face of the brep in ex(b). Because of the separation part of the phase, for each such box there is unique brep k-face to which it is associated. The box is said to lie in the orbit of this face. The orbits are processed independently of one another in the alignment phase.

Alignment works as follows. Let f be the face of the brep of dimension k, and suppose we are processing the orbit of f. For a box b in the orbit, some subface c determined to be the close subface to f. Then b is further split so that boxes that also contain c are the same size as b or larger. Once this condition holds for b, an apex is selected on the face of the brep and near the close subface. This apex will end up as a vertex of the triangulation. Once an apex is assigned, the active box is changed to a finished box. This process is carried out for every box in the orbit.

During the separation phase, the boxes do not communicate with one another. In particular, there is no rule enforced concerning the size of adjacent boxes (sometimes called a balance condition). Nonetheless, it is proved in the above papers that there is a bound on how disparate the size of two adjacent can be. (This bound depends on the sharpest angle of the input domain and on the local variation in the size and curvature control function.)

After all d+1 phases of alignment, there are no active boxes remaining and all boxes are finished. Furthermore, there are finished boxes of all possible dimensions 0 to d, and each finished box is linked to the boxes sitting on its surface via pointers. The finished boxes are now triangulated recursively; a k-dimensional box is triangulated by taking the convex hull of its apex with the simplices that triangulate its (k−1)-dimensional surface boxes.

There are many aspects of the algorithm omitted from this overview: boxes are duplicated if the intersection of ex(b) and the brep has more than one component, a table is made of all intersections between quadtree entities and model entities, etc.

The Mitchell-Vavasis algorithm is not the first to use quadtrees for mesh generation. See for instance

M. A. Yerry and M. S. Shephard, “Automatic three-dimensional mesh generation by the modified-octree technique,” Int. J. Numer. Meth. Eng. 20 (1984) 1965–1990.

The Mitchell-Vavasis work is most directly related to 2-dimensional quadtree work by

M. Bern, D. Eppstein and J. R. Gilbert, “Provably good mesh generation,” Proc. 31st IEEE Symp.Found of Computer Sci. (1990) 231–241.

In particular, the two theoretical properties of the MV algorithm described below were first established for the two-dimensional case by this paper.

Theoretical guarantees

The bound on the size of adjacent boxes mentioned above is the basis for a theoretical result concerning the Mitchell-Vavasis algorithm which is as follows. The worst-case aspect ratio of any simplex produced by the algorithm is within a constant factor of optimal for the given domain. Unfortunately, the constant factor is not known explicitly and could be extremely large. The Mitchell-Vavasis algorithm has a second theoretical property: among all possible triangulations of the domain with bounded aspect ratio, the MV algorithm produces the triangulation with the fewest number of simplices up to a constant factor. Again, the constant factor is presumably very large.

The actual implementation QMG 2.0 does not have the same mathematical guarantee as the algorithm in the papers for several reasons. First, the warping distances derived in papers (which are extremely small numbers) have not implemented; instead, some heuristics have implemented. It is not known whether these heuristics work in every case. If the heuristic fails, you will end up with a mesh with bad aspect ratio.

Roundoff error is another reason why the theoretical guarantees may fail. Roundoff error can cause the mesh generator to fail in almost any manner, including sending it into an infinite loop. If you suspect that you have a problem with roundoff, try changing the tol setting.

A third source of difficult is the geometric routines for intersecting rays with curved patches. If these routines fail, then the mesh generator may go into an infinite loop. In current research with G. Jónsson, we have developed an intersection routine that seems to be very robust.

A fourth source of difficulty is curvature in the boundary. The routine currently uses heuristics to detect sharp curves; these heuristics may fail to catch a sharp curve. In this case, an incorrect mesh or a mesh with elements of poor aspect ratio may be generated. The cause of the difficulty is that computing the maximum curvature of a parametric polynomial patch requires (expensive) solution of high-degree polynomial systems and is probably impractical for QMG.

Curvature can cause some subtler problems as well. If the curvature tolerance is fairly high and the expansion factors are small, then it is theoretically possible for the mesh generator to make a simplex with arbitrarily poor aspect ratio.

Performance of the mesh generator

The performance of a mesh generator is measured in three ways: (1) how long does it take? (2) what is the aspect ratio of the elements? (3) how many simplices are generated? These three measures are sometimes correlated for the quadtree/octree algorithm, i.e., when one is bad, so are all three. Naturally, the worst cases in terms of performance are domains with many faces and many sharp angles. The mesh generator works better when there are no sharp angles in the domain.

It works especially well in the case that many of the brep faces are aligned with the coordinate axes. A dramatic example of this axis alignment behavior is provided by test8 in the test set. This test involves meshing a cube. First, the unit cube is meshed in its natural orientation, then it is reoriented randomly and triangulated again. The running time is less than a second for the natural orientation versus about 2 seconds the rotated orientation (Tcl/Tk QMG running on a Pentium Pro). The number of vertices is 27 versus 228. The number of simplices is 48 versus 867. The aspect ratio is 2.7 versus 78.6.

To tune the performance, then you can slightly adjust the tradeoff between mesh quality and number of elements. In particular, by decreasing the curvature tolerance and increasing the expansion factors, you can sometimes improve the mesh quality (but more elements will be generated).

The GUI for the mesh generator

As the generator runs, a panel appears on the screen showing its progress. This panel can be suppressed with the showgui keyword option. The panel is updated once every 1 or 2 seconds. It displays the number of active boxes (which goes up and down as the algorithm proceeds; at the end of phase 3 separation it drops suddenly to zero), the number of finished boxes (which goes up monotonically during the alignment stages, and goes down monotonically to zero during the triangulation phase at the end), the number of vertices (which goes up monotonically during the d+1 phases and then drops once at the very end of the algorithm), and the number of simplices (which goes up monotonically during the triangulation phase).

The control panel also shows the current operation. The order of operations is: Phase 0 separation, phase 0 alignment, phase 1 separation, phase 1 alignment, phase 2 separation, phase 2 alignment, phase 3 separation (if the brep is three-dimensional), phase 3 alignment (if the brep is three-dimensional), triangulation.

During the execution of the mesh generator, a button labeled Cancel appears in the GUI; the user may press this to terminate the mesh generation. Pressing this button terminates mesh generation. Sometimes the button does not respond immediately. The cancel button is polled only when the GUI screen is updated, so if the GUI panel freezes because of some bug in the program, the cancel button won't work either.

When the mesh generation terminates, a second button labeled Dismiss appears in the panel. This button deletes the panel.

The GUI display is controlled by gm_meshgui.m in matlab and gm_meshgui.tcl in Tcl/Tk and hence is fairly easy to customize.

This documentation is written by Stephen A. Vavasis and is copyright ©1999 by Cornell University. Permission to reproduce this documentation is granted provided this notice remains attached. There is no warranty of any kind on this software or its documentation. See the accompanying file 'copyright' for a full statement of the copyright.

Stephen A. Vavasis, Computer Science Department, Cornell University, Ithaca, NY 14853, vavasis@cs.cornell.edu