v0.14.0
Loading...
Searching...
No Matches
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
using EntData = EntitiesFieldData::EntData;
template <int DIM>
using ElementsAndOps = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>;
using DomainEleOp = DomainEle::UserDataOperator;
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,
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
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:
OpDomainRhsVectorF(std::string field_name)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
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__
constexpr double a
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#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
FTensor::Index< 'i', SPACE_DIM > i
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
Data on single entity (This is passed as argument to DataOperator::doWork)
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
VectorDouble locF
local entity vector
MatrixDouble locMat
local entity block matrix
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase AssemblyDomainEleOp