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.
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".
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.
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.
If you discover a problem with the mesh generator, please save the log file when you report the problem.
If your domain itself contains a sharp angle, say a 5-degree corner, then the minimum angle requirement is not enforced in that corner.
gmchecktri
checks the output of the mesh generator.
The calling sequence is:
Matlab: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.gmchecktri(b, s {, checko {, tol}})
Tcl/Tk:
gmchecktri $b $s {$checko {$tol}}
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.
An exterior facet must be adjacent to exactly one simplex in the complex.
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.
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:
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.
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.
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.
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.
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.
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.
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:
In Tcl/Tk, there are three calling formats: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 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.gmset l [gmrefine $mesh] gmset l [gmrefine $mesh $brep] gmset l [gmrefine $mesh $brep tol]
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.
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.
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 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.
Stephen A. Vavasis, Computer Science Department, Cornell University, Ithaca, NY 14853, vavasis@cs.cornell.edu