Finite Element Method

Lecture notes for the course

‘Computational mechanics – numerical methods for fluids and solids’
(4MC10)

Ron Peerlings

2014/08/04
 Contents

1 Introduction 1

1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2 Characteristics of the finite element method . . . . . . . . . . . . . . . . . . . . . . . . 4

1.3 Objectives and contents of the course . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2 Finite element method for one-dimensional boundary value problems 7

2.1 Model boundary value problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2 Discretisation of the problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.3 Finite element shape functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.4 Assembly and solving of the linear system . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.5 Isoparametric elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3 Two and three-dimensional scalar boundary value problems 31

3.1 Model boundary value problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.2 Finite element discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.3 Two-dimensional elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

3.4 Three-dimensional elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

4 Vectorial boundary value problems 57

4.1 Model boundary value problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

4.2 Finite element discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

 1

Introduction
1.1 Background

Modern engineering design makes more and more use of advanced analysis methods to predict product
properties such as strength, dynamic behaviour, heat resistance, energetic efficiency, environmental im-
pact, etc. If reliable predictions of such properties can be obtained, more design options can be examined,
and a better final design may thus be obtained. At the same time, the need for prototyping and testing is
reduced, which may result in considerable cost savings and a shorter time-to-market. In a similar fash-
ion, analysis methods are used in the design of manufacturing processes, where they may dramatically
decrease the cost of tooling, and in the design and operation of process installations, where they assist in
controlling product properties and meeting environmental and safety standards.

Many of the problems that need to be solved in design analyses take the form of a balance of one or
several physical quantities. For example, the balances of mass and of momentum govern the flow of flu-
ids, the momentum balance governs the static and dynamic mechanical response of solids and the energy
balance governs the conduction and convection of heat. For continuous media, these balances are usually
formulated mathematically as (sets of) partial differential equations. Examples are the Navier–Stokes
equation in fluid mechanics, the equations of motion of solid continua and the convection–diffusion equa-
tion for temperature distributions. Combinations of these and other equations governing different, but
possibly coupled, physical phenomena can also be encountered. Associated with the above partial dif-
ferential equations are boundary conditions and, for instationary problems, initial conditions. Boundary
conditions describe how the physical fields within the considered domain interact with its surround-
ings. Examples are no-slip conditions and outflow pressures in fluid mechanics, kinematic restraints
and surface stresses in solid mechanics and fixed temperatures and insulation conditions in heat transfer
problems.

The combination of a set of partial differential equations and appropriate boundary conditions is called a
boundary value problem. Properly formulated – or in mathematical terms: well posed – boundary value
problems have a single, unique solution. Analytical expressions for these solutions, however, usually can
only be obtained for relatively simple equations defined on elementary geometries and therefore are of
little use for practical purposes. As a consequence, engineers facing boundary value problems on com-
plex domains must often rely on approximate solutions of these problems, which can usually be obtained
only in numerical form, i.e., for a given set of parameter values. Fortunately, the rapidly increasing power
of digital computers can be employed in generating such numerical solutions, allowing to solve relevant
problems in engineering practice with sufficient accuracy and within a reasonable amount of time.

A number of methods have been developed which, using the strengths of digital computing technology,
allow to efficiently compute numerical solutions of boundary value problems to within an arbitrary ac-
curacy (only limited by the computer hardware or software). Each of these methods uses some kind of
approximation to reduce the original, continuous equations to a discrete set of equations. For instance,
the finite difference method approximates the differential operators of the boundary value problem by
differences on a grid of points with a finite spacing. A set of algebraic equations is thus obtained, in
which the solution values in the grid points act as unknowns. Denser grids lead to more equations and a
higher accuracy of the approximate solution. Fourier methods and other spectral techniques approximate
the solution by a truncated series of orthogonal (harmonic) functions. The coefficients of these functions
become the unknowns in a finite set of algebraic equations. Here the accuracy of the solution can be
improved by increasing the number of terms after which the series is truncated. But perhaps the most
versatile and the most widely used among these methods is the finite element method, which is the subject
of this course.

3
1.2 Characteristics of the finite element method

Although originally developed in the context of solid mechanics, the finite element method has found
widespread application within the mechanics of both solids and fluids, as well as in many other fields of
applied science. Perhaps the most important reason for this popularity is its ability to deal with complex
geometries in a relatively straightforward way – a feature which is indispensable in engineering practice.
This geometrical flexibility is the result of the way in which the underlying continuum problem is dis-
cretised. The domain on which the problem is defined is described by a – possibly large – number of
volumes of elementary shape called finite elements or simply elements, see Figure 1.1 for an example.
Each of these elements contains a number of control points, or nodes, which uniquely define the posi-
tion of the element boundaries. Nodes on the boundary between adjacent elements are shared by these
elements, thus creating a connected grid of elements which we call the mesh. Within the element, the
quantities of interest are interpolated between their values in the nodes by polynomial interpolation func-
tions, the polynomial degree of which matches the number of nodes associated with the element. The
approximate solution of the problem is thus entirely controlled by the discrete set of nodal values of the
primary variables, also called the degrees of freedom. In order to determine these values, a discrete set of
equations is set up which provides the approximate solution that best satisfies – in some sense which will
be detailed later – the boundary value problem for the given discretisation. Refining this discretisation,
i.e. increasing the number of elements and nodes, results in a larger number of equations and in a more
accurate approximation of the continuous solution.

Figure 1.1: Three-dimensional finite element discretisation of a structural part of a rollercoaster;

courtesy of Femtec BV

The mathematical operations which are necessary to set up and solve the discrete set of equations for a
given finite element discretisation are relatively simple and can therefore be performed automatically and
efficiently by computer codes. The same holds to a large degree for generating the finite element mesh:
most finite element codes nowadays provide tools which automatically generate elements and nodes on
a geometry which can be defined using a graphical interface. Nevertheless, constructing geometric and

4
finite element models may take a considerable amount of time (up to months) for complex products or
structures. Computers, and in particular computer graphics, also play an important role in visualising
the results of finite element analyses in a way which is intuitive to the engineer. Here, again, most
finite element codes provide extensive graphical capabilities to present results. An example is given in
Figure 1.2.

Figure 1.2: Acoustic pressure field emitted from a gear box as computed using the finite element
method; courtesy of Paragon Numerical Engineering

1.3 Objectives and contents of the course

Operating finite element computer software and interpreting results obtained with it requires a certain
amount of knowledge and experience. Providing knowledge about a particular code and obtaining ex-
perience in using it, however, is not the main purpose of the present course. Instead, we focus on the
theory and computational techniques on which such codes are based and which are often largely hidden
from everyday users. We believe that a detailed understanding of the basic concepts of the finite element
method, its strengths and its limitations is essential in order to obtain and interpret finite element results
from any code in a responsible way. In this course we therefore start from the mathematical statement of
a number of model problems and systematically work our way towards the complete computational solu-
tion algorithm. For some problems these algorithms are implemented entirely in a simple finite element
code; for more complicated cases only relevant parts will be considered in detail.

The contents of the course and of these lecture notes can be divided roughly into three parts. In the first
part, consisting of Chapter 2, the basic concepts and theory of the finite element method are explained
for the case of a single, scalar linear differential equation. Although not very relevant from a practical
point of view, the one-dimensional case provides much insight and is easier to visualise. In the second
part (Chapter 3), the ideas developed in Chapter 2 are extended to two and three dimensions for scalar
linear partial differential equations. Vectorial linear partial differential equations, finally, are considered
in Chapter 4.

5
 2

Finite element method for one-dimensional

boundary value problems
2.1 Model boundary value problem
The basics of the finite element method are demonstrated most clearly in one dimension. This is done in
this chapter for a model problem which consists of a linear second-order differential equation combined
with one essential and one natural boundary condition. The remainder of this chapter explains the neces-
sary steps to discretise this particular problem in full detail, right up to the point where the result can be
programmed in a computer code. At the same time, however, this discussion shows the strategy that can
be followed for a wider class of problems.

The ordinary differential equation which we will use as a model equation throughout this chapter reads
 
d du
c(x) + f (x) = 0 (2.1)
dx dx
where the functions c(x) and f (x) must satisfy certain continuity requirements but are otherwise arbi-
trary. The solution u(x) clearly must be differentiable and the product c(x) du/dx in principle must be
differentiable as well. Equation (2.1) is a linear equation because it features only linear combinations
of the unknown function u(x) and its derivatives. This equation arises frequently in applied physics and
engineering. It describes for instance the longitudinal deformation of a bar (in terms of the axial dis-
placement u(x)) with stiffness c(x) as a result of a distributed force f (x). Or if u(x) is interpreted as
a temperature, it describes the stationary heat conduction in a medium with conductivity c(x) resulting
from a heat source f (x).

The differential equation (2.1) by itself does not have a unique solution. In order to obtain a well-posed
mathematical problem, it must be complemented by boundary conditions. Since we are dealing with a
second-order differential equation, two boundary conditions are needed: one at each end of the domain
on which (2.1) is valid. This domain is defined here as (0, L) and the following boundary conditions are
specified:
u(0) = g0 (2.2)
and
du
c(L) = qL (2.3)
dx x=L
where g0 and q L are constants. These conditions are of the types which are most frequently encountered.
The first condition is a so-called essential boundary condition, sometimes also referred to as a boundary
condition of the Dirichlet type. In terms of the deformation of a bar it indicates a fixed displacement,
imposed for instance by the grip of a tensile tester. If u(x) represents a temperature, boundary condition
(2.2) indicates a prescribed temperature, which might for instance model the presence of a large volume
of highly conductive material at this end of the problem domain. Condition (2.3) is a natural boundary
condition or a boundary condition of the Neumann type. For a deforming bar this condition represents a
force acting at the end of the bar. In other problems it usually represents a prescribed flux of a quantity
related to u(x) across the boundary, e.g. a heat flux in heat conduction problems.

The combination of ordinary differential equation (2.1) and boundary conditions (2.2)–(2.3) is called a
boundary value problem. It is emphasised that the present boundary conditions have been selected solely
to illustrate the treatment of the conditions which are encountered most frequently in practice. Other
combinations of the above types of boundary conditions as well as other types of boundary conditions
may also lead to well-posed boundary value problems. For instance, the natural condition at x = L may
be replaced by an essential condition similar to (2.2) if the physical problem at hand fixes the solution
u(x) at this position as well. Note, however, that imposing only natural conditions does not suffice to
uniquely define the solution, because a constant can always be added to it in such a case without violating
the differential equation or the boundary conditions.

9
2.2 Discretisation of the problem

A straightforward way to discretise the boundary value problem given by Equations (2.1)–(2.3) is to
replace the differential operators in (2.1) by quotients of finite differences between a number of sample
points. This provides a set of linear equations in terms of the sample point values. Solving the system
of equations, also taking into account the boundary conditions, results in a discrete approximation of
the solution of the boundary value problem. This method is known as the finite difference method. It
is sometimes called a direct discretisation method because it operates directly on the partial differential
equation. The class of approximation methods which we will consider in this course, on the contrary,
are indirect in the sense that they are based on an integral formulation of the boundary value problem.
Rewriting a boundary value problem in an integral form and discretising this form requires taking a
number of steps which are more or less identical for different problems. We will discuss each of these
steps in detail below for the model problem which was introduced in the previous section.

2.2.1 Weighted residuals form

First step in deriving a discretised formulation of our boundary value problem is rewriting the differential
equation in an integral form. Several methods can be used for this purpose. Here we use the weighted
residuals method. It states that the differential equation (2.1) is fully equivalent to the requirement that
the integral relation

ZL    
d du
φ(x) c(x) + φ(x) f (x) dx = 0 (2.4)
dx dx
0

holds for any φ(x). The function φ(x) is called a test function or weighting function and must satisfy
certain continuity requirements. Equation (2.4) is called the weighted residuals form of (2.1). It is
obtained by multiplying (2.1) by φ(x) and integrating on the domain (0, L).

That (2.4) is equivalent to (2.1) can be seen as follows. Suppose Equation (2.1) is not satisfied identically,
i.e. the residue
 
d du
r(x) = c(x) + f (x) (2.5)
dx dx

does not vanish for all x ∈ (0, L). A test function φ(x) can then always be constructed such that the
integral in (2.4), which can be rewritten as

ZL
φ(x)r(x) dx (2.6)
0

does not vanish. For instance, if φ(x) is set equal to r(x), this integral becomes

ZL
r 2 (x) dx (2.7)
0

which clearly will always be positive if r(x) differs from zero anywhere on the domain. The only way
this integral can vanish is therefore if r(x) = 0 for all x, which implies that (2.1) holds.

10
2.2.2 Weak form

It was mentioned already in Section 2.1 that the presence of a derivative of the product c(x) du/dx in
Equation (2.1) implies that c(x) du/dx must be differentiable. For smooth (e.g. constant) c(x) this in turn
requires u(x) to be twice differentiable. Such a strong continuity requirement would limit us considerably
in constructing approximate solutions of the boundary value problem. For this reason, Equation (2.4) is
rewritten in a form which imposes less stringent continuity requirements on u(x) before introducing an
approximate solution. The first term within the brackets in (2.4) is integrated by parts for this purpose:

ZL     L ZL
d du du dφ du
φ(x) c(x) dx = φ(x)c(x) − c(x) dx (2.8)
dx dx dx 0 dx dx
0 0

In the boundary contribution which appears, the product c(x)du/dx can be recognised as given to be
equal to q L at x = L by the natural boundary condition (2.3). At the other boundary (x = 0) we set this
product equal to −q0 :
  L
du du du
φ(x)c(x) = φ(L)c(L) − φ(0)c(0) = φ(0) q0 + φ(L) q L (2.9)
dx 0 dx x=L dx x=0

Contrary to q L , the boundary flux q0 is not known yet, since no natural boundary condition was set at
x = 0. However, this distinction is irrelevant in the next steps of our development, which are therefore
also valid for other combinations of boundary conditions. We will use the boundary data only at the end
of the process, after discretisation of the problem, at which point it will also appear that the unknown
boundary flux q0 can be retrieved as a byproduct of the approximate solution.

Now, substitution of (2.8) and (2.9) into (2.4) followed by reordering yields

ZL ZL
dφ du
c(x) dx = φ(x) f (x) dx + φ(0) q0 + φ(L) q L (2.10)
dx dx
0 0

which must again be satisfied for all possible φ(x). This equation no longer contains a second-order
derivative and thus does not require u(x) to be twice differentiable. Indeed, since the derivative du/dx
appears only in the integrand, the continuity requirements imposed upon u(x) can be relaxed even further:
u(x) must be continuous and needs to be only piecewise differentiable. Because of this less stringent
continuity requirement (2.10) is called the weak form of the differential equation (2.1). In contrast, (2.1)
is called the strong form. The price that has been paid for reducing the conditions on u(x) is a slightly
stronger continuity requirement on the test functions φ(x) compared with (2.4): these functions must
now be continuous and piecewise differentiable as well. However, this price is quite affordable because
this requirement on φ(x) is still no more stringent than that on u(x).

2.2.3 Discretisation of the approximate solution

Equation (2.10) provides the basis to generate approximate solutions u h (x) of our boundary value prob-
lem. For this purpose, we assume that u h (x) can be written as a linear combination of a finite number of
shape functions Ni (x):
n
X
h
u (x) = u i Ni (x) (2.11)
i=1

11
The shape functions are known functions of the coordinate x which can be selected more or less ar-
bitrarily. In order for the linear combination (2.11) to be unique, they should be linearly independent.
Examples of classes of functions which satisfy this requirement are polynomials of increasing degree or
harmonic functions of increasing frequency. The finite element method provides its own, special way to
generate shape functions, which we will discuss in detail in Section 2.3.

The coefficients u i in (2.11) are yet unknown. We call these coefficients the degrees of freedom. Note
that they are a finite number – n to be precise – of variables which do not depend on x. For a given
set of shape functions Ni (x), their values completely determine the approximate solution u h (x). Rather
than by a function value for each value of the continuous coordinate x, our approximate solution is now
determined by n discrete parameters u i ; we say that we have discretised the solution. The problem of
finding a solution has been reduced to finding the combination of values u i for which (2.11) provides the
best possible approximation of the exact solution u(x) in some sense.

For the following discussion it is convenient to rewrite expression (2.11) in a matrix form. Storing the
degrees of freedom u i and the shape functions Ni (x) in column matrices
h iT h iT
u = u1 u2 u3 . . . un N(x) = N1 (x) N2 (x) N3 (x) . . . Nn (x) (2.12)
˜ ˜
we can write

u h (x) = NT(x) u (2.13)

˜ ˜
We will also need the derivative of u h (x); since the coefficients in u are all constants, this derivative is
given by ˜

du h dN T
= ˜ u (2.14)
dx dx ˜

2.2.4 Galerkin approximation

The next step towards obtaining an approximate solution of our problem consists in replacing the solu-
tion u(x) in the weak form (2.10) by the approximate solution u h (x) as defined by (2.13). Since only
the derivative of u(x) appears in (2.10), this effectively means inserting expression (2.14) into (2.10),
yielding

ZL ZL
dφ dN T
c(x) ˜ dx u = φ(x) f (x) dx + φ(0) q0 + φ(L) q L (2.15)
dx dx ˜
0 0

Note that the column matrix u can be taken out of the integral on the left-hand side because it contains
only constants. ˜

Now remember that for u(x) to be the exact solution of our problem, Equation (2.10) had to be satisfied
for any test function φ(x). After replacing u(x) by u h (x) this requirement clearly cannot be satisfied
anymore. We have limited the possible shapes the approximate solution can take to those satisfying the
format (2.13). The exact solution u(x) generally will not fit into this format and therefore u h (x) cannot
satisfy (2.10) for arbitrary test functions φ(x). Thus, if Equation (2.15) is to determine our approximate
solution, we can require it to hold only for a limited class of test functions φ h (x). Given the fact that we

12
defined the approximate solution to be of the form (2.13), it is natural to require these test functions to
be of the same form:

φ h (x) = NT(x) φ (2.16)

˜ ˜
where the shape function matrix N(x) is identical to that in (2.13) and the column matrix φ contains n
˜
scalar constants φi . This way of limiting ˜
the class of test functions is known as Galerkin’s method or the
Bubnov–Galerkin method.

Taking the derivative of (2.16) yields (cf. (2.14))

dφ h dN T
= ˜ φ (2.17)
dx dx ˜
Using this expression and (2.16), Equation (2.15) can be rewritten for test functions φ(x) = φ h (x) as:

ZL ZL
T dN dN T
φ ˜ c(x) ˜ dx u = φT N(x) f (x) dx + φT N(0) q0 + φT N(L) q L (2.18)
˜ dx dx ˜ ˜ ˜ ˜ ˜ ˜ ˜
0 0

where use has been made of the fact that φ is constant. Since the test functions φ h (x) are fully determined
by the coefficients φ, requiring that (2.15)˜ is satisfied for all φ(x) = φ h (x) is equivalent to requiring that
˜ for all φ. This, in turn, is possible if and only if
Equation (2.18) holds
˜
L L
dN T
Z Z
dN
˜ c(x) ˜ dx u = N(x) f (x) dx + N(0) q0 + N(L) q L (2.19)
dx dx ˜ ˜ ˜ ˜
0 0

For a given boundary value problem and a given set of shape functions N(x), the integrals in (2.19) can
˜ multiplied by the coefficient
in principle be evaluated. On the left-hand side this yields an n × n matrix
column u. On the right-hand side we have a sum of column matrices. Equation (2.19) can thus be written
as ˜

Ku = f+ q (2.20)
¯ ˜ ˜ ˜
with
ZL
dN dN T
K= ˜ c(x) ˜ dx (2.21)
¯ dx dx
0
ZL
f= N(x) f (x) dx (2.22)
˜ ˜
0
q = N(0) q0 + N(L) q L (2.23)
˜ ˜ ˜
The matrix K is called system matrix or, particularly in a solid mechanics context, stiffness matrix. The
¯
sum f + q is usually referred to simply as the right-hand side of the system. The first terms in it, f,
˜
represents˜ the influence of the distributed source term in (2.1), f (x); the other contribution comes from˜
the two boundary terms.

One should realise that the boundary contribution q in (2.20) is not fully known. For a given set of shape
˜ without any difficulty. But in the first term, which
functions the second term in (2.23) can be determined

13
is associated with the essential boundary condition at x = 0, the flux quantity q0 is still unknown. This
means that the linear system (2.20) contains n + 1 unknowns: the n components of u plus q0 . Clearly,
one more equation is needed in addition to the n equations provided by (2.20) to be˜ able to solve for
these n + 1 unknowns. This additional equation can be retrieved from the essential boundary condition
(2.2), which we indeed have not used yet. We will require that the approximate solution u h (x) satisfies
this boundary condition exactly, which in our case means that u h (0) = g0 . Substitution of (2.13) then
provides exactly one more linear equation in terms of u:
˜
NT(0) u = g0 (2.24)
˜ ˜
If this equation is combined with the linear system (2.20) the unknowns u and q0 can in principle be
˜ element shape functions
solved for, so that the approximate solution u h (x) is fully determined. If finite
are used, this combination of (2.20) and (2.24) can be done in an elegant and straightforward way. We
will come back to this in Section 2.4, but will now first turn our attention to the finite element shape
functions.

2.3 Finite element shape functions

The way in which the model boundary value problem has been discretised in the previous section is
largely independent of the definition of the shape functions contained in N(x). The only requirements
which they must satisfy are linear independence and continuity. The finite˜ element method uses a spe-
cial class of shape functions which have certain properties that render their use in combination with
Galerkin’s method particularly attractive. These properties mainly result from the fact that relatively
simple, polynomial functions are used, each defined on only a small subregion of the problem domain.

2.3.1 Polynomial interpolation

Let us for a moment assume that the exact solution u(x) of our model boundary value problem is known
and is given by the dashed curve in Figure 2.1. A straightforward way of approximating this solution is
by a polynomial of degree p. For p = 1 we have the linear function

u h (x) = a0 + a1 x (2.25)

which clearly provides a rather poor approximation of u(x) for any values of the constants a0 and a1
(Figure 2.1). Increasing the polynomial degree to p = 2 by setting

u h (x) = a0 + a1 x + a2 x 2 (2.26)

allows to obtain a reasonable approximation of the exact solution if the coefficients a0 , a1 and a2 are
selected carefully. Further increasing the degree p results in an increasingly accurate approximation of
the exact solution, as is illustrated in Figure 2.1 by polynomials of degree four and eight. The numerical
solution is said to converge to the exact solution.

Instead of characterising the polynomial approximations by the coefficients a0 , a1 , a2 , etc. as in Equa-

tions (2.25) and (2.26), we can also define them by their values in a limited set of points, which we will
call nodes. For instance, for the linear function given by (2.25) specifying its value in two separate nodes
suffices to uniquely define the coefficients a0 and a1 , and therefore the entire function u h (x). If we put

14
 u(x)

p=1

u(x), u h (x)
p=2
p=8
p=4

0 L
x

Figure 2.1: Approximation of the exact solution u(x) by polynomials of degree one and two

these nodes at x1 = 0 and x2 = L and denote the function values in them by u 1 = u h (0) and u 2 = u h (L),
some elementary algebra shows that the coefficients a0 and a1 can be written as
u2 − u1
a0 = u 1 a1 = (2.27)
L
Using these expressions, (2.25) can be rewritten in terms of u 1 and u 2 as

u h (x) = N1 (x) u 1 + N2 (x) u 2 (2.28)

where N1 (x) and N2 (x) have been defined as

x x
N1 (x) = 1 − N2 (x) = (2.29)
L L

Expression (2.28) clearly is of the same form as (2.11). It expresses that function values u h (x) for x
between 0 and L can be determined by linear interpolation between the values u 1 at x1 = 0 and u 2 at
x2 = L. Accordingly, the shape functions N1 (x) and N2 (x) are called interpolation polynomials, in this
case of degree one. A graphical representation is given in Figure 2.2. The linear function u h (x) has been
plotted in this diagram as a solid line. The contributions from the two individual terms in (2.2) are shown
as thin dashed lines. It can easily be seen from the diagram that the sum of these contributions indeed
equals u h (x) for all x. Each of the two contributions can be obtained as the product of the nodal value
u i (i = 1, 2) and the corresponding interpolation function Ni (x). The latter have been plotted in the
diagram as thick dashed lines.

As can be seen from the diagram of Figure 2.2, or indeed directly from Equation (2.29), the interpolation
functions have the property that they equal one in the node which they correspond to and vanish in the
other node, or more briefly: Ni (x j ) = δi j , with δi j the Kronecker delta:
(
1 if i = j
δi j = (2.30)
0 if i 6 = j

This property also holds for interpolation functions of higher orders.

15
 u2
u h (x)

u h (x), Ni (x)
u1

1 N1 (x) N2 (x)

0
0 L
x

Figure 2.2: Linear approximation considered as interpolation between nodal values

An illustration of a higher-order interpolation is given in Figure 2.3, which shows the interpolation poly-
nomials for a quadratic interpolation. Since the degree of the polynomial has been increased by one
compared with the linear case, an additional coefficient must be defined and an additional node must
thus be introduced. It is common practice to define these nodes as equidistant, so that for the quadratic
case we have nodes at x1 = 0, x2 = 12 L and x3 = L. Following the same reasoning as above, the
corresponding, second-degree interpolation polynomials can be determined as

3x 2x 2 4x 4x 2 x 2x 2
N1 (x) = 1 − + 2 N2 (x) = − N3 (x) = − + 2 (2.31)
L L L L L L
These functions have been plotted as dashed curves in the diagram of Figure 2.3. Note that they indeed
are equal to one in the node with which they are associated and vanish in the two other nodes. The
quadratic approximate solution u h (x) can again be obtained by multiplying each of the interpolation
functions Ni (x) by the function value u i in the corresponding node and adding these three contributions:

u h (x) = N1 (x) u 1 + N2 (x) u 2 + N3 (x) u 3 (2.32)

Polynomial interpolations of degree p = 3 and higher can be obtained in exactly the same manner as
described above by defining p+1 equidistant nodes and deriving the interpolation polynomials associated
with them. By assembling the nodal values u i and interpolation polynomials Ni (x) in column matrices u
and N(x) according to (2.12) these interpolations can be written in the form of Equation (2.13). They can˜
then ˜be inserted directly into (2.21)–(2.23) in order to build a set of equations from which the unknowns
u can be solved. Once these nodal values have been determined, the approximate solution is entirely
˜
known. Its accuracy can be improved to an arbitrary degree by increasing the order of the interpolation
functions, and therefore the number of nodes and the number of unknowns in the linear system.

A disadvantage of using polynomial shape functions with Galerkin’s method is that the system matrix K
¯
becomes a full matrix and solving the linear system (2.20) therefore becomes expensive if a high degree
of accuracy is required. Furthermore, it can be shown that the linear system becomes ill-conditioned
very rapidly as the order of interpolation is increased. This is caused by the fact that the equidistant
interpolation used here results in non-orthogonal interpolation functions; it can be eliminated by using

16
 u3

u1
u h (x)

u h (x), Ni (x)
u2

N2 (x)
1 N1 (x) N3 (x)

0
1
0 2
L L
x
Figure 2.3: Quadratic approximation and the associated interpolation polynomials of degree two

polynomials which are orthogonal by construction. But perhaps the most serious disadvantage of using a
polynomial interpolation exhibits itself mainly in two and three dimensions, where it renders the method
rather inflexible from a geometrical point of view. In the one-dimensional case which we are considering
here, this is reflected in the fact that the order of the approximate solution is necessarily the same every-
where on the domain (0, L), even if the solution may exhibit strong variations in a certain region of this
domain (near the end right in Figure 2.1) and vary in a much smoother way in other parts (the left part in
the figure).

2.3.2 Finite element interpolation

Finite element shape functions largely avoid the disadvantages of the polynomial shape functions by
using polynomials of a relatively low degree (usually one or two), but only on a small subdomain. By
dividing the entire domain (0, L) into a number of such subdomains, which we call finite elements or
simply elements, a solution can be constructed on the entire domain. The accuracy of approximate
solutions obtained in this way can be improved by increasing the degree of the polynomials used within
the elements, but also by decreasing the size of the elements – at the same time increasing their number.
Indeed, the element size by no means needs to be uniform and can therefore be adapted to local features
of the solution. In two and three dimensions, element sizes and element shapes can be adapted to the
shape of the problem domain, which renders the finite element method particularly flexible in terms of
problem geometries.

An example of a finite element discretisation of our model problem is shown in Figure 2.4. At the
bottom of the figure the division of the domain (0, L) into m = 5 finite elements has been indicated.
The elements do not overlap, but the entire collection of elements covers the complete solution domain.
In order to emphasise that elements may have different sizes, smaller than average elements have been
used near the right end of the domain, where the solution varies more rapidly than elsewhere. Within
each element, linear interpolation functions are used. This means that each element requires two nodes.
However, these nodes are shared between adjacent elements, so that we have a total of m + 1 = 6 nodes.
This assembly of nodes and elements is usually referred to as the finite element mesh or finite element
grid.

17
 u(x), u h (x), Ni (x)
u(x)
u h (x)

Ni (x) Nn (x)

e e+1
0
x1 = 0 xi−1 xi xi+1 xn−1 xn = L
x

Figure 2.4: Finite element discretisation of the model problem using linear interpolation func-
tions

Within the element indicated by e in Figure 2.4 the approximate solution is given by linear interpolation
between nodes i − 1 and i. The interpolation functions associated with these nodes can easily be derived
as (cf. Equation (2.29))
xi − x x − xi−1
Ni−1 (x) = Ni (x) = (2.33)
xi − xi−1 xi − xi−1
for xi−1 ≤ x ≤ xi . However, if we concentrate ourselves for a moment on node i, the function value u i
in this node is also relevant for the interpolation within element e + 1, which is given by the functions
xi+1 − x x − xi
Ni (x) = Ni+1 (x) = (2.34)
xi+1 − xi xi+1 − xi
for xi ≤ x ≤ xi+1 . Node i is not connected to any of the other elements and therefore does not contribute
to them. We can express this by defining Ni (x) = 0 for x < xi−1 and for x > xi+1 . A shape function
Ni (x) associated with node i can thus been defined on the entire domain (0, L) as follows:


 0 if 0 ≤ x ≤ xi−1



x − xi−1


if xi−1 < x ≤ xi


−

Ni (x) = x i x i−1 (2.35)
 xi+1 − x
if xi < x ≤ xi+1


xi+1 − xi






0 if xi+1 < x ≤ L


This function has been plotted in Figure 2.4 as a thick dashed curve. The figure shows once more that
Ni (x) vanishes everywhere except in those elements of which node i forms a part (in this case e and
e + 1). In mathematical terms: its support consists of only the two elements which are connected to it.
Within these elements Ni (x) decreases linearly from unity at x = xi to zero at the adjacent nodes.

Finite element shape functions for the other nodes can be constructed in exactly the same manner as de-
scribed above. Note that the nodes at both ends of the domain are connected to only one element and their
shape functions therefore have a support which consists of only one element. One such shape function,

18
at the right end of the domain, has been plotted in Figure 2.4 as a thin dashed curve. An approximate
solution can now be constructed by multiplying each of these shape functions by the respective degree
of freedom u i and summation according to (2.11) or, equivalently, (2.13). The resulting function u h (x),
plotted in Figure 2.4 as a solid curve, is piecewise linear and is fully determined by the nodal values in
u. Indeed, the value of the function u h (x) in a particular node i always equals the degree of freedom u i
˜ that node because the shape function attached to that node, N (x), is always equal to one and all other
in i
shape functions vanish in it. It is obvious from the figure that a better approximation may be obtained
if the finite element mesh is refined, i.e. if the domain is divided into more elements. Indeed, it can be
proved that in the limit m → ∞ the approximate solution converges to the exact solution.

A better approximation can also be obtained by increasing the order of interpolation within the element.
A quadratic finite element interpolation can be constructed largely following the lines laid out above
for the linear case. The domain (0, L) is again divided into a number of elements, but now quadratic
interpolation polynomials are used within each element. This means that each element has three nodes:
one at each end, which can be shared with adjacent elements, and one in the centre, see Figure 2.5. We
need to distinguish between these two types of nodes when constructing the shape functions associated
with them. For a node i which is shared by two elements the shape function can be derived as


0 if 0 ≤ x ≤ xi−2



(x − xi−2 )(x − xi−1 )


if xi−2 < x ≤ xi


Ni (x) = (xi − xi−2 )(xi − xi−1 )

(2.36)
 (xi+1 − x)(xi+2 − x)
if xi < x ≤ xi+2



(x i+1 − x i )(x i+2 − x i )





0 if xi+2 < x ≤ L


If the node is situated at one of both ends of the domain, the irrelevant part of this expression again must
be dropped. For a node j at the centre of an element (Figure 2.5) we have:



 0 if 0 ≤ x ≤ x j −1


(x − x j −1 )(x j +1 − x)

N j (x) = if x j −1 < x ≤ x j +1 (2.37)
(x j − x j −1 )(x j +1 − x j )





0 if x j +1 < x ≤ L

Instances of both functions have been plotted as dashed curves in Figure 2.5. Again, the support of these
shape functions consists of only two elements, or even only one for the nodes at the centre of an element.
Using these shape functions, a continuous approximate solution can again be constructed on the basis of
a set of discrete, nodal values. An example is given in the diagram as a solid curve. As in the linear case,
this curve passes through the given values in the nodes. Between these nodes, however, it now varies
quadratically as a results of the use of quadratic polynomials within each element.

The linear as well as the quadratic shape functions which we have derived in this section are continuous,
but only piecewise continuously differentiable: within an element they can be differentiated without any
difficulty, but these derivatives exhibit jumps across the element boundaries. These continuity properties
are exactly those allowed by the weak form of our boundary value problem, see Section 2.2.2. We can
therefore use the finite element interpolations as shape functions in the discretisation of the boundary
value problem as described in Section 2.2. Indeed, finite element shape functions prove to be a very
attractive choice for a number of reasons. We already mentioned their geometric flexibility, which be-
comes even more important when we consider boundary value problems in two and three dimensions.

19
 u(x), u h (x), Ni (x)
u h (x)
u(x)

N j (x)
1
Ni (x)

0
x1 = 0 xi xj xn = L
x

Figure 2.5: Finite element discretisation of the model problem using quadratic interpolation
functions

In one dimension, the size of elements can be adapted to local features of the solution, so that a high
accuracy can be obtained where needed. Furthermore, their compact support results in system matrices
which contain a relatively large number of zeros. This property greatly facilitates the efficient solving of
the final system of linear equations. We will come back to this in the next section. The small support of
finite element shape functions also implies that shape functions associated with different nodes are rela-
tively close to orthogonal, which has a positive influence on the conditioning of the final set of equations.
Finally, their simple structure allows to largely automate the construction of a linear system of equations
for a given boundary value problem along the lines of Section 2.2. We will have a closer look at this
aspect in the next section.

2.4 Assembly and solving of the linear system

In the previous two sections, the two crucial ingredients of the finite element method have been intro-
duced: the discretisation of a boundary value problem using Galerkin’s method and a set of interpolation
functions which are particularly suitable for this purpose. This section is devoted to the integration of
these two ingredients. We will show how some of the special properties of finite element interpolation
functions can be used in Galerkin’s method to efficiently construct a linear system of equations and solve
this system.

2.4.1 Assembly of the system matrices

A key property of finite element shape functions is that they have a simple form when considered on a
single element. We can exploit this property by splitting the integration in the system matrix (2.21) and
the source term (2.22) in the right-hand side on an element-by-element basis. For the system matrix K
¯
and the linear discretisation of Figure 2.4 this implies rewriting (2.21) as
x
m Ze+1
X dN dN T
K= ˜ c(x) ˜ dx (2.38)
¯ e=1
dx dx
xe

20
If we concentrate ourselves on the contribution of the first element for a moment, the shape function
derivatives in this element can be obtained by differentiation of (2.35) as
   
dN1 1
 −
 dx   x2 − x1 
 
 dN2   1
   

   
 dx   x2 − x1 
   
 dN3   
dN   dx  
  0 
˜ = = (2.39)


dx  dN4   

 dx  
  0 

   
 dN5   

 dx  
  0 

   
 dN6   
0
dx
The zero entries at the last four positions of this column matrix are a result of the fact that the shape func-
tions of nodes not connected to the first element vanish identically on this element. Inserting expression
(2.39) into the contribution of the first element in (2.38) shows that this contribution has the following
structure:
 
k 1 −k 1 0 0 0 0
 1
−k k 1 0 0 0 0

Zx2  
dN dN T  0 0 0 0 0 0 
˜ c(x) ˜ dx =  (2.40)
 

dx dx  0 0 0 0 0 0 
x1  
 0 0 0 0 0 0
 

0 0 0 0 0 0

where
Zx2
1 1
k = c(x) dx (2.41)
(x2 − x1 )2
x1

In order to avoid storing large numbers of zeros – particularly for larger systems – it is convenient to
store the nonzero entries of (2.40) in a 2 × 2 element matrix K1 :
¯
" #
1 1
k −k
K1 = (2.42)
¯ −k 1 k 1

More generally, we can write for element e:

Zxe+1
dN e dN eT
Ke = ˜ c(x) ˜ dx (2.43)
¯ dx dx
xe

where Ne (x) contains only those shape functions which are relevant for this particular element, i.e.,
˜
those shape functions which are associated with nodes that form part of it. In the present case of one-
dimensional, linear interpolation polynomials each of these element matrices is given simply by
" #
e e
k −k
Ke = (2.44)
¯ −k e k e

21
with
Zxe+1
1
ke = c(x) dx (2.45)
(xe+1 − xe )2
xe

For higher-order interpolations and in more dimensions – as considered in the next chapters – these
expressions may become more complicated and the element matrices are larger than 2 × 2, but they can
still be evaluated relatively easily on an individual element.

Constructing the overall system matrix K now consists in adding the contributions from the different
¯
elements. Since these contributions are generally computed in the condensed form given by (2.43), they
must first be expanded to an n × n matrix by padding them with zeros in the appropriate places before
they can be added. This process of expanding and adding the element matrices is called assembly and is
sometimes denoted as
m
K=
¯ A K¯
e=1
e
(2.46)

For the example of Figure 2.4 it results in a system matrix which reads in terms of the element matrix
components k e :
 
k1 −k 1 0 0 0 0
 1 1
−k k + k 2 −k 2

0 0 0 
 
2 2 3 3
 0 −k k + k −k 0 0 
K= (2.47)
 
3 3 4 4

¯  0 0 −k k + k −k 0 
 
 0 0 0 −k 4 k 4 + k 5 −k 5 
 

0 0 0 0 −k 5 k5

Clearly, the components of the element matrix Ke end up in rows and columns e and e + 1 in this case
¯
because these rows and columns correspond to the nodes of the element. Tracing which components
should be inserted at which position in the global system matrix requires some simple book-keeping,
which lends itself well for automation. Since the element matrices Ke can be computed according to a
¯
relatively simple recipe as well, constructing the global system matrix K can be programmed relatively
¯
easily and executed in a very efficient way on a computer.

With respect to the right-hand side of the linear system (2.20), a very similar assembly procedure can be
set up for the column vector f. The non-vanishing components of the contribution of element e can be
computed as ˜

Zxe+1
e
f = Ne (x) f (x) dx (2.48)
˜ ˜
xe

Using these element contributions, the global column f can then be assembled as
˜
m
f=
˜ A f˜
e=1
e
(2.49)

The assembly operator here stores the entries of the element columns fe in the appropriate rows of the
global column f. ˜
˜
22
The other part of the right-hand side finally, i.e. the column matrix q, can be simplified using another
˜
property of the finite element interpolation functions. We saw in Section 2.3 that the shape function
associated with a particular node i equals one in that node and vanishes in all other nodes: Ni (x j ) = δi j .
Since x = 0 and x = L coincide with nodes 1 and n respectively we have for q according to (2.23):
      ˜
1 0 q0
     
0 0 0
 ..   ..  .
     
q =  .  q0 +  .  q L =  ..  (2.50)
˜       
0 0 0
   

0 1 qL
i.e. the only non-vanishing entries in q are simply the boundary fluxes q0 and q L appearing at the po-
sitions corresponding to the nodes they˜ must be applied at. Constructing q, as well as f, again requires
some elementary book-keeping and simple computations which lend themselves˜ well to˜ being executed
automatically by a computer.

2.4.2 Partitioning of the system

Now that we have established a systematic and relatively straightforward way of setting up the linear set
of equations (2.20) as a discrete replacement of the differential equation (2.1), let us have a look at how
to obtain the discrete solution u from it. As we saw at the end of Section 2.2, one more equation, given
by (2.24) must be added to the ˜system in order to be able to uniquely determine the n + 1 unknowns in it.
This additional equation states that the approximate solution u h (x) must satisfy the essential boundary
condition at x = 0. Again using the property that Ni (x j ) = δi j , it can be rewritten as
h i
1 0 0 . . . 0 u = u 1 = g0 (2.51)
˜
i.e., the first component of u must equal the value g0 prescribed by the essential boundary condition. This
result is not very surprising,˜ since we saw also in Section 2.3 that the approximate solution u h (x) passes
through the nodal values contained in u. Requiring this approximate solution to be equal to g0 at node
number one therefore boils down to setting˜ u =g .
1 0

Inserting the above result for u 1 into Equation (2.20) results in a system of n equations in terms of the n
remaining unknowns u 2 , . . . , u n and q0 :
      
K 11 K 12 . . . K 1n g0 f1 q0
      
 K 21 K 22 . . . K 2n   u 2   f 2   0 
 . = . + .  (2.52)
      
 . .. .. .. 
 .. . . .   .
.   .
.   .
. 
      
K n1 K n2 . . . K nn un fn qL
The lines which have been drawn in this matrix system indicate the separation between known quantities
and unknowns in u and q. In u, we have the known value g0 above the line, whereas the n −1 components
below the line are˜ yet unknown.
˜ ˜ The situation is exactly opposite in q, in which q – above the line – is
0
˜
the only unknown. The horizontal and vertical lines in the matrix K indicate which entries correspond to
¯
known and unknown u i and qi .

More generally, we can partition the linear system (2.20) as follows:

" #" # " # " #
Kgg Kgf ug fg q
¯ ¯ ˜ = ˜ + ˜g (2.53)
Kfg Kff uf ff qf
¯ ¯ ˜ ˜ ˜
23
where ug contains those coefficients u i which are given by essential boundary conditions and uf the
˜ degrees of freedom which are still free. On the other hand, q contains the unknown boundary
remaining ˜
g
fluxes associated with the prescribed values of ug whereas qf is fully˜ known. Note that rewriting the
˜
original system in this form may require some reordering of ˜rows and columns of K. Also note that the
¯
matrices Kgg and Kff will always be square, because an unknown qi is always associated with a prescribed
¯ ¯
u i and vice versa.

We can now use the partitioning of the system to split it into two separate matrix equations:

Kgg ug + Kgf uf = fg + qg (2.54)

¯ ˜ ¯ ˜ ˜ ˜
Kfg ug + Kff uf = ff + qf (2.55)
¯ ˜ ¯ ˜ ˜ ˜
In the second of these equations the only unknowns are the components of uf . This means that we can
rewrite this equation in the form of a standard linear system ˜

Kff uf = ff + qf − Kfg ug (2.56)

¯ ˜ ˜ ˜ ¯ ˜
in which the right-hand side is fully known and the matrix Kff is square. This system can thus be solved
¯
for the unknown coefficients uf . Combining this solution with the already known ug according to
˜ ˜
" #
u
u = ˜g (2.57)
˜ uf
˜
the entire set of degrees of freedom is known and the approximate solution u h (x) thus is fully determined.
If necessary, the unknown boundary fluxes qg can now be retrieved from the top part of the partitioned
system (Equation (2.54)) according to: ˜

qg = Kgg ug + Kgf uf − fg (2.58)

˜ ¯ ˜ ¯ ˜ ˜
In solid mechanics problems, for instance, these quantities represent reaction forces, whereas in temper-
ature problems they give the heat flux which is needed to maintain a prescribed boundary temperature.

2.4.3 Solving the condensed system

We have now reduced the problem of finding an approximate solution for the linear boundary value prob-
lem given by Equations (2.1)–(2.3) to that of finding the solution of the linear system of equations (2.56).
This linear system forms a standard problem in linear algebra and numerical algebra offers a wealth
of algorithms to solve it, the most standard of which is an LU-decomposition (sometimes also called
LU-factorisation) followed by back-substitution. This and other algorithms have been implemented in
standard software libraries and can be used more or less as a black-box solver. A detailed discussion of
the methods used by these solvers and their respective advantages and disadvantages is beyond the scope
of this course. Reference is made for this purpose to courses and textbooks on numerical linear alge-
bra. Here we will briefly discuss the specific properties of linear systems resulting from a finite element
discretisation and indicate which classes of solvers are particularly suitable for this kind of problems.

The matrix Kff which appears in (2.56) is obtained by partitioning of the bigger system matrix (2.21).
¯
For the discretisation by five linear elements considered above, for instance, this implies that the first
row and column can be deleted from expression (2.47) in order to obtain Kff . A number of properties of
¯
the resulting matrix can be noticed immediately. First of all, a considerable part of its entries are equal

24
to zero. Indeed, the density of nonzero entries becomes even much lower if we increase the number of
elements to practical values; we say that the matrix is sparse. As we have seen at the start of this section,
this sparsity is a direct consequence of the small support of finite element shape functions. Sparse solvers
take into account the sparsity of the matrix by acting solely on the nonzero entries when computing the
LU decomposition. This avoids large numbers of trivial operations (e.g. multiplications by zero) and
therefore speeds up the computation.

A further property of (2.47) – and therefore also of the resulting Kff – is that the nonzero entries are
¯
all concentrated in a band around the main diagonal of the matrix. The matrix is said to have a band
structure. Band solvers use this structure to compute the LU decomposition more efficiently than for
a matrix without a particular structure. Note that the band structure of (2.47) results from the order in
which we have stored the degrees of freedom in the global system. A different numbering of the nodes
(and thus of the degrees of freedom) would have resulted in a higher bandwidth, which in turn renders the
band solver less efficient. For this reason, algorithms have been developed which renumber the degrees
of freedom such that an optimal bandwidth is obtained. Modern sparse solvers automatically perform a
similar operation to obtain an optimal sparsity pattern. As a result of these reordering techniques, sparse
solvers have become very competitive with band solvers and seem to be gradually replacing them in
practice.

A third obvious property of (2.47) is that it is symmetric. This symmetry holds for many system matrices
obtained by finite element discretisation, but not for all. Whether a symmetric matrix is obtained depends
on the differential equation which is being solved and sometimes on the precise way in which it is
discretised. That the matrix must be symmetric for the model problem considered in this chapter is
immediately clear from (2.21). For such symmetric matrices the LU decomposition can be replaced by
the Cholesky decomposition, which can be computed more efficiently, particularly if use is also made of
the sparsity or band structure of the matrix.

In the above, we have silently assumed that a direct solver is used. Direct solvers provide a solution
which in principle can be exact – although its accuracy is limited in practice by the numerical precision
of computers. In contrast to this, iterative solvers start from an estimate of the solution and successively
improve this estimate in a number of iterations. The iteration process is stopped when a user-defined
accuracy is met. This accuracy is usually much less than that provided by direct solvers, but this is not
necessarily a problem given the fact that we are computing an approximate solution of the boundary
value problem anyway. Iterative solvers are particularly attractive for very large sparse systems because
they do not require to explicitly assemble – and store – the entire system matrix. Furthermore, they
lend themselves well to parallelisation, so that the solving process can be distributed among a number
of processors. The rate at which the iterative process converges in such methods is highly dependent on
the condition number of the matrix. Preconditioners are therefore used to improve this condition before
entering the iterative solver.

2.5 Isoparametric elements

One aspect of the discretisation to a linear system which we have not paid much attention to is the eval-
uation of the integrals which determine the system matrix and the right-hand side of this linear system.
We have seen in the previous section that the integrals in (2.21) and (2.22) can be conveniently split
into integrals on each of the elements, which can then be assembled to a global system. Most often, the
functions c(x) and f (x) will be constant or of another simple form. In such cases it may be possible
to determine analytical expressions for the integrals in (2.43) and (2.48) and evaluate these expressions

25
for each element. For slightly more complicated functions however, and for higher-order shape func-
tions, this may be difficult or even impossible – and the problem becomes even more compelling in two
and three dimensions. The so-called isoparametric transformation, combined with numerical integration
rules, provides an elegant and efficient method to deal with this difficulty. Because of its generality and
efficiency, this method is used almost exclusively in finite element codes.

2.5.1 Master element and normalised coordinates

We start our development of isoparametric elements by defining a master element which has its own, new
coordinate ξ . In terms of this coordinate, the element is defined on the interval (−1, 1). If we for the
moment consider a linear element, this means that we have nodes at ξ = ±1, see Figure 2.6. The shape
functions which we can use to construct an approximate solution can be written down immediately for
−1 ≤ ξ ≤ 1 as (cf. (2.35))
" #
1
e (1 − ξ )
N (ξ ) = 21 (2.59)
˜ 2
(1 + ξ )

It can easily be verified that these functions satisfy the usual condition Ni (ξ j ) = δi j (see also the figure).
A quadratic element can be constructed simply by adding a node at ξ = 0 and redefining the shape
functions as
 
− 12 (1 − ξ )ξ
 
Ne (ξ ) =  (1 − ξ )(1 + ξ ) (2.60)
˜ 
1

2
ξ (1 + ξ )

1 N1 (ξ ) N2 (ξ )
Ni (ξ )

0
−1 0 1
ξ

1 2

Figure 2.6: Linear master element and corresponding shape functions

The purpose of defining a master element is that we can now consider each element in the physical mesh
as being obtained by a mapping of this master element. Operations on the physical elements, such as
the evaluation of element matrices, can thus be rewritten in terms of the normalised coordinate ξ and
basically become identical for all elements. This is illustrated in Figure 2.7 for the example mesh with

26
two-node elements which we have used before. Within an element e, we can construct the approximate
solution value u h (ξ ) in a point indicated by the normalised coordinate ξ as

u h (ξ ) = NeT(ξ ) ue (2.61)
˜ ˜
where ue = [u e u e+1 ]T contains the degrees of freedom associated with this particular element.
˜

e
x1 xe xe+1 xn
x

−1 0 1
ξ

Figure 2.7: Mapping of master element to elements in mesh

2.5.2 Isoparametric transformation

The above, loose indication of a mapping between the physical coordinate x and the normalised coor-
dinate ξ does not suffice to rewrite the element matrices Ke and fe in terms of ξ . We need a properly
¯ ˜ maps ξ to x and vice versa. Such
defined coordinate transformation for this purpose, which uniquely
a unique relation is provided by the isoparametric transformation. This transformation uses the same
shape functions Ne (ξ ) which we defined above for the approximate solution to define the mapping of ξ
˜
onto x (cf. (2.61)):

x(ξ ) = NeT(ξ ) xe (2.62)

˜ ˜
where xe stands for the physical coordinates of the nodes related to this particular element. For instance,
˜
for element e in the linear mesh of Figure 2.7 we have xe = [xe xe+1 ]T . The prefix ‘iso’ in ‘isoparamet-
˜
ric transformation’ indicates that the coordinate transformation uses the same order of interpolation as
the approximate solution. Some particular applications, which we will discuss no further, require sub-
parametric or superparametric transformations, which respectively use polynomials of lower and higher
order than those of the approximate solution.

The isoparametric transformation given by (2.62) allows us to rewrite the integral which defines Ke (see
¯
Equation (2.43)) as

Z1
e dN e dN eT
K = ˜ c(x) ˜ J (ξ ) dξ (2.63)
¯ dx dx
−1

where the Jacobian

dx dN eT
J (ξ ) = = ˜ xe (2.64)
dξ dξ ˜

27
follows by differentiation of (2.62) and can be evaluated in a straightforward fashion for any ξ . For the
linear shape functions (2.59) J (ξ ) is constant, but for higher-order elements it may vary with ξ . The
shape function derivatives in (2.63) can be determined for given ξ via
dN e dN e dx dN e
˜ = ˜ = ˜ J (ξ ) (2.65)
dξ dx dξ dx
from which dividing by J (ξ ) gives
dN e dN e −1
˜ = ˜ J (ξ ) (2.66)
dx dξ

c(x) in (2.63) can be evaluated by directly inserting the transformation rule (2.62): c(x) =
The factor
c x(ξ ) . In computer codes this function is sometimes replaced by an interpolated function ch (ξ ) ac-
cording to

ch (ξ ) = NeT(ξ ) ce (2.67)
˜ ˜
where the nodal values ce are to be specified by the user. Depending on the exact form of c(x), this may
involve some degree of˜approximation. But this approximation will generally be of the same order as
that of the finite element interpolation of the solution and is therefore considered quite acceptable.

The transformation of the column matrix fe follows along the same lines as that of Ke . The result reads
˜ ¯
Z1
fe = Ne (ξ ) f (x)J (ξ ) dξ (2.68)
˜ ˜
−1


where f (x) can be evaluated as f (x) = f x(ξ ) . Alternatively, like the function c(x) above, f (x) can
be interpolated between its values in the nodes.

2.5.3 Numerical integration

If the exact functions c(x) and f (x) are used in the transformed integrals in (2.63) and (2.68), these
integrals cannot generally be evaluated analytically. But even when c(x) and f (x) are replaced by finite
element interpolations, and an exact evaluation may in principle be possible, these integrals are usually
approximated by a numerical integration rule – also referred to as a quadrature rule. These rules can
be constructed such that they give an exact result for some cases that could in principle be evaluated
analytically and an approximate value for all other cases. As long as it is sufficiently accurate, this
approximation is quite acceptable since we are constructing an approximate solution anyway.

Quadrature rules replace the integration in for instance (2.63) by a weighted sum of the integrand values
in a number of points, the so-called integration points. Perhaps the most simple approximation is that by
the one-point quadrature rule

Z1
g(ξ ) dξ ≈ 2g(0) (2.69)
−1

A graphical representation of this approximation is given in Figure 2.8. From the figure it can be im-
mediately concluded that this midpoint rule is exact for linear g(ξ ), i.e. for a polynomial of degree one.

28
 g(ξ ) −1 0 1
ξ

Figure 2.8: Graphical interpretation of the midpoint rule

More generally, we can write

Z1 s
X
g(ξ ) dξ ≈ wk g(ξk ) (2.70)
−1 k=1

where wk are the weight factors associated with the integration points ξk . It can be shown that for a given
number of integration points, their optimal positions are in the so-called Gauss points. These points are
optimal in the sense that together with the corresponding weight factors they provide an exact integration
for polynomials of the highest possible degree for a given number of points. For this reason Gauss
quadrature has become by far the most popular in finite element applications. Table 2.1 lists the Gauss
point positions and weight factors for up to three-point integration. Also listed is the polynomial degree
p which is still integrated exactly by the respective rules. Clearly, this degree and thus the accuracy
increase as the number of Gauss points increases.

s p ξk wk
1 1 0 2
√
2 3 −1/ 3 1
√
1/ 3 1
√
3 5 − 3/5 5/9
0 8/9
√
3/5 5/9

Table 2.1: Gauss point coordinates, weight factors and accuracy

Using the Gauss integration rule, the element matrix Ke can be written as
¯
s e eT
X dN
dN
Ke = wk ˜ c x(ξk ) ˜ J (ξk ) (2.71)
¯ k=1
dx ξ =ξk dx ξ =ξk

while for the contribution fe to the right-hand side we have

˜
X s
fe = wk Ne (ξk ) f x(ξk ) J (ξk )


(2.72)
˜ k=1
˜

29
 g(ξ ) −1 − √1 0 √1 1
3 3
ξ

Figure 2.9: Two-point Gauss quadrature rule

The number of Gauss point required for an exact integration of these matrices depends on the functions
c(x) and f (x). If these functions are not polynomials, an exact integration may not be possible using
Gauss quadrature – which, as indicated above, is not really a big problem given the approximate character
these integrals already have. But if we assume that c(x) are polynomials of at the most the same degree as
the interpolation polynomials, a two-point integration suffices for linear elements, whereas a three-point
integration would be appropriate for quadratic elements. For constant c(x) and f (x), these numbers can
be reduced by one.

30
 3

Two and three-dimensional scalar

boundary value problems
3.1 Model boundary value problem

In the previous chapter we have demonstrated the main concepts of the finite element method for perhaps
the simplest class of problems it can be applied to: a single, second-order, linear ordinary differential
equation accompanied by two boundary conditions. However, the main advantages of the finite element
method over other numerical solution techniques become much more pronounced for more realistic prob-
lems from our three-dimensional world. We will take the step towards two and three dimensions in this
chapter. For the time being we limit ourselves to equations which are formulated in terms of a scalar-
valued function – vector problems are the subject of the next chapter. Like in the previous chapter,
we will use a model problem to explain the steps that must be taken to obtain a finite element solution
for boundary value problems in more than one dimension. However, the techniques developed for this
specific problem are by no means limited to it, but can be applied to a wider class of equations.

The model equation which we consider in this chapter reads in its strong form
 
∇E · C(E E + f (E
x )· ∇u x) = 0 (3.1)

Compared with the one-dimensional equation (2.1), the coordinate x has been replaced by a coordinate
vector or position vector xE, which can be defined for instance in a three-dimensional Cartesian basis
{Ee1 , eE2 , eE3 } as

xE = x1 eE1 + x2 eE2 + x3 eE3 (3.2)

In two-dimensional problems the last term of the above expression can be dropped. The spatial derivative
d/dx in (2.1) has become a gradient in (3.1). In the same, three-dimensional Cartesian basis as above
the gradient operator reads

E = eE1 ∂ + eE2 ∂ + eE3 ∂

∇ (3.3)
∂ x1 ∂ x2 ∂ x3
where ∂/∂ xl denotes partial differentiation with respect to xl and the last term can again be dropped in two
dimensions. Clearly, (3.1) has become a partial differential equation instead of the ordinary differential
equation of Chapter 2.

As a consequence of the fact that directions are now associated with the spatial derivatives in (3.1), the
coefficient C(E x ) has become a second-order tensor. If this tensor is assumed to be isotropic, i.e. if it is of
the form C(E x ) = c(E x )I with I the second-order identity tensor, a slightly simpler equation is obtained:
 
∇E · c(E E + f (E
x )∇u x) = 0 (3.4)

This equation for instance governs heat conduction according to Fourier’s law. If furthermore c(x) = 1
for all xE, (3.4) reduces to the Poisson equation or, if also f (E
x ) = 0, the Laplace equation. In the
remainder of this chapter we will use the general form of Equation (3.1). The corresponding results for
the reduced equations based on (3.4) can be obtained completely analogically, or simply by substitution
of C(Ex ) = c(Ex )I.

Equation (3.1) is normally valid only on a bounded domain. In the three-dimensional case this domain is
a volume in space, whereas in two dimensions it consists of a certain area. Here and in the remainder of
these notes we will use the terminology associated with the three-dimensional situation and thus refer to
the domain as a volume V which is bounded by a surface S. The equations which are derived using this
notation, however, are equally applicable to two-dimensional problems if V is interpreted as a surface
area and S as the curve which bounds this area, see Figure 3.1 for an example. Indeed, most of the

33
 Sq

Sg

Figure 3.1: Problem domain V together with the boundary conditions applied to it

illustrations in this chapter will represent two-dimensional problems because these can be visualised
more clearly.

In order to complete the boundary value problem, boundary conditions must be specified on the surface S,
which describe the interaction of the phenomena represented by the partial differential equation with the
surroundings of V . For the type of equation we are considering here, prescribing one boundary condition
in each point of S can be shown to suffice. We will again consider two types of boundary conditions:
those of the essential (or Dirichlet) type and natural (or Neumann) conditions. The essential condition
acts on a certain part Sg ⊂ S of the boundary (Figure 3.1):

x ) = g(E
u(E x) (3.5)

with g(Ex ) a given function on Sg . The remaining part of the boundary, Sq = S \ Sg , has a natural
condition, which is of the form

nE (E
x )·C(E E = q(E
x )· ∇u x) (3.6)

with nE (E
x ) the unit outward normal vector on S and q(Ex ) a given function. In terms of heat conduction,
the first of these conditions represents a fixed temperature, whereas the second represents a given heat
flux across the boundary. In particular, the homogeneous version of (3.6), i.e. q(E x ) = 0, prohibits any
heat flow across the boundary and is therefore often referred to as an insulation condition.

3.2 Finite element discretisation

Like in the one-dimensional case, the discretisation of the boundary value problem is based on a weak
form of the differential equation. Using finite element shape functions, Galerkin’s method yields a lin-
ear system of equations in terms of nodal values of the approximate solution. We will only consider
isoparametric elements in this chapter, which means that the system matrix and right-hand side can be
assembled from element contributions, which in turn can be obtained by numerical integration on a mas-
ter element. Once the linear system has been set up, the approximate solution is obtained by partitioning
of the degrees of freedom and solving of the resulting reduced system.

34
Conceptually, the above steps are largely equivalent with those taken for the one-dimensional problem of
Chapter 2. However, the fact that we are dealing with more than one spatial dimension gives rise to some
technical differences which we will discuss in detail below for each of the steps. Most notable among
these differences is the fact that the line elements of Chapter 2 must now be replaced by elementary areas
or volumes. This opens a wealth of possible shapes and interpolations, which we will detail separately
for the two and three-dimensional cases near the end of this chapter.

3.2.1 Weak form of the partial differential equation

We will first derive the weak form of the partial differential equation (3.1) using the weighted residuals
approach. Like in the one-dimensional case, the weighted residuals equivalent of (3.1) is obtained by
multiplication by a test function φ(E
x ) and integration on the domain V :
Z h   i
φ(E E · C(E
x) ∇ E + φ(E
x )· ∇u x ) dV = 0
x ) f (E (3.7)
V

Following the same reasoning as in Chapter 2 it can easily be understood that requiring this relation to
hold for all possible φ(E
x ) is equivalent to (3.1).

In order to reduce the continuity requirements on u(E x ), the first term of Equation (3.7) is rewritten using
the product rule for the divergence and the divergence theorem as:
Z   Z   Z
φ(E E
x ) ∇ · C(E E E
x )· ∇u dV = ∇ · φ(E E
x )· ∇u dV − ∇φ
x ) C(E E ·C(E E dV
x )· ∇u
V V
Z ZV
= x ) nE (E
φ(E E dS −
x )· ∇u
x )·C(E E ·C(E
∇φ E dV
x )· ∇u (3.8)
S V

Note that these operations are equivalent with the integration by parts in the one-dimensional case.
The resulting boundary contribution, which consisted of two discrete terms in the one-dimensional case
(cf. (2.8)), now takes the form of an integral on S. If we disregard the distinction between the essential
and natural boundary parts of S for a moment and define the boundary flux

x ) = nE (E
q(E x )· C(E E
x )· ∇u (3.9)

everywhere on it, this boundary integral can be rewritten as

Z Z
x ) nE(E
φ(E x )· C(E E
x )· ∇u dS = φ(E x ) dS
x )q(E (3.10)
S S

Using the above relations, (3.7) can finally be reshaped to the weak form
Z Z Z
E
∇φ · C(E E
x )· ∇u dV = φ(E x ) dV + φ(E
x ) f (E x ) dS
x )q(E (3.11)
V V S

3.2.2 Discretisation of the problem domain

Now that we have replaced the strong form of the partial differential equation, (3.1), by the weak form
(3.11), we can start constructing an approximate solution of the boundary value problem. To this end

35
the domain V is first subdivided into a number of finite-sized elements. For isoparametric elements the
shape of these elements is derived from an elementary-shaped master element. In two dimensions square
or triangular master elements are usually used; the corresponding shapes in three dimensions are cubes
and tetrahedra. An example of a square, four-node master element is given in the left part of Figure 3.2.
Associated with it is a normalised coordinate vector ξE , with Cartesian components ξ1 and ξ2 which range
from minus one to one. Each node i of the element has a corresponding shape function Nie (ξE ) which is
defined in terms of ξE . We will detail the shape functions for this and other master elements later in this
chapter. For now it suffices to state that, like in the one-dimensional case, they are defined such that in
node j we have Nie (ξE j ) = δi j (i, j = 1, . . . , 4) and they provide some kind of interpolation – in this case
bilinear – in the interior of the element.

4 3 S
Sh
eE2

Vh
eE1

1 2
V \ Vh

Figure 3.2: Square master element (left) and quadrilateral element discretisation of a problem
domain V derived from it (right)

Unlike the master element, the elements which are derived from it need not be square, but may take on
other quadrilateral shapes which depend on the true position of the nodes associated with them. The
boundaries between these nodes, however, are governed by the isoparametric mapping

xE(ξE ) = NeT(ξE ) xEe (3.12)

˜ ˜
where the vector column xEe contains the position vectors of the element nodes. As a consequence of
˜ boundaries can only take shapes which are allowed by the shape functions.
using this mapping, element
For the bilinear element of Figure 3.2 this means that the element boundaries are straight. For higher-
order elements, some curvature of the element boundaries may be allowed, but only to the polynomial
degree of the element shape functions, e.g. a parabolic shape for quadratic elements.

The above limitation to the shape of elements implies that arbitrarily shaped domains cannot always be
completely covered by elements. An example is given on the right side of Figure 3.2. The straight edges
imposed by the bilinear master element result in the curved part of the boundary not being followed
exactly by the finite element mesh. If the domain occupied by the mesh is denoted by V h , a thin shell
V \ V h (shaded in the figure) remains which cannot be taken into account in the discretisation. Also,
since this shell will not be part of the computational problem, the boundary conditions must be applied at
S h rather than at S. Essential boundary conditions are imposed by fixing the nodal values of those nodes
which are lying on Sg – or close to it where S h does not coincide with S. This results in the element faces
between such nodes – on which values are interpolated between the nodes – to be fixed as well. We will
denote this part of the discretised boundary by Sgh . On the remaining part Sqh = S h \ Sgh the flux quantity
q(Ex ) must be prescribed based on the data for the nearby Sq .

36
Applying the boundary conditions on S h instead of S and replacing the exact domain V in the weak form
(3.11) by V h , as we will do next, clearly introduces some approximation. As the discretisation is refined
(smaller elements), however, this approximation becomes more and more accurate, so that convergence
of the numerical solution to the exact solution can nevertheless be ensured.

3.2.3 Discretisation of the weak form

If we rewrite the weak form (3.11) on the discretised domain V h rather than the exact domain V and at
the same time replace u(E
x ) and φ(E
x ) by the discretised fields

x ) = NT(E
u h (E x) u x ) = NT(E
φ h (E x) φ (3.13)
˜ ˜ ˜ ˜
we have
Z Z Z
φT ∇N E · C(E E T dV u = φT N(E
x )· ∇N x ) dV + φT N(E
x ) f (E x ) dS
x )q(E (3.14)
˜ h ˜ ˜ ˜ ˜ h ˜ ˜ h ˜
V V S

where C(E x ) are understood to be defined on the discretised domain V h and its boundary
x ) and q(E
x ), f (E
S h instead of V and S respectively. The shape functions N(E x ) must be continuous and piecewise differen-
˜
tiable; we will see later that isoparametric elements automatically satisfy these continuity requirements.
Note that for isoparametric elements the shape functions are actually defined
in terms of the local co-
E
ordinates ξ ; the column matrix N(E x ) must therefore be evaluated as N xE(ξ ) where xE(ξE ) represents the
E
isoparametric mapping on the element which point xE belongs to. The ˜same holds for the gradients of the
˜
shape functions contained in ∇N E . Note also that ∇N
E stands for a column of vectors which is defined as
˜ ˜
h iT
E = ∇N
∇N E 1 ∇N E 2 . . . ∇N E n (3.15)
˜
The inner products in (3.14) act between the components of this column matrix and the second-order
tensor C(Ex ). The result of these products is therefore a matrix of scalar functions and the integral on the
left-hand side of the equation thus yields a constant scalar matrix.

Now we must realise again that the weak form (3.11) had to hold for all φ(E x ) and that therefore the
discretised form (3.14) must hold for all possible φ. This can only be true if
˜
Ku = f+ q (3.16)
¯ ˜ ˜ ˜
where
Z
K = ∇N E · C(E E T dV
x )· ∇N (3.17)
¯ ˜ ˜
h
ZV
f = N(E x ) dV
x ) f (E (3.18)
˜ ˜
Vh
Z
q = N(E x )q(E x ) dS (3.19)
˜ ˜
Sh

We thus end up with a linear system of exactly the same form as in the one-dimensional case, but with
a slightly different definition of the matrix K and the column matrices f and q. In particular, the lat-
¯
ter contribution is defined as a boundary integral now, whereas it had a˜ very ˜simple structure in one
dimension.

37
3.2.4 Isoparametric transformation and numerical integration

The integrals in expressions (3.17)–(3.19) are not evaluated directly, but element by element and in terms
of the normalised coordinates associated with each element. The isoparametric mapping which sets the
relation between the local coordinate vector ξE and the global coordinate vector xE was given in (3.12).
Using this transformation, the gradient with respect to xE which appears in the system matrix and right-
E ξ . In a three-dimensional
hand side can be rewritten as a gradient with respect to ξE , which we denote by ∇
Cartesian basis this gradient is defined as

E ξ = eE1 ∂ + eE2 ∂ + eE3 ∂

∇ (3.20)
∂ξ1 ∂ξ2 ∂ξ3

E ξ and ∇
The relation between the operators ∇ E reads

E ξ = J (ξE )· ∇
∇ E (3.21)

where the Jacobian tensor J (ξE ) has been defined as

J (ξE ) = ∇
E ξ xE (3.22)

Inserting the isoparametric mapping (3.12) into this definition yields

E ξ NeT xE e
J (ξE ) = ∇ (3.23)
˜ ˜
which can be computed for any given ξE . Pre-multiplying Equation (3.21) by the inverse of J (ξE ) yields
the transformation of gradients in the opposite direction:
E = J −1(ξE )· ∇
∇ Eξ (3.24)

The above expression can be used to compute the gradients of the shape functions with respect to xE which
appear in the definition for the system matrix (3.17). For the relevant shape functions on element e we
have:
E e = J −1(ξE )· ∇
∇N E ξ Ne (3.25)
˜ ˜
in which ∇ E ξ Ne follows simply by differentiation of the local shape functions Ne (ξE ). The remaining factor
C(E ˜
x ) in the integrand ˜ in terms of ξE by directly
of (3.17), like in the one-dimensional case, can be written
using the mapping (3.12) on each element. The contribution of element e to the integral can now be
transformed to an integral on a cube Q = (−1, 1) × (−1, 1) × (−1, 1) as
Z Z
e E e E eT E e· C(E E eT det J (ξE ) dQ


K = ∇N · C(E x )· ∇N dV = ∇N x )· ∇N (3.26)
¯ ˜ ˜ ˜ ˜
Ve Q


where the factor det J (ξ ) appears as a consequence of the coordinate transformation (cf. the factor J (ξ )
appearing in (2.63)). This can be seen by writing the infinitesimal volume dV as
     
dV = dE x1 × dE x3 = dξE1 · J(ξE ) × dξE2 · J (ξE ) · dξE3 · J (ξE )
x2 · dE

= det J(ξE ) dξE1 × dξE2 · dξE3 = det J(ξE ) dQ

(3.27)

Evaluating the integral in (3.26) analytically is usually cumbersome, if not impossible. This is not only
due to the function C(E
x ), which may in principle be arbitrary, but also to the determinant of the Jacobian,

38
which very easily gets a complicated form for arbitrarily shaped elements. For this reason, numerical
integration is generally used, which yields for the present case (cf. (2.71)):
s
E e E eT
X
Ke = · C xE (ξEk ) · ∇N det J(ξEk )



wk ∇N (3.28)
¯ ˜ ξE =ξEk ˜ ξE =ξEk
k=1

The optimum position of the integration points ξEk and the corresponding weight factors wk depend on the
element type and are therefore discussed separately for each type of element in the following sections.

The right-hand side contribution fe which is related to the source term f (E

x ) can be written analogically
to the above as ˜
Xs
fe = wk Ne(ξEk ) f xE(ξEk ) det J (ξEk )



(3.29)
˜ k=1
˜
The contribution which stems from the influence of the boundary, q, requires a slightly different treatment
˜ implies that elements in the interior
since it involves a surface integral rather than a volume integral. This
of the domain do not contribute and only one layer of elements at the surface of the domain needs to be
considered. Furthermore, on the part of the boundary on which an essential boundary condition is acting,
Sgh , the flux q(E
x ) is yet unknown and no contribution to the column matrix q can thus be determined –
˜
indeed, these entries will be a result of the computation after partitioning and solving the reduced system,
see ahead to Section 3.2.5. What remains is the contribution from the boundary region Sqh , on which q(E x)
is known. Often q(E x ) will be equal to zero, so that the contribution from Sqh to the surface integral in
(3.19) vanishes. This is for instance the case at insulated walls in heat conduction problems.

If a non-zero q(Ex ) is defined on Sqh , the integral is evaluated numerically on an element-by-element basis.
A useful property of isoparametric element families in this respect is that the shape functions of a three-
dimensional element on each of its faces (i.e. for one component of ξE equal to plus or minus one) reduce
to the shape functions of the corresponding two-dimensional element. Similarly, the shape functions of
a two-dimensional element reduce to those of the corresponding one-dimensional element on each of its
edges. This can easily be seen for the bilinear master element of Figure 3.2, for which we stated before
that the shape functions vary linearly along its edges. Using this property, we can rewrite the contribution
of element e to the boundary integral in (3.19) as
Z Z
e e
x ) dS = Ne′(ξE ′ ) q xE(ξE ′ ) det J ′ (ξE ′ ) dQ ′



q = N (E x ) q(E (3.30)
˜ ˜ ˜
Se Q′

where the prime denotes local coordinates and shape functions associated with the lower-dimensional
boundary element. Using the appropriate quadrature rules for this surface element (or line element in
two dimensions), the integral in the above expression can be approximated by
s ′
X
e
wk′ Ne′(ξEk′ ) q xE(ξEk′ ) det J ′ (ξEk′ )



q = (3.31)
˜ k=1
˜

3.2.5 Assembly, partitioning and solving

Computing the element contributions Ke , fe and qe according to Equations (3.28), (3.29) and (3.31) is
¯
straightforward and can therefore easily be˜ done automatically
˜ by a computer. Once this has been done
for all elements, these contributions can be assembled to the global matrices:
m m m
K=
¯ A K¯
e=1
e
f=
˜ A f˜
e=1
e
q=
˜
A q˜
e=1
e
(3.32)

39
This process is in principle identical to that for the one-dimensional case. However, in general two and
three-dimensional cases one can no longer rely on a straightforward numbering of elements and nodes
like in the one-dimensional example of Section 2.4.1. The element contributions therefore generally do
not appear as compact blocks in the global matrices, but may be more dispersed. For instance, the entry
e
K 12 in the element matrix of an element e must be added to the global matrix K in row i, column j ,
¯
where i is the global node number corresponding to the first node of the element and j corresponds to
the second node. Storing all element contributions at the appropriate positions in the global system thus
requires some knowledge on which nodes are connected to which element, i.e. on the mesh connectivity.
This bookkeeping, however, can easily be dealt with by a computer code.

Similarly, the processing of essential boundary conditions by partitioning of the global system can be
done automatically if the numbers of the prescribed nodes are made available to the code. Once the
degrees of freedom have been partitioned into a set of given values ug and the remainder uf , the reduction
of the global system of equations to a solvable set of equations and ˜ equations can
˜ the solving of these
go ahead exactly as described for the one-dimensional case in Sections 2.4.2 and 2.4.3. If desired, the
unknown nodal fluxes qg associated with prescribed degrees of freedom can subsequently be determined
via (2.58). ˜

3.3 Two-dimensional elements

The remainder of this chapter is devoted to a discussion of the different kinds of elements that are the
most popular in two and three-dimensional finite element analyses. In particular we will discuss the
shape functions associated with them and the numerical integration rules which can be used to build the
element matrix and right-hand side. If these ingredients are inserted into the discussion of the previous
section, a complete discretisation strategy is obtained, which can be implemented in a computer code
straightaway.

We begin our discussion with two-dimensional element formulations because these are easier to visu-
alise. The three-dimensional counterparts, which will be discussed in Section 3.4, prove to be relatively
straightforward extensions of these two-dimensional elements.

3.3.1 Bilinear quadrilaterals

We already saw an example of a two-dimensional master element in Figure 3.2. Figure 3.3 shows this
master element again, together with a shape it could take in a physical mesh. The element has four
edges and is therefore called a quadrilateral element or simply quad. It forms a natural extension of
the one-dimensional linear isoparametric element (see Figure 2.6). The single normalised coordinate ξ
of the one-dimensional element is replaced by a vector ξE with Cartesian components ξ1 and ξ2 which
each range from minus one to plus one. A square master element is thus obtained with a node at each
corner. By convention these nodes are numbered counterclockwise, starting at the node with coordinates
ξ1 = ξ2 = −1. Note that these numbers are not the global node numbers of a mapped version of this
element, but simply serve to define the order in which we will store the shape functions associated with
this element.

The shape functions of the element can be obtained by multiplying in all possible combinations the shape

40
 ξ2
4 3 3
1 4

−1 0 1 ξ1
2
−1
1 2
1

Figure 3.3: Bilinear quadrilateral master element (left) and possible physical shape (right)

functions of the one-dimensional element, with ξ replaced by ξ1 and ξ2 respectively:

 
1
4
(1 − ξ1 )(1 − ξ2 )
1
 4 (1 + ξ1 )(1 − ξ2 )

e E
N (ξ ) =  1
  (3.33)
˜  4 (1 + ξ1 )(1 + ξ2 )

1
4
(1 − ξ1 )(1 + ξ2 )
It can easily be verified that in each node only the corresponding shape function equals one and all other
shape functions vanish. In between the nodes, (3.33) provides a bilinear interpolation: each of the shape
functions depends linearly on each of the coordinates ξ1 and ξ2 . A graphical representation of (3.33) is
given in Figure 3.4, which shows a plot of each of the first component of Ne(ξE ) versus the normalised
coordinates ξ1 and ξ2 . ˜

1 4
ξ2
N1e (ξE )

3
0
1 ξ1

Figure 3.4: Shape function N1e(ξE ) associated with the first node of the bilinear element; the shape
functions of the other nodes are obtained by rotation

It can be seen from the figure, but also directly from (3.33), that the shape functions vary linearly along
the boundaries of the master element. Since this would also hold for an adjacent bilinear element which
shares two of its nodes with our reference element, continuity of the shape functions across element
boundaries is ensured automatically – we say that the elements are compatible. As a result of this com-
patibility, continuity of the approximate solution is ensured a priori in isoparametric elements, something
that is not trivial in other classes of elements.

The fact that shape function variations are linear along the element boundaries also implies that after
mapping to the physical space described by the coordinate vector xE the boundaries of the resulting quadri-
lateral element must still be straight, since they are given by linear interpolation between the global node

41
positions xEe (cf. the discussion in Section 3.2.2). Curved boundaries can only be obtained by increasing
˜ interpolation, which we will discuss below.
the order of

The isoparametric transformation also imposes restrictions on the relative positions of the nodes in the
physical mesh. If these positions are such that the shape of the element differs too dramatically from the
square shape of the master element, the Jacobian tensor may become singular for some ξE and its inverse
thus no longer exists. This may happen for instance if the angle between two edges of the element
becomes equal to or larger than 180◦ . As a general rule, one should therefore avoid severely distorted
elements and strive for meshes in which elements are as close to their master-shape as possible.

3.3.2 Higher-order quadrilaterals: Lagrange family

Following the same strategy as for the bilinear element, expanding the one-dimensional quadratic shape
functions given by (2.60) to two dimensions results in a biquadratic element. The resulting master
element has nine nodes: eight on its boundary and one at the centre. They are usually numbered such
that the corner nodes have the same number as in the four-node element, see Figure 3.5. The shape
functions corresponding to these nodes are:

 
1
(1 − ξ1 ) ξ1 (1 − ξ2 ) ξ2
 41 
− ξ1 (1 + ξ1 )(1 − ξ2 ) ξ2

 4 
 1
ξ1 (1 + ξ1 ) ξ2 (1 + ξ2 )

 4
 1 
− (1 − ξ1 ) ξ1 ξ2 (1 + ξ2 )
 4 
Ne(ξE ) =
 1
− 2 (1− ξ1 )(1 + ξ1 )(1 − ξ2 ) ξ2  (3.34)

˜ 
 1

 2 ξ1 (1 + ξ1 )(1 − ξ2 )(1 + ξ2 ) 
 1
− ξ1 )(1 + ξ1 ) ξ2 (1 + ξ2 )

 2 (1
 
− 1 (1− ξ1 ) ξ1 (1 − ξ2 )(1 + ξ2 )
 2 
(1 − ξ1 )(1 + ξ1 )(1 − ξ2 )(1 + ξ2 )

Figure 3.6 shows a graphical representation of some of these polynomials. They vary quadratically
along the boundaries of the master element. As a result, these elements may have parabolic boundaries
in the physical mesh, the amount of curvature being determined by the position of the midside nodes
(Figure 3.5). However, these midside nodes should not be placed to far away from the centre position
between the adjacent corner nodes in order to avoid excessive element distortion and breakdown of the
isoparametric transformation.

Like the bilinear element, the biquadratic isoparametric element is compatible with itself. As a result,
continuity of the approximate solution is ensured for arbitrary meshes of biquadratic elements. Note,
however, that this does not hold for combinations of bilinear and biquadratic elements since their respec-
tive shape functions vary linearly and quadratically along the common boundary and thus are incompat-
ible.

Bicubic master elements and elements of even higher order can be generated completely analogically to
the quadratic element. The class of elements which we may thus obtain is referred to as the Lagrange
family. Orders higher than biquadratic, however, are not very often used in practice.

42
 ξ2
4 7 3 7 3
1 4

9 6
8 9 6
−1 0 1 ξ1 8
2
−1 5
1 5 2
1

Figure 3.5: Biquadratic quadrilateral element of the Lagrange family: master element (left) and
possible physical shape (right)

1 4 1 4 1
7 7

N9e (ξE )
N5e (ξE )
N1e (ξE )

8
9 3 3 3
0 0 9 0
1 6 1 6 1 6
5 5 5
2 2 2

Figure 3.6: Selected shape functions of the biquadratic element of the Lagrange family

3.3.3 Serendipity family

If we multiply the shape functions of the Lagrangian biquadratic element according to (3.34) by the
corresponding nodal degrees of freedom ue and add up all contributions of identical orders of ξ1 and ξ2 ,
the result can be written as ˜

u h (ξE ) = a00 + (a10 ξ1 + a01 ξ2 ) + a20 ξ12 + a11 ξ1 ξ2 + a02 ξ22 + a21 ξ12 ξ2 + a12 ξ1 ξ22 + a22 ξ12 ξ22 (3.35)

where the nine coefficients ai j are linear combinations of the nodal degrees of freedom and brackets
have been used to group terms of the same order. This expression shows that all polynomial terms
of orders up to two are present. The third-order terms, however, are incomplete: the combinations ξ13
and ξ23 are missing. Furthermore, there is only one term present of order four. As a consequence of
this incompleteness of the polynomial, the order of approximation of the element is not the same in all
directions. This is a situation which we generally want to avoid because the incomplete higher-order
contributions do not necessarily contribute to a higher accuracy, but nevertheless add degrees of freedom
to our system. Note that the bilinear quadrilateral element is also incomplete in this sense, because it
provides a ξ1 ξ2 term, but not the ξ12 and ξ22 terms which would accompany it in a complete quadratic
polynomial.

The so-called serendipity element family aims at providing interpolation polynomials which are as com-
plete as possible. The quadratic serendipity element has eight nodes – it lacks the centre node of the

43
Lagrange family. Its shape functions read, in the order indicated in the master element of Figure 3.7,
 1 
− 4 (1 − ξ1 ) (1 − ξ2 )(1 + ξ1 + ξ2 )
 1
− 4 (1 + ξ1 ) (1 − ξ2 )(1 − ξ1 + ξ2 )

 
− 1 (1 + ξ ) (1 + ξ )(1 − ξ − ξ )
 4 1 2 1 2 
 1
− 4 (1 − ξ1 ) (1 + ξ2 )(1 + ξ1 − ξ2 )

e E
N (ξ ) =  1
  (3.36)
˜  2 (1 − ξ12 )(1 − ξ2 )


 1 
 (1 + ξ1 )(1 − ξ 2 ) 
 2 2 
 1
 2 (1 − ξ12 )(1 + ξ2 )


1
2
(1 − ξ1 )(1 − ξ22 )

Figure 3.8 shows the two typical shapes of these functions. It can easily be verified that these shape
functions contain only two polynomial terms in excess of a complete quadratic polynomial, namely ξ12 ξ2
and ξ1 ξ22 . It thus provides the same complete polynomial degree as the quadratic Lagrange element, but
at a lower cost.

ξ2
4 7 3 7 3
1 4

6
8 6
−1 0 1 ξ1 8
2
−1 5
1 5 2
1

Figure 3.7: Quadratic quadrilateral element of the serendipity family

1 4 1 4
7 7
N5e (ξE )
N1e (ξE )

8
3 8 3
0 0
1 6 1 6
5 5
2 2

Figure 3.8: Selected shape functions of the quadratic serendipity element

Serendipity elements of orders higher than two can also be constructed with the same number of boundary
nodes as their Lagrangian counterparts, but fewer interior nodes. These elements are therefore generally
more efficient than those of the Lagrange family. A disadvantage, however, is that they are more sensitive
to distortions. Their shape should therefore be controlled properly and should not deviate too wildly from
that of the – square – master element. Serendipity elements and Lagrange elements of the same order are
compatible and can therefore be mixed in the same mesh.

44
3.3.4 Numerical integration of quadrilateral elements

The appropriate Gauss integration rules for quadrilateral elements can be obtained as a natural extension
of the one-dimensional rules. The positions of the Gauss points may be expanded into the second di-
rection and the corresponding weight factors may be multiplied. This is illustrated in Figure 3.9, which
shows the Gauss point positions in the most commonly used Gauss quadrature schemes for quadrilateral
elements. Table 3.1 lists the positions and weight factors for these schemes.

4 7 3
4 3
1 8 9 6

1 2
1 5 2

Figure 3.9: Location of the Gauss points in the one (left), four (centre) and nine-point (right)
Gauss schemes for quadrilateral elements

s p (ξ1 )k (ξ2 )k wk
1 1 0 0 4
√ √
4 3 −1/ 3 −1/ 3 1
√ √
1/ 3 −1/ 3 1
√ √
1/ 3 1/ 3 1
√ √
−1/ 3 1/ 3 1
√ √
9 5 − 3/5 − 3/5 25/81
√ √
3/5 − 3/5 25/81
√ √
3/5 3/5 25/81
√ √
− 3/5 3/5 25/81
√
0 − 3/5 40/81
√
3/5 0 40/81
√
0 3/5 40/81
√
− 3/5 0 40/81
0 0 64/81

Table 3.1: Gauss point coordinates, weight factors and accuracy for two-dimensional quadrilat-
eral elements

For the bilinear, four-node elements a four-point integration is usually considered appropriate. Indeed,
if the element is rectangular and the coefficient tensor C(Ex ) is constant, this integration rule is exact.
This combination is therefore often referred to as full integration. For non-constant C(E x ) and non-
square elements a higher-order (say nine-point) rule would in principle be more accurate, but given the
approximations made in the discretisation of the problem this higher accuracy of the integration has only
a marginal impact on the overall accuracy of the solution, whereas it more than doubles the number of
evaluations of the integrand (nine vs. four). Although even more efficient, the one-point quadrature rule,
which is sometimes referred to as reduced integration for the four-node element, should generally be

45
avoided because it may lead to singular system matrices.

For quadratic elements – both the Lagrange and serendipity versions – full integration by the nine-point
Gauss rule is a safe choice. Reduced integration, which in this case means the four-point rule, is some-
times used for reasons of efficiency, but brings about a certain risk of singularity of the resulting matrix.
If this option is used, results should be examined carefully and a comparison with fully integrated results
is recommended.

3.3.5 Triangular elements

Triangular elements, or simply triangles, are useful because of their greater geometric flexibility com-
pared to quadrilaterals. This is particularly advantageous when automatic meshing techniques are used
on arbitrary domains. Unlike quadrilaterals, triangular elements cannot be constructed as straightfor-
ward extensions of the one-dimensional element family. However, a triangular isoparametric element
family can still be built along more or less the same lines as for quadrilaterals if we use so-called triangle
coordinates as the normalised, local basis.

The definition of triangle coordinates, which are also referred to as area coordinates, is illustrated for a
general triangle in Figure 3.10. The coordinate λ1 , associated with the corner indicated by 1 in the figure,
is defined as the area of the triangle between P and the edge opposite this corner – i.e. the shaded area
in the figure – divided by the total area of the triangle. It can easily be verified that this means that λ1
increases linearly from λ1 = 0 on the edge opposite of corner 1 to λ1 = 1 in this corner. Furthermore,
lines of constant λ1 are parallel to the edge opposite corner 1, as indicated by the dashed line in the
figure. In exactly the same way two more area coordinates, λ2 and λ3 can be associated with the second
and third corners of the triangle respectively. Any point in the triangle can be uniquely defined by the
three coordinates (λ1 , λ2 , λ3 ). Note, however, that two of these coordinates in principle suffice and the
three coordinates thus are not independent. Indeed, it is immediately clear from their definition that they
must satisfy λ1 + λ2 + λ3 = 1. It can also be seen directly that the centre of gravity of the triangle has
coordinates ( 13 , 13 , 13 ).

P
2
1
constant λ1

Figure 3.10: Definition of area coordinate λ1 for a general triangle; coordinates λ2 and λ3 are
defined fully analogically

Area coordinates are useful in defining triangular elements for a number of reasons. A perhaps slightly
academic reason is that they respect the threefold symmetry which a triangle may show and do not
introduce any form of preferential direction. This is reflected in the linear master element depicted in
Figure 3.11: an equilateral triangle with edge lengths equal to one. A second, more practical reason for
using area coordinates is that the shape functions of such elements can be constructed in a straightforward
way and take on simple forms when written in terms of the area coordinates. For instance, the triangular
coordinates themselves already satisfy the requirements one would impose on the shape functions of the

46
 eE2
3 3

eE1
2
1 2
1

Figure 3.11: Linear triangular master element (left) and possible physical shape (right)

linear triangle of Figure 3.11: they vary linearly in space and equal one in the the node they are associated
with and zero in the other two nodes. We can thus immediately write down
 
λ
e E
 1
N (ξ ) = λ2 
 (3.37)
˜ 
λ3

The first of these functions has been plotted in Figure 3.12. Note that, in contrast to the bilinear shape
functions of the ‘linear’ quad, these functions are truly linear. They also vary linearly along the element
boundaries and the element is therefore compatible with the four-node quadrilateral element. Indeed,
linear triangles are sometimes inserted in otherwise quadrilateral meshes for reasons of geometrical flex-
ibility.

3
N1e (ξE )

0
1
2

Figure 3.12: Shape function N1e (ξE ) associated with the linear triangular element

The shape functions of higher-order triangular elements can be constructed in a relatively straightforward
manner in terms of area coordinates. For instance, the six-node quadratic master element shown in
Figure 3.13 has shape functions
 
λ1 (2λ1 − 1)
 λ2 (2λ2 − 1)
 
 
 λ (2λ − 1)
 3 3
Ne(ξE ) =  (3.38)


˜ 4λ1 λ2



4λ2 λ3
 

4λ3 λ1

47
Two of these functions have been plotted in Figure 3.14. It can be verified that the quadratic triangle
provides a complete quadratic polynomial. In this sense, the element is more efficient than its quadrilat-
eral counterparts. As a result of the quadratic shape functions combined with the isoparametric mapping,
edges of the element may be curved according to a parabola (Figure 3.13). This element is compatible
with the biquadratic quadrilateral elements of the Lagrange as well as the serendipity families.

eE2
3 3

5
6 5 6

eE1
2
1 4 2 4
1

Figure 3.13: Quadratic triangular master element (left) and possible physical shape (right)

1 1
3 3
N4e (ξE )
N1e (ξE )

6
5 5
6
0 0
1 1
4 4
2 2

Figure 3.14: Selected shape functions of the quadratic triangle

The evaluation of element contributions to the system matrix and the right-hand side is slightly more
involved for the triangular master elements than for the quadrilateral elements. This is due to the fact
that the isoparametric mapping requires an independent set of local coordinates. The area coordinates
clearly do not satisfy this requirement. For this reason we introduce a Cartesian basis {Ee1 , eE2 } with its
origin in the centre of gravity of the master element and eE1 parallel to the first side of the triangle (see
Figures 3.11 and 3.13). However, it is emphasised that the orientation of this system does not have an
impact on the behaviour of the element and, in particular, does not introduce any preferential direction.
Points within the triangle can now be indicated by a vector ξE = ξ1 eE1 + ξ2 eE2 . Note that the existence
of such a coordinate system was already tacitly assumed in Equations (3.37) and (3.38). We can now
express the area coordinates λ1 , λ2 , λ3 in terms of the Cartesian coordinates ξ1 and ξ2 . Simple geometric
arguments show that these relations must read
1 √1 ξ2 1 √1 ξ2 1 √2 ξ2
λ1 = 3
− ξ1 − 3
λ2 = 3
+ ξ1 − 3
λ3 = 3
+ 3
(3.39)

Inverse relations may be established by solving these expressions for ξ1 and ξ2 .

Since for the linear element the element shape functions equal the area coordinates, the local gradient
E ξ Ne , is given directly by applying the gradient operator to expressions (3.39),
of the shape functions, ∇
resulting in ˜

E ξ Ne = ∇
∇ Eξ λ (3.40)
˜ ˜
48
where the area coordinates λ1 , λ2 , λ3 have been stored in a column matrix λ(ξE ). Since relations (3.39)
E ξ λ is a column of constant vectors. For elements ˜of order two and higher, the
are linear, the gradient of ∇
˜
local gradient of the shape functions can be written as
e
∇E ξ Ne = ∂N
˜ ∇ Eξλ (3.41)
˜ ∂λ ˜
˜
where the matrix ∂Ne/∂λ is defined as
˜ ˜
e  e
∂N ∂N ∂Ne e 
∂N
˜ = ˜ ˜ ˜ (3.42)
∂λ ∂λ1 ∂λ2 ∂λ3
˜
and can be evaluated relatively easily in terms of the area coordinates.

The local shape function gradient ∇ E ξ Ne having been determined, the Jacobian J(ξE ) can be computed via
(3.23) and the element integrals which ˜ form Ke and fe can be transformed as described in Section 3.2.4,
¯ ˜ are defined on the unit triangle T rather than the
the only difference being that the transformed integrals
square Q.

3.3.6 Numerical integration of triangular elements

For triangular elements, Gauss integration schemes are conveniently formulated in terms of the area
coordinates λ1 , λ2 and λ3 . The definition of the quadrature points and weight factors is largely analogical
to the case of quadrilateral elements. However, for triangles the quadrature rules (3.28) and (3.29) are
usually reformulated as
s
E e E eT
X
e
· C xE(ξEk ) · ∇N det J (ξEk )



K =T wk ∇N (3.43)
¯ ˜ ξE =ξEk ˜ ξE =ξEk
k=1
s
X
fe = T wk Ne(ξEk ) f xE(ξEk ) det J(ξEk )



(3.44)
˜ k=1
˜
√
where T = 14 3 denotes the area of the master element. With this definition, the sum of the weight
factors wk equals one.

Table 3.2 lists these weight factors for the most common triangle schemes. The locations of the respective
integration points are listed in the table in terms of their area coordinates; these locations are also shown
in Figure 3.15. For the linear triangle, the one-point integration scheme is usually used. For the quadratic
element the three-point scheme usually suffices.

3 3

1 3 2 4 6 7 5
1 2
1 4 2
1

Figure 3.15: Location of the Gauss points in the one (left), three (centre left) four (centre right)
and nine-point (right) Gauss schemes for triangular elements

49
 s p (λ1 )k (λ2 )k (λ3 )k wk
1 1 1/3 1/3 1/3 1
3 2 1/2 1/2 0 1/3
0 1/2 1/2 1/3
1/2 0 1/2 1/3
4 3 3/5 1/5 1/5 25/48
1/5 3/5 1/5 25/48
1/5 1/5 3/5 25/48
1/3 1/3 1/3 −27/48
7 5 0.7974269853 0.1012865073 0.1012865073 0.1259391805
0.1012865073 0.7974269853 0.1012865073 0.1259391805
0.1012865073 0.1012865073 0.7974269853 0.1259391805
0.4701420641 0.4701420641 0.0597158717 0.1323941527
0.0597158717 0.4701420641 0.4701420641 0.1323941527
0.4701420641 0.0597158717 0.4701420641 0.1323941527
1/3 1/3 1/3 9/40

Table 3.2: Gauss point coordinates, weight factors and accuracy for two-dimensional triangular
elements

3.4 Three-dimensional elements

Where model simplifications which lead to a two-dimensional representation cannot be made, three-
dimensional finite element discretisations must be used. We have seen examples of such cases in Fig-
ures 1.1 and 1.2 of the introductory chapter. Although constructing three-dimensional meshes and visu-
alising the results obtained with them is more difficult than in two dimensions, the finite element method
itself in three dimensions is an almost trivial extension of the two-dimensional method. In particular,
three-dimensional elements can be obtained in a straightforward fashion by extending in the third di-
rection the arguments which led to the two-dimensional families. For this reason, we will largely limit
ourselves below to indicating how the extension from two to three dimensions can be made, without
giving much detail.

3.4.1 Hexagonal elements

The three-dimensional counterparts of the two-dimensional quadrilateral element families are hexagonal
elements, sometimes also called brick elements or bricks. The lowest-order brick element, which is

50
shown on the left in Figure 3.16, has eight nodes and trilinear shape functions

1 
8
(1 − ξ1 )(1 − ξ2 )(1 − ξ3 )
1
 8 (1 + ξ1 )(1 − ξ2 )(1 − ξ3 )

 
 1 (1 + ξ1 )(1 + ξ2 )(1 − ξ3 )
8 
1
− ξ1 )(1 + ξ2 )(1 − ξ3 )

 8 (1
Ne(ξE ) = 
 1 (1
 (3.45)
˜ 8 − ξ1 )(1 − ξ2 )(1 + ξ3 )

1 
 (1
8 + ξ1 )(1 − ξ2 )(1 + ξ3 )

1
 8 (1 + ξ1 )(1 + ξ2 )(1 + ξ3 )

1
8
(1 − ξ1 )(1 + ξ2 )(1 + ξ3 )

8 8 8
15 15
7 16 7 16 7
e3 26
5 5 14 5 14
e2 13 20 23 13 20
6 6
24 19 19
6 eE1 4 27
17 11 22 17 11
4 21 4
12 12
18 18
3 3 3
25
1 1 10 1 10
9 9
2 2 2

Figure 3.16: Hexagonal master elements: linear (left), quadratic Lagrangian (centre) and
quadratic serendipity (right)

Quadratic and higher-order elements again come in two families. On the one hand, the Lagrange family
is obtained by straightforward extension of the one-dimensional polynomials. This leads to a brick
element with 27 nodes, as depicted at the centre of Figure 3.16. The shape functions of this element
can be constructed without much difficulty and are therefore not detailed here. The quadratic serendipity
element (on the right in Figure 3.16) is generally favoured over the Lagrange element because it provides

51
fewer incomplete polynomial terms. It has only 20 nodes, with shape functions
 1 
− 8 (1 − ξ1 ) (1 − ξ2 ) (1 − ξ3 )(2 + ξ1 + ξ2 + ξ3 )
 1
− 8 (1 + ξ1 ) (1 − ξ2 ) (1 − ξ3 )(2 − ξ1 + ξ2 + ξ3 )

 
− 1 (1 + ξ ) (1 + ξ ) (1 − ξ )(2 − ξ − ξ + ξ )
 8 1 2 3 1 2 3 
 1
− 8 (1 − ξ1 ) (1 + ξ2 ) (1 − ξ3 )(2 + ξ1 − ξ2 + ξ3 )

 
− 1 (1 − ξ ) (1 − ξ ) (1 + ξ )(2 + ξ + ξ − ξ )
 8 1 2 3 1 2 3 
 1 
− (1 + ξ1 ) (1 − ξ2 ) (1 + ξ3 )(2 − ξ1 + ξ2 − ξ3 )
 8 
 1
− 8 (1 + ξ1 ) (1 + ξ2 ) (1 + ξ3 )(2 − ξ1 − ξ2 − ξ3 )

 
− 1 (1 − ξ1 ) (1 + ξ2 ) (1 + ξ3 )(2 + ξ1 − ξ2 − ξ3 )
 8 
 1 2
 4 (1 − ξ1 )(1 − ξ2 ) (1 − ξ3 )


 
 1 (1 + ξ )(1 − ξ 2 )(1 − ξ ) 
1 3
Ne(ξE ) =  41 2
(3.46)
 
 (1 − ξ12 )(1 + ξ2 ) (1 − ξ3 )

˜  4


 1 2
 4 (1 − ξ1 )(1 − ξ2 )(1 − ξ3 )


 1 
 (1 − ξ 2 )(1 − ξ2 ) (1 + ξ3 ) 
 4 1 
 1
 4 (1 + ξ1 )(1 − ξ22 )(1 + ξ3 )


 
 1 (1 − ξ 2 )(1 + ξ ) (1 + ξ ) 
 4 1 2 3 
 1
 (1 − ξ1 )(1 − ξ22 )(1 + ξ3 )


 4 
 1
 4 (1 − ξ1 ) (1 − ξ2 )(1 − ξ32 )


 1 
 (1 + ξ1 ) (1 − ξ2 )(1 − ξ 2 ) 
 4 3 
 1
 4 (1 + ξ1 ) (1 + ξ2 )(1 − ξ32 )


1
4
(1 − ξ1 ) (1 + ξ2 )(1 − ξ32 )

Table 3.3 lists the positions and weight factors which can be used in the numerical integration of the
element stiffness matrix and right-hand side. For the eight-node brick element the eight-point integra-
tion scheme should normally be used. For the quadratic elements of the Lagrange type as well as the
serendipity type, twenty-seven Gauss points are a safe choice, but the reduced eight-point scheme is also
sometimes used.

3.4.2 Tetrahedral elements

Tetrahedral elements, tetrahedra, or simply tets are the three-dimensional counterparts of triangular ele-
ments. Figure 3.17 shows the linear and quadratic master elements of this class. These master elements
consist √
of a regular tetrahedron with an edge length equal to one; their volume is therefore given by
1
T = 12 2.

The shape functions of tetrahedra are conveniently formulated in terms of volume coordinates λ1 , λ2 , λ3 ,
λ4 . These are defined fully analogically to the area coordinates for triangles. The volume coordinate λ1
of a point P is equal to the ratio of the volume of the tetrahedron which can be constructed between P
and the face opposite of the first corner and the entire volume of the tetrahedron on which it is defined,
see Figure 3.18 for a graphical representation of these two volumes. From this definition, it can easily
be seen that the plane of constant λ1 is a plane parallel to the face opposite of corner 1. The other three
volume coordinates are defined similarly. Like the area coordinates on triangles, they are not entirely
independent, but satisfy the relation λ1 + λ2 + λ3 + λ4 = 1.

52
s p (ξ1 )k (ξ2 )k (ξ3 )k wk
1 1 0 0 0 8
√ √ √
8 3 −1/ 3 −1/ 3 −1/ 3 1
√ √ √
1/ 3 −1/ 3 −1/ 3 1
√ √ √
1/ 3 1/ 3 −1/ 3 1
√ √ √
−1/ 3 1/ 3 −1/ 3 1
√ √ √
−1/ 3 −1/ 3 1/ 3 1
√ √ √
1/ 3 −1/ 3 1/ 3 1
√ √ √
1/ 3 1/ 3 1/ 3 1
√ √ √
−1/ 3 1/ 3 1/ 3 1
√ √ √
27 5 − 3/5 − 3/5 − 3/5 125/729
√ √ √
3/5 − 3/5 − 3/5 125/729
√ √ √
3/5 3/5 − 3/5 125/729
√ √ √
− 3/5 3/5 − 3/5 125/729
√ √ √
− 3/5 − 3/5 3/5 125/729
√ √ √
3/5 − 3/5 3/5 125/729
√ √ √
3/5 3/5 3/5 125/729
√ √ √
− 3/5 3/5 3/5 125/729
√ √
0 − 3/5 − 3/5 200/729
√ √
3/5 0 − 3/5 200/729
√ √
0 3/5 − 3/5 200/729
√ √
− 3/5 0 − 3/5 200/729
√ √
0 − 3/5 3/5 200/729
√ √
3/5 0 3/5 200/729
√ √
0 3/5 3/5 200/729
√ √
− 3/5 0 3/5 200/729
√ √
− 3/5 − 3/5 0 200/729
√ √
3/5 − 3/5 0 200/729
√ √
3/5 3/5 0 200/729
√ √
− 3/5 3/5 0 200/729
√
0 0 − 3/5 320/729
√
0 0 3/5 320/729
√
0 − 3/5 0 320/729
√
3/5 0 0 320/729
√
0 3/5 0 320/729
√
− 3/5 0 0 320/729
0 0 0 512/729

Table 3.3: Gauss point data for three-dimensional hexagonal elements

53
 4 4

8 10 9
3
7 3
6
1 2 1 5 2

Figure 3.17: Tetrahedral elements: linear (left) and quadratic (right)

P
3
1

Figure 3.18: Definition of volume coordinate λ1 for a general tetrahedron; coordinates λ2 , λ3

and λ4 are defined fully analogically

Like for the triangular elements, the isoparametric mapping requires that we also define an independent
coordinate system. A basis {Ee1 , eE2 , eE3 } is introduced for this purpose, with its origin in the centre of
gravity of the tetrahedron. This centre of gravity is given in volume coordinates by ( 14 , 14 , 14 , 14 ). We
define the basis vector eE1 parallel to the edge between corners 1 and 2, eE2 perpendicular to eE1 but parallel
to the face between corners 1, 2 and 3, and eE3 perpendicular to this plane in the direction of corner 4.
However, it is emphasised that the orientation of this basis does not influence the behaviour of the element
in any way. For the regular tetrahedron which forms our master element, the volume coordinates can now
be written in terms of the Cartesian components (ξ1 , ξ2 , ξ3 ) of a local position vector ξE as

1
λ1 = 1
− ξ1 − √1 ξ2 − √1 ξ3 λ2 = 4
+ ξ1 − √13 ξ2 − √1 ξ3
6
4 3 6
1
q (3.47)
λ3 = + √2 ξ2 − √1 ξ3 λ4 = 1
+ 32 ξ3
4 3 6 4

E ξ Ne can be determined via (3.41) and the isoparametric transforma-

Using these relations, the gradient ∇
tion can go ahead as described before.˜

The shape functions of the linear tetrahedral element depicted on the left in Figure 3.17 read (cf. (3.37)):
 
λ1
 
λ2 
Ne(ξE ) =   (3.48)
˜ λ 
 3
λ4

These functions again provide a truly linear interpolation between the four nodes of the element. The

54
quadratic element (on the right in the figure) has ten nodes and shape functions
 
λ1 (2λ1 − 1)
 
 λ2 (2λ2 − 1)
 
 λ3 (2λ3 − 1)
 
 
 λ (2λ − 1)
 4 4 
 
4λ1 λ2
Ne(ξE ) = 
 
 (3.49)
˜ 4λ2 λ3
 

 
4λ3 λ1 
 
4λ1 λ4
 

 
4λ2 λ4 
 
4λ3 λ4

This elements provides a quadratic interpolation within the element. It can be verified that the shape
functions of the linear element are compatible with those of the trilinear hexagonal element; these ele-
ments can thus be mixed in finite element meshes. Likewise, the quadratic tetrahedron can be combined
with the 20 and 27-node hexagons.

The numerical integration of tetrahedral elements

√ follows Equations (3.43)-(3.44), with T equal to the
1
volume of the master tetrahedron: T = 12 2. Table 3.4 lists the most common Gauss schemes. The
locations of the Gauss points are given in terms of volume coordinates. For the linear element the one-
point scheme is often used, whereas the four-point scheme is required for the quadratic element.

s p (λ1 )k (λ2 )k (λ3 )k (λ4 )k wk

1 1 1/4 1/4 1/4 1/4 1
4 2 0.58541020 0.13819660 0.13819660 0.13819660 1/4
0.13819660 0.58541020 0.13819660 0.13819660 1/4
0.13819660 0.13819660 0.58541020 0.13819660 1/4
0.13819660 0.13819660 0.13819660 0.58541020 1/4
5 3 1/2 1/6 1/6 1/6 9/20
1/6 1/2 1/6 1/6 9/20
1/6 1/6 1/2 1/6 9/20
1/6 1/6 1/6 1/2 9/20
1/4 1/4 1/4 1/4 −4/5

Table 3.4: Gauss point data for three-dimensional tetrahedral elements

55
 4

Vectorial boundary value problems

4.1 Model boundary value problem

Many boundary value problems in engineering and science cannot be formulated in terms of a single,
scalar field like we considered in the previous chapter, but contain a vector field. Most notable examples
are the displacement field in solid mechanics and the velocity field in fluid mechanics. Note that after
choosing a suitable vector basis the three (or two) components of the vector field with respect to this
basis can be also be regarded as separate, scalar fields. The vector equation which they must satisfy then
becomes a set of coupled scalar partial differential equations which must be solved simultaneously. This
is the view which has traditionally been taken when developing finite element formulations and finite
element codes. Modern software tools, however, allow to directly code in terms of vectors and tensors.
As a result, finite element formulations can be developed in a vector format, which can then be translated
almost literally into computer code. This not only leads to more transparent mathematical formulations
and computer codes, it also allows to code new problems only once for a range of coordinate systems.
Given these advantages, we will largely concentrate on vector-based descriptions here; we will come
back briefly to the equivalent scalar formulations at the end of the chapter.

The partial differential equation which we consider as a model problem in this chapter reads
 
∇E · 4C(Ex):∇E uE + fE(Ex ) = 0E (4.1)

This equation can immediately be recognised as that governing the linear elastic deformation of solids,
with uE(E x ) the elasticity tensor and fE(E
x ) the displacement field, 4C(E x ) a body force – e.g. gravity. Associ-
ated with (4.1) are essential and natural boundary conditions, which in the elasticity interpretation would
mean imposed displacements and imposed boundary tractions respectively. The essential condition reads

uE(E
x ) = gE(E
x) (4.2)

and is valid on a part Sg of the boundary S. On the remaining part Sq we have the natural condition

x )· 4C(E
nE (E E uE = qE(E
x) : ∇ x) (4.3)

4.2 Finite element discretisation

The steps which must be taken to rework the boundary value problem given by Equations (4.1)–(4.3)
into a set of linear algebraic equations are entirely equivalent with those of Chapter 3. Some technical
differences arise, however, as consequence of the vectorial character of the equation we are considering
here. For this reason, each of the steps is reiterated here with special emphasis on the differences with
the scalar case.

4.2.1 Weak form of the equation

The weighted residuals form of Equation (4.1) follows by taking the scalar product with a vector-valued
E x ) followed by integration on the domain V :
test function φ(E
Z h   i
φ(E E · 4C(E
E x )· ∇ E uE + φ(E
x) : ∇ E x )· fE(E
x ) dV = 0 (4.4)
V

59
Integration by parts of the first term in the integrand – or, more precisely, using the product rule for the
divergence and the divergence theorem – gives
Z   Z Z
E E 4 E E E φE T: 4C(E E uE dV


φ(Ex )· ∇ · C(Ex ) : ∇ uE dV = φ(E x )· qE (E
x ) dS − ∇ x):∇ (4.5)
V S V

where use has also been made of the definition

qE (E x )· 4C(E
x ) = nE (E E uE
x) : ∇ (4.6)

on the part Sg of S where this relation did not already hold because of the natural boundary condition
(4.3); it is emphasised, however, that qE (E
x ) is yet unknown on Sg . Substitution of (4.5) into (4.4) and
reordering results in the weak form
Z Z Z
E φE T : 4C(E E x )· fE(E
E uE dV = φ(E E x )· qE (E


∇ x):∇ x ) dV + φ(E x ) dS (4.7)
V V S

4.2.2 Discretisation of the weak form

Like for the scalar problem of the previous chapter, the problem domain V in (4.7) is now replaced by
a domain V h which is constructed out of a number of finite elements. Likewise, the boundary S h of V h
replaces S in (4.7). For simple geometries V and V h may be identical, but for more complex shapes this
discretisation of the problem domain involves some approximation (Figure 3.2). At the same time the
x ) is replaced by an approximate field uEh (E
vector field uE(E x ) of the form
n
X
h
x) =
uE (E x ) uEi
Ni (E (4.8)
i=1

The shape functions Ni (Ex ) are scalar-valued. They are used to interpolate the uEh(E
x )-field on the basis of
a number of nodal values uEi . Since we are dealing with a vector field, each of these nodal values takes
the shape of a vector. Like before, it is useful to write the above expression in a more compact fashion
using the matrix notation

x ) = NT(E
uEh(E x ) uE (4.9)
˜ ˜
In this expression uE is a column matrix of vectors, i.e. a column matrix which has vectors as its entries
˜ scalar entries. Apart from this minor difference, (4.9) is entirely equivalent with
rather than the usual
E x ) is replaced by the
(3.13) in the previous chapter. Following the Galerkin approach, the test function φ(E
similarly discretised field

φEh(E
x ) = NT(E E
x) φ (4.10)
˜ ˜
E also is a column matrix of vectors which contains the values of φE h(E
where φ x ) in the nodes.
˜
E
Substituting the above discretisations into the weak form and realising that the result must hold for all φ
we have ˜
Z Z Z
E 4
∇N · C(E E T
x )· ∇N dV · uE = N(E E x ) dV + N(E
x ) f (E x )E x ) dS
q (E (4.11)
˜ ˜ ˜ ˜ ˜
Vh Vh Sh

60
Inspecting the integral on the left-hand side of this equation shows that it forms an n × n matrix of
second-order tensors. If we denote this tensor matrix by K, i.e.
¯
Z
K = ∇N E · 4C(E E T dV
x )· ∇N (4.12)
¯ ˜ ˜
Vh

and furthermore define the vector columns

Z
Ef = N(E x ) fE(E
x ) dV (4.13)
˜ ˜
Vh
Z
qE = N(E x )E x ) dS
q (E (4.14)
˜ ˜
Sh

Equation (4.11) may be rewritten as

K · uE = Ef + qE (4.15)
¯ ˜ ˜ ˜

Comparing Equation (4.15) with Equations (3.16) and (2.20) shows a great deal of similarity between
the systems of equations obtained for the one-dimensional and multi-dimensional scalar boundary value
problems and the present problem in terms of a vector-valued field. The main difference is in the fact
that the unknowns of the system are now vectors rather than scalar quantities, and related to this, the fact
that the system matrix now contains second-order tensors rather than scalar values and that the product
between this matrix and the the column uE has become an inner product. The tensorial nature of the system
matrix arises in a natural way from the˜ fact that the coefficient tensor 4C(Ex ) is a fourth-order tensor in
(4.12) instead of the second-order tensor in (3.17). Similarly, the right-hand side matrices Ef and qE have
˜
become vector columns as a consequence of the vectorial character of fE(E x ) respectively.˜
x ) and qE (E

It is interesting to note that for boundary value problems which are formulated in terms of tensors exactly
the same steps as taken above lead to discrete systems in terms of tensor-valued unknowns and system
matrices which consist of higher-order tensors. For instance, a second-order tensor field will give rise to
second-order tensor column matrices and a system matrix which contains fourth-order tensors; accord-
ingly, the product in the linear system will become a double inner product. However, such problems are
not very common in engineering and we will therefore not pursue them any further.

4.2.3 Shape functions, isoparametric transformation and numerical integration

Apart from the fact that the system matrix K and the right-hand side matrices Ef and qE have got tensors and
¯ ˜ ˜ virtually identical
vectors as their entries respectively, the expressions for these matrices have remained
to those for the scalar problem, cf. Equations (3.17)–(3.19). As a result, the requirements which must
be satisfied by the shape functions contained in N(E x ) are completely identical to those arising in the
scalar problem of the previous chapter. This means ˜ that the shape functions and elements which were
introduced in Sections 3.3 and 3.4 can be adopted here without any change.

Likewise, the arguments which led to the use of the isoparametric transformation and numerical inte-
gration to evaluate the element contributions to the global system in the scalar case equally apply to the
present situation. Reasoning fully parallel to the discussion in Section 3.2.4, this leads to expressions for

61
the contributions of element e to K, Ef and qE which read (cf. Equations (3.28), (3.29), (3.31)):
¯ ˜ ˜
Xs
e E e
· C xE(ξEk ) · ∇N
4 E eT det J(ξEk )



K = wk ∇N (4.16)
¯ k=1
˜ ξE=ξEk ˜ ξE=ξEk
s
Efe = wk Ne(ξEk ) fE xE(ξEk ) det J (ξEk )
X


(4.17)
˜ k=1
˜
s ′
X
e
wk′ Ne′(ξEk′ ) qE xE(ξEk′ ) det J ′ (ξEk′ )



qE = (4.18)
˜ k=1
˜
where the primes in the latter expression again denote quantities which are defined with respect to the
relevant part of the element boundary. For triangular and tetrahedral elements, the element area or volume
must be inserted in the above expressions, cf. Equations (3.43)–(3.44). The Gauss point data given in
Tables 3.2 and 3.4 can then be used without any further change. For quadrilateral and hexagonal elements,
the data given in Tables 3.1 and 3.3 respectively can be used straightaway.

4.2.4 Assembly and partitioning

e
Having computed the element contributions Ke , Ef and qEe for each element e as described above, the
¯ ˜ ˜
global system can be assembled:
m m m
Efe
K=
¯ Ae=1
Ke
¯
Ef =
˜ e=1
A qE =
˜˜
qEe
e=1 ˜
A (4.19)

These assembly operations are identical to those for the scalar case, with the understanding that they act
on tensorial and vectorial components rather than scalar values. This also holds for the partitioning of
the resulting global system of equations, which results in a reduced set of equations which reads
Kff · uEf = Eff + qE f − Kfg · uE g (4.20)
¯ ˜ ˜ ˜ ¯ ˜
Like for the scalar case, the above operations can be executed by a computer in a relatively straightfor-
ward way. Indeed, if the programming platform which is used allows to define vectors and tensors, the
only difference with the scalar case is in the type of objects on which the operations act.

4.2.5 Vector formulation vs coupled scalar formulation

It should be mentioned, finally, that a scalar equivalent of the vector equation (4.15) can be obtained
if a vector basis is introduced right at the start of the analysis. As indicated at the beginning of this
chapter, the vector-form of the boundary value problem may then be rewritten as a coupled set of three
scalar partial differential equations (or two in a two-dimensional setting). Each of these equations can
be discretised along the lines of Chapter 3 and the resulting sets of equations can be merged into one
linear system, which will be equivalent with (4.15). Indeed, this is the perspective which is taken in
many textbooks and computer codes. However, a major disadvantage of specialising the development
– and indeed the implementation – of finite element formulations at such an early stage to a particular
coordinate basis, is that it needs to be repeated for every new basis one might want to select. By contrast,
the vector–tensor notation used in these lecture notes allows to delay the specification of a basis – or
even multiple bases – right until the final stages of development or. As a result, the development of the
method itself can be done once and for all and covers all possible bases. Furthermore, as an important
by-product, the presentation of the method becomes more transparent and is more easily recognised as
an extension of the scalar and one-dimensional cases.

62