User-specified functions in QMG

The size-control function is a string. A simple example of a legal size-control string is "(const 0.1)"; this denotes the function that returns the constant value 0.1 (i.e., the desired mesh size is 0.1). There are three places where this size-control string can be specified.

First, it can be attached as a property-value pair to a face of a brep. In this case, the property name is size_control, and the value would be the string "(const 0.1)". A size-control function attached to a face acts only on the part of the mesh near that face. The routine to add property-value pairs, including size_control, is gm_addpropval.

The second place this string is specified is as a keyword argument to gmmeshgen. The keyword name is sizecontrol (no underscore this time). Thus, the user can say (in Matlab)

m = gmmeshgen(b, 'size', '(const 0.1)');

The third place a size-control string is specified is in gmmeshgen.m or gmmeshgen.tcl; in these files the default size control is defined to be "(const 1e308)".

If multiple size control parameters are specified, then all of them act as upper bounds.

Here is an example of using "const" size control functions. Suppose the user types in Matlab:

t = gmpolygon(5);
t = gm_addpropval(t, 1, 0, 'size_control', '(const 0.05)');
m = gmmeshgen(t,'size','(const 0.3)');

gmset t [gmpolygon 5]
gmset t [gm_addpropval $t 1 0 size_control "(const 0.05)"]
gmset m [gmmeshgen $t size "(const 0.3)"]

The first line creates a regular pentagon. The second line attaches a size control property to one edge of the pentagon (the edge whose facespec is 1:0) that requests triangles of size 0.05 or smaller. The default size passed into the mesh generator is size 0.3. Therefore, the mesh generator will generate triangles with side lengths approximately 0.3, except in the neighborhood of the specified edge, where they will be size approximately 0.05. In between the two regions the mesh size will decrease smoothly. Here is how the mesh appears:

Note that for both mesh generation algorithms, the size control parameter acts as an upper bound, but the actual mesh size could be smaller. This is because the mesh generator must take into account geometric contraints imposed by the domain itself. For instance, in the above pentagon example, if we had specified "(const 20)" as the size control argument, clearly we should not expect triangles of side-length 20 as output!

Another valid size control function is an string containing a formula. In this case, the function has the form "(formula expression)". The expression is a usual arithmetic expression containing numbers, +,-,*,/,(), and so on. Special functions like sqrt, log, and so on are also OK. In addition, the expression can contain the special symbols %0, %1, %2, etc, which represent the coordinates of the point at which the size control function is being evaluated (the first coordinate is %0, the second is %1, and so on).

Here is an example of this. Suppose you type in Matlab (a similar example works in Tcl/Tk):

t = gmpolygon(5);
m = gmmeshgen(t, 'size', '(formula 1.02-%0)');

The expression is evaluated either by the Tcl/Tk expr function or the Matlab eval function. These two functions have fairly similar expression syntax; the above example works the same in both expression evaluators. One important difference is exponentiation. In Tcl/Tk you say pow(a,b) whereas in Matlab you say a^b. A second important difference is that Tcl switches to integer arithmetic if it detects all integers in an expression, in contrast to Matlab, which always uses floating-point arithmetic. Thus, you should avoid something like "%0/10" in a Tcl expression because if %0 happens to be 1, then this expression would evaluate to 0. Instead you should say "%0/10.0".

Security note: the Tcl/Tk expr and Matlab eval functions both have the ability to execute arbitrary shell commands. Therefore, if you receive a brep from an untrusted source, check the property/value pairs carefully for undesirable shell commands.

In both Matlab and Tcl, the size control string is processed by a function called gmsize_control. This function makes sure a parenthesized list is passed in, and checks that the first item in the list is the keyword either "const" or "formula". Then it substitutes the actual point coordinates (if the size control is a formula) and computes the value. This function gmsize_control is invoked from with the C++ code of the mesh generator (in file sizectl_fe.cpp; there are two versions of this source file, one for Matlab and one for Tcl.)

The size control scheme is reasonably general, but the user might want to customize it further to suit his/her own needs, for instance, by tying the size control to a posteriori error indicators. In this case the user should rewrite gmsize_control. There is an additional hook to help with customization of the size control. If the user specifies a string to the mesh generator under the keyword userdata (i.e., m = gmmeshgen(t, 'userdata', 'somestring',...)), then this string is saved by the mesh generator and passed verbatim to gmsize_control on every invocation as its "userdata" argument. Right now gm_sizecontrol does not do anything with this argument.

The other functions specified by the user---conductivity, source term, and boundary conditions---are all associated with gmfem (Matlab only). The specification of these is very similar, except there is an alternate form of specifying a formula: "vecformula". If the user specifies a vecformula, then when the formula is eval'd, %0 will have substituted a column vector of x-coordinates instead of a single x-coordinate, and similarly for %1 and so on. (The column vector for %0 will have the same number of entries as the vector for %1). The formula is expected to return a column vector of function values of the same size. The purpose of the vecformula option is to speed up processing in Matlab, since vector formulas are much faster than evaluation by loops. An example of a vector formula would be (vecformula (sqrt(%0.^2 + %1.^2))).

Conductivity is handled by gm_conductivity, which can be customized by the user as mentioned above. The default conductivity is "(const 1.0)". As above, the user can specify conductivities as property-value pairs to top level faces only or as an argument to gmfem. The property name for conductivity is "conductivity". If the user specifies both, then the property-value pair takes precedence.

Source terms are handled by gm_source. The default is "(const 0.0)". Source terms can also be property-value pairs of top level faces. As with conductivity, property value pairs take precedence over specification as an argument to gmfem. The property name for a source term is "source".

Boundary conditions are handled by gm_bdrycond. Unlike conductivity and source, there is no argument to gmfem for boundary conditions; they must be property-value pairs on the faces of the brep of dimension d-1, where d is the dimension of the brep. The property name is 'bc'. The value is the boundary condition, which has a somewhat different syntax. It is specified as an ordered pair <type string> where type is either the letter d or the letter n (for Dirichlet or Neumann) and string is a string of the kind described above, for example, (const 0.1) or a formula (formula %0+%1). The default boundary condition is <n (const 0)>, i.e., du/dn=0, which corresponds to an insulated boundary. Thus, for a 2-dimensional problem, to add a Dirichlet condition of constant 10.0 to face 1:2, you would type

t = gm_addpropval(t, 1, 2, 'bc', '<d (const 10.0)>');

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.