Finite Element Method

In the sequel the theory of the Finite Element is discussed in a compact way. This discussion is by no means comprehensive. Therefore one is invited to dive in more complete textbooks. The key contribution of this reader is that it is supported by many examples that can be easily extended and customized into efficient, production ready code. To this end, the examples are in C++14, and are specifically written such that they are mostly ‘what you see is what you get’. The entire structure is in the main-file, and not hidden somewhere in a library. To simplify your life we do use several libraries, each of which however only with a dedicated task, which can be understood, used, and checked independently of the Finite Element Method or any specific application. More specifically we use:

  • GooseMaterial

    Provides the constitutive response (and optionally the constitutive tangent) of several materials.

  • cppmat

    Provides tensor classes and operations. (The amount of tensor-operations is limited in the main program, and even non-standard, but this library is crucial to compute the material response implemented in GooseMaterial.)

  • Eigen3

    A linear algebra library. As you will notice, Eigen plays an important role in GooseFEM, and glues everything together since in the end the Finite Element Method is just a way to cast a problem into a set linear or linearized equations. Most of the efficiency of the final program will depend on the efficiency of the implementation of the linear algebra. In several examples we will simplify the structure by using dense matrices together with a simple solver which solves the resulting linear system. In reality one should always use sparse matrices combined with a more efficient solver. As you will notice, many examples need only few changes to be transformed in a production code.

Note

Compilation

Unless otherwise mentioned, the examples can be compiled as follows. Provided that pkg-config is set-up correctly one can use

clang++ `pkg-config --cflags Eigen3 cppmat GooseMaterial GooseFEM` -std=c++14 -o example example_name.cpp

(whereby clang++ can be replaced by for example g++). If one does not want to use pkg-config, one has to specify -I/path/to/library for each of the libraries.

For further development it is strongly advised to include the options -Wpedantic -Wall to get on top of mistakes. Once the code is ready, one should compile with optimization (-O3) and without assertions (-DNDEBUG). The Eigen3 documentation further recommends the option -march=native to enable vectorization optimized for you architecture.

Statics

The conceptual idea

We begin our discussion by considering a static, solid mechanics, problem. Loosely speaking the the goal is to find a deformation map, \vec{x} = \varphi(\vec{X},t), that maps a body \Omega_0 to a deformed state \Omega that satisfies equilibrium and the boundary conditions applied on \Gamma.

../_images/problem1.svg

As is the case in the illustration, the body can be subjected to two kinds of boundary conditions:

  • Essential or Dirichlet boundary conditions on \Gamma_p, whereby the displacements are prescribed.
  • Natural or Neumann boundary conditions on \Gamma_u, whereby the tractions are prescribed. (Whereby traction-free is also perfectly acceptable.)

In practice, we are not explicitly looking for the map \vec{x} = \varphi(\vec{X},t), but for the deformation gradient \bm{F}, or in fact for a displacement field \vec{u} = \vec{x} - \vec{X}. To make things a bit more explicit, the deformation gradient is defined as follows:

\vec{x} = \bm{F} \cdot \vec{X}

hence

\bm{F}
=
\frac{\partial \varphi}{\partial \vec{X}}
=
\big( \vec{\nabla}_0 \, \vec{x} \big)^T
=
\bm{I} + \big( \vec{\nabla}_0 \, \vec{u} \big)^T

Momentum balance

We start from the linear momentum balance:

\vec{\nabla} \cdot \bm{\sigma}(\vec{x}) = \vec{0}
\qquad
\vec{x} \in \Omega

where \bm{\sigma} is the Cauchy stress which depends on the new position \vec{x} and thus on the displacement \vec{u}. It has been assumed that all actions are instantaneous (no inertia) and, for simplicity, that there are no body forces. Loosely speaking the interpretation of this equation is that the sum of all forces vanishes everywhere in the domain \Omega

Note

The following nomenclature has been used

\vec{\nabla} \cdot \bm{\sigma} = \frac{ \partial \sigma_{ij} }{ \partial x_i }

The crux of the Finite Element Method is that this non-linear differential equation is solved in a weak sense. I.e.

\int\limits_\Omega
  \vec{\phi}(\vec{X}) \cdot \big[\, \vec{\nabla} \cdot \bm{\sigma}(\vec{x}) \,\big] \;
\mathrm{d}\Omega
=
0
\qquad
\forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d

where \vec{\phi} are test functions. For reasons that become obvious below, we apply integration by parts, which results in

\int\limits_\Omega
  \big[\, \vec{\nabla} \vec{\phi}(\vec{X}) \,\big] : \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\int\limits_\Omega
  \vec{\nabla}
  \cdot
  \big[\, \vec{\phi}(\vec{X}) \cdot \bm{\sigma}(\vec{x}) \,\big] \;
\mathrm{d}\Omega
\qquad
\forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d

Note

Use has been made of the following chain rule

\vec{\nabla} \cdot \big[\, \vec{\phi} \cdot \bm{\sigma}^T \,\big]
=
\big[\, \vec{\nabla} \vec{\phi} \,\big] : \bm{\sigma}^T
+
\vec{\phi} \cdot \big[\, \vec{\nabla} \cdot \bm{\sigma} \,\big]

together with the symmetry of the Cauchy stress

\bm{\sigma} = \bm{\sigma}^T

and the following nomenclature:

C = \bm{A} : \bm{B} = A_{ij} B_{ji}

The right-hand-side of this equation can be reduced to an area integral by employing Gauss’ divergence theorem. The result reads

\int\limits_\Omega
  \big[\, \vec{\nabla} \vec{\phi}(\vec{X}) \,\big] : \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\int\limits_\Gamma
  \vec{\phi}(\vec{X}) \cdot
  \underbrace{
    \vec{n}(\vec{x}) \cdot \bm{\sigma}(\vec{x})
  }_{
    \vec{t}(\vec{x})
  } \;
\mathrm{d}\Gamma
\qquad
\forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d

Note

Gauss’ divergence theorem states that

\int\limits_\Omega \vec{\nabla} \cdot \vec{a}(\vec{x}) \; \mathrm{d}\Omega
=
\int\limits_\Gamma \vec{n}(\vec{x}) \cdot \vec{a}(\vec{x}) \; \mathrm{d}\Gamma

where \vec{n} is the normal along the surface \Gamma.

Discretization

The problem is now discretized using n nodes that are connected through elements, which define the discretized domain \Omega^h_0. Shape functions N_i(\vec{X}) are used to extrapolate the nodal quantities throughout the domain \Omega^h_0 (and \Omega^h), as follows:

\vec{x}(\vec{X},t)
\approx
\vec{x}^h(\vec{X},t)
=
\sum_{i=1}^{n} N_i (\vec{X}) \; \vec{x}_i (t)
=
\underline{N}^\mathsf{T} (\vec{X}) \; \underline{\vec{x}} (t)

Following standard Galerkin

\vec{\phi}(\vec{X})
\approx
\vec{\phi}^h(\vec{X})
=
\underline{N}^\mathsf{T} (\vec{X}) \; \underline{\vec{\phi}}

Note

Applied to our problem sketch, a discretization might look like this. The nodes are clearly marked as circles. The lines connecting the nodes clearly mark the elements which are in this case three-node triangles (Tri3 in GooseFEM)

../_images/problem-discretized.svg

Applied to the balance equation we obtain

\underline{\vec{\phi}}^\mathsf{T} \cdot
\int\limits_{\Omega^h}
  \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big]
  \cdot
  \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\underline{\vec{\phi}}^\mathsf{T} \cdot
\int\limits_{\Gamma^h}
  \underline{N}(\vec{X}) \cdot
  \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\qquad
\forall \; \underline{\vec{\phi}} \in \mathbb{R}^d_n

from which the dependency on \underline{\vec{\phi}} can be dropped:

\int\limits_{\Omega^h}
  \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big]
  \cdot
  \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\int\limits_{\Gamma^h}
  \underline{N}(\vec{X}) \cdot
  \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma

This corresponds to (non-linear) set of nodal balance equations:

\underline{\vec{f}}(\vec{x})
=
\underline{\vec{t}}(\vec{x})

with:

  • Internal forces

\underline{\vec{f}}(\vec{x})
=
\int\limits_{\Omega^h}
  \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big]
  \cdot
  \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega

  • Boundary tractions

    \underline{\vec{t}}(\vec{x})
=
\int\limits_{\Gamma^h}
  \underline{N}(\vec{X}) \cdot
  \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma

    which is zero in the interior of the domain, i.e. in \Omega^h \bigcap \Gamma^h, while they can be zero or non-zero in \Gamma^h depending on the problem details.

Iterative solution – small strain

A commonly used strategy to solve the non-linear system, is the iterative Newton-Raphson scheme (see inset below). The idea is thereby to formulate an initial guess for the solution, determine possible residual forces, and use these forces to come to a better guess for the solution. This is continued until the solution has been found, i.e. when the residual vanishes.

This solution technique is discussed here in the context of small deformations, while it is later generalized. Assuming the deformations to be small allows us to assume that \Omega = \Omega_0, and thus that \nabla = \nabla_0. Also we define a strain tensor

\bm{\varepsilon}
=
\tfrac{1}{2} \left[ \nabla_0 \vec{u} + \big[\, \nabla_0 \vec{u} \,\big]^T \right]
=
\mathbb{I}_s : \big[\, \nabla_0 \vec{u} \,\big]

and use some non-linear relationship between it and the stress

\bm{\sigma} = \bm{\sigma} \big( \bm{\varepsilon} \big)

To simplify our discussion we assume the boundary tractions to be some known constant. Our nodal equilibrium equations now read

\underline{\vec{r}}(\vec{x})
=
\underline{\vec{t}}
-
\underline{\vec{f}}(\vec{x})
=
\underline{\vec{0}}

with

\underline{\vec{f}}(\vec{x})
=
\int\limits_{\Omega^h_0}
  \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big]
  \cdot
  \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega

To come to an iterative solution, we linearize as this point. This results in

\int\limits_{\Omega^h_0}
  \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big]
  \cdot
  \mathbb{K}\big(\vec{x}_{(i)}\big)
  \cdot
  \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big]^\mathsf{T} \;
\mathrm{d}\Omega
\cdot \delta \underline{\vec{x}}
=
\underline{\vec{t}}
-
\int\limits_{\Omega^h_0}
  \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big]
  \cdot
  \bm{\sigma}\big(\vec{x}_{(i)}\big) \;
\mathrm{d}\Omega

where

\mathbb{K}\big(\vec{x}_{(i)}\big)
=
\left. \frac{\partial \bm{\sigma}}{\partial \bm{\varepsilon}} \right|_{\vec{x}_{(i)}}
:
\mathbb{I}_s

where the left part is the constitutive tangent operator and the right part comes from the strain definition. Note that this right part, the symmetrization using \mathbb{I}_s, can often be omitted as many constitutive tangent operators already symmetrize.

In a shorter notation, this is our iterative update:

\underline{\underline{\mathbb{K}}}_{(i)} \cdot \delta \underline{\vec{x}}
=
\underline{\vec{t}}
-
\underline{\vec{f}}_{(i)}

with

\underline{\underline{\mathbb{K}}}_{(i)}
=
\int\limits_{\Omega^h_0}
  \big[\, \vec{\nabla}_0 \underline{N} \,\big]
  \cdot
  \mathbb{K}\big(\vec{x}_{(i)}\big)
  \cdot
  \big[\, \vec{\nabla}_0 \underline{N} \,\big]^\mathsf{T} \;
\mathrm{d}\Omega

and

\underline{\vec{f}}_{(i)}
=
\int\limits_{\Omega^h_0}
  \big[\, \vec{\nabla}_0 \underline{N} \,\big]
  \cdot
  \bm{\sigma}\big(\vec{x}_{(i)}\big) \;
\mathrm{d}\Omega

Note

This is a good point to study some examples:

  • Linear statics – small strain

    We slowly work up to an iterative scheme starting from a linear problem, written however in such a way that the step towards a non-linear problem is small.

  • Non-linear statics – small strain

    Here we employ Newton-Raphson to solve the non-linear equilibrium equation. It is easy to see that once the above examples have been understood this step is indeed trivial.

Note

Newton-Raphson in one dimension

We try to find x such that

r(x) = 0

We will make a guess for x and (hopefully) iteratively improve this guess. This iterative value is denoted using x_{(i)}. Therefore we will make use of the following Taylor expansion

r \big( x_{(i+1)} \big)
=
r \big( x_{(i)} \big)
+
\left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \delta x
+
\mathcal{O} \big( \delta x^2 \big)
\approx
0

where

\delta x = x_{(i+1)} - x_{(i)}

We now determine \delta x by neglecting higher order terms, which results in

r \big( x_{(i)} \big)
+
\left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \delta x
=
0

From which we obtain \delta x as

\delta x
=
- \left[ \left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \right]^{-1}
r \big( x_{(i)} \big)

Thereafter we set

x_{(i+1)} = x_{(i)} + \delta x

And check if we are have reached our solution within a certain accuracy \epsilon:

\left| r \big( x_{(i+1)} \big) \right| < \epsilon

If not, we repeat the above steps until we do.

The iterative scheme is well understood from the following illustration:

../_images/newton-raphson.svg

Dynamics

Momentum balance

We continue with our balance equation and add inertia an damping to it:

\rho\, \ddot{\vec{x}}
=
\vec{\nabla} \cdot
\bm{\sigma}(\vec{x})
+
\eta\, \nabla^2\dot{\vec{x}}
\qquad
\vec{x} \in \Omega

where \rho is the density and \eta the viscosity (a.k.a. the damping coefficient). The first and second time derivative of the position \vec{x} are respectively the velocity \vec{v} = \dot{\vec{x}} and the acceleration \vec{a} = \ddot{\vec{x}}.

We can generalize this as follows (which will also simplify our proceedings below)

\rho(\vec{x})\, \ddot{\vec{x}}
=
\vec{\nabla} \cdot
\big[\, \bm{\sigma}(\vec{x}) + \bm{\sigma}_{\eta}(\vec{\dot{x}} ) \,\big]
\qquad
\vec{x} \in \Omega

Note

To retrieve the original form

\bm{\sigma}_{\eta} = \eta\; \vec{\nabla} \dot{\vec{x}}

But, we can now use also other expressions. For example the damping equivalent of linear elasticity:

\bm{\sigma}_{\eta} (\vec{x}) = \mathbb{C}_{\eta} (\vec{x}) : \dot{\bm{\varepsilon}} (\vec{x})

with

\mathbb{C}_{\eta} (\vec{x})
=
\kappa (\vec{x}) \bm{I} \otimes \bm{I}
+
2 \gamma (\vec{x}) \mathbb{I}_d

where \kappa is the bulk viscosity while \gamma is the shear viscosity. Furthermore

\dot{\bm{\varepsilon}} (\vec{x})
=
\tfrac{1}{2} \big[\, \vec{\nabla} \dot{\vec{x}} + [\, \vec{\nabla} \dot{\vec{x}} \,]^T \,\big]

Our original form is retrieved when \kappa = \tfrac{2}{3} \gamma, both are independent of \vec{x}, and \dot{\vec{x}} possesses the necessary symmetries.

Like before, we will solve this equation in a weak sense

\int\limits_\Omega
  \rho(\vec{x})\; \vec{\phi}(\vec{X}) \cdot \ddot{\vec{x}} \;
\mathrm{d}\Omega
=
\int\limits_\Omega
  \vec{\phi}(\vec{X})
  \cdot
  \Big[\,
    \vec{\nabla}
    \cdot
    \big[\, \bm{\sigma}(\vec{x}) + \bm{\sigma}_{\eta}(\vec{\dot{x}} ) \,\big]
  \,\Big] \;
\mathrm{d}\Omega
\qquad
\forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d

Integration by parts results in

\int\limits_\Omega
  \rho(\vec{x})\; \vec{\phi}(\vec{X}) \cdot \ddot{\vec{x}} \;
\mathrm{d}\Omega
=
\int\limits_\Gamma
  \vec{\phi}(\vec{X}) \cdot \big[\, \vec{t}(\vec{x}) + \vec{t}_{\eta}(\vec{x}) \,\big] \;
\mathrm{d}\Gamma
-
\int\limits_\Omega
  \big[\, \vec{\nabla} \vec{\phi}(\vec{X}) \,\big]
  :
  \big[\, \bm{\sigma}(\vec{x}) + \bm{\sigma}_{\eta}(\dot{\vec{x}}) \,\big] \;
\mathrm{d}\Omega
\qquad
\forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d

Which we will discretize as before:

\underline{\vec{\phi}}^\mathsf{T} \cdot
\int\limits_\Omega
  \rho(\vec{x})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \;
\mathrm{d}\Omega \;
\underline{\ddot{\vec{x}}}
=
\underline{\vec{\phi}}^\mathsf{T} \cdot
\int\limits_\Gamma
  \underline{N}(\vec{X})\; \big[\, \vec{t}(\vec{x}) + \vec{t}_{\eta}(\vec{x}) \,\big] \;
\mathrm{d}\Gamma
-
\underline{\vec{\phi}}^\mathsf{T} \cdot
\int\limits_\Omega
  \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big]
  :
  \big[\, \bm{\sigma}(\vec{x}) + \bm{\sigma}_{\eta}(\dot{\vec{x}}) \,\big] \;
\mathrm{d}\Omega
\qquad
\forall \; \underline{\vec{\phi}} \in \mathbb{R}^d_n

Which is independent of the test functions, hence:

\int\limits_\Omega
  \rho(\vec{x})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \;
\mathrm{d}\Omega \;
\underline{\ddot{\vec{x}}}
=
\int\limits_\Gamma
  \underline{N}(\vec{X})\; \big[\, \vec{t}(\vec{x}) + \vec{t}_{\eta}(\vec{x}) \,\big] \;
\mathrm{d}\Gamma
-
\int\limits_\Omega
  \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big]
  :
  \big[\, \bm{\sigma}(\vec{x}) + \bm{\sigma}_{\eta}(\dot{\vec{x}}) \,\big] \;
\mathrm{d}\Omega

Which we can denote as follows

\underline{\underline{M}}(\vec{x})\; \underline{\ddot{\vec{x}}}
=
\underline{\vec{t}}(\vec{x})
+
\underline{\vec{t}}_{\eta}(\vec{x})
-
\underline{\vec{f}}(\vec{x})
-
\underline{\vec{f}}_{\eta}(\vec{x})

whereby we have introduced:

  • Mass matrix

    \underline{\underline{M}}(\vec{x})
=
\int\limits_\Omega
  \rho(\vec{x})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \;
\mathrm{d}\Omega

  • Boundary tractions

    \underline{\vec{t}}(\vec{x})
=
\int\limits_\Gamma
  \underline{N}(\vec{X})\; \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\qquad
\mathrm{and}
\qquad
\underline{\vec{t}}_{\eta}(\vec{x})
=
\int\limits_\Gamma
  \underline{N}(\vec{X})\; \vec{t}_{\eta}(\vec{x}) \;
\mathrm{d}\Gamma

  • Internal forces

    \underline{\vec{f}}(\vec{x})
=
\int\limits_\Omega
  \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big] : \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
\qquad
\mathrm{and}
\qquad
\underline{\vec{f}}(\vec{x})
=
\int\limits_\Omega
  \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big] : \bm{\sigma}_{\eta}(\dot{\vec{x}}) \;
\mathrm{d}\Omega

Note

In many problems it make sense to assume the mass matrix constant, as any change of volume results in an equivalent change of the density, i.e.

\int\limits_{\Omega}
  \rho(\vec{x})
\;\mathrm{d}\Omega
=
\int\limits_{\Omega_0}
  \rho(\vec{X})
\;\mathrm{d}\Omega_0

This results in the following expression for the mass matrix:

\underline{\underline{M}}(\vec{X})
=
\int\limits_{\Omega_0}
  \rho(\vec{X})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \;
\mathrm{d}\Omega_0
=
\mathrm{constant}

Time discretization

Here we will discuss several common time discretization steps. To simplify notation we will denote the velocity \vec{v} = \dot{\vec{x}} and the acceleration \vec{a} = \ddot{\vec{x}}.

Note

Most time integration schemes result is some form like

\underline{\underline{M}}\; \underline{\vec{a}}_{n+1}
=
\underline{\vec{q}}_{n}

where \underline{\vec{q}}_{n} contains the boundary tractions and internal forces, including their damping equivalents. The subscript n indicates that the variable is a known quantity, while n+1 indicates that it is an unknown quantity. To enhance computational efficiency, it may be a good option to approximate the mass matrix in such a way that it becomes diagonal. Consequently, no system has be solved to find \underline{\vec{a}}_{n+1}. One only has to invert an array of scalars. Since in addition the mass matrix is almost often assumed constant, this factorization has to be performed only once for the entire simulation.

Physically one can interpret this assumption as assuming the damping to be concentrated on the nodes.

See: Diagonal mass matrix.

Velocity Verlet with damping

  1. Compute the position at t_{n+1} = t_{n} + \Delta_t:

    \vec{x}_{n+1}
=
\vec{x}_{n} + \Delta_t \vec{v}_{n} + \tfrac{1}{2} \Delta_t^2 \vec{a}_{n}

  2. Estimate the velocity at t_{n+1} = t_{n} + \Delta_t:

\hat{\vec{v}}_{n+1}
=
\vec{v}_{n}
+
\tfrac{1}{2} \Delta_t \Big[\,
  \vec{a}_{n} + \vec{a} ( \vec{x}_{n+1} , \vec{v}_{n} + \Delta_t \vec{a}_{n} , t_{n+1} ) \,
\Big]

  1. Correct \hat{\vec{v}}_{n+1}:

    \vec{v}_{n+1}
=
\vec{v}_{n}
+
\tfrac{1}{2} \Delta_t \Big[\,
  \vec{a}_{n} + \vec{a} ( \vec{x}_{n+1} , \hat{\vec{v}}_{n+1} , t_{n+1} ) \,
\Big]

Shape functions

In the Finite Element Method a geometry is discretized using nodes. The nodes are grouped in elements which define the domain \Omega^h_0. The crux of the method is that nodal quantities, for example \vec{u}_i, are extrapolated throughout the discretized domain \Omega^h_0 using shape functions N_i (\vec{X}). Each shape function is globally supported, however in such a way that N_i (\vec{X}) \neq 0 only in the elements containing node i. It is furthermore imposed that N_i (\vec{X}_j) = \delta_{ij}, i.e. it is one in the node i, and zero in all other nodes.

For a one-dimensional problem comprising four linear elements and five nodes the shape functions are sketched below (whereby the node numbers are in color, while the element numbers are in black, in between the nodes).

../_images/shape-functions-1d.svg

From this it becomes obvious that N_i (\vec{X}) is polynomial through each of the nodes, and that \partial N_i / \partial \vec{X} is discontinuous across element boundaries. Note once more that each of the shape functions N_i (X) is globally supported, but zero outside the elements that contain the node i. For node 2, the shape function is thus:

../_images/shape-functions-1d-node-2.svg

As we can see, node 2 is only non-zero in elements 1 and 2, while it is zero everywhere else. To evaluate \vec{f}_2 we therefore only have to integrate on these elements (using Isoparametric transformation and quadrature):

\vec{f}_2
=
\int\limits_{\Omega^1}
  \big[\, \vec{\nabla} N^1_2(\vec{X}) \,\big]
  \cdot
  \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
+
\int\limits_{\Omega^2}
  \big[\, \vec{\nabla} N^2_2(\vec{X}) \,\big]
  \cdot
  \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega

By now it should be clear that the above allows us assemble \underline{f} element-by-element. For this example, graphically this corresponds to the following sum:

../_images/shape-functions-1d-element-0.svg
../_images/shape-functions-1d-element-1.svg
../_images/shape-functions-1d-element-2.svg
../_images/shape-functions-1d-element-3.svg

where the indices show that the shape functions are evaluated compared to some generic element definition (see Isoparametric transformation and quadrature).

Isoparametric transformation and quadrature

A very important concept in the Finite Element Method is the isoparametric transformation. It allows us to map an arbitrarily shaped element with volume \Omega^e onto a generic isoparametric element of constant volume Q. By using this mapping it is easy to perform numerical quadrature while even reusing an existing implementation (for example the one of GooseFEM).

../_images/isoparametric-transform.svg

The mapping between the generic domain Q and the physical domain \Omega^e is as follows

\vec{x} ( \vec{\xi} ) = \big[\, \underline{N}^{e} \,\big]^\mathsf{T} \underline{x}^e

where the column \underline{x}^e contains the real position vectors of the element nodes. In order to perform the quadrature on Q we must map also the gradient operator:

\vec{\nabla}_{\xi}\,
=
\vec{e}_i \frac{\partial}{\partial \xi_i}
=
\vec{e}_i \frac{\partial x_j(\vec{\xi})}{\partial \xi_i} \frac{\partial}{\partial x_j}
=
\vec{e}_i \frac{\partial x_j(\vec{\xi})}{\partial \xi_i} \vec{e}_j \cdot \vec{e}_k \frac{\partial}{\partial x_k}
=
\big[\, \vec{\nabla}_{\xi}\, \vec{x}(\vec{\xi}) \,\big] \cdot \vec{\nabla}
=
\bm{J}(\vec{\xi}) \cdot \vec{\nabla}

or

\vec{\nabla} = \bm{J}^{-1}(\vec{\xi}) \cdot \vec{\nabla}_{\xi}\,

with

\bm{J}(\vec{\xi})
=
\vec{\nabla}_{\xi}\, \vec{x}(\vec{\xi})
=
\big[\, \vec{\nabla}_{\xi}\, \underline{N}^{e} \,\big]^\mathsf{T} \; \underline{x}^e

Using the above:

\vec{\nabla} \underline{N}^{e}
=
\bm{J}^{-1}(\vec{\xi}) \cdot  \big[\, \vec{\nabla}_{\xi}\, \underline{N}^{e} \,\big]

We can now determine the mapping between the real and the master volume:

\mathrm{d} \Omega
=
\mathrm{d} \vec{x}_0 \times \mathrm{d} \vec{x}_1 \cdot \mathrm{d} \vec{x}_2
=
\left[ \mathrm{d} \vec{x}_0 \cdot \bm{J}(\vec{\xi}) \right] \times
\left[ \mathrm{d} \vec{x}_1 \cdot \bm{J}(\vec{\xi}) \right] \cdot
\left[ \mathrm{d} \vec{x}_2 \cdot \bm{J}(\vec{\xi}) \right]
=
\det \big( \bm{J}(\vec{\xi}) \big)\,
\mathrm{d} \vec{\xi}_0 \times \mathrm{d} \vec{\xi}_1 \cdot \mathrm{d} \vec{\xi}_2
=
\det \big( \bm{J}(\vec{\xi}) \big)\, \mathrm{d} Q

For example for the internal force this implies

\underline{\vec{f}^e}
=
\int\limits_{\Omega^e}
  \big[\, \vec{\nabla} \underline{N} \,\big]
  \cdot
  \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\int\limits_{Q}
  \big[\, \vec{\nabla} \underline{N} \,\big]
  \cdot
  \bm{\sigma}(\vec{x}) \;
  \det \big( \bm{J}(\vec{\xi}) \big) \;
\mathrm{d}Q

Numerical quadrature can be formulated (exactly) on the master element. It corresponds to taking the weighted sum of the integrand evaluated at specific quadrature points (or integration-points). Again, for our internal force:

\underline{\vec{f}^e}
=
\sum_{k}^{n_k}
w_k
\big[\, \vec{\nabla} \underline{N} \,\big]_{\vec{\xi} = \vec{\xi}_k}
\cdot
\bm{\sigma}\big(\vec{x}(\vec{\xi}_k)\big) \;
\det \big( \bm{J}(\vec{\xi}_k) \big) \;

Note

To obtain \vec{X}(\vec{\xi}), \vec{\nabla}_0, and \int\limits_{\Omega_0} . \;\mathrm{d}\Omega, simply replace \underline{x}^e with \underline{X}^e in the first equation. For this reason the same element implementation (of for example GooseFEM) can be used in small strain and finite strain (total Lagrange and updated Lagrange), proving either \underline{X}^e or \underline{X}^e + \underline{u}^e as input.

Note

The details depend on the element type. Several standard elements types are implemented in GooseFEM.