v0.13.1
MoFEM Architecture

# 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, the 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 elements to create a new capability for solving complex problems efficiently. Approximation space defines adjacency of DOFs on entities, the number of DOFs on an entity is independent on approximation base (at least for polynomials). The base on the entity is a trace of the base on the element, and opposite relation works, base on the entity is propagated to the 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 the 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.

# Basic design description

## MoFEM software ecosystem

MoFEM has not developed and will not develop all capabilities from scratch. Instead, it integrates advanced scientific computing tools for Algebra (PETSc and associated tools), Topology (MOAB and associated tools) and data structures from Boost. The resilience of this ecosystem is ensured by solid fundation of its components and by a dynamic and established group of developers backed by a constantly growing user base. MoFEM focuses its attention on complexities related to finite element technology and uses abstractions like field entity, DOF (degree of freedom), finite element and problem.

Ecosystem

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 the 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 with solvers and finally MoFEM core library managing complexities directly related to the finite element method. Each part has its design objectives and appropriate programming tools from a spectrum of solutions can be selected [34]. MoFEM makes PETSc an integral part of the code by extending PETSc via DMMOFEM interface (several other functions work directly on PETSc objects). MoAB on the other hand is adopted for internal data storage. PETSc is designed to be fast, whereas MoAB aims to be memory efficient.

MoFEM code parts

Code development, with no difference 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 implementing finite elements than the time needed to run the 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, the 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 elements for a given predefined type, but UserDataOperators operate on the finite element entities. In such a method, the order of approximation and the shape of the finite element does not influence how the user implements a problem.

Moreover, unlike in common finite element codes, 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 a given field number of DOFs or base functions on an assigned 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 elements and entities. Elements and entities share data structures avoiding unnecessary memory allocation and de-allocation. However, the motivation for this development is code elasticity (resilience to change and openness for new finite element technologies) and reusability. With relatively small code changes, large group of problems can be implemented and last but not least this approach 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;
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator("U","U",UserDataOperator::OPROWCOL,true),
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
entityLocalMatrix.resize(nb_row_dofs,nb_col_dofs,false);
// zero elements
entityLocalMatrix.clear();
// 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);
++t_base_diff_col;
++t_local_mat;
}
++t_base_diff_row;
}
}
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();
A,
nb_row,&*row_ind.begin(),
nb_col,&*col_ind.begin(),
);
}
// 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);
}
}
};
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
EntitiesFieldData::EntData EntData
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#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
FTensor::Index< 'i', SPACE_DIM > i
Definition: heat_method.cpp:38
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
double w(const double g, const double t)
double A
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element

## Advanced design of Core library and application code

Typically, internal data structures are hidden from the user and there is no need to access them directly. In the 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 the 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();
// 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();
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ H1
continuous field
Definition: definitions.h:85
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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