v0.14.0
Loading...
Searching...
No Matches
poisson_2d_nonhomogeneous.hpp
Go to the documentation of this file.
1#ifndef __POISSON2DNONHOMOGENEOUS_HPP__
2#define __POISSON2DNONHOMOGENEOUS_HPP__
3
4#include <stdlib.h>
6
9
12
14
16
18
19typedef boost::function<double(const double, const double, const double)>
21//! [OpDomainLhs]
22struct OpDomainLhs : public OpFaceEle {
23public:
24 OpDomainLhs(std::string row_field_name, std::string col_field_name)
25 : OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL) {
26 sYmm = true;
27 }
28
29 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
30 EntityType col_type, EntData &row_data,
31 EntData &col_data) {
33
34 const int nb_row_dofs = row_data.getIndices().size();
35 const int nb_col_dofs = col_data.getIndices().size();
36
37 if (nb_row_dofs && nb_col_dofs) {
38
39 locLhs.resize(nb_row_dofs, nb_col_dofs, false);
40 locLhs.clear();
41
42 // get element area
43 const double area = getMeasure();
44
45 // get number of integration points
46 const int nb_integration_points = getGaussPts().size2();
47 // get integration weights
49
50 // get derivatives of base functions on row
51 auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
52
53 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
54 for (int gg = 0; gg != nb_integration_points; gg++) {
55
56 const double a = t_w * area;
57
58 for (int rr = 0; rr != nb_row_dofs; ++rr) {
59 // get derivatives of base functions on column
60 auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
61
62 for (int cc = 0; cc != nb_col_dofs; cc++) {
63
64 locLhs(rr, cc) += t_row_diff_base(i) * t_col_diff_base(i) * a;
65
66 // move to the derivatives of the next base functions on column
67 ++t_col_diff_base;
68 }
69
70 // move to the derivatives of the next base functions on row
71 ++t_row_diff_base;
72 }
73
74 // move to the weight of the next integration point
75 ++t_w;
76 }
77
78 // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
79
80 // Fill value to local stiffness matrix ignoring boundary DOFs
81 CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
82 getKSPB(), row_data, col_data, &locLhs(0, 0), ADD_VALUES);
83
84 // Fill values of symmetric local stiffness matrix
85 if (row_side != col_side || row_type != col_type) {
86 transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
87 noalias(transLocLhs) = trans(locLhs);
88 CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
89 getKSPB(), col_data, row_data, &transLocLhs(0, 0), ADD_VALUES);
90 }
91 }
92
94 }
95
96private:
97 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
98 MatrixDouble locLhs, transLocLhs;
99};
100//! [OpDomainLhs]
101
102//! [OpDomainRhs]
103struct OpDomainRhs : public OpFaceEle {
104public:
105 OpDomainRhs(std::string field_name, ScalarFunc source_term_function)
107 sourceTermFunc(source_term_function) {}
108
109 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
111
112 const int nb_dofs = data.getIndices().size();
113
114 if (nb_dofs) {
115
116 locRhs.resize(nb_dofs, false);
117 locRhs.clear();
118
119 // get element area
120 const double area = getMeasure();
121
122 // get number of integration points
123 const int nb_integration_points = getGaussPts().size2();
124 // get integration weights
125 auto t_w = getFTensor0IntegrationWeight();
126 // get coordinates of the integration point
127 auto t_coords = getFTensor1CoordsAtGaussPts();
128
129 // get base functions
130 auto t_base = data.getFTensor0N();
131
132 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
133 for (int gg = 0; gg != nb_integration_points; gg++) {
134
135 const double a = t_w * area;
136 double body_source =
137 sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
138
139 for (int rr = 0; rr != nb_dofs; rr++) {
140
141 locRhs[rr] += t_base * body_source * a;
142
143 // move to the next base function
144 ++t_base;
145 }
146
147 // move to the weight of the next integration point
148 ++t_w;
149 // move to the coordinates of the next integration point
150 ++t_coords;
151 }
152
153 // FILL VALUES OF THE GLOBAL VECTOR ENTRIES FROM THE LOCAL ONES
154 CHKERR VecSetValues<MoFEM::EssentialBcStorage>(
155 getKSPf(), data, &*locRhs.begin(), ADD_VALUES);
156 }
157
159 }
160
161private:
163 VectorDouble locRhs;
164};
165//! [OpDomainRhs]
166
167//! [OpBoundaryLhs]
168struct OpBoundaryLhs : public OpEdgeEle {
169public:
170 OpBoundaryLhs(std::string row_field_name, std::string col_field_name)
171 : OpEdgeEle(row_field_name, col_field_name, OpEdgeEle::OPROWCOL) {
172 sYmm = true;
173 }
174
175 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
176 EntityType col_type, EntData &row_data,
177 EntData &col_data) {
179
180 const int nb_row_dofs = row_data.getIndices().size();
181 const int nb_col_dofs = col_data.getIndices().size();
182
183 if (nb_row_dofs && nb_col_dofs) {
184
185 locBoundaryLhs.resize(nb_row_dofs, nb_col_dofs, false);
186 locBoundaryLhs.clear();
187
188 // get (boundary) element length
189 const double edge = getMeasure();
190
191 // get number of integration points
192 const int nb_integration_points = getGaussPts().size2();
193 // get integration weights
194 auto t_w = getFTensor0IntegrationWeight();
195
196 // get base functions on row
197 auto t_row_base = row_data.getFTensor0N();
198
199 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
200 for (int gg = 0; gg != nb_integration_points; gg++) {
201 const double a = t_w * edge;
202
203 for (int rr = 0; rr != nb_row_dofs; ++rr) {
204 // get base functions on column
205 auto t_col_base = col_data.getFTensor0N(gg, 0);
206
207 for (int cc = 0; cc != nb_col_dofs; cc++) {
208
209 locBoundaryLhs(rr, cc) += t_row_base * t_col_base * a;
210
211 // move to the next base functions on column
212 ++t_col_base;
213 }
214 // move to the next base functions on row
215 ++t_row_base;
216 }
217
218 // move to the weight of the next integration point
219 ++t_w;
220 }
221
222 // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
223 CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
224 getKSPB(), row_data, col_data, &locBoundaryLhs(0, 0), ADD_VALUES);
225 if (row_side != col_side || row_type != col_type) {
226 transLocBoundaryLhs.resize(nb_col_dofs, nb_row_dofs, false);
227 noalias(transLocBoundaryLhs) = trans(locBoundaryLhs);
228 CHKERR MatSetValues<MoFEM::EssentialBcStorage>(
229 getKSPB(), col_data, row_data, &transLocBoundaryLhs(0, 0),
230 ADD_VALUES);
231 }
232 }
233
235 }
236
237private:
239};
240//! [OpBoundaryLhs]
241
242//! [OpBoundaryRhs]
243struct OpBoundaryRhs : public OpEdgeEle {
244public:
245 OpBoundaryRhs(std::string field_name, ScalarFunc boundary_function)
247 boundaryFunc(boundary_function) {}
248
249 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
251
252 const int nb_dofs = data.getIndices().size();
253
254 if (nb_dofs) {
255
256 locBoundaryRhs.resize(nb_dofs, false);
257 locBoundaryRhs.clear();
258
259 // get (boundary) element length
260 const double edge = getMeasure();
261
262 // get number of integration points
263 const int nb_integration_points = getGaussPts().size2();
264 // get integration weights
265 auto t_w = getFTensor0IntegrationWeight();
266 // get coordinates at integration point
267 auto t_coords = getFTensor1CoordsAtGaussPts();
268
269 // get base function
270 auto t_base = data.getFTensor0N();
271
272 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
273 for (int gg = 0; gg != nb_integration_points; gg++) {
274
275 const double a = t_w * edge;
276 double boundary_term =
277 boundaryFunc(t_coords(0), t_coords(1), t_coords(2));
278
279 for (int rr = 0; rr != nb_dofs; rr++) {
280
281 locBoundaryRhs[rr] += t_base * boundary_term * a;
282
283 // move to the next base function
284 ++t_base;
285 }
286
287 // move to the weight of the next integration point
288 ++t_w;
289 // move to the coordinates of the next integration point
290 ++t_coords;
291 }
292
293 // FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR
294 CHKERR VecSetValues<MoFEM::EssentialBcStorage>(
295 getKSPf(), data, &*locBoundaryRhs.begin(), ADD_VALUES);
296 }
297
299 }
300
301private:
303 VectorDouble locBoundaryRhs;
304};
305//! [OpBoundaryRhs]
306}; // namespace Poisson2DNonhomogeneousOperators
307
308#endif //__POISSON2DNONHOMOGENEOUS_HPP__
constexpr double a
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
boost::function< double(const double, const double, const double)> ScalarFunc
constexpr auto field_name
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 getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ 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.
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)
OpDomainLhs(std::string row_field_name, std::string col_field_name)
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.