A Multilevel Approach for Obtaining Locally Optimal Finite Element Meshes
Peter K. Jimack, Rashid Mahmood, Mark A. Walkley and Martin Berzins Computational PDE Unit, School of Computing University of Leeds, Leeds LS2 9JT, UK
Abstract In this paper we consider the adaptive ?nite element solution of a general class of variational problems using a combination of node insertion, node movement and edge swapping. The adaptive strategy that is proposed is based upon the construction of a hierarchy of locally optimal meshes starting with a coarse grid for which the location and connectivity of the nodes is optimized. This grid is then locally re?ned and the new mesh is optimized in the same manner. Results presented indicate that this approach is able to produce better meshes than those possible by more conventional adaptive strategies and in a relatively ef?cient manner.
1
Introduction
Automatic mesh generation is an important computational tool for the ?nite element analysis of a wide variety of engineering problems ranging from structural analysis through to computational ?uid dynamics for example. For many of these problems the use of unstructured grids offers many advantages over structured grids, such as permitting the straightforward triangulation of geometrically complex domains or allowing the mesh density to be adapted according to the behaviour of the solution. In this paper we are concerned mainly with the latter of these properties of unstructured grids: the natural manner in which numerous forms of mesh adaptivity are permitted. Broadly speaking mesh adaptivity algorithms may be categorised as belonging to one of two classes. The ?rst of these, generally referred to as re?nement, involves adding vertices and elements to 1
the mesh in some manner. This may be through local re?nement (e.g. [10]) or through more global remeshing (e.g. [15]) but has the general aim of increasing the number of vertices and elements in those regions of the domain where some measure of the error is unacceptably high. The second class of approach, often referred to as ?re?nement, adapts the mesh in such a way that the number of vertices and elements remains essentially unchanged. This is typically achieved through the use of node movement (e.g. [7, 8]), where the mesh is continuously deformed so as to increase the density of vertices in those regions of the domain with the highest errors, or through the use of edge swapping (e.g. [11]), where the number and position of the vertices is held ?xed but the way in which they are connected together is allowed to change. (There is also a third class of adaptive algorithm, known as ?re?nement, which involves increasing the degree of the ?nite element approximation space on a ?xed mesh, however we do not consider this approach here. See, for example, [1] or [16] for further details.) In this paper we present a new hybrid algorithm that combines the local insertion and movement of vertices with the local swapping of edges in order to attempt to obtain optimal ?nite element meshes for a general class of problem. These are variational problems which may be posed in the following form (or similar, according to the precise nature of the boundary conditions):
?
?? ?? ? ?? ?
??
?? ? ???
?
(1)
for some energy density function
?? ? ?? ? ????
?. Physically this variational form may
be used to model problems in linear and nonlinear elasticity, heat and electrical conduction, motion by mean curvature and many more (see, for example, [2],[7],[8]). Throughout the majority of this paper we restrict our attention to the twodimensional case where dimensions (?
?
?) are considered towards the end.
?, however generalizations to three
For variational problems of the form (1), the fact that the exact solution minimises the energy functional provides a natural optimality criterion for the design of computational grids using ?re?nement. Indeed, the idea of locally minimising the energy with respect to the location of the vertices of a mesh of ?xed topology has been considered by a number of authors (e.g. [5],[14]), as has the approach of locally minimising the energy with respect to the connectivity of a mesh with ?xed vertices (e.g. [11]). Accordingly, the speci?c algorithms that we use for node movement are generalizations of those used in [6] and [14], and the edge swapping is based upon [11]. Full details of these algorithms and how 2
they are combined with local re?nement are presented in the following section. In section 3 it is then demonstrated that combining the above ?re?nement and re?nement approaches in an appropriate manner allows locally optimal grids to be obtained which are better (in terms of energy minimisation) than using either strategy alone. The approach taken is to start with a very coarse mesh which is optimised using ?re?nement. This is then re?ned locally to create a new mesh with a greater number of elements and vertices which can itself be optimised. By repeating this process a number of times a hierarchy of locally optimal meshes is obtained. Since the initial mesh at each level of the hierarchy is produced by local re?nement of an optimal mesh at the previous level it follows that this typically provides a reliable starting point when optimising the new mesh. The results presented demonstrate that the proposed hybrid algorithm is able to provide a mesh which allows the solution of (1) to be approximated to an arbitrary error tolerance using substantially fewer vertices and elements than through re?nement alone. Furthermore, it is also demonstrated that, for a ?xed size of mesh, this multilevel approach invariably ?nds a better locally optimal solution than is obtained by applying ?re?nement directly to a regular starting mesh of the same ?xed size. The paper concludes by addressing possible generalizations of the technique to threedimensional problems and discussing the strengths and weaknesses of the proposed hybrid algorithm for obtaining locally optimal meshes in two dimensions.
2
Multilevel adaptivity
In this section we consider the adaptive ?nite element solution of problems of the form (1), ?rst using
?re?nement and then adding
re?nement to obtain a sequence of optimal meshes. The ?re?nement
approach is described ?rst and consists of a combination of node movement and edge swapping in order to minimize the energy functional for a given size of mesh. At this stage the mesh may be any triangulation of the domain , which is assumed to be a subset of ?? (i.e. ? we only consider piecewise linear ?nite element trial functions. 3
?
? in equation (1)), and
2.1
Locally optimal meshes
We de?ne a locally optimal mesh for the ?nite element solution of (1) to be a mesh with the following properties.
1. There exists some number
? such that if any node is displaced by any distance ?
in any
direction (subject to the constraint that a boundary node remains on the boundary and the domain is not altered) the ?nite element solution on the new mesh has an energy which is no less than the energy of the ?nite element solution on the locally optimal mesh.
2. By noting that each internal edge of a triangulation is shared by exactly two triangles then, if the union of these two triangles is a convex quadrilateral, we may obtain a modi?ed triangulation by swapping the diagonal of this quadrilateral, as shown in Figure 1. The ?nite element solution on any such modi?ed triangulation has an energy which is no less than the energy of the ?nite element solution on the locally optimal mesh.
In order to obtain such a mesh from a given starting mesh we use an approach which is based heavily on that of [14]. This approach combines node movement and edge swapping in a manner which only requires the solution of local problems in order to converge to a global solution of the full problem, (1), on a locally optimal grid. For clarity we describe the node movement and the edge swapping algorithms separately and then discuss how they may be combined. A necessary condition for the position of each node of the triangulation to be optimal is that the derivative of the energy functional with respect to each nodal position is zero. Like the approach of [14] our algorithm seeks to reduce the energy functional monotonically by moving each node in turn until the derivative with respect to the position of each node is zero. Whilst this does not guarantee with absolute certainty that we reach a local minimum (as opposed to a saddle point or even a local maximum), the presence of rounding errors combined with the downhill nature of the technique ensures that in practice any other outcome is almost impossible. The algorithm proposed consists of a number of sweeps through each of the nodes in turn until convergence is achieved. At the beginning of each sweep the gradient, with respect to the position of 4
each node, of the energy functional
?
is found (where
?? ? ?? ?
?
(2)
?
is the latest piecewise linear ?nite element solution). This is done using a slightly
different approach to that described in [14], based upon [8] instead. In [8] it is proved that if × is the position vector of node then
? ?
×
where (
?
? ? ? ?
?
? ?
·
? ?
?
(3)
?
is the usual local piecewise linear basis function at node ,
?
×
is the th component of
×
? to ?),
represents the derivative of
with respect to its ?th argument, other suf?ces represent
components of tensors, ? is the Kronecker delta and repeated suf?ces imply summation ( and
? to ?). Note that using (3) the gradients with respect to all of the vertices in the mesh may be
? to ?
assembled in a single pass of the elements. These gradients are then sorted by magnitude and the nodes visited one at a time, starting with the largest gradient. When each node is visited the direction of steepest descent,
×
? ×
(4)
is used in order to determine along which line the node should be moved. The distance that the node is moved along this line is computed using a onedimensional minimization of the energy subject to the constraint that the node should not move more than a proportion ? (
?
?
?) of the distance from
its initial position to its opposite edge (see Figure 2). Once a new position for the node has been found the value of the solution, ? say, at that node must be updated by solving the local problem
?
?? ??
?
?
?? ? ?? ?
?
(5)
? is the union of all elements which have node as a vertex and Dirichlet conditions are imposed on ? using the latest values for ? . Once this update is complete the same process is undertaken for
Here the next node in the sorted list and when each node has been visited the sweep is complete. Provided convergence has not been achieved the next sweep may then begin. Using the above approach the interior nodes could move in any direction however a slight modi?cation is required for nodes on the boundary of
?. These nodes may only be moved tangentially along
5
the boundary and even then this is subject to the constraint that the domain remains unaltered. Where this constraint is not violated the downhill direction of motion along the boundary is easily computed by projecting × from (4) onto the local tangent of the boundary. The onedimensional minimization in this direction is then completed as for any other node. On Dirichlet boundaries the updated value of ? is of course prescribed however on any other type of boundary it must be computed by solving a local problem of the same form as (5).
Once convergence with respect to the position of each node has been achieved a further reduction in the energy of the solution is sought by the use of edge swapping. Following [11] a loop through each of the internal edges of the mesh is completed and, for each edge, the local energy on the two triangles on either side is computed. The edge is then swapped in the manner illustrated in Figure 1 and the new local energy over the two triangles on either side is computed. If this energy is less than the original local energy on the quadrilateral then the edge swap is kept; otherwise it is rejected. Once the loop through each of the edges has been completed it is repeated until there is an entire pass for which no edges are swapped.
Of course the grid is unlikely to be locally optimal at this point since the edge swapping will generally cause the node locations to become suboptimal. Hence it is necessary to alternate between the node movement and the edge swapping algorithms until the whole process has converged. The downhill nature of each step in the process guarantees that this will eventually occur. Despite this guarantee however, for pragmatic reasons it is useful to be able to impose a small number of additional constraints on the allowable meshes. For example, in our implementation of the edge swapping algorithm parameters are provided for the maximum number of edges that may be connected to a single node and for the smallest angle allowed in any triangle. Similarly, for the node movement algorithm a parameter is provided to specify the smallest area allowed for any element (and any triangle which shrinks to a size below this threshold is removed from the mesh by a simple element/node deletion algorithm). Numerous other such parameters could also be included. 6
2.2
Local mesh re?nement
The main dif?culty with the ?re?nement strategy introduced in the previous subsection is that it is impossible to know a priori how many nodes or elements will be required in order to get a suf?ciently accurate ?nite element solution to any given variational problem. Even an optimal mesh with a given number of nodes may not be adequate for obtaining a solution of a desired accuracy. For this reason some form of mesh re?nement is essential. In [14] global mesh re?nement is combined with ?re?nement and it is demonstrated that this provides better solutions than the use of uniform global re?nement alone. In this work we extend these results in a number of ways. Firstly, by generalizing to systems of equations (i.e.
?
? in (1)) and
secondly, by using local (rather than uniform) re?nement in conjunction with ?re?nement. This, we demonstrate, leads to further ef?ciency gains above and beyond those observed in [14]. In addition, we also demonstrate that the hierarchical approach of starting with a coarse grid and then optimizing, re?ning, optimizing, re?ning, etc., provides a far more robust adaptive algorithm than simply re?ning ?rst and then optimizing the node positions and the mesh topology at the end. For the purposes of this twodimensional work two different local re?nement algorithms have been considered. The ?rst of these divides all triangles which are to be re?ned into four children (as used in [10] for example and illustrated in the top half of Figure 3) whilst the other divides all triangles which are to be re?ned into just two children (as used in [4] for example and shown in the bottom half of Figure 3). In each case any “hanging nodes” left at the end (again see Figure 3) are removed by bisecting the neighbouring elements and then performing local edge swapping. There are many possible ways in which the re?nement might be combined with the ?re?nement to produce a hybrid algorithm. Our experience suggests that a robust approach is to always re?ne an optimized mesh and then to interpolate the coarser solution onto the re?ned mesh as a starting point for the next level of optimization. It also appears to be advantageous to approximately optimize the nodal solution values ?rst before attempting to optimize the positions. The local re?nement itself attempts to subdivide all elements for which the local energy is greater than calculated on any single element. Typically
± of the largest energy value is chosen to be somewhere between ? and ?.
7
In the following section the performance of this hybrid algorithm is assessed using a number of
test problems. In each case comparisons are made with the approach of [14], in which ?re?nement is combined with global re?nement, and with the more conventional approach of using local re?nement on its own.
3
Numerical results
In this section we study three representative test problems in order to assess the quality of the adaptive technique that has been described. The ?rst of these is an arti?cial test case however the second problem is taken from [8] and [13] whilst the third appears in [3], [9] and [14].
3.1
Problem one
We begin by considering the simple twodimensional reactiondiffusion equation
??? · ?? ? ?
subject to the Dirichlet boundary conditions
???
?? ?? ? ?? ??
(6)
?
???
(7)
? whole domain ?. Furthermore, solving equation (6) corresponds to minimizing the energy functional ? ?
? ? ? ? ?
throughout . This boundary condition is chosen so that (7) is also the exact solution of (6) over the
·
?
?
?
(8)
and so this clearly falls into the general class de?ned by problem (1). Indeed, substituting (7) into (8) shows that the energy of the exact solution is given by
? ?? ?
for which
??
?
(9)
For the purposes of these experiments however we restrict our consideration to the single case
? ??,
? ????. ?? elements. This grid is
Initially the problem is solved on a uniform coarse grid containing just
then optimized using the ?re?nement approach of Subsection 2.1 and the total energy reduces from 8
?
? to ? ? ??
, re?ecting the fact that before optimization there were no degrees of freedom in
the boundary layer near ?? produce
?.
Following [14] this optimal grid may now be uniformly re?ned to
elements which may themselves be optimized. This leads to a solution with a total stored
energy of
? ???
. A further global re?nement and optimization then leads to a solution with a total
stored energy of in Figure 4.
? ??
on a mesh of
?? elements: this sequence of locally optimal meshes is shown ?? elements: the ?rst obtained by global re?nement of
Figure 5 illustrates two further meshes of
the initial uniform mesh and the second by optimizing this grid directly. The energies of the solutions on these meshes are
??? ?? and ? ???? respectively thus illustrating, for this example at least, the
superiority of the hierarchical approach when ?re?nement is combined with global re?nement. It is clear from these meshes that although the second grid is locally optimal it suffers from the problem that too many of the degrees of freedom, inherited from the ?rst grid, lie in a part of the domain that is far from the boundary layer. The purpose of this paper however is to propose that the hybrid algorithm should combine ?re?nement with local re?nement and Figure 6 shows a sequence of meshes obtained in this manner.
?? elements, that appears in Figure 4, whilst the second, and ??? elements respectively and were obtained by re?ning third and fourth meshes contain ?, into 2 children only those elements whose local energy exceeded ?± of the maximum local energy on any element. The total energies of the solutions on these four meshes are ? ? , ? ? ? , ? ???? and ? ?? respectively: clearly illustrating the superiority of the use of local rather than global The ?rst mesh is the same one, containing re?nement within the hybrid algorithm. To conclude our discussion of this example we illustrate the advantage of applying the hybrid approach hierarchically by contrasting it with the use of local re?nement alone, possibly followed by a single application of ?re?nement. Figure 7 shows two grids of
??
elements and two grids of
elements that were obtained in this manner (again using a threshold of ment). The total energies of the solutions on the
??
? for the local re?ne? and
element meshes (obtained by local onetofour
re?nement alone and then a single application of the mesh optimization at the end) are
?? ?
respectively, whilst the total energies of the solutions on the 9
element meshes (obtained by
local onetotwo re?nement plus a ?nal optimization) are
? ?
and
?? ?
respectively. We see
that in both cases, despite the fact that the second of each pair of meshes is locally optimal, the quality of these local optima are not as good as that obtained using the hierarchical approach. A summary of all of the computations made for this test problem is provided in Table 1. Elements 32 32 128 512 512 512 32 42 94 323 1048 1048 784 784 Energy 374.473 50.8937 50.1137 50.0158 103.630 50.2311 50.8937 50.3408 50.1010 50.0085 54.8553 50.0536 51.4939 50.0714 Description Figure 4 (top left) Figure 4 (top right) Figure 4 (bottom left) Figure 4 (bottom right) Figure 5 (left) Figure 5 (right) Figure 6 (top left) Figure 6 (top right) Figure 6 (bottom left) Figure 6 (bottom right) Figure 7 (top left) Figure 7 (top right) Figure 7 (bottom left) Figure 7 (bottom right)
Table 1: Summary of the results obtained for Problem one (the global energy minimum is
? ????).
3.2
Problem two
We now consider the more challenging problem of calculating the displacement ?eld for a twodimensional linear elastic model of an overhanging cantilever beam supporting a vertical point load at the end of the cantilever, as illustrated in Figure 8. For this problem ?
?
? and the energy functional is given by
? ?
? ? ?
? ?
??
10
?
? ??
? ×
(10)
Here, all repeated suf?ces are summed from
? to ?, C is the usual fourth order elasticity tensor (in this case corresponding to a Young’s modulus ??? and a Poisson ratio ? ???), provides
represents the traction boundary condition (in this case a
the external body forces due to gravity and
point load as illustrated in Figure 8). The left half of the lower boundary is ?xed whilst the rest of the boundary, say, is free. Unlike for the ?rst example we do not know an exact solution for this problem and so the optimal value for the stored energy is not known a priori. As before we begin by solving the problem on a uniform coarse grid, this time containing elements. This grid is then optimized using the ?re?nement algorithm to reduce the total energy from
?? ???? ? to ?? ? ????. This optimal grid may now be uniformly re?ned to produce ? elements which are also optimized, leading to a solution with a total stored energy of ?? ???? ?. One further global re?nement and optimization then leads to a solution with a total stored energy of ?? ?? on
a mesh of
???
elements. This sequence of locally optimal meshes is shown in Figure 9.
Figure 10 illustrates two further meshes of
???
elements. The ?rst of these is obtained by global
re?nement of the initial uniform mesh and the second by optimizing this grid directly. The energies of the solutions on these meshes are ?
? ?? ? and ?? ?? ?
respectively and so we again observe the
superiority of the hierarchical approach when ?re?nement is combined with global re?nement. As for the previous example, our goal is to assess the hybrid algorithm that combines ?re?nement with local re?nement hence Figure 11 shows a sequence of meshes obtained in this manner. The ?rst mesh is the same one, containing third meshes contain elements, that appears in Figure 9, whilst the second and
??
and
elements respectively and were obtained by re?ning into 2 children
?± of the maximum local energy on any element. The total energies of the solutions on these three meshes are ?? ? ????, ?? ?? ? ? and ?? ? ????
only those elements whose local energy exceeded respectively: again illustrating the superiority of the use of local rather than global re?nement within the hybrid algorithm. We again conclude our example by illustrating the advantage of applying the hybrid approach hierarchically by contrasting it with the use of local re?nement alone, possibly followed by a single application of ?re?nement. Figure 12 shows two grids of that were obtained in this manner (again using a threshold of 11 elements and two grids of
? elements
? for the local re?nement). The
total energies of the solutions on the
element meshes (obtained by local onetofour re?nement
alone and then a single application of the mesh optimization at the end) are ? respectively, whilst the total energies of the solutions on the onetotwo re?nement plus a ?nal optimization) are
? ??
and ?
?? ? ?
?? ??
? element meshes (obtained by local and ?? ? ?? respectively. As be
fore it is clear that the quality of the locally optimal meshes obtained in this manner is inferior to that of meshes obtained using the hierarchical approach. A summary of all of the computations made for this test problem is provided in Table 2. Elements 64 64 256 1024 1024 1024 64 224 455 674 674 462 462 Energy 0.201352 0.253210 0.302353 0.338964 0.306791 0.329249 0.253210 0.308351 0.363313 0.325679 0.342525 0.325879 0.342355 Description Figure 9 (top) Figure 9 (second) Figure 9 (third) Figure 9 (bottom) Figure 10 (top) Figure 10 (bottom) Figure 11 (top) Figure 11 (middle) Figure 11 (bottom) Figure 12 (top) Figure 12 (second) Figure 12 (third) Figure 12 (bottom)
Table 2: Summary of the results obtained for Problem two (the global energy minimum is unknown).
3.3
Problem three
The ?nal twodimensional problem that we consider also involves just one dependent variable (i.e.
?
? in (1)) however it features a solution which is singular at the origin.
12
The energy functional
corresponds to the Laplacian operator and is given by
? ?
disc with a
? ? ? ? ? ?
(11)
where the presence of repeated suf?ces again implies summation from to . The domain, , is the unit
? ?
?
? sector removed, as illustrated in Figure 13, and Dirichlet boundary conditions consistent
with the exact solution ?
??
×??
are applied throughout . Since the exact solution is known in in (11):
?
this case so is the true value of the global minimum of
?? ?
.
As with the previous examples the problem is ?rst solved on a coarse initial mesh, in this case with just
?
elements, which is then optimized. This locally optimal mesh is then re?ned globally and
? ? elements respectively. These ? ?, ? ? ? , meshes are shown in Figure 14 and their corresponding solutions have energies of ? ? ? ? ? and ? ? ?? .
optimized to three further levels, giving meshes of and Once again, it may be observed that the approach of optimizing the mesh at each level after global re?nement is superior to applying global re?nement alone and then optimizing the resulting mesh.
???,
? ? elements, that were obtained by this method. The energies of the solutions on these meshes are ? ? ? (uniform re?nement only) and ? ? (after
Figure 15 shows two meshes, each containing optimization), which are signi?cantly worse than for the ?nal mesh of Figure 14. To conclude this example, we now consider the application of local re?nement in our hybrid algorithm. Figure 16 shows a sequence of four meshes of
? , ?? , ?
and
??
elements respectively.
In order to contrast the solutions on these meshes with those obtained on the meshes shown in Figure 14 we have forced re?nement of each of the edges on the circular boundary so that the domains correspond to the four domains in Figure 14. Further re?nement (one element to two children) has then been permitted locally for any elements whose energy is greater than
?± of the maximum energy on any
single element. This local re?nement is executed repeatedly on each domain until it is necessary to re?ne the boundary elements again. The total energies of the solutions on the four meshes shown in Figure 16 are
?
? ? (the same mesh as in Figure 14), ? ??
,
? ?? ?? and ? ? ? ? respectively.
Again we have seen the advantage of using the hierarchical mesh optimization approach with local, rather than global, re?nement. Furthermore, when local re?nement is used on its own, even if this is followed by mesh optimization, the resulting grids are not as good. Two pairs of such grids, containing 13
??
(onetofour re?nement) and
? ?? (onetotwo re?nement) elements respectively, are illustrated
in Figure 17. For these examples the corresponding ?nite element solutions have total energies of
? ? ?? and ? ? ?? (? ? elements before and after optimization) and ? ??? and ? ? ??? (? ?? elements before and after optimization) respectively. (For the purposes of comparison, we have
arti?cially re?ned those edges on the circular boundary so as to ensure that the domains are identical to the ?nal domains in Figures 14 to 16.) A summary of all of the computations made for this test problem is provided in Table 3.
Elements 28 112 448 1792 1792 1792 28 107 255 1275 1437 1437 1413 1413
Energy 0.549242 0.434828 0.404352 0.396215 0.438164 0.405547 0.549242 0.431777 0.402413 0.395183 0.407613 0.398523 0.402199 0.398123
Description Figure 14 (top left) Figure 14 (top right) Figure 14 (bottom left) Figure 14 (bottom right) Figure 15 (left) Figure 15 (right) Figure 16 (top left) Figure 16 (top right) Figure 16 (bottom left) Figure 16 (bottom right) Figure 17 (top left) Figure 17 (top right) Figure 17 (bottom left) Figure 17 (bottom right)
Table 3: Summary of the results obtained for Problem three (the global energy minimum is
?? ?
).
14
4
4.1
Discussion
Two dimensions
The three examples of the previous section have clearly illustrated that the quality of the ?nal mesh produced when using the proposed hybrid algorithm is better, in the sense that the ?nite element solution has a lower energy, than that obtained by using either re?nement or ?re?nement alone. Furthermore it is demonstrated that combining the mesh optimization with local re?nement is superior to the global re?nement approach used in [14]. Finally, the advantage of using the hierarchical approach, whereby the intermediate level meshes are optimized, is also apparent: an excellent combination of small mesh sizes and low energies for the corresponding ?nite element solutions being achieved. When discussing the merits of our proposed algorithm it is important to note that there are some problems for which the bene?ts may not be quite so substantial as those observed in the three examples above. A common feature to each of these examples is the desirability of clustering the majority of the mesh elements in a relatively small subset of the domain. When a problem is such that the optimal mesh is more uniformly distributed across the domain the local re?nement algorithm will show little or no advantage over the global approach of [14] since almost all elements of the mesh will need to be re?ned when moving from one level to the next. This is a phenomenon that we have observed in at least one example that we have considered (the nonlinear problem used as the second test problem in [14]). Nevertheless, even in this case, our variant of the algorithm performed no worse than that used in [14].
4.2
Three dimensions
Up to this point our discussion has been restricted mainly to the solution of twodimensional problems (i.e. ?
? in (1)). We now conclude the paper by considering how the multilevel approach may also
be applied to obtain optimal tetrahedral meshes when solving problems in three dimensions. The de?nition of a locally optimal mesh in Subsection 2.1 contains two components. One is that the position of each vertex of the mesh should be locally optimal, whilst the other is that the connectivity of these vertices should also be locally optimal. The ?rst of these is quite straightforward to generalize to three dimensions however the edgeswapping part of the de?nition is more complex. This is discussed 15
in some detail in [6] for example, where a number of different stencils are used to modify the local topology of the mesh depending upon how many elements share an edge. As a general rule, if there are elements sharing a particular edge ( a way that replaces them with
?) then the union of these elements may be reconnected in
different elements. This is made even more complicated by the
? ?
fact that there are numerous alternative ways to reconnect the region in this manner, all of which need to be considered when seeking a local optimum. In [6] the objective is just to improve the geometric quality of the mesh and so it is not always necessary to consider all possible edge swaps (they never locally reconnect the mesh when for example). Moreover, because these local reconnections of
the mesh allow the possibility of introducing new elements and edges it is not entirely straightforward to guarantee the termination of an energy minimization algorithm based upon this approach. Due to these dif?culties associated with edge swapping, we restrict this initial discussion on producing locally optimal tetrahedral meshes to the problem of optimizing the node locations only. This means that we will consider a mesh to be optimal if it satis?es the ?rst of the two conditions enumerated in Subsection 2.1. The node movement part of our algorithm then generalizes simply to three dimensions. The derivatives of the energy with respect to the nodal positions may still be computed using (3) with a single loop over the elements of the mesh. This list may then be sorted and, beginning with the largest values of
×
, the nodes may be moved in turn. In each case the movement requires an (approximate) one
dimensional minimization in the direction of steepest descent (given by (4)). As in two dimensions we may also introduce arti?cial constraints on this minimization to prevent the possibility of the mesh from becoming too degenerate. Once the updated location of node has been found it is a simple matter to modify the corresponding solution value through a local solve on the patch of elements, that node. The local re?nement part of the algorithm could either be implemented by dividing each tetrahedron into two children (as in [4] for example) or into eight children (as in [12] for example). It is the latter approach that we use here, and this is illustrated in Figure 18. The removal of “hanging nodes”, which appear when neighbouring elements are at different re?nement levels, is achieved through the use of a transitional re?nement layer, as illustrated in Figure 19. 16
? , surrounding
For a simple test problem we consider the following generalization of the ?rst equation solved in the previous section.
??? · ?? ? ?
subject to the Dirichlet boundary conditions
???
?? ?? ? ?? ?? ? ?? ??
(12)
?
???
(13)
throughout . As with the twodimensional example (13) is the true solution of (12) over all of and the corresponding energy functional ((8) but with this modi?ed
?
?,
? and summation of the repeated ? ?? to yield a thin suf?ces from ? to ?) is has a minimum value given by (9). We again choose ? ????. boundary layer near ?? ? and an optimal energy
Following the approach used for testing the twodimensional algorithm in Section 3, we begin by assessing the performance of threedimensional multilevel mesh optimization when combined with global re?nement. Initially the test problem is solved on a regular coarse grid of
?
tetrahedral
elements, as illustrated in Figure 20. This mesh is then optimized and the total energy of the solution reduces from
?
?
to
??
. Three levels of uniform re?nement, each followed by optimization,
then yield solutions with energies of
? ?
?
,
??
and
?
? on meshes of ?? ?, ?
and
elements respectively. To see that this ?nal mesh is superior to one obtained without multilevel
optimization the problem is then solved on a three level uniform re?nement of the initial mesh shown in Figure 20 (with
? ?
elements therefore), to yield a solution with energy
? ?. When this mesh
is locally optimized however the energy only decreases to a value of
? ? ?. ?
element
We now demonstrate that the potential advantages of using local re?nement with the multilevel optimization also appear to apply in three dimensions. Starting with the locally optimal
grid, a sequence of three further meshes is obtained through local re?nement (again using a threshold of
?) followed by local optimization.
These meshes contain
tetrahedral elements and the corresponding solutions have energies of
? , ? ?? and ??? ?? , ? ? ? and ? ?
respectively. Finally, we demonstrate the superiority of this ?nal mesh over one obtained using only local re?nement followed by local optimization at the end. This comes from the observation that a grid of
??
elements obtained using only local re?nement yields a solution energy of 17
?
and when this is optimized the solution energy only reduces to computational results is provided in Table 4. Elements 384 3072 24576 196608 196608 196608 384 2655 16933 100866 573834 573834 Energy 104.857 59.9077 52.3988 50.7552 67.2790 52.4342 104.857 59.9024 52.3812 50.7460 54.8852 51.3324 Description
? ???
. A summary of all of these
Multilevel optimization and global re?nement.
Global re?nement followed by optimization.
Multilevel optimization and local re?nement.
Local re?nement followed by optimization.
Table 4: Summary of the results obtained for the threedimensional test problem (the global energy minimum is
? ????).
Because of the dif?culties in visualizing very large unstructured tetrahedral meshes we do not include pictures of all of the grids described above. Nevertheless, it is perhaps informative to include a couple of representative examples. Figure 21 therefore shows a mesh of
??
elements created as part
of the above sequence of local re?nements. The solution on this mesh has an energy of contrast to this, Figure 22 shows a locally optimal mesh of
??? ?. In
?
elements, created as part of the se
quence of multilevel optimizations with local re?nement. Although containing many fewer elements than the mesh in Figure 21 the solution on this mesh also has a lower energy of
??
. Despite not
appearing to be particularly smooth, the mesh in Figure 22 certainly seems to possess the excellent qualities of being both ?ne in the direction perpendicular to the boundary layer (near the face ??
?)
and quite coarse in the directions parallel with the layer. It is anticipated that the addition of a suitable 18
edgeswapping strategy will, as in two dimensions, further improve the quality of these meshes.
5
Conclusion
In this paper we have presented a technique for producing ?nite element solutions to variational problems on locally optimal meshes. The major contribution is to propose a multilevel approach which is shown to lead to better quality meshes, with fewer elements, than those obtained by using alternative techniques. Furthermore, based on the ideas presented in [14], there is no need to solve any global problems other than on an initial coarse grid. Extensive numerical results have been presented for a variety of problems in two dimensions and more provisional results have been described for a threedimensional example. All of these numerical experiments have proved to be extremely encouraging. Some additional work is still required however to turn this promising technique into ef?cient, reliable and robust generalpurpose software. For example, the use of edgeswapping has proven to be highly bene?cial in two dimensions and an approach similar to that used (in a different context) in [6] is therefore also likely to be well worth including. In addition, although global solves are not strictly necessary, there may well be ef?ciency gains to be made through the use of approximate global solves at appropriate points in the algorithm (immediately after re?nement for example): these should be investigated carefully. As another example, it is still an open question as to how accurately the mesh needs to be optimized at each level of the hierarchy before local re?nement takes place. Related to this, it also is unclear how accurately it is necessary to solve each of the onedimensional minimization problems that are encountered at each node within each sweep of the node movement algorithm. Other issues that should be considered further concern the importance of the order in which nodes and edges are visited during local optimization sweeps and the possibility of making more aggressive use of the element/node deletion algorithm that is currently only employed when elements shrink to zero.
Acknowledgements
RM gratefully acknowledges the Government of Pakistan for ?nancial support in the form of a merit scholarship. The work of MAW was funded through EPSRC research grant GR/M00077. 19
References
[1] I. Babuska, B.A. Szabo and I.N. Katz, The pVersion of the Finite Element Method, SIAM Journal on Numerical Analysis, vol.18, pp.515–545, 1981.
[2] J.M. Ball, P.K. Jimack and T. Qi, Elastostatics in the presence of a temperature distribution or inhomogeneity, Zeitschrift Fur Angewandte Mathematic Und Physik, vol.43, pp.943–973, 1992.
[3] R.E. Bank, PLTMG Users’ Guide 7.0, SIAM, Philadelphia, 1994.
[4] E. B¨ nsch, An adaptive ?nite element strategy for the 3dimensional timedependent NavierStokes a equations, Journal of Computational and Applied Mathematics, vol.36, pp.3–28, 1991.
[5] M. Delfour, G. Payre and J.P. Zol? sio, An optimal triangulation for secondorder elliptic problems, e Computer Methods in Applied Mechanics and Engineering, vol.50, pp.231–261, 1985.
[6] L.A. Freitag and C. Ollivier Gooch, Tetrahedral mesh improvement using swapping and smoothing, International Journal for Numerical Methods in Engineering, vol.40, pp.3979–4002, 1997.
[7] P.K. Jimack, A best approximation property of the moving ?nite element method, SIAM Journal on Numerical Analysis, vol.33, pp.2206–2232, 1996.
[8] P.K. Jimack, An optimal ?nite element mesh for elastostatic structural analysis problems, Computers and Structures, vol.64, pp.197–208, 1997.
[9] C. Johnson, Numerical solution of partial differential equations by the ?nite element method, Cambridge University Press, Cambridge, 1987.
[10] R. Lohner, An adaptive ?nite element scheme for transient problems in CFD, Computer Methods in Applied Mechanics and Engineering, vol.61, pp.267–281, 1987.
[11] S. Rippa and B. Schiff, Minimum energy triangulations for elliptic problems, Computer Methods in Applied Mechanics and Engineering, vol 84, pp.257–274, 1990. 20
[12] W. Speares and M. Berzins, A 3d unstructured mesh adaptation algorithm for timedependent shock dominated problems, International Journal for Numerical Methods in Fluids, vol 25, pp.81– 104, 1997. [13] B.H.V. Topping and A.I. Khan, Parallel Finite Element Computations, SaxeCoburg Publications, Edinburgh, 1996. [14] Y. Tourigny and F. Hulsemann, A new moving mesh algorithm for the ?nite element solution of variational problems, SIAM Journal on Numerical Analysis, vol. 35, pp.1416–1438, 1998. [15] N.P. Weatherill, Grid adaptation using a distribution of sources applied to inviscid compressible ?ow simulations, International Journal for Numerical Methods in Fluids, vol 19, pp.739–764, 1994. [16] O.C. Zienkiewicz, D.W. Kelly and J.P. Gago, The Hierarchical Concept in Finite Element Analysis, International Journal of Computers and Structures, vol.16, pp.53–65, 1982.
Figure 1: An illustration of the modi?cation of a mesh by the swapping of a single edge. 21
s w
Figure 2: An illustration of local node movement:
× is the direction of steepest descent for the node
motion and ? represents the maximum distance that the node may move in this direction.
Figure 3: An illustration of the re?nement of certain (shaded) elements of a mesh using onetofour subdivision (top) and onetotwo subdivision (bottom).
22
Figure 4: An initial mesh (top left) followed by a sequence of meshes obtained by ?re?nement and then combinations of global re?nement with ?re?nement.
Figure 5: A globally re?ned mesh of
?? elements and the corresponding locally optimized mesh.
23
Figure 6: A sequence of meshes obtained by ?re?nement of an initial coarse grid (top left) and then combinations of local re?nement followed by ?re?nement.
Figure 7: A pair of meshes of
??
elements obtained using local onetofour re?nement (top left) elements obtained using local onetotwo 
followed by optimization and a pair of meshes of re?nement (bottom left) followed by optimization.
24
Figure 8: An illustration of the overhanging cantilever beam with a vertical point load at the end of the cantilever.
Figure 9: An initial mesh followed by a sequence of meshes obtained by ?re?nement and then combinations of global re?nement with ?re?nement.
25
Figure 10: A globally re?ned mesh of
???
elements and the corresponding locally optimized mesh.
Figure 11: A sequence of meshes obtained by ?re?nement of an initial coarse grid and then combinations of local re?nement followed by ?re?nement.
26
Figure 12: A pair of meshes of
elements obtained using local onetofour re?nement (top) fol
lowed by optimization (second) and a pair of meshes of re?nement (third) followed by optimization (bottom).
? elements obtained using local onetotwo
O
45O
Figure 13: An illustration of the domain for the singular problem.
27
Figure 14: A sequence of meshes obtained by ?re?nement of an initial coarse grid (top left) and then combinations of global re?nement followed by ?re?nement.
Figure 15: A globally re?ned mesh of
? ? elements and the corresponding locally optimized mesh.
28
Figure 16: A sequence of meshes obtained by ?re?nement of an initial coarse grid (top left) and then combinations of local re?nement followed by ?re?nement.
Figure 17: A pair of meshes of
??
elements obtained using local onetofour re?nement (top left)
followed by optimization and a pair of meshes of re?nement (bottom left) followed by optimization.
? ?? elements obtained using local onetotwo

29
Figure 18: An illustration of the regular re?nement of a tetrahedron into eight children by bisecting each edge.
Figure 19: An illustration of the transitional re?nement of a tetrahedron when it has a single “hanging node”.
Figure 20: An illustration of an initial uniform mesh containing
?
tetrahedral elements.
30
Figure 21: An illustration of a re?nement alone.
??
element mesh for the solution of equation (12) using local 
Figure 22: An illustration of a algorithm.
?
element mesh for the solution of equation (12) using the hybrid
31