the discretization and then to perform a new run of the solution process.

A typical modi¬cation is to re¬ne or to coarsen the grid.

Error estimators may overestimate the real error measure signi¬cantly;

thus a grid adaptation procedure based on such an error estimate gener-

ates a grid that is too ¬ne, and consequently, the corresponding algebraic

problem is too large.

This e¬ect can be reduced or even avoided if the error estimator satis¬es

a two-sided inequality like (4.4). Then the ratio D2 /D1 is a measure of the

e¬ciency of the error estimator.

4.2. A Posteriori Error Estimates and Grid Adaptation 187

An error estimator · is called asymptotically exact if for an arbitrary con-

vergent sequence of approximations {uh } with u ’ uh ’ 0 the following

limit is valid:

·

’ 1.

u ’ uh

Usually, a posteriori error estimators are designed for a well-de¬ned class

of boundary or initial-boundary value problems. Within a given class of

problems, the question regarding the sensitivity of the constants D in (4.3)

or D1 , D2 in (4.4), with respect to the particular data of the problem (e.g.,

coe¬cients, inhomogeneities, geometry of the domain, grid geometry, . . . ),

arises. If this dependence of the data is not crucial, then the error estimator

is called robust within this class.

Grid Adaptation

Let us assume that the local error estimators ·K composing an e¬cient er-

ror estimator · for an approximate solution uh on some grid Th really re¬‚ect

the error on the element K and that this local error can be improved by a

re¬nement of K (e.g., following the principles of Section 4.1.5). Then the

following grid adaptation strategies can be applied until the given tolerance

µ is reached or the computer resources are exhausted.

Equidistribution strategy: The objective of the grid adaptation (re¬nement

or coarsening of elements) is to get a new grid Th new

such that the

new

local error estimators ·K for this new grid take one and the same

value for all elements K ∈ Th ; that is (cf. (4.5))

new

µ

·K ≈ for all K ∈ Th .

new new

|Th |

new

Since the number of elements of the new grid enters the right-hand

side of this criterion, the strategy is an implicit method. In practical

use, it is approximated iteratively.

Cut-o¬ strategy: Given a parameter κ ∈ (0, 1), a threshold value κ· is

de¬ned. Then the elements K with ·K > κ· will be re¬ned.

Reduction strategy: Given a parameter κ ∈ (0, 1), an auxiliary toler-

ance µ· := κ· is de¬ned. Then a couple of steps following the

equidistribution strategy with the tolerance µ· are performed.

In practice, the equidistribution strategy may perform comparatively slowly

and thus may reduce the e¬ciency of the complete solution process. The

cut-o¬ method does not allow grid coarsening. It is rather sensitive to the

choice of the parameter κ. Among all three strategies, the reduction method

represents the best compromise.

188 4. Grid Generation and A Posteriori Error Estimation

Design of A Posteriori Error Estimators

In the following, three basic principles of the design of a posteriori error

estimators will be described. In order to illustrate the underlying ideas and

to avoid unnecessary technical di¬culties, a model problem will be treated:

Consider a di¬usion-reaction equation on a polygonally bounded domain

„¦ ‚ R2 with homogeneous Dirichlet boundary conditions

’∆u + ru = f in „¦ ,

u= 0 on ‚„¦ ,

where f ∈ L2 („¦) and r ∈ C(„¦) with r(x) ≥ 0 for all x ∈ „¦. The problem

is discretized using piecewise linear, continuous ¬nite element functions as

described in Section 2.2.

Setting a(u, v) := „¦ (∇u · ∇v + ruv) dx for u, v ∈ V := H0 („¦), we have

1

the following variational (weak) formulation:

Find u ∈ V such that a(u, v) = f, v for all v ∈ V.

0

The corresponding ¬nite element method reads as follows:

Find uh ∈ Vh such that a(uh , vh ) = f, vh for all vh ∈ Vh .

0

Residual Error Estimators

Similar to the derivation of the a priori error estimate in the proof of C´a™s

e

lemma (Theorem 2.17), the V -ellipticity of a (2.43) implies that

± u ’ uh ¤ a(u ’ uh , u ’ uh ) .

2

1

Without loss of generality we may suppose u ’ uh ∈ V \ {0}, hence

1 a(u ’ uh , u ’ uh ) a(u ’ uh , v)

1

u ’ uh ¤ ¤ sup . (4.6)

1

u ’ uh 1

± ± v∈V v1

We observe that the term

a(u ’ uh , v) = a(u, v) ’ a(uh , v) = f, v ’ a(uh , v) (4.7)

0

is the residual of the variational equation; i.e., the right-hand side of in-

equality (4.6) can be interpreted as a certain norm of the variational

residual.

In a next step, the variational residual will be split into local terms

according to the given grid, and these terms are transformed by means of

integration by parts. For arbitrary v ∈ V, from (4.7) it follows that

a(u ’ uh , v) = f v dx ’ (∇uh · ∇v + ruh v) dx

K K

K∈Th

[f ’ (’∆uh + ruh )]v dx ’ ν · ∇uh v dσ

= .

K ‚K

K∈Th

4.2. A Posteriori Error Estimates and Grid Adaptation 189

The ¬rst factor in the integrals over the elements K is the classical

elementwise residual of the di¬erential equation:

rK (uh ) := [f ’ (’∆uh + ruh )] K

All quantities entering rK (uh ) are known. In the case considered here we

even have ’∆uh = 0 on K, hence rK (uh ) = [f ’ ruh ] K .

The integrals over the boundary of the elements K are further split into

a sum over the integrals along the element edges E ‚ ‚K:

ν · ∇uh v dσ = ν · ∇uh v dσ .

‚K E

E‚‚K

Since v = 0 on ‚„¦, only the integrals along edges lying in „¦ contribute to

the sum. Denoting by Eh the set of all interior edges of all elements K ∈ Th

and assigning a ¬xed unit normal νE to any of those edges, we see that in

the summation of the split boundary integrals over all K ∈ Th there occur

exactly two integrals along one and the same edge E ∈ Eh . This observation

results in the relation

ν · ∇uh v dσ = [νE · ∇uh ]E v dσ ,

‚K E

K∈Th E∈Eh

where, for a piecewise continuous function w : „¦ ’ R, the term

[w]E (x) := lim w(x + νE ) ’ lim w(x ’ νE ) , x∈E,

’+0 ’+0

denotes the jump of the function w across the edge E. If w is the normal

derivative of uh in the ¬xed direction νE , i.e., w = νE · ∇uh , then its jump

does not depend on the particular orientation of νE (see Exercise 4.6).

In summary, we have the following relation:

a(u ’ uh , v) = rK (uh )v dx ’ [νE · ∇uh ]E v dσ .

K E

K∈Th E∈Eh

Using the error equation (2.39), we obtain for an arbitrary element vh ∈ Vh

the fundamental identity

a(u ’ uh , v) = a(u ’ uh , v ’ vh )

rK (uh )(v ’ vh ) dx

=

K

K∈Th

’ [νE · ∇uh ]E (v ’ vh ) dσ ,

E

E∈Eh

which is the starting point for the construction of further estimates.

190 4. Grid Generation and A Posteriori Error Estimation

¨¨r r

r

¨¡

¢d ¢f ¢d ¢f

¡

d¡ d

¢f ¢f

¢ ¢d ¢ ¢d

d¢f d¢f

¢¢ ¢¢

¢ ¢ Kd¢ f ¢ ¢ E ¢d f

¢d ¢ ¢ d

d d¢

¢ e ¢ e

¢ e ¨ ¢ e ¨

¨ ¨

¢ e¨ ¢

e¨

Figure 4.6. The triangular neighbourhoods ∆(K) (left) and ∆(E) (right).

First we see that the Cauchy“Schwarz inequality immediately implies

a(u ’ uh , v ’ vh ) ¤ v ’ vh

rK (uh ) 0,K 0,K

K∈Th

(4.8)

[νE · ∇uh ]E v ’ vh

+ 0,E .

0,E

E∈Eh

To get this bound as small as possible, the function vh ∈ Vh is chosen such

that the element v ∈ V is approximated adequately in both spaces L2 (K)

and L2 (E). One suggestion is the use of an interpolating function according

to (2.47). However, since V ∈ C(„¦), this interpolant is not de¬ned. There-

fore other approximation procedures have to be applied. Roughly speaking,

suitable approximation principles, due to Cl´ment [52] or Scott and Zhang

e

[67], are based on taking certain local integral means. However, at this place

we cannot go further into these details and refer to the cited literature. In

fact, for our purposes it is important only that such approximations exist.

Their particular design is of minor interest.

We will formulate the relevant facts as a lemma. To do so, we need some

additional notation (see Figure 4.6):

triangular neighbourhood of a triangle K : ∆(K) := K :K ©K=… K ,

triangular neighbourhood of an edge E : ∆(E) := K.

K :K ©E=…

Thus ∆(K) consists of the union of the supports of those nodal basis

functions that are associated with the vertices of K, whereas ∆(E) is formed

by the union of those nodal basis functions that are associated with the

boundary points of E. Furthermore, the length of the edge E is denoted

by hE := |E|.

Lemma 4.2 Let a regular family (Th ) of triangulations of the domain „¦

be given. Then for any v ∈ V there exists an element Qh v ∈ Vh such that

for all triangles K ∈ Th and all edges E ∈ Eh the following estimates are

valid:

v ’ Qh v ¤ ChK |v|1,∆(K) ,

0,K

4.2. A Posteriori Error Estimates and Grid Adaptation 191

v ’ Qh v ¤C hE |v|1,∆(E) ,

0,E

where the constant C > 0 depends only on the family of triangulations.

Now, setting vh = Qh v in (4.8), the discrete Cauchy-Schwarz inequality

yields

a(u ’ uh , v) ¤ C 0,K |v|1,∆(K)

hK rK (uh )

K∈Th

hE [νE · ∇uh ]E |v|1,∆(E)

+C 0,E

E∈Eh

1/2 1/2

¤C |v|2

h2 rK (uh ) 2

K 0,K 1,∆(K)

K∈Th K∈Th

1/2 1/2

2

hE [νE · ∇uh ]E 0,E |v|2

+C .

1,∆(E)

E∈Eh E∈Eh

A detailed investigation of the two second factors shows that we can

decompose the integrals over ∆(K), ∆(E), according to

... = ... , ... = ... .

∆(K) K ∆(E) K

K ‚∆(K) K ‚∆(E)

This leads to a repeated summation of the integrals over the single elements

K. However, due to the regularity of the family of triangulations, the mul-

tiplicity of these summations is bounded independently of the particular

triangulation (see (3.93)). So we arrive at the estimates

|v|2

1,∆(K) ¤ C|v|1 , |v|2

1,∆(E) ¤ C|v|1 .

2 2

K∈Th E∈Eh

Using the inequality a + b ¤ 2(a2 + b2 ) for a, b ∈ R, we get

a(u ’ uh , v)

1/2

¤ hE [νE · ∇uh ]E |v|1 .

h2 rK (uh ) 2 2

C +

K 0,K 0,E

K∈Th E∈Eh

Finally, (4.6) yields

u ’ uh ¤ D· with · 2 := 2

·K

1

K∈Th

and

1 2

·K := h2 f ’ ruh hE [νE · ∇uh ]E

2 2

+ . (4.9)

K 0,K 0,E

2

E‚‚K\‚„¦

192 4. Grid Generation and A Posteriori Error Estimation

Here we have taken into account that in the transformation of the edge

sum

... into the double sum ...

E∈Eh K∈Th E‚‚K\‚„¦

the latter sums up every interior edge twice.

In summary, we have obtained an a posteriori error estimate of the form

(4.3). By means of re¬ned arguments it is also possible to derive lower

bounds for u ’ uh 1 . For details, we refer to the literature, for example

[35].

Error Estimation by Gradient Recovery

If we are interested in an estimate of the error u ’ uh ∈ V = H0 („¦) 1

measured in the H 1 - or energy norm · , this problem can be simpli¬ed by

means of the fact that both norms are equivalent on V to the H 1 -seminorm

1/2

|u ’ uh |1 = |∇u ’ ∇uh | dx =: ∇u ’ ∇uh

2

.

0

„¦

This is a simple consequence of the de¬nitions and the Poincar´ inequality

e

(see Theorem 2.18). That is, there exist constants C1 , C2 > 0 independent

of h such that

C1 |u ’ uh |1 ¤ u ’ uh ¤ C2 |u ’ uh |1 (4.10)

(cf. Exercise 4.8). Consequently, ∇u remains the only unknown quantity in

the error bound.

The idea of error estimation by means of gradient recovery is to replace

the unknown gradient of the weak solution u by a suitable quantity Rh uh

that is computable from the approximative solution uh at moderate ex-

pense. A popular example of such a technique is the so-called Z 2 estimate.

Here we will describe a simple version of it. Further applications can be

found in the original papers by Zienkiewicz and Zhu, e.g., [74].

Similar to the notation introduced in the preceding subsection, for a

given node a the set

∆(a) := K

K :a∈‚K

denotes the triangular neighbourhood of a (see Figure 4.7). This set co-

incides with the support of the piecewise linear, continuous basis function

associated with that node.

The gradient ∇uh of a piecewise linear, continuous ¬nite element function

uh is constant on every triangle K. This suggests that at any node a of

the triangulation Th we de¬ne the average Rh uh (a) of the values of the

gradients on those triangles sharing the vertex a:

1

∇uh |K |K| .

Rh uh (a) :=

|∆(a)|

K‚∆(a)

4.2. A Posteriori Error Estimates and Grid Adaptation 193

¢d

d

¢ ¢d

¢ ¢d

d

¢¢

a

¢d ¢

dr

¢ e

¢ e ¨ ¨

¢ e¨

Figure 4.7. The triangular neighbourhood ∆(a).

Interpolating the two components of these nodal values of Rh uh separately

in Vh , we get a recovery operator Rh : Vh ’ Vh — Vh .

Now a local error estimator can be de¬ned by the simple restriction of

the quantity · := Rh uh ’ ∇uh 0 onto a single element K:

· K := Rh uh ’ ∇uh .