The QMG mesh generator

The mesh generator takes as input a brep, that is, a boundary representation of a polyhedral object, and produces as output a triangulation of that brep. The triangulation is stored as a simplicial complex. The mesh generator introduces 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.) The calling sequence for the mesh generator is

s = gmmeshgen(b, 'size_func', 'logfilename', logfileverbosity, initial_rotation);

The return variable s is a chunk format simplicial complex. The first argument to the routine is the brep to triangulate in chunk format. The remaining arguments are all optional.

Mesh grading control

The second argument to gmmeshgen is a string containing the name of a matlab function accessible on the current matlab path. This function is the mesh grading function and is called repeatedly by the mesh generator as it executes. The mesh grading function takes as input a vector whose coordinates are some point inside or on the boundary of the brep. It should return a positive real number. This real number acts as an upper bound on the size of the boxes generated by the triangulation, which acts as an upper bound (within a constant factor) on the diameter of the simplices produced by the mesh generator. Thus, this function controls local refinement of the mesh.

The boxes generated by the mesh generator will always be no larger than the value of the size function at the evaluation point. However, the boxes could be much smaller: The mesh generator will always produce boxes small enough to fit the constraints imposed by the geometry.

The function 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. The mesh generator evaluates the function at only one point per box.

Smooth variation means that every point of the domain, the norm of the derivative of the size function should be bounded by a small constant multiple of value of the size function.

The default value of this size function is `gm_defaultsize,' which is the following function


function sz = gm_defaultsize(x);
sz = 1e308;

This size function is a do-nothing size function; the number 1e308 is essentially infinity, so this means that every box generated will be smaller than the default size function.

As an example of a more 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.

You can specify your own size-function as noted above in the calling sequence. You can also specify the size function as part of the brep itself; each face may be assigned a property-value pair where the property name is `size_control' and the value is the name of a Matlab function as above.

The order of precedence is: if a property-value pair for property size_control is specified for a face f, it gets precedence. Else the size_control argument to gmmeshgen gets precedence. If neither is specified, gm_defaultsize is used.

The property-value pair for f applies only to those boxes that are in the orbit of f. See below for the meaning of orbit. Each box has a unique orbit.

Note that the QMG mesh generator does not contain any routines for automatic or adaptive local mesh refinement, in which the level of mesh refinement is determined automatically by a posteriori or a priori information from the finite element solution. But the user can do this ``manually'' by writing a size-control function tied to some error indicator.

There are two approaches to making a very fine mesh on a domain. 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 calling gmrefine documented below. The latter method is much more efficient.

Log file

Because the mesh generator is still in an experimental stage, we have found it very handy to have log of its computations. The third argument is a string with the name of the file in which the log is written. The default is `mglog.' The fourth argument is the logfile verbosity level. 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, including the ``tolerances'' 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.

If you want to prevent a log file from being produced, call gmmeshgen with the third argument set to '/dev/null'.

Initial rotation

The mesh generator is based on a quadtree algorithm (see below) which generates boxes with axis-parallel faces. Therefore, the mesh generator will generate a completely different mesh if the object is reoriented with respect to the coordinate axes. The fifth argument is a user-specified orthogonal matrix that is applied to the initial brep, and then its inverse is applied to the resulting simplicial complex before the simplicial complex is returned. If the fifth argument is the string `r' then the mesh generator selects a random rotation. The default value for this argument is the identity matrix.

Checking the output

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

gmchecktri(b,s)

where b is the brep and s is the simplicial complex. There are no return variables. Because the mesh generator is still experimental, it is recommended that you run this utility after every single run of the mesh generator.

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.
The brep faces associated with each vertex of a simplex must form a containment chain in the brep. (In other words, in the d=3 case, if the four vertices of the tetrahedron are associated with brep faces f1, f2, f3, f4, then it must be possible to renumber these four faces so that f1 is contained in f2, f2 in f3, and f3 in f4. Some of them might be equal. For instance, if the tetrahedron is completely interior, all four vertices would be associated with the top-level face of the brep.) This is not a correctness property of triangulations in general, but it is a correctness property of our algorithm.
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.

How the 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,'' forthcoming, and
``An aspect ratio bound for triangulating a mesh cut by an affine set,'' forthcoming.

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. 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 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 can be shown 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 the sharpest local variation in the size control function.)

This bound on the size of adjacent boxes 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.

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''), box subfaces in an orbit all have a rank that is related to tolerance settings, and so on. These aspects will be in the papers.

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 mentioned above were first established for the two-dimensional case by this paper.

Quadtrees versus Delauney triangulation

A popular technique for mesh generation is Delauney triangulation. For instance, Chew's mesh generator is based on Delauney triangulation.

In two dimensions Delauney triangulations are generally regarded as superior to quadtree triangulations, particularly in terms of element quality. For instance, Chew's mesh generator achieves an aspect ratio less than 2*sqrt(3) in all cases (provided that there are no sharp angles in the boundary of the domain).

The difficulty with Delauney triangulation is that nobody knows how to extend it to dimensions higher than 2, with any kind of a bound on the aspect ratio. Furthermore, a systematic way of dealing with domain boundaries in a Delauney triangulation in three dimensions and higher is not known. (In two dimensions the technique for dealing with boundaries is called ``constrained Delauney triangulation'' and an efficient algorithm for CDT is due to Chew.) For this reason we have adopted the quadtree technique, but we hope to do research in the future on 3D Delauney approaches. There is some very recent theoretical work on this topic:

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, 1995, to appear.

Performance of the mesh generator

The performance of the 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?

The worst case I have observed is the following: On three-dimensional test cases with dozens of faces and some very sharp angles, the mesh generator ran for several hours on a Sparcstation. It produced aspect ratios up to 6e5 and over 10^5 simplices with the default size function.

Is this bad? Unfortunately, we don't have another package to compare it against. But for small examples in which we know how to generate a mesh by hand, the mesh generator is often worse than the hand-generated mesh by a factor of 50--1000, both in terms of the number of simplices and the aspect ratio. It appears that the three performance measures are all correlated with each other. If the mesh generator runs a long time, it carries out more needless refining of simplices, and produces a mesh with worse aspect ratio.

The mesh generator also appears to have problems with thrashing virtual memory even though parts of the program were extensively rewritten to improve locality of reference. The degree of thrashing depends how much RAM you have on your system. To determine whether it is thrashing, run the vmstat program at the same time as the mesh generator. If the percentage of CPU time going to the user job drops below 25% while the mesh generator is running, then the program is probably thrashing.

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 27 seconds for the natural orientation versus 145 seconds the rotated orientation. The number of vertices is 36 versus 770. The number of simplices is 80 versus 3244. And the worst-case aspect ratio is about 12 for the natural orientation versus about 400 for the random rotation.

High-dimensional mesh generation

The 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 meshg.H). We have not experimented with higher dimensions. If you want to try a higher-dimensional brep, you should change MG_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. The basic calling format is

mesh2 = gmrefine(mesh);

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.

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, there are also vertex links. In particular, there is a second calling sequence for gmrefine:

[mesh2,vtxsrc] = gmrefine(mesh);

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 without matlab

The mesh generator and accompanying routines are the only programs in the package that can run in stand-alone mode without matlab as a front end. The procedure for installing this software is described elsewhere. The Unix command to call the stand-alone mesh generator is:

meshg brepfile [-l logfile-name] [-v verbosity] [-r rotation] [-s size_function] [-o output file]

There must be space between -l and the log-file name, and similarly for the other options.

The first argument is the name of a file containing a brep in ascii format. Thus, this file should contain a string something like <brep < 3 3 ... > >. This argument can be `-' which indicates standard input. If the Matlab software is installed, this file can be created with the gm_write command.

The logfile-name and verbosity arguments are the same as the matlab calling sequence. The rotation is specified as an orthogonal matrix in row-major form enclosed in curly-braces. For example, you could say -r {0.8 0.6 -0.6 0.8} as a command argument when the brep is two-dimensional.

The size_function argument must be the name of a tcl procedure sitting in file size_control.tcl. This file is shipped with the software and initially contains the default size-control function, but you can add as many functions as you wish. You will need to learn the TCL language to write these functions. The functions take one argument, which is a TCL list of doubles, and should return a double. If the brep faces have a size_control property, this value of this property must also refer to a TCL procedure in this file.

If an output file is not specified, the default is to print the simplicial complex on standard output. The simplicial complex is written into the output file in ascii format. It can be read into the Matlab routines using gm_rd Matlab routine.

There is a Unix-level version of gmchecktri:

checktri brep-file simpcomplex-file

This routine checks that the simplicial complex in simpcomplex-file is a valid triangulation of the brep in brep-file. It makes a printout of information and any discrepancies encountered.

There is a Unix-level version of gmrefine:

refine simpcomplex [-o outputfile] [-v vtxsrc]

The file simpcomplex holds a simplicial complex. If this argument is - then standard input is used. This routine refines the complex as described above and stores the new mesh in outputfile. If no output file is specified, the new mesh is sent to standard output. The file vtxsrc gets the vertex source array (described above). If no such file is specified, then the vertex source array is discarded.

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.
In the case of degeneracy in intersecting ex(b) with the input brep (e.g., a vertex of the brep comes within a very small tolerance distance of a facet of ex(b)), the mesh generator enlarges ex(b) by a small random amount and tries again. It is possible that, from one run to the next, this enlargement is such that in one case a box is crowded but in the other case it is not, causing an additional split. This is expected to be a highly unlikely occurrence.
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 VM 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 mesh generator

As the QMG mesh generator runs, a panel appears on the screen showing its progress. This panel is displayed by the Matlab version only, not the standalone version. 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 monotonically 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 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.

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

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

meshgen.html,v 1.1.1.1 1995/05/05 01:15:54 vavasis Exp

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