v0.13.2
Loading...
Searching...
No Matches
nonlinear_poisson_2d.hpp
Go to the documentation of this file.
1#ifndef __NONLINEARPOISSON2D_HPP__
2#define __NONLINEARPOISSON2D_HPP__
3
4#include <stdlib.h>
6
9
12
14
16
18
19typedef boost::function<double(const double, const double, const double)>
21
23 // This struct is for data required for the update of Newton's iteration
24
27};
28
30public:
32 std::string row_field_name, std::string col_field_name,
33 boost::shared_ptr<DataAtGaussPoints> &common_data,
34 boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
35 : OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL),
36 commonData(common_data), boundaryMarker(boundary_marker) {
37 sYmm = false;
38 }
39
40 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
41 EntityType col_type, EntData &row_data,
42 EntData &col_data) {
44
45 const int nb_row_dofs = row_data.getIndices().size();
46 const int nb_col_dofs = col_data.getIndices().size();
47
48 if (nb_row_dofs && nb_col_dofs) {
49
50 locLhs.resize(nb_row_dofs, nb_col_dofs, false);
51 locLhs.clear();
52
53 // get element area
54 const double area = getMeasure();
55
56 // get number of integration points
57 const int nb_integration_points = getGaussPts().size2();
58 // get integration weights
60 // get solution (field value) at integration points
61 auto t_field = getFTensor0FromVec(commonData->fieldValue);
62 // get gradient of field at integration points
63 auto t_field_grad = getFTensor1FromMat<2>(commonData->fieldGrad);
64
65 // get base functions on row
66 auto t_row_base = row_data.getFTensor0N();
67 // get derivatives of base functions on row
68 auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
69
70 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
71 for (int gg = 0; gg != nb_integration_points; gg++) {
72 const double a = t_w * area;
73
74 for (int rr = 0; rr != nb_row_dofs; ++rr) {
75 // get base functions on column
76 auto t_col_base = col_data.getFTensor0N(gg, 0);
77 // get derivatives of base functions on column
78 auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
79
80 for (int cc = 0; cc != nb_col_dofs; cc++) {
81 locLhs(rr, cc) += (((1 + t_field * t_field) * t_row_diff_base(i) *
82 t_col_diff_base(i)) +
83 (2.0 * t_field * t_field_grad(i) *
84 t_row_diff_base(i) * t_col_base)) *
85 a;
86
87 // move to the next base functions on column
88 ++t_col_base;
89 // move to the derivatives of the next base function on column
90 ++t_col_diff_base;
91 }
92
93 // move to the next base function on row
94 ++t_row_base;
95 // move to the derivatives of the next base function on row
96 ++t_row_diff_base;
97 }
98
99 // move to the weight of the next integration point
100 ++t_w;
101 // move to the solution (field value) at the next integration point
102 ++t_field;
103 // move to the gradient of field value at the next integration point
104 ++t_field_grad;
105 }
106
107 // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
108
109 // store original row indices
110 auto row_indices = row_data.getIndices();
111 // mark the boundary DOFs (as -1) and fill only domain DOFs
112 if (boundaryMarker) {
113 for (int r = 0; r != row_data.getIndices().size(); ++r) {
114 if ((*boundaryMarker)[row_data.getLocalIndices()[r]]) {
115 row_data.getIndices()[r] = -1;
116 }
117 }
118 }
119 // fill value to local stiffness matrix ignoring boundary DOFs
120 CHKERR MatSetValues(getSNESB(), row_data, col_data, &locLhs(0, 0),
121 ADD_VALUES);
122 // revert back row indices to the original
123 row_data.getIndices().swap(row_indices);
124 }
125
127 }
128
129private:
130 boost::shared_ptr<DataAtGaussPoints> commonData;
131 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
133};
134
136public:
138 std::string field_name, ScalarFunc source_term_function,
139 boost::shared_ptr<DataAtGaussPoints> &common_data,
140 boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
142 sourceTermFunc(source_term_function), commonData(common_data),
143 boundaryMarker(boundary_marker) {}
144
145 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
147
148 const int nb_dofs = data.getIndices().size();
149
150 if (nb_dofs) {
151 locRhs.resize(nb_dofs, false);
152 locRhs.clear();
153
154 // get element area
155 const double area = getMeasure();
156
157 // get number of integration points
158 const int nb_integration_points = getGaussPts().size2();
159 // get integration weights
160 auto t_w = getFTensor0IntegrationWeight();
161 // get coordinates of the integration point
162 auto t_coords = getFTensor1CoordsAtGaussPts();
163 // get solution (field value) at integration point
164 auto t_field = getFTensor0FromVec(commonData->fieldValue);
165 // get gradient of field value of integration point
166 auto t_field_grad = getFTensor1FromMat<2>(commonData->fieldGrad);
167
168 // get base function
169 auto t_base = data.getFTensor0N();
170 // get derivatives of base function
171 auto t_grad_base = data.getFTensor1DiffN<2>();
172
173 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
174 for (int gg = 0; gg != nb_integration_points; gg++) {
175 const double a = t_w * area;
176 double body_source =
177 sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
178
179 // calculate the local vector
180 for (int rr = 0; rr != nb_dofs; rr++) {
181 locRhs[rr] -=
182 (t_base * body_source -
183 t_grad_base(i) * t_field_grad(i) * (1 + t_field * t_field)) *
184 a;
185
186 // move to the next base function
187 ++t_base;
188 // move to the derivatives of the next base function
189 ++t_grad_base;
190 }
191
192 // move to the weight of the next integration point
193 ++t_w;
194 // move to the coordinates of the next integration point
195 ++t_coords;
196 // move to the solution (field value) at the next integration point
197 ++t_field;
198 // move to the gradient of field value at the next integration point
199 ++t_field_grad;
200 }
201
202 // FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR
203
204 // store original row indices
205 auto row_indices = data.getIndices();
206 // mark the boundary DOFs (as -1) and fill only domain DOFs
207 if (boundaryMarker) {
208 for (int r = 0; r != data.getIndices().size(); ++r) {
209 if ((*boundaryMarker)[data.getLocalIndices()[r]]) {
210 data.getIndices()[r] = -1;
211 }
212 }
213 }
214 // fill value to local vector ignoring boundary DOFs
215 CHKERR VecSetOption(getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
216 CHKERR VecSetValues(getSNESf(), data, &*locRhs.begin(), ADD_VALUES);
217 // revert back the indices
218 data.getIndices().swap(row_indices);
219 }
220
222 }
223
224private:
226 boost::shared_ptr<DataAtGaussPoints> commonData;
227 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
229};
230
232public:
233 OpBoundaryTangentMatrix(std::string row_field_name,
234 std::string col_field_name)
235 : OpEdgeEle(row_field_name, col_field_name, OpEdgeEle::OPROWCOL) {
236 sYmm = true;
237 }
238
239 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
240 EntityType col_type, EntData &row_data,
241 EntData &col_data) {
243
244 const int nb_row_dofs = row_data.getIndices().size();
245 const int nb_col_dofs = col_data.getIndices().size();
246 // cerr << nb_row_dofs;
247
248 if (nb_row_dofs && nb_col_dofs) {
249 locBoundaryLhs.resize(nb_row_dofs, nb_col_dofs, false);
250 locBoundaryLhs.clear();
251
252 // get (boundary) element length
253 const double edge = getMeasure();
254
255 // get number of integration points
256 const int nb_integration_points = getGaussPts().size2();
257 // get integration weights
258 auto t_w = getFTensor0IntegrationWeight();
259
260 // get base function on row
261 auto t_row_base = row_data.getFTensor0N();
262
263 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
264 for (int gg = 0; gg != nb_integration_points; gg++) {
265 const double a = t_w * edge;
266
267 for (int rr = 0; rr != nb_row_dofs; ++rr) {
268 // get base function on column
269 auto t_col_base = col_data.getFTensor0N(gg, 0);
270
271 for (int cc = 0; cc != nb_col_dofs; cc++) {
272 locBoundaryLhs(rr, cc) += t_row_base * t_col_base * a;
273
274 // move to the next base function on column
275 ++t_col_base;
276 }
277
278 // move to the next base function on row
279 ++t_row_base;
280 }
281
282 // move to the weight of the next integration point
283 ++t_w;
284 }
285
286 // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
287 CHKERR MatSetValues(getSNESB(), row_data, col_data, &locBoundaryLhs(0, 0),
288 ADD_VALUES);
289
290 if (row_side != col_side || row_type != col_type) {
291 transLocBoundaryLhs.resize(nb_col_dofs, nb_row_dofs, false);
292 noalias(transLocBoundaryLhs) = trans(locBoundaryLhs);
293 CHKERR MatSetValues(getSNESB(), col_data, row_data,
294 &transLocBoundaryLhs(0, 0), ADD_VALUES);
295 }
296 // cerr << locBoundaryLhs << endl;
297 // cerr << transLocBoundaryLhs << endl;
298 }
299
301 }
302
303private:
305};
306
308public:
309 OpBoundaryResidualVector(std::string field_name, ScalarFunc boundary_function,
310 boost::shared_ptr<DataAtGaussPoints> &common_data)
312 boundaryFunc(boundary_function), commonData(common_data) {}
313
314 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
316
317 const int nb_dofs = data.getIndices().size();
318
319 if (nb_dofs) {
320
321 locBoundaryRhs.resize(nb_dofs, false);
322 locBoundaryRhs.clear();
323
324 // get (boundary) element length
325 const double edge = getMeasure();
326
327 // get number of integration points
328 const int nb_integration_points = getGaussPts().size2();
329 // get integration weights
330 auto t_w = getFTensor0IntegrationWeight();
331 // get coordinates at integration point
332 auto t_coords = getFTensor1CoordsAtGaussPts();
333 // get solution (field value) at integration point
334 auto t_field = getFTensor0FromVec(commonData->fieldValue);
335
336 // get base function
337 auto t_base = data.getFTensor0N();
338
339 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
340 for (int gg = 0; gg != nb_integration_points; gg++) {
341 const double a = t_w * edge;
342 double boundary_term =
343 boundaryFunc(t_coords(0), t_coords(1), t_coords(2));
344
345 // calculate the local vector
346 for (int rr = 0; rr != nb_dofs; rr++) {
347 locBoundaryRhs[rr] -= t_base * (boundary_term - t_field) * a;
348
349 // move to the next base function
350 ++t_base;
351 }
352
353 // move to the weight of the next integration point
354 ++t_w;
355 // move to the coordinates of the next integration point
356 ++t_coords;
357 // move to the solution (field value) at the next integration point
358 ++t_field;
359 }
360
361 // FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR
362 CHKERR VecSetValues(getSNESf(), data, &*locBoundaryRhs.begin(),
363 ADD_VALUES);
364 }
365
367 }
368
369private:
371 boost::shared_ptr<DataAtGaussPoints> commonData;
373};
374
375}; // namespace NonlinearPoissonOps
376
377#endif //__NONLINEARPOISSON2D_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
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
FTensor::Index< 'i', 2 > i
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 & getLocalIndices() const
get local indices of dofs on entity
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
OpBoundaryResidualVector(std::string field_name, ScalarFunc boundary_function, boost::shared_ptr< DataAtGaussPoints > &common_data)
boost::shared_ptr< DataAtGaussPoints > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpBoundaryTangentMatrix(std::string row_field_name, std::string col_field_name)
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.
boost::shared_ptr< DataAtGaussPoints > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
OpDomainResidualVector(std::string field_name, ScalarFunc source_term_function, boost::shared_ptr< DataAtGaussPoints > &common_data, boost::shared_ptr< std::vector< unsigned char > > boundary_marker=nullptr)
boost::shared_ptr< DataAtGaussPoints > commonData
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
OpDomainTangentMatrix(std::string row_field_name, std::string col_field_name, boost::shared_ptr< DataAtGaussPoints > &common_data, boost::shared_ptr< std::vector< unsigned char > > boundary_marker=nullptr)
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.