1 A finite-element formulation of the linear elasticity
1.1 Initial definitions
In the following, several conventions are used to facilitate the layout:
• We distinguish between an equation (using = symbol) and definition (using ≡ symbol)
• Summation over repeated indices is assumed in any term in an equation or on the right side in a
definition. Let’s take an equation:
Aijxj = bi where i = 1 . . . N, j = 1 . . . M (1)
denotes a system of N linear equations where summation over j on the left is assumed:
M
j=1
Aijxj = bi (2)
Now let’s take a definition:
Ci ≡ Dijfj where i = 1 . . . N, j = 1 . . . M (3)
denotes a definition of a vector C of size N. Each element of this vector is defined as:
Ci ≡
M
j=1
Dijfj. (4)
• Derivative is denoted by a comma notation. So
fi,j ≡
∂fi
∂xj
where i = 1 . . . N, j = 1 . . . M (5)
where f is a vectorial function (of dimension N) of M variables x1, . . . xN and here we take a partial
derivative of each component fi w.r.t. each variable xj. Therefore, if f is linear in each component,
the term fi,j denotes a matrix of size N × M.
• In all what follows in this document, we suppose that indices i, j, k are dimension-related indices. Since
we focus on 3D elasticity, each of them is running over a set {1, 2, 3}. For example xi is an element of
a vector {x1, x2, x3} (an alternative for {x, y, z}).
• A combination of all that has been said can be shown on the following example:
divσ ≡ · σ ≡ σij,j ≡
j
∂σij
∂xj
(6)
where σij is a 3×3 matrix, while σij,j is a column vector of 3 elements (since each element of σij is
derived w.r.t. xj and summation over j is applied.
1
1.2 Linear Elasticity in 3D
We recall that the linear elasticity is given by following kinetic equation:
· σ + b = 0 (7)
where · is an alternative notation for a divergence of vector field, σ is the stress tensor (3×3 matrix) defined
as
σ ≡ λtr(e)I + 2µe (8)
where λ and µ are the Lam´e coefficients and e is the linearized Green strain tensor (3×3 matrix) defined as:
e ≡
1
2
( u + u ). (9)
If written in tensor notation, the definition of the 3D linear elasticity problem is:
σij,j + bi = 0 (10)
where
σij ≡ λekkδij + 2µeij (11)
where
eij ≡
1
2
(ui,j + uj,i). (12)
Notice that only the first one is a “real equation”, the other two are just definitions of terms (stress and
strain tensors). Also notice that we are still in a “continuous realm”: e.g. the unknown variable u are
(presumably continuous) functions of x1, x2, x3. The same is true for the applied loads (bi) and so σij,j can
be regarded as a functional operator (in this case a differential operator) applied to functions.
1.3 Application of the finite element (FE) method
The finite element method reformulates a problem given by a differential operator (such as the one in Eq. 10)
as a set of algebraic equations which can be solved numerically. Basically, instead of looking for an unknown
function (e.g. a solution for an infinite number of degrees of freedom [DoFs]), we reformulate the problem so
that we solve for an unknown vector (having a finite number of DoFs) which approximates the continuous
solution given relevant interpolation functions. The (non-)linearity of the algebraic equations depends on
the (non-)linearity of the initial differential operator.
The FE method has following basic steps:
• Weak reformulation of the problem (mainly the continuous differential operator).
• Discretization of the domain and unknown variables. As a necessary input, the type of elements
together with the associated interpolation functions must be chosen here, resulting in the local system
of algebraic equations assembled for each element.
• Assembling of the global system of equations from the local element equations.
• Application of the boundary conditions.
• Numerical solution of the system.
In what follows, we apply the FE method on the 3D linear-elasticity problem defined by Eq. 10, focusing on
the first two steps.
2
1.3.1 Weak formulation
Introduction from Wikipedia:
Weak formulations are an important tool for the analysis of mathematical equations that permit the
transfer of concepts of linear algebra to solve problems in other fields such as partial differential equations.
In a weak formulation, an equation is no longer required to hold absolutely (and this is not even well defined)
and has instead weak solutions only with respect to certain ”test vectors” or ”test functions”.
In 3D linear elasticity, the weak (re-)formulation of the problem allows for redistribution of the derivatives
and incorporation of the boundary conditions.
The weak form is obtained in two steps: the problem is integrated over the domain Ω where the problem is
solved (e.g. the volume of the object undergoing the deformations) and multiplication of the test function (a
sort of a “weight” function). Intuitively, the original problem where the equality given by Eq. 10 is required to
hold at each point of the domain is replaced by a weaker condition, where only a weighted sum of integrated
terms is to be required to hold in the equality. The validity (together with corresponding limitations) of
such a weakening of the problem is shown by Lax-Milgram theorem (see corresponding literature or web).
In our case, the weak form of the original problem 10 is simply:
Ω
(σij,j + bi)widΩ = 0 (13)
where Ω is the problem domain (in 3D so we are assuming a volume integration) and wi are the test functions.
The weak form allows us to redistribute the derivatives using a corollary of divergence theorem (see the first
corollary at http://en.wikipedia.org/wiki/Divergence_theorem), resulting in:
Ω
σijwi,jdΩ =
Ω
biwidΩ +
∂Ω
tiwidΓ (14)
where a new term appears on the right: the integration is taken over the boundary of the volume (denoted
as ∂Ω, sometimes replaced by Γ) and ti are tractions applied to the surface of the object.
1.3.2 Discretization
This is the key step in the FE method: the continuous problem is approximated by discretized one. First,
the continuous domain is discretized by finite elements, denoted by index e:
Ω ≡ Ωe (15)
where Ωe is the domain of an element e.
It is supposed that the element is a geometric primitive such as triangle, quad (in 2D), tetrahedron,
hexahedron, pyramid, prism (in 3D): i.e. the element is given by a small number of points called nodes. Besides
the geometrical description of the element, interpolation (alternatively called shape, basis or blending)
functions are associated to each node of the element. The interpolation functions can be used to approximate
(interpolate) any quantity (displacement, load, stress, strain, temperature, pressure etc.) in any position
inside of the element provided the values in the nodes are given.
In the following, an interpolation function associated to an element and a node is denoted ϕen
. Since we
are in 3D, ϕne
is a scalar function of x1, x2, x3. Two basic properties are required:
• ϕne
has a local support on the element e, i.e., it is non-zero only inside the associated element e and
zero anywhere else:
loc suppe ϕne
(16)
3
• ϕne
is bounded to the node n of the element e:
ϕne
(xm
) = δnm (17)
where xm
is the location of a node m and δnm is a Kronecker symbol.
Now, the discretization can be applied to the weak form 13: first, instead of having one equation over the
volume Ω, we have a sum of terms corresponding to each element Ωe and the physical quantities (unknown
displacements, stresses, strains etc.) are now related to local elements instead of the global domain. Second,
the continuous unknown functions ue
are approximated with interpolation functions:
ue
i ≈ ϕne
Une
i (18)
where Une
is an (unknown) displacement of the node n in element e. Since the number of elements is finite
and each element has a small number of nodes, Une
is a vector of finite size where each number in the vector
corresponds to one DoF of the discretized domain. Notice, that in Eq. 18, summation over n is assumed,
whereas no summation is done over e and i (which are taken as the parameters of the definition, although
symbol ≈ is used instead of ≡). Also notice that now the derivative of the function u can be approximated
as follows:
ue
i,j ≈ ϕen
,j Uen
i ≡
∂ϕen
∂xj
Uen
i (19)
Finally, the Galerkin method is used (for the sake of simplicity, the method is not described in details, but
it’s only used as a “black-box”). This means that the same interpolation functions ϕne
(used to approximate
the unknown displacements u) are employed as the test functions wi:
wen
i ≡ ϕne
Wne
i (20)
where Wne
i are virtual displacements (note that no summation is done in this case as all the indices appear
on left).
Putting it all together, the weak form is now discretized as:
e Ωe
σijϕen
,j dΩ Wen
i =
e Ωe
biϕen
dΩ +
∂Ωe
tiϕen
dΓ Wen
i (21)
Since the equality must hold for any virtual displacements Wne
i , the above sum can be rewritten as a
system of equalities where for each element e we have:
Ωe
σe
ijϕen
,j dΩ =
Ωe
be
i ϕen
dΩ +
∂Ωe
te
i ϕen
dΓ (22)
where the discretized stress tensor σe
ij is defined as
σe
ij = λϕne
,k Une
k δij + µ(ϕen
,j Uen
i + ϕen
,i Uen
j ) (23)
where we already used the definitions Eq. 11,12 and the discretization (Eq. 18, 19).
For the sake of simplicity, we assume the tractions to be zero everywhere (i.e. te
i = 0) and we assume that
loads be
i do not depend on the position so it can be moved outside of the integral. This results in simplifying
the right-hand side of 22. Thus, after substituting 23 into 22, the discretized formulation over an element e
becomes:
Ωe
λϕne
,k Une
k δij + µ(ϕen
,j Uen
i + ϕen
,i Uen
j ) ϕen
,j dΩ = bi
Ωe
ϕne
dΩ (24)
Notice that
4
• the term on the left-hand side is linear in U and only derivatives of the interpolation functions ϕne
appear there,
• the right-hand side represent applied body loads such as the gravity. In this case, the vector b is
computed as
b = ρ [0 − g 0] (25)
where ρ is the density and g is the gravitational acceleration, here applied along −y axis.
1.3.3 Particular case: P1 elements (linear tetrahedra)
The P1 elements are the simplest elements in 3D: each element has four nodes (therefore in all what follows,
n =∈ {1, . . . 4}) and the interpolation functions are linear. The simplest linear function in 3D is defined as:
ϕ(x1, x2, x3) = a + bx1 + cx2 + dx3 (26)
In what follows, we focus on a single element e (therefore, we can drop the index e in rest of this section).
This element has four nodes with positions xn
≡ {xn
1 , xn
2 , xn
3 }, n ∈ {1...4}. Now, in order to rewrite the
Eq. 24 (which is defined generally for any element of any order) for the given P1 element, we need to
obtain four linear interpolation functions ϕn
, i.e. we need to determine the coefficients an
, bn
, cn
, dn
(from
Eq. 26). This can be done using the property of the interpolation functions defined by Eq. 17. In practical
implementations, this is achieved by constructing a nodal matrix Ve:




1 x1
1 x1
2 x1
3
1 x2
1 x2
2 x2
3
1 x3
1 x3
2 x3
3
1 x4
1 x4
2 x4
3



 (27)
and solving for linear systems (one for each node, corresponding to one vector of the right-hand side):




1 x1
1 x1
2 x1
3
1 x2
1 x2
2 x2
3
1 x3
1 x3
2 x3
3
1 x4
1 x4
2 x4
3








a
b
c
d



 =




1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1



 (28)
While the tensor notation was extremely useful in the derivation of weak forms and general equations
(such as Eq. 24), the standard matrix (or vectorial) notation is more useful for the implementation purposes.
Moreover, Voigt notation is usually employed in order to avoid higher-order tensors: for example, a 2nd order
tensor in 3D (such as the stress tensor σij) is a symmetric 3×3 matrix. In order to reduce its dimensionality,
it can be stored as a vector of size 6, using the trick depicted in Fig. 1. Similarly, the unknown vector Un
i ,
which is in fact a 4×3 matrix (for each of four nodes, we have three displacements along x1, x2, x3), can be
stored as 12×1 vector. This process can be described by introducing a mapping on the indices where the
nodal index n and the dimension-related index i (or j) are joined together; thus the new “joined” index p
would be computed as p = 3n + i (we are not going to use the joined index, it’s only to explain better what
is behind these manipulations on the implementation level).
Putting this all together, the general 3D FE-discretized linear elasticity equation (Eq. 24) for an element
e can be written as
Ke
ue
= fe
(29)
5
Figure 1: Example of the Voigt notation
where ue
is a 12×1 vector gathering all the unknowns Un
i (displacements in each degree of freedom), similarly
b is a 12×1 vector gathering the integrated loads and Ke
is a 12×12 stiffness matrix. This is possible thanks
to the fact that the formulation given by Eq. 24 is linear in U, and so all the terms on the left-hand side can
be reorganized (using the Voigt notation and the index mapping introduced above) in the matrix Ke
which
is of form:
Ke ≡
Ωe
Be DeBedΩ (30)
where Be is a geometrical matrix containing only the coefficients obtained by solving the system 28:
Be =








b1 0 0 b2 0 0 b3 0 0 b4 0 0
0 c1 0 0 c2 0 0 c3 0 0 c4 0
0 0 d1 0 0 d2 0 0 d3 0 0 d4
c1 b1 0 c2 b2 0 c3 b3 0 c4 b4 0
d1 0 b1 d2 0 b2 d3 0 b3 d4 0 b4
0 d1 c1 0 d2 c2 0 d3 c3 0 d4 c4








(31)
and De is a material matrix containing the material parameters:
De =








λ + 2µ λ λ 0 0 0
λ λ + 2µ λ 0 0 0
λ λ λ + 2µ 0 0 0
0 0 0 µ 0 0
0 0 0 0 µ 0
0 0 0 0 0 µ








(32)
The final remark concerns the integration: since only derivatives of the interpolation function appear on the
left hand side of the Eq. 24, only the coefficients b, c and d of the interpolation functions are present in the
final matrix form (spatial derivatives of the linear shape functions). Therefore, there is no term inside K
which depends on the spatial coordinates and therefore, the integration over the volume of the P1 element
6
can be replaced by the multiplication by the element volume. For linear tetrahedron, the volume can be
computed using the determinant of the vertex matrix Ve (see Eq. 27). Therefore:
Ke ≡
Ωe
Be DeBedΩ ≡
det(Ve)
6
Be DeBe (33)
However, on the right hand side, the shape functions ϕne
must be integrated over the element volume.
Since the shape functions depends on positional variables, it is not possible to replace the integration by the
multiplication by the element volume. Instead, Gauss integration rule is applied to perform the numerical
integration. I.e., the integrand is evaluated in g integration points xp and the integrated value is obtained
as a weighted sum:
Ωe
ϕne
(x)dΩ ≈
g
p=1
ωpϕne
(xp) (34)
where ωp are weights associated to each integration point. The values of points and weights can be found
for given volume element.
The same approach must be applied to the left hand side if the integrand depends on the spatial variables.
This is the case when a higher-order shape functions are used (i.e., quadratic functions).
1.4 Final remarks
The previous section gives a complete set of instructions how to implement a linear elasticity FEM. However,
before the system 29 can be solved (by direct or iterative solver), two steps are needed:
• global assembly of the system: the Eq. 33 gives a prescription for one element. In reality, such matrix
must be computed for each element and put together to form the global system. Basically, during the
global assembly is done as follows: for a node m where m denotes a global numbering in the FE mesh
of the object domain, the contributions of all the elements attached to the node m are added. This
is usually done directly during the computation of the local matrix Ke which is not stored directly,
but its elements are added to the global matrix according to the local-to-global index mapping: while
locally, the node is represented by the pair (e, n) where e is the element index and n is the local node
index in the element (e.g. for P1 n ∈ {1, 2, 3, 4}), it has also the global index ˆn = 1 . . . N where N is
the number of nodes. The mapping is then given as (n, e) → ˆn.
In 3D case employing the P1 elements, the local matrix Ke
is a 12×12 matrix composed of 3×3 blocks
representing the stiffness relations between the nodes of the element e. During the assembly, the block
placed at the position (n, m) in the local matrix is added to the global matrix on position (ˆn, ˆm) where
(n, e) → ˆn and (m, e) → ˆm.
• The essential boundary conditions such as the homogeneous Dirichlet condition must be applied using
one of the methods introduced in the lecture (penalisation, projection or Lagrange multipliers). This
step is performed on the global matrix.
It must be emphasized that linear formulation can be used only for small strains (up to 10%), otherwise,
artifacts can be observed such as increase in volume shown in Fig. 2.
Recall that the initial equations given by Eq. 10,12, 11 contain two types of linearities:
• material linearity: the material relation between the stress and strain is linear (also called StVenantKirchhoff
formulation),
7
Figure 2: Linear elasticity showing important artifacts in large deformations.
• geometrical linearity: the full Green strain tensor is linearized, i.e. in the strain tensor definition
eij =
1
2
(ui,j + uj,i + ui,juj,i) (35)
the last quadratic term is neglected. Basically, this linearlization results in the large-deformation
artifacts.
The employment of the full Green tensor with e.g. StVenant-Kirchhoff (or another such as Mooney-Rivlin)
material would result in non-linear hyperelastic formulation. However, the resulting equations would be
much more complicated, as the resulting formulation (based on the weak form and the linearization) would
be non-linear and some iterative solution method would be necessary (e.g. Newton-Raphson). Nevertheless,
in this case, Jacobian (a derivative w.r.t. spatial coordinates) must be computed using so called Fr´echet
derivatives.
An alternative approach is represented by a corotational FEM (see for instance in Nesme, M., Payan, Y.,
Faure, F.: Efficient, physically plausible finite elements. In Dingliana, J., Ganovelli, F., eds.: Eurographics
2005, Short papers, August, 2005, Trinity College, Dublin, Irlande (2005)) or C.Felippa: A systematic
approach to the element-independent corotational dynamics of finite elements, 2000.
8