v0.14.0
Loading...
Searching...
No Matches
poisson_2d_homogeneous.hpp
Go to the documentation of this file.
1/**
2 * \file poisson_2d_homogeneous.hpp
3 * \example poisson_2d_homogeneous.hpp
4 *
5 * Solution of poisson equation. Direct implementation of User Data Operators
6 * for teaching proposes.
7 *
8 * \note In practical application we suggest use form integrators to generalise
9 * and simplify code. However, here we like to expose user to ways how to
10 * implement data operator from scratch.
11 */
12
13// Define name if it has not been defined yet
14#ifndef __POISSON_2D_HOMOGENEOUS_HPP__
15#define __POISSON_2D_HOMOGENEOUS_HPP__
16
17// Use of alias for some specific functions
18// We are solving Poisson's equation in 2D so Face element is used
19using EntData = EntitiesFieldData::EntData;
20
21template <int DIM>
22using ElementsAndOps = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>;
24using DomainEleOp = DomainEle::UserDataOperator;
25
27 FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
28
29// Namespace that contains necessary UDOs, will be included in the main program
31
32// Declare FTensor index for 2D problem
34
35// For simplicity, source term f will be constant throughout the domain
36const double body_source = 5.;
37
39public:
40 OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)
41 : AssemblyDomainEleOp(row_field_name, col_field_name,
42 DomainEleOp::OPROWCOL) {}
43
44 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
45 EntitiesFieldData::EntData &col_data) {
47
48 const int nb_row_dofs = row_data.getIndices().size();
49 const int nb_col_dofs = col_data.getIndices().size();
50
51 this->locMat.resize(nb_row_dofs, nb_col_dofs, false);
52 this->locMat.clear();
53
54 // get element area
55 const double area = getMeasure();
56
57 // get number of integration points
58 const int nb_integration_points = getGaussPts().size2();
59 // get integration weights
60 auto t_w = getFTensor0IntegrationWeight();
61
62 // get derivatives of base functions on row
63 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
64
65 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
66 for (int gg = 0; gg != nb_integration_points; gg++) {
67 const double a = t_w * area;
68
69 for (int rr = 0; rr != nb_row_dofs; ++rr) {
70 // get derivatives of base functions on column
71 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
72
73 for (int cc = 0; cc != nb_col_dofs; cc++) {
74 this->locMat(rr, cc) += t_row_diff_base(i) * t_col_diff_base(i) * a;
75
76 // move to the derivatives of the next base functions on column
77 ++t_col_diff_base;
78 }
79
80 // move to the derivatives of the next base functions on row
81 ++t_row_diff_base;
82 }
83
84 // move to the weight of the next integration point
85 ++t_w;
86 }
87
89 }
90};
91
93public:
96
97 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
99
100 const int nb_dofs = data.getIndices().size();
101
102 this->locF.resize(nb_dofs, false);
103 this->locF.clear();
104
105 // get element area
106 const double area = getMeasure();
107
108 // get number of integration points
109 const int nb_integration_points = getGaussPts().size2();
110 // get integration weights
111 auto t_w = getFTensor0IntegrationWeight();
112
113 // get base function
114 auto t_base = data.getFTensor0N();
115
116 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
117 for (int gg = 0; gg != nb_integration_points; gg++) {
118 const double a = t_w * area;
119
120 for (int rr = 0; rr != nb_dofs; rr++) {
121 this->locF[rr] += t_base * body_source * a;
122
123 // move to the next base function
124 ++t_base;
125 }
126
127 // move to the weight of the next integration point
128 ++t_w;
129 }
130
132 }
133
134};
135
136}; // namespace Poisson2DHomogeneousOperators
137
138#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)
VectorDouble locF
local entity vector
MatrixDouble locMat
local entity block matrix
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpDomainLhsMatrixK(std::string row_field_name, std::string col_field_name)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)