un ’ um ’0 for n, m ’ ∞ ,

a

un ’ u ’0 for n ’ ∞ ,

a

but u ∈ V , where V is de¬ned analogously to (2.7) with d = 1.

/

In Section 3.1 we will see that a vector space V normed with · a exists

˜

such that u ∈ V and V ‚ V . Therefore, V is not complete with respect

˜ ˜

to · a ; otherwise, u ∈ V must be valid. In fact, there exists a (unique)

completion of V with respect to · a (see Appendix A.4, especially (A4.26)),

but we have to describe the new “functions” added by this process. Besides,

integration by parts must be valid such that a classical solution continues to

be also a weak solution (compare with Lemma 2.1). Therefore, the following

idea is unsuitable.

Attempt of a correct de¬nition of V :

Let V be the set of all u with the property that ‚i u exists for all x ∈ „¦

without any requirements on ‚i u in the sense of a function.

2.1. Variational Formulation 53

For instance, there exists Cantor™s function with the following properties:

f : [0, 1] ’ R, f ∈ C([0, 1]), f = 0, f is not constant, f (x) exists with

f (x) = 0 for all x ∈ [0, 1].

x

Here the fundamental theorem of calculus, f (x) = 0 f (s) ds+f (0), and

thus the principle of integration by parts, are no longer valid.

Consequently, additional conditions for ‚i u are necessary.

To prepare an adequate de¬nition of the space V, we extend the de¬nition

of derivatives by means of their action on averaging procedures. In order

to do this, we introduce the multi-index notation.

A vector ± = (±1 , . . . , ±d ) of nonnegative integers ±i ∈ {0, 1, 2, . . .} is

d

called a multi-index. The number |±| := i=1 ±i denotes the order (or

length) of ±.

For x ∈ Rd let

x± := x±1 · · · x±d . (2.15)

1 d

A shorthand notation for the di¬erential operations can be adopted by this:

For an appropriately di¬erentiable function u let

‚ ± u := ‚1 1 · · · ‚d d u .

± ±

(2.16)

We can obtain this de¬nition from (2.15) by replacing x by the symbolic

vector

T

∇ := (‚1 , . . . , ‚d )

of the ¬rst partial derivatives.

For example, if d = 2 and ± = (1, 2), then |±| = 3 and

‚3u

‚ ± u = ‚1 ‚2 u =

2

.

‚x1 ‚x2

2

Now let ± be a multi-index of length k and let u ∈ C k („¦). We then

∞

obtain for arbitrary test functions • ∈ C0 („¦) by integration by parts

‚ ± u • dx = (’1)k u ‚ ± • dx .

„¦ „¦

The boundary integrals vanish because ‚ β • = 0 on ‚„¦ for all multi-indices

β.

Therefore, we make the following de¬nition:

De¬nition 2.9 v ∈ L2 („¦) is called the weak (or generalized) derivative

∞

‚ ± u of u ∈ L2 („¦) for the multi-index ± if for all • ∈ C0 („¦),

v • dx = (’1)|±| u ‚ ± • dx .

„¦ „¦

54 2. Finite Element Method for Poisson Equation

The weak derivative is well-de¬ned because it is unique: Let v1 , v2 ∈

2

L („¦) be two weak derivatives of u. It follows that

∞

(v1 ’ v2 ) • dx = 0 for all • ∈ C0 („¦) .

„¦

∞

Since C0 („¦) is dense in L2 („¦), we can furthermore conclude that

(v1 ’ v2 ) • dx = 0 for all • ∈ L2 („¦) .

„¦

If we now choose speci¬cally • = v1 ’ v2 , we obtain

v1 ’ v2 (v1 ’ v2 ) (v1 ’ v2 ) dx = 0 ,

2

=

0

„¦

and v1 = v2 (a.e.) follows immediately. In particular, u ∈ C k („¦) has weak

¯

derivatives ‚ ± u for ± with |±| ¤ k, and the weak derivatives are identical

to the classical (pointwise) derivatives.

Also the di¬erential operators of vector calculus can be given a weak

de¬nition analogous to De¬nition 2.9. For example, for a vector ¬eld q

with components in L2 („¦), v ∈ L2 („¦) is the weak divergence v = ∇ · q if

∞

for all • ∈ C0 („¦)

v• dx = ’ q · ∇• dx .

„¦ „¦

1

The correct choice of the space V is the space H0 („¦), which will be

de¬ned below. First we de¬ne

u : „¦ ’ R u ∈ L2 („¦) , u has weak derivatives

H 1 („¦) :=

(2.17)

‚i u ∈ L2 („¦) for all i = 1, . . . , d .

A scalar product on H 1 („¦) is de¬ned by

∇u(x) · ∇v(x) dx

u, v := u(x)v(x) dx + (2.18)

1

„¦ „¦

with the norm

1/2

|u(x)| dx + |∇u(x)| dx

2 2

u := u, u = (2.19)

1 1

„¦ „¦

induced by this scalar product.

The above “temporary” de¬nition (2.7) of V takes care of the boundary

condition u = 0 on ‚„¦ by conditions for the functions. I.e. we want to

choose the basic space V analogously as:

H0 („¦) := u ∈ H 1 („¦) u = 0

1

on ‚„¦ . (2.20)

Here H 1 („¦) and H0 („¦) are special cases of so-called Sobolev spaces.

1

For „¦ ‚ Rd , d ≥ 2, H 1 („¦) may contain unbounded functions. In par-

ticular, we have to examine carefully the meaning of u|‚„¦ (‚„¦ has the

2.2. The Finite Element Method with Linear Elements 55

d-dimensional measure 0) and, in particular, u = 0 on ‚„¦. This will be

described in Section 3.1.

Exercises

2.1

(a) Consider the interval (’1, 1); prove that the function u(x) = |x| has

the generalized derivative u (x) = sign(x).

(b) Does sign(x) have a generalized derivative?

N

2.2 Let „¦ = l=1 „¦l , N ∈ N, where the bounded subdomains „¦l ‚ R2

are pairwise disjoint and possess piecewise smooth boundaries. Show that

a function u ∈ C(„¦) with u|„¦l ∈ C 1 („¦l ), 1 ¤ l ¤ N, has a weak derivative

‚i u ∈ L2 („¦), i = 1, 2, that coincides in N „¦l with the classical one.

l=1

2.3 Let V be the set of functions that are continuous and piecewise con-

tinuously di¬erentiable on [0, 1] and that satisfy the additional conditions

u(0) = u(1) = 0. Show that there exist in¬nitely many elements in V that

minimize the functional

1

2

1 ’ [u (x)]2

F (u) := dx.

0

2.2 The Finite Element Method

with Linear Elements

The weak formulation of the boundary value problem (2.1), (2.2) leads to

particular cases of the following general, here equivalent, problems:

Let V be a vector space, let a : V — V ’ R be a bilinear form, and let

b : V ’ R be a linear form.

Variational equation:

Find u ∈ V a(u, v) = b(v) for all v ∈ V .

with (2.21)

Minimization problem:

Find u ∈ V with F (u) = min F (v) v ∈ V ,

(2.22)

1

where F (v) = a(v, v) ’ b(v) .

2

The discretization approach consists in the following procedure: Replace

V by a ¬nite-dimensional subspace Vh ; i.e., solve instead of (2.21) the ¬nite-

dimensional variational equation,

¬nd uh ∈ Vh for all v ∈ Vh .

with a(uh , v) = b(v) (2.23)

56 2. Finite Element Method for Poisson Equation

This approach is called the Galerkin method. Or solve instead of (2.22) the

¬nite-dimensional minimization problem,

¬nd uh ∈ Vh F (uh ) = min F (v) v ∈ Vh .

with (2.24)

This approach is called the Ritz method.

It is clear from Lemma 2.3 and Remark 2.4 that the Galerkin method

and the Ritz method are equivalent for a positive and symmetric bilinear

form. The ¬nite-dimensional subspace Vh is called an ansatz space.

The ¬nite element method can be interpreted as a Galerkin method (and

in our example as a Ritz method, too) for an ansatz space with special

properties. In the following, these properties will be extracted by means of

the simplest example.

1

Let V be de¬ned by (2.7) or let V = H0 („¦).

The weak formulation of the boundary value problem (2.1), (2.2)

corresponds to the choice

∇u · ∇v dx ,

a(u, v) := b(v) := f v dx .

„¦ „¦

Let „¦ ‚ R2 be a domain with a polygonal boundary; i.e., the boundary

“ of „¦ consists of a ¬nite number of straight-line segments as shown in

Figure 2.2.

„¦

Figure 2.2. Domain with a polygonal boundary.

Let Th be a partition of „¦ into closed triangles K (i.e., including the

boundary ‚K) with the following properties:

(1) „¦ = ∪K∈Th K;

(2) For K, K ∈ Th , K = K ,

int (K) © int (K ) = … , (2.25)

where int (K) denotes the open triangle (without the boundary ‚K).

(3) If K = K but K ©K = …, then K ©K is either a point or a common

edge of K and K (cf. Figure 2.3).

A partition of „¦ with the properties (1), (2) is called a triangulation

of „¦. If, in addition, a partition of „¦ satis¬es property (3), it is called a

conforming triangulation (cf. Figure 2.4).

2.2. Linear Elements 57

r ¨ r

¨

rr rr d

¨

rr ¨¨ rr d

not

r

¨

¨rr d

¨

allowed: ¨

allowed:

¨ ¨

rr d

¨¨ ¨

¨

r d

¨ r ¨ d

Figure 2.3. Triangulations.

The triangles of a triangulation will be numbered K1 , . . . , KN . The

subscript h indicates the ¬neness of the triangulation, e.g.,

h := max diam (K) K ∈ Th ,

where diam (K) := sup |x ’ y| x, y ∈ K denotes the diameter of K.

Thus here h is the maximum length of the edges of all the triangles.

Sometimes, K ∈ Th is also called a (geometric) element of the partition.

The vertices of the triangles are called the nodes, and they will be

numbered

a1 , a2 , . . . , aM ,

i.e., ai = (xi , yi ), i = 1, . . . , M , where M = M1 + M2 and

∈

a1 , . . . , aM1 „¦,

(2.26)

∈

aM1 +1 , . . . , aM ‚„¦ .

This kind of arrangement of the nodes is chosen only for the sake

of simplicity of the notation and is not essential for the following

considerations.

a• a8

•

a9 d

10

£d

• £

¨

¨t d

K1 £

K2 ¨

• K5

¨ £ K6 d

t

a11 t a1 ¨

•

¨ d

t£ $•

K3 $

K4

t ¡t a2 t £ $$$ 5 a7 5

K12 t$

t¡t • • 5

¢ a3 5

K10 f

t• K11t ¢

t¡ f K8¢ K7 5

a4

t• 5

K9 f

t ¢5

a5 f•

5

f¢

a 6

Figure 2.4. A conforming triangulation with N = 12, M = 11, M1 = 3, M2 = 8.

An approximation of the boundary value problem (2.1), (2.2) with linear

¬nite elements on a given triangulation Th of „¦ is obtained if the ansatz

space Vh is de¬ned as follows:

Vh := u ∈ C(„¦) u|K ∈ P1 (K) for all K ∈ Th , u = 0 on ‚„¦ .

¯ (2.27)

58 2. Finite Element Method for Poisson Equation

Here P1 (K) denotes the set of polynomials of ¬rst degree (in 2 variables)

on K; i.e., p ∈ P1 (K) ” p(x, y) = ± + βx + γy for all (x, y) ∈ K and for

¬xed ±, β, γ ∈ R.

Since p ∈ P1 (K) is also de¬ned on the space R — R, we use the short but

inaccurate notation P1 = P1 (K); according to the context, the domain of

de¬nition will be given as R — R or as a subset of it.

We have

Vh ‚ V .

This is clear for the case of de¬nition of V by (2.7) because ‚x u|K = const,

‚y u|K = const for K ∈ Th for all u ∈ Vh . If V = H0 („¦), then this inclusion

1

is not so obvious. A proof will be given in Theorem 3.20 below.

An element u ∈ Vh is determined uniquely by the values u(ai ), i =

1, . . . , M1 (the nodal values).

In particular, the given nodal values already enforce the continuity of the

piecewise linear composed functions. Correspondingly, the homogeneous

Dirichlet boundary condition is satis¬ed if the nodal values at the boundary

nodes are set to zero.

In the following, we will demonstrate these properties by an unnecessar-

ily involved proof. The reason is that this proof will introduce all of the

considerations that will lead to analogous statements for the more general

problems of Section 3.4.

Let Xh be the larger ansatz space consisting of continuous, piecewise

linear functions but regardless of any boundary conditions, i.e.,

Xh := u ∈ C(„¦) u|K ∈ P1 (K) for all K ∈ Th .

¯

Lemma 2.10 For given values at the nodes a1 , . . . , aM , the interpolation

problem in Xh is uniquely solvable. That is, if the values u1 , . . . , uM are

given, then there exists a uniquely determined element

u ∈ Xh such that u(ai ) = ui , i = 1, . . . , M .

If uj = 0 for j = M1 + 1, . . . , M , then it is even true that

u ∈ Vh .

Proof: (1) For any arbitrary K ∈ Th we consider the local interpolation

problem:

Find p = pK ∈ P1 such that p(ai ) = ui , i = 1, 2, 3 , (2.28)

where ai , i = 1, 2, 3, denote the vertices of K, and the values ui , i = 1, 2, 3,

are given. First we show that problem (2.28) is uniquely solvable for a

particular triangle.

ˆ

A solution of (2.28) for the so-called reference element K (cf. Figure 2.5)

with the vertices a1 = (0, 0), a2 = (1, 0), a3 = (0, 1) is given by

ˆ ˆ ˆ

p(x, y) = u1 N1 (x, y) + u2 N2 (x, y) + u3 N3 (x, y)

2.2. Linear Elements 59

y

1

^

K

0 1 x

ˆ

Figure 2.5. Reference element K.

with the shape functions

1’x’y,

N1 (x, y) =

N2 (x, y) = x, (2.29)

N3 (x, y) = y.

Evidently, Ni ∈ P1 , and furthermore,

1 for i = j ,

Ni (ˆj ) = δij =

a for i, j = 1, 2, 3 ,

0 for i = j ,

and thus

3

p (ˆj ) =

a ui Ni (ˆj ) = uj

a for all j = 1, 2, 3 .

i=1

The uniqueness of the solution can be seen in the following way: If p1 , p2

satisfy the interpolation problem (2.28) for the reference element, then for

p := p1 ’ p2 ∈ P1 we have

p (ˆi ) = 0 ,

a i = 1, 2, 3 .

Here p is given in the form p(x, y) = ± + βx + γy. If we ¬x the second

variable y = 0, we obtain a polynomial function of one variable

p(x, 0) = ± + βx =: q(x) ∈ P1 (R) .

The polynomial q satis¬es q(0) = 0 = q(1), and q ≡ 0 follows by the

uniqueness of the polynomial interpolation in one variable; i.e., ± = β = 0.

Analogously, we consider

q(y) := p(0, y) = ± + γy = γy ,

and we obtain from q(1) = 0 that γ = 0 and consequently p ≡ 0.

In fact, this additional proof of uniqueness is not necessary, because the

uniqueness already follows from the solvability of the interpolation problem

because of dim P1 = 3 (compare with Section 3.3).

Now we turn to the case of a general triangle K. A general triangle K is

ˆ

mapped onto K by an a¬ne transformation (cf. Figure 2.6)

F :K ’K,

ˆ F (ˆ) = B x + d ,

x ˆ (2.30)

where B ∈ R2,2 , d ∈ R2 are such that F (ˆi ) = ai .

a

60 2. Finite Element Method for Poisson Equation

B = (b1 , b2 ) and d are determined by the vertices ai of K as follows:

a1 = F (ˆ1 ) = F (0) = d ,

a

a2 = F (ˆ2 ) = b1 + d = b1 + a1 ,

a

a3 = F (ˆ3 ) = b2 + d = b2 + a1 ;

a

i.e., b1 = a2 ’ a1 and b2 = a3 ’ a1 . The matrix B is regular because a2 ’ a1

and a3 ’ a1 are linearly independent, ensuring F (ˆi ) = ai .

a

Since

3 3

K = conv {a1 , a2 , a3 } := »i ai 0 ¤ »i ¤ 1 , »i = 1

i=1 i=1

and especially K = conv {ˆ1 , a2 , a3 }, F [K] = K follows from the fact that

ˆ ˆ

aˆˆ

the a¬ne-linear mapping F satis¬es

3 3 3

F »i ai

ˆ = »i F (ˆi ) =

a »i ai

i=1 i=1 i=1

3

for 0 ¤ »i ¤ 1, i=1 »i = 1.

ˆ

In particular, the edges (where one »i is equal to 0) of K are mapped

onto the edges of K.

a3

y

^

y

1 K

a2

a1

^

K

x

0 1 ^

x

Figure 2.6. A¬ne-linear transformation.

Analogously, the considerations can be applied to the space Rd word for

word by replacing the set of indices {1, 2, 3} by {1, . . . , d + 1}. This will be

done in Section 3.3.

The polynomial space P1 does not change under the a¬ne transforma-

tion F .

(2) We now prove that the local functions u|K can be composed

continuously:

For every K ∈ Th , let pK ∈ P1 be the unique solution of (2.28), where

the values u1 , u2 , u3 are the values ui1 , ui2 , ui3 (i1 , i2 , i3 ∈ {1, . . . , M }) that

have to be interpolated at these nodes.

Let K, K ∈ Th be two di¬erent elements that have a common edge E.

Then pK = pK on E is to be shown. This is valid because E can be

mapped onto [0, 1] — {0} by an a¬ne transformation (cf. Figure 2.7). Then

2.2. Linear Elements 61

q1 (x) = pK (x, 0) and q2 (x) := pK (x, 0) are elements of P1 (R), and they

solve the same interpolation problem at the points x = 0 and x = 1; thus

q1 ≡ q2 .

0 1

E

Figure 2.7. A¬ne-linear transformation of E on the reference element [0, 1].

Therefore, the de¬nition of u by means of

for x ∈ K ∈ Th

u(x) = pK (x) (2.31)

is unique, and this function satis¬es u ∈ C(„¦) and u ∈ Xh .

¯

(3) Finally, we will show that u = 0 on ‚„¦ for u de¬ned by (2.31) if

ui = 0 (i = M1 + 1, . . . , M ) for the boundary nodes.

The boundary ‚„¦ consists of edges of elements K ∈ Th . Let E be such

an edge; i.e., E has the vertices ai1 , ai2 with ij ∈ {M1 + 1, . . . , M }. The

given boundary values yield u(aij ) = 0 for j = 1, 2. By means of an a¬ne

transformation analogously to the above one we obtain that u|E is a poly-

nomial of ¬rst degree in one variable and that u|E vanishes at two points.

2

So u|E = 0, and the assertion follows.

The following statement is an important consequence of the unique solv-

ability of the interpolation problem in Xh irrespective of its particular

de¬nition: The interpolation conditions

•i (aj ) = δij , j = 1, . . . , M , (2.32)

uniquely determine functions •i ∈ Xh for i = 1, . . . , M. For any u ∈ Xh ,

we have

M

for x ∈ „¦ ,

u(x) = u(ai )•i (x) (2.33)

i=1

because both the left-hand side and the right-hand side functions belong

to Xh and are equal to u(ai ) at x = ai .

M

The representation u = i=1 ±i •i is unique, too, for otherwise, a func-

tion w ∈ Xh , w = 0, such that w(ai ) = 0 for all i = 1, . . . , M would

exist. Thus {•1 , . . . , •M } is a basis of Xh , especially dim Xh = M . This

basis is called a nodal basis because of (2.33). For the particular case of a

piecewise linear ansatz space on triangles, the basis functions are called

62 2. Finite Element Method for Poisson Equation

pyramidal functions because of their shape. If the set of indices is re-

stricted to {1, . . . , M1 }; i.e., we omit the basis functions corresponding to

the boundary nodes, then a basis of Vh will be obtained and dim Vh = M1 .

Summary: The function values u(ai ) at the nodes a1 , . . . , aM are the de-

grees of freedom of u ∈ Xh , and the values at the interior points a1 , . . . , aM1

are the degrees of freedom of u ∈ Vh .

The following consideration is valid for an arbitrary ansatz space Vh with

a basis {•1 , . . . , •M }. The Galerkin method (2.23) reads as follows: Find

M

uh = i=1 ξi •i ∈ Vh such that a(uh , v) = b(v) for all v ∈ Vh . Since

M

v = i=1 ·i •i for ·i ∈ R, this is equivalent to

b(•i ) for all i = 1, . . . , M ⇐’

a(uh , •i ) =

«

M

a ξj •j , •i b(•i ) for all i = 1, . . . , M ⇐’

=

j=1

M

b(•i ) for all i = 1, . . . , M ⇐’

a (•j , •i ) ξj =

j=1

Ah ξ = qh (2.34)

T

with Ah = (a(•j , •i ))ij ∈ RM,M , ξ = (ξ1 , . . . , ξM ) and q h = (b(•i ))i .

Therefore, the Galerkin method is equivalent to the system of equations

(2.34).

The considerations for deriving (2.34) show that, in the case of equiva-

lence of the Galerkin method with the Ritz method, the system of equations

(2.34) is equivalent to the minimization problem

Fh (ξ) = min Fh (·) · ∈ RM , (2.35)

where

1T

· Ah · ’ q T · .

Fh (·) = h

2

Because of the symmetry and positive de¬niteness, the equivalence of (2.34)

and (2.35) can be easily proven, and it forms the basis for the CG methods

that will be discussed in Section 5.2.

Usually, Ah is called sti¬ness matrix, and q h is called the load vector.

These names originated from mechanics. For our model problem, we have

∇•j · ∇•i dx ,

(Ah )ij = a(•j , •i ) =

„¦

(qh )i = b(•i ) = f •i dx .

„¦

By applying the ¬nite element method, we thus have to perform the

following steps:

(1) Determination of Ah , q h . This step is called assembling.

2.2. Linear Elements 63

(2) Solution of Ah ξ = q h .

If the basis functions •i have the property •i (aj ) = δij , then the solution

of system (2.34) satis¬es the relation ξi = uh (ai ), i.e., we obtain the vector

of the nodal values of the ¬nite element approximation.

Using only the properties of the bilinear form a, we obtain the following

properties of Ah :

• Ah is symmetric for an arbitrary basis {•i } because a is symmetric.

• Ah is positive de¬nite for an arbitrary basis {•i } because for u =

M

i=1 ξi •i ,

M M M

ξT Ah ξ = i,j=1 ξj a(•j , •i )ξi = j=1 ξj a •j , i=1 ξi •i

M M

=a j=1 ξj •j , i=1 ξi •i = a(u, u) > 0

(2.36)

for ξ = 0 and therefore u ≡ 0.

Here we have used only the positive de¬niteness of a.

Thus we have proven the following lemma.

Lemma 2.11 The Galerkin method (2.23) has a unique solution if a is a

symmetric, positive de¬nite bilinear form and if b is a linear form.

In fact, as we will see in Theorem 3.1, the symmetry of a is not necessary.

• For a special basis (i.e., for a speci¬c ¬nite element method), Ah is a

sparse matrix, i.e., only a few entries (Ah )ij do not vanish. Evidently,

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

(Ah )ij = 0

„¦

This can happen only if supp •i © supp •j = …, as this property is

again necessary for supp ∇•i © supp ∇•j = … because of

(supp ∇•i © supp ∇•j ) ‚ (supp •i © supp •j ) .

The basis function •i vanishes on an element that does not contain

the node ai because of the uniqueness of the solution of the local

interpolation problem. Therefore,

supp •i = K,

K∈Th

ai ∈K

cf. Figure (2.8), and thus

’ ai , aj ∈ K for some K ∈ Th ;

(Ah )ij = 0 (2.37)

i.e., ai , aj are neighbouring nodes.

64 2. Finite Element Method for Poisson Equation

If we use the piecewise linear ansatz space on triangles and if ai is an interior

node in which L elements meet, then there exist at most L nondiagonal

entries in the ith row of Ah . This number is determined only by the type

of the triangulation, and it is independent of the ¬neness h, i.e., of the

number of unknowns of the system of equations.

supp • i

ai

Figure 2.8. Support of the nodal basis function.

Example 2.12 We consider again the boundary value problem (2.1), (2.2)

on „¦ = (0, a) — (0, b) again, i.e.

’∆u =f in „¦ ,

u =0 on ‚„¦ ,

under the condition (1.4). The triangulation on which the method is based

is created by a partition of „¦ into squares with edges of length h and by

a subsequent uniform division of each square into two triangles according

to a ¬xed rule (Friedrichs“Keller triangulation). In order to do this, two

possibilities (a) and (b) (see Figures 2.9 and 2.10) exist.

(a) (b)

d d

d

d d d

d

d d

d

d d

d d d

d

d d

d

d d

d d d

d

d d

d

d d

d d d

d d d

Figure 2.9. Possibilities of Friedrichs“Keller triangulation.

In both cases, a node aZ belongs to six elements, and consequently, it

has at most six neighbours:

2.2. Linear Elements 65

aNW aN

d II d

for (a): for (b):

d Id

III aZ d

d

aW aE

d IV dVI

d Vd

d d

aS aSE

Figure 2.10. Support of the basis function.

Case (a) becomes case (b) by the transformation x ’ a ’ x, y ’ y.

This transformation leaves the di¬erential equation or the weak formula-

tion, respectively, unchanged. Thus the Galerkin method with the ansatz

space Vh according to (2.27) does not change, because P1 is invariant with

respect to the above transformation. Therefore, the discretization matrices

Ah according to (2.34) are seen to be identical by taking into account the

renumbering of the nodes by the transformation.

Thus it is su¬cient to consider only one case, say (b). A node which is

far away from the boundary has 6 neighbouring nodes in {a1 , . . . , aM1 },

a node close to the boundary has less. The entries of the matrix in the

row corresponding to aZ depend on the derivatives of the basis function

•Z as well as on the derivatives of the basis functions corresponding to the

neighbouring nodes. The values of the partial derivatives of •Z in elements

having the common vertex aZ are listed in Table 2.1, where these elements

are numbered according to Figure 2.10.

I II III IV V VI

’h ’h

1 1 1 1

‚1 •Z 0 0

h h

’h ’h

1 1 1 1

‚2 •Z 0 0

h h

Table 2.1. Derivatives of the basis functions.

Thus for the entries of the matrix in the row corresponding to aZ we

have

2 2 2

|∇•Z | dx = 2

(Ah )Z,Z=a(•Z , •Z ) = (‚1 •Z ) + (‚2 •Z ) dx,

I∪...∪VI I∪II∪III

because the integrands are equal on I and IV, on II and V, and on III and

VI. Therefore

(‚2 •Z )2 dx = 2h’2 h2 + 2h’2 h2 = 4 ,

(‚1 •Z )2 dx + 2

(Ah )Z,Z = 2

I∪III I∪II

∇•N · ∇•Z dx

(Ah )Z,N = a (•N , •Z ) =

I∪II

’h’1 h’1 dx = ’1 ,

= ‚2 •N ‚2 •Z dx =

I∪II I∪II

66 2. Finite Element Method for Poisson Equation

because ‚1 •Z = 0 on II and ‚1 •N = 0 on I. The element I for •N corre-

sponds to the element V for •Z ; i.e., ‚1 •N = 0 on I, analogously, it follows

that ‚2 •N = h’1 on I ∪ II. In the same way we get

(Ah )Z,E = (Ah )Z,W = (Ah )Z,S = ’1

as well as

(Ah )Z,NW = a (•NW , •Z ) = ‚1 •NW ‚1 •Z + ‚2 •NW ‚2 •Z dx = 0 .

II∪III

The last identity is due to ‚1 •NW = 0 on III and ‚2 •NW = 0 on III,

because the elements V and VI for •Z agree with the elements III and II

for •NW , respectively.

Analogously, we obtain for the remaining value

(Ah )Z,SE = 0 ,

such that only 5 (instead of the maximum 7) nonzero entries per row exist.

The way of assembling the sti¬ness matrix described above is called node-

based assembling. However, most of the computer programs implementing

the ¬nite element method use an element-based assembling, which will be

considered in Section 2.4.

If the nodes are numbered rowwise analogously to (1.13) and if the equa-

tions are divided by h2 , then h’2 Ah coincides with the discretization matrix

(1.14), which is known from the ¬nite di¬erence method. But here the

right-hand side is given by

h’2 (qh )i = h’2 f •i dx = h’2 f •i dx

„¦ I∪...∪VI

for aZ = ai and thus it is not identical to f (ai ), the right-hand side of the

¬nite di¬erence method.

However, if the trapezoidal rule, which is exact for g ∈ P1 , is applied to

approximate the right-hand side according to

3

1

g(x) dx ≈ vol (K) g(ai ) (2.38)

3

K i=1

for a triangle K with the vertices ai , i = 1, 2, 3 and with the area vol (K),

then

11 2 1

f •i dx ≈ h (f (aZ ) · 1 + f (aO ) · 0 + f (aN ) · 0) = h2 f (aZ ).

32 6

I

Analogous results are obtained for the other triangles, and thus

h’2 f •i dx ≈ f (aZ ) .

I∪...∪VI

In summary, we have the following result.

2.2. Linear Elements 67

Lemma 2.13 The ¬nite element method with linear ¬nite elements on a

triangulation according to Figure 2.9 and with the trapezoidal rule to ap-

proximate the right-hand side yields the same discretization as the ¬nite

di¬erence method from (1.7), (1.8).

We now return to the general formulation (2.21)“(2.24). The approach

of the Ritz method (2.24), instead of the Galerkin method (2.23), yields an

identical approximation because of the following result.

Lemma 2.14 If a is a symmetric and positive bilinear form and b is a

linear form, then the Galerkin method (2.23) and the Ritz method (2.24)

have identical solutions.

2

Proof: Apply Lemma 2.3 with Vh instead of V .

Hence the ¬nite element method is the Galerkin method (and in our

problem the Ritz method, too) for an ansatz space Vh with the following

properties:

• The coe¬cients have a local interpretation (here as nodal values).

The basis functions have a small support such that:

• the discretization matrix is sparse,

• the entries of the matrix can be assembled locally.

Finally, for the boundary value problem (2.1), (2.2) with the correspond-

ing weak formulation, we consider other ansatz spaces, which to some extent

do not have these properties:

(1) In Section 3.2.1, (3.28), we will show that mixed boundary conditions

need not be included in the ansatz space. Then we can choose the ¬-

nite dimensional polynomial space Vh = span 1, x, y, xy, x2 , y 2 , . . .

for it. But in this case, Ah is a dense matrix and ill-conditioned. Such

ansatz spaces yield the classical Ritz“Galerkin methods.

(2) Let Vh = span {•1 , . . . , •N } and let •i ≡ 0 satisfy, for some »i ,

for all v ∈ V ,

a(•i , v) = »i •i , v 0

i.e., the weak formulation of the eigenvalue problem

’∆u = »u in „¦ ,

u= 0 on ‚„¦ ,

for which eigenvalues 0 < »1 ¤ »2 ¤ . . . and corresponding eigen-

functions •i exist such that •i , •j 0 = δij (e.g., see [12, p. 335]). For

special domains „¦, (»i , •i ) can be determined explicitly, and

(Ah )ij = a(•j , •i ) = »j •j , •i = »j δij

0

68 2. Finite Element Method for Poisson Equation

is obtained. Thus Ah is a diagonal matrix, and the system of equations

Ah ξ = q h can be solved without too great expense. But this kind of

assembling is possible with acceptable costs for special cases only.

(3) The (spectral) collocation method consists in the requirement that

the equations (2.1), (2.2) be satis¬ed only at certain distinct points

xi ∈ „¦, called collocation points, for a special polynomial space Vh .

The above examples describe Galerkin methods without having the typical

properties of a ¬nite element method.

2.3 Stability and Convergence of the Finite

Element Method

We consider the general case of a variational equation of the form (2.21)

and the Galerkin method (2.23). Here let a be a bilinear form, which is not

necessarily symmetric, and let b be a linear form.

Then, if

e := u ’ uh (∈ V )

denotes the error, the important error equation

for all v ∈ Vh

a(e, v) = 0 (2.39)

is satis¬ed. To obtain this equation, it is su¬cient to consider equation

(2.21) only for v ∈ Vh ‚ V and then to subtract from the result the

Galerkin equation (2.23).

If, in addition, a is symmetric and positive de¬nite, i.e.,

a(u, u) ≥ 0 , a(u, u) = 0 ” u = 0

a(u, v) = a(v, u) ,

(i.e., a is a scalar product), then the error is orthogonal to the space Vh

with respect to the scalar product a.

Therefore, the relation (2.39) is often called the orthogonality of the error

(to the ansatz space). In general, the element uh ∈ Vh with minimal distance

to u ∈ V with respect to the induced norm · a is characterized by (2.39):

Lemma 2.15 Let Vh ‚ V be a subspace, let a be a scalar product on V ,

and let u a := a(u, u)1/2 be the norm induced by a. Then for uh ∈ Vh , it

follows that

a(u ’ uh , v) for all v ∈ Vh ”

= 0 (2.40)

u ’ uh u’v v ∈ Vh .

= min (2.41)

a a

Proof: For arbitrary but ¬xed u ∈ V , let b(v) := a(u, v) for v ∈ Vh .

Then b is a linear form on Vh , so (2.40) is a variational formulation on Vh .

2.3. Stability and Convergence 69

According to Lemma 2.14 or Lemma 2.3, this variational formulation has

the same solutions as

min F (v) v ∈ Vh

F (uh ) =

1 1

a(v, v) ’ b(v) = a(v, v) ’ a(u, v) .

with F (v) :=

2 2

Furthermore, F has the same minima as the functional

1/2 1/2

a(v, v) ’ 2a(u, v) + a(u, u)

2F (v) + a(u, u) =

1/2

a(u ’ v, u ’ v) = u’v

= ,

a

because the additional term a(u, u) is a constant. Therefore, F has the

2

same minima as (2.41).

If an approximation uh of u is to be sought exclusively in Vh , then the

element uh , determined by the Galerkin method, is the optimal choice with

respect to · a .

A general, not necessarily symmetric, bilinear form a is assumed to satisfy

the following conditions, where · denotes a norm on V :

• a is continuous with respect to · ; i.e., there exists M > 0 such that

|a(u, v)| ¤ M u for all u, v ∈ V ;

v (2.42)

• a is V -elliptic; i.e., there exists ± > 0 such that

a(u, u) ≥ ± u for u ∈ V .

2

(2.43)

If a is a scalar product, then (2.42) with M = 1 and (2.43) (as equality)

with ± = 1 are valid for the induced norm · := · a due to the

Cauchy“Schwarz inequality.

The V -ellipticity is an essential condition for the unique existence of a

solution of the variational equation (2.21) and of the boundary value prob-

lem described by it, which will be presented in more detail in Sections 3.1

and 3.2. It also implies ” without further conditions ” the stability of the

Galerkin approximation.

Lemma 2.16 The Galerkin solution uh according to (2.23) is stable in

the following sense:

1

uh ¤ b independently of h , (2.44)

±

where

|b(v)|

v∈V ,v=0

b := sup .

v

70 2. Finite Element Method for Poisson Equation

Proof: In the case uh = 0, there is nothing to prove. Otherwise, from

a(uh , v) = b(v) for all v ∈ Vh , it follows that

|b(uh )|

¤ a(uh , uh ) = b(uh ) ¤ uh ¤ b

2

± uh uh .

uh

2

Dividing this relation by ± uh , we get the assertion.

Moreover, the approximation property (2.41) holds up to a constant:

Theorem 2.17 (C´a™s lemma)

e

Assume (2.42), (2.43). Then the following error estimate for the Galerkin

solution holds:

M

u ’ uh ¤ min u ’ v v ∈ Vh . (2.45)

±

Proof: If u ’ uh = 0, then there is nothing to prove. Otherwise, let

v ∈ Vh be arbitrary. Because of the error equation (2.39) and uh ’ v ∈ Vh ,

a(u ’ uh , uh ’ v) = 0 .

Therefore, using (2.43) we have

± u ’ uh ¤ a(u ’ uh , u ’ uh ) = a(u ’ uh , u ’ uh ) + a(u ’ uh , uh ’ v)

2

= a(u ’ uh , u ’ v) .

Furthermore, by means of (2.42) we obtain

± u ’ uh ¤ a(u ’ uh , u ’ v) ¤ M u ’ uh u ’ v for arbitrary v ∈ Vh .

2

Thus the assertion follows by division by ± u ’ uh . 2

Therefore also in general, in order to get an asymptotic error estimate

in h, it is su¬cient to estimate the best approximation error of Vh , i.e.,

u’v v ∈ Vh .

min

However, this consideration is meaningful only in those cases where M/±

is not too large. Section 3.2 shows that this condition is no longer satis¬ed

for convection-dominated problems. Therefore, the Galerkin approach has

to be modi¬ed, which will be described in Chapter 9.

We want to apply the theory developed up to now to the weak formula-

tion of the boundary value problem (2.1), (2.2) with V according to (2.7)

or (2.20) and Vh according to (2.27). According to (2.4) the bilinear form

a and the linear form b read as

∇u · ∇v dx ,

a(u, v) = b(v) = f v dx .

„¦ „¦

In order to guarantee that the linear form b is well-de¬ned on V, it is su¬-

cient to assume that the right-hand side f of the boundary value problem

belongs to L2 („¦).

2.3. Stability and Convergence 71

Since a is a scalar product on V ,

1/2

|∇u| dx 2

u=u =

a

„¦

is an appropriate norm. Alternatively, the norm introduced in (2.19) for

1

V = H0 („¦) can be taken as

1/2

|u(x)| dx + |∇u(x)| dx

2 2

u = .

1

„¦ „¦

In the latter case, the question arises whether the conditions (2.42) and

(2.43) are still satis¬ed. Indeed,

|a(u, v)| ¤ u ¤u for all u, v ∈ V .

v v

a a 1 1

The ¬rst inequality follows from the Cauchy“Schwarz inequality for the

scalar product a, and the second inequality follows from the trivial estimate

1/2

|∇u(x)| dx ¤u for all u ∈ V .

2

u =

a 1

„¦

Thus a is continuous with respect to · 1 with M = 1.

The V -ellipticity of a, i.e., the property

≥± u for some ± > 0 and all u ∈ V ,

2 2

a(u, u) = u a 1

is not valid in general for V = H 1 („¦). However, in the present situation

1

of V = H0 („¦) it is valid because of the incorporation of the boundary

condition into the de¬nition of V :

Theorem 2.18 (Poincar´) Let „¦ ‚ Rn be open and bounded. Then a

e

constant C > 0 exists (depending on „¦) such that

1/2

¤C |∇u(x)| dx for all u ∈ H0 („¦) .

2 1

u 0

„¦

2

Proof: Cf. [13]. For a special case, see Exercise 2.5.

Thus (2.43) is satis¬ed, for instance with

1

±= ,

1 + C2

(see also (3.26) below) and thus in particular

¤ a(u, u) = u ¤u for all u ∈ V ,

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