v0.14.0
SCL-1: Poisson's equation (homogeneous BC)
Note
Prerequisites of this tutorial include MSH-1: Create a 2D mesh from Gmsh and MSH-2: Create a 3D mesh from Gmsh (for the 3D extension implementation)


Note
Intended learning outcome:
  • general structure of a program developed using MoFEM
  • idea of Simple Interface in MoFEM and how to use it
  • idea of Domain element in MoFEM and how to use it
  • process of implementing User Data Operators (UDOs) to calculate and assemble stiffness matrix and force vector
  • how to push the developed UDOs to the Pipeline
  • a way to handle homogeneous boundary condition in MoFEM
  • utilisation of tools to convert outputs (MOAB) and visualise them (Paraview)


Note
After finishing this tutorial, if you would like to replicate the program and practice yourself in an existing module or in your own module, you may wish to have a look at How to add a new module and program and How to compile a program.

Introduction

The problem

In this tutorial we will solve a simple Poisson's equation with homogeneous (uniform in space) boundary condition using MoFEM. The physical example of this equation is that you have a membrane that is fixed at its boundary and it experiences a uniformly distributed force on its surface, then you are asked to estimate the displacement of the membrane. The strong form of the problem as follows

\[ \begin{align} -\nabla \cdot \nabla u(\mathbf{x}) &= f \quad {\rm in} \quad { \Omega}, \\ { u}(\mathbf{x}) &= 0 \quad {\rm on} \quad \partial { \Omega}, \end{align} \]

where \( { \Omega} \) denotes the domain occupied by the body and \( \partial { \Omega} \) is the boundary of the domain. Additionally, \( f \) is the source function and \( {\bf x} \) is the position in 2D or 3D space of a point in the domain.

Let consdier a Poisson problem on a rectangle domain with homogenous boundary conditions and source function \( f \), the analytical solution can be found. However, as the first practice with MoFEM, we are solving it numerically using a finite element approach. This is done by subdividing the domain into multiple elements, building piece-wise approximations, and solving a discrete problem.

As you may have already learned from the basic finite element method, in order to approximate a field \( u \) of the problem equation. We need to derive the weak form it. The procedure to achieve it is as follows

  • First, multiplying both sides of the equation by a test function \( \delta u \) and then integrate over the domain \( \Omega \)

    \[ \begin{equation} -\int_\Omega \delta u(\nabla \cdot \nabla u ) ~d\Omega= \int_\Omega \delta u f ~d\Omega. \end{equation} \]

  • Second, apply the integration by parts on the left-hand side of the equation, we have

    \[ \begin{equation} \int_\Omega \nabla \delta u \cdot \nabla u ~d\Omega - \int_{\partial \Omega} \delta u {\bf n} \cdot \nabla u ~d\partial \Omega = \int_\Omega \delta u f ~d\Omega. \end{equation} \]

    It is noted that the test function \( \delta u \) has to satisfy the homogeneous boundary condition, i.e. \( \delta u=0 \) on \( \partial \Omega \). Substituting it to the above equation, we finally obtain the weak form of the problem

    \[ \begin{equation} \int_\Omega \nabla \delta u \cdot \nabla u ~d\Omega = \int_\Omega \delta u f ~d\Omega, \quad \forall \delta u \in H_{0}^{1}(\Omega), \end{equation} \]

    where the space for test function \( \delta u \) is \( H_{0}^{1}(\Omega):=\left\{v \in H^{1}(\Omega) \mid v=0 \text { on } \partial \Omega\right\} \).

We are now solving this equation of the weak form instead of the original equation. This equation is asking for the solution of \( u \) that is true for all test function \( \delta u \). In order to achieve it, we will approximate \( u \) following a process in finite element called discretisation which will be presented in the next part.

Discretisation

As mentioned above, instead of trying to solve the problem analytically, we will approximate \( u \) assuming its solution has the form

\[ \begin{equation} u \approx u^h = \sum_{j=0}^{N-1} \phi_j \bar{u}_j. \end{equation} \]

This expression can be interpreted as follows. Find the approximate solution \( u^h \) of \( u \) where \( u^h \) is calculated by summing the contribution of base function \( \phi_j \) with the coefficient of the contribution is \( \bar{u}_j \).

Sometimes, \( \phi_j \) is also called shape functions and \( \bar{u}_j \) called degrees of freedom (DOFs) of the problem. In the implementation process, which will be discussed later in this tutorial, MoFEM gives you values of \( \phi_j \) by default (provided that you give it some hints) and you will find the solution of \( u^h \), of course, with the help of MoFEM.

Keep in mind that, in the weak form, we have another term that also needs to be taken care of, that is the test function \( \delta u \). As we were saying, the weak form has to be satisfied with all test function \( \delta u \) meaning it has to be true with the arbitrary choice of \( \delta u_i=\phi_i \).

Substituting \( u \) and \( \delta u \) into the weak form, we have the discrete form of the problem given by

\[ \begin{equation} \int_{\Omega^e} \nabla \phi_i \cdot \nabla \left( \sum_j \bar{u}_j \phi_j \right) ~d{\Omega^e}= \int_{\Omega^e} \phi_i f ~d{\Omega^e}. \end{equation} \]

Moving \( \bar{u}_j \) outside of the parentheses and rearranging the equation, we have

\[ \begin{equation} \sum_j \left( \int_{\Omega^e} \nabla \phi_i \cdot \nabla \phi_j ~d{\Omega^e} \right) \bar{u}_j = \int_{\Omega^e} \phi_i f ~d{\Omega^e}. \end{equation} \]

Now, the problem becomes: finding the vector of coefficients (or DOFs) \( {\bf U} \) such that

\[ {\bf K u} = {\bf F}, \]

where \( {\bf K} \) and \( {\bf F} \) are the global stiffness matrix (left hand side) and global force vector (right hand side) calculated over the whole domain, respectively. \( {\bf K} \) and \( {\bf F} \) are obtained by assembling all elements (entity) in the domain, and the components of element stiffness matrix and element force vector are calculated as

\[ \begin{align} K_{ij}^e &= \int_{\Omega^e} \nabla \phi_i \cdot \nabla \phi_j ~d{\Omega^e}, \\ F_i^e &= \int_{\Omega^e} \phi_i f ~d{\Omega^e}. \end{align} \]

One thing we can notice here is that the matrix \( {\bf K} \) is symmetric which means there will be opportunity to save time later in the implementation.

We are almost at the place where we can start the implementation, but still, you may ask how computers can handle the integral terms. Of course, we will not ask computers to calculate the integrals directly in an infinite sense, instead it is done using quadrature which is in a finite sense and commonly used in finite element method. In other word, the integrals are approximated by the sum of a set of points on each element along with their weights as follows

\[ \begin{align} K_{ij}^e &= \int_{\Omega^e} \nabla \phi_i \cdot \nabla \phi_j \approx \sum_q \nabla \phi_i \left( {\bf x}_q \right) \cdot \nabla \phi_j \left( {\bf x}_q \right) W_q \left\|\mathbf{J}_{g}^{e}\right\| \label{eqn_stiffness}, \\ F_i^e &= \int_{\Omega^e} \phi_i f \approx \sum_q \phi_i \left( {\bf x}_q \right) f \left( {\bf x}_q \right) W_q \left\|\mathbf{J}_{g}^{e}\right\|, \end{align} \]

where \( {\bf x}_q \) and \( W_q \) are the location and weight of the quadrature point, respectively. Meanwhile, \( \left\|\mathbf{J}_{g}^{e}\right\| \) is the determinant of the Jacobian of the transformation of the element \( e \) from physical coordinates to parent coordinates. For triangular elements, this determinant equals to the area of the element. These equations can be interpreted as follows

The approximated value of the component at row \( i \) and column \( j \) of the element stiffness matrix \( K \) of element \( e \) is calculated by taking the summation over all quadrature points of the multiplication of

  • the first derivative of the \( i^{\rm th} \) approximation function \( \phi \) evaluated at quadrature point \( q \),
  • the first derivative of the \( j^{\rm th} \) approximation function \( \phi \) evaluated at the same quadrature point \( q \), and
  • the multiplication of the weight of that quadrature point \( q \) and the determinant of the Jacobian.

Likewise, the approximated value of the component at row \( i \) of the element force vector \( F \) of element \( e \) is calculated by taking the summation over all quadrature points of the multiplication of

  • the \( i^{\rm th} \) approximation function \( \phi \) evaluated at quadrature point \( q \),
  • the body force evaluated at the same quadrature point \( q \), and
  • the multiplication of the weight of that quadrature point \( q \) and the determinant of the Jacobian.

More details on the numerical integration can be found at the tutorial FUN-1: Integration on finite element mesh.

To this end, we have established a linear system of the primary variable \( {\bf U} \) through \( {\bf K} \) and \( {\bf F} \). In MoFEM, we will use an iterative solver (KSP) to solve for \( {\bf U} \) and apply steps to postprocess the solution and visualise it.

Implementation

An immediate question you may have regarding the implementation is how to implement the matrix \( {\bf K} \) and vector \( {\bf F} \). This is done through the implementation of the so-called User Data Operators. UDOs are the essential part of MoFEM and they are present in all finite element problems implemented in MoFEM. UDO is normally called Operator (for short) by MoFEM developers. Once UDOs are implemented, they will be pushed to the working pipeline.

Here we have introduced two important concepts in MoFEM

  • UDOs are responsible for calculation of certain things, e.g. stiffness matrix, force vector, inverse of Jacobian, field values, field gradients, etc., and
  • Pipeline is where you push UDOs to in a certain order, e.g. UDOs for calculating field values and field gradients (nonlinear problem) are pushed before UDOs calculating stiffness matrix and force vector. After being pushed to the Pipeline, those UDOs will be executed sequentially.

As a common practice, typically all the implementation of UDOs for a specific problem is put in a *.hpp file. This file will be included in the main *.cpp file which contains all the necessary classes and functions. Detailed explanation of the implementation of UDOs, how you pushed UDOs to Pipeline as well as development of classes/functions are presented below

User Data Operators

As described previously, solving the Poisson problem with homogeneous would require the computation and assembling of the matrix \( {\bf K} \) and vector \( {\bf F} \). These essential processes of computation and assembling will be handled by two different UDOs which separately deal with the matrix and vector. They are:

  1. Poisson2DHomogeneousOperators::OpDomainLhsMatrixK is responsible for calculating and assembly of the left-hand-side matrix \( {\bf K} \)
  2. Poisson2DHomogeneousOperators::OpDomainRhsVectorF is responsible to calculation and assembly of the right-hand-side vector \( {\bf F} \)

Before going into details of the implementation of the two UDOs in poisson_2d_homogeneous.hpp, let's have a look at the first few lines of code file that includes some libraries for finite element implementation. The code is based on C++ libraries with advanced aliases, and declaration/initialisation. Apart from that, the dimension(2D/3D) of the problem can easily be assihned from the EXECUTABLE_DIMENSION which will discuss later.

// Define name if it has not been defined yet
#ifndef __POISSON2DHOMOGENEOUS_HPP__
#define __POISSON2DHOMOGENEOUS_HPP__
// Include standard library and Header file for basic finite elements
// implementation
template <int DIM>
using ElementsAndOps = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>;
FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
// Namespace that contains necessary UDOs, will be included in the main program
// Declare FTensor index for a problem
// For simplicity, body source term f has been consider to be constant throughout the domain
const double body_source = 5.;
//.. ..Implementation of the UDOs below this point.. ..
}

Next, Users can also learn a general structure of any UDOs implemented in MoFEM to many extend. So, we will look in details how the following two main User Defined Operator(UDO)s of the problem can be implemented.

  1. OpDomainLhsMatrixK and
  2. OpDomainRhsVectorF

OpDomainLhsMatrixK

This UDO is responsible for the calculation and assembly of the left-hand-side matrix \( {\bf K} \)

First, let's look at the structure of this UDO to calculate matrix \( {\bf K} \).

struct OpDomainLhsMatrixK : public AssemblyDomainEleOp {
public:
OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)
: AssemblyDomainEleOp(row_field_name, col_field_name,
DomainEleOp::OPROWCOL) {}
//.. ..Implementation of iNtegrate() below this point .. ..
}
};

We have the class OpDomainLhsMatrixK that is inherited from AssemblyDomainEleOp which is the alias of the base chain of specialized template FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase. This operator also accelerate the process of assembling local matrix to global matrix. Here, the forms integrators using the DomainEleOp type is a templated parameter to assembly using the PETSc library. Here, the newly created class has two public objects OpDomainLhsMatrixK() and iNtegrate which can be accessed from outside of the class. This follows the concept of data encapsulation to hide values and objects inside the class as much as possible. We will discuss a little more details of those four objects of the class in the followings

  • OpDomainLhsMatrixK()
    OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)
    : AssemblyDomainEleOp(row_field_name, col_field_name,
    DomainEleOp::OPROWCOL) {}
    This public member function has two input arguments which are the row_field_name and col_field_name carrying the names of the field for row and column, respectively. In general, these two field names can be different.But in this particular case of Poisson's problem with homogeneous boundary condition, we have considered the same field named as "U" for both row and column.In this way, OpFaceEle::OPROWCOL because we assemble stiffness matrix that has data on both row and column. In this case, our matrix \( K \) is symmetric but if we want we could code sYmm = true inside the closing bracket{} to assemble and store data of half of the matrix. However, here we are considering the full matrix.
  • iNtegrate
    // Implementation of iNtegrate() below this point
    }
    This public member function is the core part of the UDO implementation. It decides how to calculate the local matrix components and how to assemble them to the global matrix. All mathematical derivations that we did in the previous section of this tutorial will be implemented in this function. More discussion on other essential function of it is stating below:
  • locMat
    This private member object of type MatrixDouble is used to store the results of the calculation of components of element stiffness matrix. This object is made private because only member of the class know how to use it, members from other classes have no access to this private object avoiding unpredictable consequences, i.e. errors.

Now we take a closer look on the details of the implementation of the member function iNtegrate(). Full code for this function as follows

const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
this->locMat.resize(nb_row_dofs, nb_col_dofs, false);
this->locMat.clear();
// get element area
const double area = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on row
auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
// get derivatives of base functions on column
auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
for (int cc = 0; cc != nb_col_dofs; cc++) {
this->locMat(rr, cc) += t_row_diff_base(i) * t_col_diff_base(i) * a;
// move to the derivatives of the next base functions on column
++t_col_diff_base;
}
// move to the derivatives of the next base functions on row
++t_row_diff_base;
}
// move to the weight of the next integration point
++t_w;
}
}
};

First, let's discuss the structure of this function

For specific tasks, the data structure is adapted based on the problem's dimension. Such as solving a 2D or 3D Poisson problem, the code can utilize different entities types e.g. edges, faces, volumes. For example, the entities or elements are provided accordingly using PipelineManager::ElementsAndOpsByDim<SPACE_DIM> which is aliased as ElementsAndOps.The DomainEle::UserDataOperator is used to assemble using the forms integrator AssemblyDomainEleOp. This can be seen in the code, particularly in the struct OpDomainLhsMatrixK. The function's parameters of iNtegrate are row_data and col_data are referencing to objects of type EntitiesFieldData::EntData. These parameters make ready the data related to degrees of freedom, basis function values, and field values for the rows and columns of the operation.

In MoFEM, functions are written with check and error handling using keywords like MoFEMErrorCode, MoFEMFunctionBegin, and MoFEMFunctionReturn(0). These indicate the execution of functions and the success or failure of the operation. These arguments are used to handle errors like try() or catch() in basic c++ systematically.

Now, moving into the main implementation of iNtegrate() function, we will see information about the number of DOFs on row and column are extracted from the database

const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();

Once the program is confident that the data for row and column are all valid, it starts to initialise the local stiffness matrix whose components will be calculated and assembled to the global matrix.

locMat.resize(nb_row_dofs, nb_col_dofs, false);
locMat.clear();

Then you will get the area for 2D or volume for 3D which is needed later when you integrate the function to calculate the stiffness matrix.

// get element area
const double area = getMeasure();

It is worth noting that getMeasure() is a generic function and the value you get depends on which entity type you are working with. For example, in this UDO, you are working with face entities (triangles), getMeasure() gives you face area. If you are dealing with edge (for boundary element) or volume entities (for volume domain), you would get edge length or element volume, respectively.

Next, as required for the calculation of stiffness matrix in Eq. (3), we need the number of integration points ( \( q \)) and their weights ( \( W_q \)) which can be done as follows

// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();

Here getFTensor0IntegrationWeight() is a function of MoFEM::ForcesAndSourcesCore::UserDataOperator and it returns a zeroth-order FTensor object (Tensor template library) that stores value of integration weight. We prefer to use FTensor anywhere we can because it is compact and highly efficient. In MoFEM implementation, FTensor object is normally named with the prefix t_.

Once the number of Gauss points and their weights are determined, the remaining component to calculate stiffness matrix would include the gradients of the basis function evaluated at the integration (Gauss) point, \( \nabla \phi_i, \nabla \phi_j \) ,for row and column, respectively, and the loop over all the Gauss points. They are done in the following part of the code

// get derivatives of base functions on row
auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
// get derivatives of base functions on column
auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
for (int cc = 0; cc != nb_col_dofs; cc++) {
locMat(rr, cc) += t_row_diff_base(i) * t_col_diff_base(i) * a;
// move to the derivatives of the next base functions on column
++t_col_diff_base;
}
// move to the derivatives of the next base functions on row
++t_row_diff_base;
}
// move to the weight of the next integration point
++t_w;
}

Here you can see inside the loop over the integration points, we have two other nested loops that loop over the row ( \(i\)) and column ( \(j\)) DOFs of the matrix which is ultimately calculated similarly to the way they are presented in Eq. (11)

locMat(rr, cc) += t_row_diff_base(i) * t_col_diff_base(i) * a;

where a is the intermediate variable of type double and calculated earlier as

const double a = t_w * area;

where area can be considered as the determinant of the Jacobian of the transformation from the physical (global) coordinate system to the reference (local) coordinate system. See more about the integration and Jacobian at FUN-1: Integration on finite element mesh.

It should be noticed that at the end of each loop over the row and column DOFs, t_row_diff_base and t_col_diff_base which have the values of the gradients of the basis function are moved/shifted to the next chunk of the memory which stores values of another set of gradients of basis function associated with DOFs. The same technique is applied to move t_w to the weight of the next Gauss point in the memory.

Note: You might have also recognised that in the traditional finite element implementation, you will probably have, at the same place of the code, the structure of the main nested loop like this

  • Loop over all elements
    • Loop over all integration points
      • Calculation of the local stiffness matrix components and fill them to the global matrix

In MoFEM's finite element implementation, there are some differences worth noting. While nested loops are still used, the implementation of the User Data Operator (UDO) involves a single loop over integration points. Within this loop, local stiffness matrix components are computed and then assembled into the global matrix. The loop over elements referred to as entities in the context of Hierarchical basis functions Details of how this loop is triggered is presented in solveSystem().

That concludes the description of the code segment responsible for calculating the components of the element stiffness matrix.

Figure 1: Symmetric matrix assembling (for a general 3D case).

OpDomainRhsVectorF

In this part, we will talk about the implementation of the second UDO Poisson2DHomogeneousOperators::OpDomainRhsVectorF which is responsible for the calculation and assemble of the right-hand-side force vector \( {\bf F} \)

Similar to what implemented for the \( {\bf K} \) matrix we discussed above, the implementation of the UDO for the \( {\bf F} \) vector has very similar structure with slightly different functions

struct OpDomainRhsVectorF : public AssemblyDomainEleOp {
public:
OpDomainRhsVectorF(std::string field_name)
// .. ..Implementation of iNtegrate() below this point.. ..
}
};

We have the public objects of OpDomainRhsVectorF() and iNtegrate() and private object of locF. Here the name of the main function is changed to OpDomainRhsVectorF to reflect what it does. Additionally, as we are calculating and assembling vector instead of matrix, we need only one input of field_name and we tell MoFEM that we are doing the operation of row for the vector by specifying DomainEleOp::OPROW. Similarly for function iNtegrate(), it requires all three groups of the input but only one for each group instead of two as we had for the stiffness matrix.

The implementation of this essential function of iNtegrate() for force vector calculation and assembling also has very similar structure to its matrix UDO counterpart. The code for this function is as follows

const int nb_dofs = data.getIndices().size();
this->locF.resize(nb_dofs, false);
this->locF.clear();
// get element area
const double area = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get base function
auto t_base = data.getFTensor0N();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_dofs; rr++) {
this->locF[rr] += t_base * body_source * a;
// move to the next base function
++t_base;
}
// move to the weight of the next integration point
++t_w;
}
}
};

A part from the input arguments of this function is slightly different from the UDO for stiffness matrix. In this UDO for force vector, we use the base function value to calculate the components of local vector reflecting what is presented in Eq. (4) and only need to loop over row DOFs as follows

// get base function
auto t_base = data.getFTensor0N();

and later using the constant body source (force) that is predefined at the beginning of the *.hpp to calculate the local vector

for (int rr = 0; rr != nb_dofs; rr++) {
this->locF[rr] += t_base * body_source * a;
// move to the next base function
++t_base;
}

That is all for the implementation of UDOs that are responsible for the calculation and assembling of the LHS matrix \( {\bf K} \) and RHS vector \( {\bf F} \). Having the essential components implemented, you now may ask how those UDOs are called in the main program. This is done through the process called push operator to the Pipeline and you will find out in more details later in this tutorial, assembleSystem().

Next, we will look in detail how the main program, including class and functions, is implemented.

The Poisson2DHomogeneous class

This main class Poisson2DHomogeneous contains functions and each of which is responsible for a certain task of a finite element program. Intially, the field name has named as "U" string to recognsie it well. Next the EXECUTABLE_DIMENSION is the considered for the dimension of the mesh. To generalize, the value is assined to the Space dimension of problem.

constexpr auto field_name = "U";

The public part of the class includes a constructor and function runProgram() that calls other functions to perform finite element analysis.

public:
// Declaration of the main function to run analysis
MoFEMErrorCode runProgram();

Then there are private functions doing certain tasks and can be recognised by their names

private:
// Declaration of other main functions called in runProgram()
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode setIntegrationRules();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();

It is followed by the declaration of member variables that will be used in one or some of the member functions declared above

// MoFEM interfaces
Simple *simpleInterface;
// approximation order
int oRder;
Note
For writhing program in MoFEM, we follow some coding practices for code style and naming convention described at Coding practice

Poisson2DHomogeneous()

Now is the constructor

The constructor specifies

  • mField(m_field): The MoFEM instance that is the backbone of the program

readMesh()

Now the first function that actually does some finite element task is the function responsible for reading input mesh

Apart from the codes for error handling, the main three lines of code is a standard way to read an input mesh using Simple interface. Interface in MoFEM is a set of rules through which program developers use to setup a problem. Simple interface provides simplest but also less flexible way to setup a problem. More on the interfaces can be found at MoFEM interfaces. Implementation of more advanced interfaces will be presented in later tutorials.

setupProblem()

Next is the function that is responsible for setting up the finite element problem

This function has

  • simpleInterface->addDomainField: It is important to understand that we need to add the domain field through the simpleInterface only as the field at the boundary is zero (homogeneous) for this particular example hence no need to add the boundary field. The addition of the domain field requires the followings
    • Field name ("U"),
    • Approximation space (H1 - scalar space) - function space can also be H(curl), H(div), L2 depending on physical properties of the field you are approximating, see more at FieldSpace
    • Approximation bases (AINSWORTH_BERNSTEIN_BEZIER_BASE) - base can also be AINSWORTH_LEGENDRE_BASE, AINSWORTH_LOBATTO_BASE, DEMKOWICZ_JACOBI_BASE , see more at FieldApproximationBase, and
    • Number of DOFs per shape function (1) as the current problem is a scalar-field problem. This will be 3 for vector-field problem.
  • simpleInterface->setFieldOrder: Set the polynomial order of the approximation of the field
  • simpleInterface->setUp(): Finally, do the setup of the problem.

boundaryCondition()

This function deals the boundary condition and its implementation for the current problem is as follows

auto bc_mng = mField.getInterface<BcManager>();
// Remove BCs from blockset name "BOUNDARY_CONDITION" or SETU, note that you
// can use regular expression to put list of blocksets;
CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
simpleInterface->getProblemName(), "(BOUNDARY_CONDITION|SETU)",
std::string(field_name), true);
}

To ensure the code reads the boundary conditions, it's essential to generate the mesh that includes the required attributes or values. Specifically, blocks in the mesh must be labelled with the name "BOUNDARY_CONDITION". In other cases where multiple boundary conditions are involved, you should assign distinct names using numbers, such as "BOUNDARY_CONDITION_1", "BOUNDARY_CONDITION_2", "BOUNDARY_CONDITION_3", and so on. This naming scheme helps in distinguishing and applying different boundary conditions to their respective blocks within the mesh.

Though the application procedure is pretty straightforward but the underlying process of applying these boundary conditions is not immediately obvious. Users can jump to the next section without understanding the underlying process of it.

If we consider \( {\bf u} \) as the solution vector such that \( {u} = [u_1, u_2, u_3, u_4, \ldots, u_n] \), we will essentially eliminate the DoFs associated with entities that are part of the "BOUNDARY_CONDITION" blocksets. In this context, let's say the \( \bf{u_0} \) represents the vector after removing the DOF from BC. For example, \( {u_0} = [u_1, u_2, u_3,.. .. u_\text{n-m}, 0 ,0\ldots, 0] \). However, we need to reintroduce these removed DoFs later on to make the solution admissible for the boundary values as well. To accomplish this, we will incorporate an additional vector \( {\bf u^e} \) with the current \( {\bf u_0} \), as in (14).

\[ \begin{equation} \bf {Ku}= f \end{equation} \]

\[ \begin{equation} \bf {K(u_0 + u_e)}= f \end{equation} \]

Here, \( u_e \) is the vector where we will impose those removed degrees of freedom in it. For example,

\[ \begin{equation} {\mathbf { u}_e} = \left[ \begin{array}{cccc} {0} & {0} & {0} & \ldots & {u_\text{n-m+1} } & {u_\text{n} } & \end{array} \right] \label{eq:u_e} \end{equation} \]

Here, \( u_\text{n-m+1} \) and \(u_\text{n} \) denotes DOFs the of BOUNDARY_CONDITION Blockset.

This approach will allow us to streamline the procedure and ensure that the removed degrees of freedom are properly reintegrated into the calculations. However, to simplify this cumulative process, we will address rest of the part during the assembling phase in the upcoming sub-section.

assembleSystem()

This part is about the function that is responsible for the assembling of the system of equations. In between, we will also go through the code for setting the DOFs to the BC.

At the beginning of the section assembleSystem(), recall the two important concepts in MoFEM, namely UDO and Pipeline. While how the implementation of UDOs has been shown in User Data Operators, here you will see how the implemented UDOs are pushed to the Pipeline.

Apart from pushing UDOs to the main program, assembleSystem() ()function is mainly responsible for setting operators implementing to the pipelines for KSP solver (iterative linear solver provided by PETSc). The full source code for the function is as follows:

auto pipeline_mng = mField.getInterface<PipelineManager>();
{ // Push operators to the Pipeline that is responsible for calculating LHS
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainLhsMatrixK(field_name, field_name));
}
{ // Push operators to the Pipeline that is responsible for calculating RHS
auto set_values_to_bc_dofs = [&](auto &fe) {
auto get_bc_hook = [&]() {
EssentialPreProc<TemperatureCubitBcData> hook(mField, fe, {});
return hook;
};
fe->preProcessHook = get_bc_hook();
};
// you can skip that if boundary condition is prescribing zero
auto calculate_residual_from_set_values_on_bc = [&](auto &pipeline) {
using DomainEle =
using OpInternal = FormsIntegrators<DomainEleOp>::Assembly<
auto grad_u_vals_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>(field_name,
grad_u_vals_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInternal(field_name, grad_u_vals_ptr,
[](double, double, double) constexpr { return -1; }));
};
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
pipeline_mng->getOpDomainRhsPipeline(), {H1});
set_values_to_bc_dofs(pipeline_mng->getDomainRhsFE());
calculate_residual_from_set_values_on_bc(
pipeline_mng->getOpDomainRhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpDomainRhsVectorF(field_name));
}
}
  • First, the pipeline manager (pipeline_mng) that manages two main pipelines which are:
    • getOpDomainLhsPipeline: responsible for calculations of the left hand side matrix
    • getOpDomainRhsPipeline: responsible for calculations of the right hand side vector
      auto pipeline_mng = mField.getInterface<PipelineManager>();
  • Second, pushing the operators to the Pipeline that is responsible for the calculation of the LHS matrix. This should be done in two steps. Initially, the operators to calculate the inverse of the Jacobian of the mapping from reference (parent) space to the physical space is pushed first and then the implemented UDO for stiffness matrix OpDomainLhsMatrixK is pushed to the LhsPipeline.
    { // Push operators to the Pipeline that is responsible for calculating LHS
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
    pipeline_mng->getOpDomainLhsPipeline(), {H1});
    pipeline_mng->getOpDomainLhsPipeline().push_back(
    new OpDomainLhsMatrixK(field_name, field_name));
    }
    It is to be noted that the procedure to calculate the inverse of Jacobian with the presence of the gradient of approximation function is needed but it is done automatically and user does not have to manually add codes to push any operators calculating the inverse of Jacobian.
  • Third, previously in in boundaryCondition(), we modified \( \bf{u} \) and we were waiting to impose the removed DoFs to get the admissible solution. Here, this operation is done by the lamda function 'set_values_to_bc_dofs'.

We enforce the following code snippet to set dofs. First we get the reference \( \bf{u}\) for BC. Then using 'get_bc_hook' as preprocessing of object 'fe' we set the boundary condition values to \( \bf{u_e} \).

auto set_values_to_bc_dofs = [&](auto &fe) {
auto get_bc_hook = [&]() {
EssentialPreProc<TemperatureCubitBcData> hook(mField, fe, {});
return hook;
};
fe->preProcessHook = get_bc_hook();
};

Note:For the current homogenous problem, if we consider zero boundary condition, that endsup with all the values of \( \bf {u_e} \) as zero. So then the remainings from the equation (15) is \( \bf {K u_0} = \bf {f}\). However, in case of non-zero BCs we have to consider \(\bf {K u_e}\) in the next part.

So the rest of the code should be organised as follows:

// you can skip that if boundary condition is prescribing zero
auto calculate_residual_from_set_values_on_bc = [&](auto &pipeline) {
using DomainEle =
using OpInternal = FormsIntegrators<DomainEleOp>::Assembly<
auto grad_u_vals_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>(field_name,
grad_u_vals_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInternal(field_name, grad_u_vals_ptr,
[](double, double, double) constexpr { return -1; }));
};
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
pipeline_mng->getOpDomainRhsPipeline(), {H1});
set_values_to_bc_dofs(pipeline_mng->getDomainRhsFE());
calculate_residual_from_set_values_on_bc(
pipeline_mng->getOpDomainRhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpDomainRhsVectorF(field_name));
}
}

Here, the lamda function calculate_residual_from_set_values_on_bc is to get the \(\bf {K{u_e}}\). Inside of it, the OpInternal integrator serves a dual purpose: it undertakes system of equations and assembles the matrices, effectively performing a linear multiplication. Upon pushing OpInternal to the pipeline manager, -1 is passed as it is alongside \( \bf K {u_e} \) as in equation (17).

\[ \begin{equation} \bf {K u_0} = f - \bf {K u_e} \end{equation} \]

We do not need derivative of the approximation function to calculate force vector in UDO OpDomainRhsVectorF. Consequently, the operators for recalculating inverse of Jacobian is not needed as well as it is integrated with MoFEM.

setIntegrationRules()

This function is responsible for setting the Gauss integration rules and it looks like this

auto rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
auto rule_rhs = [](int, int, int p) -> int { return p; };
auto pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
}

Here \( p \) is the polynomial order of the approximation function. For the LHS matrix, we are calculating the integral of function \(\nabla \phi_i \cdot \nabla \phi_j\) so the polynomial order of the integral is \( p - 1 \) plus \( p - 1 \) resulting \( 2 * (p - 1)\) as we see in the implementation of the integration rules for the LHS. Similarly, as we calculate the integral of \( \phi_i\) for the RHS, we use \( p \) as the integration rule of the RHS. Having the integration rules, MoFEM will automatically determine the number of integration (Gauss) points that need to be used for each entity (element). Here you can see MoFEM allows you to choose different integration rules for different operators.

solveSystem()

Having the computation of LHS and RHS is defined in the previous function. We now can actually solve the system of equations using iterative KSP solver from PETSc.

The codes for the function that is responsible for solving the systems of equations look like this

auto pipeline_mng = mField.getInterface<PipelineManager>();
auto ksp_solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(ksp_solver);
CHKERR KSPSetUp(ksp_solver);
// Create RHS and solution vectors
auto dm = simpleInterface->getDM();
auto F = createDMVector(dm);
auto D = vectorDuplicate(F);
// Solve the system
CHKERR KSPSolve(ksp_solver, F, D);
// Scatter result data on the mesh
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}

It starts first with getting the pipeline manager

auto pipeline_mng = mField.getInterface<PipelineManager>();

Then create the KSP solver using wrapped function in MoFEM (pipeline_mng->createKSP()) and setup the solver using PETSc functions

auto ksp_solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(ksp_solver);
CHKERR KSPSetUp(ksp_solver);

Next is getting the Discrete Manager (dm) which is a common object allowing things implemented in MoFEM talk to things implemented in PETSc before initilising the RHS and solution vectors

// Create RHS and solution vectors
auto dm = simpleInterface->getDM();
auto F = createDMVector(dm);

At this particular point, Discrete Manager allows the two pipelines (responsible for LHS and RHS) implemented in MoFEM to be used as the input for KSP solver implemented in PETSc. From that, the solution vector \( {\bf U} \) of the system of equations \({\bf KU=F} \) will be obtained when the KSP solver is triggered by this

// Solve the system
CHKERR KSPSolve(ksp_solver, F, D);

It is important to note that all the implementation of the UDOs presented in the previous sections was just the definition. All the calculations (all the loops to calculate matrix and vector entries implemented in UDOs) are only triggered when the line of code above calling KSPSolve() function is executed.

Lastly, the results are scattered through DM and ready to be fed to the output mesh in the next step

// Scatter result data on the mesh
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);

outputResults()

This function is solely responsible for the postprocessing of the results writing the calculated field values to the output mesh.

auto pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
post_proc_fe->getOpPtrVector(), {H1});
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues(field_name, u_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>(field_name, grad_u_ptr));
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{{"U", u_ptr}},
OpPPMap::DataMapMat{{"GRAD_U", grad_u_ptr}},
)
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile("out_result.h5m");
}

runProgram()

Finally, having all the necessary tasks implemented in the corresponding functions, we can now put them together in the last function of the main class. This function is responsible for the calling sequence which is similar to most of other finite element programs

The main() function

This main() function does not do much job apart from creating the top-level class and call the function to trigger the analysis

int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Error handling
try {
// Register MoFEM discrete manager in PETSc
DMType dm_name = "DMMOFEM";
// Create MOAB instance
moab::Core mb_instance; // mesh database
moab::Interface &moab = mb_instance; // mesh database interface
// Create MoFEM instance
MoFEM::Core core(moab); // finite element database
MoFEM::Interface &m_field = core; // finite element interface
// Run the main analysis
Poisson2DHomogeneous poisson_problem(m_field);
CHKERR poisson_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}

So you can see, at the beginning, it creates the Discrete Managers that enable information flows between MoFEM (finite element implementation), MOAB (element topology management), and PETSc (algebraic solvers). After that, it creates variable poisson_problem of type/class Poisson2DHomogeneous which is previously defined and then run the analysis by triggering the public function runProgram().

// Run the main analysis
Poisson2DHomogeneous poisson_problem(m_field);
CHKERR poisson_problem.runProgram();

Results

Run the program

In order to run the program that we have been discussing in this tutorial, you will do the following steps

  • First, go to the directory where the binary file poisson_2d_homogeneous is located. Depending on how you install MoFEM shown in this page Installation with Spack (Scripts), going to the directory would be something similar to this
    • For user version installation
      cd mofem_install/um_view/tutorials/scl-1/
    • For developer version installation
      cd mofem_install/mofem-cephas/mofem/users_modules/um-build-RelWithDebInfo-abcd1234/tutorials/scl-1
  • Second, check the parameters in the param_file.petsc. These are PETSc parameters and you should only use parameters that are needed for a particular solver, in this case KSP solver. Only the following parameters should be uncommented
    ## Linear solver
    -ksp_type fgmres
    -pc_type lu
    -pc_factor_mat_solver_type mumps
    -ksp_monitor
  • Third, in the terminal, run commands to partition the input mesh and start the analysis
    mofem_part -my_file square.h5m -output_file square_2parts.h5m -my_nparts 2 -dim 2 -adj_dim 1
    mpirun -np 2 ./poisson_2d_homogeneous -file_name square_2parts.h5m -order 4
    where the mesh of a square plate is partitioned into two parts and then the program is run using two processors (the same number of partitions) with fourth order polynomial of approximation.

Output

Once the analaysis is complete, you see all output messages printed to the terminal

  • Version of MoFEM used to run the analysis
    MoFEM version 0.11.0 (MOAB 5.2.1 Petsc Release Version 3.11.4, Sep, 28, 2019 )
    git commit id 0ac8895b15c4ad41ea5e7077c4d011f7efb50f13
  • Meshset and entity blocks defined in the input mesh
    read cubit meshset 12682136550675316738 type BLOCKSET UNKNOWNNAME msId 1 name square
    read cubit meshset 12682136550675316739 type BLOCKSET UNKNOWNNAME msId 2 name boundary
    read cubit meshset 12682136550675316740 type BLOCKSET UNKNOWNNAME msId 3 name surface
    read cubit meshset 12682136550675316741 type BLOCKSET MAT_ELASTICSET msId 100 name MAT_ELASTIC
    read cubit meshset 12682136550675316742 type BLOCKSET UNKNOWNNAME msId 110 name BOUNDARY_CONDITION
  • Domain field name, approximation space and bases, as well as rank (1 for scalar problem) as implemented in setupProblem()
    Add field U field_id 1 space H1 approximation base AINSWORTH_BERNSTEIN_BEZIER_BASE rank 1 meshset 12682136550675316745
    Add finite element dFE
  • General information about the problem including removing DOFs associated with the boundary
    Add problem SimpleProblem
    Number of dofs 1687
    Number of dofs 1659
    Finite element dFE added. Nb. of elements added 205
    Finite element dFE added. Nb. of elements added 201
    Number of adjacencies 1435
    SimpleProblem Nb. local dof 1687 by 1687 nb global dofs 3289 by 3289
    SimpleProblem Nb. local dof 1602 by 1602 nb global dofs 3289 by 3289
    FEs ghost dofs on problem SimpleProblem Nb. ghost dof 0 by 0 Nb. local dof 1687 by 1687
    FEs ghost dofs on problem SimpleProblem Nb. ghost dof 57 by 57 Nb. local dof 1602 by 1602
    removed ents on rank 0 from problem SimpleProblem dofs [ 1650 / 1650 (before 3289 / 3289) local, 0 / 0 (before 0 / 0) ghost, 3209 / 3209 (before 1687 / 1687) global]
    removed ents on rank 1 from problem SimpleProblem dofs [ 1559 / 1559 (before 3289 / 3289) local, 55 / 55 (before 57 / 57) ghost, 3209 / 3209 (before 1602 / 1602) global]
  • Convergence of the KSP iterative solver from PETSC. As shown, it is converged is one step which is expected for this simple linear problem.
    0 KSP Residual norm 1.057541430518e-01
    1 KSP Residual norm 2.720205258196e-15

Then, you also see in the directory where you run the analysis, it now has the newly created output file, namely out_result.h5m. The output can be visualised in a visualisation software. If you would like to open the output in the free software of Paraview, you would need to convert the input file to *.vtk format by running the following command line in your terminal

mbconvert out_result.h5m out_result.vtk

Then open it in Paraview and use the filter WarpByScalar, you will be able to see the deformation as below

Figure 2: Poisson homogeneous visualisation.

Discussion

As mentioned at the beginning of this tutorial, this Poisson equation with homogeneous boundary condition helps to predict the deformation of a membrane that is fixed at the boundary bearing an uniformly distributed force on its surface. In this case, the force \( f=5.0\) is hardcoded in the code.

You can test yourself how the increase/decrease in approximation order affects the number of DOFs and the analysis time. Also, you can do the same investigation but by changing the mesh density.

Regarding the implementation in MoFEM, it is important that you get the concept of developing/using UDOs to evaluate the matrices and vectors and then push them, in a certain order, to the Pipelines where calculations are done sequentially. These concepts will apply to all of the programs implemented in MoFEM.

You can find the complete code to solve 3D version of the Poisson problem with homogeneous boundary condition, using the same source code by setting variable EXECUTABLE_DIMENSION=3, see CMakeList.txt file in mofem_install/um_view/tutorials/scl-1/

tutorials_build_and_install(
poisson_3d_homogeneous ${CMAKE_CURRENT_SOURCE_DIR}/poisson_2d_homogeneous.cpp)
set_target_properties(
poisson_3d_homogeneous PROPERTIES COMPILE_FLAGS
"-DEXECUTABLE_DIMENSION=3")

To run the analysis, you will follow very similar procedure as for the 2D case using the same param_file.petsc file with slight changes in the command lines

mofem_part -my_file cube.h5m -output_file cube_2parts.h5m -my_nparts 2 -dim 3 -adj_dim 1
mpirun -np 2 ./poisson_3d_homogeneous -file_name cube_2parts.h5m -order 4 -log_quiet

Plain program

The plain program for both the implementation of the UDOs (*.hpp) and the main program (*.cpp) are as follows

Implementation of User Data Operators (*.hpp)

/**
* \file poisson_2d_homogeneous.hpp
* \example poisson_2d_homogeneous.hpp
*
* Solution of poisson equation. Direct implementation of User Data Operators
* for teaching proposes.
*
* \note In practical application we suggest use form integrators to generalise
* and simplify code. However, here we like to expose user to ways how to
* implement data operator from scratch.
*/
// Define name if it has not been defined yet
#ifndef __POISSON_2D_HOMOGENEOUS_HPP__
#define __POISSON_2D_HOMOGENEOUS_HPP__
// Use of alias for some specific functions
// We are solving Poisson's equation in 2D so Face element is used
template <int DIM>
using ElementsAndOps = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>;
FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
// Namespace that contains necessary UDOs, will be included in the main program
// Declare FTensor index for 2D problem
// For simplicity, source term f will be constant throughout the domain
const double body_source = 5.;
struct OpDomainLhsMatrixK : public AssemblyDomainEleOp {
public:
OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)
: AssemblyDomainEleOp(row_field_name, col_field_name,
DomainEleOp::OPROWCOL) {}
const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
this->locMat.resize(nb_row_dofs, nb_col_dofs, false);
this->locMat.clear();
// get element area
const double area = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on row
auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
// get derivatives of base functions on column
auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
for (int cc = 0; cc != nb_col_dofs; cc++) {
this->locMat(rr, cc) += t_row_diff_base(i) * t_col_diff_base(i) * a;
// move to the derivatives of the next base functions on column
++t_col_diff_base;
}
// move to the derivatives of the next base functions on row
++t_row_diff_base;
}
// move to the weight of the next integration point
++t_w;
}
}
};
struct OpDomainRhsVectorF : public AssemblyDomainEleOp {
public:
const int nb_dofs = data.getIndices().size();
this->locF.resize(nb_dofs, false);
this->locF.clear();
// get element area
const double area = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get base function
auto t_base = data.getFTensor0N();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_dofs; rr++) {
this->locF[rr] += t_base * body_source * a;
// move to the next base function
++t_base;
}
// move to the weight of the next integration point
++t_w;
}
}
};
}; // namespace Poisson2DHomogeneousOperators
#endif //__POISSON_2D_HOMOGENEOUS_HPP__

Implementation of the main program (*.cpp)

/**
* \file poisson_2d_homogeneous.cpp
* \example poisson_2d_homogeneous.cpp
*
* Solution of poisson equation. Direct implementation of User Data Operators
* for teaching proposes.
*
* \note In practical application we suggest use form integrators to generalise
* and simplify code. However, here we like to expose user to ways how to
* implement data operator from scratch.
*/
constexpr auto field_name = "U";
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
using namespace MoFEM;
static char help[] = "...\n\n";
public:
// Declaration of the main function to run analysis
MoFEMErrorCode runProgram();
private:
// Declaration of other main functions called in runProgram()
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode setIntegrationRules();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
// MoFEM interfaces
Simple *simpleInterface;
// Field name and approximation order
int oRder;
};
: mField(m_field) {}
//! [Read mesh]
}
//! [Read mesh]
//! [Setup problem]
int oRder = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
}
//! [Setup problem]
//! [Boundary condition]
auto bc_mng = mField.getInterface<BcManager>();
// Remove BCs from blockset name "BOUNDARY_CONDITION" or SETU, note that you
// can use regular expression to put list of blocksets;
CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
simpleInterface->getProblemName(), "(BOUNDARY_CONDITION|SETU)",
std::string(field_name), true);
}
//! [Boundary condition]
//! [Assemble system]
auto pipeline_mng = mField.getInterface<PipelineManager>();
{ // Push operators to the Pipeline that is responsible for calculating LHS
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating RHS
auto set_values_to_bc_dofs = [&](auto &fe) {
auto get_bc_hook = [&]() {
return hook;
};
fe->preProcessHook = get_bc_hook();
};
// you can skip that if boundary condition is prescribing zero
auto calculate_residual_from_set_values_on_bc = [&](auto &pipeline) {
using DomainEle =
auto grad_u_vals_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
grad_u_vals_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInternal(field_name, grad_u_vals_ptr,
[](double, double, double) constexpr { return -1; }));
};
pipeline_mng->getOpDomainRhsPipeline(), {H1});
set_values_to_bc_dofs(pipeline_mng->getDomainRhsFE());
calculate_residual_from_set_values_on_bc(
pipeline_mng->getOpDomainRhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
}
//! [Assemble system]
//! [Set integration rules]
auto rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
auto rule_rhs = [](int, int, int p) -> int { return p; };
auto pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
}
//! [Set integration rules]
//! [Solve system]
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto ksp_solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(ksp_solver);
CHKERR KSPSetUp(ksp_solver);
// Create RHS and solution vectors
auto dm = simpleInterface->getDM();
auto F = createDMVector(dm);
auto D = vectorDuplicate(F);
// Solve the system
CHKERR KSPSolve(ksp_solver, F, D);
// Scatter result data on the mesh
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Solve system]
//! [Output results]
auto pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
post_proc_fe->getOpPtrVector(), {H1});
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{{"U", u_ptr}},
OpPPMap::DataMapMat{{"GRAD_U", grad_u_ptr}},
)
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile("out_result.h5m");
}
//! [Output results]
//! [Run program]
}
//! [Run program]
//! [Main]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Error handling
try {
// Register MoFEM discrete manager in PETSc
DMType dm_name = "DMMOFEM";
// Create MOAB instance
moab::Core mb_instance; // mesh database
moab::Interface &moab = mb_instance; // mesh database interface
// Create MoFEM instance
MoFEM::Core core(moab); // finite element database
MoFEM::Interface &m_field = core; // finite element interface
// Run the main analysis
Poisson2DHomogeneous poisson_problem(m_field);
CHKERR poisson_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
//! [Main]
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
Poisson2DHomogeneousOperators::OpDomainRhsVectorF::OpDomainRhsVectorF
OpDomainRhsVectorF(std::string field_name)
Definition: poisson_2d_homogeneous.hpp:94
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
Poisson2DHomogeneous::solveSystem
MoFEMErrorCode solveSystem()
[Set integration rules]
Definition: poisson_2d_homogeneous.cpp:169
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
MoFEM::PipelineManager::getDomainRhsFE
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: PipelineManager.hpp:405
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
surface
Definition: surface.py:1
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition: FormsIntegrators.hpp:247
Poisson2DHomogeneousOperators::OpDomainRhsVectorF
Definition: poisson_2d_homogeneous.hpp:92
Poisson2DHomogeneous::oRder
int oRder
Definition: poisson_2d_homogeneous.cpp:49
Poisson2DHomogeneousOperators::OpDomainLhsMatrixK::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: poisson_2d_homogeneous.hpp:44
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
Poisson2DHomogeneous::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: poisson_2d_homogeneous.cpp:68
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProcBrokenMeshInMoabBase.hpp:700
MoFEM::PipelineManager::loopFiniteElements
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Definition: PipelineManager.cpp:19
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::PipelineManager::getOpDomainRhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Definition: PipelineManager.hpp:773
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
Poisson2DHomogeneous::assembleSystem
MoFEMErrorCode assembleSystem()
[Boundary condition]
Definition: poisson_2d_homogeneous.cpp:101
Poisson2DHomogeneousOperators::i
FTensor::Index< 'i', SPACE_DIM > i
Definition: poisson_2d_homogeneous.hpp:33
Poisson2DHomogeneous::outputResults
MoFEMErrorCode outputResults()
[Solve system]
Definition: poisson_2d_homogeneous.cpp:196
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Poisson2DHomogeneous::Poisson2DHomogeneous
Poisson2DHomogeneous(MoFEM::Interface &m_field)
Definition: poisson_2d_homogeneous.cpp:52
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
Poisson2DHomogeneous::simpleInterface
Simple * simpleInterface
Definition: poisson_2d_homogeneous.cpp:46
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
EigenMatrix::Number
FTensor::Number< N > Number
Definition: MatrixFunctionTemplate.hpp:55
UNKNOWNNAME
@ UNKNOWNNAME
Definition: definitions.h:158
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::OpBaseImpl::locF
VectorDouble locF
local entity vector
Definition: FormsIntegrators.hpp:249
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
Poisson2DHomogeneousOperators::OpDomainLhsMatrixK
Definition: poisson_2d_homogeneous.hpp:38
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1067
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
convert.type
type
Definition: convert.py:64
MoFEM::OpBaseImpl::iNtegrate
virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.
Definition: FormsIntegrators.hpp:257
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:310
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::BcScalarMeshsetType
Definition: BcManager.hpp:15
Poisson2DHomogeneousOperators::body_source
const double body_source
Definition: poisson_2d_homogeneous.hpp:36
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(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_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
MoFEM::PipelineManager::createKSP
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
Definition: PipelineManager.cpp:49
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
Poisson2DHomogeneousOperators::OpDomainLhsMatrixK::OpDomainLhsMatrixK
OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)
Definition: poisson_2d_homogeneous.hpp:40
Poisson2DHomogeneous
Definition: poisson_2d_homogeneous.cpp:27
MoFEM::PipelineManager::setDomainRhsIntegrationRule
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:530
Poisson2DHomogeneous::mField
MoFEM::Interface & mField
Definition: poisson_2d_homogeneous.cpp:45
Poisson2DHomogeneous::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
[Assemble system]
Definition: poisson_2d_homogeneous.cpp:154
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: PipelineManager.hpp:749
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:701
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', SPACE_DIM >
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
Poisson2DHomogeneous::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Setup problem]
Definition: poisson_2d_homogeneous.cpp:85
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
MoFEM::PipelineManager::setDomainLhsIntegrationRule
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:503
MoFEM::EssentialPreProc< TemperatureCubitBcData >
Specialization for TemperatureCubitBcData.
Definition: EssentialTemperatureCubitBcData.hpp:91
ElementsAndOps
Definition: child_and_parent.cpp:18
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
Poisson2DHomogeneous::runProgram
MoFEMErrorCode runProgram()
[Output results]
Definition: poisson_2d_homogeneous.cpp:242
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
Poisson2DHomogeneousOperators
Definition: poisson_2d_homogeneous.hpp:30
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
poisson_2d_homogeneous.hpp
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
Poisson2DHomogeneous::readMesh
MoFEMErrorCode readMesh()
[Read mesh]
Definition: poisson_2d_homogeneous.cpp:56
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:362
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
Poisson2DHomogeneousOperators::OpDomainRhsVectorF::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Definition: poisson_2d_homogeneous.hpp:97
EshelbianPlasticity::U
@ U
Definition: EshelbianContact.cpp:193
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
F
@ F
Definition: free_surface.cpp:394
MoFEM::PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:677
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698