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:
Provides the constitutive response (and optionally the constitutive tangent) of several materials.
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.)
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.
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
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.
We begin our discussion by considering a static, solid mechanics, problem. Loosely speaking the the goal is to find a deformation map, , that maps a body to a deformed state that satisfies equilibrium and the boundary conditions applied on .
As is the case in the illustration, the body can be subjected to two kinds of boundary conditions:
- Essential or Dirichlet boundary conditions on , whereby the displacements are prescribed.
- Natural or Neumann boundary conditions on , whereby the tractions are prescribed. (Whereby traction-free is also perfectly acceptable.)
In practice, we are not explicitly looking for the map , but for the deformation gradient , or in fact for a displacement field . To make things a bit more explicit, the deformation gradient is defined as follows:
We start from the linear momentum balance:
where is the Cauchy stress which depends on the new position and thus on the displacement . 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
The following nomenclature has been used
The crux of the Finite Element Method is that this non-linear differential equation is solved in a weak sense. I.e.
where are test functions. For reasons that become obvious below, we apply integration by parts, which results in
Use has been made of the following chain rule
together with the symmetry of the Cauchy stress
and the following nomenclature:
The right-hand-side of this equation can be reduced to an area integral by employing Gauss’ divergence theorem. The result reads
Gauss’ divergence theorem states that
where is the normal along the surface .
The problem is now discretized using nodes that are connected through elements, which define the discretized domain . Shape functions are used to extrapolate the nodal quantities throughout the domain (and ), as follows:
Following standard Galerkin
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)
Applied to the balance equation we obtain
from which the dependency on can be dropped:
This corresponds to (non-linear) set of nodal balance equations:
- Internal forces
which is zero in the interior of the domain, i.e. in , while they can be zero or non-zero in depending on the problem details.
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 , and thus that . Also we define a strain tensor
and use some non-linear relationship between it and the stress
To simplify our discussion we assume the boundary tractions to be some known constant. Our nodal equilibrium equations now read
To come to an iterative solution, we linearize as this point. This results in
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 , can often be omitted as many constitutive tangent operators already symmetrize.
In a shorter notation, this is our iterative update:
This is a good point to study some examples:
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.
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.
Newton-Raphson in one dimension
We try to find such that
We will make a guess for and (hopefully) iteratively improve this guess. This iterative value is denoted using . Therefore we will make use of the following Taylor expansion
We now determine by neglecting higher order terms, which results in
From which we obtain as
Thereafter we set
And check if we are have reached our solution within a certain accuracy :
If not, we repeat the above steps until we do.
The iterative scheme is well understood from the following illustration:
We continue with our balance equation and add inertia an damping to it:
where is the density and the viscosity (a.k.a. the damping coefficient). The first and second time derivative of the position are respectively the velocity and the acceleration .
We can generalize this as follows (which will also simplify our proceedings below)
To retrieve the original form
But, we can now use also other expressions. For example the damping equivalent of linear elasticity:
where is the bulk viscosity while is the shear viscosity. Furthermore
Our original form is retrieved when , both are independent of , and possesses the necessary symmetries.
Like before, we will solve this equation in a weak sense
Integration by parts results in
Which we will discretize as before:
Which is independent of the test functions, hence:
Which we can denote as follows
whereby we have introduced:
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.
This results in the following expression for the mass matrix:
Here we will discuss several common time discretization steps. To simplify notation we will denote the velocity and the acceleration .
Most time integration schemes result is some form like
where contains the boundary tractions and internal forces, including their damping equivalents. The subscript indicates that the variable is a known quantity, while 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 . 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¶
Compute the position at :
Estimate the velocity at :
In the Finite Element Method a geometry is discretized using nodes. The nodes are grouped in elements which define the domain . The crux of the method is that nodal quantities, for example , are extrapolated throughout the discretized domain using shape functions . Each shape function is globally supported, however in such a way that only in the elements containing node . It is furthermore imposed that , i.e. it is one in the node , 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).
From this it becomes obvious that is polynomial through each of the nodes, and that is discontinuous across element boundaries. Note once more that each of the shape functions is globally supported, but zero outside the elements that contain the node . For node 2, the shape function is thus:
As we can see, node 2 is only non-zero in elements 1 and 2, while it is zero everywhere else. To evaluate we therefore only have to integrate on these elements (using Isoparametric transformation and quadrature):
By now it should be clear that the above allows us assemble element-by-element. For this example, graphically this corresponds to the following sum:
where the indices show that the shape functions are evaluated compared to some generic element definition (see 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 onto a generic isoparametric element of constant volume . By using this mapping it is easy to perform numerical quadrature while even reusing an existing implementation (for example the one of GooseFEM).
The mapping between the generic domain and the physical domain is as follows
where the column contains the real position vectors of the element nodes. In order to perform the quadrature on we must map also the gradient operator:
Using the above:
We can now determine the mapping between the real and the master volume:
For example for the internal force this implies
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:
To obtain , , and , simply replace with 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 or as input.
The details depend on the element type. Several standard elements types are implemented in GooseFEM.