The finite element package

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. The finite element package computes approximate solutions to boundary value problems of the following form.
div (c*grad u) = 0 on D
u = g on B1
c*du/dn = h on B2
This equation is called the isotropic nonconstant-coefficient Laplace equation. The scalar field u is the unknown field. Scalar field c is a function of the domain and must be positive at every point, and is given. 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.

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.

A slightly more general problem is the isotropic nonconstant-coefficient Poisson equation in which the first equation is replaced by div (c*grad u) = f, where f is a given function. A solver for this problem has not been implemented yet but this extension to the program would be straightforward.

The call to solve the above problem is

u = gmfem(brep, scomplex, conductivity);
The simplicial complex must have been constructed from the brep using the mesh generator. The conductivity argument is optional; see below. The boundary data is not an argument to the call gmfem and must be specified via property/value pairs associated with the facets of the domain.

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 routines there is no way in the current software to visualize u.

Specifying the conductivity

The conductivity is specified either by a property of a top-level face, by the third argument to gmfem, or by the default conductivity function which is 1 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.)

If a conductivity property-value pair is specified, the property name is `conductivity' and its value is a string that is the name of a Matlab function. This function must be accessible on the matlab path when gmfem runs.

The conductivity function should take one argument; this argument is a matrix with a number of columns equal to the dimension of the brep, and a variable of rows. Each row in this matrix is a point in the interior of that top-level face. The function is expected to return a column vector with the same number of rows as its argument. This column vector should be the values of the conductivities at the input points. Thus, the conductivity function must either loop over elements in the rows of the matrix it receives, or else it must use Matlab's vector operations to construct an entire vector of conductivities at once.

Only conductivity properties associated with top-level faces are used. If a face of lower dimension has a conductivity property-value pair, it is ignored.

As mentioned above, if no property-value pair is specified for a top-level face, then the third argument to gmfem is used for conductivity, which should be a string holding the name of a matlab function as above. If no third argument is specified either, the default gm_one is used, which returns conductivity 1 everywhere.

The gmfem routine confirms which matlab routines it is using for conductivity when it starts.

The routine to add property-value pairs, including conductivity, is gm_addpropval.

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 the character d or the character n for Dirichlet or Neumann, and func is the name of a Matlab function. 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. However, if the need for this type of boundary condition arises, the user can split the face by using the gmsplit operation on the original brep. (As a side effect, the interior and some other faces will get split as well.)

A function specified as a boundary condition must be callable from matlab on the current path. It takes one argument and must return a vector of boundary values. The argument of this function is understood to be a matrix, where each row of the matrix is a vector of coordinates of points that lie on the boundary face. The function must return a column vector whose number of entries is equal to the number of rows in its calling argument, and these are the computed boundary values.

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 Dirichlet boundary conditions and will arbitrarily select one of the Dirichlet functions associated with its Dirichlet facets to get its boundary values. (Therefore, it is recommended that Dirichlet boundary conditions agree at boundaries between facets, otherwise the output of the package could be unpredictable.)

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'') which is computed by function gm_zero. The gmfem routine confirms which boundary conditions it using for each facet.

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

The transformation described in the previous paragraph for pure Neumann domains is not correct in the case of a domain with multiple connected components. The correct way to handle such a domain is to check whether each individual connected component is pure Neumann, and if so, to change one node in each such component to Dirichlet. For the current version of QMG, the user should avoid domains with multiple components in which some of the components have pure Neumann boundary data.

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 is helpful to color the faces with gmrndcolor and display the brep with gmviz or gmgeomview. The coloring of the facets can be displayed with gmshowcolor.

Boundary conditions on internal boundaries

Boundary conditions are specified on internal boundaries (both MTL and SL 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 on an internal boundary is understood as a prespecified jump on c*du/dn as the boundary crossed. The field u itself will be continuous across the internal boundary in this case as well. Thus, the default boundary condition for the package (Neumann data of 0) means that both u and c*du/dn are continuous across the boundary, which means that nodes on the internal boundary obey the original differential equation div (c*grad u) = 0.

There are other boundary conditions for internal boundaries that arise in applications; for instance, if an internal boundary represents a slit that is an insulator inside a conductor, then the boundary condition one would like is c*du/dn=0 on both sides of the internal boundary, and u is discontinuous.

This type of boundary condition is not supported in the current package. Furthermore, extending the package to handle this case is somewhat tricky. Handling this case requires making duplicate copies of all nodes in the relative interior of the internal boundary, and then renumbering the endpoints of the simplices along the internal boundary according to which copy of the duplicated node they own. This duplication procedure has not been implemented.

This type of boundary condition, however, could be approximated with the current package by putting a hole in the domain rather than an internal boundary. The hole would be long in d-1 dimensions and narrow in the last dimension to approximate a slit. There is a lower bound as to how narrow a hole could be; narrower than about 1e-6 times the domain width will confuse the geometric routines inside the mesh generator.

Numerical methods used

The method used to solve this by the package is piecewise linear finite elements. 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 conductivity is integrated on interior elements using the midpoint rule.

For a BVP in which c is constant and in which 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,'' to appear, SIAM J. Num. Anal.
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.

I have done some research on partitioning finite element meshes for efficient nested dissection ordering, but this research is not yet implemented in the current package.

QMG can generate hierarchical triangulations in support of a multigrid-type finite element solver, but there is no multigrid routine shipped with QMG.

This documentation is written by Stephen A. Vavasis and is copyright (c) 1995 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 QMG home page.

fe.html,v 1.1.1.1 1995/05/05 01:15:54 vavasis Exp

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