v0.14.0
SCL-2: Poisson's equation (non-homogeneous BC)
Note
Prerequisite of this tutorial is SCL-1: Poisson's equation (homogeneous BC)

Note
Intended learning outcome:
• idea of Boundary element in MoFEM and how to use it
• handling of non-homogeneous boundary conditions using least square approximation in MoFEM
• process of implementing additional User Data Operators (UDOs) to calculate and assemble non-homogeneous boundary condition
• practising on how to push the developed UDOs to the Pipeline

# Introduction

The focus of our current tutorial is on solving Laplace-type partial differential equations (PDEs) Poisson Problem with non-homogenous Dirichlet boundary conditions. While the SCL-2: Poisson's equation (non-homogeneous BC) has designed for solving Dirichlet problems with homogeneous boundary conditions, this tutorial guides the methods to extend this to non-homogeneous cases and plot the solution.

We are also aiming to present MoFEM approaches for handling Dirichlet problems in addition to the least square approximation while establishing your necessary foundation in MoFEM as well. Prior to that, we would like to suggest understanding the operators(i.e. OpDomainLhsMatrixK and OpDomainRhsVectorF) in SCL-1: Poisson's equation (homogeneous BC) before proceeding to this tutorial.

## Problem Statement

We are demonstrating an approach to solve a simple Poisson's problem with nonhomogenous BC and help you grasp the approach better. For this problem, we denote the domain as Ω, and Γ represents the Dirichlet boundary, as illustrated in Figure 1.

Thus, we are solving a sample problem with nonhomogeneous BCs to find $$\mathbf{u}$$. such that:

Domain Equation:

\begin{align} -\nabla \cdot\nabla u = f(\mathbf{x}) \,\,\, in \,\,\, \Omega \end{align}

Boundary condition:

\begin{align} u = g(\mathbf{x}) \,\,\, on \,\,\,\Gamma \end{align}

Where f is the source term of the body and g is denoted as boundary function. However, if g is a constant everywhere (i.e. 0), we speak of homogeneous boundary conditions but in this problem, we are considering a spatially varying function $$g(\mathbf{x})$$.

Figure 1: Nonhomogeneous problem Statement.

# Numerical Scheme

In the context of solving this partial differential equation (PDE), Galerkin finite element methods with least square approximation have been executed to enforce nonhomogeneous boundary condition. We are also utilizing MoFEM libraries conforming finite element spaces with Least-Square approximation(LSA) for the boundary by minimizing the squared differences between the approximated values and the actual, similar to the following:

$\frac{1}{2}\int_{\Gamma} (u - g(\mathbf{x}))^2 \,d\Gamma$

Remember, previously we derived a weak form for homogeneous BC in SCL-1: Poisson's equation (homogeneous BC). Following that approch, we can write similar weak form for our current problem for domain (Ω) and boundary(Γ) which also exists for $$\forall \,\delta u$$ as follows:

\begin{align} \int_{\Omega} \nabla \delta u \cdot \nabla u \,d\Omega = \int_{\Omega} \delta u \,f(\mathbf{x})\,d\Omega \quad\text{and} \end{align}

\begin{align} \int_{\Gamma} \delta u \, \,d\Gamma = \int_{\Gamma} \delta u \,g(\mathbf{x})\,d\Gamma \end{align}

## Approximation

In the context of finite element method, the problem will be solved numerically by approximating the scalar space $$\mathbf{H^1}$$ by the discrete counterparts associated with our FE mesh of the computational domain $$\Omega$$. To do that, we are utilizing two key components to solve: (i) $$N_{\alpha}$$ and $$N_{\beta}$$ are the shape functions for the field variable and for test functions respectively, (ii) $$\bar{u}_{\alpha}$$ and $$\delta\bar{u}_{\beta}$$ represent for 'n' number of DoF within that element to obtain the solution 'u' as stated below:

\begin{align} u &= \sum_{\alpha=0}^{n-1} \mathbf{N}_{\alpha} \bar{u}_\alpha \,\,\text{and,}\\ {\delta u} &= \sum_{\beta=0}^{n-1}\mathbf{{N}_{\beta}} \delta\bar{u}_{\beta} \end{align}

Further, to assemble the element stiffness matrix we are indicating $$\alpha$$ for the collums and $$\beta$$ for rows. A typical approximation can be understood well in tutorial FUN-2: Hierarchical approximation. So, to approximate the fields based on the weak forms in equation (5) and (6) can be written as follows:

$$$\sum_{\alpha=0}^{n-1} \sum_{\beta=0}^{n-1} \delta\bar{u}_{\beta} \left[\int_{\Omega} \nabla \mathbf{N}_{\alpha} \cdot \nabla \mathbf{N}_{\beta} \,d\Omega\right] \bar{u}_{\alpha} \,\,= \,\, \sum_{\beta=0}^{n-1} \delta \bar{u}_{\beta} \left[\int_{\Omega} f(\mathbf{x})\,\, \mathbf{N}_{\beta} \,d\Omega\right] \quad ;\,\forall \,\delta \bar{u}_{\beta}$$$

and

$$$\sum_{\alpha=0}^{n-1} \sum_{\beta=0}^{n-1} \delta\bar{u}_{\beta} \left[\int_{\Gamma} \mathbf{N}_{\alpha} \,\,\mathbf{N}_{\beta} \,d\Gamma\right] \bar{u}_{\alpha} \,\,= \,\,\sum_{\beta=0}^{n-1} \delta \bar{u}_{\beta} \left[\int_{\Gamma} g(\mathbf{x})\,\, \mathbf{N}_{\beta} \,d\Gamma\right]\quad ;\forall \,\delta \bar{u}_{\beta}$$$

For a better understanding the implementation of the aboves in MoFEM, we are explaining the forms in two parts: domain and boundary. Also, let's say, we are naming the individuals for respective UDOs as follows,

$\mathbf{A_k} = \int_{\Omega} \nabla \mathbf{N}_{\alpha} \cdot \nabla \mathbf{N}_{\beta} \,d\Omega, \quad \mathbf{B_f} = \int_{\Omega} f(\mathbf{x})\,\,\mathbf{N}_{\beta} \,d\Omega\,\,\ \text{ |\,For Domain},$

and

$\mathbf{C_k} = \int_{\Gamma} \mathbf{N}_{\alpha} \cdot \mathbf{N}_{\beta} \,d\Gamma, \quad \mathbf{D_f} = \int_{\Gamma} g(\mathbf{x})\,\, \mathbf{N}_{\beta} \,d\Gamma \,\,\text{ |\,For Boundary}.$

Additionally, a similar discretization and implementation approach was referred in tutorial SCL-1: Poisson's equation (homogeneous BC).

## Problem formulation

Based on the above schemes, we are considering 4 operators in MoFEM. Let's say, we are naming those as follows:

Table 1: Operators for Domain and Boundary for Relevant Expresions
Respective Operators
$$\mathbf{A_k} \text{is implemented in } \text{OpDomainLhs}$$
$$\mathbf{B_f} \text{is implemented in } \text{OpDomainRhs}$$
$$\mathbf{C_k} \text{is implemented in } \text{OpBoundaryLhs}$$
$$\mathbf{D_f} \text{is implemented in } \text{OpBoundaryRhs}$$

The numerical scheme has several similarities with SCL-1: Poisson's equation (homogeneous BC) in terms of coding in MOFEM. Therefore, we are mostly going to discuss the critical parts of it in the next section comparing with the SCL-1: Poisson's equation (homogeneous BC).

# Code Structure

Before we explore the implementation details of the User-Defined Operators (UDOs) in "poisson_2d_nonhomogeneous.hpp", let's take a look at the initial lines of code in the file. It includes MoFEM libraries inclusions for finite element functionality and alias definitions for declaration and initialization as follows:

// Define file name if it has not been defined yet
#ifndef __POISSON2DNONHOMOGENEOUS_HPP__
#define __POISSON2DNONHOMOGENEOUS_HPP__
// Include standard library and Header file for basic finite elements
// implementation
#include <stdlib.h>
// 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
// Declare FTensor index for 2D problem
//declare a Lambda function
typedef boost::function<double(const double, const double, const double)>
//.. .. Implementation of the UDOs below this point .. ..

## User Defined Operators(UDO)

MoFEM has an in-built extensive library consisting of operators and to implement those operators users need to create some new UDOs according to the problems and PDEs. First, let's look at the structure of a UDO so that users can also learn the general structure of any UDOs implemented in MoFEM.

### OpDomainLhs()

For example, we can notice that our OpDomainLhs is similar to the OpDomainLhsMatrixK UDO in SCL-1: Poisson's equation (homogeneous BC). So, let's identify the differences bellow and then the further will be discussed:

struct OpDomainLhs : public OpFaceEle {
public:
OpDomainLhs(std::string row_field_name, std::string col_field_name)
: OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL) {
sYmm = true;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data) {
const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
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
// 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++) {
locLhs(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;
}
// FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
// Fill value to local stiffness matrix ignoring boundary DOFs
CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
getKSPB(), row_data, col_data, &locLhs(0, 0), ADD_VALUES);
// Fill values of symmetric local stiffness matrix
if (row_side != col_side || row_type != col_type) {
transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transLocLhs) = trans(locLhs);
CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
getKSPB(), col_data, row_data, &transLocLhs(0, 0), ADD_VALUES);
}
}
}
private:
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
MatrixDouble locLhs, transLocLhs;
};

In the above, the class OpDomainLhs is inherited from OpFaceEle, which stands for the alias for the base class MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator. Within this class, there are two public functions: OpDomainLhs() and doWork(). Here, the doWork() does the same as iNtegrate as in SCL-1: Poisson's equation (homogeneous BC). However, the difference is that we are using a different Lambda function to calculate the body source term. Alongside, we have got two private members named locLhs and transLocLhs. This follows the concept of data encapsulation to hide values and objects inside the class as much as possible.

• locLhs
This private member object of type MatrixDouble is used to store the results of the calculation of the components of element stiffness matrix. Besides, this object is made private so that only member of the class has access to use 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 OpDomainLhs.
• CHKERR MatSetValues
The code part is as follows:
CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
getKSPB(), row_data, col_data, &locLhs(0, 0), ADD_VALUES);
// Fill values of symmetric local stiffness matrix
if (row_side != col_side || row_type != col_type) {
transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transLocLhs) = trans(locLhs);
CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
getKSPB(), col_data, row_data, &transLocLhs(0, 0), ADD_VALUES);
The function MatSetValues() is used to assemble into the global matrix. This function sets values for the global matrix using inputs like getKSPB(), row_data, col_data, and the memory address of the local stiffness matrix's and also address of the first component &locLhs(0, 0). The ADD_VALUES flag is used to add values during the assembly. Transposing is performed for the symmetric global matrix which is then solved using the iterative linear solver KSP.

### OpDomainRhs()

The characterization of OpDomainRhs is straightforward and also aligns with OpDomainRhsVectorF in the context of SCL-1 for solving Poisson's equation with homogeneous BC. However, there is a slight difference in the body source term. The accompanying code snippet illustrates in the following:

struct OpDomainRhs : public OpFaceEle {
public:
OpDomainRhs(std::string field_name, ScalarFunc source_term_function)
sourceTermFunc(source_term_function) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const int nb_dofs = data.getIndices().size();
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
// get coordinates of the integration point
auto t_coords = getFTensor1CoordsAtGaussPts();
// get base functions
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;
double body_source =
sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
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;
// move to the coordinates of the next integration point
++t_coords;
}
// FILL VALUES OF THE GLOBAL VECTOR ENTRIES FROM THE LOCAL ONES
CHKERR VecSetValues<MoFEM::EssentialBcStorage>(
getKSPf(), data, &*locRhs.begin(), ADD_VALUES);
}
}
private:
ScalarFunc sourceTermFunc;
VectorDouble locRhs;
};

In the above snippet, we can find a different body_source term which is calculated by calling a Lambda function sourceTermFunc() with coordinate(x, y, z) i.e. t_coords(0), t_coords(1), and t_coords(2) as arguments as follows:

double body_source =
sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
for (int rr = 0; rr != nb_dofs; rr++) {
locRhs[rr] += t_base * body_source * a;
// move to the next base function
++t_base;
}

### OpBoundaryLhs()

OpBoundaryLhs is focusing on the boundary( $$\Gamma$$) instead of the domain. So we're using OpEdgeEle instead of OpFaceEle. You can realize the fact easily that for 2D problems any boundary element type should be edge(1D) and for 3D problems it will be face(2D) type element.This operator is quite similar to our first operator, OpDomainLhs. Therefore we are not describing it again. Further, the code snippet is here:

struct OpBoundaryLhs : public OpEdgeEle {
public:
OpBoundaryLhs(std::string row_field_name, std::string col_field_name)
: OpEdgeEle(row_field_name, col_field_name, OpEdgeEle::OPROWCOL) {
sYmm = true;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data) {
const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
locBoundaryLhs.resize(nb_row_dofs, nb_col_dofs, false);
locBoundaryLhs.clear();
// get (boundary) element length
const double edge = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
// get base functions on row
auto t_row_base = row_data.getFTensor0N();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
// get base functions on column
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (int cc = 0; cc != nb_col_dofs; cc++) {
locBoundaryLhs(rr, cc) += t_row_base * t_col_base * a;
// move to the next base functions on column
++t_col_base;
}
// move to the next base functions on row
++t_row_base;
}
// move to the weight of the next integration point
++t_w;
}
// FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
getKSPB(), row_data, col_data, &locBoundaryLhs(0, 0), ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
transLocBoundaryLhs.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transLocBoundaryLhs) = trans(locBoundaryLhs);
CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
getKSPB(), col_data, row_data, &transLocBoundaryLhs(0, 0),
}
}
}
private:
MatrixDouble locBoundaryLhs, transLocBoundaryLhs;
};

### OpBoundaryRhs()

If we compare the OpBoundaryRhs operator, it also closely resembles the OpDomainRhsVectorF operator in its structure. Additionally, it is also similar to OpBoundaryRhs. However, we are now using boundary_term here instead of body_source as follows:

double boundary_term =
boundaryFunc(t_coords(0), t_coords(1), t_coords(2));
for (int rr = 0; rr != nb_dofs; rr++) {
locBoundaryRhs[rr] += t_base * boundary_term * a;
// move to the next base function
++t_base;
}

## Main program (*.cpp)

This main class Poisson2DNonhomogeneous contains functions each of which is responsible for a certain task as described in SCL-1. However, in this case, we need to declare two functions as a global variable of sourceTermFunction and boundaryFunction. These functions are applicable for the operators OpDomainRhs and OpBoundaryRhs respectively. In addition to that, we have declared a vector type shared pointer 'boundaryMarker'. We will use and discuss while applying the boundary condition. The code snippet is as follows:

public:
// Declaration of the main function to run analysis
private:
// Declaration of other main functions called in runProgram()
// Function to calculate the Source term
static double sourceTermFunction(const double x, const double y,
const double z) {
return 200 * sin(x * 10.) * cos(y * 10.);
// return 1;
}
// Function to calculate the Boundary term
static double boundaryFunction(const double x, const double y,
const double z) {
return sin(x * 10.) * cos(y * 10.);
// return 0;
}
// Main interfaces
Simple *simpleInterface;
// Field name and approximation order
std::string domainField;
int oRder;
// Object to mark boundary entities for the assembling of domain elements
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
// Boundary entities marked for fieldsplit (block) solver - optional
};

Next, this code defines function in MoFEM to loads and reads your mesh file.

CHKERR simpleInterface->getOptions();
}

### setupProblem()

Since we are considering both domain and boundary, we need to add the domain and boundary field in the simpleInterface to setup the problem. Besides, we have selected the polynomial base function as AINSWORTH_BERNSTEIN_BEZIER_BASE and $$\mathbf{H^1}$$ scalar space, where we are considering single(1) degree of freedom for each shape function has

int oRder = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
}

### boundaryCondition()

Next, we are proceeding to discuss the detail for setting up boundary condition. This code involves marking degrees of freedom on boundary entities, merging markers, and obtaining entities for field splitting in the following:

auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->pushMarkDOFsOnEntities(simpleInterface->getProblemName(),
"BOUNDARY_CONDITION", domainField, 0, 1,
true);
// merge markers from all blocksets "BOUNDARY_CONDITION"
boundaryMarker = bc_mng->getMergedBlocksMarker({"BOUNDARY_CONDITION"});
// get entities on blocksets "BOUNDARY_CONDITION"
bc_mng->getMergedBlocksRange({"BOUNDARY_CONDITION"});
}

The line bc_mng->pushMarkDOFsOnEntities(...) is used to collect the degrees of freedom (DOFs) of the entities only for marked block set ("BOUNDARY_CONDITION"). Then, the block sets are accumulated as in container and assembled in boundaryMarker. In addition, the boundaryEntitiesForFieldsplit is to get the range/IDs for the blocks including vertex, edges or faces. It has been done to indicate the difference between the domain and boundary entities.

### assembleSystem()

At this stage, we are all set to focus on assembling all we did before to the pipeline. However, we find some changes that happens here in codes compared with SCL-2: Poisson's equation (non-homogeneous BC). Before going to that, first we add operators to the OpDomainLhsPipeline using the AddHOOps function considering scalar space H1. Then we need to assemble to the pipeline for the domain using OpSetBc. Then push the main OpDomainLhs operator, and finally unset the assembly using OpUnSetBc. We will discuss it in detail after going through the code snippet in the following:

auto pipeline_mng = mField.getInterface<PipelineManager>();
{ // Push operators to the Pipeline that is responsible for calculating LHS of
// domain elements
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetBc(domainField, true, boundaryMarker));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainLhs(domainField, domainField));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpUnSetBc(domainField));
}
{ // Push operators to the Pipeline that is responsible for calculating RHS of
// domain elements
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetBc(domainField, true, boundaryMarker));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpDomainRhs(domainField, sourceTermFunction));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpUnSetBc(domainField));
}
{ // Push operators to the Pipeline that is responsible for calculating LHS of
// boundary elements
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpSetBc(domainField, false, boundaryMarker));
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpBoundaryLhs(domainField, domainField));
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpUnSetBc(domainField));
}
{ // Push operators to the Pipeline that is responsible for calculating RHS of
// boundary elements
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpSetBc(domainField, false, boundaryMarker));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpBoundaryRhs(domainField, boundaryFunction));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpUnSetBc(domainField));
}
}

Here, the process involves using the OpSetBc, OpUnSetBc, and UDOs to push back to the RHS and LHS pipelines for assembling. Specifically, these operations are being carried out on the pipeline for the Domain and boundary in same way i.e. getOpDomainLhsPipeline(), OpBoundaryLhsPipeline. So we are only describing for only one set of assemble processes in the following:

The OpSetBc() object is responsible for assembling only for the boundary entities or everything except the boundary. In other words, when it is flagged 'true' it focuses exclusively on the domain entities and does not account for boundary entities tagged with the 'boundaryMarker.

{ // Push operators to the Pipeline that is responsible for calculating LHS of
// domain elements
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetBc(domainField, true, boundaryMarker));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainLhs(domainField, domainField));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpUnSetBc(domainField));
}

For example, when we assemble for the boundary, first we use OpSetBc(), then add OpBoundaryLhs and OpUnSetBc to the getOpBoundaryLhsPipeline() and getOpBoundaryLhsPipeline. To be noted that it's important to set the flag as false here.

To recapitulate, when we are assembling for the domain, OpSetBc() works for skipping(marked'true' to skip boundary) the boundary entities. On the other hand, it considers only the boundaries (marked 'false' to consider) when assembling only for boundary. And then push the UDO. Then OpUnSetBc() is used to set the index to the entities as it was initially.

### setIntegrationRules()

Afterwards, the computation of LHS and RHS is defined in the setIntegrationRules. The instructions for the functions setIntegrationRules() is same as it was in SCL-1: Poisson's equation (homogeneous BC).

auto pipeline_mng = mField.getInterface<PipelineManager>();
auto domain_rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
auto domain_rule_rhs = [](int, int, int p) -> int { return 2 * (p - 1); };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(domain_rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_rhs);
auto boundary_rule_lhs = [](int, int, int p) -> int { return 2 * p; };
auto boundary_rule_rhs = [](int, int, int p) -> int { return 2 * p; };
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_rhs);
}

### solveSystem()

To solve the problem, we remember that we made our assembling systems ready using the set/unseting process in boundaryCondition() for domain and boundaries as shown in Figure: 2 where we can consider the subscripts 'D' is for domain and 'B' for boundary.

Figure 2: Solving approach

At this stage, we are not solving the system of equations directly using KSP solver as we did in SCL-1: Poisson's equation (homogeneous BC). Instead of that, we are proceeding to deploy a preconditioner to the solver named FIELDSPLIT. It splits the fields for domain and boundary as an option - yes/no if needed. The part of the code is as follows:

// Setup fieldsplit (block) solver - optional: yes/no
if (1) {
PC pc;
CHKERR KSPGetPC(ksp_solver, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
// Set up FIELDSPLIT, only when user set -pc_type fieldsplit
// Identify the index for boundary entities, remaining will be for domain
// Then split the fields for boundary and domain for solving
if (is_pcfs == PETSC_TRUE) {
IS is_domain, is_boundary;
const MoFEM::Problem *problem_ptr;
CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
problem_ptr->getName(), ROW, domainField, 0, 1, &is_boundary,
&boundaryEntitiesForFieldsplit);
CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
CHKERR ISDestroy(&is_boundary);
}
}

The algorithm for the above code is to check whether the indices in the field were marked for domain or boundary. It set up the field-split solver using PETSc libraries for boundary entities(boundaryEntitiesForFieldsplit) and the remaining will be for the domain. During the iteration process, IS considers solving the option for domain and boundary so that the index sets are then used to configure the PCFIELDSPLIT solver to solve the appropriate blockset. Finally, the code releases the memory associated with the boundary index set and goes for KSP to solve. And to accumulate all these we followed by VecGhostUpdateBegin(), VecGhostUpdateEnd() and DMoFEMMeshToLocalVector to make things ready before getting the results.

CHKERR KSPSetUp(ksp_solver);
// 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);
}

The rest of the code including outputResults() and main function does the same as SCL-2: Poisson's equation (non-homogeneous BC)

CHKERR simpleInterface->getOptions();
}

Then, run the analysis and get the newly created output file, namely out_result.h5m. Convert it to *.vtk format. Then open it in Paraview and use the filter WarpByScalar, you will be able to see the deformation as bellow:

# Results

Figure 3: Poisson nonhomogeneous visualisation. Source term: 200sin(10x)cos(10y), boundary function: sin(10x)cos(10y)

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

#ifndef __POISSON2DNONHOMOGENEOUS_HPP__
#define __POISSON2DNONHOMOGENEOUS_HPP__
#include <stdlib.h>
typedef boost::function<double(const double, const double, const double)>
//! [OpDomainLhs]
struct OpDomainLhs : public OpFaceEle {
public:
OpDomainLhs(std::string row_field_name, std::string col_field_name)
: OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL) {
sYmm = true;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data) {
const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
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
// 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++) {
locLhs(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;
}
// FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
// Fill value to local stiffness matrix ignoring boundary DOFs
CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
getKSPB(), row_data, col_data, &locLhs(0, 0), ADD_VALUES);
// Fill values of symmetric local stiffness matrix
if (row_side != col_side || row_type != col_type) {
transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transLocLhs) = trans(locLhs);
CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
getKSPB(), col_data, row_data, &transLocLhs(0, 0), ADD_VALUES);
}
}
}
private:
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
};
//! [OpDomainLhs]
//! [OpDomainRhs]
struct OpDomainRhs : public OpFaceEle {
public:
OpDomainRhs(std::string field_name, ScalarFunc source_term_function)
sourceTermFunc(source_term_function) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const int nb_dofs = data.getIndices().size();
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
// get coordinates of the integration point
auto t_coords = getFTensor1CoordsAtGaussPts();
// get base functions
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;
double body_source =
sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
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;
// move to the coordinates of the next integration point
++t_coords;
}
// FILL VALUES OF THE GLOBAL VECTOR ENTRIES FROM THE LOCAL ONES
CHKERR VecSetValues<MoFEM::EssentialBcStorage>(
getKSPf(), data, &*locRhs.begin(), ADD_VALUES);
}
}
private:
};
//! [OpDomainRhs]
//! [OpBoundaryLhs]
struct OpBoundaryLhs : public OpEdgeEle {
public:
OpBoundaryLhs(std::string row_field_name, std::string col_field_name)
: OpEdgeEle(row_field_name, col_field_name, OpEdgeEle::OPROWCOL) {
sYmm = true;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data) {
const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
locBoundaryLhs.resize(nb_row_dofs, nb_col_dofs, false);
locBoundaryLhs.clear();
// get (boundary) element length
const double edge = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
// get base functions on row
auto t_row_base = row_data.getFTensor0N();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
// get base functions on column
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (int cc = 0; cc != nb_col_dofs; cc++) {
locBoundaryLhs(rr, cc) += t_row_base * t_col_base * a;
// move to the next base functions on column
++t_col_base;
}
// move to the next base functions on row
++t_row_base;
}
// move to the weight of the next integration point
++t_w;
}
// FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
getKSPB(), row_data, col_data, &locBoundaryLhs(0, 0), ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
transLocBoundaryLhs.resize(nb_col_dofs, nb_row_dofs, false);
CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
getKSPB(), col_data, row_data, &transLocBoundaryLhs(0, 0),
}
}
}
private:
};
//! [OpBoundaryLhs]
//! [OpBoundaryRhs]
struct OpBoundaryRhs : public OpEdgeEle {
public:
OpBoundaryRhs(std::string field_name, ScalarFunc boundary_function)
boundaryFunc(boundary_function) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const int nb_dofs = data.getIndices().size();
if (nb_dofs) {
locBoundaryRhs.resize(nb_dofs, false);
locBoundaryRhs.clear();
// get (boundary) element length
const double edge = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
// get coordinates at integration point
auto t_coords = getFTensor1CoordsAtGaussPts();
// 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 * edge;
double boundary_term =
boundaryFunc(t_coords(0), t_coords(1), t_coords(2));
for (int rr = 0; rr != nb_dofs; rr++) {
locBoundaryRhs[rr] += t_base * boundary_term * a;
// move to the next base function
++t_base;
}
// move to the weight of the next integration point
++t_w;
// move to the coordinates of the next integration point
++t_coords;
}
// FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR
CHKERR VecSetValues<MoFEM::EssentialBcStorage>(
getKSPf(), data, &*locBoundaryRhs.begin(), ADD_VALUES);
}
}
private:
};
//! [OpBoundaryRhs]
}; // namespace Poisson2DNonhomogeneousOperators
#endif //__POISSON2DNONHOMOGENEOUS_HPP__

## Implementation of the main program (*.cpp)

#include <stdlib.h>
using namespace MoFEM;
constexpr int SPACE_DIM = 2;
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 setupProblem();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode setIntegrationRules();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
// Function to calculate the Source term
static double sourceTermFunction(const double x, const double y,
const double z) {
return 200 * sin(x * 10.) * cos(y * 10.);
// return 1;
}
// Function to calculate the Boundary term
static double boundaryFunction(const double x, const double y,
const double z) {
return sin(x * 10.) * cos(y * 10.);
// return 0;
}
// Main interfaces
Simple *simpleInterface;
// Field name and approximation order
std::string domainField;
int oRder;
// Object to mark boundary entities for the assembling of domain elements
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
// Boundary entities marked for fieldsplit (block) solver - optional
Range boundaryEntitiesForFieldsplit;
};
: domainField("U"), mField(m_field) {}
}
//! [Setup problem]
int oRder = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
}
//! [Setup problem]
//! [Boundary condition]
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->pushMarkDOFsOnEntities(simpleInterface->getProblemName(),
"BOUNDARY_CONDITION", domainField, 0, 1,
true);
// merge markers from all blocksets "BOUNDARY_CONDITION"
boundaryMarker = bc_mng->getMergedBlocksMarker({"BOUNDARY_CONDITION"});
// get entities on blocksets "BOUNDARY_CONDITION"
bc_mng->getMergedBlocksRange({"BOUNDARY_CONDITION"});
}
//! [Boundary condition]
//! [Assemble system]
auto pipeline_mng = mField.getInterface<PipelineManager>();
{ // Push operators to the Pipeline that is responsible for calculating LHS of
// domain elements
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating RHS of
// domain elements
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating LHS of
// boundary elements
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating RHS of
// boundary elements
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
}
}
//! [Assemble system]
//! [Set integration rules]
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto domain_rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
auto domain_rule_rhs = [](int, int, int p) -> int { return 2 * (p - 1); };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(domain_rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_rhs);
auto boundary_rule_lhs = [](int, int, int p) -> int { return 2 * p; };
auto boundary_rule_rhs = [](int, int, int p) -> int { return 2 * p; };
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_rhs);
}
//! [Set integration rules]
//! [Solve system]
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto ksp_solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(ksp_solver);
// Create RHS and solution vectors
auto dm = simpleInterface->getDM();
auto F = createDMVector(dm);
auto D = vectorDuplicate(F);
// Setup fieldsplit (block) solver - optional: yes/no
if (1) {
PC pc;
CHKERR KSPGetPC(ksp_solver, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
// Set up FIELDSPLIT, only when user set -pc_type fieldsplit
// Identify the index for boundary entities, remaining will be for domain
// Then split the fields for boundary and domain for solving
if (is_pcfs == PETSC_TRUE) {
IS is_domain, is_boundary;
const MoFEM::Problem *problem_ptr;
CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
problem_ptr->getName(), ROW, domainField, 0, 1, &is_boundary,
// CHKERR ISView(is_boundary, PETSC_VIEWER_STDOUT_SELF);
CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
CHKERR ISDestroy(&is_boundary);
}
}
CHKERR KSPSetUp(ksp_solver);
// 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();
pipeline_mng->getBoundaryLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
auto d_ptr = boost::make_shared<VectorDouble>();
auto l_ptr = boost::make_shared<VectorDouble>();
auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{domainField, d_ptr}}, {}, {}, {}));
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile("out_result.h5m");
}
}
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
Poisson2DNonhomogeneous poisson_problem(m_field);
CHKERR poisson_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
Poisson2DNonhomogeneousOperators::ScalarFunc
boost::function< double(const double, const double, const double)> ScalarFunc
Definition: poisson_2d_nonhomogeneous.hpp:20
Poisson2DNonhomogeneousOperators::OpBoundaryRhs
[OpBoundaryLhs]
Definition: poisson_2d_nonhomogeneous.hpp:243
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPf
Vec getKSPf() const
Definition: ForcesAndSourcesCore.hpp:1088
Poisson2DNonhomogeneousOperators::OpDomainRhs::sourceTermFunc
ScalarFunc sourceTermFunc
Definition: poisson_2d_nonhomogeneous.hpp:162
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
Poisson2DNonhomogeneousOperators::i
FTensor::Index< 'i', 2 > i
Definition: poisson_2d_nonhomogeneous.hpp:17
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
Poisson2DNonhomogeneousOperators::OpBoundaryLhs
[OpDomainRhs]
Definition: poisson_2d_nonhomogeneous.hpp:168
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
Poisson2DNonhomogeneous::oRder
int oRder
Definition: poisson_2d_nonhomogeneous.cpp:48
Poisson2DNonhomogeneousOperators::OpDomainRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: poisson_2d_nonhomogeneous.hpp:109
Poisson2DNonhomogeneousOperators::OpDomainLhs::OpDomainLhs
OpDomainLhs(std::string row_field_name, std::string col_field_name)
Definition: poisson_2d_nonhomogeneous.hpp:24
Poisson2DNonhomogeneousOperators::OpDomainLhs
[OpDomainLhs]
Definition: poisson_2d_nonhomogeneous.hpp:22
Poisson2DNonhomogeneousOperators::OpDomainRhs
[OpDomainLhs]
Definition: poisson_2d_nonhomogeneous.hpp:103
Poisson2DNonhomogeneous::simpleInterface
Simple * simpleInterface
Definition: poisson_2d_nonhomogeneous.cpp:44
Poisson2DNonhomogeneous::assembleSystem
MoFEMErrorCode assembleSystem()
[Boundary condition]
Definition: poisson_2d_nonhomogeneous.cpp:108
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
Poisson2DNonhomogeneousOperators::OpDomainRhs::locRhs
VectorDouble locRhs
Definition: poisson_2d_nonhomogeneous.hpp:163
Poisson2DNonhomogeneous::mField
MoFEM::Interface & mField
Definition: poisson_2d_nonhomogeneous.cpp:43
Definition: Simple.cpp:194
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
poisson_2d_nonhomogeneous.hpp
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
BasicFiniteElements.hpp
Poisson2DNonhomogeneousOperators::OpBoundaryLhs::locBoundaryLhs
MatrixDouble locBoundaryLhs
Definition: poisson_2d_nonhomogeneous.hpp:238
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Poisson2DNonhomogeneous::outputResults
MoFEMErrorCode outputResults()
[Solve system]
Definition: poisson_2d_nonhomogeneous.cpp:231
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
operator doWork is executed on FE rows &columns
Definition: ForcesAndSourcesCore.hpp:569
MoFEM::DataOperator::doWork
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.
Definition: DataOperators.hpp:40
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight
auto getFTensor0IntegrationWeight()
Get integration weights.
Definition: ForcesAndSourcesCore.hpp:1240
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::ForcesAndSourcesCore::UserDataOperator::getMeasure
double getMeasure() const
get measure of element
Definition: ForcesAndSourcesCore.hpp:1275
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1236
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
Poisson2DNonhomogeneousOperators::OpDomainLhs::doWork
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.
Definition: poisson_2d_nonhomogeneous.hpp:29
Poisson2DNonhomogeneousOperators::OpBoundaryRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: poisson_2d_nonhomogeneous.hpp:249
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:747
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
Poisson2DNonhomogeneous::Poisson2DNonhomogeneous
Poisson2DNonhomogeneous(MoFEM::Interface &m_field)
Definition: poisson_2d_nonhomogeneous.cpp:57
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
Definition: poisson_2d_nonhomogeneous.cpp:60
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Poisson2DNonhomogeneous::domainField
std::string domainField
Definition: poisson_2d_nonhomogeneous.cpp:47
Poisson2DNonhomogeneousOperators::OpBoundaryRhs::OpBoundaryRhs
OpBoundaryRhs(std::string field_name, ScalarFunc boundary_function)
Definition: poisson_2d_nonhomogeneous.hpp:245
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
double
convert.type
type
Definition: convert.py:64
Poisson2DNonhomogeneous::boundaryFunction
static double boundaryFunction(const double x, const double y, const double z)
Definition: poisson_2d_nonhomogeneous.cpp:36
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
Poisson2DHomogeneousOperators::body_source
const double body_source
Definition: poisson_2d_homogeneous.hpp:36
Poisson2DNonhomogeneousOperators::OpBoundaryLhs::doWork
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.
Definition: poisson_2d_nonhomogeneous.hpp:175
Poisson2DNonhomogeneous::sourceTermFunction
static double sourceTermFunction(const double x, const double y, const double z)
Definition: poisson_2d_nonhomogeneous.cpp:30
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
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::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
Poisson2DNonhomogeneousOperators::OpDomainRhs::OpDomainRhs
OpDomainRhs(std::string field_name, ScalarFunc source_term_function)
Definition: poisson_2d_nonhomogeneous.hpp:105
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
Poisson2DNonhomogeneous::runProgram
MoFEMErrorCode runProgram()
Definition: poisson_2d_nonhomogeneous.cpp:260
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
Poisson2DNonhomogeneous::solveSystem
MoFEMErrorCode solveSystem()
[Set integration rules]
Definition: poisson_2d_nonhomogeneous.cpp:180
Poisson2DNonhomogeneousOperators::OpBoundaryRhs::locBoundaryRhs
VectorDouble locBoundaryRhs
Definition: poisson_2d_nonhomogeneous.hpp:303
Poisson2DNonhomogeneous::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
[Assemble system]
Definition: poisson_2d_nonhomogeneous.cpp:160
Poisson2DNonhomogeneous::boundaryEntitiesForFieldsplit
Range boundaryEntitiesForFieldsplit
Definition: poisson_2d_nonhomogeneous.cpp:54
MoFEM::ForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: ForcesAndSourcesCore.hpp:482
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
Poisson2DNonhomogeneousOperators::OpDomainLhs::locLhs
MatrixDouble locLhs
Definition: poisson_2d_nonhomogeneous.hpp:98
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', 2 >
Poisson2DNonhomogeneousOperators::OpDomainLhs::transLocLhs
MatrixDouble transLocLhs
Definition: poisson_2d_nonhomogeneous.hpp:98
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:545
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator
default operator for EDGE element
Definition: EdgeElementForcesAndSourcesCore.hpp:68
Poisson2DNonhomogeneousOperators::OpDomainLhs::boundaryMarker
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: poisson_2d_nonhomogeneous.hpp:97
Poisson2DNonhomogeneousOperators::OpBoundaryRhs::boundaryFunc
ScalarFunc boundaryFunc
Definition: poisson_2d_nonhomogeneous.hpp:302
MoFEM::OpUnSetBc
Definition: FormsIntegrators.hpp:49
Range
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::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
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:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
Poisson2DNonhomogeneous::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Setup problem]
Definition: poisson_2d_nonhomogeneous.cpp:90
Poisson2DNonhomogeneous::setupProblem
MoFEMErrorCode setupProblem()
Definition: poisson_2d_nonhomogeneous.cpp:72
MoFEM::OpSetBc
Set indices on entities on finite element.
Definition: FormsIntegrators.hpp:38
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:426
Poisson2DNonhomogeneousOperators::OpBoundaryLhs::transLocBoundaryLhs
MatrixDouble transLocBoundaryLhs
Definition: poisson_2d_nonhomogeneous.hpp:238
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEMErrorCode addBoundaryField(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 boundary.
Definition: Simple.cpp:354
Poisson2DNonhomogeneous::boundaryMarker
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: poisson_2d_nonhomogeneous.cpp:51
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:683
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:389
MoFEM::DataOperator::sYmm
bool sYmm
If true assume that matrix is symmetric structure.
Definition: DataOperators.hpp:82
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:429
MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPB
Mat getKSPB() const
Definition: ForcesAndSourcesCore.hpp:1104
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor1CoordsAtGaussPts
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
Definition: ForcesAndSourcesCore.hpp:1269
F
@ F
Definition: free_surface.cpp:394
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
Poisson2DNonhomogeneousOperators::OpBoundaryLhs::OpBoundaryLhs
OpBoundaryLhs(std::string row_field_name, std::string col_field_name)
Definition: poisson_2d_nonhomogeneous.hpp:170
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698