The QMG Mesh Generator |
The mesh generator is invoked as follows in Matlab:
s = gmmeshgen(brep, opt_{1}, val_{1}, ..., ,opt_{k}, val_{k});or in Tcl/Tk:
gmset s [gmmeshgen brep opt_{1} val_{1} ... 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.
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.
The interpretation of the curvature control argument is as follows. Suppose the value is p. Let x_{1}, x_{2} be two distinct points lying on the same a boundary face F in an octree box. Let n_{1}, n_{2} be the normals at x_{1}, x_{2} if F is 2-dimensional, else let n_{1}, n_{2} be the tangents at x_{1}, x_{2} if F is 1-dimensional. Then if ||n_{1}−n_{2}||>p, the box must be split. The norm used in this comparison is like the l^{2} 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^{1} (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.
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.
If you discover a problem with the mesh generator, please save the log file when you report the problem.
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; 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.gmchecktri(b, s {, checko {, tol}})
Tcl/Tk:
gmchecktri $b $s {$checko {$tol}}
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.
gmchecktri
program also computes the worst
aspect ratio in the simplicial complex.
It also computes the globally longest edge length and globally
minimum altitude.
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:
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.
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.
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 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