Basic design

Idea behind MoFEM

Traditional finite element codes are element-centric (type of element defines approximation space and base) and therefore cannot exploit the potential of emerging approximation methods. MoFEM software utilises recent advances in finite element technology and modern data structures, enabling the efficient solution of difficult, multi-domain, multi-scale and multi-physics problems. In this code, design of data structures for approximation of field variables are independent of the specific finite element (e.g. Lagrangian, Nedelec, Rivart-Thomas) being used, such that finite element is constructed by a set of lower dimension entities on which the approximation fields are defined.

Therefore, different approximation spaces (H1, H-curl, H-div, L2) can be arbitrarily mixed in finite element to create new capability for solving complex problems efficiently. Approximation space defines adjacency of DOFs on entities, the number of DOFs on entity is independent on approximation base (at least for polynomials). The base on entity is a trace of the base on element, and opposite relation works, base on entity is extruded into element.

MoFEM data structures allow easily to enrich approximation fields or modify base functions, e.g. to resolve singularity at crack front. Applying this technology, it is effortless to construct transition elements between domain with different problem formulation and physics (e.g. from two-field mix-formulation to single-field formulation), or elements with anisotropic approximation order (e.g. in solid-shells with arbitrary high-order on shell surface and arbitrary low-order through thickness).

This approach also sets the benchmark in terms of how finite element codes are implemented, introducing a concept of user-defined data operators acting on fields that are associated with entities (vertices, edges, faces and volumes) rather on the finite element directly. Such an approach simplifies code writing, testing and validation, making the code resilient to bugs.

Before you start to look into examples and tutorials

MoFEM software eccosystem

MoFEM has not developed and will not develop all capabilities from scratch. Instead, MoFEM integrates advanced scientific computing tools for Algebra (PETSc and associated tools), Topology (MOAB and associated tools) and data structures from Boost. Resilience of this ecosystem is ensured because the underpinning components have sustainable funding, dynamic and established group of developers and significant user base. MoFEM focuses attention on complexities related to finite element technology and uses abstractions like field entity, DOF (degree of freedom), finite element and problem.


Finite element software is a complex ecosystem managing mesh and topology related complexities, sparse algebra and complications related to approximation, integration or dense tensor algebra at integration point level. Each element by itself is a standalone complicated problem. MoFEM ecosystem is composed of three parts, MoAB managing complexities related to topology, PETSc managing complexities related to sparse algebra and solvers and finally MoFEM core library managing complexities directly related to finite element method. Each part has its own design objectives and appropriate programming tools from a spectrum of solutions can be selected [32]. MoFEM makes PETSc integral part of code by extending PETSc by DMMOFEM interface (several other functions work directly on PETSc objects). MoAB from other hand is internal data storage. PETSc is designed to be fast, whereas MoAB aims to be memory efficient.

MoFEM code parts

Code development, with no different from engineering, is a compromise of different goals. The first primary design objective of MoFEM is to develop flexible and extendable code that is easy to test and bug search. Very often, we spend more time to implement finite element than the time needed to run analysis itself, flexibility and extensibility allow for rapid development of finite elements for non-standard and complex problems. Utilising PETSc and MoAB scalability is the second important design objective. Nowadays, we can have easy access to the computing clusters through cloud computing or local/regional supercomputers. Third, core library cannot compromise runtime efficiency and use of computer memory with care.

Finite element implementation and UserDataOperator

The finite element entity is composed of subentities, like nodes, edges, faces (triangle, quad, etc.) and volumes (tetrahedron, hex, prism, etc). The element is implemented by evaluating base functions, indices, data on each entity in the domain of the element. The developer does not directly implement finite element for a given predefined type, but UserDataOperators operate on finite element entities. In such method order of approximation, the shape of the finite element does not influence how user implements problem.

Moreover, unlike in common finite element code, the choice of field space and base space is independent of the type of element, it happens before we select the domain of integration and finite elements on it. Concepts like the periodic table of finite elements, or zoo of finite elements are obsolete.

We use the concept of a finite element as a domain on which integration rule is constructed and base functions are evaluated. The finite element has an auxiliary set of data, for example, a measure of the domain, i.e. volume, area or length. Some elements carry information about normal, direction depending on the dimension of the finite element domain. The base finite element interface can be seen here MoFEM::ForcesAndSourcesCore. Note that most of the derived classes do not alter finite element core behaviour, add only type-specific functionality, f.e. appropriate rule of integration for given shape and dimension of element entity. Except for special cases, differential operators could be implemented unrespectfully on the type of finite element entity. The structure and finite element interaction with user data operators are shown in the figure below

Finite element and its structure

The fundamental role in finite element implementation plays MoFEM::ForcesAndSourcesCore::UserDataOperator. The generality of this technology range of applications 0d-2d or 3d, with Nitsche's method or with H-div or H-curl spaces. The UserDataOperator are operated in three typical ways, first (OPROW/OPCOL) are usually used to evaluate the right-hand side vectors, second (OPROWCOL) to evaluate the left-hand side matrices, and third, often use to modify base functions of given space. First two are intrinsically linked to operations on particular fields, the last one is abstract and should not be used to assemble vectors and matrices. Type of operator is set by overloading function doWork. Operators are run in sequence as they are added to finite element and executed for each entity on that element. If for given field number of DOFs or base functions on given entity is zero, such entity is not evaluated.

Types of UserDataOperator

Note that the data and the algorithm operating on that data are localised and grouped by both element and entity. Elements and entities share data structure avoiding unnecessary memory allocation and deallocation. However, the motivation for this development is code elasticity (resilience to change and openness for new finite element technologies), reusability, which with relatively small code, large group of problems can be implemented and last but not least enables modularisation and testing of each element independently.

Example of implementation

Above example reflects how an engineer would formulate the problem for linear thermoelasticity, calculate physical quantities like strain, stress, temperature or flux and assemble a system of linear equations to solve it staggered method. Another way approaching this problem to think regarding differential operators. One can implement set of operators,

\[ (\nabla u ,\nabla u),\;(\nabla \cdot \tau, u),\;(u,u) \]

and example, assembly of

\[ \mathbf{A} = (\nabla u,\nabla u)_\Omega = \int_\Omega \nabla \mathbf{u} \cdot \nabla \mathbf{v} \textrm{d}\Omega \]

take form

Mat A;
OpGradGrad(Mat a):
A(a) {
MatrixDouble entityLocalMatrix;
* \brief integrate matrix
* @param row_side local number of entity on element for row of the matrix
* @param col_side local number of entity on element for col of the matrix
* @param row_type type of row entity (VERTICES/EDGE0-5/TRIANGLE0-3/TETRAHEDRON)
* @param col_type type of col entity (VERTICES/EDGE0-5/TRIANGLE0-3/TETRAHEDRON)
* @param row_data structure of data, like base functions and associated methods to access those data on rows
* @param col_data structure of data, like base functions and associated methods to access those data on rows
* @return error code
int row_side,int col_side,
EntityType row_type,EntityType col_type,
) {
// get number of dofs on rows
const int nb_row_dofs = row_data.getIndices().size();
// if no dofs, end here
if(nb_row_dofs==0) MoFEMFunctionReturnHot(0);
// get number of dofs on column
const int nb_col_dofs = col_data.getIndices().size();
// if no dofs, end here
if(nb_col_dofs==0) MoFEMFunctionReturnHot(0);
// resize entity local matrix
// zero elements
// get number of integration points
const int nb_gauss_pts = row_data.getN().size1();
// integrate
for(int gg = 0;gg!=nb_gauss_pts;gg++) {
// get integration weight scaled by volume
double w = getGaussPts()(3,gg)*getVolume();
// get derivatives of base functions on rows
auto t_base_diff_row = row_data.getFTensor1DiffN<3>();
// take first element in row "aa"
FTensor::Tensor0<double*> t_local_mat(&entityLocalMatrix(0,0),1);
// loop rows
for(int aa = 0;aa!=nb_row_dofs;aa++) {
// get derivatives of base functions on columns
auto t_base_diff_col = col_data.getFTensor1DiffN<3>(gg,0);
// loop column
for(int bb = 0;bb!=nb_col_dofs;bb++) {
t_local_mat += w*t_base_diff_row(i)*t_base_diff_col(i);
auto assemble = [A](
auto & row_ind, auto & col_ind, auto & data) {
const int nb_row = row_ind.size();
const int nb_col = col_ind.size();
// assemble global matrix
CHKERR assemble(
row_data.getIndices(), col_data.getIndices(), entityLocalMatrix);
// and its transpose
if(row_side != col_side || row_type != col_type) {
entityLocalMatrix = trans(entityLocalMatrix);
CHKERR assemble(
col_data.getIndices(), row_data.getIndices(), entityLocalMatrix);

Core library and application code

Typically, internal data structures are hidden from the user and there is no need to access them directly. In core library, no single problem or finite element in the classical sense is implemented. Core library only provides functionality for finite element implementation and application code development. For that reason, repositories of core library and user modules are in separate independent locations.

MoFEM code parts

In the figure above on first instant, a developer is interested in the implementation of user data operators and tasks in the yellow boxes. It has to define fields, which span on some entities and set approximation order. For example

MoFEM::Core core(moab);
MoFEM::Interface& m_field = core;
// define fields
// meshset consisting all entities in mesh
EntityHandle root_set = moab.get_root_set();
// add entities to field
CHKERR m_field.add_ents_to_field_by_TETs(root_set,"DISP");
CHKERR m_field.add_ents_to_field_by_TETs(root_set,"FLUX");
// set app. order (that is boring same order to all)
int order = 5;
CHKERR m_field.set_field_order(root_set,MBTET,"DISP",order);
CHKERR m_field.set_field_order(root_set,MBTET,"FLUX",order);
// build
CHKERR m_field.build_fields();

Similarly, defining finite element set implement integration domain and build internal data structures. Note, finite element structure is separated from implementation. The same element can have several implementations which work on some set of user data operators. Moreover, the same problem, defined on the same field and data structures can be implemented in separate modules. Also, several users modules exchange data by common MoFEM database, without knowledge how other modules have been implemented, that enables simplified cooperation between developers.

Internal field structure

You don't have to know how internal data structures work, like the driver of the car, does not necessarily know how a spark plug works. So if you get bored, you can skip this part.

Fields data structures

However, sometimes it is useful to know details when something goes wrong. The objective is to have an extendable and flexible data structure which can manage a wide range of approximations. At the same time, it should be easy and fast access to data structures, short set-up time and memory efficiency without redundancies. This aim is realised by use boost multi-indices, aliased shared pointers and vector sequences. See for example MoFEM::NumeredDofEntity or MoFEM::DofEntity. Data are accessed by multi-indices, for example, defined in DofsMultiIndices.hpp, FieldEntsMultiIndices.hpp or FieldMultiIndices.hpp. The similar data structures are created for finite elements and problems.

#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
@ H1
continuous field
Definition: definitions.h:177
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
Definition: VolumeElementForcesAndSourcesCore.hpp:354
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
double w(const double g, const double t)
Definition: ContactOperators.hpp:141
DataForcesAndSourcesCore::EntData EntData
Definition: quad_polynomial_approximation.cpp:27
#define CHKERR
Inline error check.
Definition: definitions.h:604
FTensor::Index< 'i', 3 >
const int order
Definition: rd_stdOperators.hpp:28
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Definition: Tensor0.hpp:17
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: DataOperators.hpp:45
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:77
Definition: definitions.h:158
double getVolume() const
element volume (linear geometry)
Definition: VolumeElementForcesAndSourcesCore.hpp:403
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
Mat A
Definition: LoopMethods.hpp:73
DEPRECATED MoFEMErrorCode add_ents_to_field_by_TETs(const EntityHandle meshset, const std::string &name, int verb=-1)
set field entities from adjacencies of tetrahedron
Definition: DeprecatedCoreInterface.cpp:530
default operator for TET element
Definition: VolumeElementForcesAndSourcesCore.hpp:45
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: DataStructures.hpp:1599
FTensor::Index< 'i', 3 > i
Definition: matrix_function.cpp:18
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
field with continuous normal traction
Definition: definitions.h:179
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1016