No Matches
VEC-0: Linear elasticity
Prerequisites of this tutorial include MSH-1: Create a 2D mesh from Gmsh and SCL-1: Poisson's equation (homogeneous BC)

Intended learning outcome:
  • general structure of a program developed using MoFEM
  • idea of Simple Interface in MoFEM and how to use it
  • implementing vector valued problems like linear elasticity
  • implementing boundary conditions specified on the part of the boundary
  • developing code that can be compiled for 2D or 3D cases
  • use of default forms integrators
  • pushing UDOs to the Pipeline
  • utilisation of tools to convert outputs (MOAB) and visualise them (Paraview)

The solution of the linear elasticity problem is presented in this tutorial. Lets consider an isotropic elastic cantilever beam with prescribed gravity load as presented in Figure 1.

Figure 1: Cantilever beam considered in this example.

Strong form

In order to compute displacement vector field \(\mathbf{u}\), every point in the beam has to satisfy balance of linear momentum and boundary conditions as follows:

\[ \begin{align} \label{eq:momentum_balance} \begin{cases} \nabla \cdot \boldsymbol \sigma \left(\mathbf u(\mathbf{x}) \right) + \mathbf b =0 & \text { in } \Omega \\ \mathbf u(\mathbf{x}) = \bar{\mathbf{u}} & \text { on } \partial\Omega^\mathbf{u} \\ \boldsymbol \sigma \cdot \mathbf n = \bar {\mathbf t} & \text { on } \partial\Omega^\mathbf{t} \end{cases} \end{align} \]

where \(\mathbf{b}\) are body forces, \( \bar{\mathbf{u}} \) are fixed displacements, \( \bar{\mathbf{t}} \) is traction vector and \(\mathbf{n}\) is the unit normal vector on the boundary and \( \boldsymbol \sigma \) is the Cauchy stress tensor. In case of linear elasticity, small strains are assumed whereby the stress measure is defined as:

\[ \begin{align} \boldsymbol \sigma =& \mathbb D : \frac{1}{2} (\nabla \mathbf u^\text{T} + \nabla \mathbf u) \end{align} \]

where \( \mathbb D \) is the 4th order elasticity tensor.


Following a standard Galerkin approach as in SCL-1: Poisson's equation (homogeneous BC), the governing equations are discretised into weak form as follows:

\[ \begin{align} \int_{\Omega} \nabla \delta \mathbf u : \boldsymbol \sigma(\mathbf u) \, \mathrm d \Omega= - \int_{\Omega} \delta \mathbf u \cdot \mathbf b \, \mathrm d \Omega + \int_{\partial{\Omega}^\mathbf t} \delta \mathbf u \cdot \bar{\mathbf t} \, \mathrm d S \quad \forall\, \delta \mathbf u \in \mathbf{H}^1_0(\Omega) \end{align} \]

To approximate displacements \( \mathbf u\) and variance of displacements \( \delta \mathbf u\), once again we will use MoFEM's hierarchical shape functions:

\[ \begin{align} \mathbf u \approx \mathbf u^{h}=\sum_{j=0}^{N-1} \phi_{j} \bar{\mathbf{u}}_{j}, \quad \quad \delta \mathbf u \approx \delta \mathbf u^{h}=\sum_{j=0}^{N-1} \phi_{j} \delta \bar{\mathbf{u}}_{j} \end{align} \]

By substituting the above definitions into the weak form, the global system \( \mathbf{K U}=\mathbf{F}\) can be assembled, where elements contributions from every finite elements are computed as follows:

\[ \begin{align} \mathbf K_{i j}^{e}&=\int_{\Omega^{e}} \nabla \phi_{i} \mathbb D \nabla \phi_{j} \, \mathrm d \Omega \\ \mathbf F_{i}^{e}&=\int_{\partial \Omega^{e}} \phi_{i} \mathbf t \, \mathrm d S + \int_{\Omega^{e}} \phi_{i} \mathbf b \, \mathrm d \Omega \end{align} \]

The Gauss quadrature is utilised to calculate the integrals.


In this section we will focus on the procedure of setting up and solving vector-valued problems like elasticity in MoFEM using MoFEM::Simple and Form Integrators, which are template UDOs used for common operations.
The code developed for this example can be compiled for both 2D and 3D cases by changing the following constant variable:

The constant variable SPACE_DIM will be used throughout the code to determine dimension of the executable. It is also possible to setup CMakeFiles.txt to compile both, 2D and 3D executables (see plastic.cpp example).

Note that based on SPACE_DIM an appropriate type of elements is defined in the following snippet:

The main program consists of a series of functions when invoking the function Example::runProblem():

MoFEMErrorCode Example::runProblem() {
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: helmholtz.cpp:106
MoFEMErrorCode assembleSystem()
[Applying essential BC]
Definition: helmholtz.cpp:154
MoFEMErrorCode readMesh()
[run problem]
Definition: helmholtz.cpp:73
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: helmholtz.cpp:299
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Definition: helmholtz.cpp:235
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:176
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:188
MoFEMErrorCode outputResults()
Definition: helmholtz.cpp:255

where Example::readMesh() executes the mesh reading procedure.

Problem setup

The function Example::setupProblem() involves adding the displacement field U on the entities of the mesh and determine the base function space, base, field rank and base order. We are adding U as both domain and a boundary field.

MoFEMErrorCode Example::setupProblem() {
Simple *simple = mField.getInterface<Simple>();
// Add field
int order = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("U", order);
CHKERR simple->setFieldOrder("GEOMETRY", 2);
CHKERR simple->setUp();
auto project_ho_geometry = [&]() {
Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
return mField.loop_dofs("GEOMETRY", ent_method);
CHKERR project_ho_geometry();
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
constexpr int SPACE_DIM
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::Interface & mField
Definition: plastic.cpp:148
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

The space chosen is H1, the base is AINSWORTH_LEGENDRE_BASE and dimension is SPACE_DIM that corresponds to a vector field with SPACE_DIM number of coefficients. Additionally, we define a command line parameter -order that specifies the order of our base, 2nd order is chosen as a default.

Common data

The next function in our main program: Example::createCommonData() defines and initialises the data passed between UDOs. In case of elasticity, we compute elasticity tensor \( \mathbb D\) using the formula below:

\[ \begin{align} \mathbb D_{ijkl} = G \left[ \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk} \right] + A (K - \frac{2}{3} G) \left[ \delta_{ij} \delta_{kl} \right] \end{align} \]

where \( K\) and \( G\) are bulk and shear modulus, respectively. The coefficient \( A\) depends on the dimension of the problem, for 3D cases and plane strain formulation it is simply \( A = 1\), whereas for plane stress it takes the following form:

\[ \begin{align} A=\frac{2 G}{K+\frac{4}{3} G} \end{align} \]

Using the FTensor::Index type these formulas can be directly implemented using Einstein's summation convention as shown below:

auto set_material_stiffness = [&]() {
double A = (SPACE_DIM == 2)
: 1;
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
t_D(i, j, k, l) =
2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
t_kd(k, l);
Kronecker Delta class symmetric.
constexpr double shear_modulus_G
Definition: elastic.cpp:60
constexpr double bulk_modulus_K
Definition: elastic.cpp:59
constexpr auto t_kd
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr AssemblyType A

Note that depending on SPACE_DIM our code will appropriately sum over all the indices and set coefficient \( A \) for plane stress or 3D case. In the next lines of the program, we define a vector that will be passed to the UDOs and specify the gravity load. For this example, we will hardcode the load to be always 1 in the Y direction \( \mathbf{ b } = [ 0,-1,0 ]\)

It is also possible to implement passing that data from command line parameter or read it from the mesh blocksets.

Finally, Example::createCommonData() we initialise containers that will be storing data between different UDOs. For elasticity problem we will be storing gradient of deformation, small strain tensor, stress tensor, elasticity tensor and body force vector.

In the last two lines we call previously defined lambda functions for calculating elasticity tensor and body force vector.

CHKERR set_material_stiffness();
CHKERR set_body_force();

Essential boundary conditions

In order to apply homogeneous Dirichlet boundary conditions for our problem, we will utilise MoFEM's utility to remove degrees of freedom from a problem.

MoFEMErrorCode Example::boundaryCondition() {
auto pip = mField.getInterface<PipelineManager>();
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3);
pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
// Essential boundary condition.
CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
auto get_pre_proc_hook = [&]() {
return EssentialPreProc<DisplacementCubitBcData>(mField,
pip->getDomainRhsFE(), {});
pip->getDomainRhsFE()->preProcessHook = get_pre_proc_hook();
pip->getBoundaryRhsFE()->preProcessHook = get_pre_proc_hook();
// Infernal forces
auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
new OpSymmetrizeTensor<SPACE_DIM>("U", mat_grad_ptr, mat_strain_ptr));
auto mat_D_ptr = boost::make_shared<MatrixDouble>();
CHKERR addMatBlockOps(pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC",
mat_D_ptr, Sev::inform);
new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
"U", mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
new OpInternalForce("U", mat_stress_ptr,
[](double, double, double) constexpr { return -1; }));
// Body forces
pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
// Add force boundary condition
pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
// Add case for mix type of BCs
pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::verbose);
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
Definition: elastic.cpp:36
MoFEMErrorCode addMatBlockOps(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev)
Definition: elastic.cpp:86

The first lambda function in the snippet above fix_disp searches for entities in a meshset with a given name. Subsequently, these entities (Range) are passed to another lambda function - remove_ents that removes them from the problem. This operation is equivalent of setting all the degrees of freedom to zero on associated entities. For the examples considered in this tutorial, the provided meshes contain only one meshset FIX_ALL that will essentially fix displacements in all directions to zero. In Figure 2 we show selected triangles on the face of triangle that are used to define the FIX_ALL meshset. Note that in case of 2D mesh, we select edge elements.

Figure 2: FIX_ALL block meshset specified on the mesh of the cantilever for 3D case.

Natural boundary conditions

For natural boundary conditions we will only specify simple gravity load on the entire domain. We have to push into our Domain RHS pipeline an operator that calculates the following integral:

\[ \begin{align} \mathbf F_{i}^{e}&=\int_{\Omega^{e}} \phi_{i} \mathbf b \, \mathrm d V \end{align} \]

This operation can be done by pushing OpBodyForce User Data Operator into RHS domain pipeline. Note that arguments of this operator are: name of the field U, container bodyForceMatPtr for vector of body force and lambda function which is defines a coefficient for our vector. The function takes three arguments, which are the coordinates of the integration points. However, in this example, we assume homogeneous body, therefore the coefficient can by simply set to 1.

The OpBodyForce is an alias of the specialisation of the form integrator template OpBaseTimesVector declared in the beginning of the file

using OpBodyForce = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 0>;
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce

The first parameter in the above template specifies the dimension of the base (we use scalar base, so set to 1), second parameter requires specification of the number of the coefficients of we field we are operating on (our vector field of displacements has SPACE_DIM coefficients). The last parameter specifies whether the passed vector is constant for all integration points. In this case bodyForceMatPtr body force vector is assumed constant, (has the same values for all gauss points). However, if the considered body is heterogeneous, then the last parameter would be set to 1.

Pushing operators to LHS domain pipeline

Another function of the main program Example::assembleSystem() is responsible for defining pipelines used to evaluate linear and bilinear forms of the system and assembly to the global matrices and vectors.

MoFEMErrorCode Example::assembleSystem() {
auto pip = mField.getInterface<PipelineManager>();
auto simple = mField.getInterface<Simple>();
pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
auto mat_D_ptr = boost::make_shared<MatrixDouble>();
CHKERR addMatBlockOps(pip->getOpDomainLhsPipeline(), "U", "MAT_ELASTIC",
mat_D_ptr, Sev::verbose);
pip->getOpDomainLhsPipeline().push_back(new OpK("U", "U", mat_D_ptr));
pip->getOpDomainLhsPipeline().push_back(new OpUnSetBc("U"));
auto integration_rule = [](int, int, int approx_order) {
return 2 * (approx_order - 1);
auto integration_rule_bc = [](int, int, int approx_order) {
return 2 * approx_order;
CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
CHKERR pip->setDomainLhsIntegrationRule(integration_rule);
CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
auto integration_rule
static constexpr int approx_order

In case of linear elasticity, after pushing boundary and source term operators (gravity load), we only have to define UDOs for LHS pipeline of the domain. In case of 2D problem we push operators (MoFEM::OpCalculateHOJacForFace) that calculate and set the inverse of Jacobian for 2D elements.

The last operator OpK integrates the bilinear form:

\[ \begin{align} \mathbf K^e_{ij} = \int_{\Omega^e} \nabla \phi_i \mathbb D \nabla \phi_j \, \mathrm d \Omega \end{align} \]

and subsequently assembles it into global system matrix \( \mathbf K \). Note that once again, a new operator does not have to be defined as OpK is simply an alias to existing operator from MoFEM's repertoire. To that operator we are passing previously defined container that stores elasticity tensor \( \mathbb D \).

using OpK = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;

The utilised OpGradSymTensorGrad form integrator template of Bilinear Form has four parameters: rank of the base, rank of row field, rank of the column field and finally, the increment of elasticity tensor, which in case of homogeneous material can be set to 0.

More info about available Form Integrators can be found in User data operators table.

Solver setup

With appropriately defined pipelines for assembling the global system of equations \( \mathbf{K} \mathbf{U} = \mathbf{F} \), the next functions Example::solveSystem() sets up PETSc solver and calculates it directly or iteratively based on specified settings.

struct SetUpSchur {
static boost::shared_ptr<SetUpSchur>
virtual MoFEMErrorCode setUp(SmartPetscObj<KSP> solver) = 0;
virtual MoFEMErrorCode preProc() = 0;
virtual MoFEMErrorCode postProc() = 0;
SetUpSchur() = default;
MoFEMErrorCode Example::solveSystem() {
auto simple = mField.getInterface<Simple>();
auto pip = mField.getInterface<PipelineManager>();
auto solver = pip->createKSP();
CHKERR KSPSetFromOptions(solver);
auto dm = simple->getDM();
auto D = smartCreateDMVector(dm);
auto F = smartVectorDuplicate(D);
auto setup_and_solve = [&]() {
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
CHKERR KSPSetUp(solver);
MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
CHKERR KSPSolve(solver, F, D);
MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
MOFEM_LOG_TAG("TIMER", "timer");
if (A == AssemblyType::SCHUR) {
CHKERR schur_ptr->setUp(solver);
pip->getDomainLhsFE()->preProcessHook = [&]() {
if (schur_ptr)
CHKERR schur_ptr->preProc();
pip->getBoundaryLhsFE()->postProcessHook = [&]() {
if (schur_ptr)
CHKERR schur_ptr->postProc();
CHKERR setup_and_solve();
} else {
CHKERR setup_and_solve();
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:511
#define MOFEM_LOG(channel, severity)
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
double D
Deprecated interface functions.
[Push operators to pipeline]
Definition: plastic.cpp:736
virtual MoFEMErrorCode setUp(KSP solver)=0
virtual MoFEMErrorCode postProc()=0
virtual MoFEMErrorCode preProc()=0
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM >, SmartPetscObj< IS >, SmartPetscObj< AO >)
Definition: plastic.cpp:1351

Note that in case of simple linear problems like elasticity considered therein the above snippet will look very similar. A detailed description of each function can be found e.g. in SCL-1: Poisson's equation (homogeneous BC).

Postprocessing pipeline

Finally, once the solver calculations are completed and solution vector is obtained, it is necessary save the data on the visualisation mesh and compute quantities of interest like strains and stresses. In function Example::outputResults() we will define pipeline to postprocess the results and save the calculated values to the output mesh.

We can reuse LHS domain pipeline, for the purpose of postprocessing. The following line resets (cleans) the operators pushed previously


In the next two lines we create a postprocessig element (object) and use its function to generate a reference mesh.

auto post_proc_fe = boost::make_shared<PostProcEle>(mField);

Next, a series of operators is pushed:

  • OpCalculateHOJacForFace computes and sets the Jacobian in case of 2D elements
  • OpCalculateVectorFieldGradient computes the gradient of the obtained displacement field \( \nabla \mathbf u^{h}=\sum_{m=0}^{N-1} \nabla \phi_m \bar{\mathbf u}_m \)
  • OpSymmetrizeTensor symmetrises the previously computed gradient to compute small strain tensor \( \boldsymbol \varepsilon = \frac{1}{2} (\nabla \mathbf u^\text{T} + \nabla \mathbf u) \)
  • OpTensorTimesSymmetricTensor is another form integrator that computes Cauchy stress: \( \boldsymbol \sigma = \mathbb D : \boldsymbol \varepsilon \)
  • Tutorial::OpPostProcElastic saves the tensors \( \boldsymbol \varepsilon \) and \( \boldsymbol \sigma \) onto the reference mesh using moab tags.

Running code and visualisation

In order to run the program that we have been discussing in this tutorial, you will do the following steps

  • First, go to the directory where the binary file named elastic is located. Depending on how you install MoFEM shown in this page Installation, going to the directory would be something similar to this
    • For user version installation
      cd mofem_install/um_view/tutorials/vec-0/
    • For developer version installation
      cd mofem_install/mofem-cephas/mofem/users_modules/um-build-RelWithDebInfo-abcd1234/tutorials/vec-0
  • Second, check the parameters in the param_file.petsc. These are PETSc parameters and you should only use parameters that are needed for a particular solver, in this case KSP solver. Only the following parameters should be uncommented
    ## Linear solver
    -ksp_type fgmres
    -pc_type lu
    -pc_factor_mat_solver_type mumps
  • Third, in the terminal, run commands to partition the input mesh and start the analysis
    ./elastic_2d -file_name beam_3D.cub -order 2
    where the mesh of a square plate is partitioned into two parts and then the program is run using two processors (the same number of partitions) with fourth order polynomial of approximation.

This will produce out_elastic.h5m file which we can convert using

mbconvert out_elastic.h5m out_elastic_3D.vtk

The file out_elastic_3D.vtk can be opened in Paraview. Using the filter WarpByScalar allows to visualise the deformation as shown in Figure 3.

We can also compute 2D case of the beam. In order to do that, change the below variable:

to the value equal to 2 and compile the code (in your build directory) using make command. You should be able then to compute 2D variant as follows

./elastic_2d -file_name beam_2D.cub -order 2

and generate vtk file using:

mbconvert out_elastic.h5m out_elastic_2D.vtk

Figure 3: Deformation of cantilever beams for 2D and 3D case. The color map represents displacements U.

Note that the beams, both for 2D and 3D cases experience spurious deformation at the free end. This is a known feature of small strain kinematic description. To mitigate this issue, a large strain geometrically nonlinear formulation can be used. The extension of the above tutorial to such nonlinear case can be found in VEC-2: Nonlinear elastic. We propose there a simple Hencky strain measure that improves the behaviour of the material experiencing large rotations.

Plain program

The full source code of the main program can be accessed below: elastic.cpp