<<

. 3
( 16)



>>

Then
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

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

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

<<

. 3
( 16)



>>