v0.14.0
Loading...
Searching...
No Matches
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:

\[ \begin{equation} \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} \end{equation} \]

and

\[ \begin{equation} \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} \end{equation} \]

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
using EntData = EntitiesFieldData::EntData;
// 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 .. ..
boost::function< double(const double, const double, const double)> ScalarFunc
Data on single entity (This is passed as argument to DataOperator::doWork)

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;
};
constexpr double a
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'i', SPACE_DIM > i
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
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.
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.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element

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;
};
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
constexpr auto field_name
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.

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),
ADD_VALUES);
}
}
}
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
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();
// 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
};
Deprecated interface functions.
static double boundaryFunction(const double x, const double y, const double z)
MoFEMErrorCode solveSystem()
[Set integration rules]
MoFEMErrorCode setIntegrationRules()
[Assemble system]
static double sourceTermFunction(const double x, const double y, const double z)
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode boundaryCondition()
[Setup problem]
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode outputResults()
[Solve system]
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode assembleSystem()
[Boundary condition]

readMesh()

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

CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

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

CHKERR simpleInterface->addBoundaryField(domainField, H1,
int oRder = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
}
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
@ H1
continuous field
Definition: definitions.h:85

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
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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);
}
static Index< 'p', 3 > p

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);
}
}
@ ROW
Definition: definitions.h:123
keeps basic data about problem

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);
}
@ F
double D

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

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)

Further requisite: https://en.wikiversity.org/wiki/Partial_differential_equations/Poisson_Equation.

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>
using EntData = EntitiesFieldData::EntData;
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),
ADD_VALUES);
}
}
}
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__
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
const double body_source
bool sYmm
If true assume that matrix is symmetric structure.
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
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.
OpBoundaryLhs(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.
OpBoundaryRhs(std::string field_name, ScalarFunc boundary_function)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
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.
OpDomainRhs(std::string field_name, ScalarFunc source_term_function)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.

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
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
// 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
};
: domainField("U"), mField(m_field) {}
//! [Read mesh]
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
//! [Read mesh]
//! [Setup problem]
CHKERR simpleInterface->addBoundaryField(domainField, H1,
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";
// 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;
}
std::string param_file
static char help[]
int main()
Definition: adol-c_atom.cpp:46
constexpr int SPACE_DIM
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
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:509
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:412
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
Poisson2DNonhomogeneous(MoFEM::Interface &m_field)