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