The Finite Element Method: A Short Introduction
What is FEM?
The Finite Element Method (FEM) introduced by engineers in late 50’s and 60’s is a numerical
techinique for solving problems which are described by Ordinary Differential Equations (ODE) /Partial
Differential Equations (PDE) with appropriate boundary/initial conditions or to solve problems that
can be formulated as a functional minimization.
Why FEM?
• Greater flexibility to model complex geometries.
• Can handle general boundary conditions.
• Variable material properties can be handled.
• Clear structure and versatility helps to construct general purpose software for applications.
• Has a solid theoretical foundation which gives added reliability and makes it possible to mathematically
analyze and estimate the error in the approximate solution.
We now describe the important steps of classical/standard finite element method. We start with
a boundary value problem in the classical formulation. After defining an appropriate variational formulation for the bounday value problem, we give a finite element formulation.
Step 1: We start with the following two point boundary value problem.

Given 0 f 0 find 0 u0 such that

0 0
(PC ) : 
 Au ≡ −(a(x)u ) + a0 (x)u = f in I = (0, 1),
u(0) = u(1) = 0.
(1)
The following assumptions on the coefficients are made:
¯ f ∈ L2 (I),
a(x) and a0 (x) are smooth functions with a(x) ≥ α0 > 0, a0 (x) ≥ 0 in I,
maxx∈I¯(|a(x)| , |a0 (x)|) ≤ M.
Step 2: Weak Formulation
For a boundary value problem defined using an operator of order 2m, admissible space is chosen as
H m (Ω) along with all essential boundary conditions i.e., boundary conditions involving derivatives of
u of order ≤ m − 1. In this case, we choose the admissible space as V = H01 (I).
Multiplying (1) by v ∈ V , integrating left hand side of (1) by parts and using the boundary conditions,
we get:

Given f in L2 (I), find u in V such that

a(u, v) = l(v)
∀v ∈ V,
(2)
where
• a(·, ·) : V × V → R is a symmetric, continuous, bilinear form defined by:
Z
a(u, v) = (a(x)u0 v 0 + a0 (x)uv) dx
∀u, v ∈ V
I
1
(3)
• l(·) : V → R is a continuous linear form defined by:
Z
l(v) = f v dx
∀v ∈ V.
(4)
I
Lemma: The weak formulation is well-posed.
Step (2)0 : Energy Minimization Formulation


Given f in L2 (I), find u in V such that




 J(u) = inf J(v)
v∈V
(5)


∀v ∈ V,
where J(v) = 12 a(v, v) − l(v)




a(·, ·) and l(·) are defined by (3) and (4) respectively.
Step 3: Galerkin Finite Element Problem
Firstly, we construct a finite dimensional subspace Vh of V .
For any positive integer N + 1, let {0 = x0 < x1 < ... < xN +1 = 1} be a partition of I¯ into subintervals
(finite elements) Ij = (xj−1 , xj ), 1 ≤ j ≤ N + 1, with length hj = xj − xj−1 . h =
max
1≤j≤N +1
hj is the
parameter associated with the mesh.
The discrete solution will be sought in Vh defined by:
¯ vh I ∈ P1 (Ij ), 1 ≤ j ≤ N + 1, vh (0) = vh (1) = 0}
Vh = {vh : vh ∈ C 0 (I),
j
(6)
• Vh is a subspace of V .
• Vh is finite dimensional, since on each subinterval Ij , vh ∈ P1 (Ij ); i.e., we have 2(N + 1) degrees of
freedom. Continuity constraints at the nodal points x1 , ..., xN and boundary conditions vh (0) = vh (1) =
0 imply dim Vh = 2(N + 1) − (N ) − 2 = N .
Remark: This is the simplest choice for Vh . Other choices are possible.
The Galerkin finite element problem (PGh ) corresponding to the weak formulation can be defined
as:

(PGh ) : 
Given f in L2 (I), find uh in Vh ⊂ V such that
a(uh , vh ) = l(vh )
∀uh ∈ Vh .
Z
Z duh dvh
+ a0 (x)uh vh dx = f vh dx
∀vh ∈ Vh .
i.e.,
a(x)
dx dx
I
I
h
The Ritz finite element problem (PR ) can be defined as:

Given f in L2 (I), find uh in Vh ⊂ V such that
(PRh ) : 
J(uh ) = inf J(vh ),
vh ∈Vh
where J(vh ) = 21 a(vh , vh ) − l(vh ).
Using Lax-Milgram lemma, (PGh ) is well-posed.
Step 4: Construction of Canonical Basis Functions for Vh .
2
(7)
(8)
dim Vh = N . The basis functions for Vh can be chosen to be the Lagrange polynomials {φih }N
i=1
which satisfy φih (xj ) = δij 1 ≤ i, j ≤ N , where

 1
δij =
 0
when i = j
when i 6= j

x − xi−1




 xi − xi−1
x − xi+1
i.e., φih (x) =

x
 i − xi+1


 0,
in [xi−1 , xi ]
in [xi , xi+1 ]
otherwise.
Note that supp φih = [xi−1 , xi+1 ].
N
N
X
X
For vh ∈ Vh , vh (x) =
βi φih (x), vh (xj ) =
βi φjh (xj ) = βj . Hence βj are nothing but values of
i=1
i=1
vh at the nodal points.
i.e., vh (x) =
N
X
vh (xi )φih (x)
(9)
i=1
The unknown function uh ∈ Vh can also be expressed as
uh =
N
X
uh (xi )φih (x) =
i=1
N
X
ui φih (x)
(10)
i=1
with (ui )1≤i≤N = (uh (xi ))1≤i≤N .
Step 5: Reduction of (PGh ) into Matrix Equations.
From (7), a(uh , vh ) = l(vh )
∀vh ∈ Vh and from (10), uh =
!
N
X
j
i
i.e., a
ui φh , φh = l(φjh )
1 ≤ j ≤ N.
N
X
ui φih
i=1
i=1
⇒
N
X
ui a(φih , φjh ) = l(φjh )
1 ≤ j ≤ N.
i=1
⇒ KU = F where K = [a(φih , φjh )]1≤i,j≤N , U = (ui )1≤i≤N , F = (l(φjh ))1≤j≤N .
Properties of K:
(1) K is sparse: a(φih , φjh ) = 0 whenever Supp φih and Supp φjh does not intersect. Since Supp φih =
[xi−1 , xi+1 ] in this case, we observe that K is tridiagonal.
(2) K is symmetric: This is because the bilinear form a(·, ·) is symmetric.
(3) K is positive-definite:
UT KU
=
N X
N
X
ui uj a(φih , φjh )
i=1 j=1


N
N
X
X
j
i
= a
ui φh ,
uj φh 
i=1
j=1
= a(uh , uh ) ≥ αkuh k2V > 0.
3
⇒ K is positive-definite.
Step 6: Construction of K, F for implementation purpose
On [xi−1 , xi+1 ],we have

x − xi−1
1




in
[x
,
x
]
in [xi−1 , xi ]
i−1
i




 xi − xi−1
 hi
i
dφ
x − xi+1
1
h
φih =
=
in [xi , xi+1 ]
−
in [xi , xi+1 ]


dx
x
−
x
h


i
i+1
i




 0
 0
otherwise.
otherwise.


x
−
x
1
i−2




in [xi−2 , xi−1 ]
in [xi−2 , xi−1 ]




 hi−1
 xi−1 − xi−2
i−1
dφ
x − xi
1
h
φhi−1 =
=
in [xi−1 , xi ]
−
in [xi−1 , xi ]


dx
x
−
x


h
i−1
i
i




 0
 0
otherwise.
elsewhere.


x − xi
1




in [xi , xi+1 ]
in [xi , xi+1 ]


 xi+1 − xi

hi+1


i+1
dφh
x − xi+2
1
φi+1
=
=
in [xi+1 , xi+2 ]
−
in [xi+1 , xi+2 ]
h


dx
x
−
x
hi+2


i+1
i+2




 0
 0
otherwise.
elsewhere.
For i = 2, ..., N − 1,
Z xi
Z xi
x − xi−1
xi − x
1
i
dx
+
a
(x)
dx
ai−1,i = a(φi−1
,
φ
)
=
−
a(x)
0
h
h
h2i
hi
hi
x
x
i−1
i−1
#
" Z
Z
xi
xi
ai−1,i =
1
h2i
ai,i =
−
xi−1
xi−1
a(φih , φih )
=
+
ai,i+1 =
a0 (x)(x − xi−1 )(xi − x)dx
a(x)dx +
a(φih , φi+1
h )
=
1
h2i+1
1
h2i
"Z
xi
Z
#
xi
2
a0 (x)(x − xi−1 ) dx
a(x)dx +
xi−1
Z
1
h2i+1
Z
−
xi−1
xi+1
Z
xi+1
a(x)dx +
xi
a0 (x)(x − xi+1 ) dx
2
xi
xi+1
Z
xi+1
a(x)dx +
xi
a0 (x)(x − xi+1 )(xi − x)dx
xi
a11 , a12 , aN N , aN −1,N are to be also computed.
For F = l(φih ) 1≤j≤N ,
Z xj+1 Z xj D
E Z xj+1
x − xj+1
x − xj−1
j
j
j
dx +
f
dx.
fj = l(φh ) = f, φh =
f φh dx =
f
hj
hj+1
xj
xj−1
xj−1
Remarks: 1. The finite element method (with constant mesh size) coincides with finite difference
method except for the fact that an average of f over (xj−1 , xj+1 ) is now used instead of point values of
f (xj ).
2. FE method is thus Galerkin method applied to a special choice of finite dimensional subspace, namely
continuous piecewise linear elements in this case.
Step 7: Error Estimates
We need the following results for obtaining the estimates ku − uh kV , ku − uh k0,I .
Cea’s Lemma: Let u and uh be the solutions of the weak formulation and (PGh ) respectively. Then
ku − uh kV ≤ C inf ku − uh kV
vh ∈Vh
4
(11)
Proof: From (2) and (7), we get ∀vh ∈ Vh ,
a(u, vh )
=
a(uh , vh )
=

l(vh ) 
⇒ a(u − uh , vh ) = 0
l(vh ) 
∀vh ∈ Vh .
(12)
(12) is called Galerkin orthogonality property, i.e., the finite element solution uh is the orthogonal
projection of the exact solution u onto Vh with respect to the inner product a(·, ·) .
αku − uh k2V
≤ a(u − uh , u − vh ) (By coercivity of a(·, ·))
= a(u − uh , u − vh ) + a(u − uh , vh − uh ) (using (12))
|
{z
}
=0
i.e., ku −
uh k2V
≤
M
ku − uh kV ku − vh kV (Boundedness of a(·, ·)).
α
In case u = uh , the result is trivial.
M
ku − vh kV ∀ vh ∈ Vh .
For u 6= uh , ku − uh kV ≤
α
M
⇒ ku − uh kV ≤ C inf ku − vh kV with C =
.
vh ∈Vh
α
Remark: uh is the best approximation of u in V w.r.t k · ka = a(·, ·)1/2 .
Interpolation Operator
¯ with v(0) = v(1) = 0, we define the interpolant operator
For a function v ∈ C 0 (Ω),
Πh : v ∈ H01 (0, 1) ,→ C 0 (0, 1) −→ Vh ,
(Πh v)(xi ) = vI (xi ) = v(xi ) 0 ≤ i ≤ (N + 1).
We are now interested to find ku − uI k0 and ku − uI k1 so that we can use this in (11) to estimate
ku − uh kV .
Approximation Properties of interpolation operator
R
1/2
1
Assume that u ∈ H 2 (I). We know |u|2 = 0 u00 (x)2 dx
is a semi-norm. By Sobolev imbedding
theorem, u ∈ C 1 (I). Applying mean-value theorem to u ∈ [xi−1 , xi ], we get ∃ξ ∈ (xi−1 , xi ) such that
u0 (ξ) =
u(xi ) − u(xi−1 )
hi
u(xi ) − u(xi−1 )
= u0 (ξ).
hi
Hence, in this interval, we have
In this interval, u0I (x) =
u0 (x) − u0I (x)
= u0 (x) − u0 (ξ)
Z x
=
u00 (t)dt
ξ
5
Applying Cauchy-Schwarz inequality, we have
Z x
1.u00 (t)dt
|u0 (x) − u0I (x)| = ξ
Z x 1/2 Z
12
≤
ξ
1/2
2
|u00 (t)| dt
ξ
!1/2
xi
Z
1/2
hi
≤
x
00
2
00
2
∀ x ∈ (xi−1 , xi )
(u (t)) dt
xi−1
0
i.e., |u (x) −
u0I (x)|
≤
!1/2
xi
Z
1/2
hi
∀ x ∈ (xi−1 , xi ).
(u (t)) dt
(13)
xi−1
To derive a bound on |u(x) − uI (x)|:
Z
u(x) − uI (x)
x
[u0 (t) − u0I (t)] dt
=
xi−1
Z
x
0
0
|u(x) − uI (x)| = (u (t) − uI (t)) dt
xi−1
Z x
≤
|u0 (t) − u0I (t)| dt
xi−1
Z
x
≤
Z
1/2
hi
i.e., |u(x) − uI (x)|
≤
00
2
(u (t)) dt
xi−1
3/2
hi
!1/2
xi
dx (using (13))
xi−1
Z
!1/2
xi
00
2
(u (t)) dt
∀ x ∈ (xi−1 , xi ).
(14)
xi−1
Squaring (13) and (14), integrating and
! summing over all elements and taking square roots, we get
Z xi
2
(u00 (t))2 dt
|u0 (x) − u0I (x)| ≤ hi
xi−1
!
Z
Z xi
n Z xi
X
2
0
0
00
2
|u (x) − uI (x)| dx =
hi
(u (t)) dt dx
I
i=1
2
xi−1
xi−1
2
⇒ |u − uI |1,I ≤ h2 |u|2,I , where h = max hi .
1≤i≤N
⇒ |u − uI |1,I ≤ h |u|2,I .
(15)
Similarly, we obtain |u − uI |0,I ≤ h2 |u|2,I from (14).
(16)
Error Estimates
(1) First we find an estimate for ku − uh kV = ku − uh k1,I .
ku − uh kV ≤ Cku − uI kV (using (11))
ku − uI k2V
= ku − uI k20,I + ku − uI k21,I
2
2
≤ h4 |u|2,I + h2 |u|2,I
=
2
2
(h4 + h2 ) |u|2,I = h2 (1 + h2 ) |u|2,I
2
⇒ ku − uI k2V ≤ 2h2 |u|2,I (using h < length(I) = 1)
... ku − uh kV ≤ Ch |u|2,I u ∈ H 2 (I) ∩ H01 (I).
6
(17)
(2) An estimate for ku − uh k0,I is found using Aubin-Nitsche’s duality argument.
Consider the following adjoint elliptic problem:
For g ∈ L2 (I), let φ be the solution of
−(a(x)φ0 )0 + a0 (x)φ = g,
x ∈ I,
(18)
φ(0) = φ(1) = 0
This problem satisfies the regularity result kφk2 ≤ Ckgk.
(19)
Multiplying (18) by e = u − uh and integrating over I, we get
a(e, φ) = (e, g)
(e, g) = a(e, φ − φI ) (using (12))
i.e., |(e, g)|
≤ CkekV kφ − φI kV
≤ ChkekV kgk0,I (using (19))
The result follows by choosing g = e.
i.e., kek20,I ≤ ChkekV kek0,I ⇒ ku − uh k0,I ≤ Ch2 |u|2,I
Remark:
(1) The analysis could have been done for a more general finite element space consisting of polynomials
of degree ≤ r − 1 where r ≥ 2. (The case discussed was for r = 2). Then, we would obtain
ku − uI k1,I ≤ Chr−1 |u|r,I and
ku − uI k0,I ≤ Chr |u|r,I under the assumption of additional regularity u ∈ H r (I).
(2) Computational Order of Convergence :
For 0 < h2 < h1 , let ku − uh1 k
≈ C(u)hα
1
ku − uh2 k ≈ C(u)hα
2
ku − uh1 k
h1
Then the order of convergence is α ≈ log
/ log( )
ku − uh2 k
h2
In the absence of exact solution, ‘u‘ may be replaced by a more refined computed solution.
7
An Elliptic Model Problem
Given f ∈ L2 (Ω), find 0 u0 such that


−∆u + u = f in Ω = (0, 1) × (0, 1) ⊂ R2 , 
(PC ) :  ∂u
|Γ = 0, Γ is the square boundary of Ω, 
∂n
• ∆u =
•
(20)
∂2u ∂2u
+ 2,
∂x2
∂y
∂u
is the normal derivative of 0 u0 in the direction of the exterior normal to Γ
∂n
• u =?? (unknown function).
Weak Formulation:
1
•Admissible space V = H (Ω) =
∂v ∂v
,
∈ L2 (Ω)
v : v ∈ L (Ω),
∂x ∂y
2
• Bilinear form a(·, ·) : V × V → R defined by
Z
Z
a(u, v) = (−∆u + u)v dΩ = (∇u.∇v + uv) dΩ
Ω
∀u, v ∈ V
(21)
(22)
Ω
is symmetric, continuous and V elliptic.
• Linear form l(·) : V → R defined by
Z
f v dΩ ∀ v ∈ V
l(v) =
(23)
Ω
is continuous.
The weak formulation corresponding to (20) can be defined as:

Given f ∈ L2 (Ω), find u ∈ V such that

a(u, v) = l(v)
∀v ∈ V.
(24)
Since a(·, ·) is a continuous, symmetric, coercive, bilinear form and l(·) is a continuous, linear form,
using Lax-Milgram lemma, the weak formulation is well-posed.
Remark: Since a(·, ·) is symmetric, we can define the energy minimization formulation also.
¯
Triangulation of Ω:
¯ ⊂ R2 , we mean a subdivision of Ω
¯ into closed trianBy a triangulation τh of a closed domain Ω
ELEM
gles/rectangles/quadrilaterals denoted by {Ti }N
such that
i=1
ELEM
¯ = ∪N
(i) Ω
Ti
i=1



 φ
(ii) Ti ∩ Tj =
common vertex



common side
(iii) Ti0 ∩ Tj0 = φ for i 6= j
8
for i 6= j
(25)
Remark: For simplicity, we shall assume that ∂Ω = Γ is a polygonal curve, i.e., Ω is a polygonal
domain. If Γ is curved, we may first approximate Γ with a polygonal curve.
We introduce the mesh parameter
h = max (diam T ),
T ∈τh
diam(T ) = longest side of T.
(26)
Galerkin Finite Element Problem:
¯ we associate a finite dimensional linear space Vh of continuous
To the triangulation τh of Ω,
piecewise polynomials of degree ≤ 1 in each triangle T of τh , i.e.,
¯ vh |T ∈ P1 (T ) ∀ T ∈ τh }
Vh = {vh : vh ∈ C 0 (Ω),
(27)
Remark: This is the simplest choice for Vh . Other choices of Vh are possible.
Note that
• Vh is finite-dimensional. (dimVh ≤ N ELEM ∗ dimP1 (T ))
• Vh is a subspace of V .
Proof: i.e., to show that vh ∈ L2 (Ω),
∂vh ∂vh
,
∈ L2 (Ω) ∀ vh ∈ Vh .
∂x ∂y
| {z }
(distributional derivatives)
...
¯ vh ∈ L2 (Ω).
vh ∈ C 0 (Ω),
0
0
vh ∈ L2 (Ω)
Z ⇒ vh ∈ Lloc (Ω) ⇒ ∃Tvh ∈ D (Ω) such that
vh φ dΩ ∀ φ ∈ D(Ω).
Tvh (φ) =
Ω
Z
Consider
T
i.e.,
X Z
T ∈τh
T
∂
(vh |T )φ dT
∂x
∂
(vh |T )φ dT
∂x
Z
Z
∂
∂φ
(vh |T φ) dT −
vh | T
dT
∂x
∂x
ZT
Z T
∂φ
=
vh |T φni ∂T −
vh | T
dT
∂x
∂T
T
Z
∂φ
= −
vh
dΩ
∂x
Ω
∂
=
(Tvh ), φ ∀ φ ∈ D(Ω)
∂x
=
i.e., the derivative of the distribution defined by Tvh is nothing but the distribution obtained from
patching up piecewise derivatives
of vh in each T ∈ τh .
X Z ∂
∂vh ∂Tvh
Denoting ∂x =
(vh |T ) dT,
= Th ∂vh i . Since vh is a polynomial in each T ∈ τh , we
∂x
∂x
∂x
T ∈τh T
!
get that
∂vh ∂vh
∂x , ∂y
(in distributional sense) ∈ L2 (Ω).
9
The Galerkin finite element problem (PGh ) can be defined as:

 Given f ∈ L2 (Ω), find u ∈ V such that
h
h
(PGh ) :
 a(uh , vh ) = l(vh ) ∀ vh ∈ Vh .
(28)
Matrix formulation: Choosing an appropriate basis {φih }i=1 for Vh (the construction of {φih }N
i=1 will
be explained with details), we can reduce (PGh ) into
KU = F
(29)
(proceeding exactly in the same way as in the 1D case),
where
• K = [a(φih , φjh )]1≤i,j≤N is a N × N matrix called the global stiffness matrix,
• F = [l(φjh )]1≤j≤N is called the global load vector
• U = [u1 , ..., uN ]t is the unknown vector.
N
X
ui φih (x, y) .
Here uh =
i=1
Construction of appropriate basis functions {φih }N
i=1 of Vh
To each global node ai (1 ≤ i ≤ 9), we associate a 0 pyramidal function0 φih (i.e., the graph of φih is a
pyramid of unit height over the triangles containing the node ai ) which vanishes in all those triangles
T not containing 0 a0i as one of its nodes, i.e.,










(i) φih (ai ) = 1, φih (aj ) = 0 ∀ i 6= j
(Height of the pyramid at ai is φih (ai ) = 1)
(ii) φih |T ∈ P1 (T ) ∀ T ∈ τh
(30)



(iii) Supp φih = ∪ai ∈T T (i.e., union of all triangles T containing ai ). 





i
0 ¯
(iv) φh ∈ C (Ω).
From (ii) and (iv), we get φih ∈ Vh .
For the explicit construction of φih , we now introduce the 0 Barycentric Co-ordinates0 for a
traingle.
Definition: The closed convex hull T defined by
T = {x : x ∈ R2 , x =
3
X
λi ai , 0 ≤ λi ≤ 1,
i=1
3
X
i=1
is called a 2-simplex in R2 with 0 30 vertices
{ai }3i=1
a
11 a12 a13 (such that D3 = a21 a22 a23 6= 0)
1
1
1 10
λi = 1}
(31)
∀(x1 , x2 ) ∈ T¯, (x1 , x2 ) 7−→ (λ1 , λ2 , λ3 ) is a one-one correspondence.
In particular, a1 7−→ (1, 0, 0), a2 7−→ (0, 1, 0), a3 7−→ (0, 0, 1).
3
X
Now x ∈ T ⇒ x = λ1 a1 + λ2 a2 + λ3 a3 s.t. 0 ≤ λi ≤ 1,
λi = 1.
i=1

i.e., 

x1

 = λ1 
x2
i.e.,
a11
a12


 + λ2 
a21
a22


 + λ3 
a31
a32


λ1 a11 + λ2 a21 + λ3 a31 = x1
λ1 a12 + λ2 a22 + λ3 a32 = x2
with
λ1 + λ2 + λ3
a
11 a21
This system has a unique solution because a12 a22
1
1
(32)
=1
a31 a32 6= 0.
1 Definition: Let T be a 2-simplex in R2 with three vertices {ai }3i=1 such that ∀ x ∈ T, x =
3
X
λi ai , 0 ≤
i=1
λi ≤ 1,
3
X
λi = 1. Then (λ1 , λ2 , λ3 ) are the barycentric co-ordinates of the point x of the 2 simplex
i=1
T ⊂ R2 .
Properties of barycentric co-ordinates:
(1) Barycentric co-ordinates of any point x ∈ T satisfies λ1 + λ2 + λ3 = 1.
(2) For any x ∈
/ T , a representation of points using barycentric co-ordinates is not
(3) To locate a point x ∈ T , we need
(i) Barycentric co-ordinates (λ1 , λ2 , λ3 ) with 0 ≤ λi ≤ 1,
3
X
λi = 1.
i=1
(ii) Cartesian co-ordinates of the vertices of T .
(4) The graphs of the barycentric co-ordinate functions can be given by:
Projecting onto the xy plane,
11
possible.
i.e., graph of λ1 is a plane with height 0 above L1 and height 1 above a1 .
...
a11 λ1 + a12 λ2 + a13 λ3
= x
a21 λ1 + a22 λ2 + a23 λ3
= y
λ1 + λ2 + λ3 = 1
Using Cramer’s rule, we obtain x a
12 a13 1
y a22 a23 = α0 + α1 x + α2 y ∈ P1 (T )
λ1 =
2(area T ) 2(area T )
1 1
1 where α0 = a12 a23 − a13 a22 , α1 = a22 − a23 , α2 = a13 − a12 .
Similarly, we can show that λ2 (x, y) ∈ P1 (T ), λ3 (x, y) ∈ P1 (T ).
(5) Barycentric co-ordinates functions λi (x, y) 1 ≤ i ≤ 3 form a canonical nodal basis for P1 (T ).
3
X
..
αi λi = 0 ⇒ αi = 0; hence they form a basis
( . λi (aj ) = δij , λi (x, y) ∈ P1 (T ) 1 ≤ i, j ≤ 3,
i=1
for P1 (T )).
p ∈ P1 (T ), p(x, y) =
3
X
p(ai )λi (x, y) ∈ P1 (T ).
i=1
(6) For the reference triangle Tb with vertices a1 = (1, 0), a2 = (1, 0) and a3 = (0, 0), the barycentric
co-ordinate functions take the following simplified form:


 λ1 (x, y) = x

(33)
λ2 (x, y) = y



λ3 (x, y) =
1−x−y
The barycentric co-ordinates in triangles will be used now to define global basis functions {φih }N
i=1 for
Vh which satisfy (i)-(iv) in (30).
¯ We will illustrate the construction of {φi }N for this
Consider the following triangulation of Ω.
h i=1
case.

 λT1 (x, y) in T
1
1
φ1h (x, y) =
 0
otherwise.
(We follow the local numbering for the vertices of a triangle. Some convention for the numbering has to
be followed consistently. In our example, we choose the vertex with right angle as a1 always (locally)
and move anti-clockwise to locate a2 and a3 in that order).

T

 λ2 1 (x, y) in T1



 λT2 (x, y) in T
2
3
2
φh (x, y) =
T3


 λ1 (x, y) in T3


 0
otherwise

T3


 λ2 (x, y) in T3
φ3h (x, y) =
λT3 4 (x, y) in T4



0
12
otherwise
and so on.
If we choose {φih }9i=1 in this way, we can verify that (i) - (iv) of (30) are satisfied.
Note that a(φih , φjh ) = 0 whenever nodes i and j donot belong to a same triangle. This makes
K = [a(φih , φjh )]1≤i,j≤N a sparse matrix. Also, K is positive definite as a(·, ·) is co-ercive.
Now, we will explain how the 0 global stiffness matrix0 K can be 0 assembled0 from 0 element stiffness
matrices0 KT and how the 0 global load vector0 F can be assembled using the element load vectors FT .
Construction of K and F : We will illustrate the construction of K and F for a simple case.
¯ = [0, 1] × [0, 1] given as:
Let τh be a triangulation of Ω
Note that, since we are solving a Neumann problem, we have no information about 0 u0 along the
boundary Γ.
The numbers 1,2,...,9 in the figure indicate the nodal points a1 , ..., a9 . As we have seen earlier, there
are 9 nodal basis functions {φih }9i=1 of Vh which can be constructed such that Supp {φih } is as small as
possible.
The matrix K = [a(φih , φjh )]9×9 will have the following structure:

×










K=










×
0
×
0
0
0
0
0
×
×
×
×
0
0
0
×
0
×
×
0
0
×
×
0
×
0
×
×
×
×
×
0
×
×
×

0 


0 


0 


0 


× 

0 


× 

×
×
9×9
Remark: In actual finite element computations, we need to store only those elements above the
diagonal (including diagonal elements and below the skyline structure, in case we are using a “SKYLINE
SOLVER”).
Let K = (kij )1≤i,j≤9 .
kij
=
=
Z "
#
∂φih ∂φjh
∂φih ∂φjh
i j
=
+
+ φh φh dΩ
∂x ∂x
∂y ∂y
Ω
"
#
X Z
∂φih ∂φjh
∂φih ∂φjh
i j
+
+ φh φh dΩ
∂x ∂x
∂y ∂y
Tij
a(φih , φjh )
Tij ∈τh
where Tij ’s are those triangles in τh which have both the nodes ai and aj .
13
For example,
a(φ1h , φ1h )
Z
=
T1
∂φ1h ∂φ1h
∂φ1h ∂φ1h
1 2
+
+ (φh ) dT1 since the node a1 is present only in
∂x ∂x
∂y ∂y
the triangle T1 .
T1
1
But in T1 , we know
Z φh = λ1 and hence
∂ T1 ∂ T1
∂ T1 ∂ T1
a(φ1h , φ1h ) =
(λ1 ) (λ1 ) +
(λ1 ) (λ1 ) + (λT1 1 )2 dxdy
∂x
∂y
∂y
T1 ∂x
Similarly,
Z ∂ T1 ∂ T1
∂ T1 ∂ T1
T1 T1
1
2
a(φh , φh ) =
(λ1 ) (λ2 ) +
(λ ) (λ ) + λ1 λ2 dxdy
∂x
∂y 1 ∂y 2
T1 ∂x
a(φ1h , φ3h ) = Z
0 asa1 and a3 donot belong to the same triangle. ∂ T1 ∂ T1
∂ T1 ∂ T1
1
4
a(φh , φh ) =
(λ1 ) (λ3 ) +
(λ ) (λ ) + λT1 1 λT3 1 dxdy
∂x
∂y 1 ∂y 3
T1 ∂x
a(φ1h , φjh ) = 0 ∀ j = 5, ..., 9.
#
2
2
Z "
∂ T2
∂ T2
T2 2
2
2
a(φh , φh ) =
(λ ) + [
(λ ) + (λ3 ) dxdy
∂x 3
∂y 3
T2
#
2
2
Z "
∂ T3
∂ T3
T3 2
(λ ) + [
(λ ) + (λ1 ) dxdy
+
∂x 1
∂y 1
T3
Thus [a(φih , φjh )] can be assembled from suitable element contributions of integrals of the kind
i
R h ∂ Tk ∂ Tk
Tk ∂
Tk
Tk Tk
∂
(λ
(λ
(λ
(λ
)
)
+
[
)
)
+
λ
λ
dT .
i
j
i
j
i
j
∂x
∂y
∂y
T ∂x
Similarly F = [l(φjh )]1≤j≤N can be also assembled by adding up suitable element contributions (from
triangles which contain the node aj ) i.e., in the implementation, we first compute the element stiffness
T
matrices (kij
) ∀ T ∈ τh and the element load vectors (f )T ∀ T ∈ τh . Then we assemble these element
matrices to obtain the 0 global stiffness matrix0 K and the global load vector 0 F 0 .
We will illustrate this 0 assembly procedure0 for our example:
For each triangle T ∈ τh , we compute the element stiffness matrix and the element load vector as
follows:

T
k11
T
k12
T
k13


f1T




T  and f T = (f T )
 T 
k23
j 1≤j≤3 =  f2 .

T
f3T
symmetric
k33
T
T
T
Local canonical basis functions in T are λ1 , λ2 , λ3

T
k T = (kij
)3x3 = 

T
k22
where,
T
kij
=
fiT =
Z ZT
∂ T ∂ T
∂
∂
(λi ) (λj ) + [ (λTi ) (λTj ) + λTi λTj dT
∂x
∂x
∂y
∂y
f λTi dT
1 ≤ i ≤ 3.
T
14
1 ≤ i, j ≤ 3.
Now,













K=











T1
k11
T1
k12
T
T
k221 +k332
T
+k113
0
T1
k13
0
0
0
0
T3
k12
T2
T1
+ k23
k23
T3
T2
+ k13
k13
0
0
0
0
T3
T4
k23
+ k23
T4
k13
0
0
0
T5
k13
0
T4
T7
k12
+ k12
T5
T6
k23
+ k23
T6
T7
k13
+ k13
0
T8
T7
+ k23
k23
T3
T4
k22
+ k33
T
T
k331 +k222
T5
+k11
T2
k12
+
T5
k12
T2
T3
k11
+k33
T4
T5
+k22
+k22
T
T
+k336 +k117
T
T
k114 +k222
T
+k338
T5
T6
k33
+ k22
T6
k12
T
T
k116 +k337
T
+k228
0


0 


0 



0 



0 


T8 

k13


0 


T8 
k12

T9
k11
9×9











F =









f1T1

f2T1 + f3T2 + f1T3




















f2T3 + f3T4
f3T1 + f2T2 + f1T5
f1T2 + f3T3 + f2T4 + f2T5 + f3T6 + f1T7
f2T7 + f3T8
f3T5 + f2T6
f1T6 + f3T7
f1T9
9×1
KU = F is now solved using direct/iterative methods. U gives the solution uh at the nodal points
(which are vertices of the triangles in the triangulation).
Error estimates:
The important steps involved are:
(1) Cea’s lemma (stated earlier in one-dimensional case).
(2) Interpolation operator.
(3) Affine mapping,(Inverse),|J(FT )|, |J(FT−1 )| estimates.
(4) Change of scale.
(5) Bramble-Hilbert lemma.
(6) Semi-norm inner product.
(7) Finite element interpolation over reference triangle.
(8) Estimates k u − uh kV , k u − uh k0,Ω (using Aubin-Nitsche duality arguments as in one-dimensional
case)
Now a detailed description of the steps are given :
15
Step (2):Interpolation operator
Define a map:
Πh : V ∩ C 0 (Ω) −→ Vh
u 7−→ uI
by
(Πh u)(x1 , x2 ) = uI (x1 , x2 ) =
N
X
−
u(→
ai )φih (x1 , x2 )
(34)
i=1
−
where {→
ai }N
i=1 are the nodes of the triangulation.
−
−
(Πh u)(→
ai ) = u(→
ai )
1 ≤ i ≤ N.
(35)
Step(3): Affine mapping
Let Tb represent reference triangle and T , any triangle of τh . Then, define a map:
FT : Tb −→ T by
b + [bT ]2×1
FT (b
x) = [BT ]2×2 x
where BT = (bi,j )2×2 and bT = (bi )2×1 are obtained using :
−
FT (abi ) = →
ai 1 ≤ i ≤ 3.
If FT (b
x, yb) = (x, y) then, in matrix form:

 
x
b

 =  11
y
b21
b12
b22


x
b
yb


+
b1
b2

(36)

−
Using: abi −→ →
ai , we get:

BT = 
b11
b12
b21
b22


 =
16
a11 − a13
a12 − a13
a21 − a23
a22 − a23



bT = 
Jacobian of the transformation FT = J(FT )
|J(FT )| = Abs from equation (36):

 
b
x
b

 =  11
yb
b21
where:
BT−1

b22

a23
will be:
∂x
∂b
x
∂y
∂b
x
−1 
b12

a13

∂x
∂b
y
∂y
∂b
y
x − b1
y − b2

1  b22
=
|BT |
−b21
b11
= Abs b21

 =
−b12
b11


b12 b22 b∗11
b∗12

b∗21
b∗22


 =
b∗11
b∗12

b∗21
b∗22

x − b1
y − b2


Jacobian of transformation FT−1 = JF −1 will be:
T
x
∂b
∂x
|J(FT −1 )| = Abs ∂b
y
∂x
∂b
x
∂b
y
∂b
y
∂y
∗
b11
= Abs b∗21
b∗12 b∗22 Thus clearly,
|J(FT−1 )| = Abs|[J(FT )]−1 |
Properties of Affine mapping FT
• FT is invertible.
• FT (Tb) = T
• FT carries vertices and sides of Tb to vertices and sides of T
We would need the following estimates of bij , b∗ij , J(FT ) and JF −1 :
T
Consider a triangle T . Let α1 , α2 , α3 denote angles of T and h1 , h2 , h3 denote the sides opposite to
these angles.
Let
α1 ≤ α2 ≤ α3 ⇒ h1 ≤ h2 ≤ h3
Let α = min{αi } = α1 and h = max{ hi } = h3
Now, h1 + h2 > h3 and h2 > h1
⇒ h2 >
17
h3
2
(37)
h2 h3 sinα1
2
h2 sinα
h2 h3 sinα1
h2 sinα
⇒
≤
≤
. Thus,
4
2
2
Area of triangle T =
h2 sinα
≤ |J(FT )| ≤ h2 sinα
2
1
h2 sinα
2
≤ |J(FT−1 )| ≤
(since|J(FT−1 )| =
h2 sinα
(38)
1
)
|J(FT )|
(39)
Again, |b11 | = |a11 − a13 | ≤ h2 ≤ h Similarly,
|bij | ≤ h ∀ 1 ≤ i, j ≤ 2
|b∗ij | =
(40)
2h
|bij |
≤ 2
|J(FT )|
h sinα
Thus,
|b∗ij | ≤
2
hsinα
(41)
Step (4) :Change of scale
(We will derive the estimates in reference triangle Tb and then do a change of scale to get back to T .)
Theorem:
Let
• Tb be the reference triangle.
• T be any generic triangle of τh which is obtained from Tb by the affine mapping FT defined above.
• u(x, y) ∈ H s (T ), u
b(b
x, yb) ∈ H s (Tb) satisfy
u
b(b
x, yb) = u(x, y), (x, y) ∈ T, FT (b
x, yb) = (x, y)
Then, there exists constants Cs , Ds > 0 (independent of u, h and α) such that:
Ds (sinα)s−1/2 hs−1 |u|s,T ≤ |b
u|s,Tb ≤ Cs (sinα)−1/2 hs−1 |u|s,T
(42)
Proof :
(for s=0): To prove that
D0 (sinα)−1/2 h−1 |u|0,T ≤ |b
u|0,Tb ≤ C0 (sinα)−1/2 h−1 |u|0,T
We know that
Z
2
Z
|u(x)| dT =
T
|b
u(b
x)|2 |J(FT )|dTb
Tb
⇒ |u|20,T = |b
u|20,Tb |J(FT )|
⇒ |J(FT )|−1 |u|20,T = |b
u|20,Tb
Using equation (39):
2|u|20,T
|u|20,T
2
≤
|b
u
|
≤
0,Tb
h2 sinα
h2 sinα
√
|u|0,T
2|u|0,T
⇒ √
≤ |b
u|0,Tb ≤ √
h sinα
h sinα
18
(43)
√
i.e., (43) holds with D0 = 1, C0 =
2.
(for s=1): To prove that ∃ D1 , C1 such that:
D1 (sinα)1/2 h0 |u|1,T ≤ |b
u|1,Tb ≤ C1 (sinα)−1/2 h0 |u|1,T
Z
|u|21,T = (u2x + u2y )dT
(44)
(45)
T
We know u(x, y) = u
b(b
x, yb) i.e.
ux = u
bxbx
bx + u
bybybx
(46)
uy = u
bxbx
by + u
bybyby
(47)
Substituting (46) and (47) in (45) :
|u|21,T =
Z
Tb
[c1 u
b2xb + c2 u
b2yb + 2c3 u
bxbu
byb]|J(FT )|dTb
where
c1 = x
b2x + x
b2y , c2 = ybx2 + yby2 , c3 = x
bx ybx + x
by x
bx yby .
Using ab ≤
a2 + b2
in (48),
2
|u|21,T
Z "
≤
c1 u
b2xb + c2 u
b2yb + 2c3
u
b2xb + u
b2yb
2
Tb
Z
=
Tb
!#
|J(FT )|dT
[(c1 + c3 )b
u2xb + (c2 + c3 )b
u2yb)]|J(FT )|dT
Z
≤ max{c1 + c3 , c2 + c3 }
Tb
(b
u2xb + u
b2yb)|J(FT )|dT
⇒ |u|21,T ≤ C |b
u|21,Tb |J(FT )|
Now, C = max{c1 + c3 , c2 + c3 }, where :
c1 + c3 ≤
16
16
, c2 + c3 ≤ 2 2 .
h2 sin2 α
h sin α
Using (38),
|u|21,T ≤
⇒
16
|b
u|2
sinα 1,Tb
1
(sinα)1/2 |u|1,T ≤ |b
u|1,Tb
4
Thus, LHS of (44) holds with D1 = 1/4.
Similarly, we can get the right hand side estimate :
|b
u|21,Tb ≤ |J(FT−1 )|C 0 |u|21,T
where:
C 0 = max{c01 + c02 , c03 + c02 }, c01 = x2xb + x2yb, c02 = yxb2 + yyb2 , c03 = xxbyxb + xybyyb.
19
(48)
From (40), C 0 ≤ 4h2 Using (39),
8
|u|2
sinα 1,T
√
8
|u|1,T
≤√
sinα
|b
u|21,Tb ≤
⇒ |b
u|1,Tb
Thus, RHS of (44) also holds.
To establish equation (44) for s > 1 we need an analog of equation (46) and equation (47) for derivatives
of order 0 s‘. Note that every partial derivative of u of order s is a polynomial of degree s w.r.t.
x
bx , x
by , ybx , yby with no term having degree less than s. Summing the squares of all sth order derivatives
a2 + b2
and using ab ≤
to eliminate cross-products, we get:
2
|u|2s,T ≤ |J(FT )|C|b
u|2s,Tb
where C = O((hsinα)−2s ).
Similarly, reversing the process we get:
|b
u|2s,Tb ≤ |J(FT−1 )|C 0 |u|2s,T ,
where C 0 = O(h2s ). Using upper bounds for |J(FT )| and |J(FT−1 )|, we get the required result.
Step (5): Bramble-Hilbert Lemma
Let
• Ω be a sufficiently smooth domain
• L : H k+1 (Ω) −→ R is a bounded linear funcional,
i.e., |L(u)| ≤ M k u kk+1,Ω ∀u ∈ H k+1 (Ω), M > 0
• L(p) = 0 ∀p ∈ Pk (Ω).
Then, there exists a constant C > 0 depending only on Ω such that:
|L(u)| ≤ CM |u|k+1,Ω ∀u ∈ H k+1 (Ω)
(49)
Step (6): Semi-norm inner product
We will define an inner product which will give rise to a semi-norm.
For u ∈ H k+1 (Ω), (k ≥ 1)
[u, v]s,Ω =
X Z
|α|=s
Dα u(x)Dα v(x)dΩ
Ω
Now, [u, u]s,Ω = |u|2s,Ω .
Step (7): Finite element interpolation over the reference triangle
Given some integer k ≥ 1 and some u
b ∈ H k+1 (Tb), consider the interpolant :
u
bI (b
x, yb) =
n
X
u
b(b
ai )φbi (b
x, yb) , (b
x, yb) ∈ Tb
i=1
20
(50)
where b
ai are nodes in the reference triangles and n denotes the number of nodes in a triangle of the
triangulation. In this case n = 3, φbi (b
x, yb) denotes the element nodal canonical basis functions in Tb.
Differentiating (50),
Dα u
bI (b
x, yb) =
n
X
u
b(b
ai )Dα φbi (b
x, yb) , (b
x, yb) ∈ Tb.
i=1
Since the standard local basis functions in Tb are polynomials, derivatives of higher order vanish. So,
∃M > 0 such that :
|Dα φbi (b
x, yb)| ≤ M ∀(b
x, yb) ∈ Tb , i = 1, 2, ..., n
Also, |b
u(b
ai )| ≤ max(bx,by)∈Tb |b
u(b
x, yb)| ≤ h k u
b kk+1,Tb ∀b
u ∈ H k+1 (Tb)
(Because H k+1 ,→ C 0 as k ≥ 1).
Here h depends on k and Tb. Hence,
|Dα u
bi (b
x, yb)| ≤ ηnM k u
b kk+1,Tb ∀b
u ∈ H k+1 (Tb).
c indeForming the Sobolov norms of u
bI (of which only a finite number are non-zero),we see that ∃M
pendent of u
b and s, such that:
cku
|b
uI |s,Tb ≤ M
b kk+1,Tb ∀b
u ∈ H k+1 (Tb) , s = 0, 1, 2, ...
(51)
Theorem:
Let
• Tb be the reference finite element.
• {φbi }ni=1 be the standard local basis functions on Tb.
• Let k be the largest integer such that Pk (Tb) ⊆ SP AN {φbi }n
i=1
b independent of u
Then, there exists a constant C,
b such that
b u|
|b
u−u
bI |s,Tb ≤ C|b
u ∈ H k+1 (Tb) s = 0, 1, ..., k + 1.
k+1,Tb ∀b
Proof :
Let 0 ≤ s ≤ k + 1. Fix vb ∈ H s (Tb).
Let L : H k+1 (Tb) −→ R be a linear functional defined by:
L(b
u) = [b
u−u
bI , vb]s,Tb
also |L(b
u)| = |[b
u−u
bI , vb]s,Tb |
≤ |b
u−u
bI |s,Tb |b
v |s,Tb
≤ (|b
u|s,Tb + |b
uI |s,Tb )(|b
v |s,Tb )
c)|b
≤ (1 + M
v |s,Tb k u
b kk+1,Tb .
So, L is a bounded linear functional on H k+1 (Tb).
If
u
b ∈ Pk (Tb) ⇒ u
bI = u
b
21
(52)
⇒ L(b
u) = [b
u−u
bI , vb]s,Tb = 0
So, by Bramble-Hilbert lemma , ∃ a contant C such that
c)|b
|L(b
u)| ≤ C(1 + M
v |s,Tb |b
u|k+1,Tb ∀b
u ∈ H k+1 (Tb).
Put vb = u
b−u
bI to get:
c)|b
|b
u−u
bI |2s,Tb ≤ C(1 + M
u−u
bI |s,Tb |b
u|k+1,Tb
b u|
|b
u−u
bI |s,Tb ≤ C|b
k+1,Tb
∀b
u ∈ H k+1 (Tb), s = 0, 1, ..., k + 1.
The global basis functions for Vh belong to C 0 (Ω). Hence, it is typical that the first order derivative
of global basis functions and of the interpolant uI , have step discontinuties across element boundaries
with the result that,in general uI does not belong to H 2 (Ω). Hence, the Sobolev-seminorm |u − uI |s
must be restricted to s = 0, 1 (In the earlier theorem,this restriction was not necessary, as the domain
here is a single element.)
Theorem:
Let
• Ω be a polygon in R2
• τh be a triangulation of Ω into closed triangles.
• The associated canonical nodal basis functions be {φi }N
i=1 .
• h, α denote the greatest edge length and the smallest angle respectively in the mesh.
• k be the largest integer for which Pk (Ω) ⊆ Vh = SP AN {φ1 , φ2 , ..., φN }.
Then,∃ a value C , independent of u and the mesh,such that
|u − uI |s ≤ C(sinα)−s hk+1−s |u|k+1 ∀u ∈ H k+1 (Ω) s = 0, 1,
(53)
where uI is the interpolant of u.
Proof :
We apply (42) (change of scale) twice and (52)(Finite element interpolation estimate) once to get (53).
Let T ∈ τh where h and α are the greatest edge length and smallest angle, respectively of that element.Then,
|u − uI |2s,T ≤ Ds−2 (sinα)1−2s h2−2s |b
u−u
bI |2s,Tb
b 2 |b
≤ Ds−2 (sinα)1−2s h2−2s C
u|2k+1,Tb
2
b 2 Ck+1
≤ Ds−2 (sinα)1−2s h2−2s C
(sinα)−1 h2k |b
u|2k+1,Tb
(using equation (42) with s=0 or 1)
2
b 2 |b
≤ Ds−2 CK+1
(sinα)−2s h2k+2−2s C
u|2k+1,Tb
Summing over all elements, we get
|u − uI |2s,Ω ≤ C 2 (sinα)−2s h2k+2−2s |u|2k+1,Ω
22
2
b 2 /min{D0 , D1 } = 4CC
b k+1 i.e.
where C 2 = Ck+1
C
|u − uI |s,Ω ≤ C(sinα)−s hk+1−s |u|k+1,Ω ∀u ∈ H k+1 (Ω) , s = 0, 1.
k u − uh kV can be computed using this. We get
k u − uh kV ≤ Ch|u|2,Ω ∀u ∈ H 2 (Ω).
23