GooseFEM¶
Introduction¶
GooseFEM provides several standard finite element simulations, some common element definitions, and some simple finite element meshes. A basic documentation is provided here, whereby one is highly urged to look into the code itself, at minimum the function decelerators and their comments.
Note
This library is free to use under the GPLv3 license. Any additions are very much appreciated, in terms of suggested functionality, code, documentation, testimonials, word of mouth advertisement, .... Bug reports or feature requests can be filed on GitHub. As always, the code comes with no guarantee. None of the developers can be held responsible for possible mistakes.
Download: .zip file  .tar.gz file.
(c  GPLv3) T.W.J. de Geus (Tom)  tom@geus.me  www.geus.me  github.com/tdegeus/GooseFEM
Contents¶
GooseFEM documentation¶
Compiling¶
Introduction¶
This module is header only. So one just has to #include <GooseFEM/GooseFEM.h>
. somewhere in the source code, and to tell the compiler where the headerfiles are. For the latter, several ways are described below.
Before proceeding, a words about optimization. Of course one should use optimization when compiling the release of the code (O2
or O3
). But it is also a good idea to switch off the assertions in the code (mostly checks on size) that facilitate easy debugging, but do cost time. Therefore, include the flag DNDEBUG
. Note that this is all C++ standard. I.e. it should be no surprise, and it always a good idea to do.
Manual compiler flags¶
GNU / Clang¶
Add the following compiler’s arguments:
I${PATH_TO_GOOSEFEM}/src std=c++14
Note
(Not recommended)
If you want to avoid separately including the header files using a compiler flag, git submodule
is a nice way to go:
Include this module as a submodule using
git submodule add https://github.com/tdegeus/GooseFEM.git
.Replace the first line of this example by
#include "GooseFEM/src/GooseFEM/GooseFEM.h"
.If you decide to manually copy the header file, you might need to modify this relative path to your liking.
Or see (Semi)Automatic compiler flags. You can also combine the git submodule
with any of the below compiling strategies.
(Semi)Automatic compiler flags¶
Install¶
To enable (semi)automatic build, one should ‘install’ GooseFEM
somewhere.
Proceed to a (temporary) build directory. For example
$ cd /path/to/GooseFEM/src/build
‘Build’
GooseFEM
$ cmake .. $ make install
(If you’ve used another build directory, change the first command to
$ cmake /path/to/GooseFEM/src
)
Proceed to a (temporary) build directory. For example
$ cd /path/to/GooseFEM/src/build
‘Build’
GooseFEM
to install it in a custom location$ mkdir /custom/install/path $ cmake .. DCMAKE_INSTALL_PREFIX:PATH=/custom/install/path $ make install
(If you’ve used another build directory, change the first command to
$ cmake /path/to/GooseFEM/src
)Add the following path to your
~/.bashrc
(or~/.zshrc
):export PKG_CONFIG_PATH=/custom/install/path/share/pkgconfig:$PKG_CONFIG_PATH
Note
(Not recommended)
If you do not wish to use CMake
for the installation, or you want to do something custom. You can of course. Follow these steps:
 Copy the file
src/GooseFEM.pc.in
toGooseFEM.pc
to some location that can be found bypkg_config
(for example by addingexport PKG_CONFIG_PATH=/path/to/GooseFEM.pc:$PKG_CONFIG_PATH
to the.bashrc
).  Modify the line
prefix=@CMAKE_INSTALL_PREFIX@
toprefix=/path/to/GooseFEM
.  Modify the line
Cflags: I${prefix}/@INCLUDE_INSTALL_DIR@
toCflags: I${prefix}/src
.  Modify the line
Version: @GOOSEFEM_VERSION_NUMBER@
to reflect the correct release version.
Compiler arguments from ‘pkgconfig’¶
Instead of I...
one can now use
`pkgconfig cflags GooseFEM` std=c++14
as compiler argument.
Compiler arguments from ‘cmake’¶
Add the following to your CMakeLists.txt
:
set(CMAKE_CXX_STANDARD 14)
find_package(PkgConfig)
pkg_check_modules(GOOSEFEM REQUIRED GooseFEM)
include_directories(${GOOSEFEM_INCLUDE_DIRS})
GooseFEM::Mesh¶
Note
Source:
#include <GooseFEM/Mesh.h>
GooseFEM::Mesh::dofs
¶
Get a sequential list of DOFnumbers for each vectorcomponent of each node.
MatS GooseFEM::Mesh::dofs ( size_t nnode , size_t ndim )
For example for 3 nodes in 2 dimensions the output is
GooseFEM::Mesh::renumber
¶
Renumber indices to lowest possible indices (returns a copy, input not modified).
MatS GooseFEM::Mesh::renumber ( const MatS &in )
For example:
is renumbered to
Reorder indices such that some items are at the beginning or the end (returns a copy, input not modified).
MatS GooseFEM::Mesh::renumber ( const MatS &in , const ColS &idx, std::string location="end" );
For example:
with
Implies that
in
is renumbered such that the 6th item becomes the onebeforelast item (), and the 4th item become the last (). The remaining items are renumbered to the lowest index while keeping the same order. The result:Consider also [
source: figures/Mesh/renumber.cpp
]
GooseFEM::Mesh::Tri3¶
Note
Source:
#include <GooseFEM/MeshTri3.h>
GooseFEM::Mesh::Quad4¶
Note
Source:
#include <GooseFEM/MeshQuad4.h>
GooseFEM::Mesh::Quad4::Regular
¶
GooseFEM::Mesh::Quad4::Regular(size_t nx, size_t ny, double h=1.);
Regular mesh of linear quadrilaterals in twodimensions. The element edges are all of the same size (by default equal to one), optional scaling can be applied afterwards. For example the mesh shown below that consists of 21 x 11 elements. In that image the element numbers are indicated with a color, and likewise for the boundary nodes.
Methods:
// A matrix with on each row a nodal coordinate:
// [ x , y ]
MatD = GooseFEM::Mesh::Quad4::Regular.coor();
// A matrix with the connectivity, with on each row to the nodes of each element
MatS = GooseFEM::Mesh::Quad4::Regular.conn();
// A list of boundary nodes
ColS = GooseFEM::Mesh::Quad4::Regular.nodesBottom();
ColS = GooseFEM::Mesh::Quad4::Regular.nodesTop();
ColS = GooseFEM::Mesh::Quad4::Regular.nodesLeft();
ColS = GooseFEM::Mesh::Quad4::Regular.nodesRight();
// A matrix with periodic node pairs on each row:
// [ independent nodes, dependent nodes ]
MatS = GooseFEM::Mesh::Quad4::Regular.nodesPeriodic();
// The node at the origin
size_t = GooseFEM::Mesh::Quad4::Regular.nodeOrigin();
// A matrix with DOFnumbers: two per node in sequential order
MatS = GooseFEM::Mesh::Quad4::Regular.dofs();
// A matrix with DOFnumbers: two per node in sequential order
// All the periodic repetitions are eliminated from the system
MatS = GooseFEM::Mesh::Quad4::Regular.dofsPeriodic();
GooseFEM::Mesh::Quad4::FineLayer
¶
Regular mesh with a fine layer of quadrilateral elements, and coarser elements above and below.
Note
The coarsening depends strongly on the desired number of elements in horizontal elements. The becomes clear from the following example:
mesh = GooseFEM::Mesh::Quad4::FineLayer(6*9 ,51); // left image : 546 elements
mesh = GooseFEM::Mesh::Quad4::FineLayer(6*9+3,51); // middle image : 703 elements
mesh = GooseFEM::Mesh::Quad4::FineLayer(6*9+1,51); // right image : 2915 elements
Methods:
// A matrix with on each row a nodal coordinate:
// [ x , y ]
MatD = GooseFEM::Mesh::Quad4::Regular.coor();
// A matrix with the connectivity, with on each row to the nodes of each element
MatS = GooseFEM::Mesh::Quad4::Regular.conn();
// A list of boundary nodes
ColS = GooseFEM::Mesh::Quad4::Regular.nodesBottom();
ColS = GooseFEM::Mesh::Quad4::Regular.nodesTop();
ColS = GooseFEM::Mesh::Quad4::Regular.nodesLeft();
ColS = GooseFEM::Mesh::Quad4::Regular.nodesRight();
// A matrix with periodic node pairs on each row:
// [ independent nodes, dependent nodes ]
MatS = GooseFEM::Mesh::Quad4::Regular.nodesPeriodic();
// The node at the origin
size_t = GooseFEM::Mesh::Quad4::Regular.nodeOrigin();
// A matrix with DOFnumbers: two per node in sequential order
MatS = GooseFEM::Mesh::Quad4::Regular.dofs();
// A matrix with DOFnumbers: two per node in sequential order
// All the periodic repetitions are eliminated from the system
MatS = GooseFEM::Mesh::Quad4::Regular.dofsPeriodic();
// A list with the element numbers of the fine elements in the center of the mesh
// (highlighted in the plot below)
ColS = GooseFEM::Mesh::Quad4::FineLayer.elementsFine();
.. image:: figures/MeshQuad4/FineLayer/example_elementsFine.svg
:width: 500px
:align: center
GooseFEM::Dynamics::Diagonal¶
Note
Source:
#include <GooseFEM/DynamicsDiagonalPeriodic.h>
#include <GooseFEM/DynamicsDiagonalSemiPeriodic.h>
#include <GooseFEM/DynamicsDiagonalLinearStrainQuad4.h>
Overview¶
Principle¶
The philosophy is to provide some structure to efficiently run a finite element simulation which remains customizable. Even more customization can be obtained by copy/pasting the source and modifying it to your need. The idea that is followed involves a hierarchy of three classes, whereby a class that is higher in hierarchy writes to some field of the class that is lower in hierarchy, runs a function, and reads from some field. In general:
Discretized system (
GooseFEM::Dynamics::Diagonal::Periodic
,GooseFEM::Dynamics::Diagonal::SemiPeriodic
).Defines the discretized system, and assembles the (inverse) mass matrix, the displacement dependent forces, and the velocity dependent forces from element arrays that are provided by the element definition.
Element definition (
GooseFEM::Dynamics::Diagonal::LinearStrain::Qaud4
)Provides the element arrays by performing numerical quadrature. At the integration point the constitutive response is probed from the quadrature point definition.
Quadrature point definition
This class is not provided, and should be provided by the user.
Example¶
A simple example is:
[source: examples/DynamicsDiagonalPeriodic/laminate/no_damping/Verlet/main.cpp
]
#include <Eigen/Eigen>
#include <cppmat/cppmat.h>
#include <GooseFEM/GooseFEM.h>
#include <GooseMaterial/AmorphousSolid/LinearStrain/Elastic/Cartesian2d.h>
// 
using MatS = GooseFEM::MatS;
using MatD = GooseFEM::MatD;
using ColD = GooseFEM::ColD;
using T2 = cppmat::cartesian2d::tensor2 <double>;
using T2s = cppmat::cartesian2d::tensor2s<double>;
namespace GM = GooseMaterial::AmorphousSolid::LinearStrain::Elastic::Cartesian2d;
// ===================================================================
class Quadrature
{
public:
T2s eps, epsdot, sig;
size_t nhard;
GM::Material hard, soft;
double Ebar, Vbar;
Quadrature(size_t nhard);
double density (size_t elem, size_t k, double V);
void stressStrain (size_t elem, size_t k, double V);
void stressStrainRate (size_t elem, size_t k, double V);
void stressStrainPost (size_t elem, size_t k, double V);
void stressStrainRatePost(size_t elem, size_t k, double V);
};
// 
Quadrature::Quadrature(size_t _nhard)
{
nhard = _nhard;
hard = GM::Material(100.,10.);
soft = GM::Material(100., 1.);
}
// 
double Quadrature::density(size_t elem, size_t k, double V)
{
return 1.0;
}
// 
void Quadrature::stressStrain(size_t elem, size_t k, double V)
{
if ( elem < nhard ) sig = hard.stress(eps);
else sig = soft.stress(eps);
}
// 
void Quadrature::stressStrainRate(size_t elem, size_t k, double V)
{
}
// 
void Quadrature::stressStrainPost(size_t elem, size_t k, double V)
{
Vbar += V;
if ( elem < nhard ) Ebar += hard.energy(eps) * V;
else Ebar += soft.energy(eps) * V;
}
// 
void Quadrature::stressStrainRatePost(size_t elem, size_t k, double V)
{
}
// ===================================================================
int main()
{
// class which provides the mesh
GooseFEM::Mesh::Quad4::Regular mesh(40,40,1.);
// class which provides the constitutive response at each quadrature point
auto quadrature = std::make_shared<Quadrature>(40*40/4);
// class which provides the response of each element
using Elem = GooseFEM::Dynamics::Diagonal::LinearStrain::Quad4<Quadrature>;
auto elem = std::make_shared<Elem>(quadrature);
// class which provides the system and an increment
GooseFEM::Dynamics::Diagonal::Periodic<Elem> sim(
elem,
mesh.coor(),
mesh.conn(),
mesh.dofsPeriodic(),
1.e2,
0.0
);
// loop over increments
for ( ... )
{
//  set displacement of fixed DOFs
...
//  compute time increment
sim.Verlet();
//  postprocess
quadrature>Ebar = 0.0;
quadrature>Vbar = 0.0;
sim.post();
...
}
return 0;
}
Pseudocode¶
What is happening inside Verlet
is evaluating the forces (and the mass matrix), and updating the displacements by solving the system. In pseudocode:
Mass matrix:
sim.computeMinv(): { for e in elements: sim>elem>xe(i,j) = ... sim>elem>ue(i,j) = ... sim>elem>computeM(e): { for k in integrationpoints: sim>elem>M(...,...) += ... * sim>elem>quad>density(e,k,V) } M(...) += sim>elem>M(i,i) }
Displacement dependent force:
sim.computeFu(): { for e in elements: sim>elem>xe(i,j) = ... sim>elem>ue(i,j) = ... sim>elem>computeFu(e): { for k in integrationpoints: sim>elem>quad>eps(i,j) = ... sim>elem>quad>stressStrain(e,k,V) sim>elem>fu(...) += ... * sim>elem>quad>sig(i,j) } Fu(...) += sim>elem>fu(i) }
Velocity dependent force:
sim.computeFv(): { for e in elements: sim>elem>xe(i,j) = ... sim>elem>ue(i,j) = ... sim>elem>ve(i,j) = ... sim>elem>computeFv(e): { for k in integrationpoints: sim>elem>quad>epsdot(i,j) = ... sim>elem>quad>stressStrainRate(e,k,V) sim>elem>fv(...) += ... * sim>elem>quad>sig(i,j) } Fv(...) += sim>elem>fu(i) }
Signature¶
From this it is clear that:
GooseFEM::Dynamics::Diagonal::Periodic
requires the following minimal signature fromGooseFEM::Dynamics::Diagonal::LinearStrain::Qaud4
:class Element { public: matrix M; // should have operator(i,j) column fu, fv; // should have operator(i) matrix xe, ue, ve; // should have operator(i,j) void computeM (size_t elem); // mass matrix < quad>density void computeFu(size_t elem); // displacement dependent forces < quad>stressStrain void computeFv(size_t elem); // displacement dependent forces < quad>stressStrainRate void post (size_t elem); // postprocess < quad>stressStrain(Rate) }
GooseFEM::Dynamics::Diagonal::LinearStrain::Qaud4
requires the minimal signature fromQuadrature
class Quadrature { public: tensor eps, epsdot, sig; // should have operator(i,j) double density (size_t elem, size_t k, double V); void stressStrain (size_t elem, size_t k, double V); void stressStrainRate (size_t elem, size_t k, double V); void stressStrainPost (size_t elem, size_t k, double V); void stressStrainRatePost(size_t elem, size_t k, double V); }
Notes for developers¶
Create a new release¶
Update the version number as follows in
src/GooseFEM/Macros.h
. The C++ and Python distributions will read from this.Upload the changes to GitHub and create a new release there (with the correct version number).
Upload the package to PyPi:
$ python3 setup.py bdist_wheel universal $ twine upload dist/*
Theory¶
Finite Element Method¶
Contents
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 mainfile, 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 tensoroperations is limited in the main program, and even nonstandard, 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.
Note
Compilation
Unless otherwise mentioned, the examples can be compiled as follows. Provided that pkgconfig
is setup correctly one can use
clang++ `pkgconfig 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 pkgconfig
, 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, , 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 tractionfree 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:
hence
Momentum balance¶
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
Note
The following nomenclature has been used
The crux of the Finite Element Method is that this nonlinear 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
Note
Use has been made of the following chain rule
together with the symmetry of the Cauchy stress
and the following nomenclature:
The righthandside of this equation can be reduced to an area integral by employing Gauss’ divergence theorem. The result reads
Note
Gauss’ divergence theorem states that
where is the normal along the surface .
Discretization¶
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
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 threenode triangles (Tri3 in GooseFEM)
Applied to the balance equation we obtain
from which the dependency on can be dropped:
This corresponds to (nonlinear) set of nodal balance equations:
with:
 Internal forces
Boundary tractions
which is zero in the interior of the domain, i.e. in , while they can be zero or nonzero in depending on the problem details.
Iterative solution – small strain¶
A commonly used strategy to solve the nonlinear system, is the iterative NewtonRaphson 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 nonlinear 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
with
To come to an iterative solution, we linearize as this point. This results in
where
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:
with
and
Note
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 nonlinear problem is small.
Nonlinear statics – small strain
Here we employ NewtonRaphson to solve the nonlinear equilibrium equation. It is easy to see that once the above examples have been understood this step is indeed trivial.
Note
NewtonRaphson 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
where
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:
Dynamics¶
Momentum balance¶
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)
Note
To retrieve the original form
But, we can now use also other expressions. For example the damping equivalent of linear elasticity:
with
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:
Mass matrix
Boundary tractions
Internal forces
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.
This results in the following expression for the mass matrix:
Time discretization¶
Here we will discuss several common time discretization steps. To simplify notation we will denote the velocity and the acceleration .
Note
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.
Note
References
Shape functions¶
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 onedimensional 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 nonzero 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 elementbyelement. 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).
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:
or
with
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 integrationpoints). Again, for our internal force:
Note
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.
Note
The details depend on the element type. Several standard elements types are implemented in GooseFEM.
Examples¶
Linear statics – small strain¶
Contents
Mixed boundary conditions¶
Problem description¶
In this first example we will look at a linear elastic, 2D plane strain problem, with imposed boundary displacements. The problem is sketched in the figure below. From which it is observed that:
 The sample is homogeneous.
 Symmetry conditions are assumed in  and direction.
 The outer (right) edge is displaced by in direction (whereby the displacement is constant in direction).
 The boundary tractions are zero everywhere, except for where the displacement is prescribed. There a nonzero reaction force may appear.
 (Not shown) We assume that , and thus . By doing this we can make the problem truly linear.
The geometry is discretized in 2D linear quadrilateral elements, which have four nodes per element. To keep things simple we use only three, equisized, elements in each direction. The mesh is shown below, whereby the element numbers have a regular font while the node numbers are italic.
At this point we focus our attention on the internal force. Thereby we first consider the constitutive response:
I.e the double contraction between the fourthorder stiffness
(with and the bulk and the shear modulus), and the linear strain
(more information in the documentation of GooseSolid). Because of the symmetries in we can simplify the constitutive expression as follows
The displacement of the final configuration, , is now decomposed some known prestrain plus an unknown update :
Note
For our simple problem, which is initially stress and strain free, we find that
Some of the above expressions could thus be simplified further, while also part of the implementation can be omitted. We keep it here to build on this problem later on.
For the internal force this implies that
Since is known we easily evaluate the original (in principle nonlinear) expression of the internal force:
For the update we use the explicit relation for the stress
We then again apply our discretization scheme to obtain
Hence:
The tractions are fully given by the boundary conditions, which also can be used to show that – or at least for the relevant components (see below). The final force balance then reads:
We continue, by writing the problem in terms of scalar degreesoffreedoms (DOFs). Each node has two DOFs, the two vector directions. For our mesh we define the DOFs as follows (whereby the cyan colored DOFnumbers correspond to the direction, while the magenta ones correspond to the direction).
By doing this, a little bit of bookkeeping allows us to write our balance equation as the following system of equations
In short
Next, we have to impose the boundary conditions. Given the fact that we impose displacements on part of the boundary here, a part of will be prescribed. More specifically, our mesh looks as follows
where the yellow nodes are prescribed DOFs while the blue ones are yet unknown. We are now ready to solve the balance equation in two steps. First we will partition the system is unknown and prescribed DOFs:
We are now ready to solve the former part:
which gives us that
and finally that
which can be reassembled as displacement vectors per node, .
Should you be interested, one can compute the reaction force, i.e. the boundary tractions there where the displacement has been prescribed. One has to do the following:
For this, one thus has to compute the new . Because this specific model is linear we can however obtain the reaction forces without having to reevaluate . Specifically
Basic implementation¶
Our first attempt of an implementation literately follows the steps above: it constructs and , which are then partitioned. The prescribed displacements are then set. Thereafter the problem is solved, and the displacements are reconstructed to nodal vector for easy postprocessing.
Prepartitioning¶
One of the things that made the previous examples not very suitable for becoming a production code is the fact that the stiffness matrix was first fully assembled and afterwards partitioned. Besides costing a lot of memory for storing the matrix twice, it might cost a lot of time since partitioning might become a costly operation in the case that sparse matrices are used. To avoid this, the system may be prepartioned. In that case we renumber the DOFs such that we end up with first all the unknown DOFs (denoted iiu
in the code), and then all the known DOFs (denoted iip
). For our example this results in:
This allows us to consider four different matrices (denoted using _uu
, _up
, _pu
, and _pp
) and two different columns (denoted using _u
and _p
) to which the internal force and the stiffness are directly assembled. The rows (and columns) of these matrices and columns follow from introducing separate indices for iiu
and iip
:
Periodic boundary conditions¶
Note
Some additional notes on the theory discussed on a simplified scalar system, for the same mesh as presented here, are included in a separate document. One is invited to study this document before continuing.
Prescribed macroscopic deformation¶
In our first example we will consider the same material and mesh as above. However, now we will assume periodicity is both spatial directions and prescribe a change in the macroscopic deformation gradient, equal to
First of all we will specify periodicity for our mesh. It applies that the following equalities hold (in terms of node numbers:
where are the microscopic fluctuations, that do not affect the macroscopic affine deformation. In terms of DOFs this is can be illustrated as follows:
where the red DOFs are said to be dependent (i.e. they directly follow from the equalities listed above). The simplest this that we can do is construct a system with only the independent DOFs (in blue above) by directly assembling to the independent DOFs. To this end we employ the following DOF numbers:
where the yellow color of the lower left corner indicates that this node is used as reference. Firstly it is used to suppress rigid body deformation. Secondly we apply the macroscopic deformation as the initial condition.
Final equilibrium is then obtained by solving
(which has the dimensions of the number of independent DOFs), and then assembling to the entire system (including the dependent nodes). This is done in the first example, whereby the resulting system is partitioned to deal with the zerodisplacement of the reference node.
Mixed macroscopic boundary conditions¶
[source: periodicvirtualbasic.cpp
]
Here we will enable mixed macroscopic boundary conditions by introduction extra DOFs for the macroscopic deformation gradient tensor, and its antagonist stress response. Since we work in two dimensions we introduce two virtual nodes, each with two DOFs:
We will now employ the following tying relation
To this end we first renumber the system to have all the dependent DOFs at the end
And then partition the system in independent and dependent DOFs
Finally, we obtain the following tying relations for the DOFs
We then employ the nodal dependency to obtain a system of the independent DOFs only:
which, after solving, we can reconstruct to the dependent DOFs using
Note
Towards a production code
Although not (yet) pursued here, it would make sense to partition the system as follows:
Accompanied with the tying relations
And use the following condensation
with the following reconstruction:
Nonlinear statics – small strain¶
Here we extend the example from Linear statics – small strain to a nonlinear constitutive response, which is however still subjected to a small deformations assumption. We treat the constitutive response as a blackbox here:
where it must be emphasized that symmetrizes. To understand more about the constitutive response, please consult the documentation of GooseSolid
Mixed boundary conditions¶
In summary, our iterative update reads
with
and
We will use this to iteratively update
From an ‘initial guess’
Note
This is a bit of particular case. We need this iteration to get things going, however typically
One should not be fooled, this does not mean to an equilibrium has been obtained.
Like before, we introduce DOFs to make the system scalar. We then need to partition the system to deal with the prescribed displacementcomponents:
From which the update of the unknown DOFs as follows
There a very important concept hidden here. Because we prescribe directly to the correct value – which we do not iteratively update – we need to set
(Which means that one does not have to compute the product for , as it will be zero.)
Note
Reaction forces
To obtain the reaction forces on the prescribed DOFs simply use that
Which should be evaluated once convergence has been reached (before that, this has no meaning).
Periodic problem¶
In some ways the periodic example is even easier than the example above. As discussed previously only leads to periodic fluctuations. In summary we begin by setting
the update
then only affects the fluctuations, and not the average. Also we do not need to worry about the above discussion on the prescribed displacements since they are always zero and traction free – as they merely suppress rigid body modes.
Note
The periodic implementation with macroscopic DOFs is identically extended to the nonlinear case as above. Here one does have care about properly adding the prescribed displacements.
Diagonal mass matrix¶
‘Theory’¶
Note
Should you still want to use the full (sparse) matrix in conjunction with a solver, please consider using the Intel compile with the Intel MKL library (which includes the Pardiso solver). Consider this example.
Because dynamic computations need so many increments, it may be worth to cut down the computational cost of a simulation. The best way to start cutting is on the solver, which is often the most costly of the entire Finite Element Program. Remember that
If we make the mass matrix diagonal, can be obtained without solving a system. Instead we merely needed the component wise inverse of the diagonal terms (which are the only remaining nonzero terms).
There are essentially four possible definitions of the mass matrix:
The consistent (or complete) mass matrix. Here the quadrature of the mass matrix is identical to that of the internal force in the system. It is the result of numerical integration of
The lumped mass matrix. In which
The lumped mass matrix with diagonal scaling:
where the constant is set such that:
A trick to accomplish this is to set
where the denominator is the trace of .
Evaluating of using a quadrature rule involving only the nodal points and thus automatically yielding a diagonal matrix. It relies on the fact that when the quadraturepoints coincide with the nodes, a local support is obtained (since the shape functions are zero in all the other nodes). The integration point volume is that part of the element’s volume that can be associated to that node.
For 2D linear triangles, which have three nodes and normally one Gauss point in the middle of the element, the number of integration points is increased to three (whose positions correspond to that of the nodes), and their weight is set to 1/3 of the normal weight. Interestingly, for this case all possible definitions of the diagonalized mass matrix are the same.
For 2D quads only the position of the integration points changes to that of the nodes.