# A posteriori error estimation techniques in practical finite element analysis

Thomas Grätsch, Klaus-Jürgen Bathe

► To cite this version:

Thomas Grätsch, Klaus-Jürgen Bathe. A posteriori error estimation techniques in practical finite element analysis. Computers & Structures, 2005, 83 (4-5), pp.235 - 265. 10.1016/j.compstruc.2004.08.011. hal-01390203

HAL Id: hal-01390203

https://hal.science/hal-01390203

Submitted on 1 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A posteriori error estimation techniques in practical finite element analysis

Thomas Grätsch, Klaus-Jürgen Bathe

Department of Mechanical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Room 3 356, Cambridge, MA 02139, USA

In this paper we review the basic concepts to obtain a posteriori error estimates for the finite element solution of an elliptic linear model problem. We give the basic ideas to establish global error estimates for the energy norm as well as goal oriented error estimates. While we show how these error estimation techniques are employed for our simple model problem, the emphasis of the paper is on assessing whether these procedures are ready for use in practical linear finite element analysis. We conclude that the actually practical error estimation techniques do not provide mathematically proven bounds on the error and need to be used with care. The more accurate estimation procedures also do not provide proven bounds that, in general, can be computed efficiently. We also briefly comment upon the state of error estimates in nonlinear and transient analyses and when mixed methods are used.

Keywords: Finite element analysis; A posteriori error estimation; Goal oriented error estimation; Dual problem; Practical procedures

## 1. Introduction

The modeling of physical phenomena arising in engineering and the sciences leads to partial differential equations in space and time, expressing the mathematical model of the problem to be solved. In general, analytical solutions of these equations do not exist, hence numerical methods such as the finite element method are employed. A major feature of numerical methods is that they involve different sources of numerical errors [1,2]. The focus of this paper is only on the discretization error which is due to the finite element (polynomial) approximation of the solution. Hence, we assume that an appropriate mathematical model has been chosen and, even for this case, we are only concerned with one specific error, namely the discretization error arising in the finite element solution of this model.

Since the late 1970s several strategies have been developed to estimate the discretization errors of finite

element solutions. Basically, there are two types of error estimation procedures available. So called a priori error estimators provide information on the asymptotic behavior of the discretization errors but are not designed to give an actual error estimate for a given mesh. In contrast, a posteriori error estimators employ the finite element solution itself to derive estimates of the actual solution errors. They are also used to steer adaptive schemes where either the mesh is locally refined (h version) or the polynomial degree is raised (p method). Most a posteriori error estimators developed prior to the mid 1990s focused on the global error in the energy norm. Then recently the theory was extended to estimate the error in particular quantities of interest. To understand the importance of this extension it must be realized that many local or global quantities of interest such as deformations, stresses, drag and lift coefficients or the heat transfer of a structural part can be obtained by applying a linear functional to the solution.
This new development is commonly referred to as *goal oriented* error estimation since the aim is to provide error estimates and error bounds for particular quantities of interest.

In the following we will use the words “error estimates” and “error bounds”. An error estimate denotes a quantity that is an approximation to the actual unknown error, whereas upper and lower error bounds are quantities that are always larger resp. smaller than the actual unknown error. Hence, error bounds can be guaranteed but still be inaccurate, whereas an error estimate should be accurate although, in general, it over or underestimates the true error.

Today, a posteriori error estimates are well developed for a large class of simple linear elliptic model problems. The crucial question is whether these procedures are also effective in practical solutions such as in the linear analysis of geometrically complex 2D, 3D and shell problems, and in the analysis of problems including nonlinear effects, time dependent loads or multi physics phenomena. From a practical point of view, there is much interest in reliable and efficient methods to estimate the error in complex analyses. In this context we shall consider reliability to mean that the error estimates can be expected to be accurate, and efficiency to mean that the computer time to obtain these estimates is small when compared to the total processing time used. In reviewing the state of a posteriori error estimators, we will conclude that efficient error estimates still need to be used with care because they are generally not based on guaranteed error bounds. In fact, nearly guaranteed error bounds are still quite expensive to evaluate for complex problems and are frequently not yet available.

The outline of the paper is as follows: We start with a summary of some basic requirements for an error estimator and give a simple schematic example to indicate what an ideal error estimator would give in practical engineering analysis. Next, we consider an elliptic linear model problem for which explicit and implicit error estimators as well as recovery based error estimators for the global energy norm are presented. We then introduce for the same model problem goal oriented error estimators which can be used to estimate the error in the calculation of arbitrary quantities of interest. For both, the global energy and goal oriented error estimators, we give some concluding remarks on the application of these methods in the solution of complex but linear problems. Further more, we briefly comment upon error estimation techniques in nonlinear and time dependent analyses as well as in the use of mixed methods. Then we present some demonstrative numerical solutions. Finally, we conclude with some general suggestions on future developments and open issues.

The theory of error estimation is large and this paper only presents some experiences. Indeed, we only aim to provide the reader with the main ideas of a posteriori

error estimation and a state of the art assessment regarding the use of these procedures when applied to practical problems. For more details on finite element error estimation procedures, and additional views and approaches, see for example Verfürth [3], Estep et al. [4], Ainsworth and Oden [5], Babuška and Strouboulis [6], and the references therein.

## 2. Requirements for an error estimator

The main purpose of any a posteriori error estimator is to provide an estimate and ideally bounds for the solution error in a specified norm or in a functional of interest if the problem data and the finite element solution are available. Some characteristics of an effective error estimator include:

* The error estimate should be accurate in the sense that the predicted error is close to the actual (unknown) error.

* The error estimate should be asymptotically correct in the sense that with increasing mesh density the error estimate should tend to zero at the same rate as the actual error.

* Ideally, the error estimator should yield guaranteed and sharp upper and lower bounds of the actual error.

* The error estimator should be computationally simple, with the error estimate (and bounds) inexpensive to compute when measured on the total computations of the analysis.

* The error estimator should be robust with regard to a wide range of applications, including nonlinear analysis.

* An implementation of the error estimator should be possible to steer an adaptive refinement process with the error estimate used to optimize the mesh with respect to the goal of the computation.

Of course, an ideal error estimator that meets all these requirements is not yet available. But even for linear problems, it is in general not possible to provide inexpensively computable and guaranteed error bounds which would be of much practical interest. To illustrate the requirements, we consider the static analysis of a plate with a hole acted upon by a longitudinal tensile load (see Fig. 1) where one fourth of the plate is shown. Ideally, the application of error estimation procedures should result with little computational cost in an accurate error estimate, based ideally on sharp and guaranteed maximum absolute error bounds, for the dis placements and the stresses at every point of the structure (i.e. for every quantity of interest) which we call the “ideal error estimator solution” in finite element analysis. With the ideal error estimator solution available, the analyst could refine the mesh in those regions
<figure><img src="image_2.png" /><figcaption>Fig. 1. The ideal error estimator solution in finite element analysis: guaranteed maximum absolute error bounds at every point of the structure. Not available yet.</figcaption></figure>

where the error is large in order to obtain an efficient mesh and be certain about the accuracy of the solution for any mesh used.

Of course, a key requirement for this error estimator to be useful in engineering practice is that the computational cost of the error estimate (and ideally the error bounds) must be much smaller than the added computational cost to simply use a very fine mesh.

Today, even when considering only linear analysis, we are far from the ideal error estimator solution, since either error bounds for quantities of interest are guaranteed but expensive to compute or (mostly) not computable at all, or they are computable but not guaranteed. Also, the efficient computation of accurate and nearly guaranteed error bounds achieved in the solution of simple model problems may actually be very expensive for complex engineering problems. Note that global error estimates for the energy norm consider only the error in the global energy norm and do not provide any local information.

In engineering practice, error measures are not yet much used and there are only a few contributions that address the application of effective error estimation procedures in the solution of complex problems of practical interest. We mention these contributions below.

### 3. Model problem

The elliptic linear model problem is Poisson's equation on a bounded, polyhedral and two dimensional domain main $\Omega \in \mathbb{R}^2$ with Lipschitz boundary $\Gamma$ $\Gamma_D \cup \Gamma_N$ where $\Gamma_D$ and $\Gamma_N$ are the Dirichlet and Neumann

boundaries, respectively. The boundary value problem consists of finding the solution $u$ that satisfies

$$ \begin{array}{l} \Delta u = f \quad \text{on } \Omega; \quad u = 0 \quad \text{on } \Gamma_D; \\ \mathbf{n} \cdot \nabla u = g \quad \text{on } \Gamma_N \end{array} \qquad (1) $$

The data are assumed to be sufficiently smooth, that is, $f \in L^2(\Omega)$, $g \in L^2(\Gamma_N)$ and $\mathbf{n}$ is the unit outward normal vector to $\Gamma$. An equivalent formulation of the boundary value problem is the variational formulation seeking $u \in V$ such that

$$ a(u, v) = l(v) \quad \forall v \in V \qquad (2) $$

where the trial and test space $V$ is the usual Sobolev space of functions from $H^1(\Omega)$ whose trace vanishes on the Dirichlet part of the boundary

$$ V = \{v \in H^1(\Omega) : v = 0 \text{ on } \Gamma_D\} \qquad (3) $$

The form $a(u, v)$ is assumed to be a $V$ coercive bilinear form on $V(\Omega) \times V(\Omega)$ and the linear functional $l(v)$ is an element of the dual space $V'(\Omega)$

$$ \begin{aligned} a(u, v) &= \int_{\Omega} \nabla u \cdot \nabla v \, d\Omega \\ l(v) &= \int_{\Omega} f v \, d\Omega + \int_{\Gamma_N} g v \, ds \end{aligned} \qquad (4) $$

Associated with the bilinear form is the energy norm defined by

$$ \|v\|_{E} = \sqrt{a(v, v)} \qquad (5) $$
As is well known, the existence and uniqueness of the variational solution is provided by the Lax Milgram theorem which requires the bilinear form to be bounded and elliptic on $V(\Omega) \times V(\Omega)$

$$|a(v, w)| \le M \|v\|_V \|w\|_V \quad \forall v, w \in V \times V \quad (6)$$

$$a(v, v) \ge \alpha \|v\|_V^2 \quad \forall v \in V \quad (7)$$

where $M$ and $\alpha$ are positive constants independent of $v$ and $w$.

While we consider here a model problem with a scalar for solution, the concepts given in the paper are, of course, also applicable to the general linear elasticity problem by using the appropriate vectors of solution variables and the corresponding solution spaces. Indeed, in the example solutions in Section 9, we only consider more general elasticity problems.

### 3.1. Finite element approximation

The finite element formulation is based on the Bubnov Galerkin procedure where $V_h \subset V$ is assumed to be a finite element subspace consisting of cellwise poly nomial functions of order $p$ over the finite element partition $\mathcal{T}_h$. The partition, or simply mesh, formed by the union of all elements, is assumed to coincide exactly with the domain $\Omega$ and any two elements are either disjoint or share a common edge. Of course, the boundary of each element is also assumed to be Lipschitz continuous. Then the finite element approximation means to find a function $u_h \in V_h$ such that

$$a(u_h, v_h) = l(v_h) \quad (8)$$

for all test functions $v_h \in V_h \subset V$. The error of the finite element approximation denoted by $e_h$ and $u_h$ satisfies the error representation

$$\begin{aligned} a(e_h, v) &= a(u, v) \quad a(u_h, v) \\ &= l(v) \quad a(u_h, v) = \mathcal{R}_h(v) \quad \forall v \in V \end{aligned} \quad (9)$$

which is the basis for a large class of error estimators using the energy norm. Here, $\mathcal{R}_h(\cdot)$ is called the residual functional or the weak residual. If in (9) the choice of test functions is restricted to the finite element space, the fundamental Galerkin orthogonality condition follows [1, 7, 8]

$$a(e_h, v_h) = \mathcal{R}_h(v_h) = 0 \quad \forall v_h \in V_h \quad (10)$$

Assuming that the bilinear form is positive definite it follows that

$$\| \mathcal{R}_h \|_{V'} = \sup_{v \in V(\Omega)} \frac{| \mathcal{R}_h(v) |}{\| v \|_E} = \| e_h \|_E \quad (11)$$

where $\| \mathcal{R}_h \|_{V'}$ denotes the norm of the residual in the dual space $V'(\Omega)$.

### 3.2. A priori error estimates

A priori error estimates provide useful information on the asymptotic behavior of the approximation. The most important property of any conforming finite element formulation based on a symmetric bilinear form is the optimality condition

$$\|u\|_{u_h} = \min_{v_h \in V_h} \|u\|_{V} \quad \forall v_h \in V_h \quad (12)$$

which states that there is no “better” approximation in the finite element space $V_h$ (i.e., that is closer to the exact solution) than the finite element solution itself, if the distance is measured in the energy norm [1]. For both symmetric and nonsymmetric variational problems there is another optimality condition called Céa’s Lemma which asserts that for elliptic bilinear forms it holds

$$\|u\|_{u_h} \le \frac{M}{\alpha} \inf_{v_h \in V_h} \|u\|_{V} \quad (13)$$

where $M$ and $\alpha$ are the constants defined in (6) and (7). Basically, Céa’s Lemma asserts that the error of the finite element solution measured in the $V$ norm is of the same order as the interpolation error. Choosing $V \subset H^1$ and employing interpolation estimates it turns out that the error measured in the $H^1$ norm is of the order $O(h^p)$

$$\|u\|_{H^1(\Omega)} \le c h^p \|u\|_{H^{p+1}(\Omega)} \quad (14)$$

where $c$ is a stability and interpolation constant which does not depend on the actual Ansatz (interpolation) space $^1$ and $h$ denotes the maximum of all element sizes. Furthermore, we have for the error in the $L^2$ norm

$$\|u\|_{L^2(\Omega)} \le c h^{p+1} \|u\|_{H^{p+1}(\Omega)} \quad (15)$$

which means that the convergence rate for the solution itself is $O(h^{p+1})$ [1].

## 4. Global error estimates for the energy norm

In this section, we present various error estimators for the global error in the energy norm of the above mentioned elliptic model problem when a specific (not very fine) mesh has been used. Of course, it is always possible to solve problems (8) or (9) very accurately using a very fine mesh and then the (almost) exact error can be calculated. However, instead, we want to estimate the discretization error while not knowing the exact solution.

In order to evaluate the accuracy of any error estimate, we use the effectivity index defined by

<footnote n="1">¹ Throughout this paper we use $c$ as a generic constant, which thus may take different values in successive occurrences, even in the same equation.</footnote>
$$ \lambda = \frac{E_h}{\|e_h\|_E} \qquad (16) $$

with $E_h$ denoting an estimate for the error in the energy norm. The effectivity index represents the degree of over or underestimation and should be ideally close to 1.0.

Error estimators that are based directly on the finite element approximation and the data of the problem are usually referred to as explicit error estimators. In contrast, implicit error estimators require the solution of auxiliary local boundary value problems. Hence, explicit error estimators in general require less computational effort than implicit schemes, but they involve compro mises in robustness and accuracy. A third class of error estimators are the recovery based error estimators. The main idea of these error estimators is to smooth the gradients of the solution and compare the unsmoothed and the smoothed gradients in order to assess the solution error.

## 4.1. Explicit error estimators

Explicit error estimators involve a direct computation of the interior element residuals and the jumps at the element boundaries to find an estimate for the error in the energy norm, see the fundamental work of Babuška and Rheinboldt [9,11], Babuška and Miller [10] and Kelly et al. [12]. The starting point is the error representation

$$ a(e_h, v) = l(v) \quad a(u_h, v) \quad \forall v \in V \qquad (17) $$

which holds true for arbitrarily chosen test functions $v \in V$. If the domain integral is split into the contributions from each element, (17) can be rewritten for our model problem as

$$ a(e_h, v) = \sum_{K \in \mathcal{F}_h} \left\{ \int_K f v \, d\Omega + \int_{\partial K \cap \Gamma_N} g v \, ds \right. \\ \left. \qquad \int_K \nabla u_h \cdot \nabla v \, d\Omega \right\} \quad \forall v \in V \qquad (18) $$

where $K$ denotes the volume of an element in $\mathcal{F}_h$ and $\partial K$ denotes its boundary. Applying integration by parts to the last term in (18) and rearranging terms leads to

$$ a(e_h, v) = \sum_{K \in \mathcal{F}_h} \int_K R v \, d\Omega + \sum_{\gamma \in \partial \mathcal{F}_h} \int_{\gamma} J v \, ds \quad \forall v \in V \qquad (19) $$

where $R$ is the interior element residual

$$ R = f + \Delta u_h \quad \text{in } K \qquad (20) $$

and $J$ is the jump of the gradient across the element edge $\gamma$

$$ J = \begin{cases} (\mathbf{n} \cdot \nabla u_h + \mathbf{n}' \cdot \nabla u_h') & \text{if } \gamma \not\subseteq \Gamma \\ g & \text{if } \gamma \subseteq \Gamma_N \\ 0 & \text{if } \gamma \subseteq \Gamma_D \end{cases} \qquad (21) $$

where on interelement edges, $\gamma \not\subseteq \Gamma$, the edge $\gamma$ separates elements $K$ and $K'$.

Next, we utilize the Galerkin orthogonality condition (10) to introduce the interpolant $\mathcal{I}_h v$ into (19) which results in

$$ a(e_h, v) = \sum_{K \in \mathcal{F}_h} \int_K R(v \quad \mathcal{I}_h v) \, d\Omega \\ + \sum_{\gamma \in \partial \mathcal{F}_h} \int_{\gamma} J(v \quad \mathcal{I}_h v) \, ds \quad \forall v \in V \qquad (22) $$

Applying the Cauchy Schwarz inequality elementwise yields

$$ a(e_h, v) \le \sum_{K \in \mathcal{F}_h} \|R\|_{L^2(\mathcal{K})} \|v\|_{\mathcal{I}_h v} \|L^2(\mathcal{K})\| \\ + \sum_{\gamma \in \partial \mathcal{F}_h} \|J\|_{L^2(\gamma)} \|v\|_{\mathcal{I}_h v} \|L^2(\gamma)\| \qquad (23) $$

According to results of interpolation theory we have

$$ \|v \quad \mathcal{I}_h v\|_{L^2(\mathcal{K})} \le c h_K \|v\|_{H^1(\bar{K})} \qquad (24) $$

$$ \|v \quad \mathcal{I}_h v\|_{L^2(\partial \mathcal{K})} \le c \sqrt{h_K} \|v\|_{H^1(\bar{K})} \qquad (25) $$

where $h_K$ is the diameter of the element $K$. The symbol $\bar{K}$ denotes the subdomain of elements sharing a common edge with $K$ and $c$ is an unknown interpolation constant which depends for our model problem on the shape of the elements. Using these estimates in (23) leads to

$$ a(e_h, v) \le c \|v\|_{H^1(\Omega)} \left\{ \sum_{K \in \mathcal{F}_h} h_K^2 \|R\|_{L^2(\mathcal{K})}^2 + \sum_{\gamma \in \partial \mathcal{F}_h} h_K \|J\|_{L^2(\gamma)}^2 \right\}^{1/2} \qquad (26) $$

Employing the inequality $\|v\|_{H^1(\Omega)} \le c\|v\|_E$ and substitution $e_h$ in place of $v$ yields the final error bound

$$ \|e_h\|_E^2 \le c \left\{ \sum_{K \in \mathcal{F}_h} h_K^2 \|R\|_{L^2(\mathcal{K})}^2 + \sum_{\gamma \in \partial \mathcal{F}_h} h_K \|J\|_{L^2(\gamma)}^2 \right\} \qquad (27) $$

where apart of the constant $c$ all quantities on the right hand side of (27) can be calculated explicitly. In practice (27) is regrouped as

$$ \|e_h\|_E^2 \le \sum_{K \in \mathcal{F}_h} \left\{ c_1 h_K^2 \|R\|_{L^2(\mathcal{K})}^2 + c_2 h_K \|J\|_{L^2(\partial \mathcal{K})}^2 \right\} \qquad (28) $$

where the constant $c$ is split into contributions $c_1$ and $c_2$ corresponding to the element residual and the jump terms, respectively. We also note that on interelement boundaries $\partial K \not\subseteq \Gamma$ the jump $J$ is multiplied by the factor 1/2 to distribute the error equally onto the two elements sharing the common edge.

The expression in (28) directly leads to a local error indicator $\eta_K$ defined by

$$ \|e_h\|_E^2 \le (E_h)^2 = \sum_{K \in \mathcal{F}_h} \eta_K^2 \\ \text{with } \eta_K^2 = c_1 h_K^2 \|R\|_{L^2(\mathcal{K})}^2 + c_2 h_K \|J\|_{L^2(\partial \mathcal{K})}^2 \qquad (29) $$
Of course, the constants $c_1$ and $c_2$ in (29) are in general unknown. However, research has been conducted to evaluate these constants for specific problems considering more general cases than our model problem but the values obtained are related to worst case scenarios and the error bound is generally not sharp. For instance, Johnson and Hansbo solve in [13] an additional eigenvalue problem in order to obtain computable values for the constants. In case of shell elements, the residual terms in (28) consist of bending, membrane and shear parts so that different weighting factors for each of these parts have to be applied (see e.g. Ref. [14]) and specific care is necessary regarding the “locking” phenomenon [1]. For general linear elastic analysis, the constants are unknown and cannot be calculated with a reasonable computational effort. However, we should point out that using (23) instead of (22) has significant consequences regarding the accuracy of any subsequent error estimate or bound, since the cancellation of errors over the domain is lost and usually results in severe inaccuracy of the error estimate and bound. Therefore, it is hardly an important issue to obtain the best possible constants in (27) (29).

Although the exact constants are not known, the error indicators defined in (29) with approximate constants might be used for driving mesh adaptivity. Given the estimate $E_h$ for the error in the energy norm, the relative error is calculated from

$$e_{\text{rel}} = \frac{\|e_h\|_E}{\|u\|_E} \approx \frac{E_h}{\|u_h\|_E} \qquad (30)$$

One strategy is to enrich the finite element space if the relative error $e_{\text{rel}}$ exceeds a specified tolerance, say $\gamma_{\text{tol}} = 0.01$. The mesh is considered optimal when the contribution of each element to the total error is about the same so that the following ratio indicates if a single element is to be refined or not

$$\xi_K = \frac{m \cdot \eta_K}{\gamma_{\text{tol}} \cdot \|u_h\|_E} = \begin{cases} > 1 & \text{refine element} \\ < 1 & \text{enlarge element} \\ = 1 & \text{no change} \end{cases} \qquad (31)$$

where $m$ is the number of elements used. To avoid over refinement due to overestimation, it is commonly used practice that per refinement step only some percentage of those elements which exceed the ratio $\xi_K > 1$ is refined.

## 4.2. Implicit error estimators

Implicit error estimators involve the solution of auxiliary boundary value problems whose solution yields an approximation to the actual error. The interest in implicit schemes stems from the fact that in explicit schemes the whole information for the total error is obtained only from the given solution, when it might be possible

to obtain more accurate information on the error by solving additional auxiliary problems.

The boundary value problems to be solved are local, which means that they are posed either on a small patch of elements (subdomain residual method) or even only on one single element (element residual method). In general, a drawback of the subdomain residual method can be that solving the local problems is rather expensive, since each element is considered several times. On the other hand, the element residual method needs to approximate the prescribed Neumann boundary data on each single element.

### 4.2.1. Element residual method

In the element residual method (see Bank and Weiser [15]), we define by $e_K = u - u_h$ the local error on a single element $K$ that satisfies for our model problem the variational problem

$$a(e_K, v)_K = \int_K R v \, d\Omega + \int_{\partial K} \left( \frac{\partial u}{\partial n} - \frac{\partial u_h}{\partial n} \right) v \, ds \quad \forall v \in V_K \qquad (32)$$

In (32) the subscript $K$ denotes the restriction of the bilinear form to a single element where the trial and test space is defined by

$$V_K = \{v \in H^1(K) : v = 0 \text{ on } \partial K \cap \Gamma_D\} \qquad (33)$$

In order to impose the correct boundary conditions we have to check whether the boundary $\partial K$ of the single element intersects a portion of the boundary $\Gamma$ of the domain $\Omega$. On the Dirichlet part of the boundary the contribution to the local error is zero. Clearly, on the Neumann part of the global boundary, the true flux in (32) equals the prescribed data $g$. If the element boundary matches an interelement boundary, the true flux, of course, is unknown. Therefore, an approximation to the true flux is taken from the finite element solution

$$\frac{\partial u}{\partial n} \approx \frac{\partial \bar{u}_h}{\partial n} = \frac{1}{2} \mathbf{n} \cdot (\nabla u_h + \nabla u'_h) \qquad (34)$$

where $\nabla u'_h$ denotes the gradient in element $K'$ sharing a common edge with element $K$. Hence, the idea is to aver age the discontinuous finite element flux at interelement boundaries to find an approximation to the true flux. Then the local problem means to find a function $\zeta_K \in V_K$ that satisfies

$$a(\xi_K, v)_K = \int_K R v \, d\Omega + \int_{\partial K} \left( \frac{\partial \bar{u}_h}{\partial n} - \frac{\partial u_h}{\partial n} \right) v \, ds \quad \forall v \in V_K \qquad (35)$$

With the solutions $\xi_K$ of the single elements known, the error is estimated using

$$\|e_h\|_E^2 \approx (E_h)^2 = \sum_{K \in \mathcal{F}_h} \eta_K^2 \quad \text{with} \quad \eta_K^2 = \|\xi_K\|_E^2 \qquad (36)$$
Unfortunately, the existence and uniqueness of the variational problem (35) is not guaranteed due to the possible incompatibility of the prescribed Neumann data. To overcome this drawback several techniques have been proposed. For instance, the problem is reformulated over a subspace $\tilde{V}_K \subset V_K$ for which the local bilinear form $a(\cdot, \cdot)_K$ is elliptic. Another approach employs equal ibrated boundary data so that the local problems remain well posed and therefore the consistency of the estimator is recovered (see e.g. Refs. [5,6]). It should be noted that, in this approach, an upper error bound can easily be obtained using the Cauchy Schwarz inequality. However, the upper estimate of the error is only guaranteed if the local problems are computed exactly. If not, the estimator could underestimate the error [6], but in [16] it is shown that the error arising in the approximation of the local problems can be estimated with explicit estimation schemes.

### 4.2.2. Subdomain residual method

The basic idea of the subdomain residual method is to decompose the global residual Eq. (17) into a number of local problems on small element patches with homo geneous Dirichlet boundary conditions [9,11,17]. Re cently, new versions of the subdomain residual method have been proposed which are more flexible in the choice of the boundary conditions (see Refs. [18 21]). In the approach of Prudhomme et al. [21], the subdomain residual method starts with the fundamental error representation (17) and utilizes the partition of unity property of the shape functions. Consider that in our model problem for a mesh of $n$ nodes, we use $n$ element patches, each node defining as its patch the elements coupling into the node. Let $\varphi_i(x)$ be the shape (or inter polation) function corresponding to node $i$. Then of course

$$ \sum_{i=1}^{n} \varphi_i(x) = 1 \qquad (37) $$

where we consider all boundary conditions removed. Inserting (37) into (17) leads to

$$ a(e_h, v) = \mathcal{R}_h \int v \sum_{i=1}^{n} \varphi_i = \sum_{i=1}^{n} \mathcal{R}_h(v\varphi_i) \quad \forall v \in V \quad (38) $$

We introduce on each patch $\omega_i$ supp $\varphi_i$ the “weighted” bilinear form

$$ a_{\varphi_i}(u, v) = \int_{\omega_i} \varphi_i(\nabla u \cdot \nabla v) \, d\Omega \quad (39) $$

associated with the norm $\|v\|_{\varphi_i} = \sqrt{a_{\varphi_i}(v, v)}$. Then on each patch $\omega_i$ the following local Neumann problems are considered:

Find $\zeta_i \in W$ such that

$$ a_{\varphi_i}(\zeta_i, \psi) = \mathcal{R}_h(\psi \varphi_i) \quad \forall \psi \in W \quad (40) $$

where $W$ is a space of functions on $\omega_i$ defined by

$$ W = \left\{ \begin{array}{l} \psi \in H^1(\omega_i) : \displaystyle\int_{\omega_i} \psi \varphi_i \, d\Omega = 0 \\ \text{and} \\ \displaystyle\int_{\omega_i} |\nabla \psi|^2 \varphi_i \, d\Omega < \infty \end{array} \right\} \qquad (41) $$

for an interior node $i$ and a node $i$ located on the boundary $\Gamma_N$ and

$$ W = \left\{ \begin{array}{l} \psi \in H^1(\omega_i) : \psi = 0 \quad \text{on } \partial\omega_i \cap \Gamma_D \\ \text{and} \\ \displaystyle\int_{\omega_i} |\nabla \psi|^2 \varphi_i \, d\Omega < \infty \end{array} \right\} \qquad (42) $$

for a node $i$ located on the boundary $\Gamma_D$. The reason for this particular choice of the space $W$ is that rigid body modes are this way eliminated, otherwise the solution of some local problems in (40) would only be defined up to a constant. Clearly, after assembling the local problems we have

$$ a(e_h, v) = \sum_{i=1}^{n} \mathcal{R}_h(v\varphi_i) = \sum_{i=1}^{n} a_{\varphi_i}(\zeta_i, v) \quad \forall v \in V \qquad (43) $$

Taking $v$ $e_h$ and using the Cauchy Schwarz inequality yields the guaranteed upper error bound

$$ \|e_h\|_E = \sqrt{a(e_h, e_h)} \le \sqrt{\sum_{i=1}^{n} a_{\varphi_i}(\zeta_i, \zeta_i)} \qquad (44) $$

Of course, the exact solutions $\zeta_i$ of the local problems are unknown. Therefore, we seek approximations $\zeta_i^h$ on each patch yielding the upper error bound [21]

$$ \|e_h\|_E \le \sqrt{\sum_{i=1}^{n} a_{\varphi_i}(\zeta_i^h, \zeta_i^h) + 2 \inf_{v \in V^{p+q}} \|u\|_E \|v\|_E} \quad \forall v \in V^{p+q} \qquad (45) $$

where $p$ is the degree of the polynomials used in the approximation of the original problem and $q$ is the additional degree of the polynomials used for solving the discrete local problems. We note that for $q=0$ this error estimator is not sharp since it doubles the bound which is already provided by the optimality condition. However, the second term on the right hand side of (45) should decrease as the polynomial degree used to solve the local problems increases. In practice we might simply employ

$$ E_h = \sqrt{\sum_{i=1}^{n} a_{\varphi_i}(\zeta_i^h, \zeta_i^h)} \qquad (46) $$

Using the approximations $\zeta_i^h$, it is also possible to obtain a lower bound on the error (see Ref. [21]).

The essential feature of this error estimator is that it localizes both the residual and the bilinear form, and therefore allows to pose Neumann boundary conditions for the local problems on interior patches. The effect of the flux jumps between elements is implicitly taken into
account since the local problems on each patch necessar-
ily include the interelement edges. The error estimator is
free of any constants and does not require flux equilibra-
tion. In [21] the error estimator was tested in some
numerical experiments. However, for practical analysis,
the error estimate established in (46) is, of course, expen-
sive to compute.

4.3. Recovery based error estimators

Recovery based error estimators make use of the fact
that the gradient of the finite element solution is in gen-
eral discontinuous across the interelement boundaries.
Of course, a simple way to visualize the error across
the element boundaries is to use iso bands of stresses
as proposed by Sussman and Bathe [22], see also [1],
which are compared with smoothed values. Here, the
underlying idea is to post process the gradient and to
find an estimate for the true error by comparing the
post processed gradient and the nonpost processed gra-
dient of the approximation. In particular, let M'[u_h] de
note an “improved” approximation to the gradient,
then the a posteriori error estimator is taken to be

$$
(E_h)^2 = \int_{\Omega} | \mathcal{M}[u_h] | \nabla u_h |^2 d\Omega \quad (47)
$$

It is remarkable that this rather heuristic approach gives
surprisingly good results. Under suitable conditions, the
post processed gradient M'[u_h] is labeled superconver-
gent because the finite element non post processed
approximation is closer to the smoothed gradient than
to the exact gradient. This means that for linear interpo-
lation functions we have

$$
d_1 = \| \cdot \mathcal{M}[u_h] \ \nabla u_h \|_{L^2(\Omega)} \le c h^2 \quad (48)
$$

whereas the distance between the exact gradient and the
finite element gradient behaves like

$$
d_2 = \| \nabla u \quad \nabla u_h \|_{L^2(\Omega)} \le c h \quad (49)
$$

But there are situations in which this insight is not of
practical value. Suppose that the finite element gradient
is an inaccurate approximation of the exact gradient,
then the “superconvergence” relation $d_1 < d_2$ does not
help to estimate the true error (see Fig. 2). <fnref n="2" />

A much researched recovery based error estimator
was proposed by Zienkiewicz and Zhu [23], who sug-
gested to post process the discontinuous gradient in
terms of the interpolation functions $\varphi_i \in V_h$ as follows

<figure><img src="image_3.png" /><figcaption>Fig. 2. One dimensional schematic example for the superconvergence property of the recovered gradient.</figcaption></figure>

$$
\nabla u_h^* = \sum_{i=1}^{n} (\nabla u_h^*)_i \varphi_i \quad (50)
$$

Here, the unknown nodal values $(\nabla u_h^*)_i$ are determined
by a standard $L^2$ projection

$$
\int_{\Omega} \varphi_j (\nabla u_h^* \quad \nabla u_h) \, d\Omega = 0 \quad j = 1, \dots, n \quad (51)
$$

This projection results in a linear system for determining
the nodal values

$$
\sum_{i=1}^{n} \int_{\Omega} \varphi_j \varphi_i d\Omega (\nabla u_h^*)_i = \int_{\Omega} \varphi_j \nabla u_h d\Omega \quad j = 1, \dots, n \quad (52)
$$

which is often solved only approximately by diagonaliz-
ing the coefficient matrix. Then, to obtain an error esti-
mator, the improved gradient ∇u_h\* is used instead of the
true gradient

$$
\|e_h\|_E^2 \approx a(e_h^*, e_h^*) = \int_{\Omega} (\nabla u_h^* \quad \nabla u_h)^2 d\Omega \quad (53)
$$

In practical analysis the error estimate is calculated
elementwise

$$
\|e_h\|_E^2 \approx (E_h)^2 = \sum_{K \in \mathcal{F}_h} \eta_K^2 \\
\text{with } \eta_K^2 = \|\nabla u_h^* \nabla u_h\|_{L^2(K)}^2 \quad (54)
$$

This error representation is clearly close to the proce-
dure proposed by Sussman and Bathe [22].

The crucial questions are whether the improved gra-
dients provide a better solution than the originally calcu-
lated gradients and whether the error estimate can be
used to obtain a bound on the true error. Carstensen
and Funken [24] have shown that under certain smooth
ness conditions this error estimate is asymptotically
exact up to higher order terms. Numerical examples
usually show that this error estimator is effective in
smooth problems when using interpolation functions
of degree p 1.

<footnote n="2">² In this context, the attribute “superconvergent” is adopted
from the literature, but it is somewhat misleading since it does
not provide any information about the distance between the
exact solution and the post processed finite element
approximation.</footnote>
Refs. [25,26] show that the approach of Zienkiewicz and Zhu can lead to improved estimates if the interpolation is not using the nodes but so called “superconvergent points” in element patches of the domain. This approach is called the superconvergent patch recovery technique. Actually, in [6] it is shown that sometimes superconvergence can also be obtained without using superconvergent points. Further improvements of recovery based methods can be found in [27] and an analysis of these methods for use on unstructured meshes is presented in [28].

However, of course, the Zienkiewicz Zhu algorithm is not effective in the presence of material discontinuities since in a patch based algorithm these effects are smoothed out. Therefore, Hiller and Bathe [29] propose an element based error estimator that uses higher order accuracy points to recover the strain field. For one dimensional test cases this error estimator is found to be highly accurate compared to the superconvergent patch recovery technique. However, the extension of this approach to 2D and 3D problems is not obvious.

Another weak point of the Zienkiewicz Zhu algorithm is the implicit assumption that oscillations indicate errors and that smooth stresses mean accurate stresses. Hence, the method might break down if the non post processed stresses are smooth. In the following problem [5]

$$EAu''(x) = \mu \sin(2^m \pi x); \quad u(0) = u(1) = 0; \quad m > 0 \qquad (55)$$

a bar is stretched and compressed by an oscillating longitudinal force with an amplitude $\mu$. Suppose this problem is solved on a uniform mesh of $2^n$ linear elements ($n \le m$) with the nodes located at the points

$$x_k = k/2^n; \quad k = 0, 1, \dots, 2^n \qquad (56)$$

Then the piecewise linear finite element solution is zero, because we know that it interpolates the exact solution at the nodes where the exact solution is zero. For example, let $EA = 1$, then we have

$$u(x) = \frac{\mu}{4^m \pi^2} \sin(2^m \pi x) \qquad (57)$$

Hence, the finite element normal forces as well as the recovered normal forces are zero and infinitely smooth and the Zienkiewicz Zhu error estimator predicts zero error, while in truth the exact normal force is far from zero and oscillates as shown in Fig. 3. On the other hand, the error estimator of Hiller and Bathe [29] yields a maximal amplitude of the normal force equal to 15.2789 kN for the choice $N' = 2$ and 12.4528 kN for $N' = 4$ compared to the exact value of 12.5664 kN, see Ref. [29] for the definition of $N'$. Hence, using the estimator of Ref. [29] the error is correctly detected.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Oscillating longitudinal displacements $u(x)$ produce an oscillating normal force $N(x)$ but the finite element solution $u_h(x)$ is zero.</figcaption></figure>

## 4.4. Concluding remarks

In order to assess the presented techniques, we have to realize that there are two different goals using error estimation procedures. The first goal is to estimate the actual error in a suitable norm and ideally provide useful actual error bounds. Note that upper and lower error bounds give of course always an error estimate which is, however, not very useful if the difference in the bounds is too large. An effective error estimate as sesses the true error accurately but may not give a guar anteed error bound. Actually, all the above mentioned error bounds when evaluated in general analyses are not guaranteed and must therefore be used with caution.

The second goal is to steer an adaptive scheme to obtain meshes which are optimal with respect to the aim of the computation. Then while the re meshing or refinement algorithms can be expensive the error estimation need not to be very accurate and inexpensive, explicit schemes with approximate constants can be used to provide criteria for the adaptive mesh refinement. Recovery based error estimators are inexpensive too and also provide quite useful error estimates for this purpose. While there is of course no need to use an actual error bound to establish improved and refined meshes, it is clear that in the final step of mesh and solution acceptance an actual error bound is also desired.

From a practical point of view, recovery based error estimators are efficient compared to other methods and can be generally used. Since all the error estimators do not provide in general guaranteed bounds, the estimation technique most useful is probably the method that works efficiently in general analyses (including nonlinear analyses) and for general finite element discretizations, and provides sufficient or reasonable accuracy in the error estimation.
## 5. Goal-oriented error estimates

In finite element analysis it is frequently the case that the analyst is more interested in certain output data of the finite element approximation than in the global energy norm. In order to find an estimate for the error in the output data pertaining to a specific quantity, or to find at least an effective mesh to accurately solve for this quantity, error estimators for the energy norm are not useful. Hence, more recently so called goal oriented error estimates were developed, which estimate the error in individual quantities of interest using duality tech niques [30 41]. Let $Q(u)$ denote such a (linear) quantity of interest, as for example the mean value of the $x$ com ponent of the gradient in our model problem over a (not necessarily) small patch $\Omega_{\varepsilon}$,

$$Q(u) = |\Omega_{\varepsilon}|^{-1} \int_{\Omega_{\varepsilon}} \nabla_x u \, d\Omega \quad (58)$$

The key for estimating the error in such quantities is the formulation of an auxiliary problem, which is the dual problem to the primal problem actually considered, and which filters out the necessary information for an accurate estimate for the error in the quantity of interest. In the following we recall the basic framework of these procedures and then briefly discuss specific applications using this framework.

### 5.1. The basic framework

In order to develop a general framework, it is convenient to regard the quantity of interest as a bounded, linear functional $Q: V \to \mathbb{R}$ in the dual space $V'$ associated with the norm

$$\|Q\|_{V'} = \sup_{v \in V} \frac{|Q(v)|}{\|v\|_V} \quad (59)$$

Of course, for the bilinear form considered in our model problem, we have a unique $z$ which is the solution of the variational problem

$$a(z, v) = Q(v) \quad \forall v \in V \quad (60)$$

The solution of this problem is referred to as the dual solution and can be interpreted as the generalized Green's function, or the influence function, related to the functional $Q(v)$.

The objective in the following is to find an estimate for the error

$$Q(e_h) = Q(u) \quad Q(u_h) \quad (61)$$

where $Q(u_h)$ denotes the finite element approximation of the quantity of interest and $Q(u)$ is the exact value. Setting $v = e_h$ in (60) gives the exact error representation

$$Q(e_h) = a(z, e_h) = l(z) \quad a(u_h, z) = \mathcal{R}_h(z) \quad (62)$$

Employing the Galerkin orthogonality condition (10) with $v_h$, $z_h$ yields the starting point for goal oriented error estimators

$$Q(e_h) = a(z \quad z_h, e_h) = \mathcal{R}_h(z \quad z_h) \quad (63)$$

In (63) the finite element approximation $z_h$ of the dual solution is obtained by solving the discrete problem

$$a(z_h, v_h) = Q(v_h) \quad \forall v_h \in V_h \subset V \quad (64)$$

Before focusing on the error estimation part in detail, there are three points that we want to mention. Firstly, to solve the discrete dual problem (64), we have to apply equivalent nodal forces $f_i$ defined by the right hand side of (64) and the definition of the functional as in (58). Thus, for (58), for a certain interpolation function $\varphi_i \in V_h$ we have

$$f_i = Q(\varphi_i) = |\Omega_{\varepsilon}|^{-1} \int_{\Omega_{\varepsilon}} \nabla_x \varphi_i \, d\Omega \quad (65)$$

Secondly, using that the finite element spaces for the primal problem and the dual problem are the same, we have

$$\begin{aligned} Q(u_h) &= a(z_h, u_h) = a(u_h, z_h) = l(z_h) \\ &= \int_{\Omega} f z_h \, d\Omega + \int_{\Gamma_N} g z_h \, ds \end{aligned} \quad (66)$$

which means that $Q(u_h)$ can be evaluated by forming the scalar product between the data $f$ and $g$ of the primal problem and the finite element solution $z_h$ of the dual problem [42]. We note that formula (66) expresses a fundamental symmetry condition, since there are always two ways to calculate the quantity of interest $Q(u_h)$:

1. Calculate the finite element solution $u_h$ as usual; post process the gradient by differentiation of the shape functions and perform the necessary integrations.

2. Calculate the influence function $z_h$ and form the scalar product between the data of the primal problem and $z_h$.

Clearly, the more accurate the approximation $z_h$ for the influence function $z$, the more accurate is the calculation of the quantity of interest itself. Note that the well known “unit dummy load” method in structural mechanics is a special case of Eq. (66). Furthermore, the use of (66) in many mixed finite element methods is direct since (66) only uses the symmetry of the bilinear form.

Thirdly, it should be pointed out that the dual approach is very general in that even global quantities like the mean value of the solution can be defined as a quantity of interest

$$Q(u) = |\Omega|^{-1} \int_{\Omega} u \, d\Omega \quad (67)$$
which yields the following variational problem seeking the dual solution

$$a(z, v) = Q(v) = |\Omega|^{-1}(1, v) \quad \forall v \in V \qquad (68)$$

Hence, the solution of the dual problem equals the solution of the primal problem for a constant domain load $f$ 1 up to the factor $|\Omega|^{-1}$.

## 5.2. Specific applications using the basic framework

In the following, we discuss briefly specific applications using the basic framework presented in Section 5.1 in order to estimate the error in a quantity of interest.

### 5.2.1. Energy norm based estimates

There are several strategies for goal oriented error estimation based on energy norm estimates of the primal problem and the dual problem. For instance, using the error representation (63) and employing the Cauchy Schwarz inequality, the following upper error bound is obtained

$$|Q(e_h)| = |a(z \ z_h, u \ u_h)| \le \|z \ z_h\|_E \|u \ u_h\|_E \quad (69)$$

Hence, the error in the quantity of interest is bounded by the error in the energy norm of the primal problem weighted with the error in the energy norm of the dual problem. It follows that any of the error estimators derived in Section 4 may be used to estimate the error in the evaluation of the quantity of interest. Moreover, the estimate (69) shows that the convergence rate of $|Q(e_h)|$ is larger than the convergence rate of the error in the energy norm, which means that any linear functional acting on the finite element space is superconvergent compared to the convergence in the energy norm. Of course, the actual convergence rate of the functional depends on the regularity of the solutions of the primal problem and the dual problem.

In practice, the Cauchy Schwarz inequality is applied elementwise in (69), which leads to

$$|Q(e_h)| \le \sum_{K \in \mathcal{F}_h} \{\|z \ z_h\|_{E(K)} \|u \ u_h\|_{E(K)}\} \quad (70)$$

However, goal oriented error estimators based on energy norm estimates as in (69) or (70) eliminate the can cellation of errors over the domain, and therefore produce significant error overestimation in general. In [43] it is shown that this overestimation gets even worse by increasing the polynomial degree of the interpolation functions used in a p version of the finite element method.

Different techniques with the aim to provide tighter error bounds have been proposed. For instance, Prudhomme and Oden [44,45] employ the parallelogram identity to problems with symmetric bilinear forms leading to upper and lower bounds on the error. In this ap

proach, using $\epsilon_h$ z $z_h$, the starting point is the error representation (see also [46])

$$Q(e_h) = a(e_h, \epsilon_h) = \frac{1}{4} \left\| se_h + \frac{\epsilon_h}{s} \right\|_E^2 \quad \frac{1}{4} \left\| se_h \frac{\epsilon_h}{s} \right\|_E^2 \quad (71)$$

where $s \in \mathbb{R}$ is a scaling factor chosen to be

$$s = \sqrt{\frac{\|\epsilon_h\|_E}{\|e_h\|_E}} \qquad (72)$$

Let $\eta_{\text{low}}^{+}$, $\eta_{\text{upp}}^{+}$, $\eta_{\text{low}}^{-}$, $\eta_{\text{upp}}^{-}$ denote error estimators that satisfy

$$\eta_{\text{low}}^{+} \le \|se_h + \frac{\epsilon_h}{s}\|_{E} \le \eta_{\text{upp}}^{+} \qquad (73)$$

$$\eta_{\text{low}}^{-} \le \|se_h - \frac{\epsilon_h}{s}\|_{E} \le \eta_{\text{upp}}^{-} \qquad (74)$$

then (71) implies

$$\frac{1}{4} (\eta_{\text{low}}^{+})^2 \frac{1}{4} (\eta_{\text{upp}}^{-})^2 \le Q(e_h) \\ \le \frac{1}{4} (\eta_{\text{upp}}^{+})^2 \frac{1}{4} (\eta_{\text{low}}^{-})^2 \qquad (75)$$

which provides lower and upper bounds for the error $Q(e_h)$. In [45] these error bounds have numerically been found to be accurate for simple model problems when the error estimators defined in (73) and (74) are calculated using implicit methods.

In this context, Stein and Ohnimus [47] note that by neglecting the lower error bound $\eta_{\text{low}}^{-}$, the upper error bound in (75) reduces to (69), which demonstrates that (75) gives a better error bound than (69). Hence, the lower bound error estimation is essential to obtain an efficient upper error bound for goal oriented error estimation. Of course, (75) is also providing better error bounds than

$$\frac{1}{4} (\eta_{\text{upp}}^{-})^2 \le Q(e_h) \le \frac{1}{4} (\eta_{\text{upp}}^{+})^2 \qquad (76)$$

but these bounds are less expensive to evaluate [37,44].

To steer an adaptive refinement process using energy norm based estimates, we can use the upper error bound defined in (70) in combination with explicit error estimates for the energy norm error yielding

$$|Q(e_h)| \le \eta_h = \sum_{K \in \mathcal{F}_h} \omega_K \eta_K \qquad (77)$$

where the local contributions to the error can be calculated from

$$\omega_K = \left( c_1 h_K^2 \|R_z\|_{L^2(\hat{K})}^2 + c_2 h_K \|J_z\|_{L^2(\hat{K})}^2 \right)^{1/2} \qquad (78)$$

$$\eta_K = \left( c_1 h_K^2 \|R\|_{L^2(\hat{K})}^2 + c_2 h_K \|J\|_{L^2(\hat{K})}^2 \right)^{1/2} \qquad (79)$$

In (78) and (79) $c_1, c_2$ are the unknown constants mentioned already with (29) and the $R_z$ and $J_z$ denote the
element residuals and jumps of the gradient correspond
ing to the dual problem. Given the estimate for the
quantity of interest based on (77), or any other goal ori-
ented error estimator, the relative error in the goal quan
tity is calculated from

$$
e_{\text{rel}} = \frac{|Q(e_h)|}{|Q(u)|} \approx \frac{\eta_h}{|Q(u_h)|} \qquad (80)
$$

As in the case of error estimates for the energy norm, one strategy is to enrich the finite element space if the relative error $e_{\text{rel}}$ exceeds a specified tolerance $\gamma_{\text{tol}}$. For equal distribution of errors across elements, the following ratio is used to indicate whether a single element is to be refined or not

$$
\xi_K = \frac{m \cdot (\omega_K \eta_K)}{\gamma_{\text{tol}} \cdot |Q(u_h)|} = \begin{cases} > 1 & \text{refine element} \\ < 1 & \text{enlarge element} \\ = 1 & \text{no change} \end{cases} \quad (81)
$$

Since we are considering the error bound in (70) and hence overrefinement due to overestimation is still an is sue, in practice, it may be best if per refinement step only some percentage of those elements which exceed the ratio $\xi_K > 1$ is refined. In order to assess the accuracy of the error estimate for the quantity of interest, we use the effectivity index

$$
\lambda = \frac{\eta_h}{|Q(e_h)|} \qquad (82)
$$

5.2.2. The dual weighted residual method

Rannacher and co workers (see [48] and the refer
ences therein), proposed the dual weighted residual
method, which is based on an exact representation of
the error in the quantity of interest as follows

$$
\begin{align}
Q(e_h) &= a(e_h, z \quad v_h) \nonumber \\
&= \sum_{K \in \mathcal{F}_h} \left\{ \int_K R(z \quad v_h) \, d\Omega + \int_{\partial K} J(z \quad v_h) \, ds \right\} \tag{83}
\end{align}
$$

where R and J are the element residuals and jumps of the
primal problem defined in (20) and (21), respectively <fnref n="3" />
and $v_h \in V_h$. To drive an adaptive mesh refinement
process, the following element error indicators can be
used

$$
\eta_K = |(R, z \quad v_h)_K + (J, z \quad v_h)_{\partial K}| \quad (84)
$$

which yield the guaranteed upper error bound

$$
|Q(e_h)| \le \sum_{K \in \mathcal{F}_h} \eta_K \qquad (85)
$$

There are several strategies for evaluating the unknown
dual solution in (83). One possibility is to solve the dual
problem by using a higher order method where, for
instance, biquadratic functions are used instead of
bilinear functions. This yields the approximate error
representation

$$
Q(e_h) \approx \sum_{K \in \mathcal{F}_h} \left\{ \int_K R(z_h^{(2)} \quad \mathcal{I}_h z_h^{(2)}) \, d\Omega + \int_{\partial K} J(z_h^{(2)} \quad \mathcal{I}_h z_h^{(2)}) \, ds \right\} \quad (86)
$$

where the upper index (2) denotes the solution obtained
with a higher order method and $\mathcal{I}_h$ is the interpolant on
$V_h$. Instead of using this rather expensive approach, an
other technique consists of approximating the dual solu
tion by employing higher order interpolation functions
defined over single patches of the domain, resulting in
the error representation

$$
Q(e_h) \approx \sum_{K \in \mathcal{F}_h} \left\{ \int_K R(\mathcal{I}_{2h} z_h \quad z_h) \, d\Omega + \int_{\partial K} J(\mathcal{I}_{2h} z_h \quad z_h) \, ds \right\} \quad (87)
$$

where $\mathcal{I}_{2h}$ denotes the interpolant obtained with higher order interpolation functions. Both strategies lead to useful error indicators and can even lead to effectivity indices close to 1.0 depending on the problem consid- ered. However, error estimates using (86) or (87) may underestimate the actual error since the dual problem is only solved approximately. In contrast, error bounds based on (85) may significantly overestimate the true error, since all the error indicators contribute positive terms to the estimate.

Note that instead of applying integration by parts to
obtain the element residuals and jump terms in (83), we
could also leave the residual in weak form yielding for
our model problem the error representation

$$
Q(e_h) = \sum_{K \in \mathcal{F}_h} \left\{ \int_K l(z \quad v_h) \, d\Omega \quad \int_K \nabla u_h \cdot \nabla (z \quad v_h) \, d\Omega \right\} \quad (88)
$$

Compared to (83), error estimates based on (88) are
accurate as well and the derivation of local error indica
tors according to (84) is straightforward.

The dual weighted residual method has been success
fully applied in the solution of various problem, however
the extension to nonstandard problems such as fluid
structure interactions is still to be achieved.

5.2.3. Direct use of influence functions

Instead of calculating element based residuals, we
can also directly employ influence functions and (66).
This approach was used by Grätsch and Bathe in shell
analyses [49]. Using a continuum mechanics shell formu-
lation [1,50], the error in a linear quantity of interest can
be written as

<footnote n="3">³ In (83) we assume that on interelement boundaries $\partial K \not\subseteq \Gamma$
the jumps $J$ are multiplied by $1/2$ to distribute the error equally
onto the two elements sharing the common edge.</footnote>
$$Q(\mathbf{E}_h) = \int_{\Omega_{\text{ref}}} \mathbf{F} \cdot \mathbf{Z}_{\text{ref}} \, dV \quad \int_{\Omega_h} \mathbf{F} \cdot \mathbf{Z}_h \, dV \qquad (89)$$

where $\mathbf{E}_h$ denotes the error in the shell finite element approximation of the primal problem, $\mathbf{Z}_{\text{ref}}$ is the exact influence function corresponding to a 3D reference do main $\Omega_{\text{ref}}$, $\mathbf{Z}_h$ denotes the finite element solution for $\mathbf{Z}_{\text{ref}}$ obtained on a 3D computational domain $\Omega_h$, and $\mathbf{F}$ are body loads. Assuming $\mathbf{Z}_h$ is a MITC4 shell finite element solution [51,52], a useful strategy is to replace the un known dual solution with an approximation $\mathbf{Z}_h^{(\text{MITC9})}$ corresponding to a MITC9 shell finite element formula tion such that

$$Q(\mathbf{E}_h) \approx \int_{\Omega_{h0}} \mathbf{F} \cdot \mathbf{Z}_h^{(\text{MITC9})} \, dV \quad \int_{\Omega_h} \mathbf{F} \cdot \mathbf{Z}_h \, dV \qquad (90)$$

where $\Omega_{h0}$ denotes the 3D reference mesh using MITC9 elements. In [49] it is found in some numerical studies that the error representation (90) leads to accurate error estimates and to effectivity indices close to 1.0 (since (90) takes into account the cancellation effect of the error over the domain). The error estimate includes the error due to the approximation in the geometry of the shell structure.

For the purpose of driving mesh adaptivity we can also obtain local error indicators by considering

$$Q(\mathbf{E}_h) \approx \sum_{K \in \mathcal{F}_h} \left\{ \int_K \mathbf{F} \cdot \mathbf{Z}_h^{(\text{MITC9})} \, dV \quad A_h(\mathbf{Z}_h, \mathbf{U}_h)_K \right\} \qquad (91)$$

in which $A_h(\cdot, \cdot)$ denotes the bilinear form of the MITC4 shell variational problem and $\mathbf{U}_h$ is the finite element solution. An example using the error representation (90) is given in Section 9.1.

### 5.2.4. The Green's function decomposition method

An important step in the use of goal oriented error estimates is to seek an accurate solution of the dual problem. Considering individual output data such as point values, the presented procedures face the problem that in $H^1(\Omega)$ functions are not continuous (i.e. have meaningful point values) in two or three dimensions. Although an approximation of the Green's function can still be obtained, this problem can be circumvented either by using regularized functionals, as proposed for example in [44] or [48], or, when available, using analyt ical functions (see Kelly et al. [53,54]) or Grätsch and Hartmann [55]. The idea is to approximate the Green's function by splitting it into a regular part and a well known fundamental solution, and for this reason the method can be referred to as the Green's function decomposition method. Hence, the method can only be used for the rather restricted case where the analytical fundamental solution is available. Fundamental solu tions in finite element discretizations have of course been used in various applications, and notably in fracture mechanics.

Consider our model problem (see Section 3) and let $\mathbf{u}(\mathbf{x})$ denote the point quantity of interest at the point $\mathbf{x} \in \Omega$. Then using the corresponding Green's function and the problem data we have

$$\mathbf{u}(\mathbf{x}) = \int_{\Omega} G_0(\mathbf{y}, \mathbf{x}) f(\mathbf{y}) \, d\Omega + \int_{\Gamma_N} G_0(\mathbf{y}, \mathbf{x}) g(\mathbf{y}) \, ds \qquad (92)$$

where $G_0(\mathbf{y}, \mathbf{x})$ denotes the Green's function, which in this case is the displacement at the point $\mathbf{x}$ due to a point load at the point $\mathbf{y}$. The next step is to decompose the Green's function into a regular part $u_R(\mathbf{y}, \mathbf{x})$ and the corresponding fundamental solution $F_0(\mathbf{y}, \mathbf{x})$

$$G_0(\mathbf{y}, \mathbf{x}) = u_R(\mathbf{y}, \mathbf{x}) + F_0(\mathbf{y}, \mathbf{x}) \qquad (93)$$

where the fundamental solution is given by

$$F_0(\mathbf{y}, \mathbf{x}) = \frac{1}{2\pi} \ln r \quad \text{with} \quad r = |\mathbf{y} - \mathbf{x}| \qquad (94)$$

The regular part is the solution of the following bound ary value problem

$$\begin{gathered} \Delta u_R = 0 \text{ in } \Omega; \quad u_R = F_0 \text{ on } \Gamma_D; \\ g_R = \mathbf{n} \cdot \nabla F_0 \text{ on } \Gamma_N \end{gathered} \qquad (95)$$

Hence the boundary conditions to obtain $u_R$ are given by (the negative of) the fundamental solution on the boundary. The finite element approximation $u_R^h$ to $u_R$ is obtained by solving:

Find $u_R^h \in V_h$ such that

$$a(u_R^h, v_h) = l(v_h) \quad \forall v_h \in V_{0,h} \qquad (96)$$

$$u_R^h = F_0 \quad \text{on } \Gamma_D \qquad (97)$$

where $V_{0,h} = H_0^1(\Omega) \cap V_h(\Omega)$ and the linear functional is defined by

$$l(v_h) = \int_{\Gamma_N} g_R v_h \, ds \qquad (98)$$

Denoting the so constructed Green's function by $G_{0,h}(\mathbf{y}, \mathbf{x}) = u_R^h(\mathbf{y}, \mathbf{x}) + F_0(\mathbf{y}, \mathbf{x})$, the integral representa tion of the improved displacement becomes

$$u_h(\mathbf{x}) = \int_{\Omega} G_{0,h}(\mathbf{y}, \mathbf{x}) f(\mathbf{y}) \, d\Omega + \int_{\Gamma_N} G_{0,h}(\mathbf{y}, \mathbf{x}) g(\mathbf{y}) \, ds \qquad (99)$$

To derive an a posteriori upper error bound for the im proved quantity we proceed as follows. Let $e_h(\mathbf{x})$ $u_h(\mathbf{x})$ $u_h(\mathbf{x})$ denote the local error. Then we have by subtracting (99) from (92)

$$e_h(\mathbf{x}) = \int_{\Omega} (G_0 - G_{0,h}) f \, d\Omega + \int_{\Gamma_N} (G_0 - G_{0,h}) g \, ds \qquad (100)$$

Since the exact Green's function $G_0$ and the approxi mated Green's function $G_{0,h}$ both contain the same fun damental solution we have as well
$$ 
\begin{aligned}
e_h(\mathbf{x}) &= \int_{\Omega} (u_R \quad u_R^h) f \, d\Omega + \int_{\Gamma_N} (u_R \quad u_R^h) g \, ds \\
&= a(u_R \quad u_R^h, u)
\end{aligned}
\quad (101) $$

Using the Galerkin orthogonality to introduce the finite element solution and applying the Cauchy Schwarz inequality as usual yields the error bound

$$ | e_h(\mathbf{x}) | \le \| u_R \quad u_R^h \|_E \| u \quad u_h \|_E \quad (102) $$

which states that the pointwise error is bounded by the error in the energy norm of the primal problem weighted with the error in the energy norm of the corresponding regular part problem. Comparing this result with (69), we realize that, to use (102), we can employ the usual error estimation techniques of goal oriented error esti mates since the error bounds in (69) and (102) only differ by the weighting factors. Clearly, the method is also clo sely related to the dual weighted residual technique. The only difference is of course that in the Green’s function decomposition method we actually construct an approx imation of the Green’s function for evaluating the local quantity with this function, while in other goal oriented error estimates we simply differentiate the shape func tions. Of course, the method can also be used for point wise stresses; the only difference is that appropriate fundamental solutions need to be employed, which we can expect to yield high solution accuracy.

Indeed, in some numerical experiments the Green’s function decomposition method gave excellent results in that high convergence rates in the required solutions were observed (see Ref. [55], and also Section 9.3). But, as mentioned already, the method is rather restric tive in that the fundamental solution needs to be avail able (see (97) and (98)). Additional considerations arise also, for example, if the load is not smooth, the material is not homogeneous [56], or the quantity of interest is a point value on the boundary [53].

### 5.2.5. Exact bounds approach

The “exact bounds approach” proposed by Peraire and co workers [57] can be used to obtain guaranteed upper and lower bounds for quantities of interest. The basic idea in this approach is to use a displacement based finite element method to obtain a guaranteed lower bound on the exact strain energy and to employ a hybrid (stress based) finite element method using com plementary energy principles to obtain a guaranteed upper bound. Recasting the problem in terms of the dual problem corresponding to the functional of interest, guaranteed bounds can be obtained for the quantity of interest using the solutions of local problems.

In contrast to other techniques presented in this paper, this approach is fundamentally different since it provides results to actually bound the true error in the finite element solution. However, from a practical point of view it is noted that in this approach the primal prob

lem and the dual problem have to be solved two times each, namely for the displacement method and the stress method. In addition, local problems have to be analyzed to obtain the error bounds, which could make the proce dure costly in actual applications. Also, for some types of problems, notably (almost) incompressible linear and nonlinear problems and shell problems, mixed (or hybrid) formulations must be used to circumvent the “locking” phenomenon [1], and any approach to esti mate errors needs to take this fact into account. So far the method has been presented for coercive linear model problems including (compressible) elasticity and advec tion diffusion reaction problems [58].

### 5.2.6. Reduced basis output bounds approach

We also want to mention some achievements which do not exactly fit into the basic framework of Section 5.1 but which also provide estimates for quantities of interest.

The “reduced basis output bounds approach” pro posed by Patera and co workers [59] addresses the prac tical case in which the quantity of interest needs to be computed for a certain number of parameters that are the values of specific variables describing some input data such as physical properties or geometry data.

The underlying idea of this approach is to solve the problem for a sample of the parameters on a very fine finite element mesh. Then, having these solutions, the solution corresponding to any other configuration with in the parameter set is obtained by some interpolation technique. In other words, the unknown solution is ex pressed as a linear combination of the basis functions of a reduced basis approximation space, which are the finite element solutions of the sample problems.

In practice, we generally do not know how to choose the size of the parameter space and it certainly depends on the desired accuracy, the selected quantities of inter est, and the particular problem analyzed. In order to use a parameter space that is neither inefficient nor unac ceptably uncertain, a posteriori error bounds for the quantity of interest are computed which, e.g., are based on the use of a richer approximation space regarding the number of parameters chosen in the sample. Clearly, the whole strategy of this approach relies on the fact that the solutions of the parametrized problems do not vary randomly with the parameters, but in fact vary in a smooth fashion.

## 5.3. Global or pollution error versus local error

In finite element analysis we may ask whether it is sufficient to just use a graded and reasonably refined mesh only around the point of interest. We will address this question somewhat in the context of the global error which has also been labeled the pollution error [6].
The global (or pollution) error can be defined as follows: Let $\mathcal{I}_h u$ denote the interpolant of the exact solution, then we can write

$$e_h = (u - \mathcal{I}_h u) + (\mathcal{I}_h u - u_h) = e_{h,\text{loc}} + e_{h,\text{glob}} \quad (103)$$

where $e_{h,\text{loc}}$ and $e_{h,\text{glob}}$ denote, respectively, the local error and the global error. Of course, in trivial 1D problems the global error is zero if the finite element solution interpolates the exact solution at the nodes. However, in general, considering a structural analysis (elliptic) problem, the geometry, the loads, and boundary conditions have a global influence on the solution, and $e_{h,\text{glob}}$ can not be expected to be zero. Indeed, we need to expect that this part of the error is large when there are stress singularities or high stress gradients in the solution that are not sufficiently resolved by the finite element solution.

To illustrate the global error numerically, we consider a bar of length $l$ with varying cross section governed by the equation

$$(EA(x)u'(x))' = f(x) \quad \text{for } 0 < x < l, \quad u(0) = u(l) = 0 \quad (104)$$

This problem is similar to the problem studied by Baška and Strouboulis [6], but we give further results. Here, we choose for the cross sectional area

$$A(x) = A_1 + x \frac{A_2}{l} \quad A_1 = t \cdot A_2 \quad (105)$$

to simulate a singular solution at $x=0$ if $t$ approaches zero, so that the solution is expected to be polluted. The exact solution of this model problem for different values of $t$ and chosen data is shown in Fig. 4. For our numerical studies we choose $t=0.01$, which, with the data in Fig. 4, yields the exact solution

$$u(x) = 1.01x + 0.219341 \ln(1 + 99x) \quad (106)$$

For the approximation we use eight linear finite elements on a uniform mesh with exact integration of the stiffness matrix for each element of varying cross sectional area.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Exact solutions for the bar with varying cross section $A_1 = t \cdot A_2$ ($E=1, l=1, A_2=1, f=1$).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Global and local errors for $t=0.01$ ($E=1, l=1, A_2=1, f=1$). Eight element solution.</figcaption></figure>

The results in Fig. 5 show that only in the first element the local error is larger than the global error and else where the local error becomes negligible; thus, the global error dominates in all other elements. Finally, in Fig. 6 we give only the global error with respect to the factor $t$. As expected, the global error decreases if the singularity in the exact solution decreases. Of course, for values of about $t<0.1$ the problem solved is not of much practical relevance.

To exemplify the global error in a 2D problem numerically we consider a plate in plane stress conditions (see Fig. 7). The mesh consists of bilinear 4 node elements. We compare the results for the given mesh with respect to a very fine finite element solution representing an almost exact reference solution. Fig. 8 shows the absolute values of the global error at different $y$ levels. Using these absolute values to evaluate percentage values, the maximum percent error in the $u_x$ displacement is ~6% at the point A ($x=4.5/y=4.5$) and ~24% in the $u_y$ direction, and the percentage error decreases considerably similar as in the one dimensional problem considered in Fig. 5 with increasing distance from the fixed corner

<figure><img src="image_7.png" /><figcaption>Fig. 6. Global error with varying cross section $A_1 = t \cdot A_2$ ($E=1, l=1, A_2=1, f=1$). Eight element solution.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Plate in plane stress condition to study the global error.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Global error for the plate in plane stress condition: (a) error in $u_x$ displacement, (b) error in $u_y$ displacement.</figcaption></figure>

towards the point B. At the point B the percentage error is ~2% in both displacements, $u_x$ and $u_y$. Also, the global error would be smaller if the corner were not sharp as in a good practical design. For further details on this part of the error in 2D analysis see Ref. [6].

In practice, frequently the errors can successfully be controlled by refining the mesh in the neighborhood of the high stress gradients. Hence, to obtain an accurate approximation to the quantity of interest, we need to refine the mesh mostly around the domain of interest and at singularities. Indeed, in practice, the areas of high stress gradients are usually also the areas of interest and hence, natural mesh refinements in those areas frequently control the global error. By employing goal oriented error estimation procedures the proper balance between the local error and the global error is always maintained and the mesh is refined to be most beneficial

to the accuracy of the quantity of interest. Of course, this approach has the disadvantage that the mesh refinement is specific to only the one chosen goal quantity, whereas, in practice, the stresses over the complete analysis domain are frequently required.

### 5.4. Concluding remarks

Goal oriented error estimation techniques focus on establishing accurate and computationally inexpensive error estimates for a quantity of interest. We conclude that any error estimator for the energy norm may be applied to goal oriented error estimation procedures if (69) is used to bound the error and the experiences regarding these error estimators are applicable (see Section 4).
Considering the dual weighted residual method, the error estimate may over or underestimate the actual error if the approximate error representations (86) or (87) are used. Error bounds as in (70) and (85) lead gen- erally to an overestimation of the true error, since the cancellation effect of the error is neglected [60]. More over, the upper bound (85) is, strictly, not guaranteed if the unknown dual solution is approximated. However, in practice, the main issue of this error bound is the overestimation rather than providing a not guaranteed upper bound estimate.

When considering complex problems of practical interest, we can expect that the error bounds in (75) may lead to a rather large distance between the lower and upper bounds. In this regard, some results for 3D problems were recently published [61] but additional experiences need be obtained.

6. Nonlinear problems

In solving nonlinear problems by a Galerkin finite element method, we basically solve a series of linear problems in using an iteration scheme such as the New- ton Raphson method. Any of the error estimators derived in Sections 4 and 5 may therefore also in prin- ciple be employed to estimate the error in nonlinear problems. In the following we explain the basic ideas for a simple model problem and emphasize goal ori- ented error estimates.

As a nonlinear model problem we choose a modifica-
tion of the earlier model considered (see Section 3),
namely the diffusion reaction equation <fnref n="4" />

$$
\Delta u \quad u^3 = f \quad \text{on } \Omega; \quad u = 0 \quad \text{on } \Gamma \tag{107}
$$

The standard variational form is to find $u \in V$ such that

$$
a(u; \psi) = l(\psi) \quad \forall \psi \in V \quad (108)
$$

where *u* is the exact (weak) solution, *a*(:,.) is a semi lin-
ear form that is linear only in the second argument, *l*(.) is
a linear functional and *V* is the appropriate Sobolev
space. The semi linear form in our model case is defined
by

$$
a(u; \psi) = (\nabla u, \nabla \psi) \quad (u^3, \psi) \quad (109)
$$

The corresponding finite element approximation uses fi-
nite dimensional subspaces $V_h \subset V$ to determine $u_h \in V_h$
such that

$$
a(u_h; \psi_h) = l(\psi_h) \quad \forall \psi_h \in V_h \quad (110)
$$

assuming that *u*<sub>*h*</sub> is the exact solution of (110). The error
*e*<sub>*h*</sub> u *u*<sub>*h*</sub> of the finite element approximation satisfies
the error representation

$$
\mathcal{R}_h(\psi) = l(\psi) \quad a(u_h; \psi) \quad \forall \psi \in V \quad (111)
$$

from which we obtain the Galerkin orthogonality

$$
\mathcal{R}_h(\psi_h) = l(\psi_h) \quad a(u_h; \psi_h) = 0 \quad \forall \psi_h \in V_h \quad (112)
$$

In order to linearize the error representation (111) with
respect to e<sub>h</sub>, we consider

$$
\begin{align}
\mathcal{R}_h(\psi) &= a(u; \psi) \quad a(u_h; \psi) \tag{113} \\
&= (\nabla u, \nabla \psi) \quad (u^3, \psi) \quad (\nabla u_h, \nabla \psi) + (u_h^3, \psi) \tag{114} \\
&= (\nabla e_h, \nabla \psi) \quad (u^3 \quad u_h^3, \psi) \tag{115} \\
&= (\nabla e_h, \nabla \psi) \quad (e_h^3 + 3u_h^2 e_h + 3u_h e_h^2, \psi) \tag{116}
\end{align}
$$

and provided that the finite element solution $u_h$ is suffi-
ciently close to the exact solution $u$, we obtain

$$
\mathcal{R}_h(\psi) \approx (\nabla e_h, \nabla \psi) \quad (3u_h^2 e_h, \psi) \quad (117)
$$

This error representation corresponds to the assumption
that the error $e_h$ $u$ $u_h$ can be estimated by the solu-
tion of the linearized problem utilizing the displacement
$w$ as an approximation to $u$ that satisfies

$$
a'(u_h; w - u_h, \psi) = \sum_{K \in \mathcal{F}_h} \left\{ \int_K R \psi \, d\Omega + \int_{\partial K} J \psi \, ds \right\} \quad \forall \psi \in V \quad (118)
$$

to which we can apply any of the error estimators pre-
sented in Section 4. In (118) we use the Gateaux deriva-
tive of the semi linear form with respect to the first
argument defined by

$$
\begin{align}
a'(u; \varphi, \psi) &= \left[ \frac{d}{d\varepsilon} a(u + \varepsilon\varphi; \psi) \right]_{\varepsilon=0} \nonumber \\
&= (\nabla\varphi, \nabla\psi) \quad (3u^2\varphi, \psi) \tag{119}
\end{align}
$$

which is, clearly, the usual tangent form used in stan-
dard Newton iteration schemes [1].

To estimate the error in an individual (not necessarily
linear) quantity of interest we distinguish between the
error in the quantity

$$
Q(u) \quad Q(u_h) = a(z; u) \quad a(z_h; u_h) \tag{120}
$$

and the quantity of the error

$$
Q(e_h) = a(z; e_h) = a(z; u) \quad a(z; u_h) \quad (121)
$$

But the solution *z* in (120) and (121) cannot be inter-
preted as an influence function in the original meaning
of linear analyses. In linear analysis, of course, (120)
and (121) express the same relationships and hence it
is prudent to linearize in order to estimate the error in
the quantity of interest.

We therefore consider for a specified *u* the dual prob-
lem: calculate *z* ∈ *V* such that

$$
a'(u; z, \psi) = Q'(u; \psi) \quad \forall \psi \in V \quad (122)
$$

<footnote n="4">⁴ We employ this problem since there are no difficulties in proving the uniqueness and stability of the solution (see Ref. [62]), hence we can straightforwardly focus on the error estimate.</footnote>
in which the functional is assumed to be Gateaux differ
entiable, that is

$$
Q'(u; \psi) = \left[ \frac{d}{d\epsilon} Q(u + \epsilon\psi) \right]_{\epsilon=0} \quad (123)
$$

We then may seek the error in the goal quantity, with
$\psi$ w $u_h$,

$$
a'(u_h; z, w \quad u_h) = Q'(u_h; w \quad u_h) \quad (124)
$$

where $Q'(u_h; w \quad u_h)$ is a linear approximation of the true error in the quantity of interest at the linearization point. Invoking the standard Galerkin orthogonality yields

$$
Q'(u_h; w \quad u_h) = a'(u_h; w \quad u_h, z \quad z_h) \quad (125)
$$

$$
= \sum_{K \in \mathcal{X}_h} \left\{ \int_K R(z \quad z_h) \, d\Omega \right. \\
\left. \qquad + \int_{\partial K} J(z \quad z_h) \, ds \right\} \quad (126)
$$

which can be employed in the dual weighted residual
method or, as well, in the Green's function decomposi-
tion method in case of point quantities. Clearly, employ
ing the Cauchy Schwarz inequality gives the error
bound

$$
| Q'(u_h; w \quad u_h) | \le \|z \quad z_h\|_E \|w \quad u_h\|_E \quad (127)
$$

It should be emphasized that as in the linear case we sim-
ply solve an additional load case at the linearization
point in order to find an estimate for the error in a quan
tity of interest.

It is interesting to note that the linearized error rep
resentation (124) corresponds to a linearization of
(121) with respect to $e_h$ if we proceed as in (113) (117)
for $\psi$ $z$ and assume that, as in the linear case, the
semi linear form $a(\cdot, \cdot)$ is symmetric.

A slightly different technique to estimate the error
in a (nonlinear) target functional was proposed in
[33,34,63]. This approach is based on the secant form

$$
a_T(u, u_h; \varphi, \psi) = \int_0^1 a'(u_h + s e_h; \varphi, \psi) \, ds \quad (128)
$$

which is still bilinear in $\varphi$ and $\psi$. Using the secant form, the error $e_h$ u u_h is the solution of the variational problem

$$
a_T(u, u_h; e_h, \psi) = I(\psi) \quad a(u_h; \psi) = \mathcal{R}_h(\psi) \quad \forall \psi \in V \quad (129)
$$

Here the exact solution $u$ on the left hand side in (129) is
unknown. One approach is to linearize about the finite
element approximation, and then of course we obtain
$a'(u_h; w \quad u_h, \psi)$ as used above.

To estimate the error in a functional of interest, we
similarly define the secant form of the functional

$$
Q_T(u, u_h; \psi) = \int_0^1 Q'(u_h + s e_h; \psi) \, ds \quad (130)
$$

where $Q_T$ is a linear functional in $\psi$. Here again $u$ is un
known and the straight linearization about $u_h$ of
$Q_T(u, u_h; \psi)$ results of course in obtaining $Q'(u_h; \psi)$.
However, it is possible to use enhanced solutions instead
of $u_h$ [4,63] and then a (somewhat) better approximation
to the error may be obtained.

For estimating the error in inelastic analysis, notably
elastoplasticity, several techniques have been presented,
see for example [64 68] and the references therein. In
general, in addition to the spatial approximation the
topic of this paper also the errors resulting in the inte-
gration of the stresses and the evaluation of the
constitutive relations need be considered. The estimation
of the complete solution errors in general plasticity solu-
tions is therefore much more difficult, and true error
bounds are not yet available. However, the techniques
reviewed above can directly be employed, as part of
the procedures to estimate the complete error in inelastic
analyses.

7. Time-dependent problems

The solution of time dependent problems is of special
interest in many engineering applications. From a nu-
merical point of view, in addition to the usual approxi-
mation in space, also a time discretization scheme has
to be applied. Our model transient problem is obtained
from the original model problem in Section 3, and is the
hyperbolic wave equation

$$
\ddot{u}(\mathbf{x}, t) - \beta \Delta u(\mathbf{x}, t) = \bar{f}(\mathbf{x}, t) \text{ in } \Omega_T = \Omega \times I; \quad I = (O, T) \tag{131}
$$

with $\beta$ a physical constant and the following initial and
boundary conditions

$$
u(\mathbf{x}, 0) = u_0; \quad \dot{u}(\mathbf{x}, 0) = \dot{u}_0 \text{ on } \Omega; \quad u(\mathbf{x}, t) = 0 \text{ on } \Gamma \tag{132}
$$

These equations govern the dynamic motion of a
prestressed membrane with $\beta \ N/\rho$, $\bar{f} = f/\rho$ where $N$
denotes the prestressing force, $\rho$ the mass density and $f$
is the load per unit area.

The commonly used finite element solution approach
to solve hyperbolic differential equations is to approxi-
mate the problem with finite elements in space and to
apply a finite difference scheme in time. The variational
formulation of (131) and (132) is:

Find a function $u(\mathbf{x}, t)$ that is for all $t \in (O, T)$ in the trial space $V = H_0^1(\Omega)$ such that

$$
(\ddot{u}, v) + a(u, v) = I(v) \quad \forall v \in V \quad (133)
$$
Then the solution is optimal in the sense

$$
\|\boldsymbol{\sigma} \quad \boldsymbol{\sigma}_h\|_{W} + \|u \quad u_h\|_{V} \\
\le c \left\{ \inf_{\tau_h \in W_h} \|\boldsymbol{\sigma} \quad \boldsymbol{\tau}_h\|_{W} + \inf_{v_h \in V_h} \|u \quad v_h\|_{V} \right\} \quad (150)
$$

Clearly, it is possible to formulate the mixed variational problem as follows

$$
\begin{align}
A(\mathbf{w}, \boldsymbol{\varphi}) &= (\boldsymbol{\sigma}, \boldsymbol{\tau}) + (u, \operatorname{div} \boldsymbol{\tau}) + (\operatorname{div} \boldsymbol{\sigma}, v) \nonumber \\
&= l(v) \quad \forall \boldsymbol{\varphi} \in V \times W \tag{151}
\end{align}
$$

with the symmetric bilinear form $A(\cdot, \cdot)$. Then for the
mixed finite element approximation we have

$$
A(\mathbf{w}_h, \varphi_h) = l(v_h) \quad (152)
$$

which directly leads to the Galerkin orthogonality for
the error **e**<sub>h</sub> **w** **w**<sub>h</sub>

$$
A(\mathbf{e}_h, \varphi_h) = (\boldsymbol{\sigma} \quad \boldsymbol{\sigma}_h, \boldsymbol{\tau}_h) + (u \quad u_h, \operatorname{div} \boldsymbol{\tau}_h) + (\operatorname{div}(\boldsymbol{\sigma} \quad \boldsymbol{\sigma}_h), v_h) \\
0 \quad \forall \varphi_h \in V_h \times W_h \quad (153)
$$

To derive an a posteriori error estimator for the mixed
problem we proceed as in Section 4.1. Thus, we obtain
the following error representation

$$
A(\mathbf{e}_h, \boldsymbol{\varphi}) = l(v) \quad A(\mathbf{w}_h, \boldsymbol{\varphi}) = \mathcal{R}_h(\boldsymbol{\varphi}) \quad (154)
$$

to hold true for each pair $\varphi$ $\{v, \tau\} \in V \times W$. If the
domain integral is split into its contributions from each
element, (154) can be rewritten as

$$
A(\mathbf{e}_h, \varphi) = \sum_{K \in \mathcal{F}_h} \left\{ \int_K f v \, d\Omega \quad \int_K \boldsymbol{\sigma}_h \cdot \tau \, d\Omega \right. \\
\left. \qquad \int_K u_h \, \operatorname{div} \tau \, d\Omega \quad \int_K \operatorname{div} \boldsymbol{\sigma}_h v \, d\Omega \right\}
$$

Applying integration by parts to the third integral and
assuming the finite element solution $u_h$ is sufficiently
smooth on $K$ yields

$$
\int_K u_h \operatorname{div} \tau \, d\Omega = \int_{\partial K} u_h \mathbf{n} \cdot \tau \, ds \quad \int_K \nabla u_h \cdot \tau \, d\Omega \quad (155)
$$

The boundary term does not vanish after summing over
all elements because the finite element solution $u_h$ is only
in $L^2(\Omega)$. Yet the expression $\mathbf{n} \cdot \tau$ is continuous at the ele-
ment boundaries since $\tau \in H^{\text{div}}(\Omega)$. Rearranging terms
yields

$$
\begin{equation}
\begin{split}
A(\mathbf{e}_h, \boldsymbol{\varphi}) = {}& \sum_{K \in \mathcal{F}_h} \left\{ \int_K R v \, d\Omega + \int_K \zeta \cdot \tau \, d\Omega \right\} \\
& + \sum_{\gamma \in \partial \mathcal{F}_h} \int_{\gamma} J \mathbf{n} \cdot \tau \, ds \quad \forall \boldsymbol{\varphi} \in V \times W
\end{split}
\tag{156}
\end{equation}
$$

where the element residuals are now defined by

$$
R = \operatorname{div} \boldsymbol{\sigma}_h \quad f \quad \text{in } K \qquad (157)
$$

$$
\zeta = \sigma_h + \nabla u_h \quad \text{in } K \qquad (158)
$$

and $J$ is the jump of the finite element solution $u_h$ across
the element edge $\gamma$ of elements $K$ and $K'$

$$
J = \begin{cases} (u_h, u'_h) & \text{if } \gamma \not\subseteq \Gamma \\ 0 & \text{if } \gamma \subset \Gamma \end{cases} \tag{159}
$$

Note that the normal vector **n** and the jumps change sign
if the orientation of the edge is reversed.

For driving an adaptive mesh iteration scheme, the
following local error indicators might now be used for
this model problem

$$
\eta_K^2 = h_K^2 \|R\|_{L^2(\kappa)}^2 + \|\zeta\|_{L^2(\kappa)}^2 + h_K \|J\|_{L^2(\partial\kappa)}^2 \quad (160)
$$

where on interelement boundaries ∂K ⊈ Γ the jump J
is multiplied by 1/2 to distribute the error equally onto
the two elements sharing the common edge.

Similar error indicators can be established for other
mixed formulations (see e.g. Ref. [75]), but the difficulty
in complex practical engineering analysis is to truly mea
sure the error in appropriate norms and establish error
bounds. Some discussion in this regard is given in Refs.
[50,76].

To obtain error estimates for the solution of a linear
quantity of interest $Q(\mathbf{w})$, we use a slightly different pro-
cedure than in the standard approach. Assuming that
the functional is well defined on the solution space
$V \times W$, the corresponding dual problem is:

Find $\mathbf{z} \{z, \pi\} \in V \times W$ such that

$$
\begin{align}
Q(\boldsymbol{\varphi}) &= (\boldsymbol{\pi}, \boldsymbol{\tau}) + (z, \operatorname{div} \boldsymbol{\tau}) + (\operatorname{div} \boldsymbol{\pi}, v) \nonumber \\
&= A(\mathbf{z}, \boldsymbol{\varphi}) \quad \forall \boldsymbol{\varphi} \in V \times W \tag{161}
\end{align}
$$

where **z** is the dual solution and ∇z **π** **0**. Taking
**φ** **e**<sub>h</sub> and utilizing the Galerkin orthogonality (153)
yields the error identity

$$
\begin{align}
Q(\mathbf{e}_h) = & ((\boldsymbol{\pi} \quad \boldsymbol{\pi}_h), (\boldsymbol{\sigma} \quad \boldsymbol{\sigma}_h)) + (z \quad z_h, \operatorname{div}(\boldsymbol{\sigma} \quad \boldsymbol{\sigma}_h)) \nonumber \\
& + (\operatorname{div}(\boldsymbol{\pi} \quad \boldsymbol{\pi}_h), u \quad u_h) \nonumber \\
= & A(\mathbf{z} \quad \mathbf{z}_h, \mathbf{w} \quad \mathbf{w}_h) \tag{162}
\end{align}
$$

To find an estimate for this quantity, we could write

$$
\begin{align}
|Q(\mathbf{e}_h)| &= |A(\mathbf{z} \quad \mathbf{z}_h, \mathbf{w} \quad \mathbf{w}_h)| \nonumber \\
&\le \|\mathbf{z} \quad \mathbf{z}_h\| \|\mathbf{w} \quad \mathbf{w}_h\| \tag{163}
\end{align}
$$

and apply some error estimators for the norm $||.||$. How
ever, note that we have $||\varphi|| \ne \sqrt{A(\varphi, \varphi)}$ since the bilin
ear form is not positive definite. Alternatively, we might
follow the suggestions of Rannacher and Suttmeier in
[35,68], who propose a method that is directly related
to the dual weighted residual method described in Sec
tion 5.2.2.

9. Numerical examples

In this section, we give some example solutions to
illustrate the use of error measures. Although a scalar
<figure><img src="image_10.png" /><figcaption>Fig. 9. Goal oriented error measures for a shell structure: problem description.</figcaption></figure>

model problem is considered in this paper, we consider in our numerical examples 2D linear elasticity and shell problems. We do not report upon actual practical engineering analyses but only give illustrative examples. Also, we are looking in some cases at very small errors, smaller than needed in practice, but do so because we want to study the convergence behavior of the procedure used.

## 9.1. Analysis of a shell

We consider a cylindrical shell that is partially clamped and loaded by a vertical line load as shown in Fig. 9. First, we study the results of a finite element computation on a reasonably fine mesh consisting of 20 × 20 MITC9 elements [51,52]. Fig. 10 shows that high stress concentrations are present in two regions near the clamped boundary (corresponding to the membrane stresses) and at the tip of the structure where the loading boundary conditions change (corresponding to the bending moment). Hence, using uniform MITC4 element meshes [52] and employing goal oriented error estimates based on (90), we want to evaluate the following quantities of interest (see also Fig. 11):

$$Q_1(\mathbf{U}) = |\Omega_1|^{-1} \int_{\Omega_1} \sigma_{ss}(\mathbf{U}) \, d\Omega \quad (164)$$

$$Q_2(\mathbf{U}) = |\Omega_2|^{-1} \int_{\Omega_2} \sigma_{ss}(\mathbf{U}) \, d\Omega \quad (165)$$

$$Q_3(\mathbf{U}) = |\Omega_3|^{-1} \int_{\Omega_3} z \sigma_{rr}(\mathbf{U}) \, d\Omega \quad (166)$$

These quantities are evaluated using the local Cartesian coordinate system $(r, s, z)$ (see Fig. 9). The reference solution was obtained using a uniform mesh of 100 × 100 MITC9 elements (with 201,000 degrees of

freedom) for the complete structure. For measuring the accuracy of our error estimate we use the effectivity index, which is the ratio between the estimated error using (90) and the calculated error using the reference solution defined in (89). As seen in Fig. 12 for every quantity of interest the estimated relative percentage error decreases quickly and the corresponding effectivity indices are close to 1.0.

Fig. 13 shows the predicted pointwise accuracy of the influence function for the quantity of interest defined in (165) when using the 20 × 20 mesh. We consider the absolute percentage value of the approximate error of the influence function normalized with the quantity of interest, $(Z_h^{(M9)} Z_h)/Q_2(\mathbf{U})$, and we define a tolerance on the absolute value, e.g. tol = 1.0%. Then the grey regions indicate that the error in the quantity of interest is smaller than the tolerance if the load is applied there, while the white areas correspond to errors larger than the tolerance. For additional shell analyses using this approach of error estimation see Ref. [49].

## 9.2. Analysis of a frame structure

Fig. 14 shows the frame structure considered. We want to employ adaptive mesh refinement based on energy norm estimates and the goal oriented techniques. For the analysis we use the 4 node displacement based bilinear element. The special focus is on the integrated shear stress in the two sections A A at level y = 4.0, so that the quantity of interest is:

$$Q(u_h) = \int_0^2 \sigma_{xy}^h \, dx + \int_{14}^{16} \sigma_{xy}^h \, dx \quad (167)$$

Local equilibrium asserts that the exact value is 144 which is equal to the volume above the cross section times the load. Clearly, the finite element solution yields
<figure><img src="image_11.png" /><figcaption>Fig. 10. Stresses of the shell structure under line load: (a) $\sigma_{ss}$ stress field and (b) $\sigma_{rr}$ stress field.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Goal oriented error measures for a shell structure: locations of three quantities of interest. 20 × 20 mesh shown.</figcaption></figure>

only an approximation for the integrated stresses be cause the stresses in the element interiors and the fluxes at the boundaries are only approximations to the true values. Note that, of course, the finite element nodal forces in the cross section A A are in exact equilibrium with the applied load [1].

For the solution we use two approaches: first, the refinement is based on the norms given in (28) and (29) and, second, the refinement is based on the goal ori

ented error indicators given in (77) (79). In both ap proaches we use the constants <fnref n="5" /> [13]

$$c_1 = \frac{0.16}{\lambda + 5\mu} \qquad c_2 = \frac{1.44}{\lambda + 5\mu} \qquad \lambda, \mu = \text{Lamé constants} \tag{168}$$

and the hanging node concept, in which the unknown displacements at a hanging node (a node not shared by all elements surrounding it) are eliminated by inter polation using the neighboring nodes. Also, in the ap proaches, the refinement criteria defined in (31) and (81) are used, respectively, and the 30% of the elements with $\xi_K > 1$ (and corresponding to the largest values of $\xi_K$) are refined.

Fig. 15 compares the results obtained using these two solution approaches with the results obtained by simply using a uniform mesh refinement. In all these solutions, the quantity of interest is obtained by differentiation of the finite element displacement field to evaluate the stres ses. Fig. 15 shows that the goal oriented approach leads,

<footnote n="5">⁵ The numerator in the definition of the constant $c_2$ in (168) is changed compared to the value of 1.21 given in [13] since we are using quadrilaterals instead of triangles, for which the constants originally were derived.</footnote>
<figure><img src="image_13.png" /><figcaption>Fig. 12. Results of goal oriented error estimation for the shell structure shown in Fig. 9: (a) estimated (absolute) relative errors in the quantities of interest and (b) corresponding effectivity indices. The numbers of degrees of freedom refer to the uniform MITC4 element meshes used for the complete shell.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Predicted absolute pointwise accuracy of the influence function normalized with the quantity of interest defined in (165): body loads in the *z* direction in the grey region yield an error in the quantity of interest smaller than a tolerance of (a) tol = 1.0%, (b) tol = 0.8%, (c) tol = 0.6%, (d) tol = 0.4%, (e) tol = 0.2%, (f) tol = 0.1%.</figcaption></figure>

as expected, to the smallest error in the quantity of interest. The initial and final meshes are shown in Fig. 16. The refinement based on (28) and (29) did not refine in the cross section but at the singularities of the structure, whereas the density of the goal oriented mesh obtained using (77) (79) is the highest in the cross section. Actually, for the range of degrees of freedom considered,

there is only a slight improvement in the quantity of interest employing the energy norm for refinement of the mesh.

The stress error indicator in ADINA plotted in Fig. 16(a) is obtained by taking the difference between the maximum and minimum von Mises stress at the nodes, normalized to the maximum value [22]. These results
<figure><img src="image_15.png" /><figcaption>Fig. 14. Frame under horizontal body load to study the error in the integrated shear stresses in the two sections A A.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Comparison of different adaptive refinements to approximate the integrated shear stresses in the two sections A A.</figcaption></figure>

show that, as expected, the error indicator for the von Mises stress obtained with ADINA corresponds quite well with the energy norm refinement, so that, indeed, this error indicator would have provided a good guide for driving a mesh refinement “by hand” without any usage of an error estimator for the energy norm.

Comparing the computational effort we need to keep in mind that for the goal oriented approach two finite element solutions per mesh are necessary. However, the solution of the dual problem differs only in a modi fied load vector from the analysis of the primal problem; that is, the stiffness matrix remains the same.

9.3. Analysis of a plate in plane stress

Next, we study the example described in Fig. 17. The focus is on the stresses $\sigma_{xx}$ and $\sigma_{yy}$ at the point A. The

<figure><img src="image_17.png" /><figcaption>Fig. 16. Results obtained in the analysis of the frame in Fig. 15: (a) Initial mesh and error indicators obtained with ADINA, (b) final mesh obtained with the energy norm control and (c) final mesh obtained with the goal oriented control.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 17. Plate in plane stress condition to study the error in the stresses at the point A: (a) problem data, (b) $\sigma_{xx}$ stress field, (c) $\sigma_{yy}$ stress field and (d) $\sigma_{xy}$ stress field.</figcaption></figure>

<table><caption>Table 1
Stresses at A employing the Green's function decomposition approach using uniform meshes</caption><thead><tr><th>h</th><th>DOF</th><th>σ<sub>xx</sub></th><th>|e<sub>rel</sub>| [%]</th><th>σ<sub>yy</sub></th><th>|e<sub>rel</sub>| [%]</th></tr></thead><tbody><tr><td>1/5</td><td>32</td><td>53.8018</td><td>2.931 · 10<sup>1</sup></td><td>33.9038</td><td>1.055 · 10<sup>0</sup></td></tr><tr><td>1/10</td><td>162</td><td>53.6705</td><td>4.884 · 10<sup>2</sup></td><td>33.6288</td><td>2.349 · 10<sup>1</sup></td></tr><tr><td>1/15</td><td>392</td><td>53.6540</td><td>1.808 · 10<sup>2</sup></td><td>33.5831</td><td>9.866 · 10<sup>2</sup></td></tr><tr><td>1/20</td><td>722</td><td>53.6492</td><td>9.134 · 10<sup>3</sup></td><td>33.5681</td><td>5.395 · 10<sup>2</sup></td></tr><tr><td>1/25</td><td>1152</td><td>53.6473</td><td>5.592 · 10<sup>3</sup></td><td>33.5614</td><td>3.398 · 10<sup>2</sup></td></tr><tr><td>1/32</td><td>1922</td><td>53.6460</td><td>3.169 · 10<sup>3</sup></td><td>33.5569</td><td>2.057 · 10<sup>2</sup></td></tr><tr><td>1/40</td><td>3042</td><td>53.6454</td><td>2.051 · 10<sup>3</sup></td><td>33.5544</td><td>1.311 · 10<sup>2</sup></td></tr><tr><td>1/50</td><td>4802</td><td>53.6450</td><td>1.305 · 10<sup>3</sup></td><td>33.5528</td><td>8.346 · 10<sup>3</sup></td></tr><tr><td>Reference value</td><td></td><td>53.6443</td><td></td><td>33.5500</td><td></td></tr></tbody></table>

reference stress fields are obtained using a very fine mesh. In all analyses the 4 node displacement based bilinear element is used.

First, we give some results obtained employing the Green's function decomposition method using uniform meshes. As a result of this uniform refinement shown in Table 1, the error in the stresses decreases reasonably fast, and indeed highly accurate results are achieved even on coarse meshes.

It is interesting to compare the results obtained with the Green's function decomposition method given in Table 1 with the results obtained using different mesh refinement techniques (see Fig. 18(a)). First, we calculate finite element approximations on uniformly refined meshes. Next, the refinement is steered using the global energy norm control based on the explicit error estima tor in (28) and (29) with the constants given in (168), and in a third solution the refinement is obtained using the
<figure><img src="image_19.png" /><figcaption>Fig. 18. Comparison of results obtained using different mesh refinement techniques for the stresses at point A: (a) relative error in stresses, (b) related error bounds (GFD are the Green's function decomposition method results in Table 1).</figcaption></figure>

energy norm based goal oriented strategy based on (77) using also the constants in (168). For the energy norm control, the refinement criterion in (31) is used, while for the energy norm based goal oriented strategy the refinement criterion given in (81) is employed. In both approaches 25% of the elements with $\xi_K > 1$ (and corresponding to the largest values of $\xi_K$) are refined.

In the goal oriented error estimate, to solve the dual problem in (64) and the corresponding nodal forces in (65), we consider integrated stresses over a small circular domain and take numerically the limit value as the radius of the domain tends to zero. For evaluating the stress quantity on a given mesh, we simply differentiate the finite element displacement field, as usual in finite element analysis.

As can be seen in Fig. 18, the stresses obtained with the uniformly refined meshes, the global energy norm control, and the goal oriented approach are somewhat erratic since they depend strongly on the current mesh design. In contrast, the stresses given in Table 1 show convergence rates of algebraic type on a doubly logarithmic scale, and the results are highly accurate compared with the global energy norm approach and the energy norm based goal oriented approach.

In addition, we calculate the error bounds for the stresses obtained with the Green's function decomposition approach on uniform meshes based on (102) (see Fig. 18(b)). Here, we use the explicit error estimate de

fined in (28) for the energy norm errors used in (102) with the constants in (168). As seen, these error bounds produce an envelope which contains the exact solution, and which is quite narrow although the estimates for the local errors are based on explicit error estimators.

Next, we employ the energy norm based goal oriented approach and the energy norm based Green's function decomposition approach, starting with a coarse mesh of 5 × 5 elements. We steer in each case the refinement process using (81) and the local error indicators given in (78) and (79), for the primal, the dual and the regular part problems. The refinement is performed using $\xi_K$ as above.

In Fig. 19 the final meshes obtained are displayed. As expected, the refinement of the energy norm based goal oriented procedure concentrates on the region near the point of interest. The energy norm based Green's function decomposition approach can neglect this region because the fundamental solution is already the quasi optimal choice for the region. Instead, the mesh is refined almost uniformly, which in this example is sufficient to approximate the regular part accurately.

While the Green's function decomposition approach is clearly very effective in this illustrative example, as already pointed out above, the method is however rather restrictive in that only analysis cases can be considered for which the fundamental solution is available.
<figure><img src="image_20.png" /><figcaption>Fig. 19. Final meshes to approximate the stress $\sigma_{xx}$ at point A: (a) goal oriented refinement, (b) refinement using the Green's function decomposition approach.</figcaption></figure>

## 10. Conclusions

In this paper we reviewed some basic a posteriori error estimation techniques which broadly can be classi-fied into global error estimators for the energy norm and goal oriented error estimators to provide error estimates and error bounds for linear quantities of interest. We also discussed the case when the goal quantity is a point value, which normally poses a difficulty since the dual solution is not in the solution space.

The crucial issues of any error estimator relate to questions of reliability, accuracy and computational cost where, clearly, the definition what is an admissible cost always depends on the purpose of the computation [1]. As pointed out, there may be different reasons for using an error estimator. To only obtain an indication of the error or to drive a mesh adaptivity scheme with respect to the goal of the computation, explicit error indicators might be sufficient and these are generally quite inexpensive to use. On the other hand, to actually bound the error almost guaranteed in a suitable norm is at present only possible for certain problems, and then very expensive in practical analysis. However, it is also clear that the more accurate estimators have to be more expensive because a nearly exact estimate is close to finding the exact solution. The question is then whether it is not more effective in practice to simply use a very fine mesh.

A major point is that, in essence, error bounds are either guaranteed but, in practice for complex problems, hardly computable or they are computable but not guaranteed. Another major point is that while, in engineering practice, the analysis of shell structures constitutes a large percentage of all linear analyses, there are only some contributions that address bounding the error in suitable norms [139, 50, 61, 76] and goal oriented procedures need still be further explored, in particular using mixed finite element discretizations [49].

Since sharp and effective error estimators are not yet available for many practical analyses, we are left with the common advice:

* The mesh should be reasonably fine and, ideally, solutions would be obtained for a coarse mesh and a finer mesh, for comparison purposes. Also, simple visual criteria to approximately assess the error might be used such as the stress iso bands in ADINA.

* At all high stress gradients the mesh should be sufficiently fine. The results of the analysis are frequently most accurate when the relative error is uniform over the complete analysis domain.

* Integrated (averaged) quantities are usually more accurate than point values.

* In order to obtain highly accurate local quantities of interest, the mesh should of course be reasonably fine around the region of interest (and also in general in the areas of high stress gradients as mentioned already).

Future research work should address the development of actually implementable and practically useful im-proved error estimators that are applicable for a large class of problems. From an engineering point of view, inexpensive to calculate and guaranteed bounds on the error at every point of the structure, in the sense of our ideal error estimator solution (see Fig. 1), would be very valuable. However, such error bound solution will likely be very difficult to achieve without a significant computational expense. In engineering practice, the calculation of the error measure should not be more expensive than the added expense to simply run a very fine mesh.

A quite promising approach is to use goal oriented error measures in order to establish a coarse but still appropriate mesh in computationally intensive finite ele-
ment solutions, notably in multi physics and multi scale
analyses involving optimization. Here the premise is that
an integrated quantity might be calculated with suffi
cient accuracy using a well chosen coarse mesh in part
of the domain. For example, considering a fluid flow
structural interaction analysis, a coarse mesh represent
ing the fluid might be sufficient to calculate the total
force and moment on the structure [81]. And the compu
tational expense to establish and use the coarse mesh of
the fluid might be much less than the expense to use a
very fine fluid mesh, in particular, if a structural optimi
zation is required.

Hence, while the theory of error estimation has pro
vided much valuable insight into the finite element solu
tion process, many of the proposed techniques are at
present only valuable to a limited extent in engineering
practice.

Of course, throughout the paper we assumed that an
appropriate mathematical model has been chosen and
we only focused on the discretization errors arising in
the finite element solution of this model (see Section
1). In practical engineering analysis, the errors arising
due to an inappropriate mathematical model can natu
rally be much more significant than the error we have
discussed in this paper [1,82].

Acknowledgement

We would like to thank Donald Estep, Colorado
State University, Ivo Babuška, University of Texas at
Austin, Zhimin Zhang, Wayne State University, Jaime
Peraire, MIT, Anthony Patera, MIT, and Slimane Adje
rid, Virginia Tech, for their valuable comments on this
paper.

The work of Thomas Grätsch was supported by the
German Research Foundation (DFG) under contract
GR 1894/2 1. We are grateful for this support.

References

[1] Bathe KJ. Finite element procedures. Englewood Cliffs, NJ: Prentice Hall; 1996.

[2] Szabó B, Babuška I. Finite element analysis. New York, USA: John Wiley & Sons Inc; 1991.

[3] Verfurth R. A review of a posteriori error estimation and adaptive mesh refinement techniques. Chichester: John Wiley & Sons; 1996.

[4] Estep D, Larson M, Williams R. Estimating the error of numerical solutions of systems of reaction diffusion equations. American Mathematical Society, No. 696; 2000.

[5] Ainsworth M, Oden JT. A posteriori error estimation in finite element analysis. New York: John Wiley & Sons; 2000.

[6] Babuška I, Strouboulis T. The finite element method and its reliability. Oxford: Oxford Science Publications; 2001.

[7] Brenner SC, Scott LR. The mathematical theory of finite element methods. Texts in applied mathematics Bd. 15. New York: Springer Verlag; 1994.

[8] Johnson C. Numerical solution of partial differential equations by the finite element method. Cambridge: Cam bridge University Press; 1987.

[9] Babuška I, Rheinboldt WC. A posteriori error estimates for the finite element method. Int J Numer Methods Engng 1978;12:1597 615.

[10] Babuška I, Miller AD. A feedback finite element method with a posteriori error estimation. Part I: The finite element method and some basic properties of the a posteriori error estimator. Comput Methods Appl Mech Engng 1987;61: 1 40.

[11] Babuška I, Rheinboldt WC. Analysis of optimal finite element meshes in $R^1$. Math Comput 1978;33:435 63.

[12] Kelly DW, Gago OC, Zienkiewicz OC, Babuška I. A posteriori error analysis and adaptive processes in the finite element method: Part I error analysis. Int J Numer Methods Engng 1983;19:1593 619.

[13] Johnson C, Hansbo P. Adaptive finite element methods in Computational Mechanics. Comput Methods Appl Mech Engng 1992;101:143 81.

[14] Stein E, Seifert B, Ohnimus S, Carstensen C. Adaptive finite element analysis of geometrically non linear plates and shells, especially buckling. Int J Numer Methods Engng 1994;37:2631 55.

[15] Bank RE, Weiser A. Some a posteriori error estimators for elliptic partial differential equations. Math Comput 1985;44:283 301.

[16] Babuška I, Strouboulis T, Gangaraj SK. Guaranteed computable bounds for the exact error in the finite element solution. Part I: One dimensional model problem. Comput Methods Appl Mech Engng 1999;176:51 79.

[17] Babuška I, Rheinboldt WC. Error estimates for adaptive finite element computations. SIAM J Numer Anal 1978;15: 736 54.

[18] Carstensen C, Funken S. Fully reliable localized error control in the FEM. SIAM J Sci Comput 2000;21:1465 84.

[19] Machiels L, Maday Y, Patera A. A flux free nodal Neumann subproblem approach to output bounds for partial differential equations. CR Acad. Sci. Paris, vol. 330. Serie I; 2000. p. 249 54.

[20] Morin P, Nochetto RH, Siebert KG. Local problems on stars: a posteriori error estimators, convergence, and performance. Math Comput 2003;72:1067 97.

[21] Prudhomme S, Nobile F, Oden JT. Analysis of a subdo main based error estimator for finite element approxima tions of elliptic problems. TICAM Report 02 34. The University of Texas at Austin, 2002. Numer Meth Part D and E, in press.

[22] Sussman T, Bathe KJ. Studies of finite element proce dures stress band plots and the evaluation of finite element meshes. Engng Comput 1986;3:178 91.

[23] Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptive procedure for practical engineering analysis. Int J Numer Methods Engng 1987;24:337 57.

[24] Carstensen C, Funken S. Averaging technique for FE a posteriori error control in elasticity. Part I: conform ing FEM. Comput Methods Appl Mech Engng 2001;190: 2483 98.
[25] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique. Int J Numer Methods Engng 1992; 33:1331 64.

[26] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity. Int J Numer Methods Engng 1992;33:1365 82.

[27] Zhang Z, Naga A. A new finite element gradient recovery method: superconvergence property. SIAM J Sci Comput, in press.

[28] Xu J, Zhang Z. Analysis of recovery type a posteriori error estimators for mildly structured grids. Math Comput 2004;73:1139 52.

[29] Hiller JF, Bathe KJ. On higher order accuracy points in isoparametric finite element analysis and an application to error assessment. Comput Struct 2001;79:1275 85.

[30] Babuška I, Miller AD. The post processing approach in the finite element method, I: calculations of displacements, stresses and other higher derivatives of the displacement. Int J Numer Methods Engng 1984;20:1085 109.

[31] Eriksson K, Johnson C. Adaptive finite element methods for parabolic problems I: a linear model problem. SIAM J Numer Anal 1991;28:43 77.

[32] Estep D. A posteriori error bounds and global error control for approximations of ordinary differential equations. SIAM J Numer Anal 1995;32:1 48.

[33] Eriksson K, Estep D, Hansbo P, Johnson C. Introduction to adaptive methods for differential equations. Acta Numer 1995;105 58.

[34] Becker R, Rannacher R. A feed back approach to error control in finite element methods: basic analysis and examples. East West J Numer Math 1996;4: 237 64.

[35] Rannacher R, Suttmeier FT. A feed back approach to error control in finite element methods: application to linear elasticity. Comput Mech 1997;19:434 46.

[36] Giles MB, Larson MG, Levenstam JM, Suli E. Adaptive error control for finite element approximations of the lift and drag coefficients in viscous flow. Report NA 97/06. Oxford University Computing Laboratory; 1997.

[37] Paraschivoiu M, Peraire J, Patera A. A posteriori finite element bounds for linear functional outputs of elliptic partial differential equations. Comput Methods Appl Mech Engng 1997;150:289 312.

[38] Paraschivoiu M, Patera A. A hierarchical duality approach to bounds for the outputs of partial differential equations. Comput Methods Appl Mech Engng 1998;158:389 407.

[39] Çirak F, Ramm E. A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem. Comput Methods Appl Mech Engng 1998; 156:351 62.

[40] Pierce NA, Giles MB. Adjoint recovery of superconvergent functionals from PDE approximations. Siam Rev 2000;42: 247 64.

[41] Giles MB, Suli E. Adjoint methods for PDEs: a posteriori error analysis and postprocessing by duality. Acta Numer 2002;145 236.

[42] Tottenham H. Basic principles. In: Tottenham H, Brebbia C, editors. Finite element techniques in structural mechanics. Southampton: Southampton University Press; 1970.

[43] Adjerid S, Salim M. Even odd goal oriented a posteriori error estimation for elliptic problems. Appl Numer Math, in press.

[44] Prudhomme S, Oden JT. On goal oriented error estimation for elliptic problems: application to the control of point wise errors. Comput Methods Appl Mech Engng 1999;176:313 31.

[45] Oden JT, Prudhomme S. Goal oriented error estimation and adaptivity for the finite element method. Comput Math Appl 2001;41:735 56.

[46] Babuška I, Miller AD. The post processing approach in the finite element method, III: a posteriori error estimates and adaptive mesh selection. Int J Numer Methods Engng 1984;20:2311 24.

[47] Ohnimus S, Stein E, Walhorn E. Local error estimates of FEM for displacements and stresses in linear elasticity by solving local Neumann problems. Int J Numer Methods Engng 2001;52:727 46.

[48] Bangerth W, Rannacher R. Adaptive finite element methods for differential equations. Basel: Birkhauser; 2003.

[49] Gratsch T, Bathe KJ. Influence functions and goal oriented error estimation for finite element analysis of shell structures. Int J Numer Methods Engng, submitted for publication.

[50] Chapelle D, Bathe KJ. The finite element analysis of shells fundamentals. Berlin: Springer Verlag; 2003.

[51] Bathe KJ, Lee PS, Hiller JF. Towards improving the MITC9 shell element. Comput Struct 2003;81: 1085 97.

[52] Bathe KJ, Iosilevich A, Chapelle D. An evaluation of the MITC shell elements. Comput Struct 2000;75:1 30.

[53] Wang S, Sloan IH, Kelly DW. Computable error bounds for pointwise derivatives of a Neumann problem. IMA J Numer Anal 1998;18:251 71.

[54] Cao D, Kelly DW, Sloan IH. Local error bounds for post-processed finite element calculations. Int J Numer Methods Engng 1999;45:1085 98.

[55] Gratsch T, Hartmann F. Finite element recovery techniques for local quantities of linear problems using fundamental solutions. Comput Mech 2003;33:15 21.

[56] Hartmann F, Katz C. Structural mechanics with finite elements. Berlin: Springer Verlag; 2004.

[57] Sauer Budge AM, Bonet J, Huerta A, Peraire J. Computing bounds for linear functionals of exact solutions to Poisson's equation. SIAM J Numer Anal, in press.

[58] Sauer Budge AM, Peraire J. Computing bounds for linear functionals of exact solutions to the advection diffusion reaction equation. SIAM J Sci Comput, in press.

[59] Prud'homme C, Rovas DV, Veroy K, Machiels L, Maday Y, Patera A, Torinici G. Reliable real time solution of parametrized partial differential equations: reduced basis output bound methods. ASME 2002;124:70 80.

[60] Estep D, Holst M, Larson M. Generalized Green's functions and the effective domain of influence. SIAM J Sci Comput, in press.

[61] Prudhomme S, Oden JT, Westermann T, Bass J, Botkin M. Practical methods for a posteriori error estimation in engineering applications. Int J Numer Methods Engng 2003;56:1193 224.
[62] Caloz G, Rappaz J. Numerical analysis for nonlinear and bifurcation problems, Part 2. Ciarlet PG, Lions JL, editors. Handbook of numerical analysis, Vol. V. Amsterdam: Elsevier; 1997. p. 487 638.

[63] Larsson F, Hansbo P, Runesson K. Strategies for computing goal oriented a posteriori error measures in non linear elasticity. Int J Numer Methods Engng 2002;55: 879 94.

[64] Çirak F, Ramm E. A posteriori error estimation and adaptivity for elastoplasticity using the reciprocal theorem. Int J Numer Methods Engng 2000;47:379 93.

[65] Larsson F, Runesson K, Hansbo P. Time finite elements and error computation for (visco) plasticity with hardening or softening. Int J Numer Methods Engng 2003;56: 2213 32.

[66] Lee NS, Bathe KJ. Error indicators and adaptive reme shing in large deformation finite element analysis. Finite Elem Anal Des 1994;16:99 139.

[67] Gallimard L, Ladèveze P, Pelle JP. Error estimation and adaptivity in elastoplasticity. Int J Numer Methods Engng 1996;39:189 217.

[68] Rannacher R, Suttmeier FT. A posteriori error estimation and mesh adaption for finite element models in elasto-plasticity. Comput Methods Appl Mech Engng 1999;176: 333 61.

[69] Suli E, Houston P. Finite element methods for hyperbolic problems: a posteriori error analyses and adaptivity. Invited Lecture at State of the Art in Numerical Analysis Conference. York, 1 4 April, 1996. Oxford University Computing Laboratory, Report 96/06, 1996.

[70] Schleupen A, Ramm E. Local and global error estimation in linear structural dynamics. Comp Struct 2000;76:741 756.

[71] Wiberg NE, Li XD. A postprocessed error estimate and an adaptive procedure for the semidiscrete finite element

method in dynamic analysis. Int J Numer Methods Engng 1994;37:3585 603.

[72] Brezzi F, Fortin M. Mixed and hybrid finite elements methods. New York: Springer Verlag; 1991.

[73] Bathe KJ. The inf sup condition and its evaluation for mixed finite element methods. Comput Struct 2001;79: 243 52. p. 971.

[74] Brezzi F, Bathe KJ. A discourse on the stability conditions for mixed finite element formulations. Comput Methods Appl Mech Engng 1990;82:27 57.

[75] Braess D, Klaas O, Niekamp R, Stein E, Wobschal F. Error indicators for mixed finite elements in 2 dimensional linear elasticity. Comput Methods Appl Mech Engng 1995;127:345 56.

[76] Hiller JF, Bathe KJ. Measuring convergence of mixed finite element discretizations: An application to shell structures. Comput Struct 2003;81:639 54.

[77] Huerta A, Rodriguez Ferran A, Diez P, Sarrate J. Adaptive finite element strategies based on error assessment. Int J Numer Methods Engng 1999;46:1803 18.

[78] Yazdani AA, Riggs HR, Tessler A. Stress recovery and error estimation for shell structures. Int J Numer Methods Engng 2000;47:1825 40.

[79] Duster A, Broker H, Rank E. The p version of the finite element method for three dimensional curved thin walled structures. Int J Numer Methods Engng 2001;52: 673 703.

[80] Oden JT, Prudhomme S, Westermann T, Bass J, Botkin M. Error estimation of eigenfrequencies for elasticity and shell problems. Math Mod Methods Appl S 2003;13: 323 44.

[81] Gratsch T, Bathe KJ. Goal oriented error estimation in the analysis of fluid flows with structural interactions, in preparation.

[82] Bucalem M, Bathe KJ. Hierarchical modeling in finite element analysis. Springer Verlag, in preparation.