v0.14.0
Loading...
Searching...
No Matches
rd_stdOperators.hpp
Go to the documentation of this file.
1#ifndef __STDOPERATORS_HPP__
2#define __STDOPERATORS_HPP__
3
4#include <stdlib.h>
6
7namespace StdRDOperators {
8
10using OpEle = FaceElementForcesAndSourcesCore::UserDataOperator;
11
14
15using EntData = DataForcesAndSourcesCore::EntData;
16
17// const double alpha = 0.01;
18// const double gma = 0.002;
19// const double b = 0.15;
20// const double c = 8.00;
21// const double mu1 = 0.20;
22// const double mu2 = 0.30;
23
24const double B = 0.0;
26const double natural_bc_values = 0.0;
27const double essential_bc_values = 0.0;
28const int order = 2;
29// const int dim = 3;
31
32struct BlockData {
34 MatrixDouble coef;
35 VectorDouble rate;
36
38
39 double B0; // species mobility
40
41 BlockData() : B0(2e-3) {
42 coef.resize(3, 3, false);
43 rate.resize(3, false);
44 coef.clear();
45 rate.clear();
46 coef(0, 0) = 1.0;
47 coef(0, 1) = 0.0;
48 coef(0, 2) = 0.0;
49 coef(1, 0) = 0.0;
50 coef(1, 1) = 1.0;
51 coef(1, 2) = 0.0;
52 coef(2, 0) = 0.0;
53 coef(2, 1) = 0.0;
54 coef(2, 2) = 1.0;
55
56 for (int i = 0; i < 3; ++i) {
57 rate[i] = 1.0;
58 }
59 }
60};
61
63 MatrixDouble grads; ///< Gradients of field "u" at integration points
64 VectorDouble values; ///< Values of field "u" at integration points
65 VectorDouble
66 dot_values; ///< Rate of values of field "u" at integration points
67 VectorDouble slow_values;
68
69 MatrixDouble invJac; ///< Inverse of element jacobian
70
72};
73
75 OpAssembleMass(std::string fieldu, SmartPetscObj<Mat> m)
76 : OpEle(fieldu, fieldu, OpEle::OPROWCOL), M(m)
77
78 {
79 sYmm = true;
80 }
81 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
82 EntityType col_type, EntData &row_data,
83 EntData &col_data) {
85 const int nb_row_dofs = row_data.getIndices().size();
86 const int nb_col_dofs = col_data.getIndices().size();
87 if (nb_row_dofs && nb_col_dofs) {
88 const int nb_integration_pts = getGaussPts().size2();
89 mat.resize(nb_row_dofs, nb_col_dofs, false);
90 mat.clear();
91 auto t_row_base = row_data.getFTensor0N();
93 const double vol = getMeasure();
94 for (int gg = 0; gg != nb_integration_pts; ++gg) {
95 const double a = t_w * vol;
96 for (int rr = 0; rr != nb_row_dofs; ++rr) {
97 auto t_col_base = col_data.getFTensor0N(gg, 0);
98 for (int cc = 0; cc != nb_col_dofs; ++cc) {
99 mat(rr, cc) += a * t_row_base * t_col_base;
100 ++t_col_base;
101 }
102 ++t_row_base;
103 }
104 ++t_w;
105 }
106 CHKERR MatSetValues(M, row_data, col_data, &mat(0, 0), ADD_VALUES);
107 if (row_side != col_side || row_type != col_type) {
108 transMat.resize(nb_col_dofs, nb_row_dofs, false);
109 noalias(transMat) = trans(mat);
110 CHKERR MatSetValues(M, col_data, row_data, &transMat(0, 0), ADD_VALUES);
111 }
112 }
114 }
115
116private:
117 MatrixDouble mat, transMat;
119};
120
121struct OpComputeSlowValue : public OpEle {
122 OpComputeSlowValue(std::string mass_field,
123 boost::shared_ptr<PreviousData> &data1,
124 boost::shared_ptr<PreviousData> &data2,
125 boost::shared_ptr<PreviousData> &data3,
126 std::map<int, BlockData> &block_map)
127 : OpEle(mass_field, OpEle::OPROW), commonData1(data1), commonData2(data2),
128 commonData3(data3), massField(mass_field), setOfBlock(block_map) {}
129 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
131 boost::shared_ptr<VectorDouble> slow_value_ptr1(commonData1,
132 &commonData1->slow_values);
133 boost::shared_ptr<VectorDouble> slow_value_ptr2(commonData2,
134 &commonData2->slow_values);
135 boost::shared_ptr<VectorDouble> slow_value_ptr3(commonData3,
136 &commonData3->slow_values);
137
138 VectorDouble &vec1 = *slow_value_ptr1;
139 VectorDouble &vec2 = *slow_value_ptr2;
140 VectorDouble &vec3 = *slow_value_ptr3;
141 const int nb_integration_pts = getGaussPts().size2();
142 if (type == MBVERTEX) {
143 vec1.resize(nb_integration_pts, false);
144 vec2.resize(nb_integration_pts, false);
145 vec3.resize(nb_integration_pts, false);
146 vec1.clear();
147 vec2.clear();
148 vec3.clear();
149 }
150 const int nb_dofs = data.getIndices().size();
151
152 if (nb_dofs) {
153 auto find_block_data = [&]() {
155 BlockData *block_raw_ptr = nullptr;
156 for (auto &m : setOfBlock) {
157 if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
158 block_raw_ptr = &m.second;
159 break;
160 }
161 }
162 return block_raw_ptr;
163 };
164
165 auto block_data_ptr = find_block_data();
166 if (!block_data_ptr)
167 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
168
169 auto &block_data = *block_data_ptr;
170
171 const int nb_integration_pts = getGaussPts().size2();
172
173 auto t_slow_values1 = getFTensor0FromVec(vec1);
174 auto t_slow_values2 = getFTensor0FromVec(vec2);
175 auto t_slow_values3 = getFTensor0FromVec(vec3);
176
177 auto t_mass_values1 = getFTensor0FromVec(commonData1->values);
178 auto t_mass_values2 = getFTensor0FromVec(commonData2->values);
179 auto t_mass_values3 = getFTensor0FromVec(commonData3->values);
180
181 // cout << "r1 : " << block_data.rate[0] << endl;
182
183 for (int gg = 0; gg != nb_integration_pts; ++gg) {
184 t_slow_values1 = block_data.rate[0] * t_mass_values1 *
185 (1.0 - block_data.coef(0, 0) * t_mass_values1 -
186 block_data.coef(0, 1) * t_mass_values2 -
187 block_data.coef(0, 2) * t_mass_values3);
188 t_slow_values2 = block_data.rate[1] * t_mass_values2 *
189 (1.0 - block_data.coef(1, 0) * t_mass_values1 -
190 block_data.coef(1, 1) * t_mass_values2 -
191 block_data.coef(1, 2) * t_mass_values3);
192
193 t_slow_values3 = block_data.rate[2] * t_mass_values3 *
194 (1.0 - block_data.coef(2, 0) * t_mass_values1 -
195 block_data.coef(2, 1) * t_mass_values2 -
196 block_data.coef(2, 2) * t_mass_values3);
197 ++t_slow_values1;
198 ++t_slow_values2;
199 ++t_slow_values3;
200
201 ++t_mass_values1;
202 ++t_mass_values2;
203 ++t_mass_values3;
204 }
205 }
207 }
208
209private:
210 std::string massField;
211 boost::shared_ptr<PreviousData> commonData1;
212 boost::shared_ptr<PreviousData> commonData2;
213 boost::shared_ptr<PreviousData> commonData3;
214 std::map<int, BlockData> setOfBlock;
215};
216
218 typedef boost::function<double(const double, const double, const double)>
220 typedef boost::function<FTensor::Tensor1<double, 3>(
221 const double, const double, const double)>
223 OpAssembleSlowRhs(std::string field,
224 boost::shared_ptr<PreviousData> &data,
225 FVal exact_value,
226 FVal exact_dot,
227 FVal exact_lap)
228 : OpEle(field, OpEle::OPROW), commonData(data), exactValue(exact_value),
229 exactDot(exact_dot), exactLap(exact_lap) {}
230
231 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
233
234 const int nb_dofs = data.getIndices().size();
235 if (nb_dofs) {
236 vecF.resize(nb_dofs, false);
237 vecF.clear();
238 const int nb_integration_pts = getGaussPts().size2();
239
240 auto t_slow_value = getFTensor0FromVec(commonData->slow_values);
241
242 auto t_base = data.getFTensor0N();
243 auto t_w = getFTensor0IntegrationWeight();
244 const double vol = getMeasure();
245
246 // cout << "measure : " << getMeasure() << endl;
247 const double ct = getFEMethod()->ts_t;
248 auto t_coords = getFTensor1CoordsAtGaussPts();
249 for (int gg = 0; gg != nb_integration_pts; ++gg) {
250 const double a = vol * t_w;
251
252 double u_dot = exactDot(t_coords(NX), t_coords(NY), ct);
253 double u_lap = exactLap(t_coords(NX), t_coords(NY), ct);
254
255 // double f = u_dot - u_lap;
256 const double f = t_slow_value;
257
258 for (int rr = 0; rr != nb_dofs; ++rr) {
259 const double b = a * f * t_base;
260 vecF[rr] += b;
261 ++t_base;
262 }
263 ++t_slow_value;
264 ++t_w;
265 ++t_coords;
266 }
267 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
268 PETSC_TRUE);
269 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
270 ADD_VALUES);
271 }
273 }
274
275private:
276 boost::shared_ptr<PreviousData> commonData;
277 VectorDouble vecF;
278
282
286};
287
288
289template <int DIM> struct OpAssembleStiffRhs : OpEle {
290 typedef boost::function<double(const double, const double, const double)>
292 typedef boost::function<FTensor::Tensor1<double, 3>(
293 const double, const double, const double)>
295
296 OpAssembleStiffRhs(std::string field, boost::shared_ptr<PreviousData> &data,
297 std::map<int, BlockData> &block_map, FVal exact_value,
298 FVal exact_dot, FVal exact_lap)
299 : OpEle(field, OpEle::OPROW), commonData(data), setOfBlock(block_map),
300 exactValue(exact_value), exactDot(exact_dot), exactLap(exact_lap) {}
301
302 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
304 const int nb_dofs = data.getIndices().size();
305 if (nb_dofs) {
306 auto find_block_data = [&]() {
308 BlockData *block_raw_ptr = nullptr;
309 for (auto &m : setOfBlock) {
310 if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
311 block_raw_ptr = &m.second;
312 break;
313 }
314 }
315 return block_raw_ptr;
316 };
317
318 auto block_data_ptr = find_block_data();
319 if (!block_data_ptr)
320 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
321
322 auto &block_data = *block_data_ptr;
323
324 vecF.resize(nb_dofs, false);
325 vecF.clear();
326 const int nb_integration_pts = getGaussPts().size2();
327
328 auto t_dot_val = getFTensor0FromVec(commonData->dot_values);
329 auto t_val = getFTensor0FromVec(commonData->values);
330 auto t_grad = getFTensor1FromMat<DIM>(commonData->grads);
331
332 auto t_base = data.getFTensor0N();
333 auto t_diff_base = data.getFTensor1DiffN<DIM>();
334 auto t_w = getFTensor0IntegrationWeight();
335
336 FTensor::Index<'i', DIM> i;
337
338 // cout << "B0 : " << block_data.B0 << endl;
339
340 const double vol = getMeasure();
341
342 const double ct = getFEMethod()->ts_t;
343 auto t_coords = getFTensor1CoordsAtGaussPts();
344
345 for (int gg = 0; gg != nb_integration_pts; ++gg) {
346 const double a = vol * t_w;
347
348 double u_dot = exactDot(t_coords(NX), t_coords(NY), ct);
349 // double u_lap = exactLap(t_coords(NX), t_coords(NY), ct);
350
351 // double f = u_dot - u_lap;
352 double u_lap =
353 -block_data.B0 * exactLap(t_coords(NX), t_coords(NY), ct);
354 // cout << "B0 : " << block_data.B0 << endl;
355 // double f = u_dot + u_lap;
356 double f = 0.0;
357
358 for (int rr = 0; rr != nb_dofs; ++rr) {
359 vecF[rr] +=
360 a * (t_base * t_dot_val +
361 (block_data.B0 + B * t_val) * t_diff_base(i) * t_grad(i) -
362 t_base * f);
363 ++t_diff_base;
364 ++t_base;
365 }
366 ++t_dot_val;
367 ++t_grad;
368 ++t_val;
369 ++t_w;
370 ++t_coords;
371 }
372 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
373 PETSC_TRUE);
374 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
375 ADD_VALUES);
376 }
378 }
379
380private:
381 boost::shared_ptr<PreviousData> commonData;
382 std::map<int, BlockData> setOfBlock;
383 VectorDouble vecF;
384 std::string field;
385
389
393};
394
395template <int DIM> struct OpAssembleStiffLhs : OpEle {
396
397 OpAssembleStiffLhs(std::string fieldu, boost::shared_ptr<PreviousData> &data,
398 std::map<int, BlockData> &block_map)
399 : OpEle(fieldu, fieldu, OpEle::OPROWCOL), commonData(data),
400 setOfBlock(block_map) {
401 sYmm = true;
402 }
403 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
404 EntityType col_type, EntData &row_data,
405 EntData &col_data) {
407
408 const int nb_row_dofs = row_data.getIndices().size();
409 const int nb_col_dofs = col_data.getIndices().size();
410 // cerr << "In doWork() : (row, col) = (" << nb_row_dofs << ", " <<
411 // nb_col_dofs << ")" << endl;
412 if (nb_row_dofs && nb_col_dofs) {
413 auto find_block_data = [&]() {
415 BlockData *block_raw_ptr = nullptr;
416 for (auto &m : setOfBlock) {
417 if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
418 block_raw_ptr = &m.second;
419 break;
420 }
421 }
422 return block_raw_ptr;
423 };
424
425 auto block_data_ptr = find_block_data();
426 if (!block_data_ptr)
427 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
428
429 auto &block_data = *block_data_ptr;
430
431 mat.resize(nb_row_dofs, nb_col_dofs, false);
432 mat.clear();
433 const int nb_integration_pts = getGaussPts().size2();
434 auto t_row_base = row_data.getFTensor0N();
435
436 auto t_val = getFTensor0FromVec(commonData->values);
437 auto t_grad = getFTensor1FromMat<DIM>(commonData->grads);
438
439 auto t_row_diff_base = row_data.getFTensor1DiffN<DIM>();
440 auto t_w = getFTensor0IntegrationWeight();
441 const double ts_a = getFEMethod()->ts_a;
442 const double vol = getMeasure();
443
444 FTensor::Index<'i', DIM> i;
445
446 // cout << "B0 : " << block_data.B0 << endl;
447
448 for (int gg = 0; gg != nb_integration_pts; ++gg) {
449 const double a = vol * t_w;
450
451 for (int rr = 0; rr != nb_row_dofs; ++rr) {
452 auto t_col_base = col_data.getFTensor0N(gg, 0);
453 auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
454 for (int cc = 0; cc != nb_col_dofs; ++cc) {
455
456 mat(rr, cc) +=
457 a * (t_row_base * t_col_base * ts_a +
458 (block_data.B0 + B * t_val) * t_row_diff_base(i) *
459 t_col_diff_base(i) +
460 B * t_col_base * t_grad(i) * t_row_diff_base(i));
461
462 ++t_col_base;
463 ++t_col_diff_base;
464 }
465 ++t_row_base;
466 ++t_row_diff_base;
467 }
468 ++t_w;
469 ++t_val;
470 ++t_grad;
471 }
472 CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
473 ADD_VALUES);
474 if (row_side != col_side || row_type != col_type) {
475 transMat.resize(nb_col_dofs, nb_row_dofs, false);
476 noalias(transMat) = trans(mat);
477 CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
478 &transMat(0, 0), ADD_VALUES);
479 }
480 }
482 }
483
484 private :
485 boost::shared_ptr<PreviousData> commonData;
486 std::map<int, BlockData> setOfBlock;
487 MatrixDouble mat, transMat;
488};
489
491{
493 : OpBoundaryEle(mass_field, OpBoundaryEle::OPROW),
495 cerr << "OpAssembleNaturalBCRhsTau()" << endl;
496 }
497
498 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
500 const int nb_dofs = data.getIndices().size();
501
502 if (nb_dofs) {
503 EntityHandle row_side_ent = getFEEntityHandle();
504 bool is_natural =
505 (natural_bd_ents.find(row_side_ent) != natural_bd_ents.end());
506 if (is_natural) {
507 // cerr << "In NaturalBCRhsTau..." << endl;
508 vecF.resize(nb_dofs, false);
509 vecF.clear();
510 const int nb_integration_pts = getGaussPts().size2();
511 auto t_row_base = data.getFTensor0N();
512
513 auto dir = getDirection();
514 FTensor::Tensor1<double, 3> t_normal(-dir[1], dir[0], dir[2]);
515 FTensor::Index<'i', 3> i;
516 auto t_w = getFTensor0IntegrationWeight();
517 const double pi = 3.141592653589793;
518 const double ct = getFEMethod()->ts_t;
519 for (int gg = 0; gg != nb_integration_pts; ++gg) {
520 const double a = t_w;
521 double x = getCoordsAtGaussPts()(gg, 0);
522 double y = getCoordsAtGaussPts()(gg, 1);
523
524 double mm = -10 * 8 * pi * cos(2 * pi * x) * sin(2 * pi * y) *
525 sin(2 * pi * ct);
526 double nn = -10 * 8 * pi * sin(2 * pi * x) * cos(2 * pi * y) *
527 sin(2 * pi * ct);
528
529 FTensor::Tensor1<double, 3> t_bd_val(mm, nn, 0.0);
530 double h = t_bd_val(i) * t_normal(i);
531 for (int rr = 0; rr != nb_dofs; ++rr) {
532 vecF[rr] += t_row_base * h * a;
533 ++t_row_base;
534 }
535 ++t_w;
536 }
537 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
538 PETSC_TRUE);
539 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
540 ADD_VALUES);
541 }
542 }
544 }
545
546private:
547 VectorDouble vecF;
549};
550
551struct OpError : public OpEle {
552 typedef boost::function<double(const double, const double, const double)>
554 typedef boost::function<FTensor::Tensor1<double, 3>(
555 const double, const double, const double)>
557 double &eRror;
558 OpError(FVal exact_value, FVal exact_lap, FGrad exact_grad,
559 boost::shared_ptr<PreviousData> &prev_data,
560 std::map<int, BlockData> &block_map, double &err)
561 : OpEle("ERROR", OpEle::OPROW), exactVal(exact_value),
562 exactLap(exact_lap), exactGrad(exact_grad), prevData(prev_data),
563 setOfBlock(block_map), eRror(err) {}
564 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
566 const int nb_dofs = data.getFieldData().size();
567 // cout << "nb_error_dofs : " << nb_dofs << endl;
568 if (nb_dofs) {
569
570 auto find_block_data = [&]() {
572 BlockData *block_raw_ptr = nullptr;
573 for (auto &m : setOfBlock) {
574 if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
575 block_raw_ptr = &m.second;
576 break;
577 }
578 }
579 return block_raw_ptr;
580 };
581
582 auto block_data_ptr = find_block_data();
583 if (!block_data_ptr)
584 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
585
586 auto &block_data = *block_data_ptr;
587
588 auto t_value = getFTensor0FromVec(prevData->values);
589 auto t_grad = getFTensor1FromMat<2>(prevData->grads);
590 // cout << "t_grad : " << t_grad << endl;
591 // auto t_grad = getFTensor1FromMat<3>(prevData->grads);
592 data.getFieldData().clear();
593 const double vol = getMeasure();
594 const int nb_integration_pts = getGaussPts().size2();
595 auto t_w = getFTensor0IntegrationWeight();
596 double dt;
597 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
598 double ct = getFEMethod()->ts_t - dt;
599 auto t_coords = getFTensor1CoordsAtGaussPts();
600
601 FTensor::Index<'j', 2> j;
602 // cout << "B0 : " << block_data.B0 << endl;
603
604 for (int gg = 0; gg != nb_integration_pts; ++gg) {
605 const double a = vol * t_w;
606
607 double mass_exact = exactVal(t_coords(NX), t_coords(NY), ct);
608 double flux_lap =
609 -block_data.B0 * exactLap(t_coords(NX), t_coords(NY), ct);
610 auto flux_exact = exactGrad(t_coords(NX), t_coords(NY), ct);
611
612 // cout << "grad_exact : " << flux_exact << endl;
613 // cout << "grad_value : " << t_grad << endl;
614 // cout << "--------------------" << endl;
615
616 double flux_error =
617 pow(block_data.B0, 2) * (pow(flux_exact(0) - t_grad(0), 2) +
618 pow(flux_exact(1) - t_grad(1), 2));
619
620 double local_error = pow(mass_exact - t_value, 2); // + flux_error;
621
622 data.getFieldData()[0] += a * local_error;
623 eRror += a * local_error;
624
625 ++t_w;
626 ++t_value;
627 ++t_grad;
628 ++t_coords;
629 }
630
631 data.getFieldDofs()[0]->getFieldData() = data.getFieldData()[0];
632 }
634 }
635
636private:
640 boost::shared_ptr<PreviousData> prevData;
641 std::map<int, BlockData> setOfBlock;
642
645};
646
647struct Monitor : public FEMethod {
648 double &eRror;
649 Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj<DM> &dm,
650 boost::shared_ptr<PostProcFaceOnRefinedMesh> &post_proc, double &err)
651 : cOmm(comm), rAnk(rank), dM(dm), postProc(post_proc), eRror(err){};
656 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every_nth_step",
657 &save_every_nth_step, PETSC_NULL);
658 if (ts_step % save_every_nth_step == 0) {
660 CHKERR postProc->writeFile(
661 "out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
662 }
663 Vec error_per_proc;
664 CHKERR VecCreateMPI(cOmm, 1, PETSC_DECIDE, &error_per_proc);
665 auto get_global_error = [&]() {
667 CHKERR VecSetValue(error_per_proc, rAnk, eRror, INSERT_VALUES);
669 };
670 CHKERR get_global_error();
671 CHKERR VecAssemblyBegin(error_per_proc);
672 CHKERR VecAssemblyEnd(error_per_proc);
673 double error_sum;
674 CHKERR VecSum(error_per_proc, &error_sum);
675 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Error : %3.4e \n", error_sum);
676 eRror = 0;
678 }
679
680private:
682 boost::shared_ptr<PostProcFaceOnRefinedMesh> postProc;
683 MPI_Comm cOmm;
684 const int rAnk;
685};
686
687}; // namespace StdRDOperators
688
689#endif
constexpr double a
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
double dt
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
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.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'i', 3 > i
const double essential_bc_values
const double natural_bc_values
double h
bool sYmm
If true assume that matrix is symmetric structure.
structure for User Loop Methods on finite elements
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
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
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
intrusive_ptr for managing petsc objects
PetscReal ts_t
time
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
PetscInt ts_step
time step number
MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProc
MoFEMErrorCode operator()()
function is run for every finite element
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcFaceOnRefinedMesh > &post_proc, double &err)
MoFEMErrorCode postProcess()
function is run at the end of loop
OpAssembleMass(std::string fieldu, SmartPetscObj< Mat > m)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleNaturalBCRhs(std::string mass_field, Range &natural_bd_ents)
boost::shared_ptr< PreviousData > commonData
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> FGrad
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleSlowRhs(std::string field, boost::shared_ptr< PreviousData > &data, FVal exact_value, FVal exact_dot, FVal exact_lap)
boost::function< double(const double, const double, const double)> FVal
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleStiffLhs(std::string fieldu, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
boost::shared_ptr< PreviousData > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> FGrad
OpAssembleStiffRhs(std::string field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map, FVal exact_value, FVal exact_dot, FVal exact_lap)
boost::function< double(const double, const double, const double)> FVal
std::map< int, BlockData > setOfBlock
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpComputeSlowValue(std::string mass_field, boost::shared_ptr< PreviousData > &data1, boost::shared_ptr< PreviousData > &data2, boost::shared_ptr< PreviousData > &data3, std::map< int, BlockData > &block_map)
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData2
boost::shared_ptr< PreviousData > commonData3
boost::shared_ptr< PreviousData > commonData1
FTensor::Number< 1 > NY
boost::shared_ptr< PreviousData > prevData
std::map< int, BlockData > setOfBlock
OpError(FVal exact_value, FVal exact_lap, FGrad exact_grad, boost::shared_ptr< PreviousData > &prev_data, std::map< int, BlockData > &block_map, double &err)
boost::function< double(const double, const double, const double)> FVal
FTensor::Number< 0 > NX
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> FGrad
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
VectorDouble dot_values
Rate of values of field "u" at integration points.
MatrixDouble invJac
Inverse of element jacobian.
VectorDouble values
Values of field "u" at integration points.
MatrixDouble grads
Gradients of field "u" at integration points.