Jens Hugger
The Density Method for the recovery and representation of Finite Element Meshes (Research supported by SNF (The Danish Natural Science Research Council) grant 11-9030).
The goals for the NFEARS project -- Providing computationally cheap finite element solutions, within a given tolerance from the exact solution to nonlinear, parameterized boundary value problems -- and the means to reach the goals, are described in detail. Also it is explained how one particular subgoal of the project has been realized. This subgoal was to construct a theory for finding near optimal finite element meshes for the given class of problems, and representing such meshes in a fashion suitable for the realization of the general NFEARS goals.
The NFEARS project was instigated in the mid 1980's. It has two major goals:
The nonlinear, parameterized boundary value problems make a highly relevant problem class: The reaction of physical systems, and mathematical and engineering models of such, to various external quantities like forces or radiation, is a basic problem in applied mathematics. The mathematical models of such problems are often nonlinear boundary value problems with parameters describing the external quantities. Normally interest is centered around stable systems. It must be certified that a given system remains stable, i.e. that the qualitative behavior does not change, within a pre specified range of the parameters. In some cases however the bifurcation points, where the stability properties of the system change, is the main interest. As an example consider the ``snap-through'' problem where a system goes from one stable state to a qualitatively very different stable or unstable state, as the result of a slight change of the applied forces. Generally this situation is unwanted but in special cases it may be useful. (Think about the small metal frogs that were part of the Japanese industrial miracle of the last half of the 20th century. They had build in small metal plates that when forced to snap-through made a sound somewhat related to that of a frog). For an introduction to the numerical analysis of parameterized nonlinear equations see [1]. For the mathematical theory of bifurcation theory a large amount of literature is available. See for example [2] and [3].
Since the introduction of undetermined parameters makes the problem under determined for numerical methods, parameterized problems are generally solved by separately considering a series of individual parameter or load cases. This obviously makes the numerical solution of parameterized problems a costly affair. On the other hand, the fact that the individual problems in the series of load cases are closely related, and the solutions generally qualitatively the same (except when bifurcation points are encountered), also opens up the possibility of reusing numerically obtained information to simplify the solution process: Assume that a finite element mesh could be described by a small number of real numbers, from here on called Mesh Description Parameters to distinguish them from the External Parameters describing the load case. Assume that there is the same number of mesh description parameters, and that they have the same physical significance for all meshes for a given parameterized problem. Assume further that a few load cases have been considered, say with external parameters on a straight line in the external parameter space, and that the optimal mesh description parameters have been found. Then it is obvious that a simple extrapolation could give a very good indication of the mesh description parameters and thus of the optimal mesh for a nearby load case on the same line in external parameter space. (See [4] for more details). Alternatively, since mesh modification is one of the most expensive parts of obtaining a finite element solution, assume that the mesh has been fixed, i.e. the mesh modification parameters have been held constant, for the previous few load cases. Then extrapolating the (estimated) error to the new load case will reveal whether the error tolerance will be broken by keeping the mesh for another step. These two possibilities would allow for the automatic, adaptive solution process described in algorithm 1 below. (Finer details such as stopping criteria, and how to behave if a bifurcation point or edge is passed, are omitted). In step 1 and 7 storage is provided for previous load cases. Apart from the external parameters defining the load case, we also store the estimated current error (but only for the previous load cases with the same current computational mesh, and the (near) optimal ideal mesh (which is cheap to compute, but expensive to implement). In step 2 and 6 the load case and the computational mesh is selected. If possible, the current mesh is kept from step to step. Only if the extrapolated error from the previous load cases becomes too large is an extrapolation of the ideal meshes of the previous load cases used. Step 3-5 constitutes the main solution loop. (Solve, estimate error, and compute ideal mesh).
Table: Algorithm for the adaptive solution of parameterized boundary
value problems with the finite element method.
Another anticipated situation is that, where during a sequence of solutions of ``less interesting cases,'' where the error tolerance is kept rather high, suddenly an interesting load case is found. It would be interesting in such a situation to be able to increase the precision by decreasing the tolerance, without having to start from scratch. This would be possible if the mesh description parameters were independent of the error tolerance. Then a separate parameter would be needed to take into account the error tolerance.
The two major theoretical chores of the NFEARS project is then
. For
more details see [4], [5], [6], [7], and
[8]. It relies on an error estimator for the local error in
elements of the finite element mesh. The approach to provide a
posteriori local and global error estimation in a wide variety of
norms for the finite element method will not be elaborated here. Some
progress connected to the NFEARS project is described in the papers
[9], [10], [11], [12], [13], [14],
and [15]. Any other approach to local a posteriori error
estimation in arbitrary norms would be applicable however. For the
practical implementation of algorithm 1 the reader is
referred to the NFEARS manuals [16] and [17].
The density method approach falls naturally in four parts: First a theoretical method is derived to represent meshes, and in particular optimal meshes, by functions called Mesh-Size Functions. Then the optimal mesh-size functions are replaced by density functions having the same functionality, together with the property of being independent of the error tolerance. Then a computational version of the method based on local a posteriori error estimation is given. Finally the recovered mesh-size function is discretized by approximating it with the interpolant from a finite dimensional function space, thus giving a few parameter representation. Here we shall concentrate on the first two issues and only comment briefly on the latter two.
For an n-dimensional computational domain 
, we seek
mesh-size functions 
describing a finite
element mesh, in the sense that for any point P in , an
element 
containing P will have approximately the size
. To simplify the
presentation, the treatment will be restricted to two dimensional,
square, computational domains with meshes 
consisting of square elements 
of side length 
and
fixed polynomial degree p. The only meshes allowed are those
obtained from a reference mesh consisting of the entire computational
domain , by iteratively dividing any square element into
four equivalent squares. The introduction of shadow nodes allows for
local mesh refinement for this particular case. Note however that the
density theory applies to much more general cases. The restriction is
solely for ease of exposition. The precise formulation for element
generation can now be based on element area instead of side length,
and we have chosen to construct meshes from a single
mesh-size function h, based on
To ensure the feasibility of this approach we shall require of h that
In practice we keep on refining until all elements satisfy
while
where 
is the
parent element consisting of the union of the four elements
that were created when was created. This approach is
consistent with equation (1) asymptotically, with the added
assumption that 
as 
.
Assuming also that the mesh-size function h is non-negative and
smooth in the interior of all possible elements we get from
equation (1) that 
for some
point 
in , stating that since we consider only square
elements, there exists a point in any element such
that the side length of is 
. (The
smoothness assumption on h implies that singularities may occur only
in points or on edges positioned on fixed element boundaries).
Knowing now how to construct finite element meshes from mesh-size
functions, we turn to the problem of constructing optimal meshes, or
rather Optimal Mesh-Size Functions leading to optimal meshes.
Since optimal meshes has minimal global error for a given number of
elements we need a formula for the global error. We shall measure
error in the standard Sobolev norms 
. Let
be the exact solution, 
the finite element
approximation on the mesh 
, and 
the approximation error for the problem under consideration. We shall
assume that the following local error equations hold:
This error equation is proved for the Lagrange and Taylor
interpolation cases in [14]. Here 
is a point value in
of a known function 
of the 
st
derivatives of , which will be assumed to exist, be bounded,
and be continuous except possibly for a few points and curves in
. Also 
stands for Higher Order Terms in
the element side length. We shall consider only the asymptotical case
where N is so large that we can neglect such higher order terms.
Now by the addition rules for the Sobolev norms, and the smoothness of
and h
Up to higher order terms that we are willing to neglect, the optimal mesh with N elements is then described by
The minimization problem (5) is easily solved with Lagrange multiplier techniques, giving
Note that the Mesh Density Function
can be used for the construction of meshes as well as the mesh-size function (using equation (1)), and that it has the required property of being independent of the error tolerance, by being independent of N, the number of elements.
To determine the Optimal Number of Elements necessary to meet a given error tolerance 
,
by equation (4) we simply (up to negligible higher order
terms) need to solve the minimization problem
The solution is easily seen to be
With the mesh density procedure we obtain the following results with the optimal mesh for the special case considered (neglecting higher order terms):
and
where 
, and C is the constant appearing in equation (6). Also
where 
is the actual number of elements generated by the mesh
construction algorithm, when 
is entered as the wanted number
of elements. For the proofs of these results see the literature
mentioned at the end of section .
Having a local error estimator
for 
for some known mesh , can be approximated by a
piecewise constant function over the elements in using
equation (3). This in turns leads to an
approximation of the optimal number of elements and the mesh density
function. To get a mesh independent representation for the recovered
mesh density function, we interpolate this piecewise constant function
with a smoother function from a fixed finite dimensional function
space. The details are described in the literature mentioned in
section .
