QMG project: mesh generators available

QMG currently provides two mesh generators, a Delaunay-based mesh generator that handles meshes of intrinsic dimension up to 2, and the quadtree/octree mesh generator.

Both mesh generators takes as input a brep, that is, a boundary representation of a polyhedral object, and produce as output a triangulation of that brep. The triangulation is stored as a simplicial complex. Both mesh generators introduce so-called ``Steiner points,'' that is, triangulation vertices that are not necessarily vertices of the original brep. (Indeed, in three dimensions and higher it is not possible in general to triangulate a polyhedral object without introducing Steiner points.) Both mesh generators have a common "front end" routine gmmeshgen, which has the following calling sequence (in Matlab):

s = gmmeshgen(brep, opt1, val1, opt2, val2, ..., optk, valk);

or in Tcl/Tk:

gmset s [gmmeshgen $brep $opt1 $val1 ... $optk $valk]

The return variable s is a chunk format simplicial complex. The first argument to the routine is the brep to triangulate, also in chunk format. The remaining arguments are keyword-value pairs. In other words, opt1,..., optk are strings that specify a keyword, and val1,...,valk are the values for those keyword options. The keywords may be abbreviated. For instance, "dim" may be shortened to "d".

The first three options we describe are common to both mesh generators and are as follows.

The dimension to mesh

The keyword is "dim" tells the mesh generator which dimension of the brep to mesh. In general, with mesh generation there are three different dimensions involved: the embedded dimension of the brep d, the intrinsic dimension of the brep g, and the dimension to mesh m. These dimensions must satisfy the inequalities m<=g<=d.

For instance, if you want mesh on the 2D-faces (the surface) of a 3D brep, you would want m=2. This would be specified in matlab via gmmeshgen(b,'dim',2) or in Tcl/Tk as gmmeshgen $b dim 2.

If you do not specify a dimension, the default is the intrinsic dimension of the brep g.

The frontend program gmmeshgen determines which of the two algorithms (Delaunay or quadtree) to invoke based on the values of m, g, d. The Delaunay mesh generator can handle only those cases when m<=2. The quadtree/octree mesh generator can handle only those cases when m=g=d and 2<=d<=3. (The upper limit of 3 can be removed by slight modifications to the source code; see below.) Thus, the frontend determines which case your call falls into and invokes the correct mesh generator.

Note that not all the cases are covered; for instance, there is no mesh generator in QMG for triangulating the 3-faces of a 4-dimensional brep.

Note also that there is one case where both generators could be applied, namely, the case m=g=d=2, triangulating the area inside a polygon. In this case, gmmeshgen defaults to using the Delaunay mesh generator. You can override this choice by explicitly specifying which algorithm to use via the keyword "algorithm". The value taken by this keyword should be either q or d, for "quadtree" or "Delaunay".

Mesh grading control

Two other keyword arguments to the mesh generator are userdata and size_control, which are for mesh grading control. The mesh grading control function determines an upper bound on how big the mesh elements should be at every point. As the mesh generators execute, they repeatedly call a function gmsize_control (either the m-file or tcl file). The arguments to gmsize_control are a point in the domain, the userdata string, and the size control function. The default behavior is an upper bound of 1e308 on the mesh size. This is a huge number that means essentially that there is no upper bound on the mesh size.

The generators interpret the upper bound in slightly different ways. The Delaunay mesh generator uses the number as an upper bound on its longest triangle side. The quadtree/octree mesh generator uses the number as an upper bound on its largest box size.

Click here for the details on how the mesh grading control works, along with examples.

If you specify a nonconstant size-control function, then it should vary smoothly over the domain, because sharp variations may be missed by the mesh generator. Furthermore, if there is a very sharp variation and the mesh generator does catch it, then it could generate a simplex with very bad aspect ratio.

Note: If you want the mesh generator to use very small triangles in some very local area in the domain, you should accomplish this by putting an internal boundary in the domain where you want the small triangles, and tagging it with a size_control property. This will ensure that the mesh generators do not overlook the sudden variation.

As an example of a nonconstant interesting size-control function, take a look at test3 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.

Suppose you would like a very fine uniform mesh on a domain. There are two approaches: The first is to have the size-control function return a small constant value. The second is to generate a coarse mesh (by using the default size-control function, or using a fairly large constant value in the size-control function), and then to call gmrefine documented below. The latter method is preferable because it is much more efficient.

Tol keyword

The tol keyword controls the tolerance used by the mesh generator; see the writeup on tolerances.

The preceding keyword options apply to both the quadtree/octree mesh generator and the Delaunay mesh generator. The following keyword apply either to one or the other.

Keyword options for the quadtree/octree mesh generator

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. An allowable keyword is "logfile" which specifies the name of the file in which the log is written. The default is mglog. If you don't want a log, set "logfile" to "/dev/null" (under Unix -- I don't know the corresponding setting for Windows).

Verbosity

The keyword "verbosity" determines how information to log. This is a nonnegative integer; larger integers mean that more information is logged. For 3D test cases, if the verbosity level is 4 or higher, then there is a risk that it will fill up your disk. The default value is 2. The logfile begins with a preamble that confirms the calling arguments to the routine. Then there is a section of log providing information about each face of the brep, associated with the face. Next, it logs the quadtree generation process. For verbosity 3 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 an argument. If showgui is set to 1 (the default), then the quadtree/octree mesh generator puts up an informational display as it runs to update you on its progress. This GUI is described below.

Expansion factors

The expansion factors are a decreasing sequence of d+1 numbers. All must be positive except the last, which may be zero. 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, although they may have to be decreased if QMG makes an error. See below.

Keyword options for the Delaunay mesh generator

Minimum angle

The parameter minang specifies the minimum acceptable angle (in degrees) that the Delaunay mesh generator will accept. If you set it 0, then the mesh generator can produce arbitrarily bad triangles. The default is 14.4. If you set it higher than 14.4, the Delaunay mesh generator may go into an infinite loop.

If your domain itself contains a sharp angle, say a 5-degree corner, then the minimum angle requirement is not enforced in that corner.

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. If checko=0, then correct orientation of simplices is not checked. The default is checko=1.

We believe that this utility will verify the correctness of the mesh generator output but we do not have a proof of this conjecture. 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 vertex must lie in the face of the brep with which it is labeled, and furthermore, it must not lie in a proper subface of that face.
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=1).
Let a ``facet'' refer to a (d-1)-dimensional face of a simplex in the triangulation. A facet is exterior if all of its vertices lie on a boundary of the brep (this definition is more complicated in the presence of internal boundaries). Otherwise it is interior.
An exterior facet must be adjacent to exactly one simplex in the complex.
An interior facet 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.

If the gmchecktri program reports that a simplex has the wrong orientation or overlaps on the interior with another simplex, and you are using the quadtree/octree mesh generator, you may be able to fix the problem by decreasing the expansion factors.

If the gmchecktri program reports that a vertex of the brep does not appear in the simplicial complex, then this may mean that your original brep had a dangling face, i.e., a face not adjacent to any top-level face.

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.

The current version is based on the following two papers by Mitchell and Vavasis:

``Quality mesh generation in higher dimensions.'' Click here for a preprint in gzip'd postscript format.
``An aspect ratio bound for triangulating a mesh cut by an affine set,'' Proceedings of the ACM Computational Geometry Conference, 1996. Click here for a preprint in compressed postscript format.

The algorithm is a quadtree 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. The subdivision continues recursively on subboxes. (In the case d=3 the quadtree is known as an octree.) 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 box with the same center as b but a larger diameter (expanded by an expansion factor). 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 f' in ex(b) such that f' 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.

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 comes first, then 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 b' determined to be the ``close subface'' to f. Then b is further split so that boxes that also contain b' 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 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. The triangulation procedure has an optimization to reduce the number of simplices and vertices: under certain conditions, the apex of a box may be replaced by the apex of one of its 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 have more than one connected component (and thus, the separation phases have an expensive geometric connected-component test in the ``inner loop''), ex(b) is recomputed at the beginning of each phase, 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.

How the Delaunay mesh generator works

The Delaunay mesh generator computes an initial triangulation of each 2D face using a non-optimal but reasonably robust O(n^2) algorithm. Then it refines the triangulation until the size bound is met and all angles are at least as big as the minimum requested angle. Refinement is by splitting the longest edge of a triangle. However, if the triangle has an angle greater than 120 degrees, then we move to the neighbor and split the longest edge of the neighbor (unless that neighbor has an angle greater than 120, and so on). Special rules apply to edges on the boundary or adjacent to the boundary. After every edge is split, the algorithm carries out "flipping" to restore the mesh to a Delaunay triangultaion.

This algorithm is similar to an algorithm proposed by Paul Chew of Cornell, except that Chew inserts new points at Delaunay centers instead of longest edges. Longest edge splitting has been proposed by W. Mitchell of NIST and also by M. Rivara of University of Chile. However, the algorithm in QMG is not exactly the same as any of these previous algorithms.

The termination proof for this algorithm has only been partially completed, so you may run into a case where the algorithm gets stuck in an infinite loop. I would like to hear about such a case.

If you use the default size control and you set the minang keyword option to 0, then this mesh generator will not introduce any Steiner points, i.e., every vertex of the mesh it computes will also be a vertex of the input brep. This is how the routine delpt works.

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 QMG1.1 does not have the same mathematical guarantee as the algorithm in the papers for two 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 that has simplices with the wrong orientation or that overlap. You can force the mesh generator to use smaller warping distances by decreasing the expansion factors. I would like to hear about a case where the heuristic fails.

The second reason the theoretical guarantees may fail is because of roundoff error. Roundoff error can cause the mesh generator to fail in almost any manner, including infinite loops. If you suspect that you have a problem with roundoff, try changing the tol setting.

Quadtrees versus Delaunay triangulation

In two dimensions the default mesh generator invoked by gmmeshgen is the Delaunay mesh generator because our experiments show that it is always superior in terms of performance.

Delaunay triangulation extends to three dimensions, but there is no known way to carry out three-dimensional Delaunay triangulation with aspect ratio bounds. Paul Chew at Cornell has some very recent results on this problem not yet written up.

There also some recent work stating that aspect ratio bounds may not be necessary for three dimensions; see for example

G. L. Miller, D. Talmor, S.-H. Teng, and N. Walkington, ``A Delauney Based Numerical Method for Three Dimensions: generation, formulation, and partition, Proc. ACM Symp. Theory of Computing, 683-692, 1995.

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 appeared to be 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 if there are no sharp angles in the domain. It works better for convex regions.

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 test2 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 2 seconds for the natural orientation versus 74 seconds the rotated orientation (Tcl/Tk QMG running on a Pentium). The number of vertices is 27 versus 328. The number of simplices is 48 versus 1305. And the worst-case aspect ratio is about 2.5 for the natural orientation versus about 87 for the random rotation.

High-dimensional mesh generation

The quadtree/octree mesh generation algorithm is written to work for any dimension. However, for efficiency purposes, certain arrays are statically allocated and so the maximum dimension is hard-coded into the program as 3 (in file rect.h). We have not experimented with higher dimensions. If you want to try a higher-dimensional brep, you should change QMG_MAXDIM and related constants in this header file and recompile the program.

Mesh refinement

The QMG package includes a routine gmrefine that refines an existing simplicial complex. In 2D it subdivides every triangle into four smaller triangles, In 3D it subdivides every tetrahedron into eight smaller tetrahedra. Every vertex in the new simplicial complex is either a copy of a vertex from the old complex or lies at the midpoint of two vertices in the old complex. There are several calling formats in Matlab:


mesh2 = gmrefine(mesh);
mesh2 = gmrefine(mesh, brep);
mesh2 = gmrefine(mesh, brep, tol);
[mesh2,vtxsrc] = gmrefine(mesh);
[mesh2,vtxsrc] = gmrefine(mesh, brep);
[mesh2,vtxsrc] = gmrefine(mesh, brep, tol);

In Tcl/Tk, there are three calling formats:


gmset l [gmrefine $mesh]
gmset l [gmrefine $mesh $brep]
gmset l [gmrefine $mesh $brep tol]

In the Tcl/Tk calling sequence, return variable l is a list of two items, the first being a mesh (call it mesh2) and the second being vtxsrc described below.

The new mesh mesh2 will have exactly 2^d times the number of simplices in mesh. It will have approximately 2^d times the number of vertices.

If a brep is specified as an argument, it should go together with mesh argument (in other words, the mesh argument should have been generated by gmmeshgen from the brep). If a brep is not specified, then the mesh coming out of gmrefine will not have any sourcespec information.

This routine can be used to generate fine meshes much faster than gmmeshgen. (In other words, the fast way to make a fine mesh is to generate a coarse mesh with gmmeshgen and then apply gmrefine to the coarse mesh.)

This routine can also be used in connection with a hierarchical finite element method such as multigrid because gmrefine provides information about the origin of each simplex and vertex in the refined mesh in terms of the original mesh. In particular, the numbering of the simplices in mesh2 is set up so that its simplices numbered (i-1)*2^d+1 ... i*2^d are a subdivision of simplex i in mesh (assuming 1-based numbering).

Furthermore, gmrefine can also provide links for vertices. This is the optional second return value vtxsrc. Suppose mesh has n1 vertices and mesh2 has n2 vertices. Then vtxsrc is assigned an n2-by-2 matrix. The i-th row of this matrix is the origin of the i-th vertex of mesh2. This origin is specified by the two integers of the row, say a(i) and b(i), which lie between 0 and n1-1. These integers mean that vertex i (1-based numbering) of mesh2 is the midpoint of vertices a(i) and b(i) (0-based numbering). If a(i)=b(i) then this vertex is a copy of a(i).

If you generate several levels of refinement by calling gmrefine iteratively on a 2D mesh, then at each level the maximum aspect ratio is exactly the same as at the previous level of refinement because all four subtriangles that subdivide an original triangle are similar to the original triangle. In 3D this property doesn't hold. But the following weaker property holds in dimensions d=3 and higher. There is a small constant C (depending on the dimension d) such that if A is the aspect ratio of a simplex T in the original triangulation, then no matter how many times the mesh is refined, the aspect ratio of all the subsimplices of T at any refinement level have aspect ratio bounded by C*A.

Note that although gmrefine has been written in a way that supports multigrid, there are no multigrid routines in QMG.

Using the mesh generator as a standalone application

The mesh generator and checktri routines can run in "standalone" mode under Unix. This means that they can be executed as Unix shell commands. This is accomplished via shell scripts that start QMG as a Tcl/Tk program and feed it a command to read a brep from a file, generate a mesh, write the mesh out to a file, and finish. Thus, in order use the standalone Unix commands, you must first install QMG with the Tcl/Tk scripting option.

The standalone commands are in directory $QMG_ROOT/stdalone. To use these commands, the executable "qmg" (produced by running make to build QMG with a Tcl/Tk front end) must be on your path. In addition, you must set the other relevant environment variables QMG_ROOT, QMG_DATA, TCL_LIBRARY, and TK_LIBRARY according to the installation instructions.

The Unix command for mesh generation has the following form:

meshg brepfile meshfile opt1 val1 opt2 val2 ... optk valk

The file brepfile holds a brep stored in Ascii format and the meshfile stores the resulting mesh in Ascii format. If brepfile is a `-', then the brep is read from standard input. If meshfile is a `-', then the mesh is written to standard output.

There is a Unix-level version of gmchecktri:

checktri brepfile meshfile { checko { tol } }

This routine checks that the simplicial complex in meshfile is a valid triangulation of the brep in brepfile (both in Ascii format). It makes a printout of information and any discrepancies encountered. The mesh file may be a hyphen indicating that the mesh should be read from standard input. The checko and tol arguments may be omitted.

Is the mesh generator deterministic?

The mesh generator is not completely deterministic, i.e., it will not return exactly the same simplicial complex if you run it twice in a row with the same arguments. Here are the sources of nondeterminism:

To generate nullspace bases in matrix factorizations, the mesh generator picks a random matrix and then runs modified Gram-Schmidt on the random matrix. This means that the nullspace bases vary from run to run. The numerical results should be the same, except the last digit or two of positions of the vertices may vary from one run to the next.
To determine whether a point is inside a brep, the mesh generator shoots rays in random directions. However, the containment test should return the same direction regardless of which way the ray is headed, so this should not affect the output.
The mesh generator in several places sorts pointers to boxes according to their memory addresses before processing the boxes. It was found that this sorting improves virtual memory behavior. This introduces nondeterminism into the algorithm because the memory address order is dependent on the state of the operating system's memory allocator at the start of the program.
As far as I know, the order that the boxes are processed does not affect the triangulation, but it does alter the numbering of vertices and simplices. Thus, you should expect two consecutive runs to yield different numberings of vertices and simplices.

The GUI for the quadtree/octree mesh generator

As the quadtree/octree mesh 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-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. If this button is pressed, all the memory allocated by the mesh generator is lost (i.e., it remains allocated but inaccessible) until Matlab or Tcl/Tk is exited and re-entered. 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. The cancel button in Matlab running under Windows works only if the "Enable background processing" menu item has been selected from the "Options" menu in the menubar over the Matlab command window.

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_mggui.m in matlab and gm_meshgui.tcl in Tcl/Tk.

This documentation is written by Stephen A. Vavasis and is copyright (c) 1996 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.

Back to the QMG1.1 home page.

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