# Locally Optimal Unstructured Finite Element Meshes in Three Dimensions

Locally Optimal Unstructured Finite Element Meshes in Three Dimensions

Rashid Mahmood and Peter K. Jimack School of Computing University of Leeds, LS2 9JT, UK.

Abstract

This paper investigates the adaptive ?nite element solution of a general class of variational problems in three dimensions using a combination of node movement, edge swapping, face swapping and node insertion. The adaptive strategy proposed is a generalization of previous work in two dimensions and is based upon the construction of a hierarchy of locally optimal meshes. Results presented, both for a single equation and a system of coupled equations, suggest that this approach is able to produce better meshes of tetrahedra than those obtained by more conventional adaptive strategies and in a relatively ef?cient manner. Keywords: ?nite elements, variational problems, mesh optimization, tetrahedral elements, node movement, edge swapping, node insertion.

1 Introduction

In this paper we present an extension of our previous work on mesh optimization, presented in [7], from two space dimensions to three. The approach that we follow is to consider the adaptive ?nite element solution of a general class of variational problems using a combination of node movement, edge swapping, face swapping and node insertion. The particular adaptive scheme that is used is based upon the construction of a hierarchy of locally optimal tetrahedral 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. The class of problem that we consider in this work may be posed in the following form (or similar, according to the precise nature of the boundary conditions):

?

?? ?? ? ?? ?

1

??

?? ? ???

?

(1)

for some energy density function ?? ? ?? ? ???? ?. Here ? is the dimension of the problem and ? is the dimension of the dependent variable ?. 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. Throughout this . paper we restrict our attention to the three-dimensional case where ? For variational problems of the form (1), the fact that the exact solution minimizes the energy functional provides a natural optimality criterion for the design of computational grids using ? -re?nement (de?ned here to include both node relocation and mesh reconnection). 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. [2],[14]), as has the approach of locally minimising the energy with respect to the connectivity of a mesh with ?xed vertices (e.g. [12]). All of this work has been undertaken in only two space dimensions however and, to our knowledge, this is the ?rst work in which mesh optimization with respect to the solution energy has been attempted for unstructured tetrahedral meshes in three space dimensions. The algorithm that we use 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 each node, of the energy functional

?

?

?? ? ?? ?

?

(2)

is found (where ? is the latest piecewise linear ?nite element solution). When each node is visited the direction of steepest decent 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 one-dimensional constrained minimization of (2), and once this new position for the node has been found the value of the solution at that node is updated by solving a local problem. Once this update is complete the same process is undertaken for the next node and when each node has been visited the sweep is complete. Provided convergence has not been achieved the next sweep may then begin. 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 and face swapping. In three dimensions there are a large number of different ways in which the local connectivity of the nodes may be altered, see for example [3, 5, 8, 9]. In this work we use the same edge and face swapping stencils as [3, 4], whose work is restricted to improving the geometric quality of the mesh rather than minimizing energy as we do here. Of course the positions of the nodes are likely to be no longer locally optimal at this point due to the edge/face swapping. Hence it is necessary to alternate between the node movement and the swapping algorithms until the whole process has converged (at least approximately). At this stage we allow the application of local mesh re?nement to obtain a new mesh at the next level which must itself now be optimized. The process is complete when either a desired accuracy has been obtained or a max2

Stop = false repeat repeat undertake node optimization undertake connectivity optimization until converged if (accuracy satisfactory) or (maximum mesh size reached) then Stop = true else re?ne mesh solve discrete problem on new mesh end if until Stop Figure 1: Overview of proposed mesh optimization algorithm for the ?nite element solution of (1). imum number of nodes or elements has been reached. Figure 1 illustrates the overall algorithm proposed.

2 Node Movement

A necessary condition for the position of each node of the tetrahedral mesh to be optimal is that the derivative of the energy functional with respect to each nodal position is zero. Like the approaches of [7, 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 a local minimum (as opposed to a saddle point or a local maximum) is reached, the presence of rounding errors combined with the downhill nature of the technique ensures that in practice any other outcome is almost impossible. As indicated above the node optimization phase of the overall algorithm in Figure 1 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 each node, of the energy functional (2) is found. This is done using the same approach as described in [7], based upon [6]. In [6] it is proved that if × is the position vector of node then

×

?

?

? ? ?

?

? ?

·

?

?

?

(3)

where ? is the usual local piecewise linear basis function at node , × is the th to ?), ? represents the derivative of with respect to its ?th component of × (

?

3

argument, other suf?ces represent components of tensors, ? is the Kronecker delta and repeated suf?ces imply summation ( to ? and to ?). Note that using (3) the gradients with respect to all of the vertices in the mesh may be 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 one-dimensional minimization of the energy subject to the constraint that the node should not move more than ? ) of the distance from its initial position to its nearest a proportion ? ( neighbour. 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

?

?

?

?? ??

?

is the union of all elements which have node as a vertex and Dirichlet Here conditions are imposed on using the latest values for ? . All nodes in the sorted list (based upon the magnitude of the gradient in (4)) are updated in this way in turn in order to complete a single sweep of the node optimization step. A number of sweeps are generally taken in order to converge, at least approximately, to an optimal solution. Using the above approach the interior nodes may 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 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 one-dimensional 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).

?

?

?? ? ?? ?

?

(5)

?

?

3 Optimizing Connectivity

In three dimensions tetrahedral mesh connectivities may be altered either by undertaking so-called edge swaps or face swaps. In this work we make use of both of these techniques by exploiting their implementation within the GRUMMP software package, described in [3, 4]. This software seeks to optimize three-dimensional mesh connectivity based upon geometric criteria such as angle conditions and similar qualitative mesh quality measures. Since the source code is publicly available it is possible to modify this in order to undertake optimization of the mesh connectivity based upon 4

our own criteria: speci?cally minimization of the energy functional (2) on the patches of elements surrounding an edge or a face respectively. The two algorithms used for edge and face swapping are now brie?y described.

3.1 Edge Swapping

Edge swapping in three dimensions is not really a swap but a removal of an edge followed by its replacement by one, two or many edges depending upon how many elements surround that edge (see Figure 2 for example). Edge swapping recon?gures the tetrahedra incident on an edge of the mesh by removing that edge and replacing these tetrahedra by ? new tetrahedra. As an example, consider an initial con?guration with ?ve tetrahedra incident to an edge. The left side of Figure 2 shows ?ve tetrahedra incident to an edge OP and the right side shows one possible recon?guration of this sub-mesh into six tetrahedra. This new con?guration is speci?ed by de?ning three “equatorial triangles”, i.e. which are not incident on either of vertices ? and ? . In Figure 2 these triangles are , and . There are four other possible con?gurations for this case (each corresponding to a different set of equatorial triangles), which can be obtained by rotating the interior triangle in Figure 2. As ? tetrahedra, when edge swapping replaces original tetrahedra into more elements are created than are removed. For all of the ?gures in this section solid lines are used to show the front view of the diagram, lines with dashes show the back of the diagram and dotted lines are used in the interior of the convex hull of the points.

?

??

??

?

?

O

1

5 4

T1 = O12P T2 = O1P5 T3 = O2P3 T4 = OP34 T5 = OP45

O

5:6

1

5 4

T1 = O145 T2 = O124 T3 = O234 T4 = P145 T5 = P124 T6 = P234

2

3

2

3

P

P

Figure 2: Edge swapping for 5 tetrahedra to 6, where edge tetrahedra.

?? is surrounded by 5

In addition, the number of possible ways that elements can be reconnected after deleting an edge increases with and is given by

?

?? ? ? ? ? ?? ? ? ??

(6)

(see [5]). When this gives the ?ve possibilities noted in the previous paragraph. However, as grows the number of possible con?gurations grows very rapidly and 5

so, following [3, 4], only edges with are considered as candidates for edge swapping. The possible con?gurations for are shown diagrammatically in Figure 3, where equatorial triangles are shown along with the number of unique rotations for each con?guration. An optimization method therefore has to search through a large number of connectivity permutations for large in order to determine which recon?guration of the original tetrahedra has the lowest energy. For this it is necessary to compute the energy for each tetrahedron in each con?guration. Fortunately, when is large, the number of unique tetrahedra is much smaller than the number of con?gurations times the number of tetrahedra since many tetrahedra appear in more than one con?guration. This is shown in Table 1 (taken from [3]) and means that the cost of performing a local mesh optimization is not quite as high as (6) initially suggests.

O O

4:4

5:6

2

P P

5

O

6:8

6

P

3

3

2

O

7 : 10

7

P

7

7

7

7

7

Figure 3: Equatorial triangles after swapping edge OP, surrounded by 4,5,6 and 7 tetrahedra, including the number of unique rotations for each con?guration shown.

3.2 Face Swapping

Face swapping is cheaper to execute, although possibly more complicated to implement, than edge swapping in three dimensions. It is based upon the possible con?gu6

Tets before 4 5 6 7

Tets after 4 6 8 10

Con?gurations 2 5 14 42

Tets ? con?gs Unique tets 8 8 30 20 112 40 420 70

Table 1: Number of unique tetrahedra and possible con?gurations for edge swapping (taken from [3]).

rations of sets of ?ve distinct non-coplanar points [8, 10] (since each interior face in a tetrahedral mesh separates two tetrahedra, which contain a total of ?ve points between them). These ?ve vertices may be connected to form two, three or four tetrahedra as shown in Figures 4 and 5. The most common con?guration to arise is con?guration in Figure 4, but the others can all occur depending on the geometry of the points ?.

?

O

O

T1 = ABCO T2 = ABCD

T1 = ABDO T2 = BCDO T3 = AOCD

O

T1 = ABCD T2 = ABCO T3 = ABOD T4 = BDCO

B

A B

C

2:3

A B

C

A

C

D D D

Configuration 1A

Configuration 1B

Configuration 2

Figure 4: Possible con?gurations of ?ve points where no four of the ?ve points are coplanar. In the two con?gurations shown in Figure 4, no four of the ?ve points are coplanar. In con?guration the point is in the interior of the convex hull formed by the points , , and ? . Face swapping is possible only if the ?ve points are in con?guration where the triangulation may change from the to or vice versa. In the three con?gurations shown in Figure 5, the points , , and are coplanar, and in con?guration points , and are colinear. Face swapping is possible only if the ?ve points are in con?guration , where con?gurations and are interchangeable. Unlike with edge swapping, where many possible recon?gurations are possible, if a face swap is possible (con?gurations and in Figures 4 and 5 respectively) then only two possible choices need to be compared. This allows a simple and quick comparison to ?nd the one with the lower energy. Details of the way in which the face swapping can be implemented in practice can be found in [9, 10]. In [3, 4] face swapping is the primary algorithm for reconnecting the mesh and edge swapping is

?

?

?

?

?

?

?

7

O

O

T1 = ACBO T2 = ADCO

T1 = ADBO T2 = BDCO

A

2:2

B C

A B C

D

D

Configuration 3A

Configuration 3B

O

O

T1 = ADBO T2 = BDCO

T1 = ADBO T2 = ABCO T3 = BDCO

A B C

A C B

D

D

Configuration 4

Configuration 5

Figure 5: Possible con?gurations of ?ve points where four of the ?ve points are coplanar. used as a supplement to it. The edge swapping routines are also used as part of a separate procedure speci?cally designed to remove poor quality tetrahedra but we do not make use of this procedure in this work since we are motivated only by energy reduction.

4 Node Insertion

The main dif?culty with the node movement and edge/face swapping strategies above 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 this work we use the regular re?nement algorithm implemented in [13]. This divides each tetrahedral element that is to be re?ned into eight children by introducing nodes at the mid-points of each edge. Each new node is then connected to the other 8

two new nodes lying on each face as illustrated in Figure 6. The three new edges on each face may be seen to cut off four child elements at the corners of the parent tetrahedron, leaving an octahedron at the centre. This may be divided into four more child tetrahedra by adding a further edge (LJ in Figure 6) connecting two opposite vertices. The choice of which internal diagonal to insert is important: the approach used in [13] is to choose the longest one but other approaches are possible (see, for example, [11]).

O O

T1 = OIJN T2 = T3 = T4 = T5 = T6 = T7 = T8 = IKAL JBKM NMLC JMKL LIJN LIJK LMJN

1 0

1:8

C I

N

1 0

C

1 0J

1 0

L A B A

1 0

M

1 0

K

B

Figure 6: Regular re?nement of a tetrahedron into 8 child tetrahedra, by bisecting all of the edges. For the results that are presented in the following section both global and local re?nement examples are included. In the former case the regular re?nement algorithm alone is suf?cient however, when local mesh re?nement is used, an additional re?nement scheme is required to deal with the hanging nodes that are left on an unre?ned element which has one or more neighbour that has been re?ned. In [13] these cases are dealt with through the use of a number of so-called green re?nement stencils which deal with elements that have one or more hanging node.

5 Numerical Results

In this section we consider two example problems of the form (1). The ?rst of these is a single equation (i.e. ? ), and the second of these is a system for which ? ? .

?

?

5.1 Problem One

For an initial test problem we consider the following equation:

??? · ?? ? ?

???

???

9

?? ?? ? ?? ?? ? ?? ??

(7)

subject to the Dirichlet boundary conditions

?

(8)

throughout . This is chosen so that (8) is the exact solution of (7) throughout . Hence, for any given value of the analytic solution, and therefore the true energy is chosen and the optimal value for minimum, are both known (in this case the energy is ). Following the approach used in [7] for testing the two-dimensional algorithm, we begin by assessing the performance of three-dimensional multilevel mesh optimization when combined with global -re?nement. Initially the test problem is solved on a regular coarse grid of tetrahedral elements. This mesh is then optimized locally using node movement and edge/face swapping and the total energy of the solution reduces from to . However the number of elements increases from to due to the application of edge/face swapping. Three levels of uniform re?nement, each followed by optimization, then yield solutions with energies of , and on meshes of , and elements respectively. For each of these three levels the number of elements increased by slightly more than a factor of eight due to the edge/face swapping. 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, (with elements therefore), to yield a solution with energy . When this mesh is optimized however the energy only decreases to a value of , with an increase in the number of elements to due to edge/face swapping. We now demonstrate the potential advantages of using local re?nement with the element grid, a semultilevel optimization. Starting with the locally optimal quence of three further meshes is obtained through local -re?nement (by re?ning those elements whose local energy exceeded of the maximum local energy on any element) followed by local optimization. These meshes contain , and tetrahedral elements and the corresponding solutions have energies of , 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 elements obtained using only local -re?nement the observation that a grid of yields a solution energy of and, when this is optimized, the solution energy only reduces to . A summary of all of these computational results is provided in Table 2 and an illustration of the meshes obtained using multilevel optimization with local -re?nement is given in Figure 7.

?

?

? ???? ? ? ?

? ??

? ? ? ? ???? ? ???

? ????

? ? ???

???? ? ?

???

? ?

?

? ?? ?±

? ?? ?

?

???? ? ? ?? ? ? ???? ?

? ?? ? ?

? ? ??

???? ? ???? ? ? ?

5.2 Problem Two

The second problem that we consider involves the calculation of the displacement ?eld for a three dimensional linear elastic model of an overhanging cantilever beam with domain

?

?? ? ?? ?

?

10

?

?

??

?

?

Elements 384 407 3330 27346 220769 196608 197070 407 2931 18741 110170 232140 233506

Energy 378.62763 62.113265 51.223148 50.200687 50.048211 67.278957 52.338504 62.113263 51.226773 50.200292 50.043149 54.813215 51.443760

Description Initial mesh. 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 2: Summary of the results obtained for the ?rst test problem (the global energy minimum is ).

? ????

The bottom half of the beam is ?xed as illustrated by the shaded region in Figure 8 and the energy functional is given by,

? ?

? ? ?

? ?

??

Here, all repeated suf?ces are summed from to , C is the usual fourth order elasticity tensor, chosen to correspond to an isotropic material with a non-dimensionalized Young’s modulus and a Poisson ratio , provides the external body forces due to gravity. The small value of Poisson’s ratio is chosen to ensure that the beam deforms signi?cantly under its own weight. This makes the problem suitable for mesh adaptivity. As before we begin by solving the problem on a uniform coarse mesh, this time containing elements. This mesh is then optimized using the node movement and edge/face swapping algorithms to reduce the total energy from ? to ? . For this particular mesh the edge/face swapping keeps the number of elements same. Three levels of uniform re?nement, each followed by mesh optimization, are undertaken. This produces meshes with , and elements and solutions with energies of ? ,? and ? respectively. and elements. The ?rst of these is We consider two further meshes of obtained by global re?nement of the initial uniform mesh and the second by optimizing this mesh 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 with the previous example, our goal is to assess the hybrid algorithm that combines ? -re?nement with local -re?nement hence we now consider a sequence of

? ?

?

? ?

(9)

???

? ???

??

? ??

?? ?

?? ? ? ?? ? ??

?

?? ? ? ?? ? ??

? ? ???

? ? ??

11

meshes obtained in this manner. The ?rst mesh is the same optimized mesh, containing elements, used as the basis for the global re?nement results. The energy . Four further locally optimal meshes are of the solution on this mesh is ? then obtained, each time via the use of local re?nement (of those elements whose local energy exceeds of the maximum local energy on any element) followed by mesh optimization. These meshes contain , , and elements and yield solutions with energies of ? ,? ,? and ? respectively. 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. We re?ne locally the initial mesh of elements in ?ve levels to achieve a mesh of elements (again using a threshold of for the local re?nement). The total energy of the solution on this mesh is ? . The mesh is then optimized to reduce the total stored energy to ? , with an increased number of elements, , due to edge/face swapping. As before 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 3 and an illustration of the meshes obtained using multilevel optimization with local -re?nement is given in Figure 9.

??

? ??

?±

? ? ??

??

? ? ??

?? ? ? ?? ?

? ? ???

??

?± ? ? ?? ? ? ???

???

???

Elements 192 192 1548 12415 99349 98304 98370 192 958 4529 15315 48403 132698 132958

Energy -0.168295 -0.208546 -0.26773 -0.280849 -0.285704 -0.272196 -0.283207 -0.208546 -0.252279 -0.267699 -0.281052 -0.286102 -0.278015 -0.284321

Description Initial mesh 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 3: Summary of the results obtained for Problem Two (the global energy minimum is unknown).

12

6 Discussion

The two examples of the previous section have clearly illustrated that the quality of the ?nal mesh produced when using the proposed algorithm is better, in the sense that the ?nite element solution has a lower energy, than that obtained by either -re?nement or ? -re?nement alone. Furthermore it is demonstrated that combining the mesh optimization with local -re?nement is superior to combining it with global -re?nement. Finally, the advantage of using the hierarchical approach, whereby intermediate level mesh are optimized, is also apparent: an excellent combination of small mesh sizes and low energies for the corresponding ?nite element solutions being achieved. It should be noted that, although quite complex to implement in -d, the edge/face swapping component of the hybrid algorithm is crucial. This may be demonstrated, for example, by contrasting the results of Table 2 with those obtained for the same test problem but without the connectivity optimization step included in Figure 1. Such modi?ed results are presented in Table 4 and clearly demonstrate the limitations of the adaptive algorithm when edge/face swapping is neglected.

?

Elements 384 384 3072 24576 196608 196608 196608 384 2655 16933 100866 573834 573834

Energy 378.62763 104.85725 59.907732 52.398871 50.755212 67.279033 52.434265 104.85704 59.902412 52.381223 50.746025 54.885230 51.332477

Description Initial mesh. 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 ?rst test problem without edge/face ). swapping (the global energy minimum is

? ????

To conclude this paper we observe that only two numerical examples have been included here and that further work is likely to be required to ensure the robustness of the proposed algorithm for a wide variety of application problems. In particular, it is likely that the mesh re?nement technique used here will be sub-optimal for problems with highly anisotropic solutions, which may well bene?t from a more anisotropic -d re?nement algorithm, such as [1] for example. It is also possible that different criteria could be used for deciding which elements should be locally re?ned (e.g. based upon energy gradients rather than energy values) in order to enhance the technique further.

?

13

Nevertheless, the provisional implementation and results presented here suggest that this approach has signi?cant potential and that further research is indeed likely to be fruitful.

References

[1] E. B¨ nsch, “An adaptive ?nite element strategy for the 3-dimensional timea dependent Navier-Stokes equations”, Journal of Computational and Applied Mathematics, 36, 3–28, 1991. [2] M. Delfour, G. Payre and J.-P. Zol? sio, “An optimal triangulation for seconde order elliptic problems”, Computer Methods in Applied Mechanics and Engineering, 50, 231–261, 1985. [3] L.A. Freitag and C. Ollivier Gooch, “Tetrahedral mesh improvement using swapping and smoothing”, International Journal for Numerical Methods in Engineering, 40, 3979–4002, 1997. [4] C. Ollivier Gooch, GRUMMP users’ guide 7.0”, University of British Columbia, 1999. [5] E.B. de l’Isle and P.-L. George, “Optimization of tetrahedral meshes”, in Modeling, Mesh Generation and Adaptive Numerical Methods for Partial Differential Equations (Springer, Berlin), eds. I.Babuska et al., 97–127, 1995. [6] P.K. Jimack, “An optimal ?nite element mesh for elastostatic structural analysis problems”, Computers and Structures, 64, 197–208, 1997. [7] P.K. Jimack and R. Mahmood, “Multilevel Approach for Obtaining Locally Optimal Finite Element Meshes”, in Developments in Engineering Computational Technology, ed. B.H.V. Topping (Civil-Comp Press), 191–197, 2000. [8] B. Joe, “Three-dimensional triangulations from local transformations”, SIAM Journal of Scienti?c and Statistical Computing, 10, 718–741, 1989. [9] B. Joe, “Construction of three-dimensional improved quality triangulations using local transformations”, SIAM Journal of Scienti?c Computing, 16, 1292– 1307, 1995. [10] C.L. Lawson, “Properties of n-dimensional triangulations”, Computer Aided Geometric Design, 3, 231–246, 1986. [11] M.E. Ong, “Uniform re?nement of tetrahedron”, Siam Journal of Scienti?c Computing, 15, 1994. [12] S. Rippa and B. Schiff, “Minimum energy triangulations for elliptic problems”, Computer Methods in Applied Mechanics and Engineering, 84, 257–274, 1990. [13] W. Speares and M. Berzins, “A 3-d unstructured mesh adaptation algorithm for time-dependent shock dominated problems”, International Journal for Numerical Methods in Fluids, 25, 81–104, 1997. [14] Y. Tourigny and F. Hulsemann, “A new moving mesh algorithm for the ?nite element solution of variational problems”, SIAM Journal on Numerical Analysis, 35, 1416–1438, 1998.

14

Figure 7: An initial locally optimised mesh (top left) followed by a sequence of meshes obtained by combinations of local -re?nement with ? -re?nement for the ?rst test problem.

1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000

Figure 8: An illustration of the overhanging cantilever beam

15

Figure 9: A sequence of meshes obtained by combinations of local -re?nement with ?-re?nement for the second test problem.

16