way no condition about the geometric quality of the family of triangulations

arises. Of course, these results must be combined with estimates for the

158 3. Finite Element Methods for Linear Elliptic Problems

approximation error of Vh , for which, in particular, both approaches of

Section 3.4 (either regularity or maximum angle condition) are admissible.

For the sake of simplicity, we restrict our attention in the following to the

case of the polynomial ansatz space P = Pk (K). More general results of

similar type, in particular for triangulations with the cuboid element and

ˆ ˆ

P = Qk (K) as reference element, are summarized in [9, p. 207].

We recall the notation and the relations introduced in (3.114), (3.115)

for the local errors. In the following theorems we make use of the Sobolev

spaces W∞ on „¦ and on K with the norms · l,∞ and · l,∞,K , respectively,

l

and the seminorms | · |l,∞ and | · |l,∞,K , respectively. The essential local

assertion is the following:

Theorem 3.41 Suppose k ∈ N and P = Pk (K) and the quadrature

ˆ ˆ

formula is exact for P2k’2 (K):

ˆ

for all v ∈ P2k’2 (K) .

ˆv ˆ

E(ˆ) = 0 ˆ (3.128)

Then there exist some constant C > 0 independent of h > 0 and K ∈ Th

such that for l ∈ {1, k} the following estimates are given:

|EK (apq)| ¤ Chl a

(1) p q

k,∞,K l’1,K 0,K

K

for a ∈ W∞ (K) , p, q ∈ Pk’1 (K) ,

k

|EK (cpq)| ¤ Chl c

(2) p q

k,∞,K l’1,K 1,K

K

for c ∈ W∞ (K) , p ∈ Pk’1 (K) , q ∈ Pk (K) ,

k

|EK (rpq)| ¤ Chl r

(3) p q

k,∞,K l,K 1,K

K

for r ∈ W∞ (K) , p, q ∈ Pk (K) ,

k

|EK (f q)| ¤ Chk f 1/2

(4) k,∞,K vol (K) q 1,K

K

for f ∈ W∞ (K) , q ∈ Pk (K) .

k

The (unnecessarily varied) notation of the coe¬cients already indicates

the ¬eld of application of the respective estimate. The smoothness assump-

tion concerning the coe¬cients in (1)“(3) can be weakened to some extent.

We prove only assertion (1). However, a direct application of this proof to

assertions (2)“(4) leads to a loss of convergence rate (or higher exactness

conditions for the quadrature). Here, quite technical considerations includ-

ing the insertion of projections are necessary, which can be found to some

extent in [9, pp. 201“203]. In the following proof we intensively make use

of the fact that all norms are equivalent on the “¬xed” ¬nite-dimensional

ansatz space Pk (K). The assumption (3.128) is equivalent to the same con-

ˆ

dition on a general element. However, the formulation already indicates an

assumption that is also su¬cient in more general cases.

3.6. Convergence Rate Results in Case of Quadrature and Interpolation 159

Proof of Theorem 3.41, (1): We consider a general element K ∈ Th

and mappings a ∈ W∞ (K), p, q ∈ Pk’1 (K) on it and, moreover, mappings

k

a ∈ W∞ (K), p, q ∈ Pk’1 (K) de¬ned according to (3.75). First, the proof

kˆ ˆ

ˆ ˆˆ

is done for l = k. On the reference element K, for v ∈ W∞ (K) and q ∈

ˆ kˆ

ˆ ˆ

Pk’1 (K)), we have

ˆ

R

ωl (ˆq )(ˆl ) ¤ C v q

v q dˆ ’ ¤C v

ˆ vˆ

E(ˆq ) = ˆˆ x ˆ vˆ b ˆˆ ˆ q

ˆ ,

ˆ ˆ ˆ

∞,K ∞,K ∞,K

ˆ

K l=1

kˆ ˆ

where the continuity of the embedding of W∞ (K) in C(K) is used (see [8,

p. 181]). Therefore, by the equivalence of · ∞,K and · 0,K on Pk’1 (K),

ˆ

ˆ ˆ

it follows that

E(ˆq ) ¤ C v k,∞,K q 0,K .

ˆ vˆ ˆ ˆˆ ˆ

If a ¬xed q ∈ Pk’1 (K) is chosen, then a linear continuous functional G is

ˆ

ˆ

de¬ned on W∞ (K) by v ’ E(ˆq ) that has the following properties:

kˆ ˆ vˆ

ˆ

G ¤ C q ˆ and G(ˆ) = 0 for v ∈ Pk’1 (K) ˆ

ˆ v ˆ

0,K

by virtue of (3.128).

The Bramble“Hilbert lemma (Theorem 3.24) implies

E(ˆq ) ¤ C |ˆ|k,∞,K q 0,K .

ˆ vˆ v ˆˆ ˆ

According to the assertion we now choose

v = ap for a ∈ W k,∞ (K) , p ∈ Pk’1 (K) ,

ˆ ˆ

ˆ ˆˆ ˆ ˆ

and we have to estimate |ˆp|k,∞K (thanks to the Bramble“Hilbert lemma

aˆ ˆ

not ap k,∞,K ). The Leibniz rule for the di¬erentiation of products implies

ˆˆ ˆ

the estimate

k

|ˆp|k,∞,K ¤ C |ˆ|k’j,∞,K |ˆ|j,∞,K .

aˆ a ˆp (3.129)

ˆ ˆ

j=0

ˆ

Here the constant C depends only on k, but not on the domain K.

Since p ∈ Pk’1 (K), the last term of the sum in (3.129) can be omitted.

ˆ

ˆ

Therefore, we have obtained the following estimate holding for a ∈ W∞ (K),

kˆ

ˆ

p, q ∈ Pk’1 (K):

ˆ

ˆˆ

k’1

¤ |ˆ|k’j,∞,K |ˆ|j,∞,K

ˆ a ˆˆ

E(ˆpq ) C a ˆp q

ˆ

ˆ ˆ

0,K

j=0

(3.130)

k’1

¤ |ˆ|k’j,∞,K |ˆ|j,K

C a ˆp ˆ q

ˆ .

ˆ

0,K

j=0

The last estimate uses the equivalence of · ∞ and · 0 on Pk’1 (K). ˆ

ˆ

We suppose that the transformation F of K to the general element K

has, as usual, the linear part B. The ¬rst transformation step yields the

160 3. Finite Element Methods for Linear Elliptic Problems

factor | det(B)| according to (3.115), and for the backtransformation it

follows from Theorem 3.26 and Theorem 3.27 that

C hk’j |a|k’j,∞,K ,

|ˆ|k’j,∞,K ¤

a ˆ K

C hj | det(B)|’1/2 |p|j,K ,

|ˆ|j,K ¤ (3.131)

pˆ K

C | det(B)|’1/2 q

¤

q

ˆ ˆ 0,K

0,K

for 0 ¤ j ¤ k ’ 1. Here a, p, q are the mappings a, p, q (back)transformed

ˆˆˆ

according to (3.75). Substituting these estimates into (3.130) therefore

yields

k’1

EK (apq) ¤ C |a|k’j,∞,K |p|j,K

hk q 0,K

K

j=0

and from this, assertion (1) follows for l = k.

If l = 1, we modify the proof as follows. Again, in (3.130) we estimate

by using the equivalence of norms:

k’1

¤C |ˆ|k’j,∞,K p

ˆ a ˆˆ

E(ˆpq ) a ˆˆ q

ˆ

ˆ ˆ

j,∞,K 0,K

j=0

k’1

¤C |ˆ|k’j,∞,K

a p

ˆ q

ˆ .

ˆ ˆ ˆ

0,K 0,K

j=0

The ¬rst and the third estimates of (3.131) remain applicable; the second

estimate is replaced with the third such that we have

k’1

EK (apq) ¤ C hK |a|k’j,∞,K ,p q

0,K 0,K

j=0

since the lowest hK -power arises for j = k ’ 1. This estimate yields the

2

assertion (1) for l = 1 .

Finally, we can now verify the assumptions of Theorem 3.38 with the

following result:

Theorem 3.42 Consider a family of a¬ne equivalent Lagrange ¬nite el-

ement discretizations in Rd , d ¤ 3, with P = Pk for some k ∈ N as local

ansatz space. Suppose that the family of triangulations is regular or sat-

is¬es the maximum angle condition in the case of triangles with k = 1.

Suppose that the applied quadrature formulas are exact for P2k’2 . Let the

function w satisfying the Dirichlet boundary condition and let the solution

u of the boundary value problem (3.12), (3.18)“(3.20) (with g3 = 0) belong

to H k+1 („¦).

¯

Then there exist some constants C > 0, h > 0 independent of u and w

such that for the ¬nite element approximation uh +wh according to (3.105),

¯

3.6. Convergence Rate Results in Case of Quadrature and Interpolation 161

¯

(3.121), it follows for h ¤ h that

d

u + w ’ (¯h + wh ) ¤Ch |u|k+1 + |w|k+1 +

k

u kij k,∞

1

i,j=1

d

+ ci +r u +w +f .

k,∞ k,∞ k+1 k+1 k,∞

i=1

¯

¯

Proof: According to (3.108), we aim at estimating u+w’(uh + Ih (w)) 1 ,

¯

where uh satis¬es (3.122).

By virtue of Theorem 3.29 or Theorem 3.35 (set formally “3 = …) we

have

w ’ Ih (w) ¤ Chk |w|k+1 .

¯ (3.132)

1

For the bilinear form ah de¬ned in (3.120), it follows from Theorem 3.41

for v, w ∈ Vh and l ∈ {0, k} that

d

a(v, w) ’ ah (v, w) ¤ EK (kij ‚j (v|K )‚i (w|K )) (3.133)

i,j=1

K∈Th

d

+ EK (ci ‚i (v|K )w) + EK (rvw)

i=1

d d

¤C hl kij + ci +r

k,∞,K k,∞,K k,∞,K

K

i,j=1 i=1

K∈Th

—v w

l,K 1,K

d d

¤ Ch l

kij + ci +r

k,∞ k,∞ k,∞

i,j=1 i=1

1/2

— 2

v w 1,

l,K

K∈Th

by estimating the · k,∞,K -norms in terms of norms on the domain „¦ and

then applying the Cauchy“Schwarz inequality with “index” K ∈ Th .

From this we obtain for l = 1 an estimate of the form

|a(v, w) ’ ah (v, w)| ¤ Ch v w

1 1

such that the estimate required in Lemma 3.39 holds (with C(h) = C · h).

Therefore, there exists some ¯ > 0 such that ah is uniformly Vh -elliptic

h

¯ Hence, the estimate (3.126) is applicable, and the ¬rst addend,

for h ¤ h.

the approximation error, behaves as asserted according to Theorem 3.29 or

Theorem 3.35 (again, choose v = Ih (u) for the comparison element).

In order to estimate the consistency error of ah , a comparison element

v ∈ Vh has to be found for which the corresponding part of the norm in

162 3. Finite Element Methods for Linear Elliptic Problems

(3.133) is uniformly bounded. This is satis¬ed for v = Ih (u), since

1/2 1/2

¤ u ’ Ih (u)

2 2

Ih (u) u +

k

k,K k,K

K∈Th K∈Th

¤ + Ch|u|k+1 ¤ u

u k k+1

due to Theorem 3.29 or Theorem 3.35.

Hence, the consistency error in a behaves as asserted according to (3.133),

so that only the consistency error of l has to be investigated: We have

l ’ lh = b ’ bh ’ a(w, ·) + ah (Ih (w), ·) ,

¯

where bh is de¬ned in (3.120).

If v ∈ Vh , then

a(w, v)’ah (Ih (w), v) ¤ a(w, v)’a(Ih (w), v) + a(Ih (w), v)’ah (Ih (w), v) .

¯ ¯ ¯ ¯

For the ¬rst addend the continuity of a implies

a(w, v) ’ a(Ih (w), v) ¤ C w ’ Ih (w)

¯ ¯ v ,

1

1

so that the corresponding consistency error part behaves like w ’ Ih (w) 1 ,

¯

which has already been estimated in (3.132). The second addend just cor-

responds to the estimate used for the consistency error in a (here, the

¯

di¬erence between Ih and Ih is irrelevant), so that the same contribution to

the convergence rate, now with u k+1 replaced by w k+1 , arises. Finally,

Theorem 3.41, (4) yields for v ∈ Vh ,

|b(v) ’ b(vh )| ¤ |EK (f v)| ¤ C hk vol (K)1/2 f v

k,∞,K 1,K

K

K∈Th K∈Th

¤ C hk |„¦|1/2 f v

k,∞ 1

by proceeding as in (3.133). This implies the last part of the asserted

2

estimate.

If the uniform Vh -ellipticity of ah is ensured in a di¬erent way (perhaps

by Lemma 3.40), one can dispense with the smallness assumption about h.

If estimates as given in Theorem 3.41 are also available for other types of

elements, then triangulations consisting of combinations of various elements

can also be considered.

3.7. The Condition Number of Finite Element Matrices 163

3.7 The Condition Number of Finite Element

Matrices

The stability of solution algorithms for linear systems of equations as de-

scribed in Section 2.5 depends on the condition number of the system

matrix (see [28, Chapter 1]). The condition number also plays an impor-

tant role for the convergence behavior of iterative methods, which will be

discussed in Chapter 5. Therefore, in this section we shall estimate the

spectral condition number (see Appendix A.3) of the sti¬ness matrix

A = a(•j , •i ) (3.134)

i,j=1,...,M

and also of the mass matrix (see (7.45))

B= •j , •i , (3.135)

0 i,j=1,...,M

which is of importance for time-dependent problems. Again, we consider a

¬nite element discretization in the general form of Section 3.4 restricted to

Lagrange elements. In order to simplify the notation, we assume the a¬ne

equivalence of all elements. Further we suppose that

• the family (Th )h of triangulations is regular.

We assume that the variational formulation of the boundary value problem

leads to a bilinear form a that is V -elliptic and continuous on V ‚ H 1 („¦).

As a modi¬cation of de¬nition (1.18), let the following norm (which is

also induced by a scalar product) be de¬ned in the ansatz space Vh =

span{•1 , . . . , •M }:

1/2

|v(ai )|

hd 2

v := . (3.136)

0,h K

ai ∈K

K∈Th

Here, a1 , . . . , aM denote the nodes of the degrees of freedom, where in order

to simplify the notation, M instead of M1 is used for the number of degrees

of freedom. The norm properties follow directly from the corresponding

properties of | · |2 except for the de¬niteness. But the de¬niteness follows

from the uniqueness of the interpolation problem in Vh with respect to

degrees of freedom ai .

Theorem 3.43 (1) There exist constants C1 , C2 > 0 independent of h

such that for v ∈ Vh :

¤v ¤ C2 v

C1 v .

0 0,h 0

(2) There exists a constant C > 0 independent of h such that for v ∈ Vh ,

’1

¤C

v min hK v 0.

1

K∈Th

164 3. Finite Element Methods for Linear Elliptic Problems

Proof: As already known from Sections 3.4 and 3.6, the proof is done

locally in K ∈ Th and there transformed to the reference element K by ˆ

means of F (ˆ) = B x + d.

x ˆ

ˆ

Ad (1): All norms are equivalent on the local ansatz space P , thus also

· 0,K and the Euclidean norm in the degrees of freedom. Hence, there

ˆ

exist some C1 , C2 > 0 such that for v ∈ P ,

ˆˆ ˆˆ

1/2

L

¤ |ˆ(ˆi )|2 ¤ C2 v

ˆˆ ˆˆ

C1 v va .

ˆ ˆ

0,K 0,K

i=1

ˆ

Here, a1 , . . . , aL are the degrees of freedom in K. Due to (3.50) we have

ˆ ˆ

vol (K) = vol (K) | det(B)| ,

ˆ

and according to the de¬nition of hK and the regularity of the family (Th )h ,

˜

there exist constants Ci > 0 independent of h such that

C1 hd ¤ C3 ¤ | det(B)| ¤ C2 hd .

˜ ˜ ˜

d

K K K

By the transformation rule we thus obtain for v ∈ PK , the ansatz space on

K, that

1/2

L

1/2

C1 | det(B)|1/2 v ¤ C2 hd |ˆ(ˆi )|2

ˆ ˆ ˜

C1 v = ˆ va

0,K ˆ K

0,K

i=1

1/2

1/2 L

1/2

˜ 1/2 hd |v(ai )|2 |ˆ(ˆi )|2

˜K

= C2 hd

= C2 va

K

ai ∈K i=1

1/2 1/2

C2 | det(B)|’1/2 v

¤ ˜ ˆˆ ˜ ˆ

C2 hd = C2 hd

C2 v ˆ 0,K

K K

0,K

˜ 1/2 ˆ ˜ ’1/2 v

¤ C2 C2 C1 0,K .

This implies assertion (1).

Ad (2): Arguing as before, now using the equivalence of · 1,K andˆ

· 0,K in P , it follows by virtue of (3.86) for v ∈ PK (with the generic

ˆ

ˆ

constant C) that

¤ C h’1 v

¤ C | det(B)|1/2 B ’1 ¤ C B ’1

v v

ˆ v

ˆ

1,K 2 2 0,K 0,K

K

0,K

by Theorem 3.27 and the regularity of (Th )h , and from this, the assertion

2

(2).

In order to make the norm · 0,h comparable with the (weighted)

Euclidean norm we assume in the following:

• There exists a constant CA > 0 independent of h

such that for every node of Th , the number of elements (3.137)

to which this node belongs is bounded by CA .

3.7. Condition Number of Finite Element Matrices 165

This condition is (partly) redundant: For d = 2 and triangular elements,

the condition follows from the uniform lower bound (3.93) for the smallest

angle as an implication of the regularity. Note that the condition need not

be satis¬ed if only the maximum angle condition is required.

In general, if C ∈ RM,M is a matrix with real eigenvalues »1 ¤ · · · ¤

»M and an orthonormal basis of eigenvectors ξ 1 , . . . , ξ M , for instance a

symmetric matrix, then it follows for ξ ∈ RM \ {0} that

ξT Cξ

»1 ¤ T ¤ »M , (3.138)

ξξ

and the bounds are assumed for ξ = ξ 1 and ξ = ξ M .

Theorem 3.44 There exists a constant C > 0 independent of h such that

we have

d

h

κ(B) ¤ C

min hK

K∈Th

for the spectral condition number of the mass matrix B according to (3.135).

Proof: κ(B) = »M /»1 must be determined. For arbitrary ξ ∈ RM \ {0}

we have

2

ξ T Bξ ξT Bξ v 0,h

= ,

v2

ξT ξ ξT ξ

0,h

M

where v := i=1 ξi •i ∈ Vh . By virtue of ξ T Bξ = v, v 0 , the ¬rst factor on

the right-hand side is uniformly bounded from above and below according

to Theorem 3.43. Further, by (3.137) and ξ = (v(a1 ), . . . , v(aM ))T it follows

that

min hd |ξ|2 ¤ v ¤ CA hd |ξ|2 ,

2

K 0,h

K∈Th

and, thus the second factor is estimated from above and below. This leads

to estimates of the type

»1 ≥ C1 min hd , »M ¤ C2 hd ,

K

K∈Th

2

and from this, the assertion follows.

Therefore, if the family of triangulations (Th )h is quasi-uniform in the

sense that there exists a constant C > 0 independent of h such that

h ¤ C hK for all K ∈ Th , (3.139)

then κ(B) is uniformly bounded.

In order to be able to argue analogously for the sti¬ness matrix, we

assume that we stay close to the symmetric case:

166 3. Finite Element Methods for Linear Elliptic Problems

Theorem 3.45 Suppose the sti¬ness matrix A (3.134) admits real eigen-

values and a basis of eigenvectors. Then there exists a constant C > 0

independent of h such that the following estimates for the spectral condition

number κ hold:

’2

’1

A) ¤

κ(B C min hK ,

K∈Th

’2

κ(A) ¤ C min hK κ(B) .

K∈Th

Proof: With the notation of (3.138), we proceed analogously to the proof

of Theorem 3.44. Since

ξ T Aξ ξ T Aξ ξ T Bξ

=T ,

ξT ξ ξ Bξ ξ T ξ

it su¬ces to bound the ¬rst factor on the right-hand side from above and

below. This also yields a result for the eigenvalues of B ’1 A, since we have

for the variable · := B 1/2 ξ,

· T B ’1/2 AB ’1/2 ·

ξ T Aξ

= ,

ξ T Bξ ·T ·

and the matrix B ’1/2 AB ’1/2 possesses the same eigenvalues as B ’1 A by

virtue of B ’1/2 (B ’1/2 AB ’1/2 )B 1/2 = B ’1 A. Here, B 1/2 is the symmet-

ric positive de¬nite matrix that satis¬es B 1/2 B 1/2 = B, and B ’1/2 is its

inverse.

Since ξT Aξ/ξT Bξ = a(v, v)/ v, v 0 and

≥ ≥± v

2 2

a(v, v) ±v ,

1 0

’2 (3.140)

¤ ¤C

2 2

a(v, v) Cv min hK v 0,

1

K∈Th

with a generic constant C > 0 (the last estimate is due to Theorem 3.43,

2), it follows that

’2

ξT Aξ

a(v, v) a(v, v)

±¤ ¤C

=T = min hK , (3.141)

v, v 0 v, v 0

ξ Bξ K∈Th

2

and from this the assertion.

The analysis of the eigenvalues of the model problem in Example 2.12

shows that the above-given estimates are not too pessimistic.

3.8. General Domains and Isoparametric Elements 167

3.8 General Domains and Isoparametric Elements

All elements considered so far are bounded by straight lines or plane sur-

faces. Therefore, only polyhedral domains can be decomposed exactly by

means of a triangulation. Depending on the application, domains with a

curved boundary may appear. With the available elements the obvious way

of dealing with such domains is the following (in the two-dimensional case):

for elements K that are close to the boundary put only the nodes of one

edge on the boundary ‚„¦. This implies an approximation error for the

domain, for „¦h := K∈Th K, there holds in general neither „¦ ‚ „¦h nor

„¦h ‚ „¦ (see Figure 3.14).

B

Figure 3.14. „¦ and „¦h .

As the simplest example, we consider homogeneous Dirichlet boundary

conditions, thus V = H0 („¦), on a convex domain for which therefore „¦h ‚

1

„¦ is satis¬ed. If an ansatz space Vh is introduced as in Section 3.3, then

functions de¬ned on „¦h are generated. Therefore, these functions must be

extended to „¦ in such a way that they vanish on ‚„¦, and consequently, for

the generated function space Vh , Vh ‚ V . This is supposed to be done by

˜˜

adding the domains B whose boundary consists of a boundary part of some

element K ∈ Th close to the boundary and a subset of ‚„¦ to the set of

elements with the ansatz space P (B) = {0}. C´a™s lemma (Theorem 2.17)

e

can still be applied, so that for an error estimate in · 1 the question of

how to choose a comparison element v ∈ Vh arises. The ansatz v = Ih (u),

˜ ˜

˜

where Ih (u) denotes the interpolation on „¦h extended by 0 on the domains

B, is admissible only for the (multi-)linear ansatz: Only in this case are all

nodes of an edge “close to the boundary” located on ‚„¦ and therefore have

homogeneous degrees of freedom, so that the continuity on these edges is

ensured. For the present, let us restrict our attention to this case, so that

u ’ Ih (u) 1 has to be estimated where u is the solution of the boundary

˜

value problem.

The techniques of Section 3.4 can be applied to all K ∈ Th , and by the

conditions assumed there about the triangulation, this yields

u ’ uh ¤ u ’ Ih (u)

C +u

1 1,„¦h 1,„¦\„¦h

¤ C h|u|2,„¦h + u .

1,„¦\„¦h

168 3. Finite Element Methods for Linear Elliptic Problems

If ‚„¦ ∈ C 2 , then we have the estimate

¤ Ch u

u 2,„¦

1,„¦\„¦h

for the new error part due to the approximation of the domain, and thus

the convergence rate is preserved. Already for a quadratic ansatz this is no

longer satis¬ed, where only

u ’ uh ¤ Ch3/2 u

1 3

holds instead of the order O(h2 ) of Theorem 3.29 (see [31, pp. 194 ¬]).

One may expect that this decrease of the approximation quality arises only

locally close to the boundary, however, one may also try to obtain a better

approximation of the domain by using curved elements. Such elements can

ˆ ˆˆ

be de¬ned on the basis of the reference elements (K, P , Σ) of Lagrange type

introduced in Section 3.3 if a general element is obtained from this one by

an isoparametric transformation; that is, choose an

F ∈ (P )d

ˆ (3.142)

that is injective and then

P := p —¦ F ’1 p ∈ P , Σ := F (ˆ) a ∈ Σ .

ˆ ˆˆ aˆˆ

K := F (K) , ˆ

Since the bijectivity of F : K ’ K is ensured by requirement, a ¬nite

ˆ

element is thus de¬ned in terms of (3.58). By virtue of the unique solvability

of the interpolation problem, F can be de¬ned by prescribing a1 , . . . , aL ,

L = |Σ|, and requiring

ˆ

F (ˆi ) = ai ,

a i = 1, . . . , L .

However, this does not in general ensure the injectivity. Since, on the other

hand, in the grid generation process elements are created by de¬ning the

nodes (see Section 4.1), geometric conditions about their positions that

characterize the injectivity of F are desirable. A typical curved element

that can be used for the approximation of the boundary can be generated

on the basis of the unit simplex with P = P2 (K) (see Figure 3.15).

ˆ ˆ

â3 a3

a1 3

â2 3 2

F ∈ P2 (K)

ˆ

â1 3 a2 3

E

a1

a2

â1 â2

â1 2 a1 2

Figure 3.15. Isoparametric element: quadratic ansatz on triangle.

Elements with, in general, one curved edge and otherwise straight edges

thus are suggested for the problem of boundary approximation. They are

3.8. General Domains and Isoparametric Elements 169

combined with a¬ne “quadratic triangles” in the interior of the domain.

Subparametric elements can be generated analogously to the isoparametric

elements if (the components of) the transformations in (3.142) are restricted

to some subspace PT ‚ P . If PT = P1 (K), we again obtain the a¬ne

ˆ ˆ ˆ ˆ

equivalent elements.

However, isoparametric elements are also important if, for instance, the

unit square or cube is supposed to be the reference element. Only the

isoparametric transformation allows for “general” quadrilaterals and hex-

ahedra, respectively, which are preferable in anisotropic cases (for instance

in generalization of Figure 3.11) to simplices due to their adaptability to

ˆ ˆ ˆ

local coordinates. In what follows, let K = [0, 1]d , P = Q1 (K).

In general, since also a ¬nite element (in Rd’1 ) is de¬ned for every face

S of K with P |S and Σ|S , the “faces” of K, that is, F [S], are already

ˆ ˆ ˆˆ ˆˆ ˆ

uniquely de¬ned by the related nodes.

Consequently, if d = 2, the edges of the general quadrilateral are straight

lines (see Figure 3.16), but if d = 3, we have to expect curved surfaces

(hyperbolic paraboloids) for a general hexahedron.

a3

â4 â3 a4

ˆ2

F ∈ Q1 (K)

E

a1 a2

â2

â1

Figure 3.16. Isoparametric element: bilinear ansatz on rectangle.

A geometric characterization of the injectivity of F is still unknown (to

our knowledge) for d = 3, but it can be easily derived for d = 2: Let the

nodes a1 , a2 , a3 , a4 be numbered counterclockwise and suppose that they

are not on a straight line, and thus (by rearranging) T = conv (a1 , a2 , a4 )

forms a triangle such that

2 vol (T ) = det(B) > 0 .

Here FT (ˆ) = B x + d is the a¬ne-linear mapping that maps the refer-

x ˆ

’1

ence triangle conv (ˆ1 , a2 , a4 ) bijectively to T . If a3 := FT (a3 ), then the

aˆˆ ˜

˜

quadrilateral K with the vertices a1 , a2 , a3 , a4 is mapped bijectively to K

ˆˆ˜ˆ

by FT .

The transformation F can be decomposed into

F = FT —¦ FQ ,

2

where FQ ∈ Q1 (K)

ˆ denotes the mapping de¬ned by

FQ (ˆi ) = ai ,

a ˆ i = 1, 2, 4 , FQ (ˆ3 ) = a3

a ˜

170 3. Finite Element Methods for Linear Elliptic Problems

(see Figure 3.17).

x2

^

x2

^

a4

1 1

â4 â3 â4

ã3

a3

FT

FQ ˜

^

K K K

â2

â1 â1 a1

â2

x1 a2

^

x1

^

1

1

F

Figure 3.17. Decomposition of the bilinear isoparametric mapping.

Therefore, the bijectivity of F is equivalent to the bijectivity of FQ .

We characterize a “uniform” bijectivity which is de¬ned by

det (DF (ˆ1 , x2 )) = 0 for the functional matrix DF (ˆ1 , x2 ):

xˆ xˆ

Theorem 3.46 Suppose Q is a quadrilateral with the vertices a1 ,. . .,a4

(numbered counterclockwise). Then,

for all (ˆ1 , x2 ) ∈ [0, 1]2 ⇐’

det (DF (ˆ1 , x2 )) = 0

xˆ xˆ

for all (ˆ1 , x2 ) ∈ [0, 1]2 ⇐’

det (DF (ˆ1 , x2 )) > 0

xˆ xˆ

Q is convex and does not degenerate into a triangle or straight line .

Proof: By virtue of

det (DF (ˆ1 , x2 )) = det(B) det (DFQ (ˆ1 , x2 ))

xˆ xˆ

and det(B) > 0, F can be replaced with FQ in the assertion. Since

a3,1 ’ 1

x1

ˆ ˜

FQ (ˆ1 , x2 ) =

xˆ + x1 x2 ,

ˆˆ

a3,2 ’ 1

x2

ˆ ˜

it follows by some simple calculations that

det (DFQ (ˆ1 , x2 )) = 1 + (˜3,2 ’ 1)ˆ1 + (˜3,1 ’ 1)ˆ2

xˆ a x a x

is an a¬ne-linear mapping because the quadratic parts just cancel each

other. This mapping assumes its extrema on [0, 1]2 at the 4 vertices, where

we have the following values:

(1, 1) : a3,1 + a3,2 ’ 1 .

(0, 0) : 1 , (1, 0) : a3,2 ,

˜ (0, 1) : a3,1 ,

˜ ˜ ˜

A uniform sign is thus obtained if and only if the function is everywhere

positive. This is the case if and only if

a3,1 , a3,2 , a3,1 + a3,2 ’ 1 > 0 ,

˜ ˜ ˜ ˜

˜

which just characterizes the convexity and the nondegeneration of K. By

2

the transformation FT this also holds for K.

3.9. The Maximum Principle for Finite Element Methods 171

According to this theorem it is not allowed that a quadrilateral degener-

ates into a triangle (now with linear ansatz). But a more careful analysis [55]

shows that this does not a¬ect negatively the quality of the approximation.

In general, for isoparametric elements we have the following:

From the point of view of implementation, only slight modi¬cations have

to be made: In the integrals (3.111), (3.111) transformed to the reference

element or their approximation by quadrature (3.120), | det B| has to be

replaced with |det (DF (ˆ))| (in the integrand).

x

The analysis of the order of convergence can be done along the same

lines as in Section 3.4 (and 3.6), however, the transformation rules for the

integrals become more complex (see [9, pp. 237 ¬.]).

3.9 The Maximum Principle for Finite Element

Methods

In this section maximum and comparision principles that have been intro-

duced for the ¬nite di¬erence method are outlined for the ¬nite element

method.

In the case of two-dimensional domains „¦ the situation has been well

investigated for linear elliptic boundary value problems of second order

and linear elements. For higher-dimensional problems (d > 2) as well as

other types of elements, the corresponding assumptions are much more

complex, or there does not necessarily exist any maximum principle.

From now on, let „¦ ‚ R2 be a polygonally bounded domain and let Xh

denote the ¬nite element space of continuous, piecewise linear functions

for a conforming triangulation Th of „¦ where the function values in the

nodes on the Dirichlet boundary “3 are included in the degrees of freedom.

First, we consider the discretization developed for the Poisson equation

’∆u = f with f ∈ L2 („¦). The algebraization of the method is done

according to the scheme described in Section 2.4.3. According to this, ¬rst

all nodes inside „¦ and on “1 and “2 are numbered consecutively from 1

to a number M1 . The nodal values uh (ar ) for r = 1, . . . , M1 are arranged

in the vector uh . Then, the nodes that belong to the Dirichlet boundary

are numbered from M1 + 1 to some number M1 + M2 , the corresponding

ˆ ˆ

nodal values generate the vector uh . The combination of uh and uh gives

uh ∈ RM , M = M + M .

˜

the vector of all nodal values uh = uh 1 2

ˆ

This leads to a linear system of equations of the form (1.31) described

in Section 1.4:

Ah uh = ’Ah uh + f

ˆˆ

with Ah ∈ RM1 ,M1 , Ah ∈ RM1 ,M2 , uh , f ∈ RM1 and uh ∈ RM2 .

ˆ ˆ

Recalling the support properties of the basis functions •i , •j ∈ Xh ,

˜

we obtain for a general element of the (extended) sti¬ness matrix Ah :=

172 3. Finite Element Methods for Linear Elliptic Problems

Ah Ah ∈ RM1 ,M following the relation

ˆ

∇•j · ∇•i dx = ∇•j · ∇•i dx .

˜

(Ah )ij =

supp •i ©supp •j

„¦

Therefore, if i = j, the actual domain of integration consists of at most

two triangles. Hence, for the present it is reasonable to consider only one

triangle as the domain of integration .

Lemma 3.47 Suppose Th is a conforming triangulation of „¦. Then for

an arbitrary triangle K ∈ Th with the vertices ai , aj (i = j), the following

relation holds:

1

∇•j · ∇•i dx = ’ cot ±K , ij

2

K

where ±K denotes the interior angle of K that is opposite to the edge with

ij

the boundary points ai , aj .

Proof: Suppose the triangle K has the vertices ai , aj , ak (see Figure 3.18).

On the edge opposite to the point aj , we have

•j ≡ 0 .

Therefore, ∇•j has the direction of a normal vector to this edge and ” by

considering in which direction •j increases ” the orientation opposite to

the outward normal vector νki , that is,

∇•j = ’ |∇•j | νki |νki | = 1 .

with (3.143)

. aj

ν jk

. hj

±K

ij

ak

.

ν ki ai

Figure 3.18. Notation for the proof of Lemma 3.47.

In order to calculate |∇•j | we use the following: From (3.143) we obtain

|∇•j | = ’∇•j · νki ;

that is, we have to compute a directional derivative. By virtue of •j (aj ) = 1,

we have

0’1 1

∇•j · νki = =’ ,

hj hj

3.9. Maximum Principle for Finite Element Methods 173

where hj denotes the height of K with respect to the edge opposite aj .

Thus we have obtained the relation

1

∇•j = ’ νki .

hj

Hence we have

cos ±K

νki · νjk ij

∇•j · ∇•i = =’ .

hj hi hj hi

Since

2 |K| = hj |ak ’ ai | = hi |aj ’ ak | = |ak ’ ai | |aj ’ ak | sin ±K ,

ij

we obtain

cos ±K 1 1

ij

∇•j · ∇•i = ’ |ak ’ ai | |aj ’ ak | = ’ cot ±K ,

4 |K|2 |K|

ij

2

2

so that the assertion follows by integration.

Corollary 3.48 If K and K are two triangles of Th which have a common

edge spanned by the nodes ai , aj , then

1 sin(±K + ±K )

ij ij

∇•j · ∇•i dx = ’

˜

(Ah )ij = .

K )(sin ±K )

2 (sin ±ij

K∪K ij

Proof: The formula follows from the addition theorem for the cotangent

2

function.

Lemma 3.47 and Corollary 3.48 are the basis for the proof of the as-

˜

sumption (1.32)* in the case of the extended system matrix Ah . Indeed,

additional assumptions about the triangulation Th are necessary:

Angle condition: For any two triangles of Th with a common edge, the

sum of the interior angles opposite to this edge does not exceed the

value π. If a triangle has an edge on the boundary part “1 or “2 , then

the angle opposite this edge must not be obtuse.

Connectivity condition: For every pair of nodes both belonging to „¦ ∪

“1 ∪ “2 there exists a polygonal line between these two nodes such

that the polygonal line consists only of triangle edges whose boundary

points also belong to „¦ ∪ “1 ∪ “2 (see Figure 3.19).

Discussion of assumption (1.32)*: The proof of (1), (2), (5), (6)* is rather

elementary. For the “diagonal elements,”

|∇•r |2 dx = |∇•r |2 dx > 0 ,

(Ah )rr = r = 1, . . . , M1 ,

„¦ K

K‚supp •r

which already is (1). Checking the sign conditions (2) and (5) for the

˜

“nondiagonal elements” of Ah requires the analysis of two cases:

174 3. Finite Element Methods for Linear Elliptic Problems

Figure 3.19. Example of a nonconnected triangulation (“3 = ‚„¦).

(i) For r = 1, . . . , M1 and s = 1, . . . , M with r = s, there exist two

triangles that have the common vertices ar , as .

(ii) There exists only one triangle that has ar as well as as as vertices.

In case (i), Corollary 3.48 can be applied, since if K, K just denote the two

triangles with a common edge spanned by ar , as , then 0 < ±K + ±K ¤ π

rs rs

and thus (Ah )rs ¤ 0, r = s. In case (ii), Lemma 3.47, due to the part of

˜

the angle condition that refers to the boundary triangles, can be applied

directly yielding the assertion.

Further, since M •s = 1 in „¦, we obtain

s=1

M M M

∇•s · ∇•r dx = ∇ · ∇•r dx = 0 .

˜

(Ah )rs = •s

„¦ „¦

s=1 s=1 s=1

This is (6)*.

The sign condition in (3) now follows from (6)* and (5), since we have

M1 M M

(Ah )rs ’ (Ah )rs ≥ 0 .

˜ ˆ

(Ah )rs = (3.144)

s=1 s=1 s=M1 +1

=0

The di¬cult part of the proof of (3) consists in showing that at least one of

these inequalities (3.144) is satis¬ed strictly. This is equivalent to the fact

ˆ

that at least one element (Ah )rs , r = 1, . . . , M1 and s = M1 + 1, . . . , M ,

is negative, which can be shown in terms of an indirect proof by using

Lemma 3.47 and Corollary 3.48, but is not done here in order to save

space. Simultaneously, this also proves the condition (7).

The remaining condition (4)* is proved similarly. First, due to the con-

nectivity condition, the existence of geometric connections between pairs of

nodes by polygonal lines consisting of edges is obvious. It is more di¬cult

to prove that under all possible connections there exists one along which

3.9. Maximum Principle for Finite Element Methods 175

the corresponding matrix elements do not vanish. This can be done by the

same technique of proof as used in the second part of (3), which, however,

is not presented here.

If the angle condition given above is replaced with a stronger angle con-

dition in which stretched and right angles are excluded, then the proof of

(3) and (4)* becomes trivial.

Recalling the relations

max uh (x) = max (˜ h )r

u

r∈{1,...,M }

x∈„¦

and

max uh (x) = max (ˆ h )r ,

u

x∈“3 r∈{M1 +1,...,M}

which hold for linear elements, the following result can be derived from

Theorem 1.10.

Theorem 3.49 If the triangulation Th satis¬es the angle condition and

the connectivity condition, then we have the following estimate for the ¬nite

element solution uh of the Poisson equation in the space of linear elements

for a nonpositive right-hand side f ∈ L2 („¦):

max uh (x) ¤ max uh (x) .

x∈“3

x∈„¦

Finally, we make two remarks concerning the case of more general

di¬erential equations.

If an equation with a variable scalar di¬usion coe¬cient k : „¦ ’ R is con-

sidered instead of the Poisson equation, then the relation in Corollary 3.48

loses its purely geometric character. Even if the di¬usion coe¬cient is

supposed to be elementwise constant, the data-dependent relation

1

(Ah )ij = ’

˜ kK cot ±K + kK cot ±K

ij ij

2

would arise, where kK and kK denote the constant restriction of k to the

triangles K and K , respectively. The case of matrix-valued coe¬cients

K : „¦ ’ Rd,d is even more problematic.

The second remark concerns di¬erential expressions that also contain

lower-order terms, that is, convective and reactive parts. If the di¬usive

term ’∇ · (K∇u) can be discretized in such a way that a maximum

principle holds, then this maximum principle is preserved if the discretiza-

tion of the other terms leads to matrices whose “diagonal elements” are

nonnegative and whose “nondiagonal elements” are nonpositive. These ma-

trix properties are much simpler than the conditions (1.32) and (1.32)*.

However, satisfying these properties causes di¬culties in special cases,

e.g., for convection-dominated equations (see Chapter 9), unless additional

restrictive assumptions are made or special discretization schemes are used.

4

Grid Generation and A Posteriori

Error Estimation

4.1 Grid Generation

As one of the ¬rst steps, the implementation of the ¬nite element method

(and also of the ¬nite volume method as described in Chapter 6) requires

a “geometric discretization” of the domain „¦.

This part of a ¬nite element program is usually included in the so-called

preprocessor (see also Section 2.4.1). In general, a ¬nite element program

consists further of the intrinsic kernel (assembling of the ¬nite-dimensional

system of algebraic equations, rearrangement of data (if necessary), solution

of the algebraic problem) and the postprocessor (editing of the results, ex-

traction of intermediate results, preparation for graphic output, a posteriori

error estimation).

4.1.1 Classi¬cation of Grids

Grids can be grouped according to di¬erent criteria: One criterion considers

the geometric shape of the elements (triangles, quadrilaterals, tetrahedra,

hexahedra, prisms, pyramids; possibly with curved boundaries). A further

criterion distinguishes the logical structure of the grid (structured or un-

structured grids). Beside these rough classes, in practice one can ¬nd a large

number of variants combining grids of di¬erent classes (combined grids).

A structured grid in the strict sense is characterized by a regular arrange-

ment of the grid points (nodes), that is, the connectivity pattern between

neighbouring nodes is identical everywhere in the interior of the grid. The

4.1. Grid Generation 177

only exceptions of that pattern may occur near the boundary of the domain

„¦.

Typical examples of structured grids are rectangular Cartesian two- or

three-dimensional grids as they are also used within the framework of the

¬nite di¬erence methods described in Chapter 1 (see, e.g., Figure 1.1).

A structured grid in the wider sense is obtained by the application of a

piecewise smooth bijective transformation to some “reference grid”, which

is a structured grid in the strict sense. Grids of this type are also called

logically structured, because only the logical structure of the connectivity

pattern is ¬xed in the interior of the grid. However, the edges or faces of

the geometric elements of a logically structured grid are not necessarily

straight or even.

Logically structured grids have the advantage of simple implementation,

because the pattern already de¬nes the neighbours of a given node. Fur-

thermore, there exist e¬cient methods for the solution of the algebraic

system resulting from the discretization, including parallelized resolution

algorithms.

In contrast to structured grids, unstructured grids do not have a self-

repeating node pattern. Moreover, elements of di¬erent geometric type can

be combined in unstructured grids.

Unstructured grids are suitable tools for the modelling of complex ge-

ometries of „¦ and for the adjustment of the grid to the numerical solution

(local grid adaptation).

In the subsequent sections, a survey of a few methods for generating

unstructured grids will be given. Methods to produce structured grids can

be found, for instance, in the books [23] or [33].

4.1.2 Generation of Simplicial Grids

A simplicial grid consists of triangles (in two dimensions) or tetrahedra (in

three dimensions). To generate simplicial grids, the following three types

of methods are widely used:

• overlay methods,

• Delaunay triangulations,

• advancing front methods.

Overlay Methods

The methods of this type start with a structured grid (the overlay grid)

that covers the whole domain. After that, this basic grid is modi¬ed near

the boundary to ¬t to the domain geometry. The so-called quadtree (in

two dimensions) or octree technique (in three dimensions) forms a typical

example of an overlay method, where the overlay grid is a relatively coarse

rectangular Cartesian two- or three-dimensional grid. The substantial part

178 4. Grid Generation and A Posteriori Error Estimation

of the algorithm consists of ¬tting routines for those parts of the starting

grid that are located near the boundary and of simplicial subdivisions of

the obtained geometric elements. The ¬tting procedures perform recursive

subdivisions of the boundary rectangles or rectangular parallelepipeds in

such a way that at the end every geometric element contains at most one

geometry de¬ning point (i.e., a vertex of „¦ or a point of ‚„¦, where the

type of boundary conditions changes). Finally, the so-called smoothing step,

which optimizes the grid with respect to a certain regularity criterion, can

be supplemented; see Section 4.1.4.

Typically, grids generated by overlay methods are close to structured

grids in the interior of the domain. Near the boundary, they lose the

structure. Further details can be found in the references [68] and [72].

Delaunay Triangulations

The core algorithm of these methods generates, for a given cloud of isolated

points (nodes), a triangulation of their convex hull. Therefore, a grid gen-

erator based on this principle has to include a procedure for the generation

of this point set (for example, the points resulting from an overlay method)

as well as certain ¬tting procedures (to cover, for example, nonconvex

domains, too).

The Delaunay triangulation of the convex hull of a given point set in

R is characterized by the following property (empty sphere criterion, Fig-

d

ure 4.1): Any open d-ball, the boundary of which contains d + 1 points

from the given set, does not contain any other points from that set. The

triangulation can be generated from the so-called Voronoi tesselation of Rd

for the given point set. In two dimensions, this procedure is described in

Chapter 6, which deals with ¬nite volume methods (Section 6.2.1). How-

Figure 4.1. Empty sphere criterion in two dimensions (d = 2).

4.1. Grid Generation 179

ever, practical algorithms ([48] or [71]) apply the empty sphere criterion

more directly.

The interesting theoretical properties of Delaunay triangulations are one

of the reasons for the “popularity” of this method. In two dimensions, the

so-called max-min-angle property is valid: Among all triangulations of the

convex hull G of a given point set, the Delaunay triangulation maximizes

the minimal interior angle over all triangles. In the case d = 3, this nice

property does not remain true. In contrast, even badly shaped elements (the

so-called sliver elements) may occur. A further important property of a two-

dimensional Delaunay triangulation is that the sum of two angles opposite

an interior edge is not more than π. For example, such a requirement is a

part of the angle condition formulated in Section 3.9.

Advancing Front Methods

The idea of these methods, which are also known in the literature (see, e.g.,

[50], [56], [60], [62]) as moving front methods, is to generate a triangulation

recursively from a discretization of the current boundary. The methods

start with a partition of the boundary of G0 := „¦. For d = 2, this “initial

front” is a polygonal line, whereas in d = 3 it is a triangulation of a curved

surface (the so-called “2.5-dimensional triangulation”). The method con-

sists of an iteration of the following general step (Figure 4.2): An element

of the current front (i.e., a straight-line segment or a triangle) is taken

and then, either generating a new inner point or taking an already existing

point, a new simplex Kj that belongs to Gj’1 is de¬ned. After the data of

the new simplex are saved, the simplex is deleted from Gj’1 . In this way,

a smaller domain Gj with a new boundary ‚Gj (a new “current front”)

results. The general step is repeated until the current front is empty. Of-

ten, the grid generation process is supplemented by the so-called smoothing

step; see Section 4.1.4.

¢d d

d d

;

¢ Kj d d

Gj’1 Gj

¢

¡ ¡

¢

d d

¡ ¡

d d

d d

Figure 4.2. Step j of the advancing front method: The new simplex Kj is deleted

from the domain Gj’1 .

180 4. Grid Generation and A Posteriori Error Estimation

4.1.3 Generation of Quadrilateral and Hexahedral Grids

Grids consisting of quadrilaterals or hexahedra can also be generated by

means of overlay methods (e.g., [66]) or advancing front methods (e.g.,

[46], [47]). An interesting application of simplicial advancing front meth-

ods in the two-dimensional case is given in the paper [73]. The method is

based on the simple fact that any two triangles sharing a common edge

form a quadrilateral. Obviously, a necessary condition for the success of

the method is that the triangulation should consist of an even number of

triangles. Unfortunately, the generalization of the method to the three-

dimensional situation is di¬cult, because a comparatively large number of

adjacent tetrahedra should be united to form a hexahedron.

Multiblock Methods

The basic idea of these methods is to partition the domain into a small num-

ber of subdomains (“blocks”) of simple shape (quadrilaterals, hexahedra,

as well as triangles, tetrahedra, prisms, pyramids, etc.) and then generate

structured or logically structured grids in the individual subdomains (see,

e.g., [23], [33]).

In multiblock grids, special attention has to be devoted to the treatment

of common boundaries of adjacent blocks. Unless special discretization

methods such as, for example, the so-called mortar ¬nite element method

(cf. [45]) are used in this situation, there may be a con¬‚ict between certain

compatibility conditions at the common block interfaces (to ensure, e.g.,

the continuity of the ¬nite element functions across the interfaces) on the

one hand and the output directives of an error estimation procedure that

may advise to re¬ne a block-internal grid locally on the other hand.

Hierarchically Structured Grids

These grids are a further, hybrid variant of structured and unstructured

grids, though not yet very widespread. Starting with a logically structured

grid, hierarchically structured grids are generated by a further logically

structured re¬nement of certain subdomains. As in multiblock methods,

the interfaces between blocks of di¬erent re¬nement degrees have to be

treated carefully.

Combined Grids

Especially in three-dimensional situations, the generation of “purely” hexa-

hedral grids may be very di¬cult for complicated geometries of the domain.

Therefore, the so-called combined grids that consist of hexahedral grids in

geometrically simple subdomains and tetrahedral, prismatic, pyramidal,

etc. grids in more critical subregions are used.

Chimera Grids

These grids are also called overset grids (see, e.g., [51]). In contrast to the

multiblock grids described above, here the domain is covered by a compar-

4.1. Grid Generation 181

atively small number of domains of simple shape, and then structured or

logically structured grids are generated on the individual domains. That is,

a certain overlapping of the blocks and thus of the subgrids is admitted.

4.1.4 Grid Optimization

Many grid generation codes include “smoothing algorithms” that optimize

the grid with respect to certain regularity criteria. In the so-called r-method

(relocation method) the nodes are slightly moved, keeping the logical struc-

ture (connectivities) of the grid ¬xed. Another approach is to improve the

grid connectivities themselves.

A typical example for r-methods is given by the so-called Laplacian

smoothing (or barycentric smoothing), where any inner grid point is moved

into the barycentre of its neighbours (see [50]). A local weighting of selected

neighbours can also be used (weighted barycentric smoothing). From a for-

mal point of view, the application of the Laplacian smoothing corresponds

to the solution of a system of linear algebraic equations that is obtained

from the equations of the arithmetic (or weighted) average of the nodes.

The matrix of this system is large but sparse. The structure of this matrix

is very similar to the one that results from a ¬nite volume discretization

of the Poisson equation as described in Section 6.2 (see the correspond-

ing special case of (6.9)). In general, there is no need to solve this system

exactly. Typically, only one to three steps of a simple iterative solver (as

presented in Section 5.1) are performed. When the domain is almost con-

vex, Laplacian smoothing will produce good results. It is also clear that for

strongly nonconvex domains or other special situations, the method may

produce invalid grids.

Among the methods to optimize the grid connectivities, the so-called

2:1-rule and, in the two-dimensional case, the edge swap (or diagonal swap,

[59]) are well known. The 2:1-rule is used within the quadtree or octree

method to reduce the di¬erence of the re¬nement levels between neighbour-

ing quadrilaterals or hexahedra to one by means of additional re¬nement

steps; see Figure 4.3.

In the edge swap method, a triangular grid is improved. Since any two

triangles sharing an edge form a convex quadrilateral, the method decides

which of the two diagonals of the quadrilateral optimizes a given criterion.

If the optimal diagonal does not coincide with the common edge, the other

con¬guration will be taken; i.e., the edge will be swapped.

Finally, it should be mentioned that there exist grid optimization

methods that delete nodes or even complete elements from the grid.

4.1.5 Grid Re¬nement

A typical grid re¬nement algorithm for a triangular grid, the so-called

red/green re¬nement, has previously been introduced in Section 2.4.1. A

182 4. Grid Generation and A Posteriori Error Estimation

dq dq

q q

e e

;

eq eq

d d

dqr dqr

r r

rq rq

t t

t t

Figure 4.3. 2:1-rule.

further class of methods is based on bisection, that is, a triangle is divided

by the median of an edge. A method of bisection is characterized by the

number of bisections used within one re¬nement step (stage number of the

method of bisection) and by the criterion of how to select the edge where

the new node is to be located. A popular strategy is to take the longest of

the three edges. The general (recursive) re¬nement step for some triangle

K is of the following form:

(i) Find the longest edge of K and insert the median connecting the

midpoint of that edge with the opposite vertex.

(ii) If the resulting new node is not a vertex of an already existing triangle

or is not a boundary point of the domain „¦, then the adjacent triangle

that shares the re¬ned edge has to be divided, too.

Since the longest edge of the adjacent triangle need not coincide with the

common edge, the application of this scheme leads to a nonconforming

triangulation, in general. To obtain a conforming triangulation, all new

nodes resulting from substep (i) have to be detected, and then certain

closure rules have to be applied.

The red/green re¬nement as well as the method of bisection can be gener-

alized to the three-dimensional case. However, since the number of di¬erent

con¬gurations is signi¬cantly larger than in the case d = 2, only a few

illustrative examples will be given.

The red/green re¬nement of a tetrahedron K (see Figure 4.4) yields a

partition of K into eight subtetrahedra with the following properties: All

vertices of the subtetrahedra coincide either with vertices or with edge mid-

points of K. At all the faces of K, the two-dimensional red/green re¬nement

scheme can be observed.

In addition to the di¬culties arising in the two-dimensional situation,

the (one-stage) bisection applied to the longest edge of a tetrahedron also

yields faces that violate the conformity conditions. Therefore, the closure

rules are rather complicated, and in practice, multistage (often three-stage)

Exercises 183

Figure 4.4. Representation of the red/green re¬nement of a tetrahedron.

methods of bisection are used to circumvent these di¬culties (see Figure

4.5).