## The QMG Finite Element Solver |

div (This equation is called an isotropic linear second-order elliptic boundary value problem with nonconstant coefficients. The scalar fieldcgradu) = −fonD

u=gonB_{1}

c·du/dn=honB_{2}

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

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 tou= gmfem(brep, scomplex {, conductivity {, source {, userdata}}});

`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.

`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.

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`

.

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:This routine takes a brep whose regions have internal boundaries. Some subset of those internal boundaries is specified byor:newbrep= gmdouble(brep, facelist);

[Tcl/Tk:newbrep, newmesh] = gmdouble(brep, facelist, mesh);

gmsetor:newbrep[gmdoublebrep facelist]

gmset {newbrep newmesh} [gmdoublebrep 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.

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 ^{T}DA* and

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

The advantage of assembling the matrix *K* as
*A ^{T}DA* 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