<< . .

. 25
( : 59)

. . >>

â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

<< . .

. 25
( : 59)

. . >>