v0.14.0
poisson_2d_homogeneous.hpp

Solution of poisson equation. Direct implementation of User Data Operators for teaching proposes.

Note
In practical application we suggest use form integrators to generalise and simplify code. However, here we like to expose user to ways how to implement data operator from scratch.
/**
* \file poisson_2d_homogeneous.hpp
* \example poisson_2d_homogeneous.hpp
*
* Solution of poisson equation. Direct implementation of User Data Operators
* for teaching proposes.
*
* \note In practical application we suggest use form integrators to generalise
* and simplify code. However, here we like to expose user to ways how to
* implement data operator from scratch.
*/
// Define name if it has not been defined yet
#ifndef __POISSON_2D_HOMOGENEOUS_HPP__
#define __POISSON_2D_HOMOGENEOUS_HPP__
// Use of alias for some specific functions
// We are solving Poisson's equation in 2D so Face element is used
template <int DIM>
using ElementsAndOps = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>;
FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
// Namespace that contains necessary UDOs, will be included in the main program
// Declare FTensor index for 2D problem
// For simplicity, source term f will be constant throughout the domain
const double body_source = 5.;
struct OpDomainLhsMatrixK : public AssemblyDomainEleOp {
public:
OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)
: AssemblyDomainEleOp(row_field_name, col_field_name,
DomainEleOp::OPROWCOL) {}
const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
this->locMat.resize(nb_row_dofs, nb_col_dofs, false);
this->locMat.clear();
// get element area
const double area = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on row
auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
// get derivatives of base functions on column
auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
for (int cc = 0; cc != nb_col_dofs; cc++) {
this->locMat(rr, cc) += t_row_diff_base(i) * t_col_diff_base(i) * a;
// move to the derivatives of the next base functions on column
++t_col_diff_base;
}
// move to the derivatives of the next base functions on row
++t_row_diff_base;
}
// move to the weight of the next integration point
++t_w;
}
}
};
struct OpDomainRhsVectorF : public AssemblyDomainEleOp {
public:
: AssemblyDomainEleOp(field_name, field_name, DomainEleOp::OPROW) {}
const int nb_dofs = data.getIndices().size();
this->locF.resize(nb_dofs, false);
this->locF.clear();
// get element area
const double area = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get base function
auto t_base = data.getFTensor0N();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_dofs; rr++) {
this->locF[rr] += t_base * body_source * a;
// move to the next base function
++t_base;
}
// move to the weight of the next integration point
++t_w;
}
}
};
}; // namespace Poisson2DHomogeneousOperators
#endif //__POISSON_2D_HOMOGENEOUS_HPP__
Poisson2DHomogeneousOperators::OpDomainRhsVectorF::OpDomainRhsVectorF
OpDomainRhsVectorF(std::string field_name)
Definition: poisson_2d_homogeneous.hpp:94
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition: FormsIntegrators.hpp:247
Poisson2DHomogeneousOperators::OpDomainLhsMatrixK::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: poisson_2d_homogeneous.hpp:44
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
Poisson2DHomogeneousOperators::i
FTensor::Index< 'i', SPACE_DIM > i
Definition: poisson_2d_homogeneous.hpp:33
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
MoFEM::OpBaseImpl::locF
VectorDouble locF
local entity vector
Definition: FormsIntegrators.hpp:249
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
a
constexpr double a
Definition: approx_sphere.cpp:30
Poisson2DHomogeneousOperators::body_source
const double body_source
Definition: poisson_2d_homogeneous.hpp:36
Poisson2DHomogeneousOperators::OpDomainLhsMatrixK::OpDomainLhsMatrixK
OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)
Definition: poisson_2d_homogeneous.hpp:40
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', SPACE_DIM >
EntData
EntitiesFieldData::EntData EntData
Definition: poisson_2d_homogeneous.hpp:19
ElementsAndOps
Definition: child_and_parent.cpp:18
DomainEleOp
Poisson2DHomogeneousOperators
Definition: poisson_2d_homogeneous.hpp:30
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
AssemblyDomainEleOp
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::OpBase AssemblyDomainEleOp
Definition: poisson_2d_homogeneous.hpp:27
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
Poisson2DHomogeneousOperators::OpDomainRhsVectorF::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Definition: poisson_2d_homogeneous.hpp:97
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346