v0.14.0
Loading...
Searching...
No Matches
RDOperators.hpp
Go to the documentation of this file.
1#ifndef __RDOPERATORS_HPP__
2#define __RDOPERATORS_HPP__
3
4#include <stdlib.h>
6
8
9using FaceEle = MoFEM::FaceElementForcesAndSourcesCoreSwitch<0>;
10
11using BoundaryEle = MoFEM::EdgeElementForcesAndSourcesCoreSwitch<0>;
12
15
16using EntData = DataForcesAndSourcesCore::EntData;
17
18
19
20const double B = 0.0;
21const double B_epsilon = 0.0;
22
24// const int order = 3; ///< approximation order
25const double init_value = 1.0;
26const double essen_value = 0;
27const double natu_value = 0;
28// const int dim = 3;
30
34
37
39
42
44 jac.resize(2, 2, false);
45 inv_jac.resize(2, 2, false);
46 }
47};
48
49struct BlockData {
51 double a11, a12, a13, a21, a22, a23, a31, a32, a33;
52
53 double r1, r2, r3;
54
56
57 double B0; // species mobility
58
60 : a11(1), a12(0), a13(0),
61 a21(0), a22(1), a23(0),
62 a31(0), a32(0), a33(1),
63 B0(2e-3), r1(1), r2(1), r3(1) {}
64};
65
66double compute_init_val(const double x, const double y, const double z) {
67 return 0.0;
68}
69
70double compute_essen_bc(const double x, const double y, const double z) {
71 return 0.0;
72 }
73
74 double compute_natu_bc(const double x, const double y, const double z){
75 return 0.0;
76 }
77
78
79
81 OpComputeSlowValue(std::string mass_field,
82 boost::shared_ptr<PreviousData> &data1,
83 boost::shared_ptr<PreviousData> &data2,
84 boost::shared_ptr<PreviousData> &data3,
85 std::map<int, BlockData> &block_map)
86 : OpFaceEle(mass_field, OpFaceEle::OPROW), commonData1(data1),
87 commonData2(data2), commonData3(data3), massField(mass_field),
88 setOfBlock(block_map) {}
89 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
91 boost::shared_ptr<VectorDouble> slow_value_ptr1(
92 commonData1, &commonData1->slow_values);
93 boost::shared_ptr<VectorDouble> slow_value_ptr2(
94 commonData2, &commonData2->slow_values);
95 boost::shared_ptr<VectorDouble> slow_value_ptr3(
96 commonData3, &commonData3->slow_values);
97
98 VectorDouble &vec1 = *slow_value_ptr1;
99 VectorDouble &vec2 = *slow_value_ptr2;
100 VectorDouble &vec3 = *slow_value_ptr3;
101 const int nb_integration_pts = getGaussPts().size2();
102 if (type == MBVERTEX) {
103 vec1.resize(nb_integration_pts, false);
104 vec2.resize(nb_integration_pts, false);
105 vec3.resize(nb_integration_pts, false);
106 vec1.clear();
107 vec2.clear();
108 vec3.clear();
109 }
110 const int nb_dofs = data.getIndices().size();
111
112 if (nb_dofs) {
113 auto find_block_data = [&]() {
115 BlockData *block_raw_ptr = nullptr;
116 for (auto &m : setOfBlock) {
117 if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
118 block_raw_ptr = &m.second;
119 break;
120 }
121 }
122 return block_raw_ptr;
123 };
124
125 auto block_data_ptr = find_block_data();
126 if (!block_data_ptr)
127 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
128
129 auto &block_data = *block_data_ptr;
130
131 const int nb_integration_pts = getGaussPts().size2();
132
133 auto t_slow_values1 = getFTensor0FromVec(vec1);
134 auto t_slow_values2 = getFTensor0FromVec(vec2);
135 auto t_slow_values3 = getFTensor0FromVec(vec3);
136
137 auto t_mass_values1 = getFTensor0FromVec(commonData1->mass_values);
138 auto t_mass_values2 = getFTensor0FromVec(commonData2->mass_values);
139 auto t_mass_values3 = getFTensor0FromVec(commonData3->mass_values);
140 // cout << "r1 : " << block_data.r1 << endl;
141 for (int gg = 0; gg != nb_integration_pts; ++gg) {
142 t_slow_values1 = block_data.r1 * t_mass_values1 *
143 (1.0 - block_data.a11 * t_mass_values1 -
144 block_data.a12 * t_mass_values2 -
145 block_data.a13 * t_mass_values3);
146 t_slow_values2 = block_data.r2 * t_mass_values2 *
147 (1.0 - block_data.a21 * t_mass_values1 -
148 block_data.a22 * t_mass_values2 -
149 block_data.a23 * t_mass_values3);
150
151 t_slow_values3 = block_data.r3 * t_mass_values3 *
152 (1.0 - block_data.a31 * t_mass_values1 -
153 block_data.a32 * t_mass_values2 -
154 block_data.a33 * t_mass_values3);
155 ++t_slow_values1;
156 ++t_slow_values2;
157 ++t_slow_values3;
158
159 ++t_mass_values1;
160 ++t_mass_values2;
161 ++t_mass_values3;
162 }
163 }
165 }
166
167 private:
168 std::string massField;
169 boost::shared_ptr<PreviousData> commonData1;
170 boost::shared_ptr<PreviousData> commonData2;
171 boost::shared_ptr<PreviousData> commonData3;
172 std::map<int, BlockData> setOfBlock;
173};
174
176 OpEssentialBC(const std::string &flux_field, Range &essential_bd_ents)
177 : OpBoundaryEle(flux_field, OpBoundaryEle::OPROW),
179
180 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
182 int nb_dofs = data.getIndices().size();
183 if (nb_dofs) {
185 bool is_essential =
186 (essential_bd_ents.find(fe_ent) != essential_bd_ents.end());
187 if (is_essential) {
188 int nb_gauss_pts = getGaussPts().size2();
189 int size2 = data.getN().size2();
190 if (3 * nb_dofs != static_cast<int>(data.getN().size2()))
191 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
192 "wrong number of dofs");
193 nN.resize(nb_dofs, nb_dofs, false);
194 nF.resize(nb_dofs, false);
195 nN.clear();
196 nF.clear();
197
198 auto t_row_tau = data.getFTensor1N<3>();
199
200 auto dir = getDirection();
201 double len = sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]);
202
203 FTensor::Tensor1<double, 3> t_normal(-dir[1] / len, dir[0] / len,
204 dir[2] / len);
205
206 auto t_w = getFTensor0IntegrationWeight();
207 const double vol = getMeasure();
208
209 for (int gg = 0; gg < nb_gauss_pts; gg++) {
210 const double a = t_w * vol;
211 for (int rr = 0; rr != nb_dofs; rr++) {
212 auto t_col_tau = data.getFTensor1N<3>(gg, 0);
213 nF[rr] += a * essen_value * t_row_tau(i) * t_normal(i);
214 for (int cc = 0; cc != nb_dofs; cc++) {
215 nN(rr, cc) += a * (t_row_tau(i) * t_normal(i)) *
216 (t_col_tau(i) * t_normal(i));
217 ++t_col_tau;
218 }
219 ++t_row_tau;
220 }
221 ++t_w;
222 }
223
225 cholesky_solve(nN, nF, ublas::lower());
226
227 for (auto &dof : data.getFieldDofs()) {
228 dof->getFieldData() = nF[dof->getEntDofIdx()];
229 }
230 }
231 }
233 }
234
235private:
239};
240
242 typedef boost::function<double(const double)> FVal;
243 typedef boost::function<double(const double, const double, const double)> ExactFunVal;
244 OpSkeletonSource(const std::string &mass_field,
245 FVal skeleton_fun,
246 ExactFunVal smooth_fun,
247 Range &internal_edge_ents)
248 : OpBoundaryEle(mass_field, OpBoundaryEle::OPROW),
249 internalEdges(internal_edge_ents),
250 smoothFun(smooth_fun),
251 skeletonFun(skeleton_fun) {}
252
253
254
255 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
257 int nb_dofs = data.getIndices().size();
258 if (nb_dofs) {
260 bool is_intEdge =
261 (internalEdges.find(fe_ent) != internalEdges.end());
262 if (is_intEdge) {
263 int nb_gauss_pts = getGaussPts().size2();
264 int size2 = data.getN().size2();
265 if (3 * nb_dofs != static_cast<int>(data.getN().size2()))
266 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
267 "wrong number of dofs");
268 vecF.resize(nb_dofs, false);
269
270 vecF.clear();
271
272 auto t_row_v_base = data.getFTensor0N();
273
274 auto t_w = getFTensor0IntegrationWeight();
275 const double vol = getMeasure();
276 double dt;
277 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
278 double ct = getFEMethod()->ts_t - dt;
279 auto t_coords = getFTensor1CoordsAtGaussPts();
280
281 for (int gg = 0; gg < nb_gauss_pts; gg++) {
282 const double a = t_w * vol;
283 double trace_val = -2.0 * skeletonFun(t_coords(NX)) *
284 smoothFun(t_coords(NX), t_coords(NY), ct);
285 for (int rr = 0; rr != nb_dofs; ++rr) {
286
287 vecF[rr] += a * trace_val * t_row_v_base;
288 }
289 ++t_row_v_base;
290 }
291 ++t_w;
292 }
293 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
294 PETSC_TRUE);
295 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
296 ADD_VALUES);
297 }
299 }
300
301private:
306
310};
311
312// Assembly of system mass matrix
313// //***********************************************
314
315// Mass matrix corresponding to the flux equation.
316// 01. Note that it is an identity matrix
317
318struct OpInitialMass : public OpFaceEle {
319 OpInitialMass(const std::string &mass_field, Range &inner_surface)
320 : OpFaceEle(mass_field, OpFaceEle::OPROW), innerSurface(inner_surface) {
321 }
325 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
327 int nb_dofs = data.getFieldData().size();
328 if (nb_dofs) {
330 bool is_inner_side = (innerSurface.find(fe_ent) != innerSurface.end());
331 if (is_inner_side) {
332 int nb_gauss_pts = getGaussPts().size2();
333 if (nb_dofs != static_cast<int>(data.getN().size2()))
334 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
335 "wrong number of dofs");
336 nN.resize(nb_dofs, nb_dofs, false);
337 nF.resize(nb_dofs, false);
338 nN.clear();
339 nF.clear();
340
341 auto t_row_mass = data.getFTensor0N();
342 auto t_w = getFTensor0IntegrationWeight();
343 const double vol = getMeasure();
344
345 for (int gg = 0; gg < nb_gauss_pts; gg++) {
346 const double a = t_w * vol;
347 // double r = ((double) rand() / (RAND_MAX));
348 for (int rr = 0; rr != nb_dofs; rr++) {
349 auto t_col_mass = data.getFTensor0N(gg, 0);
350 nF[rr] += a * init_value * t_row_mass;
351 for (int cc = 0; cc != nb_dofs; cc++) {
352 nN(rr, cc) += a * t_row_mass * t_col_mass;
353 ++t_col_mass;
354 }
355 ++t_row_mass;
356 }
357 ++t_w;
358 }
359
361 cholesky_solve(nN, nF, ublas::lower());
362
363 for (auto &dof : data.getFieldDofs()) {
364 dof->getFieldData() = nF[dof->getEntDofIdx()];
365
366 // this is only to check
367 // data.getFieldData()[dof->getEntDofIdx()] = nF[dof->getEntDofIdx()];
368 }
369 }
370 }
372 }
373};
374
375// Assembly of RHS for explicit (slow)
376// part//**************************************
377
378// 2. RHS for explicit part of the mass balance equation
380{
381 typedef boost::function<double(const double, const double, const double)>
383 typedef boost::function<FTensor::Tensor1<double, 3>(
384 const double, const double, const double)>
386 OpAssembleSlowRhsV(std::string mass_field,
387 boost::shared_ptr<PreviousData> &common_data,
388 FVal exact_value,
389 FVal exact_dot,
390 FVal exact_lap
391 )
392 : OpFaceEle(mass_field, OpFaceEle::OPROW)
393 , commonData(common_data)
394 , exactValue(exact_value)
395 , exactDot(exact_dot)
396 , exactLap(exact_lap)
397 {}
398
399 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
401 // cerr << "In OpAssembleSlowRhsV...." << endl;
402 const int nb_dofs = data.getIndices().size();
403 if (nb_dofs) {
404 // cerr << "In SlowRhsV..." << endl;
405 if (nb_dofs != static_cast<int>(data.getN().size2()))
406 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
407 "wrong number of dofs");
408 vecF.resize(nb_dofs, false);
409 mat.resize(nb_dofs, nb_dofs, false);
410 vecF.clear();
411 mat.clear();
412 const int nb_integration_pts = getGaussPts().size2();
413 auto t_mass_value = getFTensor0FromVec(commonData->mass_values);
414 auto t_slow_value = getFTensor0FromVec(commonData->slow_values);
415 auto t_row_v_base = data.getFTensor0N();
416 auto t_w = getFTensor0IntegrationWeight();
417 const double vol = getMeasure();
418
419 const double ct = getFEMethod()->ts_t - 0.01;
420 auto t_coords = getFTensor1CoordsAtGaussPts();
421 for (int gg = 0; gg != nb_integration_pts; ++gg) {
422 const double a = vol * t_w;
423
424 double u_dot = exactDot(t_coords(NX), t_coords(NY), ct);
425 double u_lap = exactLap(t_coords(NX), t_coords(NY), ct);
426
427 // double f = u_dot - u_lap;
428 // double f = 0;
429 for (int rr = 0; rr != nb_dofs; ++rr) {
430 auto t_col_v_base = data.getFTensor0N(gg, 0);
431 vecF[rr] += a * t_slow_value * t_row_v_base;
432 // vecF[rr] += a * f * t_row_v_base;
433 for (int cc = 0; cc != nb_dofs; ++cc) {
434 mat(rr, cc) += a * t_row_v_base * t_col_v_base;
435 ++t_col_v_base;
436 }
437 ++t_row_v_base;
438 }
439 ++t_mass_value;
440 ++t_slow_value;
441 ++t_w;
442 ++t_coords;
443 }
445 cholesky_solve(mat, vecF, ublas::lower());
446
447 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
448 PETSC_TRUE);
449 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
450 ADD_VALUES);
451 }
453 }
454
455private:
456 boost::shared_ptr<PreviousData> commonData;
462
466};
467
468// 5. RHS contribution of the natural boundary condition
470{
472 : OpBoundaryEle(flux_field, OpBoundaryEle::OPROW),
474
475 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
477 const int nb_dofs = data.getIndices().size();
478
479 if (nb_dofs) {
480 EntityHandle row_side_ent = data.getFieldDofs()[0]->getEnt();
481
482 bool is_natural =
483 (natural_bd_ents.find(row_side_ent) != natural_bd_ents.end());
484 if (is_natural) {
485 // cerr << "In NaturalBCRhsTau..." << endl;
486 vecF.resize(nb_dofs, false);
487 vecF.clear();
488 const int nb_integration_pts = getGaussPts().size2();
489 auto t_tau_base = data.getFTensor1N<3>();
490
491 auto dir = getDirection();
492 FTensor::Tensor1<double, 3> t_normal(-dir[1], dir[0], dir[2]);
493
494 auto t_w = getFTensor0IntegrationWeight();
495
496 for (int gg = 0; gg != nb_integration_pts; ++gg) {
497 const double a = t_w;
498 for (int rr = 0; rr != nb_dofs; ++rr) {
499 vecF[rr] += (t_tau_base(i) * t_normal(i) * natu_value) * a;
500 ++t_tau_base;
501 }
502 ++t_w;
503 }
504 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
505 PETSC_TRUE);
506 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
507 ADD_VALUES);
508 }
509 }
511 }
512
513private:
516};
517
518// Assembly of RHS for the implicit (stiff) part excluding the essential
519// boundary //**********************************
520// 3. Assembly of F_tau excluding the essential boundary condition
521template <int dim>
523{
524 OpAssembleStiffRhsTau(std::string flux_field,
525 boost::shared_ptr<PreviousData> &data,
526 std::map<int, BlockData> &block_map)
527 : OpFaceEle(flux_field, OpFaceEle::OPROW), commonData(data),
528 setOfBlock(block_map) {}
529
530 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
532
533 const int nb_dofs = data.getIndices().size();
534 if (nb_dofs) {
535 auto find_block_data = [&]() {
537 BlockData *block_raw_ptr = nullptr;
538 for (auto &m : setOfBlock) {
539 if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
540 block_raw_ptr = &m.second;
541 break;
542 }
543 }
544 return block_raw_ptr;
545 };
546
547 auto block_data_ptr = find_block_data();
548 if (!block_data_ptr)
549 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
550 auto &block_data = *block_data_ptr;
551
552 vecF.resize(nb_dofs, false);
553 vecF.clear();
554
555 const int nb_integration_pts = getGaussPts().size2();
556 auto t_flux_value = getFTensor1FromMat<3>(commonData->flux_values);
557 auto t_mass_value = getFTensor0FromVec(commonData->mass_values);
558 auto t_tau_base = data.getFTensor1N<3>();
559
560 auto t_tau_grad = data.getFTensor2DiffN<3, 2>();
561
562 auto t_w = getFTensor0IntegrationWeight();
563 const double vol = getMeasure();
564
565 for (int gg = 0; gg < nb_integration_pts; ++gg) {
566
567 const double K = B_epsilon + (block_data.B0 + B * t_mass_value);
568 const double K_inv = 1. / K;
569 const double a = vol * t_w;
570 for (int rr = 0; rr < nb_dofs; ++rr) {
571 double div_base = t_tau_grad(0, 0) + t_tau_grad(1, 1);
572 vecF[rr] += (K_inv * t_tau_base(i) * t_flux_value(i) -
573 div_base * t_mass_value) *
574 a;
575 ++t_tau_base;
576 ++t_tau_grad;
577 }
578 ++t_flux_value;
579 ++t_mass_value;
580 ++t_w;
581 }
582 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
583 PETSC_TRUE);
584 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
585 ADD_VALUES);
586 }
588 }
589
590private:
591 boost::shared_ptr<PreviousData> commonData;
593 std::map<int, BlockData> setOfBlock;
594};
595// 4. Assembly of F_v
596template <int dim>
598{
599 typedef boost::function<double(const double, const double, const double)>
601 OpAssembleStiffRhsV(std::string flux_field,
602 boost::shared_ptr<PreviousData> &data,
603 std::map<int, BlockData> &block_map,
604 FVal exact_value,
605 FVal exact_dot,
606 FVal exact_lap)
607 : OpFaceEle(flux_field, OpFaceEle::OPROW)
608 , commonData(data)
609 , setOfBlock(block_map)
610 , exactValue(exact_value)
611 , exactDot(exact_dot)
612 , exactLap(exact_lap)
613 {}
614
615 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
617 const int nb_dofs = data.getIndices().size();
618 // cerr << "In StiffRhsV ..." << endl;
619 if (nb_dofs) {
620 auto find_block_data = [&]() {
622 BlockData *block_raw_ptr = nullptr;
623 for (auto &m : setOfBlock) {
624 if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
625 block_raw_ptr = &m.second;
626 break;
627 }
628 }
629 return block_raw_ptr;
630 };
631
632 auto block_data_ptr = find_block_data();
633 if (!block_data_ptr)
634 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
635 auto &block_data = *block_data_ptr;
636
637
638 vecF.resize(nb_dofs, false);
639 vecF.clear();
640 const int nb_integration_pts = getGaussPts().size2();
641 auto t_mass_dot = getFTensor0FromVec(commonData->mass_dots);
642 auto t_flux_div = getFTensor0FromVec(commonData->flux_divs);
643 auto t_row_v_base = data.getFTensor0N();
644 auto t_w = getFTensor0IntegrationWeight();
645 const double vol = getMeasure();
646
647 const double ct = getFEMethod()->ts_t;
648 auto t_coords = getFTensor1CoordsAtGaussPts();
649 for (int gg = 0; gg < nb_integration_pts; ++gg) {
650 const double a = vol * t_w;
651 double u_dot = exactDot(t_coords(NX), t_coords(NY), ct);
652 double u_lap = exactLap(t_coords(NX), t_coords(NY), ct);
653
654 double f = u_dot - block_data.B0 * u_lap;
655 for (int rr = 0; rr < nb_dofs; ++rr) {
656 vecF[rr] += (t_row_v_base * (t_mass_dot + t_flux_div)) * a;
657 ++t_row_v_base;
658 }
659 ++t_mass_dot;
660 ++t_flux_div;
661 ++t_w;
662 ++t_coords;
663 }
664 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
665 PETSC_TRUE);
666 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
667 ADD_VALUES);
668 }
670 }
671
672private:
673 boost::shared_ptr<PreviousData> commonData;
675 std::map<int, BlockData> setOfBlock;
676
680
683};
684
685// Tangent operator
686// //**********************************************
687// 7. Tangent assembly for F_tautau excluding the essential boundary condition
688template <int dim>
689struct OpAssembleLhsTauTau : OpFaceEle // A_TauTau_1
690{
691 OpAssembleLhsTauTau(std::string flux_field,
692 boost::shared_ptr<PreviousData> &commonData,
693 std::map<int, BlockData> &block_map)
694 : OpFaceEle(flux_field, flux_field, OpFaceEle::OPROWCOL),
695 setOfBlock(block_map), commonData(commonData)
696 {
697 sYmm = true;
698 }
699
700 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
701 EntityType col_type, EntData &row_data,
702 EntData &col_data) {
704 const int nb_row_dofs = row_data.getIndices().size();
705 const int nb_col_dofs = col_data.getIndices().size();
706
707 if (nb_row_dofs && nb_col_dofs) {
708 auto find_block_data = [&]() {
710 BlockData *block_raw_ptr = nullptr;
711 for (auto &m : setOfBlock) {
712 if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
713 block_raw_ptr = &m.second;
714 break;
715 }
716 }
717 return block_raw_ptr;
718 };
719
720 auto block_data_ptr = find_block_data();
721 if (!block_data_ptr)
722 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
723 auto &block_data = *block_data_ptr;
724
725 mat.resize(nb_row_dofs, nb_col_dofs, false);
726 mat.clear();
727 const int nb_integration_pts = getGaussPts().size2();
728 auto t_mass_value = getFTensor0FromVec(commonData->mass_values);
729
730 auto t_row_tau_base = row_data.getFTensor1N<3>();
731
732 auto t_w = getFTensor0IntegrationWeight();
733 const double vol = getMeasure();
734
735 for (int gg = 0; gg != nb_integration_pts; ++gg) {
736 const double a = vol * t_w;
737 const double K = B_epsilon + (block_data.B0 + B * t_mass_value);
738 const double K_inv = 1. / K;
739 for (int rr = 0; rr != nb_row_dofs; ++rr) {
740 auto t_col_tau_base = col_data.getFTensor1N<3>(gg, 0);
741 for (int cc = 0; cc != nb_col_dofs; ++cc) {
742 mat(rr, cc) += (K_inv * t_row_tau_base(i) * t_col_tau_base(i)) * a;
743 ++t_col_tau_base;
744 }
745 ++t_row_tau_base;
746 }
747 ++t_mass_value;
748 ++t_w;
749 }
750 CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
751 ADD_VALUES);
752 if (row_side != col_side || row_type != col_type) {
753 transMat.resize(nb_col_dofs, nb_row_dofs, false);
754 noalias(transMat) = trans(mat);
755 CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
756 &transMat(0, 0), ADD_VALUES);
757 }
758 }
760 }
761
762private:
763 boost::shared_ptr<PreviousData> commonData;
766 std::map<int, BlockData> setOfBlock;
767};
768
769// 9. Assembly of tangent for F_tau_v excluding the essential bc
770template <int dim>
771struct OpAssembleLhsTauV : OpFaceEle // E_TauV
772{
773 OpAssembleLhsTauV(std::string flux_field, std::string mass_field,
774 boost::shared_ptr<PreviousData> &data,
775 std::map<int, BlockData> &block_map)
776 : OpFaceEle(flux_field, mass_field, OpFaceEle::OPROWCOL),
777 commonData(data), setOfBlock(block_map)
778 {
779 sYmm = false;
780 }
781
782 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
783 EntityType col_type, EntData &row_data,
784 EntData &col_data) {
786 const int nb_row_dofs = row_data.getIndices().size();
787 const int nb_col_dofs = col_data.getIndices().size();
788
789 if (nb_row_dofs && nb_col_dofs) {
790 auto find_block_data = [&]() {
792 BlockData *block_raw_ptr = nullptr;
793 for (auto &m : setOfBlock) {
794 if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
795 block_raw_ptr = &m.second;
796 break;
797 }
798 }
799 return block_raw_ptr;
800 };
801
802 auto block_data_ptr = find_block_data();
803 if (!block_data_ptr)
804 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
805 auto &block_data = *block_data_ptr;
806 mat.resize(nb_row_dofs, nb_col_dofs, false);
807 mat.clear();
808 const int nb_integration_pts = getGaussPts().size2();
809 auto t_w = getFTensor0IntegrationWeight();
810 auto t_row_tau_base = row_data.getFTensor1N<3>();
811
812 auto t_row_tau_grad = row_data.getFTensor2DiffN<3, 2>();
813 auto t_mass_value = getFTensor0FromVec(commonData->mass_values);
814 auto t_flux_value = getFTensor1FromMat<3>(commonData->flux_values);
815 const double vol = getMeasure();
816
817 for (int gg = 0; gg != nb_integration_pts; ++gg) {
818 const double a = vol * t_w;
819 const double K = B_epsilon + (block_data.B0 + B * t_mass_value);
820 const double K_inv = 1. / K;
821 const double K_diff = B;
822
823 for (int rr = 0; rr != nb_row_dofs; ++rr) {
824 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
825 for (int cc = 0; cc != nb_col_dofs; ++cc) {
826 double div_row_base = t_row_tau_grad(0, 0) + t_row_tau_grad(1, 1);
827 mat(rr, cc) += (-(t_row_tau_base(i) * t_flux_value(i) * K_inv *
828 K_inv * K_diff * t_col_v_base) -
829 (div_row_base * t_col_v_base)) *
830 a;
831 ++t_col_v_base;
832 }
833 ++t_row_tau_base;
834 ++t_row_tau_grad;
835 }
836 ++t_w;
837 ++t_mass_value;
838 ++t_flux_value;
839 }
840 CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
841 ADD_VALUES);
842 }
844 }
845
846private:
847 boost::shared_ptr<PreviousData> commonData;
849 std::map<int, BlockData> setOfBlock;
850};
851
852// 10. Assembly of tangent for F_v_tau
853struct OpAssembleLhsVTau : OpFaceEle // C_VTau
854{
855 OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
856 : OpFaceEle(mass_field, flux_field, OpFaceEle::OPROWCOL)
857 {
858 sYmm = false;
859 }
860
861 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
862 EntityType col_type, EntData &row_data,
863 EntData &col_data) {
865 const int nb_row_dofs = row_data.getIndices().size();
866 const int nb_col_dofs = col_data.getIndices().size();
867
868 if (nb_row_dofs && nb_col_dofs) {
869 mat.resize(nb_row_dofs, nb_col_dofs, false);
870 mat.clear();
871 const int nb_integration_pts = getGaussPts().size2();
872 auto t_w = getFTensor0IntegrationWeight();
873 auto t_row_v_base = row_data.getFTensor0N();
874 const double vol = getMeasure();
875
876 for (int gg = 0; gg != nb_integration_pts; ++gg) {
877 const double a = vol * t_w;
878 for (int rr = 0; rr != nb_row_dofs; ++rr) {
879 auto t_col_tau_grad = col_data.getFTensor2DiffN<3, 2>(gg, 0);
880 for (int cc = 0; cc != nb_col_dofs; ++cc) {
881 double div_col_base = t_col_tau_grad(0, 0) + t_col_tau_grad(1, 1);
882 mat(rr, cc) += (t_row_v_base * div_col_base) * a;
883 ++t_col_tau_grad;
884 }
885 ++t_row_v_base;
886 }
887 ++t_w;
888 }
889 CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
890 ADD_VALUES);
891 }
893 }
894
895private:
897};
898
899// 11. Assembly of tangent for F_v_v
901{
902 OpAssembleLhsVV(std::string mass_field)
903 : OpFaceEle(mass_field, mass_field, OpFaceEle::OPROWCOL)
904 {
905 sYmm = true;
906 }
907
908 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
909 EntityType col_type, EntData &row_data,
910 EntData &col_data) {
912
913 const int nb_row_dofs = row_data.getIndices().size();
914 const int nb_col_dofs = col_data.getIndices().size();
915 if (nb_row_dofs && nb_col_dofs) {
916
917 mat.resize(nb_row_dofs, nb_col_dofs, false);
918 mat.clear();
919 const int nb_integration_pts = getGaussPts().size2();
920
921 auto t_row_v_base = row_data.getFTensor0N();
922
923 auto t_w = getFTensor0IntegrationWeight();
924 const double ts_a = getFEMethod()->ts_a;
925 const double vol = getMeasure();
926
927 for (int gg = 0; gg != nb_integration_pts; ++gg) {
928 const double a = vol * t_w;
929
930 for (int rr = 0; rr != nb_row_dofs; ++rr) {
931 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
932 for (int cc = 0; cc != nb_col_dofs; ++cc) {
933 mat(rr, cc) += (ts_a * t_row_v_base * t_col_v_base) * a;
934
935 ++t_col_v_base;
936 }
937 ++t_row_v_base;
938 }
939 ++t_w;
940 }
941 CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
942 ADD_VALUES);
943 if (row_side != col_side || row_type != col_type) {
944 transMat.resize(nb_col_dofs, nb_row_dofs, false);
945 noalias(transMat) = trans(mat);
946 CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
947 &transMat(0, 0), ADD_VALUES);
948 }
949 }
951 }
952
953private:
955};
956
957struct OpError : public OpFaceEle {
958 typedef boost::function<double(const double, const double, const double)>
960 typedef boost::function<FTensor::Tensor1<double, 3>(
961 const double, const double, const double)>
963 double &eRror;
964 OpError(FVal exact_value,
965 FVal exact_lap, FGrad exact_grad,
966 boost::shared_ptr<PreviousData> &prev_data,
967 std::map<int, BlockData> &block_map,
968 double &err)
969 : OpFaceEle("ERROR", OpFaceEle::OPROW)
970 , exactVal(exact_value)
971 , exactLap(exact_lap)
972 , exactGrad(exact_grad)
973 , prevData(prev_data)
974 , setOfBlock(block_map)
975 , eRror(err)
976 {}
977 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
979 const int nb_dofs = data.getFieldData().size();
980 // cout << "nb_error_dofs : " << nb_dofs << endl;
981 if (nb_dofs) {
982 auto find_block_data = [&]() {
984 BlockData *block_raw_ptr = nullptr;
985 for (auto &m : setOfBlock) {
986 if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
987 block_raw_ptr = &m.second;
988 break;
989 }
990 }
991 return block_raw_ptr;
992 };
993
994 auto block_data_ptr = find_block_data();
995 if (!block_data_ptr)
996 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
997
998 auto &block_data = *block_data_ptr;
999
1000
1001
1002 auto t_flux_value = getFTensor1FromMat<3>(prevData->flux_values);
1003 // auto t_mass_dot = getFTensor0FromVec(prevData->mass_dots);
1004 auto t_mass_value = getFTensor0FromVec(prevData->mass_values);
1005 // auto t_flux_div = getFTensor0FromVec(prevData->flux_divs);
1006 data.getFieldData().clear();
1007 const double vol = getMeasure();
1008 const int nb_integration_pts = getGaussPts().size2();
1009 auto t_w = getFTensor0IntegrationWeight();
1010 double dt;
1011 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
1012 double ct = getFEMethod()->ts_t - dt;
1013 auto t_coords = getFTensor1CoordsAtGaussPts();
1014
1015 FTensor::Tensor1<double, 3> t_exact_flux, t_flux_error;
1016
1017 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1018 const double a = vol * t_w;
1019 double mass_exact = exactVal(t_coords(NX), t_coords(NY), ct);
1020 // double flux_lap = - block_data.B0 * exactLap(t_coords(NX), t_coords(NY), ct);
1021 t_exact_flux(i) = - block_data.B0 * exactGrad(t_coords(NX), t_coords(NY), ct)(i);
1022 t_flux_error(0) = t_flux_value(0) - t_exact_flux(0);
1023 t_flux_error(1) = t_flux_value(1) - t_exact_flux(1);
1024 t_flux_error(2) = t_flux_value(2) - t_exact_flux(2);
1025 double local_error = pow(mass_exact - t_mass_value, 2) + t_flux_error(i) * t_flux_error(i);
1026 // cout << "flux_div : " << t_flux_div << " flux_exact : " << flux_exact << endl;
1027 data.getFieldData()[0] += a * local_error;
1028 eRror += a * local_error;
1029
1030 ++t_w;
1031 ++t_mass_value;
1032 // ++t_flux_div;
1033 ++t_flux_value;
1034 // ++t_mass_dot;
1035 ++t_coords;
1036 }
1037
1038 data.getFieldDofs()[0]->getFieldData() = data.getFieldData()[0];
1039 }
1041 }
1042
1043 private:
1047 boost::shared_ptr<PreviousData> prevData;
1048 std::map<int, BlockData> setOfBlock;
1049
1052};
1053
1054// struct ExactMass : public OpFaceEle {
1055// typedef boost::function<double(const double, const double, const double)>
1056// FVal;
1057
1058// ExacMass(FVal exact_value)
1059// : OpFaceEle("EXACT_M", OpFaceEle::OPROWCOL), exactVal(exact_value),
1060// exactLap(exact_lap), exactGrad(exact_grad), prevData(prev_data),
1061// setOfBlock(block_map), eRror(err) {}
1062
1063// }
1064
1065struct Monitor : public FEMethod {
1066 double &eRror;
1067 Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj<DM> &dm,
1068 boost::shared_ptr<PostProcFaceOnRefinedMesh> &post_proc, double &err)
1069 : cOmm(comm)
1070 , rAnk(rank)
1071 , dM(dm)
1072 , postProc(post_proc)
1073 , eRror(err)
1074 {};
1079 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every_nth_step",
1080 &save_every_nth_step, PETSC_NULL);
1081 if (ts_step % save_every_nth_step == 0) {
1083 CHKERR postProc->writeFile(
1084 "out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
1085 }
1086 Vec error_per_proc;
1087 CHKERR VecCreateMPI(cOmm, 1, PETSC_DECIDE, &error_per_proc);
1088 auto get_global_error = [&]() {
1090 CHKERR VecSetValue(error_per_proc, rAnk, eRror, INSERT_VALUES);
1092 };
1093 CHKERR get_global_error();
1094 CHKERR VecAssemblyBegin(error_per_proc);
1095 CHKERR VecAssemblyEnd(error_per_proc);
1096 double error_sum;
1097 CHKERR VecSum(error_per_proc, &error_sum);
1098 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Error : %3.4e \n",
1099 error_sum);
1100 eRror = 0;
1101 // PetscPrintf(PETSC_COMM_SELF, "global_error : %3.4e\n", eRror);
1102 // eRror = 0;
1104 }
1105
1106private:
1108
1109 boost::shared_ptr<PostProcFaceOnRefinedMesh> postProc;
1110 MPI_Comm cOmm;
1111 const int rAnk;
1112};
1113
1114}; // namespace ReactionDiffusion
1115
1116#endif //__RDOPERATORS_HPP__
constexpr double a
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#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
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
Definition: heat_method.cpp:26
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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.
Definition: Templates.hpp:135
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
FTensor::Index< 'i', 3 > i
Definition: RDOperators.hpp:29
const double B_epsilon
Definition: RDOperators.hpp:21
double compute_natu_bc(const double x, const double y, const double z)
Definition: RDOperators.hpp:74
MoFEM::EdgeElementForcesAndSourcesCoreSwitch< 0 > BoundaryEle
Definition: RDOperators.hpp:11
double compute_init_val(const double x, const double y, const double z)
Definition: RDOperators.hpp:66
MoFEM::FaceElementForcesAndSourcesCoreSwitch< 0 > FaceEle
Definition: RDOperators.hpp:9
const double essen_value
Definition: RDOperators.hpp:26
DataForcesAndSourcesCore::EntData EntData
Definition: RDOperators.hpp:16
const double natu_value
Definition: RDOperators.hpp:27
const double init_value
Definition: RDOperators.hpp:25
double compute_essen_bc(const double x, const double y, const double z)
Definition: RDOperators.hpp:70
const double B
Definition: RDOperators.hpp:20
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.
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
SmartPetscObj< DM > dM
MoFEMErrorCode postProcess()
function is run at the end of loop
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)
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProc
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsTauTau(std::string flux_field, boost::shared_ptr< PreviousData > &commonData, std::map< int, BlockData > &block_map)
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)
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData
OpAssembleLhsTauV(std::string flux_field, std::string mass_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsVV(std::string mass_field)
OpAssembleNaturalBCRhsTau(std::string flux_field, Range &natural_bd_ents)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonData
boost::function< double(const double, const double, const double)> FVal
OpAssembleSlowRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &common_data, FVal exact_value, FVal exact_dot, FVal exact_lap)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> FGrad
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleStiffRhsTau(std::string flux_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData
boost::shared_ptr< PreviousData > commonData
OpAssembleStiffRhsV(std::string flux_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map, FVal exact_value, FVal exact_dot, FVal exact_lap)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
std::map< int, BlockData > setOfBlock
boost::function< double(const double, const double, const double)> FVal
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: RDOperators.hpp:89
boost::shared_ptr< PreviousData > commonData2
boost::shared_ptr< PreviousData > commonData1
std::map< int, BlockData > setOfBlock
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)
Definition: RDOperators.hpp:81
boost::shared_ptr< PreviousData > commonData3
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> FGrad
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > prevData
std::map< int, BlockData > setOfBlock
FTensor::Number< 0 > NX
boost::function< double(const double, const double, const double)> FVal
FTensor::Number< 1 > NY
OpError(FVal exact_value, FVal exact_lap, FGrad exact_grad, boost::shared_ptr< PreviousData > &prev_data, std::map< int, BlockData > &block_map, double &err)
OpEssentialBC(const std::string &flux_field, Range &essential_bd_ents)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpInitialMass(const std::string &mass_field, Range &inner_surface)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::function< double(const double, const double, const double)> ExactFunVal
boost::function< double(const double)> FVal
OpSkeletonSource(const std::string &mass_field, FVal skeleton_fun, ExactFunVal smooth_fun, Range &internal_edge_ents)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)