FUN-0: Hello world

Table of Contents

• Hello world
• Code dissection
  • Initialization
  • Adding fields
  • Set-up
  • Allocating finite element instances
  • Iterating finite elements
  • Implementation of user data operators
• Running the program
• Exercises

In this tutorial, basic ideas about MoFEM::Simple interface, MoFEM::PipelineManager, and programming with user data operators (UDO) are presented. This tutorial aims at providing a big picture of the sequence of how MoFEM works in a simple and typical problem. Before proceeding to the next parts, readers are strongly encouraged to look at the MoFEM Architecture to have fundamental ideas of MoFEM and its ecosystem.

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 user data operators (UDO), see Figure 1. In this figure, there are operators in green, finite element instances which are presented in blue, and red shows finite element objects. Typical operators including OpRow, OpRowCol and OpVolume are mainly used to construct the force vector (the right hand side) and the stiffness matrix and to assemble the system, respectively. Meanwhile, the finite element instances include domain_fe, boundary_fe, and skeleton_fe. Operators can be set to volume or face elements on one of the finite element instances.

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. User 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
 *
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
static char help[] = "...\n\n";
 
#define HelloFunctionBegin                                                     \
  MoFEMFunctionBegin;                                                          \
  MOFEM_LOG_CHANNEL("SYNC");                                                   \
  MOFEM_LOG_FUNCTION();                                                        \
  MOFEM_LOG_TAG("WORLD", "HelloWorld");
 
struct OpRow : public ForcesAndSourcesCore::UserDataOperator {
  OpRow(const std::string &field_name)
      : ForcesAndSourcesCore::UserDataOperator(field_name, field_name, OPROW) {}
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
    HelloFunctionBegin;
    if (type == MBVERTEX) {
      // get number of evaluated element in the loop
      MOFEM_LOG("SYNC", Sev::inform) << "**** " << getNinTheLoop() << " ****";
      MOFEM_LOG("SYNC", Sev::inform) << "**** Operators ****";
    }
    MOFEM_LOG("SYNC", Sev::inform)
        << "Hello Operator OpRow:"
        << " field name " << rowFieldName << " side " << side << " type "
        << CN::EntityTypeName(type) << " nb dofs on entity "
        << data.getIndices().size();
    MoFEMFunctionReturn(0);
  }
};
 
struct OpRowCol : public ForcesAndSourcesCore::UserDataOperator {
  OpRowCol(const std::string row_field, const std::string col_field,
           const bool symm)
      : ForcesAndSourcesCore::UserDataOperator(row_field, col_field, OPROWCOL,
                                               symm) {}
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        EntitiesFieldData::EntData &row_data,
                        EntitiesFieldData::EntData &col_data) {
    HelloFunctionBegin;
    MOFEM_LOG("SYNC", Sev::inform)
        << "Hello Operator OpRowCol:"
        << " row field name " << rowFieldName << " row side " << row_side
        << " row type " << CN::EntityTypeName(row_type)
        << " nb dofs on row entity " << row_data.getIndices().size() << " : "
        << " col field name " << colFieldName << " col side " << col_side
        << " col type " << CN::EntityTypeName(col_type)
        << " nb dofs on col entity " << col_data.getIndices().size();
    MoFEMFunctionReturn(0);
  }
};
 
struct OpVolume : public VolumeElementForcesAndSourcesCore::UserDataOperator {
  OpVolume(const std::string &field_name)
      : VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,
                                                            field_name, OPROW) {
  }
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
    HelloFunctionBegin;
    if (type == MBVERTEX) {
      MOFEM_LOG("SYNC", Sev::inform) << "Hello Operator OpVolume:"
                                     << " volume " << getVolume();
    }
    MoFEMFunctionReturn(0);
  }
};
 
struct OpFace : public FaceElementForcesAndSourcesCore::UserDataOperator {
  OpFace(const std::string &field_name)
      : FaceElementForcesAndSourcesCore::UserDataOperator(field_name, OPROW) {}
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
    HelloFunctionBegin;
    if (type == MBVERTEX) {
      MOFEM_LOG("SYNC", Sev::inform) << "Hello Operator OpFace:"
                                     << " normal " << getNormal();
    }
    MoFEMFunctionReturn(0);
  }
};
 
struct OpFaceSide : public FaceElementForcesAndSourcesCore::UserDataOperator {
  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> &feSidePtr;
  OpFaceSide(
      const std::string &field_name,
      boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> &fe_side_ptr)
      : FaceElementForcesAndSourcesCore::UserDataOperator(field_name, OPROW),
        feSidePtr(fe_side_ptr) {}
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
 
    HelloFunctionBegin;
    if (type == MBVERTEX) {
      MOFEM_LOG("SYNC", Sev::inform) << "Hello Operator OpSideFace";
      CHKERR loopSideVolumes("dFE", *feSidePtr);
    }
    MoFEMFunctionReturn(0);
  }
};
 
struct OpVolumeSide
    : public VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator {
  OpVolumeSide(const std::string &field_name)
      : VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator(
            field_name, field_name, OPROW) {}
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
    HelloFunctionBegin;
    if (type == MBVERTEX) {
      MOFEM_LOG("SYNC", Sev::inform)
          << "Hello Operator OpVolumeSide:"
          << " volume " << getVolume() << " normal " << getNormal();
    }
    MoFEMFunctionReturn(0);
  }
};
 
int main(int argc, char *argv[]) {
 
  // initialize petsc
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    // Register DM Manager
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    // 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;
    CHKERR m_field.getInterface(simple);
 
    // get options from command line
    CHKERR simple->getOptions();
    // load mesh file
    CHKERR simple->loadFile();
    // add fields
    CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 3);
    CHKERR simple->addBoundaryField("L", H1, AINSWORTH_LEGENDRE_BASE, 3);
    CHKERR simple->addSkeletonField("S", H1, AINSWORTH_LEGENDRE_BASE, 3);
    // set fields order
    CHKERR simple->setFieldOrder("U", 4);
    CHKERR simple->setFieldOrder("L", 3);
    CHKERR simple->setFieldOrder("S", 3);
 
    // setup problem
    CHKERR simple->setUp();
 
    auto pipeline_mng = m_field.getInterface<PipelineManager>();
 
    //! [set operator to the volume element instance]
    pipeline_mng->getOpDomainRhsPipeline().push_back(new OpRow("U"));
    pipeline_mng->getOpDomainRhsPipeline().push_back(new OpVolume("U"));
    pipeline_mng->getOpDomainLhsPipeline().push_back(
        new OpRowCol("U", "U", true));
    //! [set operator to the volume element instance]
 
    //! [set operator to the face element instance]
    pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpRow("L"));
    pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpFace("L"));
    pipeline_mng->getOpBoundaryLhsPipeline().push_back(
        new OpRowCol("U", "L", false));
    //! [set operator to the face element instance]
 
    //! [create skeleton side element and push operator to it]
    boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> side_fe(
        new VolumeElementForcesAndSourcesCoreOnSide(m_field));
    side_fe->getOpPtrVector().push_back(new OpVolumeSide("U"));
    //! [create skeleton side element and push operator to it]
 
    //! [set operator to the face element on skeleton instance]
    pipeline_mng->getOpSkeletonRhsPipeline().push_back(new OpRow("S"));
    pipeline_mng->getOpSkeletonRhsPipeline().push_back(
        new OpFaceSide("S", side_fe));
    //! [set operator to the face element on skeleton instance]
 
    //! [executing finite elements]
    CHKERR pipeline_mng->loopFiniteElements();
    //! [executing finite elements] 
 
    MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
  }
  CATCH_ERRORS;
 
  // finish work cleaning memory, getting statistics, etc.
  MoFEM::Core::Finalize();
 
  return 0;
}

Code dissection

While the details of different parts of the code are discussed later in this section, the main structure of the code is as follows:

1. Initialization: PETSc, MoFEM, and MoAB are initialized
2. Adding fields: Three groups of fields are added including fields defined in the domain, fields defined on boundary, and fields defined on the skeleton
3. Setting-up: Building of fields, finite elements and problems, and all other data structures
4. Allocating finite element instances: Finite element instances are created including domain_fe, boundary_fe, and skeleton_fe
5. Accessing discrete manager: Get the discrete manager which is the interface developed for PETSc to manage complexities associated with topology (mesh) and algebra
6. Iterating finite elements: Make iteration over finite elements and run sequences of UDO for each of the element. For demonstation, operator names and entity information are printed out in this part
7. Implementation of UDO: Define different user data operators

Initialization

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";
CHKERR DMRegister_MoFEM(dm_name); 
// 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;
CHKERR m_field.getInterface(simple); 
// get options from command line
CHKERR simple->getOptions(); 
// load mesh file
CHKERR simple->loadFile(); 

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 visualization of the meshed object is given in Figure 2.

The indication that MoFEM database has been initialized is that the MoFEM version and git commit id are shown on the terminal as below.

[0] <inform> MoFEM version 0.10.0 (MOAB 5.1.0 Petsc Release Version 3.11.3, Jun, 26, 2019 )
[0] <inform> git commit id 8e79d6cc5f02e48a97287c6663c53b1e871766f4
[0] <inform> Local time: 2020-Sep-17 15:55:3
[0] <inform> UTC time: 2020-Sep-17 14:55:3

It should be noted that the message can be different depending on the time of installation and update.

Figure 2. Meshed object

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.

CHKERR simple->addDomainField("U",H1,AINSWORTH_LEGENDRE_BASE,3); 
CHKERR simple->addBoundaryField("L",H1,AINSWORTH_LEGENDRE_BASE,3); 
CHKERR simple->addSkeletonField("S",H1,AINSWORTH_LEGENDRE_BASE,3); 

Next, we set approximation order to those fields

// set fields order
CHKERR simple->setFieldOrder("U",4); 
CHKERR simple->setFieldOrder("L",3); 
CHKERR simple->setFieldOrder("S",3); 

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. FUN-2: Hierarchical approximation provides more details on the different entities used in MoFEM.

When one runs the program later in Running the program, the output of this part of the code is printed on terminal as follows

[0] <inform> [FieldCore] Add field U field_id 1 space H1 approximation base AINSWORTH_LEGENDRE_BASE rank 3 meshset 12682136550675316767
[0] <inform> [FieldCore] Add field L field_id 2 space H1 approximation base AINSWORTH_LEGENDRE_BASE rank 3 meshset 12682136550675316768
[0] <inform> [FieldCore] Add field S field_id 4 space H1 approximation base AINSWORTH_LEGENDRE_BASE rank 3 meshset 12682136550675316769

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.

Set-up

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

CHKERR simple->setUp(); 

Similarly, when one runs the code, the following outputs are expected on the terminal

• Three finite elements structures are added to the database including domain finite element structure, boundary finite element structure and skeleton boundary element structure. Those elements are added by default via MoFEM::Simple interface.

add finite element: dFE
add finite element: bFE
add finite element: sFE

• Next problem is added to the database. This is the default name for MoFEM::Simple interface

add problem: SimpleProblem

• In the following lines, degrees of freedom (DOFs) are constructed for fields U, L and S added in Adding fields. The number of DOFs on each field entity depends on the approximation order and approximation space which have also been defined in Adding fields.

  [0] <verbose> [FieldCore] Build field U
  [0] <verbose> [FieldCore] Nb. of dofs (vertices) 27 (inactive 0)
  [0] <verbose> [FieldCore] Nb. of dofs (edge) 234 (inactive 0)
  [0] <verbose> [FieldCore] Nb. of dofs (triangles) 270 (inactive 0)
  [0] <verbose> [FieldCore] Nb. of dofs (tetrahedra) 36 (inactive 0)
  [0] <verbose> [FieldCore] Nb. added dofs 567 (number of inactive dofs 0 )
  [0] <verbose> [FieldCore] Build field L
  [0] <verbose> [FieldCore] Nb. of dofs (vertices) 24 (inactive 0)
  [0] <verbose> [FieldCore] Nb. of dofs (edge) 108 (inactive 0)
  [0] <verbose> [FieldCore] Nb. of dofs (triangles) 36 (inactive 0)
  [0] <verbose> [FieldCore] Nb. added dofs 168 (number of inactive dofs 0 )
  [0] <verbose> [FieldCore] Build field S
  [0] <verbose> [FieldCore] Nb. of dofs (vertices) 27 (inactive 0)
  [0] <verbose> [FieldCore] Nb. of dofs (edge) 156 (inactive 0)
  [0] <verbose> [FieldCore] Nb. of dofs (triangles) 90 (inactive 0)
  [0] <verbose> [FieldCore] Nb. added dofs 273 (number of inactive dofs 0 )
  [0] <inform> [FieldCore] Number of dofs 1008
  

• The following lines indicate domain, boundary and skeleton finite element structures have been succesfully constructed.

  [0] <verbose> [FECore] Build Finite Elements dFE
  [0] <inform> [FECore] Finite element dFE added. Nb. of elements added 12
  [0] <verbose> [FECore] Build Finite Elements bFE
  [0] <inform> [FECore] Finite element bFE added. Nb. of elements added 12
  [0] <verbose> [FECore] Build Finite Elements sFE
  [0] <inform> [FECore] Finite element sFE added. Nb. of elements added 30
  [0] <inform> [FECore] Number of adjacencies 1374
  

• Finally, the problem data structure are constructed, in particular knowing finite element DOFs adjacencies and enumeration of DOFs, matrices and vectors can be created.

  [0] <inform> [ProblemsManager] SimpleProblem Nb. local dof 1008 by 1008 nb global dofs 1008 by 1008
  [0] <verbose> [ProblemsManager] SimpleProblem nb. elems 54
  [0] <inform> [ProblemsManager]  FEs ghost dofs on problem SimpleProblem Nb. ghost dof 0 by 0 Nb. local dof 1008 by 1008
  

Allocating finite element instances

Now we push user data operators (UDOs) to pipelines. We have domain pipelines, for evaluating left and the right-hand side, i.e. matrices and vectors. Also, we have pipelines for evaluation elements on boundary and skeleton.

• Domain UDOs instances are created and pushed as follows

      pipeline_mng->getOpDomainRhsPipeline().push_back(new OpRow("U"));
      pipeline_mng->getOpDomainRhsPipeline().push_back(new OpVolume("U"));
      pipeline_mng->getOpDomainLhsPipeline().push_back(
          new OpRowCol("U", "U", true));

• Boundary UDOs instances are created and pushed as follows

      pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpRow("L"));
      pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpFace("L"));
      pipeline_mng->getOpBoundaryLhsPipeline().push_back(
          new OpRowCol("U", "L", false));

• Skeleton UDOs instances are created and pushed as follows

      pipeline_mng->getOpSkeletonRhsPipeline().push_back(new OpRow("S"));
      pipeline_mng->getOpSkeletonRhsPipeline().push_back(
          new OpFaceSide("S", side_fe));

• We explicitly create finite element instance for skeleton, and pushe UDO to it,

      boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> side_fe(
          new VolumeElementForcesAndSourcesCoreOnSide(m_field));
      side_fe->getOpPtrVector().push_back(new OpVolumeSide("U"));

  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 it is realised in a way that the generic element instance is iterated over adjacent entities to faces. We start with allocation memory for "side" finite element as before

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

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 user data operator for each of them.

    CHKERR pipeline_mng->loopFiniteElements();

[0] <inform> **** 0 ****
[0] <inform> **** Operators ****
[0] <inform> Hello Operator OpRow: field name U side 0 type Vertex nb dofs on entity 12
[0] <inform> Hello Operator OpRow: field name U side 0 type Edge nb dofs on entity 9
...
[0] <inform> Hello Operator OpRow: field name U side 0 type Triangle nb dofs on entity 9
[0] <inform> Hello Operator OpRow: field name U side 1 type Triangle nb dofs on entity 9
[0] <inform> Hello Operator OpRow: field name U side 2 type Triangle nb dofs on entity 9
[0] <inform> Hello Operator OpRow: field name U side 3 type Triangle nb dofs on entity 9
[0] <inform> Hello Operator OpRow: field name U side 0 type Tetrahedra nb dofs on entity 3
[0] <inform> 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
[0] <inform> 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 Edge nb dofs on col entity9
...

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 as in MoFEM we always start counting from zero.

Note that 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 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 OpSideFace, 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, user 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 user data operators

Now we focus attention on the implementation of user data operators including OpRow, OpRowCol, OpVolume, etc. which have been used previously. The first operator has the structure

struct OpRow: public MoFEM::ForcesAndSourcesCore::UserDataOperator {
 OpRow(const std::string& field_name):
 ForcesAndSourcesCore::UserDataOperator(field_name,field_name,OPROW) {
 }
 MoFEMErrorCode doWork(int side,EntityType type,EntitiesFieldData::EntData &data);
};

This user 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,EntitiesFieldData::EntData &data);

which as arguments take entity side number (local entity number on the finite element), entity type (e.g. MBVERTEX, MBEDGE, MBTET) and reference to structure MoFEM::EntitiesFieldData::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 user data operator is implemented here

struct OpRowCol: public MoFEM::ForcesAndSourcesCore::UserDataOperator {
 OpRowCol(
   const std::string row_field,
   const std::string col_field,
   const bool symm
 ):
 ForcesAndSourcesCore::UserDataOperator(row_field,col_field,OPROWCOL,symm) {
 }
 MoFEMErrorCode doWork(
   int row_side,int col_side,
   EntityType row_type,EntityType col_type,
   EntitiesFieldData::EntData &row_data,
   EntitiesFieldData::EntData &col_data
 );
};

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 user data operator class can be used, e.g.

struct OpVolume: public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator;
struct OpFace: public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator;
struct OpFaceSide: public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator;

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

struct OpVolumeSide: public MoFEM::VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator;

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

Running the program

In order to run the program, one should first go to the directory where the problem is located, compile the code and then run the executable file. This can be done as follows

cd mofem_install/um/build/tutorials/fun-0
make -j2
./hello_world

Exercises

• Exercise 1: Add OpRow and OpRowCol operator to side_fe
• Exercise 2: Print indices in OpRow operator.
• Exercise 3: Print base functions and derivative of base functions on OpRow.
• Exercise 4: Change space of field "U", "L" and "S", how that changes a number of DOFs on entities?
• Exercise 5: Change order of field "U", "L" and "S", how that changes a number of DOFs on entities?
• Exercise 6: Modify code and calculate a volume of the body. See simple.cpp for example.