v0.13.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.

In this practical example, that actually solves something using MoFEM, we will solve a simple Poisson's equation in 2D with homogeneous boundary condition (zero boundary values). 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 \( {x-y} \) space of a point in the domain.

For a Poisson problem, on rectangle domains, 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 the field \( u \) of the problem equation, we need to derive the weak form. 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.

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 KU} = {\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 order 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.

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

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

- Poisson2DHomogeneousOperators::OpDomainLhsMatrixK is responsible for the calculation and assembly of the left-hand-side matrix \( {\bf K} \)
- Poisson2DHomogeneousOperators::OpDomainRhsVectorF is responsible for the 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 of the file that includes some libraries for finite element implementation, some aliases, and declaration/initialisation as commented below.

// 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

#include <stdlib.h>

#include <BasicFiniteElements.hpp>

// Use of alias for some specific functions

// We are solving Poisson's equation in 2D so Face element is used

// Namespace that contains necessary UDOs, will be included in the main program

namespace Poisson2DHomogeneousOperators {

// Declare FTensor index for 2D problem

// For simplicity, source term f will be constant throughout the domain

// Implementation of the UDOs below this point

}

ForcesAndSourcesCore::UserDataOperator UserDataOperator

FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore

Face finite element default.

const double body_source

Data on single entity (This is passed as argument to DataOperator::doWork)

default operator for TRI element

Now, we will look in details how the two main operators Poisson2DHomogeneousOperators::OpDomainLhsMatrixK and Poisson2DHomogeneousOperators::OpDomainRhsVectorF are implemented.

The first UDO Poisson2DHomogeneousOperators::OpDomainLhsMatrixK is responsible for the calculation and assembly of the left-hand-side matrix \( {\bf K} \)

First, let's look at the structure of the UDO to calculate matrix \( {\bf K} \), also the general structure of all of the UDOs implemented in MoFEM

public:

OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)

sYmm = true;

}

EntityType col_type, EntData &row_data,

EntData &col_data) {

// Implementation of doWork() below this point

}

private:

MatrixDouble locLhs, transLocLhs;

};

virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)

Operator for bi-linear form, usually to calculate values on left hand side.

We have the class `OpDomainLhsMatrixK`

that is inherited from `OpFaceEle`

which is the alias of the base class MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator. This newly created class has two public objects `OpDomainLhsMatrixK()`

and `doWork()`

which can be accessed from outside of the class, and two private objects `locLhs`

and `transLocLhs`

of type `MatrixDouble`

which can be accessed only from inside of the class for calculation and assembly of local stiffness matrix. 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()`

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. In this particular case of Poisson's problem with homogeneous boundary condition, we will see later, they take the same field name of "U" for both row and column. OpFaceEle::OPROWCOL because we assemble stiffness matrix that has data on both row and column. The booleanOpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)sYmm = true;}`sYmm = true`

is set since we know matrix \( K \) is symmetric and we can just assemble and store data of half of the matrix.`doWork()`

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 of the mathematical derivations that we did in the previous section of this tutorial will be implemented in this function. More discussion on this essential function will be given below.`locLhs`

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.`transLocLhs`

This private member object also of type MatrixDouble is used to assemble the transpose of`locLhs`

. This object is used only in function`doWork()`

so it is also declared as a private object of the class`OpDomainLhsMatrixK`

.

Now we take a closer look on the details of the implementation of the member function `doWork()`

which is the most essential part of the UDO implementation. Full code for this function as follows

EntityType col_type, EntData &row_data,

EntData &col_data) {

if (nb_row_dofs && nb_col_dofs) {

locLhs.resize(nb_row_dofs, nb_col_dofs, false);

locLhs.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<2>();

// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX

for (int gg = 0; gg != nb_integration_points; gg++) {

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++) {

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

}

// FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX

CHKERR MatSetValues(getKSPB(), row_data, col_data, &locLhs(0, 0),

ADD_VALUES);

if (row_side != col_side || row_type != col_type) {

transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);

noalias(transLocLhs) = trans(locLhs);

CHKERR MatSetValues(getKSPB(), col_data, row_data, &transLocLhs(0, 0),

ADD_VALUES);

}

}

}

#define MoFEMFunctionBegin

First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...

#define MoFEMFunctionReturn(a)

Last executable line of each PETSc function used for error handling. Replaces return()

constexpr auto MatSetValues

FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)

Get derivatives of base functions.

const VectorInt & getIndices() const

Get global indices of dofs on entity.

First, let's talk about the structure of this function

You will notice that all functions written in MoFEM will include error handling as we can see here with keywords: `MoFEMErrorCode`

, `MoFEMFunctionBegin;`

, and `MoFEMFunctionReturn(0);`

Regarding arguments of this `doWork()`

function, it has six arguments and all of them are input arguments. The values of these input arguments are provided by MoFEM based on the type of the elements we are solving (edge, face, volume). In this particular case, as we solve 2D Poisson problem using triangular mesh, we use MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator (alias `OpFaceEle`

) as seen previously at the first lines of the code for the struct `OpDomainLhsMatrixK`

. And the data structure of the entities (elements) are provided accordingly thanks to EntitiesFieldData::EntData (alias `EntData`

).

The arguments are divided in three groups

`row_side`

and`col_side`

takes integer values of {0, 1, 2}`row_type`

and`col_type`

takes value of {MBVERTEX, MBEDGE, MBTRI} (see more Hierarchical basis functions)`&row_data`

and`&col_data`

are two important pointers that have information of the DOFs, and basis function values, field values, field gradients at integration (Gauss) points.

Moving into the main implementation of `doWork()`

function, we will see information about the number of DOFs on row and column are extracted from the database

Then there is an `if statement`

to make sure that data for row and column are both valid.

if (nb_row_dofs && nb_col_dofs) {

// Main implementation goes here

}

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

locLhs.resize(nb_row_dofs, nb_col_dofs, false);

locLhs.clear();

Then you will get the area 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++) {

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++) {

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

where `a`

is the intermediate variable of type `double`

and calculated earlier as

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

- Loop over all integration points

However, finite element implementation in MoFEM is slightly different as you may have already noticed. MoFEM still has the nested loops; however, in the implementation of UDO, you have only one loop over the integration points inside which you calculate the local stiffness matrix components and assemble them to the global matrix. The loop over the elements (or `entities`

in the context of Hierarchical basis functions used in MoFEM) is done in the main program, outside of the UDO. Details of how this loop is triggered is presented in solveSystem().

That was the description for the part of the code responsible for the calculation of the components of the element stiffness matrix. Now we will talk about the last part of the current UDO which is about the assembling/filling from local (element) stiffness matrix to global stiffness matrix handled by PETSc. This is done by the following lines of code

// FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX

CHKERR MatSetValues(getKSPB(), row_data, col_data, &locLhs(0, 0),

ADD_VALUES);

This process of assembling matrix is done using the function `MatSetValues()`

which set the values for the global matrix and later will be solved using PETSc. This function requires the inputs of `getKSPB()`

, `row_data`

, `col_data`

, the memory address where the first components of the local stiffness matrix `&locLhs(0, 0)`

is stored, and a flag to add value `ADD_VALUES`

. Here `getKSPB()`

is the method to set the value of LHS matrix that will be solved by the iterative solver namely KSP in PETSc.

And finally, as global \( {\bf K} \) matrix is also symmetric, we transpose everything that has already assembled to form the full matrix of \( {\bf K} \). This process is implemented and visualised as follows

if (row_side != col_side || row_type != col_type) {

transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);

noalias(transLocLhs) = trans(locLhs);

CHKERR MatSetValues(getKSPB(), col_data, row_data, &transLocLhs(0, 0),

ADD_VALUES);

}

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

public:

OpDomainRhsVectorF(std::string field_name)

// Implementation of doWork() below this point

}

private:

VectorDouble locRhs;

};

We have the public objects of `OpDomainRhsVectorF()`

and `doWork()`

and private object of `locRhs`

. 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 `OpFaceEle::OPROW`

. Similarly for function `doWork()`

, 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 `doWork()`

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

if (nb_dofs) {

locRhs.resize(nb_dofs, false);

locRhs.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++) {

for (int rr = 0; rr != nb_dofs; rr++) {

locRhs[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;

}

// FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR

// Ignoring DOFs on boundary (index -1)

CHKERR VecSetOption(getKSPf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);

CHKERR VecSetValues(getKSPf(), data, &locRhs(0), ADD_VALUES);

}

}

constexpr auto VecSetValues

FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)

Get base function as Tensor0.

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++) {

locRhs[rr] += t_base * body_source * a;

// move to the next base function

++t_base;

}

Similarly to the implementation of matrix UDO, the final part of the code involves filling/assembling the calculated values from local force vector to the global one. This procedure for force vector is slightly different from the matrix. It is done in two steps. First, we need to tell the program to ignore the assembly of DOFs associated with boundary entities (marked as -1, show where and how to mark those DOFs) and assemble only DOFs associated with the domain entities.

CHKERR VecSetOption(getKSPf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);

And then assemble the local force vector to the global one using

CHKERR VecSetValues(getKSPf(), data, &locRhs(0), ADD_VALUES);

This is very similar to the previous case for matrix, but now `getKSPf()`

is used to set the value of the RHS vector later solved by KSP iterative solver in PETSc.

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.

This main class Poisson2DHomogeneous contains functions and each of which is responsible for a certain task of a finite element program.

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

public:

Poisson2DHomogeneous(MoFEM::Interface &m_field);

// Declaration of the main function to run analysis

MoFEMErrorCode runProgram();

Deprecated interface functions.

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

MoFEM::Interface &mField;

Simple *simpleInterface;

// Field name and approximation order

std::string domainField;

int oRder;

- Note
- For writhing program in MoFEM, we follow some coding practices for code style and naming convention described at Coding practice

Now is the constructor

: domainField("U"), mField(m_field) {}

Poisson2DHomogeneous(MoFEM::Interface &m_field)

The constructor specifies

`domainField("U")`

: The domain field name as U (the field that we need to solve for solution)`mField(m_field)`

: The MoFEM instance that is the backbone of the program

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

CHKERR simpleInterface->getOptions();

CHKERR simpleInterface->loadFile();

}

MoFEMErrorCode getInterface(IFACE *&iface) const

Get interface refernce to pointer of interface.

Simple * simpleInterface

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.

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

int oRder = 3;

CHKERR simpleInterface->setUp();

}

PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)

This function has

`simpleInterface->addDomainField`

: Add the domain field through the`simpleInterface`

. We need to add the domain field 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 (
`domainField`

), - 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.

- Field name (
`simpleInterface->setFieldOrder`

: Set the polynomial order of the approximation of the field`simpleInterface->setUp()`

: Finally, do the setup of the problem.

This function helps to deal with the boundary condition and its implementation for the current problem is as follows

// Get boundary edges marked in block named "BOUNDARY_CONDITION"

Range boundary_entities;

std::string entity_name = it->getName();

if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {

boundary_entities, true);

}

}

// Add vertices to boundary entities

Range boundary_vertices;

boundary_vertices, true);

boundary_entities.merge(boundary_vertices);

// Remove DOFs as homogeneous boundary condition is used

simpleInterface->getProblemName(), domainField, boundary_entities);

}

#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)

Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.

virtual moab::Interface & get_moab()=0

There is two parts of the implementation for the boundary condition for the current problem of Poisson's equation with homogeneous boundary condition

- First, identify all the entities, from the input mesh, to which you would like to apply boundary condition. This is done in two steps
- Getting the edges (up to faces for 3D problems) in the block that has been marked with a certain name, e.g.
*BOUNDARY_CONDITION*, in the input mesh. So essentially, when creating the mesh, you select the edges where you want to apply boundary condition, you put those edges into one block and name that block making sure that the name in the mesh and in the code are matching. See more on how to create the mesh having that properties from Gmsh at MSH-1: Create a 2D mesh from Gmsh. The part of the code that is responsible for this task is as follows// Get boundary edges marked in block named "BOUNDARY_CONDITION"Range boundary_entities;std::string entity_name = it->getName();if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {boundary_entities, true);}} - Including the vertices that are shared between the marked edges and the interior edges and faces to the set of entities on which boundary condition applies. We may normally forget this step but it is essential to have those vertices included, making sure the implementation of boundary condition is correct and avoids errors in results. The inclusion is done in this part of the code

- Getting the edges (up to faces for 3D problems) in the block that has been marked with a certain name, e.g.
- Second, set the value of the boundary condition to the identified entities, i.e. (homogeneous) zero value boundary condition in this case. This process is simply done using the
`removeDofsOnEntities()`

function like thisOf course, applying boundary condition for non-homogeneous boundary condition will not be this easy. We will show ways how to do it in MoFEM in other tutorials.// Remove DOFs as homogeneous boundary condition is usedCHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(simpleInterface->getProblemName(), domainField, boundary_entities);

This part is about the function that is responsible for the assembling of the system of equations. As a reminder, at the beginning of the section Implementation, we were mentioning *two important concepts* in MoFEM, namely **UDO** and **Pipiline**. 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, this function is also responsible for setting operators for KSP solver (iterative linear solver provided by PETSc) from the implemented pipelines. The full source code for the function is as follows

auto det_ptr = boost::make_shared<VectorDouble>();

auto jac_ptr = boost::make_shared<MatrixDouble>();

auto inv_jac_ptr = boost::make_shared<MatrixDouble>();

{ // Push operators to the Pipeline that is responsible for calculating LHS

pipeline_mng->getOpDomainLhsPipeline().push_back(

new OpCalculateHOJacForFace(jac_ptr));

pipeline_mng->getOpDomainLhsPipeline().push_back(

new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));

pipeline_mng->getOpDomainLhsPipeline().push_back(

new OpSetInvJacH1ForFace(inv_jac_ptr));

pipeline_mng->getOpDomainLhsPipeline().push_back(

}

{ // Push operators to the Pipeline that is responsible for calculating RHS

pipeline_mng->getOpDomainRhsPipeline().push_back(

new OpDomainRhsVectorF(domainField));

}

}

OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace

- First, the pipeline manager (
`pipeline_mng`

) that manages two main pipelines are created. These two main pipelines are`LhsPipeline`

: responsible for calculations of the left hand side matrix`RhsPipeline`

: responsible for calculations of the right hand side vectorauto pipeline_mng = mField.getInterface<PipelineManager>();

- Second, pushing the operators to the Pipeline that is responsible for the calculation of the LHS matrix. This is 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.It is noted that the procedure to push operators calculating the inverse of Jacobian with the presence of the gradient of approximation function is needed only for 2D problem. For 3D problem, it is done automatically and user does not have to manually add codes to push operators calculating the inverse of Jacobian. This is because MoFEM was initially designed to solve 3D problems and it made some adjustments along the way for solving 2D problems.{ // Push operators to the Pipeline that is responsible for calculating LHSauto det_ptr = boost::make_shared<VectorDouble>();auto jac_ptr = boost::make_shared<MatrixDouble>();auto inv_jac_ptr = boost::make_shared<MatrixDouble>();pipeline_mng->getOpDomainLhsPipeline().push_back(new OpCalculateHOJacForFace(jac_ptr));pipeline_mng->getOpDomainLhsPipeline().push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));pipeline_mng->getOpDomainLhsPipeline().push_back(new OpDomainLhsMatrixK(domainField, domainField));} - Third, pushing the implemented UDO to Pipeline that is responsible for the calculation of the RHS vector Here you can see, we do not need derivative of the approximation function to calculate force vector in UDO{ // Push operators to the Pipeline that is responsible for calculating RHSpipeline_mng->getOpDomainRhsPipeline().push_back(new OpDomainRhsVectorF(domainField));}
`OpDomainRhsVectorF`

. Consequently, the operators for calculating inverse of Jacobian is not needed as well.

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

integration rules

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.

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 ksp_solver = pipeline_mng->createKSP();

CHKERR KSPSetFromOptions(ksp_solver);

CHKERR KSPSetUp(ksp_solver);

// Create RHS and solution vectors

auto dm = simpleInterface->getDM();

auto F = smartCreateDMVector(dm);

// Solve the system

// Scatter result data on the mesh

}

PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)

set local (or ghosted) vector values on mesh for partition only

SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)

Create duplicate vector of smart vector.

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 = smartCreateDMVector(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

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

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

pipeline_mng->getDomainLhsFE().reset();

auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);

post_proc_fe->generateReferenceElementMesh();

post_proc_fe->addFieldValuesPostProc(domainField);

pipeline_mng->getDomainRhsFE() = post_proc_fe;

CHKERR pipeline_mng->loopFiniteElements();

CHKERR post_proc_fe->writeFile("out_result.h5m");

}

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

}

This `main()`

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

// Initialisation of MoFEM/PETSc and MOAB data structures

// Error handling

try {

// Register MoFEM discrete manager in PETSc

DMType dm_name = "DMMOFEM";

CHKERR DMRegister_MoFEM(dm_name);

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

}

PetscErrorCode DMRegister_MoFEM(const char sname[])

Register MoFEM problem.

static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])

Initializes the MoFEM database PETSc, MOAB and MPI.

static MoFEMErrorCode Finalize()

Checks for options to be called at the conclusion of the program.

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

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, 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

- For user version installation
- 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 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.

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
- 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 12682136550675316745Add finite element dFE
- General information about the problem including removing DOFs associated with the boundary Add problem SimpleProblemNumber of dofs 1687Number of dofs 1659Finite element dFE added. Nb. of elements added 205Finite element dFE added. Nb. of elements added 201Number of adjacencies 1435SimpleProblem Nb. local dof 1687 by 1687 nb global dofs 3289 by 3289SimpleProblem Nb. local dof 1602 by 1602 nb global dofs 3289 by 3289FEs ghost dofs on problem SimpleProblem Nb. ghost dof 0 by 0 Nb. local dof 1687 by 1687FEs ghost dofs on problem SimpleProblem Nb. ghost dof 57 by 57 Nb. local dof 1602 by 1602removed 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-011 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

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.

With the implementation for 2D problem done, you can effortlessly extend it to solve 3D problem with few changes as follows

- Change all
`Face`

to`Volume`

(`.cpp`

and`.hpp`

files)- Using
`PostProcVolumeOnRefinedMesh`

instead of`PostProcFaceOnRefinedMesh`

- Using
`VolumeElementForcesAndSourcesCore`

instead of`FaceElementForcesAndSourcesCore`

- Using
- Change FTensore objects that have two components to three (
`.hpp`

file)- Using
*FTensor::Index<'i', 3> i;*instead of*FTensor::Index<'i', 2> i;* - Using
*getFTensor1DiffN<3>*instead of*getFTensor1DiffN<2>*

- Using
- No need to specifically calculate inverse of Jacobian in 3D case, it is done automatically behind the curtain (
`.cpp`

file)auto det_ptr = boost::make_shared<VectorDouble>();auto jac_ptr = boost::make_shared<MatrixDouble>();auto inv_jac_ptr = boost::make_shared<MatrixDouble>();pipelineLhs->getOpPtrVector().push_back(new OpCalculateHOJacForFace(jac_ptr));pipelineLhs->getOpPtrVector().push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr)); - [Optional] Objects' name changed for readability (
`.cpp`

and`.hpp`

file)- Using
`Poisson3DHomogeneous`

instead of`Poisson2DHomogeneous`

- Using
`Poisson3DHomogeneousOperators`

instead of`Poisson2DHomogeneousOperators`

- Using variable
`volume`

instead of`area`

- Using

You can find the complete code to solve 3D version of the Poisson problem with homogeneous boundary condition in poisson_3d_homogeneous.hpp and poisson_3d_homogeneous.cpp located in the same path with the source codes for 2D problem `mofem_install/um_view/tutorials/scl-1/`

.

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

The plain program for both the implementation of the UDOs (`*`

.hpp) and the main program (`*`

.cpp) are as follows

// 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

#include <stdlib.h>

#include <BasicFiniteElements.hpp>

// Use of alias for some specific functions

// We are solving Poisson's equation in 2D so Face element is used

// Namespace that contains necessary UDOs, will be included in the main program

namespace Poisson2DHomogeneousOperators {

// Declare FTensor index for 2D problem

// For simplicity, source term f will be constant throughout the domain

public:

OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)

sYmm = true;

}

EntityType col_type, EntData &row_data,

EntData &col_data) {

if (nb_row_dofs && nb_col_dofs) {

locLhs.resize(nb_row_dofs, nb_col_dofs, false);

locLhs.clear();

// get element area

// get number of integration points

// get integration weights

auto t_w = getFTensor0IntegrationWeight();

// 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++) {

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++) {

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

}

// FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX

ADD_VALUES);

if (row_side != col_side || row_type != col_type) {

transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);

noalias(transLocLhs) = trans(locLhs);

ADD_VALUES);

}

}

}

private:

};

public:

OpDomainRhsVectorF(std::string field_name)

if (nb_dofs) {

locRhs.resize(nb_dofs, false);

locRhs.clear();

// get element area

// get number of integration points

// 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++) {

for (int rr = 0; rr != nb_dofs; rr++) {

// move to the next base function

++t_base;

}

// move to the weight of the next integration point

++t_w;

}

// FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR

// Ignoring DOFs on boundary (index -1)

}

}

private:

};

}; // namespace Poisson2DHomogeneousOperators

#endif //__POISSON2DHOMOGENEOUS_HPP__

Vec getKSPf() const

auto getFTensor0IntegrationWeight()

Get integration weights.

@ OPROWCOL

Mat getKSPB() const

MatrixDouble & getGaussPts()

matrix of integration (Gauss) points for Volume Element

MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)

Operator for bi-linear form, usually to calculate values on left hand side.

MatrixDouble transLocLhs

MatrixDouble locLhs

OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)

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

Operator for linear form, usually to calculate values on right hand side.

VectorDouble locRhs

OpDomainRhsVectorF(std::string field_name)

#include <stdlib.h>

#include <BasicFiniteElements.hpp>

#include <poisson_2d_homogeneous.hpp>

using namespace MoFEM;

using namespace Poisson2DHomogeneousOperators;

struct Poisson2DHomogeneous {

public:

Poisson2DHomogeneous(MoFEM::Interface &m_field);

// 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

MoFEM::Interface &mField;

Simple *simpleInterface;

// Field name and approximation order

std::string domainField;

int oRder;

};

: domainField("U"), mField(m_field) {}

//! [Read mesh]

CHKERR simpleInterface->getOptions();

CHKERR simpleInterface->loadFile();

}

//! [Read mesh]

//! [Setup problem]

int oRder = 3;

CHKERR simpleInterface->setUp();

}

//! [Setup problem]

//! [Boundary condition]

// Get boundary edges marked in block named "BOUNDARY_CONDITION"

Range boundary_entities;

std::string entity_name = it->getName();

if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {

boundary_entities, true);

}

}

// Add vertices to boundary entities

Range boundary_vertices;

boundary_vertices, true);

boundary_entities.merge(boundary_vertices);

// Remove DOFs as homogeneous boundary condition is used

simpleInterface->getProblemName(), domainField, boundary_entities);

}

//! [Boundary condition]

//! [Assemble system]

auto det_ptr = boost::make_shared<VectorDouble>();

auto jac_ptr = boost::make_shared<MatrixDouble>();

auto inv_jac_ptr = boost::make_shared<MatrixDouble>();

{ // Push operators to the Pipeline that is responsible for calculating LHS

pipeline_mng->getOpDomainLhsPipeline().push_back(

new OpCalculateHOJacForFace(jac_ptr));

pipeline_mng->getOpDomainLhsPipeline().push_back(

new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));

pipeline_mng->getOpDomainLhsPipeline().push_back(

new OpSetInvJacH1ForFace(inv_jac_ptr));

pipeline_mng->getOpDomainLhsPipeline().push_back(

}

{ // Push operators to the Pipeline that is responsible for calculating RHS

pipeline_mng->getOpDomainRhsPipeline().push_back(

new OpDomainRhsVectorF(domainField));

}

}

//! [Assemble system]

//! [Set integration rules]

CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);

CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);

}

//! [Set integration rules]

//! [Solve system]

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 = smartCreateDMVector(dm);

// Solve the system

// Scatter result data on the mesh

}

//! [Solve system]

//! [Output results]

pipeline_mng->getDomainLhsFE().reset();

auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);

post_proc_fe->generateReferenceElementMesh();

post_proc_fe->addFieldValuesPostProc(domainField);

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]

// Initialisation of MoFEM/PETSc and MOAB data structures

// Error handling

try {

// Register MoFEM discrete manager in PETSc

DMType dm_name = "DMMOFEM";

CHKERR DMRegister_MoFEM(dm_name);

// 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]

Generated on Thu May 5 2022 22:31:06 for MoFEM by Doxygen 1.9.1 and hosted at