v0.14.0 |

VEC-0: Linear elasticity

- Note
- Prerequisites of this tutorial include MSH-1: Create a 2D mesh from Gmsh and SCL-1: Poisson's equation (homogeneous BC)
- High order geometry and schurs complement are present in the current code but outside the scope of this tutorial refer to ADV-0: Plastic problem

- Note
- 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.

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 to setup and solve vector-valued problems such as 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 assigning the input parameter *EXECUTABLE_DIMENSION* in the *CMakeFiles.txt* prior to compilation and is used throughout the code to define the constant variable *SPACE_DIM* which is used as a template paramenter to assign the dimension of the executable:

Based on *SPACE_DIM* appropriate element types are defined in the following snippet for example, *DomainEle* in 2D are defined as faces while in 3D they are defined as volumes:

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

where Example::readMesh() executes the mesh reading procedure, refer to MSH-1: Create a 2D mesh from Gmsh and MSH-2: Create a 3D mesh from Gmsh for an introduction to mesh generation.

The first function, Example::setupProblem() is used to setup the problem and involves adding the displacement field *U* on the entities of the mesh and determine the base function space, base, field rank and base order. For this problem We are adding displacement field *U* as both domain and a boundary field.

The space selected for displacement field U is `H1`

and dimension is *SPACE_DIM* that corresponds to a vector field with *SPACE_DIM* number of coefficients. Additionally, we define command line parameters *-base* and *-order* that specify the base and order of our base respectively with `AINSWORTH_LEGENDRE_BASE`

and 2nd order selected by default.

The next function, Example::boundaryCondition() is used to apply homogeneous Dirichlet boundary conditions for our problem by utilise MoFEM's utility to remove degrees of freedom from a problem through the *BcManager* interface.

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 triangles that are used to define the *FIX_ALL* meshset. While for the 2D case, we select edge elements instead.

The next function of the main program, Example::assembleSystem() is responsible for defining pipelines used to evaluate linear and bilinear forms of the system (refer to User data operators table) and assembly of the global matrices and vectors.

The first part of this function, sets and pushes the integration rule for the finite element method to the pipeline.

The second part of this function, first calculates the local elasticity tensor through `addMatBlockOps`

function depending on the `MAT_ELASTIC`

block defined in the mesh file, then pushes the LHS stiffness matrix using alias `OpK`

to the LHS pipeline.

The local elasticty tensor \( \mathbb D\) is calculate 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, within the Example::addMatBlockOps() function:

- Note
- 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.

The 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 \).

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.

The next part of the Example::assembleSystem() function calculates the internal forces for the system which is required to enforce the boundary conditions (defined as \( \mathbf K \mathbf u_e \) refer to SCL-1: Poisson's equation (homogeneous BC) for an in depth description). For elastic problems, this requires calculating the displacement gradient \( \nabla \mathbf u \), strain \( \boldsymbol \varepsilon \) and stress \( \boldsymbol \sigma \) of the system by pushing the relevant operatiors to the RHS pipeline.

This is achieved by pushing a series of operators as follows:

`OpCalculateVectorFieldGradient`

computes the gradient of the 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 \)`OpInternalForce`

is an alias of the FormsIntegrators`OpGradTimesSymTensor`

which computes \( (\mathbf K \mathbf u_e)_j = \nabla \phi_i \mathbf \sigma_{ij} \)

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 utilsing the `DomanRhsBc::AddFluxToPipeline`

alias which applies a flux User Data Operator into RHS pipeline as shown below.

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.

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).

Finally, once the solver calculations are completed and solution vector is obtained, it is necessary to 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 for each pipeline.

In the next few lines we create various postprocessing objects including the reference mesh and lists of elements to postprocess.

Then we define a series of functions including, `calculate_stress_ops`

, `post_proc_domain`

and `post_proc_boundary`

which push the operators required to output the displacement \( \mathbf{u} \), current position \( \mathbf{x} \), strain \( \boldsymbol{\varepsilon} \) and stress \( \boldsymbol{\sigma} \) for the domain and boundary.

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 with Spack (Scripts), 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

- For user version installation
- 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-ksp_monitor - Third, in the terminal, run commands to partition the input mesh and start the analysis where a second order polynomial of approximation is selected../elastic_3d -file_name beam_3D.cub -order 2

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 `WarpByVector`

allows to visualise the deformation as shown in Figure 3.

We can also compute 2D case of the beam. Following the same process as that for 3D with the commands:

./elastic_2d -file_name beam_2D.cub -order 2

and generate vtk file using:

mbconvert out_elastic.h5m out_elastic_2D.vtk

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.

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

Generated by Doxygen 1.8.17 and hosted at