The example shows how to solve the linear elastic problem.

/** \file simple_elasticity.cpp
* \example simple_elasticity.cpp
The example shows how to solve the linear elastic problem.
/* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* License for more details.
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <BasicFiniteElements.hpp>
using namespace boost::numeric;
using namespace MoFEM;
static char help[] = "-order approximation order\n"
// Finite element stiffness sub-matrix K_ij
// Elastic stiffness tensor (4th rank tensor with minor and major symmetry)
// Young's modulus
double yOung;
// Poisson's ratio
double pOisson;
OpK(bool symm = true)
symm) {
// Evaluation of the elastic stiffness tensor, D
// hardcoded choice of elastic parameters
pOisson = 0.1;
yOung = 10;
// coefficient used in intermediate calculation
const double coefficient = yOung / ((1 + pOisson) * (1 - 2 * pOisson));
tD(i, j, k, l) = 0.;
tD(0, 0, 0, 0) = 1 - pOisson;
tD(1, 1, 1, 1) = 1 - pOisson;
tD(2, 2, 2, 2) = 1 - pOisson;
tD(0, 1, 0, 1) = 0.5 * (1 - 2 * pOisson);
tD(0, 2, 0, 2) = 0.5 * (1 - 2 * pOisson);
tD(1, 2, 1, 2) = 0.5 * (1 - 2 * pOisson);
tD(0, 0, 1, 1) = pOisson;
tD(1, 1, 0, 0) = pOisson;
tD(0, 0, 2, 2) = pOisson;
tD(2, 2, 0, 0) = pOisson;
tD(1, 1, 2, 2) = pOisson;
tD(2, 2, 1, 1) = pOisson;
tD(i, j, k, l) *= coefficient;
* \brief Do calculations for give operator
* @param row_side row side number (local number) of entity on element
* @param col_side column side number (local number) of entity on element
* @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
* @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
* @param row_data data for row
* @param col_data data for column
* @return error code
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
// get number of dofs on row
nbRows = row_data.getIndices().size();
// if no dofs on row, exit that work, nothing to do here
if (!nbRows)
// get number of dofs on column
nbCols = col_data.getIndices().size();
// if no dofs on Columbia, exit nothing to do here
if (!nbCols)
// K_ij matrix will have 3 times the number of degrees of freedom of the
// i-th entity set (nbRows)
// and 3 times the number of degrees of freedom of the j-th entity set
// (nbCols)
K.resize(nbRows, nbCols, false);
// get number of integration points
nbIntegrationPts = getGaussPts().size2();
// check if entity block is on matrix diagonal
if (row_side == col_side && row_type == col_type) {
isDiag = true;
} else {
isDiag = false;
// integrate local matrix for entity block
CHKERR iNtegrate(row_data, col_data);
// assemble local matrix
CHKERR aSsemble(row_data, col_data);
int nbRows; ///< number of dofs on rows
int nbCols; ///< number if dof on column
int nbIntegrationPts; ///< number of integration points
bool isDiag; ///< true if this block is on diagonal
* \brief Integrate B^T D B operator
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
iNtegrate(DataForcesAndSourcesCore::EntData &row_data,
// get sub-block (3x3) of local stiffens matrix, here represented by second
// order tensor
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
&m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
&m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
// get element volume
double vol = getVolume();
// get intergration weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
// iterate over integration points
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
// calculate scalar weight times element volume
const double a = t_w * vol;
// iterate over row base functions
for (int rr = 0; rr != nbRows / 3; ++rr) {
// get sub matrix for the row
auto t_m = get_tensor2(K, 3 * rr, 0);
// get derivatives of base functions for columns
auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
// iterate column base functions
for (int cc = 0; cc != nbCols / 3;++cc) {
// integrate block local stiffens matrix
t_m(i, k) +=
a * (tD(i, j, k, l) * (t_row_diff_base(j) * t_col_diff_base(l)));
// move to next column base function
// move to next block of local stiffens matrix
// move to next row base function
// move to next integration weight
* \brief Assemble local entity block matrix
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
// get pointer to first global index on row
const int *row_indices = &*row_data.getIndices().data().begin();
// get pointer to first global index on column
const int *col_indices = &*col_data.getIndices().data().begin();
Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
: getFEMethod()->snes_B;
// assemble local matrix
CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
&*K.data().begin(), ADD_VALUES);
if (!isDiag && sYmm) {
// if not diagonal term and since global matrix is symmetric assemble
// transpose term.
K = trans(K);
CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
&*K.data().begin(), ADD_VALUES);
double pressureVal;
OpPressure(const double pressure_val = 1)
pressureVal(pressure_val) {}
// vector used to store force vector for each degree of freedom
MoFEMErrorCode doWork(int side, EntityType type,
// check that the faces have associated degrees of freedom
const int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
// size of force vector associated to the entity
// set equal to the number of degrees of freedom of associated with the
// entity
nF.resize(nb_dofs, false);
// get number of gauss points
const int nb_gauss_pts = data.getN().size1();
// create a 3d vector to be used as the normal to the face with length equal
// to the face area
auto t_normal = getFTensor1Normal();
// get intergration weights
auto t_w = getFTensor0IntegrationWeight();
// vector of base functions
auto t_base = data.getFTensor0N();
// loop over all gauss points of the face
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
// weight of gg gauss point
double w = 0.5 * t_w;
// create a vector t_nf whose pointer points an array of 3 pointers
// pointing to nF memory location of components
for (int bb = 0; bb != nb_dofs / 3; ++bb) {
// scale the three components of t_normal and pass them to the t_nf
// (hence to nF)
t_nf(i) += (w * pressureVal * t_base) * t_normal(i);
// move the pointer to next element of t_nf
// move to next base function
// move to next integration weight
// add computed values of pressure in the global right hand side vector
CHKERR VecSetValues(getFEMethod()->ksp_f, nb_dofs, &data.getIndices()[0],
&nF[0], ADD_VALUES);
Range fixFaces, fixNodes, fixSecondNode;
ApplyDirichletBc(const Range &fix_faces, const Range &fix_nodes,
const Range &fix_second_node)
: MoFEM::FEMethod(), fixFaces(fix_faces), fixNodes(fix_nodes),
fixSecondNode(fix_second_node) {
// constructor
MoFEMErrorCode postProcess() {
std::set<int> set_fix_dofs;
for (_IT_NUMEREDDOF_ROW_FOR_LOOP_(problemPtr, dit)) {
if (dit->get()->getDofCoeffIdx() == 2) {
if (fixFaces.find(dit->get()->getEnt()) != fixFaces.end()) {
if (fixSecondNode.find(dit->get()->getEnt()) != fixSecondNode.end()) {
if (dit->get()->getDofCoeffIdx() == 1) {
if (fixNodes.find(dit->get()->getEnt()) != fixNodes.end()) {
std::vector<int> fix_dofs(set_fix_dofs.size());
std::copy(set_fix_dofs.begin(), set_fix_dofs.end(), fix_dofs.begin());
CHKERR VecAssemblyBegin(ksp_f);
CHKERR VecAssemblyEnd(ksp_f);
Vec x;
CHKERR VecDuplicate(ksp_f, &x);
CHKERR VecZeroEntries(x);
CHKERR MatZeroRowsColumns(ksp_B, fix_dofs.size(), &fix_dofs[0], 1, x,
CHKERR VecDestroy(&x);
* \brief Set integration rule to volume elements
* This rule is used to integrate \f$\nabla v \cdot \nabla u\f$, thus
* if the approximation field and the testing field are polynomials of order "p",
* then the rule for the exact integration is 2*(p-1).
* Integration rule is order of polynomial which is calculated exactly. Finite
* element selects integration method based on return of this function.
struct VolRule {
int operator()(int, int, int p) const {
return 2 * (p - 1);
* \brief Set integration rule to boundary elements
* This rule is used to integrate the work of external forces on a face,
* i.e. \f$f \cdot v\f$, where f is the traction vector and v is the test
* vector function. The current problem features a Neumann bc with
* a pre-defined constant pressure. Therefore, if the test field is
* represented by polynomials of order "p", then the rule for the exact
* integration is also p.
* Integration rule is order of polynomial which is calculated exactly. Finite
* element selects integration method based on return of this function.
struct FaceRule {
int operator()(int, int, int p) const {
return p;
int main(int argc, char *argv[]) {
const string default_options = "-ksp_type fgmres \n"
"-pc_type lu \n"
"-pc_factor_mat_solver_package mumps \n"
string param_file = "param_file.petsc";
if (!static_cast<bool>(ifstream(param_file))) {
std::ofstream file(param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file << default_options;
MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
// Create mesh database
moab::Core mb_instance; // create database
moab::Interface &moab = mb_instance; // create interface to database
try {
// Create MoFEM database and link it to MoAB
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
// Get command line options
int order = 3; // default approximation order
PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "SimpleElasticProblem",
// Set approximation order
CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
// Set testing (used by CTest)
CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
&flg_test, PETSC_NULL);
ierr = PetscOptionsEnd();
Simple *simple_interface = m_field.getInterface<MoFEM::Simple>();
CHKERR simple_interface->getOptions();
CHKERR simple_interface->loadFile();
Range fix_faces, pressure_faces, fix_nodes, fix_second_node;
enum MyBcTypes {
EntityHandle meshset = bit->getMeshset();
int id = bit->getMeshsetId();
if (id == FIX_BRICK_FACES) { // brick-faces
CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 2,
fix_faces, true);
Range adj_ents;
CHKERR m_field.get_moab().get_adjacencies(fix_faces, 0, false, adj_ents,
CHKERR m_field.get_moab().get_adjacencies(fix_faces, 1, false, adj_ents,
} else if (id == FIX_NODES) { // node(s)
CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 0,
fix_nodes, true);
} else if (id == BRICK_PRESSURE_FACES) { // brick pressure faces
CHKERR m_field.get_moab().get_entities_by_dimension(
meshset, 2, pressure_faces, true);
} else if (id ==
FIX_NODES_Y_DIR) { // restrained second node in y direction
CHKERR m_field.get_moab().get_entities_by_dimension(
meshset, 0, fix_second_node, true);
} else {
CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
CHKERR simple_interface->setFieldOrder("U", order);
CHKERR simple_interface->defineFiniteElements();
// Add pressure element
CHKERR m_field.add_finite_element("PRESSURE");
CHKERR m_field.modify_finite_element_add_field_row("PRESSURE", "U");
CHKERR m_field.modify_finite_element_add_field_col("PRESSURE", "U");
CHKERR m_field.modify_finite_element_add_field_data("PRESSURE", "U");
CHKERR simple_interface->defineProblem();
DM dm;
CHKERR simple_interface->getDM(&dm);
CHKERR simple_interface->buildFields();
CHKERR simple_interface->buildFiniteElements();
CHKERR m_field.add_ents_to_finite_element_by_dim(pressure_faces, 2,
CHKERR m_field.build_finite_elements("PRESSURE", &pressure_faces);
CHKERR simple_interface->buildProblem();
// Create elements instances
boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
// Set integration rule to elements instances
elastic_fe->getRuleHook = VolRule();
pressure_fe->getRuleHook = FaceRule();
// Add operators to element instances
// push operators to elastic_fe
elastic_fe->getOpPtrVector().push_back(new OpK());
// push operators to pressure_fe
pressure_fe->getOpPtrVector().push_back(new OpPressure());
boost::shared_ptr<FEMethod> fix_dofs_fe(
new ApplyDirichletBc(fix_faces, fix_nodes, fix_second_node));
boost::shared_ptr<FEMethod> null_fe;
// Set operators for KSP solver
dm, simple_interface->getDomainFEName(), elastic_fe, null_fe, null_fe);
CHKERR DMMoFEMKSPSetComputeRHS(dm, "PRESSURE", pressure_fe, null_fe,
// initialise matrix A used as the global stiffness matrix
Mat A;
// initialise left hand side vector x and right hand side vector f
Vec x, f;
// allocate memory handled by MoFEM discrete manager for matrix A
CHKERR DMCreateMatrix(dm, &A);
// allocate memory handled by MoFEM discrete manager for vector x
CHKERR DMCreateGlobalVector(dm, &x);
// allocate memory handled by MoFEM discrete manager for vector f of the
// same size as x
CHKERR VecDuplicate(x, &f);
// precondition matrix A according to fix_dofs_fe and elastic_fe finite
// elements
elastic_fe->ksp_B = A;
fix_dofs_fe->ksp_B = A;
// precondition the right hand side vector f according to fix_dofs_fe and
// elastic_fe finite elements
fix_dofs_fe->ksp_f = f;
pressure_fe->ksp_f = f;
CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
CHKERR DMoFEMLoopFiniteElements(dm, "PRESSURE", pressure_fe);
// This is done because only post processor is implemented in the
// ApplyDirichletBc struct
CHKERR DMoFEMPostProcessFiniteElements(dm, fix_dofs_fe.get());
// make available a KSP solver
KSP solver;
// make the solver available for parallel computing by determining its MPI
// communicator
// making available running all options available for KSP solver in running
// command
CHKERR KSPSetFromOptions(solver);
// set A matrix with preconditioner
CHKERR KSPSetOperators(solver, A, A);
// set up the solver data strucure for the iterative solver
CHKERR KSPSetUp(solver);
// solve the system of linear equations
CHKERR KSPSolve(solver, f, x);
// destroy solver no needed any more
CHKERR KSPDestroy(&solver);
// make vector x available for parallel computations for visualization
// context
// save solution in vector x on mesh
// Set up post-processor. It is some generic implementation of finite
// element
PostProcVolumeOnRefinedMesh post_proc(m_field);
// Add operators to the elements, starting with some generic
CHKERR post_proc.generateReferenceElementMesh();
CHKERR post_proc.addFieldValuesPostProc("U");
CHKERR post_proc.addFieldValuesGradientPostProc("U");
dm, simple_interface->getDomainFEName().c_str(), &post_proc);
// write output
CHKERR post_proc.writeFile("out.h5m");
if (flg_test == PETSC_TRUE) {
const double x_vec_norm_const = 0.4;
// Check norm_1 value
double norm_check;
// Takes maximal element of the vector, that should be maximal
// displacement at the end of the bar
CHKERR VecNorm(x, NORM_INFINITY, &norm_check);
if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const > 1.e-10) {
"test failed (nrm 2 %6.4e)", norm_check);
// free memory handled by mofem discrete manager for A, x and f
CHKERR MatDestroy(&A);
CHKERR VecDestroy(&x);
CHKERR VecDestroy(&f);
// free memory allocated for mofem discrete manager
CHKERR DMDestroy(&dm);
// This is a good reference for the future
// finish work cleaning memory, getting statistics, etc
return 0;
static Index< 'm', 3 > m
static Index< 'p', 3 > p
static Index< 'K', 3 > K
const std::string default_options
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Catch errors.
Definition: definitions.h:441
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
@ H1
continuous field
Definition: definitions.h:177
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:552
Definition: definitions.h:217
Definition: definitions.h:132
Definition: definitions.h:123
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
Definition: elastic.cpp:34
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:982
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:425
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:478
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
Definition: DMMMoFEM.cpp:556
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:457
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
Definition: DMMMoFEM.cpp:597
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
use with loops to iterate row DOFs
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1129
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
double w(const double g, const double t)
constexpr auto VecSetValues
constexpr auto MatSetValues
const double r
rate factor
DataForcesAndSourcesCore::EntData EntData
int main(int argc, char *argv[])
static char help[]
Set integration rule to boundary elements.
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:77
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
Deprecated interface functions.
structure for User Loop Methods on finite elements
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
Post processing.
Set integration rule.