v0.13.2
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;
// 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 OpFaceEle {
public:
OpDomainLhsMatrixK(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
CHKERR MatSetValues(getKSPB(), row_data, col_data, &locLhs(0, 0),
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transLocLhs) = trans(locLhs);
CHKERR MatSetValues(getKSPB(), col_data, row_data, &transLocLhs(0, 0),
ADD_VALUES);
}
}
}
private:
MatrixDouble locLhs, transLocLhs;
};
struct OpDomainRhsVectorF : public OpFaceEle {
public:
OpDomainRhsVectorF(std::string field_name)
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 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++) {
locRhs[rr] += t_base * body_source * a;
// move to the next base function
++t_base;
}
// move to the weight of the next integration point
++t_w;
}
// FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR
// Ignoring DOFs on boundary (index -1)
CHKERR VecSetOption(getKSPf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
CHKERR VecSetValues(getKSPf(), data, &locRhs(0), ADD_VALUES);
}
}
private:
VectorDouble locRhs;
};
}; // namespace Poisson2DHomogeneousOperators
#endif //__POISSON_2D_HOMOGENEOUS_HPP__
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
constexpr auto field_name
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator OpFaceEle
bool sYmm
If true assume that matrix is symmetric structure.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.