The QMG Finite Element Solver

div (c grad u) = −f on D
u = g on B₁
c·du/dn = h on B₂

This equation arises in many physical problems including electrostatics, certain kinds of fluid flow, thermodynamics, chemistry, and astrophysics. The finite element package currently handles breps of dimension 2 and 3 and handles only linear elements. The remainder of this page will assume that the reader is at least somewhat familiar with boundary value problems and the finite element method.

The call to solve the above problem is

u = gmfem(brep, scomplex {, conductivity {, source {, userdata}}});

The brep and simplicial complex must both be full-dimensional, i.e., they must have their embedded dimensional equal to their intrinsic dimension.

The return value u is a zba, in other words, subscripting starts at 0. This is because node numbering also starts at 0.

Boundary conditions, conductivity and source terms are all associated with the brep rather than the mesh (unlike some other finite element software packages). This approach has the advantage that the user can change the boundary conditions, source terms, and so on after the mesh has been generated. The disadvantage is that the finite element solver needs to have access to the brep as well as the mesh to look up this data.

The return value represents the solution field u to the BVP. The elements of this vector are in correspondence with the nodes of the mesh: each of vector u is the the computed value of the solution field at the corresponding node. (For Dirichlet boundary nodes, the value of u is copied directly from the Dirichlet boundary data.)

For 2D problems, the solution may be plotted with gmplot. For 3D problems, the boundary values of the solution may be plotted on the boundary of the domain with the same function. In this case the solution u must first be restricted to the boundary using the gmboundary function.

Specifying the conductivity & source terms

The source term is specified either as a value associated with a region (under the property name source) by the fourth argument to gmfem, or by the default source-term function which is 0 everywhere (Laplace's equation). The order of precedence is as given in the previous sentence (i.e., the property-value pair gets highest precedence, etc.) The user-specified source is a string. This string is passed to gm_source by gmfem. See the page on user-specified functions.

A conductivity or source property attached to a face of lower dimension than the top-level is ignored.

Boundary conditions

Dirichlet boundary data means u is specified, whereas Neumann data means c·du/dn is specified. Thus, the same type of boundary condition must be applied to an entire facet, i.e., it is not possible to have a Dirichlet boundary condition on half a facet and a Neumann condition on the other half. The desired effect can be obtained by splitting the brep facet into two subfacets.

Note that boundary conditions are specified only on facets. Boundary conditions for faces of lower dimension (e.g. an edge of a 3D brep) are inferred by gmfem from the facets that contain them. In particular, a mesh vertex on a low-dimensional face will take on the Dirichlet boundary condition equal to the average of the value computed from those facets, assuming at least one facet is Dirichlet.

When every facet adjacent to a low-dimensional face has Neumann data, the nodes on that face will also take on Neumann conditions. Since Neumann data is integrated, such a vertex will combine all the Neumann integrations from the facets that contain it and so local consistency is not needed.

The default boundary condition for faces with no specified bc property is Neumann condition of 0 (i.e., the face is “insulated”).

If all the facets of the domain have Neumann conditions, then the gmfem routine will change one node of the mesh (chosen arbitrarily) to a Dirichlet condition of u=0 to make the problem well-posed. It will also verify that the overall integral of the Neumann boundary conditions are zero. The routine will solve the problem anyway even if the integral is nonzero, but in this case it may issue a warning.

For a disconnected domain, gmfem carries out the transformation described in the previous paragraph on each connected component of the domain individually.

In order to figure out which facet is which (for assigning boundary conditions), it may be helpful to color the faces with gmrndcolor and display the brep with gmviz. The mapping of facet indices to colors can be displayed with gmshowcolor.

Boundary conditions on internal boundaries

A Neumann condition of 0 on an internal boundary (the default) means that the internal boundary does not act as a boundary at all, i.e., the nodes on that boundary are treated like interior nodes.

A nonzero Neumann condition on an internal boundary is treated like a prespecified jump in c·du/dn, whereas u remains continuous across the boundary.

Sometimes it is desirable to have separate boundary conditions on the two sides of an internal boundary. For this purpose, a procedure for doubling internal boundaries is available: gmdouble. It has two calling sequences:

Matlab: newbrep = gmdouble(brep, facelist);
or: [newbrep, newmesh] = gmdouble(brep, facelist, mesh);
Tcl/Tk: gmset newbrep [gmdouble brep facelist]
or: gmset {newbrep newmesh} [gmdouble brep facelist mesh];

These must be repeated-boundary internal boundaries. The routine replaces all the specified internal boundaries with two copies. It may also duplicate lower dimensional faces that interconnect the specified internal boundaries. The resulting brep is the return argument newbrep.

The facelist array can optionally be an asterisk (specified as {'*'} in matlab, since the argument must still be a cell array of strings). In this case, the routine will duplicate all repeated-boundary internal boundaries.

The second calling format takes a mesh as the third argument. This argument is a mesh generated for the original brep by gmmeshgen. The gmdouble routine returns another mesh in which mesh nodes are duplicated if they sit on faces of the brep that were duplicated. The simplices in the mesh are adjusted accordingly. The resulting mesh is returned in argument newmesh.

You can now attach different boundary conditions to the two copies of the doubled faces. To help you figure out which of the two doubled faces points in which direction, you can use the orientation information of the faces.

Note that you should not run the mesh generator on the brep produced by gmdouble because it is an invalid brep (there are two coincident faces). But it is OK—in fact, recommended—that you run gmchecktri on the new brep-mesh pair. Thus, for a domain with doubled internal boundaries, first create them as single-layer internal boundaries, then generate a mesh, then execute gmdouble, then execute gmchecktri, and then finally execute gmfem.

The nodes and simplices of the new mesh are in correspondence with the nodes and simplices of the previous mesh according to the following rules. The order of simplices associated with each face is preserved. The nodes that are duplicated occur together. Duplicates of (global) node number i in the new mesh are numbered i, i+(n+1), i+2(n+1), etc., where n is the maximum global node number in the original mesh. This system prevents any numbering conflicts.

The gmdouble routine uses a tolerance for casting rays to determine which side of the internal boundary contains mesh points. There is no way in the current calling sequence to specify the tolerance, so gmdouble always uses the default tolerance.

Numerical methods used

For a BVP in which c is constant, f=0, the solution u is linear, and the domain is polyhedral, the approximations made in the last paragraph are exact, so the package should return the exact value of u up to rounding error for any combination Dirichlet and Neumann boundary conditions.

The stiffness matrix and Dirichlet boundary data are assembled by forming triple products K=A^TDA and f₀=A^TDb, where A is matrix that depends only on the mesh, and D, a diagonal matrix, depends on the conductivity field. The vector b is depends on the boundary data, source term, and the mesh.

The matrix K is symmetric, positive definite, and sparse. The number of rows and columns of K is equal to the number of nodes in the mesh that are not Dirichlet boundary nodes.

This method for assembling the matrix is nonstandard but is mathematically equivalent to the standard method. (The standard method involves forming ``element stiffness matrices'' and then adding them to get K.) The method is described in detail in section 3 of

S. Vavasis, ``Stable finite elements for problems with wild coefficients,'' SIAM J. Num. Anal. 33 (1996) 890.

The advantage of assembling the matrix K as A^TDA is that all the operations vectorize well in Matlab. There are no explicit loops over elements or nodes anywhere in the Matlab finite element m-files so it performs well.

The linear system Ku=f is solved via the Matlab backslash command applied to the sparse matrix K. By default this is sparse symmetric minimum degree ordering.

The ordering used by gmfem by default is the ordering that emerges from the mesh generator together with the setting of spparms (which defaults to minimum degree).

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.