â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