div (c*grad u) = -f on DThis 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 polyhedral 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.
u = g on B1
c*du/dn = h on B2
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 only, but the underlying methods can be extended to any dimension.
The remainder of this document 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 {,tol}}}});
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.
Note that in QMG, boundary conditions, conductivity and source terms are all associated with the brep rather than the mesh (unlike some other finite element software packages). Each mesh node has a source facespec field that stores the brep face containing the node. (See the description of simplicial complexes for more information.) The nodes use this field to look up boundary conditions from the brep. 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 u is a Matlab vector with the solution field u to the BVP. The elements of this vector are in correspondence with the nodes of the mesh, and the entry of u is the computed value of the solution field corresponding to that node. (For Dirichlet boundary nodes, the value of u is copied directly from the Dirichlet boundary data.) The routine does not return the value of the flow vector field c*grad u but this would be a straightforward extension.
For 2D problems, the solution u 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 specified conductivity as a string. This string is passed
to gm_conductivity by gmfem.
The details on defining conductivity
functions are here.
For conductivity, the default is "(const 1.0)".
If a face of lower dimension than the top-level has a conductivity property-value pair, it is ignored.
The routine to add property-value pairs, including conductivity, is gm_addpropval.
Dirichlet boundary data means u is specified, whereas Neumann data means c*du/dn is specified. Notice that this system imposes the restriction that 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. If this type of boundary condition is needed, the user should split the facet into subfaces.
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, if a mesh vertex on a low-dimensional face is adjacent
to one or more facets with Dirichlet boundary conditions, then that
mesh vertex will take on the Dirichlet boundary condition that is the
average of the value computed from those facets.
If every facet adjacent to a low-dimensional face has Neumann data, then 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 no local consistency is needed.
If a facet has no boundary condition specified, the default 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. (If the integral
is not zero,
the routine will go ahead and try to solve the problem anyway.)
If the domain is disconnected. then gmfem carries out the transformation described in the previous paragraph on each connected component of the domain individually.
The routine to add property-value pairs, including boundary conditions (bc), is gm_addpropval. To use this routine you must figure out which face of the brep has which index number. 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. This case generally does not arise in practice.
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: gmdoubleib. Its calling sequence is:
Matlab:In Tcl/Tk the return argument l is a 2-element list whose first element is brep2 and whose second element is a mesh2. This routine takes a brep whose top level face or faces have internal boundaries. Some subset of those internal boundaries is specified by faceindexlist; this is a list of integers specifying the indices of the faces of dimension d-1 that are internal boundaries to be doubled. These must be SL internal boundaries; there is no way to double MTL internal boundaries. The routine replaces all the requested 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 brep2.[brep2,mesh2] = gmdoubleib(brep, mesh, faceindexlist {, tol});Tcl/Tk:
gmset l [gmdoubleib $brep $mesh $faceindexlist {$tol}]
The mesh argument to this routine is a mesh generated for the original brep by gmmeshgen. The gmdoubleib routine goes through the nodes of the mesh and duplicates the ones that sit on faces that were duplicated. Then for each simplex it determines which of the duplicated nodes go with that simplex. The resulting mesh the return argument mesh2.
You can now attach boundary conditions to the doubled faces. To help you figure out which of the two doubled faces points in which direction, the newly created doubled internal boundaries are tagged with a property "normal" whose value is the outward-pointing normal of each new face.
The routine will fail (because it can't compute consistent outward normals) if your specified internal boundary faces are arranged in some kind of Moebius strip configuration.
Note that you should not run the mesh generator, gmviz, or any of the geometric modeling routines on the brep produced by gmdoubleib 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 and mesh. Thus, if you have a domain with doubled internal boundaries, you should first create them as single-layer internal boundaries, then generate a mesh, then execute gmdoubleib, then gmchecktri, and then finally gmfem.
For a BVP in which c is constant, f=0, and the solution u is linear, 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'*D*A and f_0=A'*D*b, 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 Dirichlet and on the mesh, f_0 is the contribution to the right-hand side of the assembled linear system Ku=f arising from the Dirichlet data.
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, because we are still searching for ways to make that method more efficient.
The advantage of assembling the matrix K as A'*D*A is that all the operations vectorize well in Matlab. There are no explicit loops over elements or nodes anywhere in the Matlab finite element package, 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).
If you would like to try other orderings, you can explore
the mesh partitioning toolbox,
Chaco or Metis, and then
reorder the matrix with gmrenumber.
Comments in
gmfem explain how you can use these.
QMG can generate hierarchical triangulations in support of a multigrid-type finite element solver, but there is no multigrid routine shipped with QMG. See the gmrefine description.
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.
Stephen A. Vavasis, Computer Science Department, Cornell University, Ithaca, NY 14853, vavasis@cs.cornell.edu