Hello world

Table of Contents

In this tutorial basic ideas about MoFEM::Simple interface and programming with user data operators (UDO) are presented.

MoFEM is designed for quick development of finite element application for multi-physics problems. We designed MoFEM to enable decomposition of the complex problem into small sub-problems, wherein each individual sub-problem can be independently tested and debugged or used in different context. That is achieved by pipelines of users data operators (UDO), see figure 1.

In MoFEM, there is a clear separation of declaration, definition and implementation of a field, finite element and problem. Such approach allows to use the same implementation in different problems, e.g. use the same implementation of an elastic finite element in a mechanical problem and thermo-mechanical problem without the need to introduce changes into a code and maximising code reusability. On the other hand for the same problem declaration and definition, one can test various formulations and implementations. In this example problem declaration and definition is managed by MoFEM::Simple interface, in which we focus attention only on field declaration and implementation of finite elements, in particular, UDO.

Figure 1. Users data operators

Hello world

* \file hello_world.cpp
* \ingroup mofem_simple_interface
* \example hello_world.cpp
* Prints basic information about users data operator.
* See more details in \ref hello_world_tut1
/* This file is part of MoFEM.
* 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/>. */
using namespace MoFEM;
static char help[] = "...\n\n";
static map<EntityType, std::string> type_name;
OpRow(const std::string &field_name)
: ForcesAndSourcesCore::UserDataOperator(field_name, field_name, OPROW) {}
MoFEMErrorCode doWork(int side, EntityType type,
if (type == MBVERTEX) {
// get number of evaluated element in the loop
std::cout << std::endl
<< "**** " << getNinTheLoop() << " **** " << std::endl;
std::cout << "**** Operators **** " << std::endl;
std::cout << "Hello Operator OpRow:"
<< " field name " << rowFieldName << " side " << side << " type "
<< type_name[type] << " nb dofs on entity "
<< data.getIndices().size() << std::endl;
OpRowCol(const std::string row_field, const std::string col_field,
const bool symm)
: ForcesAndSourcesCore::UserDataOperator(row_field, col_field, OPROWCOL,
symm) {}
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
std::cout << "Hello Operator OpRowCol:"
<< " row field name " << rowFieldName << " row side " << row_side
<< " row type " << type_name[row_type] << " nb dofs on row entity"
<< row_data.getIndices().size() << " : "
<< " col field name " << colFieldName << " col side " << col_side
<< " col type " << type_name[col_type] << " nb dofs on col entity"
<< col_data.getIndices().size() << std::endl;
OpVolume(const std::string &field_name)
field_name, OPROW) {
MoFEMErrorCode doWork(int side, EntityType type,
if (type == MBVERTEX) {
std::cout << "Hello Operator OpVolume:"
<< " volume " << getVolume() << endl;
OpFace(const std::string &field_name)
MoFEMErrorCode doWork(int side, EntityType type,
if (type == MBVERTEX) {
std::cout << "Hello Operator OpFace:"
<< " normal " << getNormal() << endl;
boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> &feSidePtr;
const std::string &field_name,
boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> &fe_side_ptr)
feSidePtr(fe_side_ptr) {}
MoFEMErrorCode doWork(int side, EntityType type,
if (type == MBVERTEX) {
std::cout << "Hello Operator OpSideFace" << endl;
CHKERR loopSideVolumes("dFE", *feSidePtr);
OpVolumeSide(const std::string &field_name)
field_name, field_name, OPROW) {}
MoFEMErrorCode doWork(int side, EntityType type,
if (type == MBVERTEX) {
std::cout << "Hello Operator OpVolumeSide:"
<< " volume " << getVolume() << " normal " << getNormal()
<< endl;
int main(int argc, char *argv[]) {
// initialize petsc
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
type_name[MBVERTEX] = "Vertex";
type_name[MBEDGE] = "Edge";
type_name[MBTRI] = "Triangle";
type_name[MBTET] = "Tetrahedra";
// Register DM Manager
DMType dm_name = "DMMOFEM";
// Create MoAB database
moab::Core moab_core;
moab::Interface &moab = moab_core;
// Create MoFEM database and link it to MoAB
MoFEM::Core mofem_core(moab);
MoFEM::Interface &m_field = mofem_core;
// Simple interface
Simple *simple_interface;
CHKERR m_field.getInterface(simple_interface);
// get options from command line
CHKERR simple_interface->getOptions();
// load mesh file
CHKERR simple_interface->loadFile();
// add fields
// set fields order
CHKERR simple_interface->setFieldOrder("U", 4);
CHKERR simple_interface->setFieldOrder("L", 3);
CHKERR simple_interface->setFieldOrder("S", 3);
// setup problem
CHKERR simple_interface->setUp();
// create elements instances
boost::shared_ptr<ForcesAndSourcesCore> domain_fe(
boost::shared_ptr<ForcesAndSourcesCore> boundary_fe(
boost::shared_ptr<ForcesAndSourcesCore> skeleton_fe(
boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> side_fe(
// set operator to the volume element instance
domain_fe->getOpPtrVector().push_back(new OpRow("U"));
domain_fe->getOpPtrVector().push_back(new OpRowCol("U", "U", true));
domain_fe->getOpPtrVector().push_back(new OpVolume("U"));
// set operator to the face element instance
boundary_fe->getOpPtrVector().push_back(new OpRow("L"));
boundary_fe->getOpPtrVector().push_back(new OpRowCol("U", "L", false));
boundary_fe->getOpPtrVector().push_back(new OpFace("L"));
// set operator to the face element on skeleton instance
skeleton_fe->getOpPtrVector().push_back(new OpRow("S"));
skeleton_fe->getOpPtrVector().push_back(new OpFaceSide("S", side_fe));
// set operator to the volume on side of the skeleton face
side_fe->getOpPtrVector().push_back(new OpVolumeSide("U"));
DM dm;
// get dm
CHKERR simple_interface->getDM(&dm);
// iterate domain elements and execute element instance with operator on
// mesh entities
// iterate boundary elements and execute element instance with operator on
// mesh entities
// iterate skeleton elements and execute element instance with operator on
// mesh entities
// destroy dm
CHKERR DMDestroy(&dm);
// finish work cleaning memory, getting statistics, etc.
return 0;

Code dissection


We initialize PETSc internal variables and register MoFEM discrete manager in PETSc. Next MoAB instance is created and assigned to it. Similarly, MoFEM instance is created and assigned. MoAB and MoFEM interfaces are abstract classes used to access data in database.

// Initialise PETSc
PetscInitialize(&argc,&argv,(char *)0,help);
// Register DM Manager
DMType dm_name = "DMMOFEM";
ierr = DMRegister_MoFEM(dm_name); CHKERRG(ierr);
// Create MoAB database
moab::Core moab_core;
moab::Interface& moab = moab_core;
// Create MoFEM database and link it to MoAB
MoFEM::Core mofem_core(moab);
MoFEM::Interface& m_field = mofem_core;

Next, we get access to database by MoFEM::Simple interface

Simple *simple_interface;
ierr = m_field.getInterface(simple_interface); CHKERRG(ierr);
// get options from command line
ierr = simple_interface->getOptions(); CHKERRG(ierr);
// load mesh file
ierr = simple_interface->loadFile(); CHKERRG(ierr);

and get options from the command line and load mesh file. Default mesh file has name mesh.h5m. Particular name of file can be given in command line using option -file_name my_file.h5m

The only indication that MoFEM database has been initialized is the following on the terminal

lib version 0.5.76
git commit id 253b1dd833ada49fe779e7f5fb60e52236ba51a9

where current MoFEM version and commit id is printed.

Adding fields

We add fields to the database. In MoFEM::Simple interface, three groups of fields can be declared. Fields defined in the domain, fields defined on boundary and fields defined on the skeleton. The same field can be defined on domain, boundary and skeleton. Fields are declared by giving its name, approximation space, base and number of coefficients. See for details here MoFEM::Simple::addDomainField.

ierr = simple_interface->addDomainField("U",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRG(ierr);
ierr = simple_interface->addBoundaryField("L",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRG(ierr);
ierr = simple_interface->addSkeletonField("S",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRG(ierr);

Next, we set approximation order to those fields

// set fields order
ierr = simple_interface->setFieldOrder("U",4); CHKERRG(ierr);
ierr = simple_interface->setFieldOrder("L",3); CHKERRG(ierr);
ierr = simple_interface->setFieldOrder("S",3); CHKERRG(ierr);

For more details see MoFEM::Simple::setFieldOrder. If needed function MoFEM::Simple::setFieldOrder can be supplemented by the additional parameter to set order to particular field entities, enabling heterogeneous approximation base.

As result the following is printed on terminal

add: name U BitFieldId 1 bit number 1 space H1 approximation base AINSWORTH_LEGENDRE_BASE rank 3 meshset 12682136550675316768
add: name L BitFieldId 2 bit number 2 space H1 approximation base AINSWORTH_LEGENDRE_BASE rank 3 meshset 12682136550675316769
add: name S BitFieldId 4 bit number 3 space H1 approximation base AINSWORTH_LEGENDRE_BASE rank 3 meshset 12682136550675316770

where at the end entity handle to meshset is printed. This meshset consisting all entities to which field is set. The effect of the setting of approximation order will be viable when fields are constructed during set-up stage.


The fields, finite elements and problems and all other data structures are build with the following code

ierr = simple_interface->setUp(); CHKERRG(ierr);

As results following lines are printed on the terminal

add finite element: dFE
add finite element: bFE
add finite element: sFE
add problem: SimpleProblem

Allocating finite element instances

Once the problem is defined and set up, we create finite element instances. We are using generic implementation of finite elements which is tailored to the dimension of finite element entities. Note that in the following code we use MoFEM::VolumeElementForcesAndSourcesCore for domain elements and MoFEM::FaceElementForcesAndSourcesCore for boundary element and skeleton elements, since we use 3d mesh, which boundary and skeleton are 2d entities emerged in 3d space

boost::shared_ptr<ForcesAndSourcesCore> domain_fe(new VolumeElementForcesAndSourcesCore(m_field));
boost::shared_ptr<ForcesAndSourcesCore> boundary_fe(new FaceElementForcesAndSourcesCore(m_field));
boost::shared_ptr<ForcesAndSourcesCore> skeleton_fe(new FaceElementForcesAndSourcesCore(m_field));

Once finite elements objects are created we add operators to them according to figure 1. Starting from domain_fe

domain_fe->getOpPtrVector().push_back(new OpRow("U"));
domain_fe->getOpPtrVector().push_back(new OpRowCol("U","U",true));
domain_fe->getOpPtrVector().push_back(new OpVolume("U"));

next boundary_fe

boundary_fe->getOpPtrVector().push_back(new OpRow("L"));
boundary_fe->getOpPtrVector().push_back(new OpRowCol("U","L",false));
boundary_fe->getOpPtrVector().push_back(new OpFace("L"));

and finally skeleton_fe

skeleton_fe->getOpPtrVector().push_back(new OpRow("S"));

To be more explicit, we can write for example for domain_fe, we have

boost::ptr_vector<UserDataOperator>& sequence_of_operators = domain_fe->getOpPtrVector();
sequence_of_operators.push_back(new OpRow("U"));
sequence_of_operators.push_back(new OpRowCol("U","U",true));
sequence_of_operators.push_back(new OpVolume("U"));

Note that boost::ptr_vector (see link) is container of pointers, once vector is destroyed memory pointed by pointers is released.

We have to discuss integration over the skeleton. This functionality is used for a class of Petrov-Discontinuous Galerkin methods, or when Nitsche method is applied. Those methods involve of integration of traces of approximation functions on faces. That enforces the evaluation of derivatives of base functions on finite elements adjacent to the face. In MoFEM that is realised in a way, that generic element instance is iterated over adjacent entities to faces. We start with allocation memory for "side" finite element as before code boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> side_fe(new VolumeElementForcesAndSourcesCoreOnSide(m_field)); Note that this time finite element class is of type MoFEM::VolumeElementForcesAndSourcesCoreOnSide, since this element has integration points which are associated with face adjacent face. Next, to that element we add user data operators, in this particular case only one

side_fe->getOpPtrVector().push_back(new OpVolumeSide("U"));

Note that pointer to finite element instance is passed in the constructor of user data operator for skeleton element, what is discussed in detail in next section. Finally, we add operator to skeleton face element itself

skeleton_fe->getOpPtrVector().push_back(new OpFaceSide("S",side_fe));

Discrete manager

Following line gives access to the discrete manager (DM). DM is the interface developed for PETSc to manage complexities associated with topology (mesh) and algebra and several implementation of this interface are available in native PETSc library. Here we are using the particular implementation developed for MoFEM. To get interface the following lines of the code are needed

DM dm;
ierr = simple_interface->getDM(&dm); CHKERRG(ierr);

When DM interface is no longer needed, it needs to be destroyed by user in the standard way well-known from PETSc interface

ierr = DMDestroy(&dm); CHKERRG(ierr);
In MoFEM implementation, in comparison to other finite element codes, we do not wrap PETSc library functions but extend PETSc interface with features native to MoFEM. You can use DM interface in the same way how you use native PETSc interface, see Solving the Poisson equation it is shown how to set up KSP solver and get the solution. That philosophy applies to MoAB functions as well.

Iterating finite elements

Finally, we get to the point when we can put our machine in motion, we iterate over finite elements and run sequences of users data operator for each of them. We start with domain_fe

ierr = DMoFEMLoopFiniteElements(dm,simple_interface->getDomainFEName(),domain_fe); CHKERRG(ierr);

as results we get

**** 0 ****
**** Operators ****
Hello Operator OpRow: field name U side 0 type Vertex nb dofs on entity 12
Hello Operator OpRow: field name U side 0 type Tetrahedra nb dofs on entity 3
Hello Operator OpRowCol: row field name U row side 0 row type Vertex nb dofs on row entity12 :  col field name U col side 0 col type Vertex nb dofs on col entity12
Hello Operator OpRowCol: row field name U row side 0 row type Tetrahedra nb dofs on row entity3 :  col field name U col side 0 col type Tetrahedra nb dofs on col entity3
Hello Operator OpVolume: volume 0.0834851

**** 2 ****
**** Operators ****

**** 11 ****
**** Operators ****

where dots are added for the abbreviation of output. Note that operators are called in the order we pushed them to finite element operators vector. Since we have twelve volume (Tetrahedra) elements, iteration ends on eleven, since in MoFEM we always start counting from zero.

To iterate over boundary element, analogically we do,

ierr = DMoFEMLoopFiniteElements(dm,simple_interface->getBoundaryFEName(),boundary_fe); CHKERRG(ierr);

and we get similar output to the one shown before, with one difference being that the last operator does not print volume of the element but that is normal since entity of boundary finite element in this particular case is the triangle.

The same procedure is applied to iterate over skeleton finite elements entities

ierr = DMoFEMLoopFiniteElements(dm,simple_interface->getSkeletonFEName(),skeleton_fe); CHKERRG(ierr);

We have thirty skeleton elements and output looks as follows

**** 0 ****
**** Operators ****
Hello Operator OpRow: field name S side 0 type Vertex nb dofs on entity 9
Hello Operator OpRow: field name S side 0 type Edge nb dofs on entity 6
Hello Operator OpRow: field name S side 1 type Edge nb dofs on entity 6
Hello Operator OpRow: field name S side 2 type Edge nb dofs on entity 6
Hello Operator OpRow: field name S side 0 type Triangle nb dofs on entity 3
Hello Operator OpSideFace
Hello Operator OpVolumeSide: volume 0.0782402 normal [3](0,0,1)

**** 13 ****
**** Operators ****
Hello Operator OpRow: field name S side 0 type Vertex nb dofs on entity 9
Hello Operator OpRow: field name S side 0 type Edge nb dofs on entity 6
Hello Operator OpRow: field name S side 1 type Edge nb dofs on entity 0
Hello Operator OpRow: field name S side 2 type Edge nb dofs on entity 0
Hello Operator OpRow: field name S side 0 type Triangle nb dofs on entity 0
Hello Operator OpSideFace
Hello Operator OpVolumeSide: volume 0.0834851 normal [3](0,0.530559,-0.50091)
Hello Operator OpVolumeSide: volume 0.0884264 normal [3](0,0.530559,-0.50091)

Note that first operator is OpRow, the second operator is OpFaceSide, this operator prints its name and runs integration over adjacent to given face elements, which is side_fe. Once this element is run for each adjacent finite element entity, users data operators are run on it, i.e. OpVolumeSide which prints volume of the adjacent entity and normal of the face. Note that first element has only one run of OpVolumeSide, since skeleton finite element "0" is on the body boundary, while skeleton finite element "13" is in body volume and it has two volume neighbours.

Implementation of users data operators

Now we focus attention on the implementation of users data operator. The first operator has the structure

OpRow(const std::string& field_name):
ForcesAndSourcesCore::UserDataOperator(field_name,field_name,OPROW) {
MoFEMErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data);

This users data operator class is derived from MoFEM::ForcesAndSourcesCore::UserDataOperator which can be used with any type of entity. It is the type of OPROW, which indicates that it only iterates lower dimension entities on the element. On each lower entity overload method is called

MoFEMErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data);

which as arguments takes entity side number (local entity number on the finite element), entity type (e.g. MBVERTEX, MBEDGE, MBTET) and reference to structure MoFEM::DataForcesAndSourcesCore::EntData, which keeps information on DOFs, approximation on given entity. This type of entity is usually used to integrate the right-hand side vector.

Another type of users data operator is implemented here

const std::string row_field,
const std::string col_field,
const bool symm
ForcesAndSourcesCore::UserDataOperator(row_field,col_field,OPROWCOL,symm) {
int row_side,int col_side,
EntityType row_type,EntityType col_type,

This user data operator is of type OPROWCOL, which takes an additional parameter in constructor, i.e. symm, which is used to set symmetry of operator. Operator of this type iterates over all unique pairs of entities. If a symmetric operator pair is set of two elements (i.e. entities), thus order of entities is not important. If an operator is not symmetric, then pairs are the sequence of two elements and all variations of entities pairs are considered. This type of operator is used to integrate matrices. Note that this time function is overloaded, which takes as argument data for rows and columns, respectively.

Performing calculations on entity of specific dimension additional data like volume, normal need to be attained, for such case derived users data operator class can be used, e.g.

and for case of operator working on adjacent to face volume entity

where from members of this class information about face normal and adjacent entity volume can be accessed.

Running the program

The easiest installation is with Docker containers; please see youtube video how to run installation link. You can also look here Installation with Docker (Linux, Mac OS X and some versions of Windows). Open terminal and run Docker container

docker run --rm=true -it --volumes-from mofem_build -v $HOME/mofem-cephas/mofem:/mofem -v $HOME:$HOME -e HOSTHOME=$HOME mofem_build /bin/bash

You are now in Docker container, go to directory where Poisson problem is located, and compile code

cd /mofem_build/um/basic_finite_elements/hello_world/
make -j2

It can take some time when you compile it for the first time since all essential functionality is compiled

Once code is compiled, we can run test