2 2 2

±u (2.46)

1 a 1

72 2. Finite Element Method for Poisson Equation

i.e., the norms · 1 and · a are equivalent on V = H0 („¦) and therefore

1

they generate the same convergence concept:

uh ’ u with respect to · ” uh ’ u ’0

1 1

” uh ’ u ’ 0 ” uh ’ u with respect to · .

a a

In summary the estimate (2.45) holds for · = · 1 with the constant

1/±.

Because of the Cauchy“Schwarz inequality for the scalar product on

2

L („¦) and

b(v) = f (x)v(x) dx ,

„¦

i.e., |b(v)| ¤ f 0 v 0 ¤ f 0 v 1 , and thus b ¤ f 0 , the stability

estimate (2.44) for a right-hand side f ∈ L2 („¦) takes the particular form

1

¤

uh f 0.

1

±

Up to now, our considerations have been independent of the special form

of Vh . Now we make use of the choice of Vh according to (2.27). In order

to obtain an estimate of the approximation error of Vh , it is su¬cient to

estimate the term u ’ v for some special element v ∈ Vh . For this element

¯ ¯

v ∈ Vh , we choose the interpolant Ih (u), where

¯

Ih : u ∈ C(„¦) u = 0 on ‚„¦ ’ Vh ,

¯

(2.47)

u ’ Ih (u) with Ih (u)(ai ) = u(ai ) .

This interpolant exists and is unique (Lemma 2.10). Obviously,

min u ’ v 1 v ∈ Vh ¤ u ’ Ih (u) 1 for u ∈ C(„¦) and u = 0 on ‚„¦ .

¯

If the weak solution u possesses weak derivatives of second order, then for

¯ ¯

certain su¬ciently ¬ne triangulations Th , i.e., 0 < h ¤ h for some h > 0,

an estimate of the type

u ’ Ih (u) ¤ Ch (2.48)

1

holds, where C depends on u but is independent of h (cf. (3.88)). The

proof of this estimate will be explained in Section 3.4, where also su¬cient

conditions on the family of triangulations (Th )h will be speci¬ed.

Exercises

1

x2 u v dx for arbitrary u, v ∈ H0 (0, 1).

1

2.4 Let a(u, v) := 0

(a) Show that there is no constant C1 > 0 such that the inequality

1

2

a(u, u) ≥ C1 for all u ∈ H0 (0, 1)

1

(u ) dx

0

2.3. Stability and Convergence 73

is valid.

N

(b) Now let Th := {(xi’1 , xi )}i=1 , N ∈ N, be an equidistant partition of

(0, 1) with the parameter h = 1/N and Vh := span {•i }N ’1 , where

i=1

±

(x ’ xi’1 )/h in (xi’1 , xi ) ,

(xi+1 ’ x)/h in (xi , xi+1 ) ,

•i (x) :=

0 otherwise .

Does there exist a constant C2 > 0 with

1

2

a(uh , uh ) ≥ C2 for all uh ∈ Vh ?

(uh ) dx

0

2.5

(a) For „¦ := (±, β) — (γ, δ) and V according to (2.7), prove the inequality

of Poincar´: There exists a positive constant C with

e

¤C u for all u ∈ V .

u 0 a

x

Hint: Start with the relation u(x, y) = ‚x u(s, y) ds .

±

(b) For „¦ := (±, β) and v ∈ C([±, β]) with a piecewise continuous

derivative v and v(γ) = 0 for some γ ∈ [±, β], show that

¤ (β ’ ±) v

v .

0 0

2.6 Let „¦ := (0, 1)—(0, 1). Given f ∈ C(„¦), discretize the boundary value

problem ’∆u = f in „¦, u = 0 on ‚„¦, by means of the usual ¬ve-point

di¬erence stencil as well as by means of the ¬nite element method with

linear elements. A quadratic grid as well as the corresponding Friedrichs“

Keller triangulation will be used.

Prove the following stability estimates for the matrix of the linear system

of equations:

1 1

(a) A’1 , (b) A’1 2 ¤ , (c) A’1 0 ¤ 1 ,

¤

∞

h h h

8 16

where · ∞ , · 2 denote the maximum row sum norm and the spectral

norm of a matrix, respectively, and A’1 0 := supvh ∈Vh vh 2 / vh 2 with

0 a

h

vh 2 := „¦ |∇vh |2 dx.

a

Comment: The constant in (c) is not optimal.

2.7 Let „¦ be a domain with polygonal boundary and let Th be a conform-

ing triangulation of „¦. The nodes ai of the triangulation are enumerated

from 1 to M.

Let the triangulation satisfy the following assumption: There exist

constants C1 , C2 > 0 such that for all triangles K ∈ Th the relation

C1 h2 ¤ vol (K) ¤ C2 h2

74 2. Finite Element Method for Poisson Equation

is satis¬ed. h denotes the maximum of the diameters of all elements of Th .

(a) Show the equivalence of the following norms for uh ∈ Vh in the space

Vh of continuous, piecewise linear functions over „¦ :

M

1/2 1/2

|uh | dx

2

u2 (ai )

uh := , uh := h .

0 0,h h

„¦ i=1

(b) Consider the special case „¦ := (0, 1)—(0, 1) with the Friedrichs“Keller

triangulation as well as the subspace Vh © H0 („¦) and ¬nd “as good

1

as possible” constants in the corresponding equivalence estimate.

2.4 The Implementation of the Finite Element

Method: Part 1

In this section we will consider some aspects of the implementation of

the ¬nite element method using linear ansatz functions on triangles for

the model boundary value problem (1.1), (1.2) on a polygonally bounded

domain „¦ ‚ R2 . The case of inhomogeneous Dirichlet boundary conditions

will be treated also to a certain extent as far as it is possible up to now.

2.4.1 Preprocessor

The main task of the preprocessor is to determine the triangulation.

An input ¬le might have the following format:

Let the number of variables (including also the boundary nodes for

Dirichlet boundary conditions) be M. We generate the following list:

x-coordinate of node 1 y-coordinate of node 1

... ...

x-coordinate of node M y-coordinate of node M

Let the number of (triangular) elements be N. These elements will be

listed in the element-node table. Here, every element is characterized by the

indices of the nodes corresponding to this element in a well-de¬ned order

(e.g., counterclockwise); cf. Figure 2.11.

11

7 10

4

Figure 2.11. Element no. 10 with nodes nos. 4, 11, 7.

2.4. The Implementation of the Finite Element Method: Part 1 75

For example, the 10th row of the element-node table contains the entry

4 11 7

Usually, a triangulation is generated by a triangulation algorithm. A short

overview on methods for the grid generation will be given in Section 4.1.

One of the simplest versions of a grid generation algorithm has the following

structure (cf. Figure 2.12):

Figure 2.12. Re¬nement by quartering.

Prescribe a coarse triangulation (according to the above format) and

re¬ne this triangulation (repeatedly) by subdividing a triangle into 4 con-

gruent triangles by connecting the midpoints of the edges with straight

lines.

If this uniform re¬nement is done globally, i.e., for all triangles of the

coarse grid, then triangles are created that have the same interior angles as

the elements of the coarse triangulation. Thus the quality of the triangu-

lation, indicated, for example, by the ratios of the diameters of an element

and of its inscribed circle (see De¬nition 3.28), does not change. However,

if the subdivision is performed only locally, the resulting triangulation is

no longer admissible, in general. Such an inadmissible triangulation can be

corrected by bisection of the corresponding neighbouring (unre¬ned) tri-

angles. But this implies that some of the interior angles are bisected and

consequently, the quality of the triangulation becomes poorer if the bisec-

tion step is performed too frequently. The following algorithm circumvents

the depicted problem. It is due to R. Bank and is implemented, for example,

in the PLTMG code (see [4]).

76 2. Finite Element Method for Poisson Equation

A Possible Re¬nement Algorithm

Let a (uniform) triangulation T be given (e.g., by repeated uniform re¬ne-

ment of a coarse triangulation). The edges of this triangulation are called

red edges.

(1) Subdivide the edges according to a certain local re¬nement criterion

(introduction of new nodes) by successive bisection (cf. Figure 2.13).

. . . .

. .

.

. . .

. . . .

Figure 2.13. New nodes on edges.

(2) If a triangle K ∈ T has on its edges in addition to the vertices two

or more nodes, then subdivide K into four congruent triangles.

Iterate over step 2 (cf. Figure 2.14).

(3) Subdivide the triangles with nodes at the midpoints of the edges into

2 triangles by bisection. This step introduces the so-called green edges.

(4) If the re¬nement is to be continued, ¬rst remove the green edges.

2.4.2 Assembling

Denote by •1 , . . . , •M the global basis functions. Then the sti¬ness matrix

Ah has the following entries:

N

(m)

∇•j · ∇•i dx =

(Ah )ij = Aij

„¦ m=1

with

(m)

∇•j · ∇•i dx .

Aij =

Km

Let a1 , . . . , aM denote the nodes of the triangulation. Because of the

implication

(m)

= 0 ’ a i , aj ∈ K m

Aij

(m)

(cf. (2.37)), the element Km yields nonzero contributions for Aij only if

ai , aj ∈ Km at best. Such nonzero contributions are called element entries

of Ah . They add up to the entries of Ah .

2.4. The Implementation of the Finite Element Method: Part 1 77

.

. ..

.

.

.

. . ..

.

. ..

.

.

.

. . ..

.

. ..

.

.

.

. . ..

: green edges

Figure 2.14. Two re¬nement sequences.

In Example 2.12 we explained a node-based assembling of the sti¬ness

matrix. In contrast to this and on the basis of the above observations, in

the following we will perform an element-based assembling of the sti¬ness

matrix.

To assemble the entries of A(m) , we will start from a local numbering

(cf. Figure 2.15) of the nodes by assigning the local numbers 1, 2, 3 to the

global node numbers r1 , r2 , r3 (numbered counterclockwise). In contrast to

the usual notation adopted in this book, here indices of vectors according to

the local numbering are included in parentheses and written as superscripts.

78 2. Finite Element Method for Poisson Equation

r3 3

Km

1

r1

r2 2

Figure 2.15. Global (left) and local numbering.

Thus in fact, we generate

˜(m)

A(m) as Aij .

ri rj

i,j=1,2,3 i,j=1,2,3

To do this, we ¬rst perform a transformation of Km onto some reference

element and then we evaluate the integral on this element exactly.

Hence the entry of the element sti¬ness matrix reads as

˜(m) ∇•rj · ∇•ri dx .

Aij =

Km

ˆ

The reference element K is transformed onto the global element Km by

means of the relation F (ˆ) = B x + d, therefore

x ˆ

Dx u(F (ˆ)) = Dx u(F (ˆ)) Dx F (ˆ) = Dx u(F (ˆ)) B ,

x x x x

ˆ ˆ

where Dx u denotes the row vector (‚1 u, ‚2 u), i.e., the corresponding dif-

ferential operator. Using the more standard notation in terms of gradients

and taking into consideration the relation B ’T := (B ’1 )T , we obtain

∇x u (F (ˆ)) = B ’T ∇x (u (F (ˆ)))

x x (2.49)

ˆ

and thus

˜(m) ∇x •rj (F (ˆ)) · ∇x •ri (F (ˆ)) |det(DF (ˆ))| dˆ

Aij = x x x x

ˆ

K

B ’T ∇x •rj (F (ˆ)) · B ’T ∇x (•ri (F (ˆ))) |det(B)| dˆ

= x x x

ˆ ˆ

ˆ

K

B ’T ∇x •rj (ˆ) · B ’T ∇x •ri (ˆ) |det(B)| dˆ

= ˆˆ x ˆˆ x x (2.50)

ˆ

K

B ’T ∇x Nj (ˆ) · B ’T ∇x Ni (ˆ) |det(B)| dˆ ,

= x x x

ˆ ˆ

ˆ

K

where the transformed basis functions •ri , •(ˆ) := •(F (ˆ)) coincide with

ˆ ˆx x

ˆ i.e., with the shape functions Ni :

the local basis functions on K,

•ri (ˆ) = Ni (ˆ) for x ∈ K .

ˆˆ

ˆx x

The shape functions Ni have been de¬ned in (2.29) (where (x, y) there must

be replaced by (ˆ1 , x2 ) here) for the standard reference element de¬ned

xˆ

there.

2.4. The Implementation of the Finite Element Method: Part 1 79

’1

T

Introducing the matrix C := B ’1 B ’1 = BT B , we can write

˜(m) C ∇x Nj (ˆ) · ∇x Ni (ˆ) |det(B)| dˆ .

Aij = x x x (2.51)

ˆ ˆ

ˆ

K

Denoting the matrix B by B = b(1) , b(2) , then it follows that

’1

b(1) · b(1) b(1) · b(2) b(2) · b(2) ’b(1) · b(2)

1

C= =

b(1) · b(2) b(2) · b(2) ’b(1) · b(2) b(1) · b(1)

det(B)2

because det(B T B) = det(B)2 . The previous considerations can be eas-

ily extended to the computation of the sti¬ness matrices of more general