83_

INTERNATIONAL JOURNAL, FOR NUMERICAL METHODS IN ENGINEERING, VOL. 19, 162 1-1656 A POSTERIORI ERROR ANALYSIS AND ADAPTIVE PROCESSES IN THE FINITE ELEMENT METHOD: PART 11-ADAPTIVE MESH REFINEMENT J. P. DE S. R. GAGO, D. W. KELLY AND 0. C. ZIENKIEWICZ Department o f Civil Engineering, University College of Swansea, U.K. I . BABUSKA Institute for Physical Science and Technology, University of Maryland, U.S.A. SUMMARY This is a paper presented in two parts dealing respectively with error analysis and adaptive processes applied to finite element calculations. Part I contains the basic theory and methods of deriving error estimates for second-order problems. Part I1 of the paper deals with the strategy for adaptive refinement and concentrates again on the p-convergent methods. It is shown that an extremely high rate of convergence is reached in practical problems using such procedures. Applications to realistic stress analysis and potential problems are presented. INTRODUCTION The analysis described in Part I of this paper had as its objective the estimation of the error in a given finite element solution and the identification of the regions for refinement given that the present solution is in some way deficient. For this purpose the concept of error estimators and indicators was developed, the estimators giving a relatively accurate measure of the error in a given mesh and the indicators a measure of the change in energy associated with the introduction of a new variable. In this second part of the paper we shall be concerned with the important engineering problem of achieving a solution of a given accuracy with the minimum total cost. Much effort has been put in the determination of optimal meshes in the sense that the error was minimized for a given number of degrees-of-freedom (see McNeice and M a r ~ a l , ~ ~ and Felippa31.”). With the cost of this optimization being prohibitive, general guidelines, frequently used in practice, were established to obtain efficient meshes (see Oliveira,sl T ~ r c k e , ~ ~ , ~ ~ Melosh et al.44-46 and Shephard et ~ l . ~ ~ , ~ ’ ) . However, the basic problem of achieving a specific accuracy cannot be so addressed and we propose a different path. We start from an arbitrarily defined mesh, using if available any previous experience in solving similar problems, and then estimate the error in the solution. If this error estimate exceeds our error tolerance we will proceed with refinement according to the error indicators. Questions can immediately be raised in connection with the initial starting mesh. At one extreme a sufficiently refined mesh can be presented such that the error tolerance is immedi- ately satisfied and no further refinement is needed. At the other, a very crude mesh can 0029-5981/83/111621-36$03.60 @ 1983 by John Wiley & Sons, Ltd. Received June 1982 Revised September 1982 1622 J. P. DE S. R. GAG0 ET AL. be used as a starting point and a large number of iterations be necessary to achieve the pre-specified error tolerance. The path to follow will obviously depend on the relative costs of salaries and computer runs. At the present stage a reasonable compromise is to use the advantages of previous experience, but in the future a more completely automatic path may be desired. Another question can be raised in relation to the number of degrees-of-freedom introduced at each step. In general we shall choose to add a set with the larger indicators, but the process could be used by introducing a single refinement at a time. The final choice has to be a function of the cost since the more restrictive we are, the more re-solution steps will be necessary to achieve a given accuracy. The important fact is that a process of this nature followed by evaluation of the error at each stage will lead to an arbitrarily reduced error (apart from computer limitations). The process will be stopped when a satisfactory estimate is achieved. By refining the higher indicators a more uniform distribution of indicators occurs at each stage. In the limit we can argue heuristically that the subdivision introduced by a uniform indicator distribution is optimal in the sense of achieving the best solution for the now given number of degrees-of-freedom. This has been proved precisely when subdivision of linear elements in one-dimensional problems is considered (see the fourth section of this paper), but for other types of discretization the optimality is only intuitively evident. However, the proof of optimality is not central to the adaptive process since the process will be stopped when a solution of a prescribed accuracy has been achieved economically. We include this proof since this property reinforces and justifies the power of the proposed scheme, and contributes to its understanding. It should also be realized that this error equilibration is associated with a smoothing of the different singularity effects over the mesh. This implies rates of convergence not influenced by the singularities. Thus optimality will be achieved in the sense that each element will behave at its full approximation power and, if there is not a restriction of polynomial order, arbitrarily high rates of convergence can be obtained. In the previous discussion we have not specified the type of refinement because the statements are completely general. There are, however, basically two types of refinement-the h and p versions of the finite element method. The h version of the finite element method is defined as a refinement by subdivision of the already existing finite elements into elements of the same type, and the p version as the refinement by adding successive higher order finite element trial functions. The rates of convergence for the h and p versions in terms of the number of degrees-of-freedom (non- adaptive situations) have been identified and quantified in Babuska and Szabo. l6 Obviously, when refining an already existing finite element mesh there are other possibilities apart from h and p. For example, intuitively, a combination of strategies may provide better convergence than the pure h or pure p version defined above, and the u priori error analysis in Babuska and Dorr4 proves that this combination does in fact provide better rates of convergence. Another alternative strategy consists in the introduction of a special shape function known to model approximately the behaviour of the exact solution in the expansion of the unknown variable. This is not new and was shown to work in where a ‘local-global’ formulation of the finite element method was presented. Here the name ‘global’ derives from the fact that the special shape functions can be based over the whole domain, and although in Mote47 this concept was not used in relation to adaptive refinement schemes, the extension to this case is immediate. A POSTERIORI ERROR ANALYSIS-PART I1 1623 However, the concept of adding shape functions modelling the behaviour of the exact solution is not restricted to global formulations. This class of formulation also includes all the special-purpose finite elements such as the ones used in fracture mechanics, where either shape functions modelling the expected behaviour of the singularity are introduced or node positioning is such that by isoparametric mapping the same effect may be achieved. The results presented in this paper have been obtained through single h- or p-convergence programs. These programs are called adaptive since each step depends on information provided by the previous ones, i.e. the error indicators and the error estimators. We remark that, although closely related, these two measures are independent of one another. This means that we can have an adaptive program without any error estimation, and vice versa. But the combination of both techniques provides the best approach. The programs are called self- adaptive in the sense that no user interaction is necessary to activate the adaptive process and, quite generally, the p-convergence program shows a better capacity to obtain the desired accuracy with the minimum degrees-of-freedom. This is a consequence of better rates of convergence introduced by p-refinement. We have chosen the number of degrees-of-freedom to characterize the rates of convergence since it allows a direct comparison between h- and p-refinement processes. Also, as this number is related to memory requirements and the dimension of the system of equations being solved, it also provides a measure of the cost associated with a given solution. Note, however, that cost should not on every occasion be the conditioning factor since a more expensive process might be appropriate if we want to ‘fully’ optimize the number of degrees-of- freedom. Finally, we want to point out that the initial mesh can be created either by a ‘normal’ mesh generator or by the self-adaptive program itself by specifying a certain number of a priori h-convergence refinements on the given initial subdivision. For a very thorough review of the different types of mesh generators available, we refer to Buel and Bush2’ and Thacker.” We do, however, defend that there is no reason for an elaborate a priori mesh generation since only an a posteriori mesh generator can define a mesh that is directly related to the solution of the given problem. In p-refinement results it will be shown that a mesh designed taking into consideration experience in solving similar problems is rewarding on the basis that the computer has to spend less effort to achieve a given accuracy and avoids high stress oscillations near singularities. This is related to the better rates of convergence presented by combinations of h- and p-refinements. Before implementing an adaptive algorithm we must define our a posteriori error measures. We will discuss them in the following two sections. We must also consider what type of refinement we are going to allow the program to perform. We will discuss this in the fourth section. Following this, we will consider one-dimensional examples and then the two- dimensional case. Finally, we will draw some conclusions. CRITERION FOR REFINEMENT The essential question to be considered in this section is the identification of the adaptive principles on which the decision for refinement will be based. The common principle to both h- and p-refinements is to add those degrees-of-freedom whose contribution to the increase of global accuracy is maximal. This approach has clear intuitive reasoning. This principle is closely related to the equilibration of the error indicators mentioned earlier, as it will obviously lead to this result. 1624 J. P. DE S. R. G A G 0 ET AL In Part I of this paper we named the first of the a posteriori error measures an error indicator. The main idea behind the error indication definition was to consider that our present solution G may be improved by the addition of a hierarchic mode Nh. There is no restriction whether this hierarchic mode is going to be a global or a local mode, or if it is going to correspond to an h- or a p-refinement. In this case we assume the error in our solution to be where Nh is the possible new refinement and Lib, the related variable. Knowing that the energy of the error is given by its internal product with the residual forces (r) corresponding to a given solution, our prediction of the energy absorbed by this hierarchic mode will be 35.36.42.89 where ch will be defined in the following. As we are not interested in solving the total problem in order to evaluate the i?h coefficient, we simply admit it to be directly proportional to the corresponding force and inversely proportional to the corresponding stiffness: 35,36,42*89 This results in We will then compare the available Nh and decide on the ones that give the greatest error decrease. This decision has been based on a parameter y and we refine all the degrees-of- freedom corresponding to error indicators 7 * such that 77 * 2 Y v m a x ( 5 ) This parameter y controls the speed of convergence and we have used y = 0.25 and y = 0-5 in our p-convergence examples. Although affected by the choice of y, the rate of convergence will in all cases be optimal in the sense of cancelling the singularity effects. The h-version does not have a hierarchic structure. We have used an alternative set of indicators obtained by keeping track of the error estimators given by successive stages of refinement. If the error for the original element of size ho is E: and when halved is E ? , then assuming that we can write e =ch" (6) where h is the parameter representing the mesh size at a given stage, we have (7) (8) 2 & F = c h z F = c (ho/2)" 2 A POSTERIORI ERROR ANALYSIS-PART I1 1625 The above results allow us to estimate the effect of a further subdivision by computing ( E : ) ' -- 2 F ; =c(ho/4)" = - ch 0" E F Now, ordering the present error estimates by magnitude, we can refine all the elements which have error estimators over the largest prediction. TERMINATION OF THE ADAPTIVE PROCESS The termination of the adaptive algorithms has to bring under consideration: 1. The achieved accuracy-here we impose that the error norm has to be below the specified tolerance. 2. Storage limitation and time-it is obvious that we cannot achieve an accuracy which is too high. The derivation of the proposed error estimator was based on an eigenvalue analysis and we have adopted an estimate of the form (9) e i 2 = c f I r 2 d n f C , l J 2 d T where c, is a function of the type of problem and type of element we are using, and r and J the domain and boundary residuals, respectively. This error measure is not exact since several approximations are made during its derivation. It does, however, present a good estimation of the error in the solution and we refer to Part I of this paper for results on non-adaptive meshes and the following examples for results in the adaptive case. We also note that adaptive schemes can be used without an error estimate. The process must then be terminated by specifying the maximum number of re-solutions (or degrees-of- freedom) and the accuracy of the final solution decided on an ad hoc basis. OPTIMAL MESHES The problem of optimal meshes has been addressed in the Introduction and it was argued that equilibration of the error is associated with the optimal nodal location. We present a proof of this for the one-dimensional situation in Appendix I. In the two-dimensional section we will show experimentally that this proof is also valid for two-dimensional problems. Having identified the optimal mesh, we add in Appendix I1 a small section on the equivalence of both indications for refinement presented in the previous section, and in Appendix I11 the proof that optimal meshes are stable. This last property implies that for meshes near the optimal we have already a very good error estimate. In practice, special elements are often used to face difficulty with the singular behaviour. This approach is usually effective, although use of only one or few elements without overlap does not solve the accuracy in general. Another question arises of how to judge the quality of the meshes constructed in adaptive mode. We have essentially two possibilities: 1. On a benchmark problem, to compute by the minimization process (very expensive) the best accuracy achievable with given degrees-of-freedom and then compare this best error with 1626 J . P. DE S. R. GAG0 ET AL. the error actually achieved by the adaptive procedure. Although this is theoretically advan- tageous, in practice it is not easy to implement. 2. Check the rates of convergence. For example, using linear elements, then the maximal accuracy for the smooth solution in the energy norm is N-’ /* , N being degrees-of-freedom. For solutions having singularity, as in the corner or crack problem, it can be experimentally shown that with proper refinement the rate N-”’ is achievable. Of course, we also have to count the cost. It seems to be achievable and reasonable to ask that the total cumulative time does not exceed the 2 x time of the solution of the final mesh which gives the required accuracy. The optimal mesh strongly depends on the norm used and the elements being used. Consider the problem -4” = f 4(-1) = 4(0) = 0, f = 1.0 (10) where the exact solution is a parabola and 4’ is a linear function. For linear elements we are dealing with‘the approximation of a linear function by constants. Defining the position of the nodal point by z and using the symmetry we have, from Figure 1, an error in the energy norm in the form [z2/4+(1 - ~ ) ~ / 4 ] ” ~ Figure 1. Linear finite element approximation to quadratic solution (derivatives) Minimizing its square we get 22 -2(1 - z ) = 0 2 = 1/2 So the uniform mesh is optimal and the error indicators will be equal, as they should. Let US now take into consideration Now, considering (see Appendix I) f =fb) (11) that the optimal mesh corresponds to the equilibration of the estimators A POSTERIORI ERROR ANALYSIS-PART I1 1627 For this problem, using linear elements the domain residual becomes and the estimate can be written for an arbitrary element as h2 x i + h 2 E . =- f(x)2dx 1 12 Jxi It is possible to see that having decided on a given number of nodal points this approach allows us to determine the best nodal positions. We are now in a position to comment on the strain energy difference concept that has been frequently adopted in practice, with the objective of obtaining optimal partitions. In this case the exact strain energy is of the form and if the refinement is to be performed in the regions of high strain energy gradients, as adopted for example in T ~ r c k e ~ ~ or S h e ~ h a r d , ~ ~ we refine according to We can conclude that this criterion takes no account of previous refinements and does not provide a method to decide upon the relative distribution of elements to be introduced, but only gives an indication of the area of refinement. We also remark the important fact that the residual is a function of the type of element, implying a variation of the indicator for refinement according to the element type. This does not happen with the exact strain energy variation that remains constant for any element type. The difficulty with the adaptive criterion for the h-p version is well illustrated, because in switching from linear elements to quadratic ones the optimal mesh completely changes the optimal nodal positions. TYPE OF REFINEMENT We must now address the problem of the type of elements being used. In the case of the p-version the hierarchic approach is the best one for many reasons. For the h-version the situation is not so clear, despite the many advantages of the hierarchic approach. The solutions presented here rely on two programs developed for this research which, although self-adaptive, are essentially different from other algorithms such as FEARS14 and COMETX (Szabo et ul.72’74 and Basu et ul.”). The type of elements used in the h-adaptive program is the simple isoparametric four-noded bilinear element with a 2 X 2 Gauss integration rule to determine the respective stiffness matrix. Whenever the adaptive strategy implies a subdivision, as in Figure 2(a), the irregular node (A) is eliminated by introducing in the stiffness matrix of the two irregular elements the information that the final displacement of (A) should be a linear function of the displacement of the neighbouring nodes (B) and (C). The p-adaptive program is based on elements of the type in Figure 2(b), with the Gauss integration rule corresponding to the degree of the highest s