# The QMG Finite Element Solver

The QMG finite element package is intended to be a demonstration of the use of the mesh generator (as well as sparse-matrix operations in Matlab) and is by no means a full-scale finite element package. These routines are available only in Matlab, not in Tcl/Tk. The finite element package computes approximate solutions to boundary value problems of the following form.
div (c grad u) = −f on D
u = g on B1
c·du/dn = h on B2
This equation is called an isotropic linear second-order elliptic boundary value problem with nonconstant coefficients. The scalar field u is the unknown field. Scalar field c is a function of the domain called the conductivity and must be positive at every point, and is given. Scalar field f is a function of the domain called the source term. Domain D is a brep. The boundary of domain D is denoted B, and B is partitioned into two subsets B1 and B2. On B1 we have Dirichlet data given by function g, and on B2 we have Neumann data given by function h.

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 simplicial complex must have been constructed from the brep using the mesh generator. The conductivity and source arguments are optional; see below. Boundary conditions are not arguments to `gmfem`; they are specified using property/value pairs associated with the facets of the domain.

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 conductivity is specified either as a value associated with a region (under property name conductivity), by the third argument to `gmfem`, or by the default conductivity function which is 1 everywhere (Poisson'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 specifies conductivity as a string. This string is passed to `gm_conductivity` by `gmfem`. See the page on user-defined functions for details on defining conductivity functions.

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

Boundary conditions are specified on facets, that is, faces of the brep of dimension d−1 (assuming the dimension of the brep is d). They are specified by a property-value pair associated with the face. The property name is bc. The value is an ordered pair of the form (type func) where type is either d for Dirichlet or n for Neumann, and func is a string defining the function. See the page on user-defined functions.

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

Boundary conditions are specified on internal boundaries (both multiple-region and repeated-boundary internal boundaries) in the same way as on ordinary boundaries. An internal boundary with a Dirichlet condition will force u to take on the prespecified value on that boundary, and hence c·du/dn will probably be discontinuous.

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];
This routine takes a brep whose regions have internal boundaries. Some subset of those internal boundaries is specified by facelist by name; this is a list of strings specifying the names of faces of dimension d−1 that are internal boundaries to be doubled. In matlab, this argument is a cell-array of strings. In Tcl/Tk, it is a list of strings.

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

The BVP is discretized using piecewise linear basis functions. The Dirichlet data is linearly interpolated from nodal values, and the Neumann data is integrated using a low-order quadrature rule that evaluates the Neumann data at mesh node points only. The source data is integrated using a low-order quadrature rule that evaluates the source function at interior and boundary nodes. The conductivity is integrated on interior elements using the midpoint rule, and hence the conductivity function is evaluated only in the interior.

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=ATDA and f0=ATDb, 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.
It should be noted that the main result of that paper, which is an algorithm for accurately solving the finite element linear system Ku=f for BVP's in which c varies by orders of magnitude over the domain, has not been implemented in the current package.

The advantage of assembling the matrix K as ATDA 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.

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