. 7
( 16)


both transformation steps cancel each other (which will happen), in this
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,

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

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

|EK (cpq)| ¤ Chl c
(2) p q
k,∞,K l’1,K 1,K

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

|EK (rpq)| ¤ Chl r
(3) p q
k,∞,K l,K 1,K

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

|EK (f q)| ¤ Chk f 1/2
(4) k,∞,K vol (K) q 1,K

for f ∈ W∞ (K) , q ∈ Pk (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

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
ω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 ˆ

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
|ˆp|k,∞,K ¤ C |ˆ|k’j,∞,K |ˆ|j,∞,K .
aˆ a ˆp (3.129)
ˆ ˆ

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),

p, q ∈ Pk’1 (K):
¤ |ˆ|k’j,∞,K |ˆ|j,∞,K
ˆ a ˆˆ
E(ˆpq ) C a ˆp q
ˆ ˆ
¤ |ˆ|k’j,∞,K |ˆ|j,K
C a ˆp ˆ q
ˆ .

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
ˆ ˆ 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
EK (apq) ¤ C |a|k’j,∞,K |p|j,K
hk q 0,K

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:
¤C |ˆ|k’j,∞,K p
ˆ a ˆˆ
E(ˆpq ) a ˆˆ q
ˆ ˆ
j,∞,K 0,K
¤C |ˆ|k’j,∞,K
a p
ˆ q
ˆ .
ˆ ˆ ˆ
0,K 0,K

The ¬rst and the third estimates of (3.131) remain applicable; the second
estimate is replaced with the third such that we have
EK (apq) ¤ C hK |a|k’j,∞,K ,p q
0,K 0,K

since the lowest hK -power arises for j = k ’ 1. This estimate yields the
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
u + w ’ (¯h + wh ) ¤Ch |u|k+1 + |w|k+1 +
u kij k,∞
+ ci +r u +w +f .
k,∞ k,∞ k+1 k+1 k,∞

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
w ’ Ih (w) ¤ Chk |w|k+1 .
¯ (3.132)

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
a(v, w) ’ ah (v, w) ¤ EK (kij ‚j (v|K )‚i (w|K )) (3.133)
+ EK (ci ‚i (v|K )w) + EK (rvw)
d d
¤C hl kij + ci +r
k,∞,K k,∞,K k,∞,K
i,j=1 i=1
—v w
l,K 1,K
d d
¤ Ch l
kij + ci +r
k,∞ k,∞ k,∞
i,j=1 i=1

— 2
v w 1,

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
¯ 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∈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 ,

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∈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

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
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)

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 }:

|v(ai )|
hd 2
v := . (3.136)
0,h K
ai ∈K

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 ,
v min hK v 0.
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 ,
ˆˆ ˆˆ
¤ |ˆ(ˆi )|2 ¤ C2 v
ˆˆ ˆˆ
C1 v va .
ˆ ˆ
0,K 0,K

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 .
˜ ˜ ˜

By the transformation rule we thus obtain for v ∈ PK , the ansatz space on
K, that
C1 | det(B)|1/2 v ¤ C2 hd |ˆ(ˆi )|2
ˆ ˆ ˜
C1 v = ˆ va
0,K ˆ K
1/2 L
˜ 1/2 hd |v(ai )|2 |ˆ(ˆi )|2
= C2 hd
= C2 va
ai ∈K i=1
1/2 1/2
C2 | det(B)|’1/2 v
¤ ˜ ˆˆ ˜ ˆ
C2 hd = C2 hd
C2 v ˆ 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

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

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
κ(B) ¤ C
min hK

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
ξ T Bξ ξT Bξ v 0,h
= ,
ξT ξ ξT ξ
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
min hd |ξ|2 ¤ v ¤ CA hd |ξ|2 ,
K 0,h

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 ,

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:
A) ¤
κ(B C min hK ,
κ(A) ¤ C min hK κ(B) .

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
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,

with a generic constant C > 0 (the last estimate is due to Theorem 3.43,
2), it follows that
ξT Aξ
a(v, v) a(v, v)
±¤ ¤C
=T = min hK , (3.141)
v, v 0 v, v 0
ξ Bξ K∈Th

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).


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 ‚

„¦ 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)
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 .
168 3. Finite Element Methods for Linear Elliptic Problems

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

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
â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.

â4 â3 a4
F ∈ Q1 (K)

a1 a2

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 ˆ
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 ,
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).


1 1
â4 â3 â4
FQ ˜
â1 â1 a1
x1 a2

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
ˆ ˜
FQ (ˆ1 , x2 ) =
xˆ + x1 x2 ,
a3,2 ’ 1
ˆ ˜
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
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).
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
In this section maximum and comparision principles that have been intro-
duced for the ¬nite di¬erence method are outlined for the ¬nite element
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:
∇•j · ∇•i dx = ’ cot ±K , ij

where ±K denotes the interior angle of K that is opposite to the edge with
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

ν 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
∇•j = ’ νki .
Hence we have
cos ±K
νki · νjk ij
∇•j · ∇•i = =’ .
hj hi hj hi
2 |K| = hj |ak ’ ai | = hi |aj ’ ak | = |ak ’ ai | |aj ’ ak | sin ±K ,

we obtain
cos ±K 1 1
∇•j · ∇•i = ’ |ak ’ ai | |aj ’ ak | = ’ cot ±K ,
4 |K|2 |K|
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

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 · ∇•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


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
r∈{1,...,M }

max uh (x) = max (ˆ h )r ,
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) .

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
(Ah )ij = ’
˜ kK cot ±K + kK cot ±K
ij ij
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.
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
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-

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


. 7
( 16)