v0.14.0
Loading...
Searching...
No Matches
level_set.cpp
Go to the documentation of this file.
1/**
2 * @file level_set.cpp
3 * @example level_set.cpp
4 * @brief Implementation DG upwind method for advection/level set problem
5 * @date 2022-12-15
6 *
7 * @copyright Copyright (c) 2022
8 *
9 */
10
11#include <MoFEM.hpp>
12using namespace MoFEM;
13
14static char help[] = "...\n\n";
15
16//! [Define dimension]
18constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
19constexpr IntegrationType I =
20 IntegrationType::GAUSS; //< selected integration type
21
30
33using FaceSideEleOp = FaceSideEle::UserDataOperator;
34
36
38constexpr size_t potential_velocity_field_dim = SPACE_DIM == 2 ? 1 : 3;
39
40constexpr bool debug = false;
41constexpr int nb_levels = 3; //< number of refinement levels
42
43constexpr int start_bit =
44 nb_levels + 1; //< first refinement level for computational mesh
45
46constexpr int current_bit =
47 2 * start_bit + 1; ///< dofs bit used to do calculations
48constexpr int skeleton_bit = 2 * start_bit + 2; ///< skeleton elemets bit
49constexpr int aggregate_bit =
50 2 * start_bit + 3; ///< all bits for advection problem
51constexpr int projection_bit =
52 2 * start_bit + 4; //< bit from which data are projected
54 2 * start_bit + 5; ///< all bits for projection problem
55
56struct LevelSet {
57
58 LevelSet(MoFEM::Interface &m_field) : mField(m_field) {}
59
61
62private:
63 using VecSideArray = std::array<VectorDouble, 2>;
64 using MatSideArray = std::array<MatrixDouble, 2>;
65
66 /**
67 * @brief data structure carrying information on skeleton on both sides.
68 *
69 */
70 struct SideData {
71 // data for skeleton computation
72 std::array<EntityHandle, 2> feSideHandle;
73 std::array<VectorInt, 2>
74 indicesRowSideMap; ///< indices on rows for left hand-side
75 std::array<VectorInt, 2>
76 indicesColSideMap; ///< indices on columns for left hand-side
77 std::array<MatrixDouble, 2> rowBaseSideMap; // base functions on rows
78 std::array<MatrixDouble, 2> colBaseSideMap; // base function on columns
79 std::array<int, 2> senseMap; // orientation of local element edge/face in
80 // respect to global orientation of edge/face
81
82 VecSideArray lVec; //< Values of level set field
83 MatSideArray velMat; //< Values of velocity field
84
85 int currentFESide; ///< current side counter
86 };
87
88 /**
89 * @brief advection velocity field
90 *
91 * \note in current implementation is assumed that advection field has zero
92 * normal component.
93 *
94 * \note function define a vector velocity potential field, curl of potential
95 * field gives velocity, thus velocity is divergence free.
96 *
97 * @tparam SPACE_DIM
98 * @param x
99 * @param y
100 * @param z
101 * @return auto
102 */
103 template <int SPACE_DIM>
104 static double get_velocity_potential(double x, double y, double z);
105
106 /**
107 * @brief inital level set, i.e. advected filed
108 *
109 * @param x
110 * @param y
111 * @param z
112 * @return double
113 */
114 static double get_level_set(const double x, const double y, const double z);
115
116 /**
117 * @brief read mesh
118 *
119 * @return MoFEMErrorCode
120 */
122
123 /**
124 * @brief create fields, and set approximation order
125 *
126 * @return MoFEMErrorCode
127 */
129
130 /**
131 * @brief push operators to integrate operators on domain
132 *
133 * @return MoFEMErrorCode
134 */
136
137 /**
138 * @brief evaluate error
139 *
140 * @return MoFEMErrorCode
141 */
142 std::tuple<double, Tag> evaluateError();
143
144 /**
145 * @brief Get operator calculating velocity on coarse mesh
146 *
147 * @param vel_ptr
148 * @return DomainEleOp*
149 */
151 getZeroLevelVelOp(boost::shared_ptr<MatrixDouble> vel_ptr);
152
153 /**
154 * @brief create side element to assemble data from sides
155 *
156 * @param side_data_ptr
157 * @return boost::shared_ptr<FaceSideEle>
158 */
159 boost::shared_ptr<FaceSideEle>
160 getSideFE(boost::shared_ptr<SideData> side_data_ptr);
161
162 /**
163 * @brief push operator to integrate on skeleton
164 *
165 * @return MoFEMErrorCode
166 */
168
169 /**
170 * @brief test integration side elements
171 *
172 * Check consistency between volume and skeleton integral.
173 *
174 * @return MoFEMErrorCode
175 */
177
178 /**
179 * @brief test consistency between tangent matrix and the right hand side
180 * vectors
181 *
182 * @return MoFEMErrorCode
183 */
185
186 /**
187 * @brief initialise field set
188 *
189 * @param level_fun
190 * @return MoFEMErrorCode
191 */
193 boost::function<double(double, double, double)> level_fun =
195
196 /**
197 * @brief initialise potential velocity field
198 *
199 * @param vel_fun
200 * @return MoFEMErrorCode
201 */
203 boost::function<double(double, double, double)> vel_fun =
204 get_velocity_potential<SPACE_DIM>);
205
206 /**
207 * @brief dg level set projection
208 *
209 * @param prj_bit
210 * @param mesh_bit
211 * @return MoFEMErrorCode
212 */
213 MoFEMErrorCode dgProjection(const int prj_bit = projection_bit);
214
215 /**
216 * @brief solve advection problem
217 *
218 * @return * MoFEMErrorCode
219 */
221
222 /**
223 * @brief Wrapper executing stages while mesh refinement
224 */
226 WrapperClass() = default;
227 virtual MoFEMErrorCode setBits(LevelSet &level_set, int l) = 0;
228 virtual MoFEMErrorCode runCalcs(LevelSet &level_set, int l) = 0;
229 virtual MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l) = 0;
230 virtual double getThreshold(const double max) = 0;
231 };
232
233 /**
234 * @brief Used to execute inital mesh approximation while mesh refinement
235 *
236 */
238
239 WrapperClassInitalSolution(boost::shared_ptr<double> max_ptr)
240 : WrapperClass(), maxPtr(max_ptr) {}
241
242 MoFEMErrorCode setBits(LevelSet &level_set, int l) {
244 auto simple = level_set.mField.getInterface<Simple>();
245 simple->getBitRefLevel() =
247 simple->getBitRefLevelMask() = BitRefLevel().set();
248 simple->reSetUp(true);
250 };
251
256 }
257
259 auto bit_mng = level_set.mField.getInterface<BitRefManager>();
260 auto set_bit = [](auto l) { return BitRefLevel().set(l); };
262 Range level;
263 CHKERR bit_mng->getEntitiesByRefLevel(set_bit(start_bit + l),
264 BitRefLevel().set(), level);
266 ->synchroniseEntities(level);
267 CHKERR bit_mng->setNthBitRefLevel(current_bit, false);
268 CHKERR bit_mng->setNthBitRefLevel(level, current_bit, true);
269 CHKERR bit_mng->setNthBitRefLevel(level, aggregate_bit, true);
271 }
272
273 double getThreshold(const double max) {
274 *maxPtr = std::max(*maxPtr, max);
275 return 0.05 * (*maxPtr);
276 }
277
278 private:
279 boost::shared_ptr<double> maxPtr;
280 };
281
282 /**
283 * @brief Use peculated errors on all levels while mesh projection
284 *
285 */
287 WrapperClassErrorProjection(boost::shared_ptr<double> max_ptr)
288 : maxPtr(max_ptr) {}
289
290 MoFEMErrorCode setBits(LevelSet &level_set, int l) { return 0; };
291 MoFEMErrorCode runCalcs(LevelSet &level_set, int l) { return 0; }
293 auto bit_mng = level_set.mField.getInterface<BitRefManager>();
294 auto set_bit = [](auto l) { return BitRefLevel().set(l); };
296 Range level;
297 CHKERR bit_mng->getEntitiesByRefLevel(set_bit(start_bit + l),
298 BitRefLevel().set(), level);
300 ->synchroniseEntities(level);
301 CHKERR bit_mng->setNthBitRefLevel(current_bit, false);
302 CHKERR bit_mng->setNthBitRefLevel(level, current_bit, true);
303 CHKERR bit_mng->setNthBitRefLevel(level, aggregate_bit, true);
305 }
306 double getThreshold(const double max) { return 0.05 * (*maxPtr); }
307
308 private:
309 boost::shared_ptr<double> maxPtr;
310
311 };
312
314
315 struct OpRhsDomain; ///< integrate volume operators on rhs
316 struct OpLhsDomain; ///< integrate volume operator on lhs
317 struct OpRhsSkeleton; ///< integrate skeleton operators on rhs
318 struct OpLhsSkeleton; ///< integrate skeleton operators on khs
319
320 // Main interfaces
322
325 using OpMassLL =
327 using OpSourceL =
332 I>::OpSource<potential_velocity_field_dim, potential_velocity_field_dim>;
334 I>::OpBaseTimesScalar<1>;
335
338
340
341private:
342 boost::shared_ptr<double> maxPtr;
343
344};
345
346template <>
347double LevelSet::get_velocity_potential<2>(double x, double y, double z) {
348 return (x * x - 0.25) * (y * y - 0.25);
349}
350
351double LevelSet::get_level_set(const double x, const double y, const double z) {
352 constexpr double xc = 0.1;
353 constexpr double yc = 0.;
354 constexpr double zc = 0.;
355 constexpr double r = 0.2;
356 return std::sqrt(pow(x - xc, 2) + pow(y - yc, 2) + pow(z - zc, 2)) - r;
357}
358
363
364 if constexpr (debug) {
366 CHKERR testOp();
367 }
368
370
371 maxPtr = boost::make_shared<double>(0);
373
375 simple->getBitRefLevel() = BitRefLevel().set(skeleton_bit) |
377 simple->getBitRefLevelMask() = BitRefLevel().set();
378 simple->reSetUp(true);
379
381
383}
384
386
387 OpRhsDomain(const std::string field_name,
388 boost::shared_ptr<VectorDouble> l_ptr,
389 boost::shared_ptr<VectorDouble> l_dot_ptr,
390 boost::shared_ptr<MatrixDouble> vel_ptr);
392
393private:
394 boost::shared_ptr<VectorDouble> lPtr;
395 boost::shared_ptr<VectorDouble> lDotPtr;
396 boost::shared_ptr<MatrixDouble> velPtr;
397};
398
400 OpLhsDomain(const std::string field_name,
401 boost::shared_ptr<MatrixDouble> vel_ptr);
402 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
403
404private:
405 boost::shared_ptr<MatrixDouble> velPtr;
406};
407
409
410 OpRhsSkeleton(boost::shared_ptr<SideData> side_data_ptr,
411 boost::shared_ptr<FaceSideEle> side_fe_ptr);
412 MoFEMErrorCode doWork(int side, EntityType type,
414
415private:
416 boost::shared_ptr<SideData> sideDataPtr;
417 boost::shared_ptr<FaceSideEle>
418 sideFEPtr; ///< pointer to element to get data on edge/face sides
419
421};
422
424 OpLhsSkeleton(boost::shared_ptr<SideData> side_data_ptr,
425 boost::shared_ptr<FaceSideEle> side_fe_ptr);
426 MoFEMErrorCode doWork(int side, EntityType type,
428
429private:
430 boost::shared_ptr<SideData> sideDataPtr;
431 boost::shared_ptr<FaceSideEle>
432 sideFEPtr; ///< pointer to element to get data on edge/face sides
433
435};
436
437int main(int argc, char *argv[]) {
438
439 // Initialisation of MoFEM/PETSc and MOAB data structures
440 const char param_file[] = "param_file.petsc";
442
443 try {
444
445 // Create MoAB database
446 moab::Core moab_core;
447 moab::Interface &moab = moab_core;
448
449 // Create MoFEM database and link it to MoAB
450 MoFEM::Core mofem_core(moab);
451 MoFEM::Interface &m_field = mofem_core;
452
453 // Register DM Manager
454 DMType dm_name = "DMMOFEM";
455 CHKERR DMRegister_MoFEM(dm_name);
456
457 // Add logging channel for example
458 auto core_log = logging::core::get();
459 core_log->add_sink(
461 LogManager::setLog("LevelSet");
462 MOFEM_LOG_TAG("LevelSet", "LevelSet");
463
464 LevelSet level_set(m_field);
465 CHKERR level_set.runProblem();
466 }
468
469 // finish work cleaning memory, getting statistics, etc.
471
472 return 0;
473}
474
478 // get options from command line
479 CHKERR simple->getOptions();
480
481 // Only L2 field is set in this example. Two lines bellow forces simple
482 // interface to creat lower dimension (edge) elements, despite that fact that
483 // there is no field spanning on such elements. We need them for DG method.
484 simple->getAddSkeletonFE() = true;
485 simple->getAddBoundaryFE() = true;
486
487 // load mesh file
488 simple->getBitRefLevel() = BitRefLevel();
489 CHKERR simple->loadFile();
490
491 // Initialise bit ref levels
492 auto set_problem_bit = [&]() {
494 auto bit0 = BitRefLevel().set(start_bit);
495 BitRefLevel start_mask;
496 for (auto s = 0; s != start_bit; ++s)
497 start_mask[s] = true;
498
499 auto bit_mng = mField.getInterface<BitRefManager>();
500
501 Range level0;
502 CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(0),
503 BitRefLevel().set(), level0);
504
505 CHKERR bit_mng->setNthBitRefLevel(level0, current_bit, true);
506 CHKERR bit_mng->setNthBitRefLevel(level0, aggregate_bit, true);
507 CHKERR bit_mng->setNthBitRefLevel(level0, skeleton_bit, true);
508
509 // Set bits to build adjacencies between parents and children. That is used
510 // by simple interface.
511 simple->getBitAdjEnt() = BitRefLevel().set();
512 simple->getBitAdjParent() = BitRefLevel().set();
513 simple->getBitRefLevel() = BitRefLevel().set(current_bit);
514 simple->getBitRefLevelMask() = BitRefLevel().set();
515
516#ifndef NDEBUG
517 if constexpr (debug) {
518 auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
519 CHKERR bit_mng->writeBitLevelByDim(
520 BitRefLevel().set(0), BitRefLevel().set(), SPACE_DIM,
521 (proc_str + "level_base.vtk").c_str(), "VTK", "");
522 }
523#endif
524
526 };
527
528 CHKERR set_problem_bit();
529
531}
532
536 // Scalar fields and vector field is tested. Add more fields, i.e. vector
537 // field if needed.
538 CHKERR simple->addDomainField("L", L2, AINSWORTH_LEGENDRE_BASE, 1);
539 CHKERR simple->addDomainField("V", potential_velocity_space,
541
542 // set fields order, i.e. for most first cases order is sufficient.
543 CHKERR simple->setFieldOrder("L", 4);
544 CHKERR simple->setFieldOrder("V", 4);
545
546 // setup problem
547 CHKERR simple->setUp();
548
550}
551
555
557 boost::shared_ptr<VectorDouble> l_ptr,
558 boost::shared_ptr<VectorDouble> l_dot_ptr,
559 boost::shared_ptr<MatrixDouble> vel_ptr)
561 lPtr(l_ptr), lDotPtr(l_dot_ptr), velPtr(vel_ptr) {}
562
564 boost::shared_ptr<MatrixDouble> vel_ptr)
566 AssemblyDomainEleOp::OPROWCOL),
567 velPtr(vel_ptr) {
568 this->sYmm = false;
569}
570
572 boost::shared_ptr<SideData> side_data_ptr,
573 boost::shared_ptr<FaceSideEle> side_fe_ptr)
575 sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr) {}
576
578 boost::shared_ptr<SideData> side_data_ptr,
579 boost::shared_ptr<FaceSideEle> side_fe_ptr)
581 sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr) {}
582
585
586 const auto nb_int_points = getGaussPts().size2();
587 const auto nb_dofs = data.getIndices().size();
588 const auto nb_base_func = data.getN().size2();
589
590 auto t_l = getFTensor0FromVec(*lPtr);
591 auto t_l_dot = getFTensor0FromVec(*lDotPtr);
592 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
593
594 auto t_base = data.getFTensor0N();
595 auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
596
597 auto t_w = getFTensor0IntegrationWeight();
598 for (auto gg = 0; gg != nb_int_points; ++gg) {
599 const auto alpha = t_w * getMeasure();
600 auto res0 = alpha * t_l_dot;
602 t_res1(i) = (alpha * t_l) * t_vel(i);
603 ++t_w;
604 ++t_l;
605 ++t_l_dot;
606 ++t_vel;
607
608 auto &nf = this->locF;
609
610 int rr = 0;
611 for (; rr != nb_dofs; ++rr) {
612 nf(rr) += res0 * t_base - t_res1(i) * t_diff_base(i);
613 ++t_base;
614 ++t_diff_base;
615 }
616 for (; rr < nb_base_func; ++rr) {
617 ++t_base;
618 ++t_diff_base;
619 }
620 }
621
623}
624
626 EntData &col_data) {
628
629 const auto nb_int_points = getGaussPts().size2();
630 const auto nb_base_func = row_data.getN().size2();
631 const auto nb_row_dofs = row_data.getIndices().size();
632 const auto nb_col_dofs = col_data.getIndices().size();
633
634 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
635
636 auto t_row_base = row_data.getFTensor0N();
637 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
638
639 auto t_w = getFTensor0IntegrationWeight();
640 for (auto gg = 0; gg != nb_int_points; ++gg) {
641 const auto alpha = t_w * getMeasure();
642 const auto beta = alpha * getTSa();
643 ++t_w;
644
645 auto &mat = this->locMat;
646
647 int rr = 0;
648 for (; rr != nb_row_dofs; ++rr) {
649 auto t_col_base = col_data.getFTensor0N(gg, 0);
650 for (int cc = 0; cc != nb_col_dofs; ++cc) {
651 mat(rr, cc) +=
652 (beta * t_row_base - alpha * (t_row_diff_base(i) * t_vel(i))) *
653 t_col_base;
654 ++t_col_base;
655 }
656 ++t_row_base;
657 ++t_row_diff_base;
658 }
659 for (; rr < nb_base_func; ++rr) {
660 ++t_row_base;
661 ++t_row_diff_base;
662 }
663
664 ++t_vel;
665 }
666
668}
669
674
675 // Collect data from side domain elements
676 CHKERR loopSideFaces("dFE", *sideFEPtr);
677 const auto in_the_loop =
678 sideFEPtr->nInTheLoop; // return number of elements on the side
679
680 auto not_side = [](auto s) {
681 return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
682 };
683
684 auto get_ntensor = [](auto &base_mat) {
686 &*base_mat.data().begin());
687 };
688
689 if (in_the_loop > 0) {
690
691 // get normal of the face or edge
692 auto t_normal = getFTensor1Normal();
693 const auto nb_gauss_pts = getGaussPts().size2();
694
695 for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
696
697 // gent number of DOFs on the right side.
698 const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
699
700 if (nb_rows) {
701
702 resSkelton.resize(nb_rows, false);
703 resSkelton.clear();
704
705 // get orientation of the local element edge
706 const auto opposite_s0 = not_side(s0);
707 const auto sense_row = sideDataPtr->senseMap[s0];
708#ifndef NDEBUG
709 const auto opposite_sense_row = sideDataPtr->senseMap[opposite_s0];
710 if (sense_row * opposite_sense_row > 0)
711 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
712 "Should be opposite sign");
713#endif
714
715 // iterate the side cols
716 const auto nb_row_base_functions =
717 sideDataPtr->rowBaseSideMap[s0].size2();
718
719 auto t_w = getFTensor0IntegrationWeight();
720 auto arr_t_l =
721 make_array(getFTensor0FromVec(sideDataPtr->lVec[LEFT_SIDE]),
722 getFTensor0FromVec(sideDataPtr->lVec[RIGHT_SIDE]));
723 auto arr_t_vel = make_array(
724 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[LEFT_SIDE]),
725 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[RIGHT_SIDE]));
726
727 auto next = [&]() {
728 for (auto &t_l : arr_t_l)
729 ++t_l;
730 for (auto &t_vel : arr_t_vel)
731 ++t_vel;
732 };
733
734#ifndef NDEBUG
735 if (nb_gauss_pts != sideDataPtr->rowBaseSideMap[s0].size1())
736 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
737 "Inconsistent number of DOFs");
738#endif
739
740 auto t_row_base = get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
741 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
743 t_vel(i) = (arr_t_vel[LEFT_SIDE](i) + arr_t_vel[RIGHT_SIDE](i)) / 2.;
744 const auto dot = sense_row * (t_normal(i) * t_vel(i));
745 const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
746 const auto l_upwind = arr_t_l[l_upwind_side];
747 const auto res = t_w * dot * l_upwind;
748 next();
749 ++t_w;
750 auto rr = 0;
751 for (; rr != nb_rows; ++rr) {
752 resSkelton[rr] += t_row_base * res;
753 ++t_row_base;
754 }
755 for (; rr < nb_row_base_functions; ++rr) {
756 ++t_row_base;
757 }
758 }
759 // assemble local operator vector to global vector
760 CHKERR ::VecSetValues(getTSf(),
761 sideDataPtr->indicesRowSideMap[s0].size(),
762 &*sideDataPtr->indicesRowSideMap[s0].begin(),
763 &*resSkelton.begin(), ADD_VALUES);
764 }
765 }
766 }
767
769}
770
775
776 // Collect data from side domain elements
777 CHKERR loopSideFaces("dFE", *sideFEPtr);
778 const auto in_the_loop =
779 sideFEPtr->nInTheLoop; // return number of elements on the side
780
781 auto not_side = [](auto s) {
782 return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
783 };
784
785 auto get_ntensor = [](auto &base_mat) {
787 &*base_mat.data().begin());
788 };
789
790 if (in_the_loop > 0) {
791
792 // get normal of the face or edge
793 auto t_normal = getFTensor1Normal();
794 const auto nb_gauss_pts = getGaussPts().size2();
795
796 for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
797
798 // gent number of DOFs on the right side.
799 const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
800
801 if (nb_rows) {
802
803 // get orientation of the local element edge
804 const auto opposite_s0 = not_side(s0);
805 const auto sense_row = sideDataPtr->senseMap[s0];
806
807 // iterate the side cols
808 const auto nb_row_base_functions =
809 sideDataPtr->rowBaseSideMap[s0].size2();
810
811 for (auto s1 : {LEFT_SIDE, RIGHT_SIDE}) {
812
813 // gent number of DOFs on the right side.
814 const auto nb_cols = sideDataPtr->indicesColSideMap[s1].size();
815 const auto sense_col = sideDataPtr->senseMap[s1];
816
817 // resize local element matrix
818 matSkeleton.resize(nb_rows, nb_cols, false);
819 matSkeleton.clear();
820
821 auto t_w = getFTensor0IntegrationWeight();
822 auto arr_t_vel = make_array(
823 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[LEFT_SIDE]),
824 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[RIGHT_SIDE]));
825
826 auto next = [&]() {
827 for (auto &t_vel : arr_t_vel)
828 ++t_vel;
829 };
830
831 auto t_row_base = get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
832 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
834 t_vel(i) =
835 (arr_t_vel[LEFT_SIDE](i) + arr_t_vel[RIGHT_SIDE](i)) / 2.;
836 const auto dot = sense_row * (t_normal(i) * t_vel(i));
837 const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
838 const auto sense_upwind = sideDataPtr->senseMap[l_upwind_side];
839 auto res = t_w * dot; // * sense_row * sense_upwind;
840 next();
841 ++t_w;
842 auto rr = 0;
843 if (s1 == l_upwind_side) {
844 for (; rr != nb_rows; ++rr) {
845 auto get_ntensor = [](auto &base_mat, auto gg, auto bb) {
846 double *ptr = &base_mat(gg, bb);
848 };
849 auto t_col_base =
850 get_ntensor(sideDataPtr->colBaseSideMap[s1], gg, 0);
851 const auto res_row = res * t_row_base;
852 ++t_row_base;
853 // iterate columns
854 for (size_t cc = 0; cc != nb_cols; ++cc) {
855 matSkeleton(rr, cc) += res_row * t_col_base;
856 ++t_col_base;
857 }
858 }
859 }
860 for (; rr < nb_row_base_functions; ++rr) {
861 ++t_row_base;
862 }
863 }
864 // assemble system
865 CHKERR ::MatSetValues(getTSB(),
866 sideDataPtr->indicesRowSideMap[s0].size(),
867 &*sideDataPtr->indicesRowSideMap[s0].begin(),
868 sideDataPtr->indicesColSideMap[s1].size(),
869 &*sideDataPtr->indicesColSideMap[s1].begin(),
870 &*matSkeleton.data().begin(), ADD_VALUES);
871 }
872 }
873 }
874 }
876}
877
879LevelSet::getZeroLevelVelOp(boost::shared_ptr<MatrixDouble> vel_ptr) {
880 auto get_parent_vel_this = [&]() {
881 auto parent_fe_ptr = boost::make_shared<DomianParentEle>(mField);
883 parent_fe_ptr->getOpPtrVector(), {potential_velocity_space});
884 parent_fe_ptr->getOpPtrVector().push_back(
886 "V", vel_ptr));
887 return parent_fe_ptr;
888 };
889
890 auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr) {
891 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
892 for (int l = 0; l <= nb_levels; ++l)
893 parents_elems_ptr_vec.emplace_back(
894 boost::make_shared<DomianParentEle>(mField));
895 for (auto l = 1; l <= nb_levels; ++l) {
896 parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
897 new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
898 BitRefLevel().set(0).flip(), this_fe_ptr,
899 BitRefLevel().set(0), BitRefLevel().set()));
900 }
901 return parents_elems_ptr_vec[0];
902 };
903
904 auto this_fe_ptr = get_parent_vel_this();
905 auto parent_fe_ptr = get_parents_vel_fe_ptr(this_fe_ptr);
906 return new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
907 BitRefLevel().set(0).flip(), this_fe_ptr,
908 BitRefLevel().set(0), BitRefLevel().set());
909}
910
913 auto pip = mField.getInterface<PipelineManager>(); // get interface to
914 // pipeline manager
915
916 pip->getOpDomainLhsPipeline().clear();
917 pip->getOpDomainRhsPipeline().clear();
918
919 pip->setDomainLhsIntegrationRule([](int, int, int o) { return 3 * o; });
920 pip->setDomainRhsIntegrationRule([](int, int, int o) { return 3 * o; });
921
922 pip->getDomainLhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
923 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
925 };
926 pip->getDomainRhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
927 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
929 };
930
931 auto l_ptr = boost::make_shared<VectorDouble>();
932 auto l_dot_ptr = boost::make_shared<VectorDouble>();
933 auto vel_ptr = boost::make_shared<MatrixDouble>();
934
936 pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
937 pip->getOpDomainRhsPipeline().push_back(
938 new OpCalculateScalarFieldValues("L", l_ptr));
939 pip->getOpDomainRhsPipeline().push_back(
940 new OpCalculateScalarFieldValuesDot("L", l_dot_ptr));
941 pip->getOpDomainRhsPipeline().push_back(getZeroLevelVelOp(vel_ptr));
942 pip->getOpDomainRhsPipeline().push_back(
943 new OpRhsDomain("L", l_ptr, l_dot_ptr, vel_ptr));
944
946 pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
947 pip->getOpDomainLhsPipeline().push_back(getZeroLevelVelOp(vel_ptr));
948 pip->getOpDomainLhsPipeline().push_back(new OpLhsDomain("L", vel_ptr));
949
951}
952
953boost::shared_ptr<FaceSideEle>
954LevelSet::getSideFE(boost::shared_ptr<SideData> side_data_ptr) {
955
957
958 auto l_ptr = boost::make_shared<VectorDouble>();
959 auto vel_ptr = boost::make_shared<MatrixDouble>();
960
961 struct OpSideData : public FaceSideEleOp {
962 OpSideData(boost::shared_ptr<SideData> side_data_ptr)
963 : FaceSideEleOp("L", "L", FaceSideEleOp::OPROWCOL),
964 sideDataPtr(side_data_ptr) {
965 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
966 for (auto t = moab::CN::TypeDimensionMap[SPACE_DIM].first;
967 t <= moab::CN::TypeDimensionMap[SPACE_DIM].second; ++t)
968 doEntities[t] = true;
969 sYmm = false;
970 }
971
972 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
973 EntityType col_type, EntData &row_data,
974 EntData &col_data) {
976 if ((CN::Dimension(row_type) == SPACE_DIM) &&
977 (CN::Dimension(col_type) == SPACE_DIM)) {
978
979 auto reset = [&](auto nb_in_loop) {
980 sideDataPtr->feSideHandle[nb_in_loop] = 0;
981 sideDataPtr->indicesRowSideMap[nb_in_loop].clear();
982 sideDataPtr->indicesColSideMap[nb_in_loop].clear();
983 sideDataPtr->rowBaseSideMap[nb_in_loop].clear();
984 sideDataPtr->colBaseSideMap[nb_in_loop].clear();
985 sideDataPtr->senseMap[nb_in_loop] = 0;
986 };
987
988 const auto nb_in_loop = getFEMethod()->nInTheLoop;
989 if (nb_in_loop == 0)
990 for (auto s : {0, 1})
991 reset(s);
992
993 sideDataPtr->currentFESide = nb_in_loop;
994 sideDataPtr->senseMap[nb_in_loop] = getSkeletonSense();
995
996 } else {
997 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Should not happen");
998 }
999
1001 };
1002
1003 private:
1004 boost::shared_ptr<SideData> sideDataPtr;
1005 };
1006
1007 struct OpSideDataOnParent : public DomainEleOp {
1008
1009 OpSideDataOnParent(boost::shared_ptr<SideData> side_data_ptr,
1010 boost::shared_ptr<VectorDouble> l_ptr,
1011 boost::shared_ptr<MatrixDouble> vel_ptr)
1013 sideDataPtr(side_data_ptr), lPtr(l_ptr), velPtr(vel_ptr) {
1014 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
1015 for (auto t = moab::CN::TypeDimensionMap[SPACE_DIM].first;
1016 t <= moab::CN::TypeDimensionMap[SPACE_DIM].second; ++t)
1017 doEntities[t] = true;
1018 sYmm = false;
1019 }
1020
1021 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
1022 EntityType col_type, EntData &row_data,
1023 EntData &col_data) {
1025
1026 if ((CN::Dimension(row_type) == SPACE_DIM) &&
1027 (CN::Dimension(col_type) == SPACE_DIM)) {
1028 const auto nb_in_loop = sideDataPtr->currentFESide;
1029 sideDataPtr->feSideHandle[nb_in_loop] = getFEEntityHandle();
1030 sideDataPtr->indicesRowSideMap[nb_in_loop] = row_data.getIndices();
1031 sideDataPtr->indicesColSideMap[nb_in_loop] = col_data.getIndices();
1032 sideDataPtr->rowBaseSideMap[nb_in_loop] = row_data.getN();
1033 sideDataPtr->colBaseSideMap[nb_in_loop] = col_data.getN();
1034 (sideDataPtr->lVec)[nb_in_loop] = *lPtr;
1035 (sideDataPtr->velMat)[nb_in_loop] = *velPtr;
1036
1037#ifndef NDEBUG
1038 if ((sideDataPtr->lVec)[nb_in_loop].size() !=
1039 (sideDataPtr->velMat)[nb_in_loop].size2())
1040 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1041 "Wrong number of integaration pts %d != %d",
1042 (sideDataPtr->lVec)[nb_in_loop].size(),
1043 (sideDataPtr->velMat)[nb_in_loop].size2());
1044 if ((sideDataPtr->velMat)[nb_in_loop].size1() != SPACE_DIM)
1045 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1046 "Wrong size of velocity vector size = %d",
1047 (sideDataPtr->velMat)[nb_in_loop].size1());
1048#endif
1049
1050 if (!nb_in_loop) {
1051 (sideDataPtr->lVec)[1] = sideDataPtr->lVec[0];
1052 (sideDataPtr->velMat)[1] = (sideDataPtr->velMat)[0];
1053 } else {
1054#ifndef NDEBUG
1055 if (sideDataPtr->rowBaseSideMap[0].size1() !=
1056 sideDataPtr->rowBaseSideMap[1].size1()) {
1057 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1058 "Wrong number of integration pt %d != %d",
1059 sideDataPtr->rowBaseSideMap[0].size1(),
1060 sideDataPtr->rowBaseSideMap[1].size1());
1061 }
1062 if (sideDataPtr->colBaseSideMap[0].size1() !=
1063 sideDataPtr->colBaseSideMap[1].size1()) {
1064 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1065 "Wrong number of integration pt");
1066 }
1067#endif
1068 }
1069
1070 } else {
1071 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Should not happen");
1072 }
1073
1075 };
1076
1077 private:
1078 boost::shared_ptr<SideData> sideDataPtr;
1079 boost::shared_ptr<VectorDouble> lPtr;
1080 boost::shared_ptr<MatrixDouble> velPtr;
1081 };
1082
1083 // Calculate fields on param mesh bit element
1084 auto get_parent_this = [&]() {
1085 auto parent_fe_ptr = boost::make_shared<DomianParentEle>(mField);
1087 parent_fe_ptr->getOpPtrVector(), {potential_velocity_space, L2});
1088 parent_fe_ptr->getOpPtrVector().push_back(
1089 new OpCalculateScalarFieldValues("L", l_ptr));
1090 parent_fe_ptr->getOpPtrVector().push_back(
1091 new OpSideDataOnParent(side_data_ptr, l_ptr, vel_ptr));
1092 return parent_fe_ptr;
1093 };
1094
1095 auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
1096 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
1097 for (int l = 0; l <= nb_levels; ++l)
1098 parents_elems_ptr_vec.emplace_back(
1099 boost::make_shared<DomianParentEle>(mField));
1100 for (auto l = 1; l <= nb_levels; ++l) {
1101 parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
1102 new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
1103 BitRefLevel().set(current_bit).flip(), this_fe_ptr,
1104 BitRefLevel().set(current_bit), BitRefLevel().set()));
1105 }
1106 return parents_elems_ptr_vec[0];
1107 };
1108
1109 // Create aliased shared pointers, all elements are destroyed if side_fe_ptr
1110 // is destroyed
1111 auto get_side_fe_ptr = [&]() {
1112 auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
1113
1114 auto this_fe_ptr = get_parent_this();
1115 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
1116
1117 side_fe_ptr->getOpPtrVector().push_back(new OpSideData(side_data_ptr));
1118 side_fe_ptr->getOpPtrVector().push_back(getZeroLevelVelOp(vel_ptr));
1119 side_fe_ptr->getOpPtrVector().push_back(
1120 new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
1121 BitRefLevel().set(current_bit).flip(), this_fe_ptr,
1122 BitRefLevel().set(current_bit), BitRefLevel().set()));
1123
1124 return side_fe_ptr;
1125 };
1126
1127 return get_side_fe_ptr();
1128};
1129
1132 auto pip = mField.getInterface<PipelineManager>(); // get interface to
1133
1134 pip->getOpSkeletonLhsPipeline().clear();
1135 pip->getOpSkeletonRhsPipeline().clear();
1136
1137 pip->setSkeletonLhsIntegrationRule([](int, int, int o) { return 18; });
1138 pip->setSkeletonRhsIntegrationRule([](int, int, int o) { return 18; });
1139
1140 pip->getSkeletonLhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
1141 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1142 skeleton_bit);
1143 };
1144 pip->getSkeletonRhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
1145 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1146 skeleton_bit);
1147 };
1148
1149 auto side_data_ptr = boost::make_shared<SideData>();
1150 auto side_fe_ptr = getSideFE(side_data_ptr);
1151
1152 pip->getOpSkeletonRhsPipeline().push_back(
1153 new OpRhsSkeleton(side_data_ptr, side_fe_ptr));
1154 pip->getOpSkeletonLhsPipeline().push_back(
1155 new OpLhsSkeleton(side_data_ptr, side_fe_ptr));
1156
1158}
1159
1160std::tuple<double, Tag> LevelSet::evaluateError() {
1161
1162 struct OpErrorSkel : BoundaryEleOp {
1163
1164 OpErrorSkel(boost::shared_ptr<FaceSideEle> side_fe_ptr,
1165 boost::shared_ptr<SideData> side_data_ptr,
1166 SmartPetscObj<Vec> error_sum_ptr, Tag th_error)
1168 sideFEPtr(side_fe_ptr), sideDataPtr(side_data_ptr),
1169 errorSumPtr(error_sum_ptr), thError(th_error) {}
1170
1171 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1173
1174 // Collect data from side domain elements
1175 CHKERR loopSideFaces("dFE", *sideFEPtr);
1176 const auto in_the_loop =
1177 sideFEPtr->nInTheLoop; // return number of elements on the side
1178
1179 auto not_side = [](auto s) {
1180 return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
1181 };
1182
1183 auto nb_gauss_pts = getGaussPts().size2();
1184
1185 for (auto s : {LEFT_SIDE, RIGHT_SIDE}) {
1186
1187 auto arr_t_l =
1188 make_array(getFTensor0FromVec(sideDataPtr->lVec[LEFT_SIDE]),
1189 getFTensor0FromVec(sideDataPtr->lVec[RIGHT_SIDE]));
1190 auto arr_t_vel = make_array(
1191 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[LEFT_SIDE]),
1192 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[RIGHT_SIDE]));
1193
1194 auto next = [&]() {
1195 for (auto &t_l : arr_t_l)
1196 ++t_l;
1197 for (auto &t_vel : arr_t_vel)
1198 ++t_vel;
1199 };
1200
1201 double e = 0;
1202 auto t_w = getFTensor0IntegrationWeight();
1203 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1204 e += t_w * getMeasure() *
1205 pow(arr_t_l[LEFT_SIDE] - arr_t_l[RIGHT_SIDE], 2);
1206 next();
1207 ++t_w;
1208 }
1209 e = std::sqrt(e);
1210
1211 moab::Interface &moab =
1212 getNumeredEntFiniteElementPtr()->getBasicDataPtr()->moab;
1213 const void *tags_ptr[2];
1214 CHKERR moab.tag_get_by_ptr(thError, sideDataPtr->feSideHandle.data(), 2,
1215 tags_ptr);
1216 for (auto ff : {0, 1}) {
1217 *((double *)tags_ptr[ff]) += e;
1218 }
1219 CHKERR VecSetValue(errorSumPtr, 0, e, ADD_VALUES);
1220 };
1221
1223 }
1224
1225 private:
1226 boost::shared_ptr<FaceSideEle> sideFEPtr;
1227 boost::shared_ptr<SideData> sideDataPtr;
1228 SmartPetscObj<Vec> errorSumPtr;
1229 Tag thError;
1230 };
1231
1232 auto simple = mField.getInterface<Simple>();
1233
1234 auto error_sum_ptr = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
1235 Tag th_error;
1236 double def_val = 0;
1237 CHKERR mField.get_moab().tag_get_handle("Error", 1, MB_TYPE_DOUBLE, th_error,
1238 MB_TAG_CREAT | MB_TAG_SPARSE,
1239 &def_val);
1240
1241 auto clear_tags = [&]() {
1243 Range fe_ents;
1244 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, fe_ents);
1245 double zero;
1246 CHKERR mField.get_moab().tag_clear_data(th_error, fe_ents, &zero);
1248 };
1249
1250 auto evaluate_error = [&]() {
1252 auto skel_fe = boost::make_shared<BoundaryEle>(mField);
1253 skel_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
1254 auto side_data_ptr = boost::make_shared<SideData>();
1255 auto side_fe_ptr = getSideFE(side_data_ptr);
1256 skel_fe->getOpPtrVector().push_back(
1257 new OpErrorSkel(side_fe_ptr, side_data_ptr, error_sum_ptr, th_error));
1258 auto simple = mField.getInterface<Simple>();
1259
1260 skel_fe->exeTestHook = [&](FEMethod *fe_ptr) {
1261 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1262 skeleton_bit);
1263 };
1264
1266 simple->getSkeletonFEName(), skel_fe);
1267
1269 };
1270
1271 auto assemble_and_sum = [](auto vec) {
1272 CHK_THROW_MESSAGE(VecAssemblyBegin(vec), "assemble");
1273 CHK_THROW_MESSAGE(VecAssemblyEnd(vec), "assemble");
1274 double sum;
1275 CHK_THROW_MESSAGE(VecSum(vec, &sum), "assemble");
1276 return sum;
1277 };
1278
1279 auto propagate_error_to_parents = [&]() {
1281
1282 auto &moab = mField.get_moab();
1283 auto fe_ptr = boost::make_shared<FEMethod>();
1284 fe_ptr->exeTestHook = [&](FEMethod *fe_ptr) {
1285 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1286 current_bit);
1287 };
1288
1289 fe_ptr->preProcessHook = []() { return 0; };
1290 fe_ptr->postProcessHook = []() { return 0; };
1291 fe_ptr->operatorHook = [&]() {
1293
1294 auto fe_ent = fe_ptr->numeredEntFiniteElementPtr->getEnt();
1295 auto parent = fe_ptr->numeredEntFiniteElementPtr->getParentEnt();
1296 auto th_parent = fe_ptr->numeredEntFiniteElementPtr->getBasicDataPtr()
1297 ->th_RefParentHandle;
1298
1299 double error;
1300 CHKERR moab.tag_get_data(th_error, &fe_ent, 1, &error);
1301
1302 boost::function<MoFEMErrorCode(EntityHandle, double)> add_error =
1303 [&](auto fe_ent, auto error) {
1305 double *e_ptr;
1306 CHKERR moab.tag_get_by_ptr(th_error, &fe_ent, 1,
1307 (const void **)&e_ptr);
1308 (*e_ptr) += error;
1309
1310 EntityHandle parent;
1311 CHKERR moab.tag_get_data(th_parent, &fe_ent, 1, &parent);
1312 if (parent != fe_ent && parent)
1313 CHKERR add_error(parent, *e_ptr);
1314
1316 };
1317
1318 CHKERR add_error(parent, error);
1319
1321 };
1322
1323 CHKERR DMoFEMLoopFiniteElements(simple->getDM(), simple->getDomainFEName(),
1324 fe_ptr);
1325
1327 };
1328
1329 CHK_THROW_MESSAGE(clear_tags(), "clear error tags");
1330 CHK_THROW_MESSAGE(evaluate_error(), "evaluate error");
1331 CHK_THROW_MESSAGE(propagate_error_to_parents(), "propagate error");
1332
1333 return std::make_tuple(assemble_and_sum(error_sum_ptr), th_error);
1334}
1335
1336/**
1337 * @brief test side element
1338 *
1339 * Check consistency between volume and skeleton integral
1340 *
1341 * @return MoFEMErrorCode
1342 */
1345
1346 /**
1347 * @brief calculate volume
1348 *
1349 */
1350 struct DivergenceVol : public DomainEleOp {
1351 DivergenceVol(boost::shared_ptr<VectorDouble> l_ptr,
1352 boost::shared_ptr<MatrixDouble> vel_ptr,
1353 SmartPetscObj<Vec> div_vec)
1354 : DomainEleOp("L", DomainEleOp::OPROW), lPtr(l_ptr), velPtr(vel_ptr),
1355 divVec(div_vec) {}
1356 MoFEMErrorCode doWork(int side, EntityType type,
1359 const auto nb_dofs = data.getIndices().size();
1360 if (nb_dofs) {
1361 const auto nb_gauss_pts = getGaussPts().size2();
1362 const auto t_w = getFTensor0IntegrationWeight();
1363 auto t_diff = data.getFTensor1DiffN<SPACE_DIM>();
1364 auto t_l = getFTensor0FromVec(*lPtr);
1365 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
1366 double div = 0;
1367 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1368 for (int rr = 0; rr != nb_dofs; ++rr) {
1369 div += getMeasure() * t_w * t_l * (t_diff(i) * t_vel(i));
1370 ++t_diff;
1371 }
1372 ++t_w;
1373 ++t_l;
1374 ++t_vel;
1375 }
1376 CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
1377 }
1379 }
1380
1381 private:
1382 boost::shared_ptr<VectorDouble> lPtr;
1383 boost::shared_ptr<MatrixDouble> velPtr;
1384 SmartPetscObj<Vec> divVec;
1385 };
1386
1387 /**
1388 * @brief calculate skeleton integral
1389 *
1390 */
1391 struct DivergenceSkeleton : public BoundaryEleOp {
1392 DivergenceSkeleton(boost::shared_ptr<SideData> side_data_ptr,
1393 boost::shared_ptr<FaceSideEle> side_fe_ptr,
1394 SmartPetscObj<Vec> div_vec)
1396 sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr), divVec(div_vec) {}
1397 MoFEMErrorCode doWork(int side, EntityType type,
1400
1401 auto get_ntensor = [](auto &base_mat) {
1403 &*base_mat.data().begin());
1404 };
1405
1406 auto not_side = [](auto s) {
1407 return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
1408 };
1409
1410 // Collect data from side domain elements
1411 CHKERR loopSideFaces("dFE", *sideFEPtr);
1412 const auto in_the_loop =
1413 sideFEPtr->nInTheLoop; // return number of elements on the side
1414
1415 auto t_normal = getFTensor1Normal();
1416 const auto nb_gauss_pts = getGaussPts().size2();
1417 for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
1418 const auto nb_dofs = sideDataPtr->indicesRowSideMap[s0].size();
1419 if (nb_dofs) {
1420 auto t_base = get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
1421 auto nb_row_base_functions = sideDataPtr->rowBaseSideMap[s0].size2();
1422 auto side_sense = sideDataPtr->senseMap[s0];
1423 auto opposite_s0 = not_side(s0);
1424
1425 auto arr_t_l =
1426 make_array(getFTensor0FromVec(sideDataPtr->lVec[LEFT_SIDE]),
1427 getFTensor0FromVec(sideDataPtr->lVec[RIGHT_SIDE]));
1428 auto arr_t_vel = make_array(
1429 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[LEFT_SIDE]),
1430 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[RIGHT_SIDE]));
1431
1432 auto next = [&]() {
1433 for (auto &t_l : arr_t_l)
1434 ++t_l;
1435 for (auto &t_vel : arr_t_vel)
1436 ++t_vel;
1437 };
1438
1439 double div = 0;
1440
1441 auto t_w = getFTensor0IntegrationWeight();
1442 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1444 t_vel(i) =
1445 (arr_t_vel[LEFT_SIDE](i) + arr_t_vel[RIGHT_SIDE](i)) / 2.;
1446 const auto dot = (t_normal(i) * t_vel(i)) * side_sense;
1447 const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
1448 const auto l_upwind =
1449 arr_t_l[l_upwind_side]; //< assume that field is continues,
1450 // initialisation field has to be smooth
1451 // and exactly approximated by approx
1452 // base
1453 auto res = t_w * l_upwind * dot;
1454 ++t_w;
1455 next();
1456 int rr = 0;
1457 for (; rr != nb_dofs; ++rr) {
1458 div += t_base * res;
1459 ++t_base;
1460 }
1461 for (; rr < nb_row_base_functions; ++rr) {
1462 ++t_base;
1463 }
1464 }
1465 CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
1466 }
1467 if (!in_the_loop)
1468 break;
1469 }
1470
1472 }
1473
1474 private:
1475 boost::shared_ptr<SideData> sideDataPtr;
1476 boost::shared_ptr<FaceSideEle> sideFEPtr;
1477 boost::shared_ptr<MatrixDouble> velPtr;
1478 SmartPetscObj<Vec> divVec;
1479 };
1480
1481 auto vol_fe = boost::make_shared<DomainEle>(mField);
1482 auto skel_fe = boost::make_shared<BoundaryEle>(mField);
1483
1484 vol_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
1485 skel_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
1486
1487 auto div_vol_vec = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
1488 auto div_skel_vec = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
1489
1490 auto l_ptr = boost::make_shared<VectorDouble>();
1491 auto vel_ptr = boost::make_shared<MatrixDouble>();
1492 auto side_data_ptr = boost::make_shared<SideData>();
1493 auto side_fe_ptr = getSideFE(side_data_ptr);
1494
1496 vol_fe->getOpPtrVector(), {potential_velocity_space, L2});
1497 vol_fe->getOpPtrVector().push_back(
1498 new OpCalculateScalarFieldValues("L", l_ptr));
1499 vol_fe->getOpPtrVector().push_back(getZeroLevelVelOp(vel_ptr));
1500 vol_fe->getOpPtrVector().push_back(
1501 new DivergenceVol(l_ptr, vel_ptr, div_vol_vec));
1502
1503 skel_fe->getOpPtrVector().push_back(
1504 new DivergenceSkeleton(side_data_ptr, side_fe_ptr, div_skel_vec));
1505
1506 auto simple = mField.getInterface<Simple>();
1507 auto dm = simple->getDM();
1508
1509 /**
1510 * Set up problem such that gradient of level set field is orthogonal to
1511 * velocity field. Then volume and skeleton integral should yield the same
1512 * value.
1513 */
1514
1516 [](double x, double y, double) { return x - y; });
1518 [](double x, double y, double) { return x - y; });
1519
1520 vol_fe->exeTestHook = [&](FEMethod *fe_ptr) {
1521 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1522 current_bit);
1523 };
1524 skel_fe->exeTestHook = [&](FEMethod *fe_ptr) {
1525 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1526 skeleton_bit);
1527 };
1528
1529 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), vol_fe);
1530 CHKERR DMoFEMLoopFiniteElements(dm, simple->getSkeletonFEName(), skel_fe);
1531 CHKERR DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(), skel_fe);
1532
1533 auto assemble_and_sum = [](auto vec) {
1534 CHK_THROW_MESSAGE(VecAssemblyBegin(vec), "assemble");
1535 CHK_THROW_MESSAGE(VecAssemblyEnd(vec), "assemble");
1536 double sum;
1537 CHK_THROW_MESSAGE(VecSum(vec, &sum), "assemble");
1538 return sum;
1539 };
1540
1541 auto div_vol = assemble_and_sum(div_vol_vec);
1542 auto div_skel = assemble_and_sum(div_skel_vec);
1543
1544 auto eps = std::abs((div_vol - div_skel) / (div_vol + div_skel));
1545
1546 MOFEM_LOG("WORLD", Sev::inform) << "Testing divergence volume: " << div_vol;
1547 MOFEM_LOG("WORLD", Sev::inform) << "Testing divergence skeleton: " << div_skel
1548 << " relative difference: " << eps;
1549
1550 constexpr double eps_err = 1e-6;
1551 if (eps > eps_err)
1552 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1553 "No consistency between skeleton integral and volume integral");
1554
1556};
1557
1560
1561 // get operators tester
1562 auto simple = mField.getInterface<Simple>();
1563 auto opt = mField.getInterface<OperatorsTester>(); // get interface to
1564 // OperatorsTester
1565 auto pip = mField.getInterface<PipelineManager>(); // get interface to
1566 // pipeline manager
1567
1570
1571 auto post_proc = [&](auto dm, auto f_res, auto out_name) {
1573 auto post_proc_fe =
1574 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
1575
1577
1578 auto l_vec = boost::make_shared<VectorDouble>();
1579 post_proc_fe->getOpPtrVector().push_back(
1580 new OpCalculateScalarFieldValues("L", l_vec, f_res));
1581
1582 post_proc_fe->getOpPtrVector().push_back(
1583
1584 new OpPPMap(
1585
1586 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1587
1588 {{"L", l_vec}},
1589
1590 {},
1591
1592 {}, {})
1593
1594 );
1595
1596 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
1597 post_proc_fe);
1598 post_proc_fe->writeFile(out_name);
1600 };
1601
1602 constexpr double eps = 1e-4;
1603
1604 auto x =
1605 opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}, {"V", {-1, 1}}});
1606 auto dot_x = opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}});
1607 auto diff_x = opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}});
1608
1609 auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline,
1610 auto rhs_pipeline) {
1612
1613 auto diff_res = opt->checkCentralFiniteDifference(
1614 simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
1615 SmartPetscObj<Vec>(), diff_x, 0, 1, eps);
1616
1617 if constexpr (debug) {
1618 // Example how to plot direction in direction diff_x. If instead
1619 // directionalCentralFiniteDifference(...) diff_res is used, then error
1620 // on directive is plotted.
1621 CHKERR post_proc(simple->getDM(), diff_res, "tangent_op_error.h5m");
1622 }
1623
1624 // Calculate norm of difference between directive calculated from finite
1625 // difference, and tangent matrix.
1626 double fnorm;
1627 CHKERR VecNorm(diff_res, NORM_2, &fnorm);
1628 MOFEM_LOG_C("LevelSet", Sev::inform,
1629 "Test consistency of tangent matrix %3.4e", fnorm);
1630
1631 constexpr double err = 1e-9;
1632 if (fnorm > err)
1633 SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1634 "Norm of directional derivative too large err = %3.4e", fnorm);
1635
1637 };
1638
1639 CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
1640 pip->getDomainRhsFE());
1641 CHKERR test_domain_ops(simple->getSkeletonFEName(), pip->getSkeletonLhsFE(),
1642 pip->getSkeletonRhsFE());
1643
1645};
1646
1648 boost::function<double(double, double, double)> level_fun) {
1650
1651 // get operators tester
1652 auto simple = mField.getInterface<Simple>();
1653 auto pip = mField.getInterface<PipelineManager>(); // get interface to
1654 // pipeline manager
1655 auto prb_mng = mField.getInterface<ProblemsManager>();
1656
1657 boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(mField);
1658 boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(mField);
1659 auto swap_fe = [&]() {
1660 lhs_fe.swap(pip->getDomainLhsFE());
1661 rhs_fe.swap(pip->getDomainRhsFE());
1662 };
1663 swap_fe();
1664
1665 pip->setDomainLhsIntegrationRule([](int, int, int o) { return 3 * o; });
1666 pip->setDomainRhsIntegrationRule([](int, int, int o) { return 3 * o; });
1667
1668 auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
1669 CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "LEVELSET_POJECTION");
1670 CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
1671 CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
1672 CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
1673 CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
1674 CHKERR DMSetUp(sub_dm);
1675
1676 BitRefLevel remove_mask = BitRefLevel().set(current_bit);
1677 remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
1678 CHKERR prb_mng->removeDofsOnEntities("LEVELSET_POJECTION", "L",
1679 BitRefLevel().set(), remove_mask);
1680 auto test_bit_ele = [&](FEMethod *fe_ptr) {
1681 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1682 current_bit);
1683 };
1684 pip->getDomainLhsFE()->exeTestHook = test_bit_ele;
1685 pip->getDomainRhsFE()->exeTestHook = test_bit_ele;
1686
1688 pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
1690 pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
1691 pip->getOpDomainLhsPipeline().push_back(new OpMassLL("L", "L"));
1692 pip->getOpDomainRhsPipeline().push_back(new OpSourceL("L", level_fun));
1693
1694 CHKERR mField.getInterface<FieldBlas>()->setField(0, "L");
1695
1696 auto ksp = pip->createKSP(sub_dm);
1697 CHKERR KSPSetDM(ksp, sub_dm);
1698 CHKERR KSPSetFromOptions(ksp);
1699 CHKERR KSPSetUp(ksp);
1700
1701 auto L = createDMVector(sub_dm);
1702 auto F = vectorDuplicate(L);
1703
1704 CHKERR KSPSolve(ksp, F, L);
1705 CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
1706 CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
1707 CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
1708
1709 auto [error, th_error] = evaluateError();
1710 MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
1711#ifndef NDEBUG
1712 auto fe_meshset =
1713 mField.get_finite_element_meshset(simple->getDomainFEName());
1714 std::vector<Tag> tags{th_error};
1715 CHKERR mField.get_moab().write_file("error.h5m", "MOAB",
1716 "PARALLEL=WRITE_PART", &fe_meshset, 1,
1717 &*tags.begin(), tags.size());
1718#endif
1719
1720 auto post_proc = [&](auto dm, auto out_name, auto th_error) {
1722 auto post_proc_fe =
1723 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
1724 post_proc_fe->setTagsToTransfer({th_error});
1725 post_proc_fe->exeTestHook = test_bit_ele;
1726
1728
1729 auto l_vec = boost::make_shared<VectorDouble>();
1730 auto l_grad_mat = boost::make_shared<MatrixDouble>();
1732 post_proc_fe->getOpPtrVector(), {potential_velocity_space, L2});
1733 post_proc_fe->getOpPtrVector().push_back(
1734 new OpCalculateScalarFieldValues("L", l_vec));
1735 post_proc_fe->getOpPtrVector().push_back(
1736 new OpCalculateScalarFieldGradient<SPACE_DIM>("L", l_grad_mat));
1737
1738 post_proc_fe->getOpPtrVector().push_back(
1739
1740 new OpPPMap(
1741
1742 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1743
1744 {{"L", l_vec}},
1745
1746 {{"GradL", l_grad_mat}},
1747
1748 {}, {})
1749
1750 );
1751
1752 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
1753 post_proc_fe);
1754 post_proc_fe->writeFile(out_name);
1756 };
1757
1758 if constexpr (debug)
1759 CHKERR post_proc(sub_dm, "initial_level_set.h5m", th_error);
1760
1761 swap_fe();
1762
1764}
1765
1767 boost::function<double(double, double, double)> vel_fun) {
1769
1770 // get operators tester
1771 auto simple = mField.getInterface<Simple>();
1772 auto pip = mField.getInterface<PipelineManager>(); // get interface to
1773 // pipeline manager
1774 auto prb_mng = mField.getInterface<ProblemsManager>();
1775
1776 boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(mField);
1777 boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(mField);
1778 auto swap_fe = [&]() {
1779 lhs_fe.swap(pip->getDomainLhsFE());
1780 rhs_fe.swap(pip->getDomainRhsFE());
1781 };
1782 swap_fe();
1783
1784 pip->setDomainLhsIntegrationRule([](int, int, int o) { return 3 * o; });
1785 pip->setDomainRhsIntegrationRule([](int, int, int o) { return 3 * o; });
1786
1787 auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
1788 CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "VELOCITY_PROJECTION");
1789 CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
1790 CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
1791 CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
1792 CHKERR DMMoFEMAddSubFieldRow(sub_dm, "V");
1793 CHKERR DMSetUp(sub_dm);
1794
1795 // Velocities are calculated only on corse mesh
1796 BitRefLevel remove_mask = BitRefLevel().set(0);
1797 remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
1798 CHKERR prb_mng->removeDofsOnEntities("VELOCITY_PROJECTION", "V",
1799 BitRefLevel().set(), remove_mask);
1800
1801 auto test_bit = [&](FEMethod *fe_ptr) {
1802 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(0);
1803 };
1804 pip->getDomainLhsFE()->exeTestHook = test_bit;
1805 pip->getDomainRhsFE()->exeTestHook = test_bit;
1806
1808 pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
1810 pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
1811
1812 pip->getOpDomainLhsPipeline().push_back(new OpMassVV("V", "V"));
1813 pip->getOpDomainRhsPipeline().push_back(new OpSourceV("V", vel_fun));
1814
1815 auto ksp = pip->createKSP(sub_dm);
1816 CHKERR KSPSetDM(ksp, sub_dm);
1817 CHKERR KSPSetFromOptions(ksp);
1818 CHKERR KSPSetUp(ksp);
1819
1820 auto L = createDMVector(sub_dm);
1821 auto F = vectorDuplicate(L);
1822
1823 CHKERR KSPSolve(ksp, F, L);
1824 CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
1825 CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
1826 CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
1827
1828 auto post_proc = [&](auto dm, auto out_name) {
1830 auto post_proc_fe =
1831 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
1832 post_proc_fe->exeTestHook = test_bit;
1833
1835 post_proc_fe->getOpPtrVector(), {potential_velocity_space});
1836
1838
1839 if constexpr (SPACE_DIM == 2) {
1840 auto l_vec = boost::make_shared<VectorDouble>();
1841 auto potential_vec = boost::make_shared<VectorDouble>();
1842 auto velocity_mat = boost::make_shared<MatrixDouble>();
1843
1844 post_proc_fe->getOpPtrVector().push_back(
1845 new OpCalculateScalarFieldValues("V", potential_vec));
1846 post_proc_fe->getOpPtrVector().push_back(
1848 SPACE_DIM>("V", velocity_mat));
1849
1850 post_proc_fe->getOpPtrVector().push_back(
1851
1852 new OpPPMap(
1853
1854 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1855
1856 {{"VelocityPotential", potential_vec}},
1857
1858 {{"Velocity", velocity_mat}},
1859
1860 {}, {})
1861
1862 );
1863
1864 } else {
1865 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1866 "3d case not implemented");
1867 }
1868
1869 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
1870 post_proc_fe);
1871 post_proc_fe->writeFile(out_name);
1873 };
1874
1875 if constexpr (debug)
1876 CHKERR post_proc(sub_dm, "initial_velocity_potential.h5m");
1877
1878 swap_fe();
1879
1881}
1882
1885
1888
1889 // get operators tester
1890 auto simple = mField.getInterface<Simple>();
1891 auto pip = mField.getInterface<PipelineManager>(); // get interface to
1892 auto prb_mng = mField.getInterface<ProblemsManager>();
1893
1896
1897 auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
1898 CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "ADVECTION");
1899 CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
1900 CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
1901 CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
1902 CHKERR DMMoFEMAddElement(sub_dm, simple->getSkeletonFEName());
1903 CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
1904 CHKERR DMSetUp(sub_dm);
1905
1906 BitRefLevel remove_mask = BitRefLevel().set(current_bit);
1907 remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
1908 CHKERR prb_mng->removeDofsOnEntities("ADVECTION", "L", BitRefLevel().set(),
1909 remove_mask);
1910
1911 auto add_post_proc_fe = [&]() {
1912 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
1913
1914 Tag th_error;
1915 double def_val = 0;
1916 CHKERR mField.get_moab().tag_get_handle(
1917 "Error", 1, MB_TYPE_DOUBLE, th_error, MB_TAG_CREAT | MB_TAG_SPARSE,
1918 &def_val);
1919 post_proc_fe->setTagsToTransfer({th_error});
1920
1921 post_proc_fe->exeTestHook = [&](FEMethod *fe_ptr) {
1922 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1923 current_bit);
1924 };
1925
1927 auto l_vec = boost::make_shared<VectorDouble>();
1928 auto vel_ptr = boost::make_shared<MatrixDouble>();
1929
1931 post_proc_fe->getOpPtrVector(), {H1});
1932 post_proc_fe->getOpPtrVector().push_back(
1933 new OpCalculateScalarFieldValues("L", l_vec));
1934 post_proc_fe->getOpPtrVector().push_back(getZeroLevelVelOp(vel_ptr));
1935
1936 post_proc_fe->getOpPtrVector().push_back(
1937
1938 new OpPPMap(
1939
1940 post_proc_fe->getPostProcMesh(),
1941
1942 post_proc_fe->getMapGaussPts(),
1943
1944 {{"L", l_vec}},
1945
1946 {{"V", vel_ptr}},
1947
1948 {}, {})
1949
1950 );
1951 return post_proc_fe;
1952 };
1953
1954 auto post_proc_fe = add_post_proc_fe();
1955
1956 auto set_time_monitor = [&](auto dm, auto ts) {
1957 auto monitor_ptr = boost::make_shared<FEMethod>();
1958
1959 monitor_ptr->preProcessHook = []() { return 0; };
1960 monitor_ptr->operatorHook = []() { return 0; };
1961 monitor_ptr->postProcessHook = [&]() {
1963
1964 if (!post_proc_fe)
1965 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1966 "Null pointer for post proc element");
1967
1968 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
1969 post_proc_fe);
1970 CHKERR post_proc_fe->writeFile(
1971 "level_set_" +
1972 boost::lexical_cast<std::string>(monitor_ptr->ts_step) + ".h5m");
1974 };
1975
1976 boost::shared_ptr<FEMethod> null;
1977 DMMoFEMTSSetMonitor(sub_dm, ts, simple->getDomainFEName(), monitor_ptr,
1978 null, null);
1979
1980 return monitor_ptr;
1981 };
1982
1983 auto ts = pip->createTSIM(sub_dm);
1984
1985 auto set_solution = [&](auto ts) {
1987 auto D = createDMVector(sub_dm);
1988 CHKERR DMoFEMMeshToLocalVector(sub_dm, D, INSERT_VALUES, SCATTER_FORWARD);
1989 CHKERR TSSetSolution(ts, D);
1991 };
1992 CHKERR set_solution(ts);
1993
1994 auto monitor_pt = set_time_monitor(sub_dm, ts);
1995 CHKERR TSSetFromOptions(ts);
1996
1997 auto B = createDMMatrix(sub_dm);
1998 CHKERR TSSetIJacobian(ts, B, B, TsSetIJacobian, nullptr);
1999 level_set_raw_ptr = this;
2000
2001 CHKERR TSSetUp(ts);
2002
2003 auto ts_pre_step = [](TS ts) {
2004 auto &m_field = level_set_raw_ptr->mField;
2005 auto simple = m_field.getInterface<Simple>();
2006 auto bit_mng = m_field.getInterface<BitRefManager>();
2008
2009 auto [error, th_error] = level_set_raw_ptr->evaluateError();
2010 MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
2011
2012 auto get_norm = [&](auto x) {
2013 double nrm;
2014 CHKERR VecNorm(x, NORM_2, &nrm);
2015 return nrm;
2016 };
2017
2018 auto set_solution = [&](auto ts) {
2020 DM dm;
2021 CHKERR TSGetDM(ts, &dm);
2022 auto prb_ptr = getProblemPtr(dm);
2023
2024 auto x = createDMVector(dm);
2025 CHKERR DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_FORWARD);
2026 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
2027 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
2028
2029 MOFEM_LOG("LevelSet", Sev::inform)
2030 << "Problem " << prb_ptr->getName() << " solution vector norm "
2031 << get_norm(x);
2032 CHKERR TSSetSolution(ts, x);
2033
2035 };
2036
2037 auto refine_and_project = [&](auto ts) {
2039
2041 WrapperClassErrorProjection(level_set_raw_ptr->maxPtr));
2042 simple->getBitRefLevel() = BitRefLevel().set(skeleton_bit) |
2043 BitRefLevel().set(aggregate_bit) |
2045
2046 simple->reSetUp(true);
2047 DM dm;
2048 CHKERR TSGetDM(ts, &dm);
2050
2051 BitRefLevel remove_mask = BitRefLevel().set(current_bit);
2052 remove_mask
2053 .flip(); // DOFs which are not on bit_domain_ele should be removed
2055 ->removeDofsOnEntities("ADVECTION", "L", BitRefLevel().set(),
2056 remove_mask);
2057
2059 };
2060
2061 auto ts_reset_theta = [&](auto ts) {
2063 DM dm;
2064 CHKERR TSGetDM(ts, &dm);
2065
2066 // FIXME: Look into vec-5 how to transfer internal theta method variables
2067
2068 CHKERR TSReset(ts);
2069 CHKERR TSSetUp(ts);
2070
2072 CHKERR set_solution(ts);
2073
2074 auto B = createDMMatrix(dm);
2075 CHKERR TSSetIJacobian(ts, B, B, TsSetIJacobian, nullptr);
2076
2078 };
2079
2080 CHKERR refine_and_project(ts);
2081 CHKERR ts_reset_theta(ts);
2082
2084 };
2085
2086 auto ts_post_step = [](TS ts) { return 0; };
2087
2088 CHKERR TSSetPreStep(ts, ts_pre_step);
2089 CHKERR TSSetPostStep(ts, ts_post_step);
2090
2091 CHKERR TSSolve(ts, NULL);
2092
2094}
2095
2098
2099 auto simple = mField.getInterface<Simple>();
2100 auto bit_mng = mField.getInterface<BitRefManager>();
2101 ParallelComm *pcomm =
2102 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2103 auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
2104
2105 auto set_bit = [](auto l) { return BitRefLevel().set(l); };
2106
2107 auto save_range = [&](const std::string name, const Range &r) {
2109 auto meshset_ptr = get_temp_meshset_ptr(mField.get_moab());
2110 CHKERR mField.get_moab().add_entities(*meshset_ptr, r);
2111 CHKERR mField.get_moab().write_file(name.c_str(), "VTK", "",
2112 meshset_ptr->get_ptr(), 1);
2114 };
2115
2116 // select domain elements to refine by threshold
2117 auto get_refined_elements_meshset = [&](auto bit, auto mask) {
2118 Range fe_ents;
2119 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit, mask, SPACE_DIM, fe_ents);
2120
2121 Tag th_error;
2122 CHK_MOAB_THROW(mField.get_moab().tag_get_handle("Error", th_error),
2123 "get error handle");
2124 std::vector<double> errors(fe_ents.size());
2126 mField.get_moab().tag_get_data(th_error, fe_ents, &*errors.begin()),
2127 "get tag data");
2128 auto it = std::max_element(errors.begin(), errors.end());
2129 double max;
2130 MPI_Allreduce(&*it, &max, 1, MPI_DOUBLE, MPI_MAX, mField.get_comm());
2131 MOFEM_LOG("LevelSet", Sev::inform) << "Max error: " << max;
2132 auto threshold = wp.getThreshold(max);
2133
2134 std::vector<EntityHandle> fe_to_refine;
2135 fe_to_refine.reserve(fe_ents.size());
2136
2137 auto fe_it = fe_ents.begin();
2138 auto error_it = errors.begin();
2139 for (auto i = 0; i != fe_ents.size(); ++i) {
2140 if (*error_it > threshold) {
2141 fe_to_refine.push_back(*fe_it);
2142 }
2143 ++fe_it;
2144 ++error_it;
2145 }
2146
2147 Range ents;
2148 ents.insert_list(fe_to_refine.begin(), fe_to_refine.end());
2149 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2150 ents, nullptr, NOISY);
2151
2152 auto get_neighbours_by_bridge_vertices = [&](auto &&ents) {
2153 Range verts;
2154 CHKERR mField.get_moab().get_connectivity(ents, verts, true);
2155 CHKERR mField.get_moab().get_adjacencies(verts, SPACE_DIM, false, ents,
2156 moab::Interface::UNION);
2157 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ents);
2158 return ents;
2159 };
2160
2161 ents = get_neighbours_by_bridge_vertices(ents);
2162
2163#ifndef NDEBUG
2164 if (debug) {
2165 auto meshset_ptr = get_temp_meshset_ptr(mField.get_moab());
2166 CHK_MOAB_THROW(mField.get_moab().add_entities(*meshset_ptr, ents),
2167 "add entities to meshset");
2168 CHKERR mField.get_moab().write_file(
2169 (proc_str + "_fe_to_refine.vtk").c_str(), "VTK", "",
2170 meshset_ptr->get_ptr(), 1);
2171 }
2172#endif
2173
2174 return ents;
2175 };
2176
2177 // refine elements, and set bit ref level
2178 auto refine_mesh = [&](auto l, auto &&fe_to_refine) {
2179 Skinner skin(&mField.get_moab());
2181
2182 // get entities in "l-1" level
2183 Range level_ents;
2184 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2185 set_bit(start_bit + l - 1), BitRefLevel().set(), SPACE_DIM, level_ents);
2186 // select entities to refine
2187 fe_to_refine = intersect(level_ents, fe_to_refine);
2188 // select entities not to refine
2189 level_ents = subtract(level_ents, fe_to_refine);
2190
2191 // for entities to refine get children, i.e. redlined entities
2192 Range fe_to_refine_children;
2193 bit_mng->updateRangeByChildren(fe_to_refine, fe_to_refine_children);
2194 // add entities to to level "l"
2195 fe_to_refine_children =
2196 fe_to_refine_children.subset_by_dimension(SPACE_DIM);
2197 level_ents.merge(fe_to_refine_children);
2198
2199 auto fix_neighbour_level = [&](auto ll) {
2201 // filter entities on level ll
2202 auto level_ll = level_ents;
2203 CHKERR bit_mng->filterEntitiesByRefLevel(set_bit(ll), BitRefLevel().set(),
2204 level_ll);
2205 // find skin of ll level
2206 Range skin_edges;
2207 CHKERR skin.find_skin(0, level_ll, false, skin_edges);
2208 // get parents of skin of level ll
2209 Range skin_parents;
2210 for (auto lll = 0; lll <= ll; ++lll) {
2211 CHKERR bit_mng->updateRangeByParent(skin_edges, skin_parents);
2212 }
2213 // filter parents on level ll - 1
2214 BitRefLevel bad_bit;
2215 for (auto lll = 0; lll <= ll - 2; ++lll) {
2216 bad_bit[lll] = true;
2217 }
2218 // get adjacents to parents
2219 Range skin_adj_ents;
2220 CHKERR mField.get_moab().get_adjacencies(skin_parents, SPACE_DIM, false,
2221 skin_adj_ents,
2222 moab::Interface::UNION);
2223 CHKERR bit_mng->filterEntitiesByRefLevel(bad_bit, BitRefLevel().set(),
2224 skin_adj_ents);
2225 skin_adj_ents = intersect(skin_adj_ents, level_ents);
2226 if (!skin_adj_ents.empty()) {
2227 level_ents = subtract(level_ents, skin_adj_ents);
2228 Range skin_adj_ents_children;
2229 bit_mng->updateRangeByChildren(skin_adj_ents, skin_adj_ents_children);
2230 level_ents.merge(skin_adj_ents_children);
2231 }
2233 };
2234
2235 CHKERR fix_neighbour_level(l);
2236
2237 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2238 level_ents);
2239
2240 // get lower dimension entities for level "l"
2241 for (auto d = 0; d != SPACE_DIM; ++d) {
2242 if (d == 0) {
2243 CHKERR mField.get_moab().get_connectivity(
2244 level_ents.subset_by_dimension(SPACE_DIM), level_ents, true);
2245 } else {
2246 CHKERR mField.get_moab().get_adjacencies(
2247 level_ents.subset_by_dimension(SPACE_DIM), d, false, level_ents,
2248 moab::Interface::UNION);
2249 }
2250 }
2251 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2252 level_ents);
2253
2254 // set bit ref level to level entities
2255 CHKERR bit_mng->setNthBitRefLevel(start_bit + l, false);
2256 CHKERR bit_mng->setNthBitRefLevel(level_ents, start_bit + l, true);
2257
2258#ifndef NDEBUG
2259 auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
2260 CHKERR bit_mng->writeBitLevelByDim(
2261 set_bit(start_bit + l), BitRefLevel().set(), SPACE_DIM,
2262 (boost::lexical_cast<std::string>(l) + "_" + proc_str + "_ref_mesh.vtk")
2263 .c_str(),
2264 "VTK", "");
2265#endif
2266
2268 };
2269
2270 // set skeleton
2271 auto set_skelton_bit = [&](auto l) {
2273
2274 // get entities of dim-1 on level "l"
2275 Range level_edges;
2276 CHKERR bit_mng->getEntitiesByDimAndRefLevel(set_bit(start_bit + l),
2277 BitRefLevel().set(),
2278 SPACE_DIM - 1, level_edges);
2279
2280 // get parent of entities of level "l"
2281 Range level_edges_parents;
2282 CHKERR bit_mng->updateRangeByParent(level_edges, level_edges_parents);
2283 level_edges_parents =
2284 level_edges_parents.subset_by_dimension(SPACE_DIM - 1);
2285 CHKERR bit_mng->filterEntitiesByRefLevel(
2286 set_bit(start_bit + l), BitRefLevel().set(), level_edges_parents);
2287
2288 // skeleton entities which do not have parents
2289 auto parent_skeleton = intersect(level_edges, level_edges_parents);
2290 auto skeleton = subtract(level_edges, level_edges_parents);
2291
2292 // add adjacent domain entities
2293 CHKERR mField.get_moab().get_adjacencies(unite(parent_skeleton, skeleton),
2294 SPACE_DIM, false, skeleton,
2295 moab::Interface::UNION);
2296
2297 // set levels
2298 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(skeleton);
2299 CHKERR bit_mng->setNthBitRefLevel(skeleton_bit, false);
2300 CHKERR bit_mng->setNthBitRefLevel(skeleton, skeleton_bit, true);
2301
2302#ifndef NDEBUG
2303 CHKERR bit_mng->writeBitLevel(
2304 set_bit(skeleton_bit), BitRefLevel().set(),
2305 (boost::lexical_cast<std::string>(l) + "_" + proc_str + "_skeleton.vtk")
2306 .c_str(),
2307 "VTK", "");
2308#endif
2310 };
2311
2312 // Reset bit sand set old current and aggregate bits as projection bits
2313 Range level0_current;
2314 CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(current_bit),
2315 BitRefLevel().set(), level0_current);
2316
2317 Range level0_aggregate;
2318 CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(aggregate_bit),
2319 BitRefLevel().set(), level0_aggregate);
2320
2321 BitRefLevel start_mask;
2322 for (auto s = 0; s != start_bit; ++s)
2323 start_mask[s] = true;
2324 CHKERR bit_mng->lambdaBitRefLevel(
2325 [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
2326 CHKERR bit_mng->setNthBitRefLevel(level0_current, projection_bit, true);
2327 CHKERR bit_mng->setNthBitRefLevel(level0_aggregate, aggregate_projection_bit,
2328 true);
2329
2330 // Set zero bit ref level
2331 Range level0;
2332 CHKERR bit_mng->getEntitiesByRefLevel(set_bit(0), BitRefLevel().set(),
2333 level0);
2334 CHKERR bit_mng->setNthBitRefLevel(level0, start_bit, true);
2335 CHKERR bit_mng->setNthBitRefLevel(level0, current_bit, true);
2336 CHKERR bit_mng->setNthBitRefLevel(level0, aggregate_bit, true);
2337 CHKERR bit_mng->setNthBitRefLevel(level0, skeleton_bit, true);
2338
2339 CHKERR wp.setBits(*this, 0);
2340 CHKERR wp.runCalcs(*this, 0);
2341 for (auto l = 0; l != nb_levels; ++l) {
2342 MOFEM_LOG("WORLD", Sev::inform) << "Process level: " << l;
2343 CHKERR refine_mesh(l + 1, get_refined_elements_meshset(
2344 set_bit(start_bit + l), BitRefLevel().set()));
2345 CHKERR set_skelton_bit(l + 1);
2346 CHKERR wp.setAggregateBit(*this, l + 1);
2347 CHKERR wp.setBits(*this, l + 1);
2348 CHKERR wp.runCalcs(*this, l + 1);
2349 }
2350
2352}
2353
2356
2357 // get operators tester
2358 auto simple = mField.getInterface<Simple>();
2359 auto pip = mField.getInterface<PipelineManager>(); // get interface to
2360 auto bit_mng = mField.getInterface<BitRefManager>();
2361 auto prb_mng = mField.getInterface<ProblemsManager>();
2362
2363 auto lhs_fe = boost::make_shared<DomainEle>(mField);
2364 auto rhs_fe_prj = boost::make_shared<DomainEle>(mField);
2365 auto rhs_fe_current = boost::make_shared<DomainEle>(mField);
2366
2367 lhs_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
2368 rhs_fe_prj->getRuleHook = [](int, int, int o) { return 3 * o; };
2369 rhs_fe_current->getRuleHook = [](int, int, int o) { return 3 * o; };
2370
2371 auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
2372 CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "DG_PROJECTION");
2373 CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
2374 CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
2375 CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
2376 CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
2377 CHKERR DMSetUp(sub_dm);
2378
2379 Range current_ents;
2380 CHKERR bit_mng->getEntitiesByDimAndRefLevel(BitRefLevel().set(current_bit),
2381 BitRefLevel().set(), SPACE_DIM,
2382 current_ents);
2383 Range prj_ents;
2384 CHKERR bit_mng->getEntitiesByDimAndRefLevel(BitRefLevel().set(projection_bit),
2385 BitRefLevel().set(), SPACE_DIM,
2386 prj_ents);
2387 for (auto l = 0; l != nb_levels; ++l) {
2388 CHKERR bit_mng->updateRangeByParent(prj_ents, prj_ents);
2389 }
2390 current_ents = subtract(current_ents, prj_ents);
2391
2392 auto test_mesh_bit = [&](FEMethod *fe_ptr) {
2393 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2394 current_bit);
2395 };
2396 auto test_prj_bit = [&](FEMethod *fe_ptr) {
2397 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2399 };
2400 auto test_current_bit = [&](FEMethod *fe_ptr) {
2401 return current_ents.find(fe_ptr->getFEEntityHandle()) != current_ents.end();
2402 };
2403
2404 lhs_fe->exeTestHook = test_mesh_bit;
2405 rhs_fe_prj->exeTestHook = test_prj_bit;
2406 rhs_fe_current->exeTestHook = test_current_bit;
2407
2408 BitRefLevel remove_mask = BitRefLevel().set(current_bit);
2409 remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
2410 CHKERR prb_mng->removeDofsOnEntities(
2411 "DG_PROJECTION", "L", BitRefLevel().set(), remove_mask, nullptr, 0,
2412 MAX_DOFS_ON_ENTITY, 0, 100, NOISY, true);
2413
2415 lhs_fe->getOpPtrVector(), {potential_velocity_space, L2});
2416 lhs_fe->getOpPtrVector().push_back(new OpMassLL("L", "L"));
2417
2418 auto l_vec = boost::make_shared<VectorDouble>();
2419
2420 // This assumes that projection mesh is refined, current mesh is coarsened.
2421 auto set_prj_from_child = [&](auto rhs_fe_prj) {
2423 rhs_fe_prj->getOpPtrVector(), {potential_velocity_space, L2});
2424
2425 // Evaluate field value on projection mesh
2426 rhs_fe_prj->getOpPtrVector().push_back(
2427 new OpCalculateScalarFieldValues("L", l_vec));
2428
2429 // This element is used to assemble
2430 auto get_parent_this = [&]() {
2431 auto fe_parent_this = boost::make_shared<DomianParentEle>(mField);
2432 fe_parent_this->getOpPtrVector().push_back(
2433 new OpScalarFieldL("L", l_vec));
2434 return fe_parent_this;
2435 };
2436
2437 // Create levels of parent elements, until current element is reached, and
2438 // then assemble.
2439 auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
2440 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
2441 for (int l = 0; l <= nb_levels; ++l)
2442 parents_elems_ptr_vec.emplace_back(
2443 boost::make_shared<DomianParentEle>(mField));
2444 for (auto l = 1; l <= nb_levels; ++l) {
2445 parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
2446 new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
2447 BitRefLevel().set(current_bit).flip(), this_fe_ptr,
2448 BitRefLevel().set(current_bit),
2449 BitRefLevel().set()));
2450 }
2451 return parents_elems_ptr_vec[0];
2452 };
2453
2454 auto this_fe_ptr = get_parent_this();
2455 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
2456 rhs_fe_prj->getOpPtrVector().push_back(
2457 new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
2458 BitRefLevel().set(current_bit).flip(), this_fe_ptr,
2459 BitRefLevel().set(current_bit), BitRefLevel().set()));
2460 };
2461
2462 // This assumed that current mesh is refined, and projection mesh is coarser
2463 auto set_prj_from_parent = [&](auto rhs_fe_current) {
2464
2465 // Evaluate field on coarser element
2466 auto get_parent_this = [&]() {
2467 auto fe_parent_this = boost::make_shared<DomianParentEle>(mField);
2468 fe_parent_this->getOpPtrVector().push_back(
2469 new OpCalculateScalarFieldValues("L", l_vec));
2470 return fe_parent_this;
2471 };
2472
2473 // Create stack of evaluation on parent elements
2474 auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
2475 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
2476 for (int l = 0; l <= nb_levels; ++l)
2477 parents_elems_ptr_vec.emplace_back(
2478 boost::make_shared<DomianParentEle>(mField));
2479 for (auto l = 1; l <= nb_levels; ++l) {
2480 parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
2481 new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
2482 BitRefLevel().set(projection_bit).flip(),
2483 this_fe_ptr, BitRefLevel().set(projection_bit),
2484 BitRefLevel().set()));
2485 }
2486 return parents_elems_ptr_vec[0];
2487 };
2488
2489 auto this_fe_ptr = get_parent_this();
2490 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
2491
2492 auto reset_op_ptr = new DomainEleOp(NOSPACE, DomainEleOp::OPSPACE);
2493 reset_op_ptr->doWorkRhsHook = [&](DataOperator *op_ptr, int, EntityType,
2494 EntData &) {
2495 l_vec->resize(static_cast<DomainEleOp *>(op_ptr)->getGaussPts().size2(),
2496 false);
2497 l_vec->clear();
2498 return 0;
2499 };
2500 rhs_fe_current->getOpPtrVector().push_back(reset_op_ptr);
2501 rhs_fe_current->getOpPtrVector().push_back(
2502 new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
2503 BitRefLevel().set(projection_bit).flip(), this_fe_ptr,
2504 BitRefLevel(), BitRefLevel()));
2505
2506 // At the end assemble of current finite element
2507 rhs_fe_current->getOpPtrVector().push_back(new OpScalarFieldL("L", l_vec));
2508 };
2509
2510 set_prj_from_child(rhs_fe_prj);
2511 set_prj_from_parent(rhs_fe_current);
2512
2513 boost::shared_ptr<FEMethod> null_fe;
2514 getDMKspCtx(sub_dm)->clearLoops();
2515 CHKERR DMMoFEMKSPSetComputeOperators(sub_dm, simple->getDomainFEName(),
2516 lhs_fe, null_fe, null_fe);
2517 CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(), rhs_fe_prj,
2518 null_fe, null_fe);
2519 CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(),
2520 rhs_fe_current, null_fe, null_fe);
2521 auto ksp = MoFEM::createKSP(mField.get_comm());
2522 CHKERR KSPSetDM(ksp, sub_dm);
2523
2524 CHKERR KSPSetDM(ksp, sub_dm);
2525 CHKERR KSPSetFromOptions(ksp);
2526 CHKERR KSPSetUp(ksp);
2527
2528 auto L = createDMVector(sub_dm);
2529 auto F = vectorDuplicate(L);
2530
2531 CHKERR KSPSolve(ksp, F, L);
2532 CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
2533 CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
2534 CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
2535
2536 auto [error, th_error] = evaluateError();
2537 MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
2538
2539 auto post_proc = [&](auto dm, auto out_name, auto th_error) {
2541 auto post_proc_fe =
2542 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
2543 post_proc_fe->setTagsToTransfer({th_error});
2544 post_proc_fe->exeTestHook = test_mesh_bit;
2545
2547
2548 auto l_vec = boost::make_shared<VectorDouble>();
2549 auto l_grad_mat = boost::make_shared<MatrixDouble>();
2551 post_proc_fe->getOpPtrVector(), {potential_velocity_space, L2});
2552 post_proc_fe->getOpPtrVector().push_back(
2553 new OpCalculateScalarFieldValues("L", l_vec));
2554 post_proc_fe->getOpPtrVector().push_back(
2555 new OpCalculateScalarFieldGradient<SPACE_DIM>("L", l_grad_mat));
2556
2557 post_proc_fe->getOpPtrVector().push_back(
2558
2559 new OpPPMap(
2560
2561 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2562
2563 {{"L", l_vec}},
2564
2565 {{"GradL", l_grad_mat}},
2566
2567 {}, {})
2568
2569 );
2570
2571 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
2572 post_proc_fe);
2573 post_proc_fe->writeFile(out_name);
2575 };
2576
2577 if constexpr (debug)
2578 CHKERR post_proc(sub_dm, "dg_projection.h5m", th_error);
2579
2581}
static Index< 'o', 3 > o
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
static const double eps
constexpr int SPACE_DIM
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
BoundaryEle::UserDataOperator BoundaryEleOp
@ NOISY
Definition: definitions.h:211
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
auto get_ntensor(T &base_mat)
ElementsAndOps< SPACE_DIM >::FaceSideEle FaceSideEle
@ F
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:219
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:442
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
Definition: DMMoFEM.cpp:623
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
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:988
PetscErrorCode DMSubDMSetUp_MoFEM(DM subdm)
Definition: DMMoFEM.cpp:1314
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
Definition: DMMoFEM.cpp:664
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpSkeletonLhsPipeline()
Get the Op Skeleton Lhs Pipeline object.
virtual EntityHandle get_finite_element_meshset(const std::string name) const =0
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
double D
constexpr int skeleton_bit
skeleton elemets bit
Definition: level_set.cpp:48
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomianParentEle DomianParentEle
Definition: level_set.cpp:25
FTensor::Index< 'j', SPACE_DIM > j
Definition: level_set.cpp:553
constexpr int current_bit
dofs bit used to do calculations
Definition: level_set.cpp:46
constexpr FieldSpace potential_velocity_space
Definition: level_set.cpp:37
static char help[]
Definition: level_set.cpp:14
FTensor::Index< 'k', SPACE_DIM > k
Definition: level_set.cpp:554
constexpr int nb_levels
Definition: level_set.cpp:41
FTensor::Index< 'i', SPACE_DIM > i
Definition: level_set.cpp:552
constexpr int SPACE_DIM
[Define dimension]
Definition: level_set.cpp:17
constexpr bool debug
Definition: level_set.cpp:40
constexpr int start_bit
Definition: level_set.cpp:43
FaceSideEle::UserDataOperator FaceSideEleOp
Definition: level_set.cpp:33
constexpr int aggregate_projection_bit
all bits for projection problem
Definition: level_set.cpp:53
constexpr IntegrationType I
Definition: level_set.cpp:19
DomainEle::UserDataOperator DomainEleOp
Definition: level_set.cpp:31
constexpr AssemblyType A
Definition: level_set.cpp:18
LevelSet * level_set_raw_ptr
Definition: level_set.cpp:1883
constexpr int aggregate_bit
all bits for advection problem
Definition: level_set.cpp:49
constexpr int projection_bit
Definition: level_set.cpp:51
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1884
constexpr size_t potential_velocity_field_dim
Definition: level_set.cpp:38
FTensor::Index< 'l', 3 > l
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto createKSP(MPI_Comm comm)
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1042
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobina in TS solver.
Definition: TsCtx.cpp:131
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:424
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
Definition: DMMoFEM.hpp:1017
constexpr auto make_array(Arg &&...arg)
Create Array.
Definition: Templates.hpp:1871
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1767
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition: DMMoFEM.hpp:976
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
OpLhsDomain(const std::string field_name, boost::shared_ptr< MatrixDouble > vel_ptr)
Definition: level_set.cpp:563
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.
Definition: level_set.cpp:625
boost::shared_ptr< MatrixDouble > velPtr
Definition: level_set.cpp:405
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
Definition: level_set.cpp:432
boost::shared_ptr< SideData > sideDataPtr
Definition: level_set.cpp:430
MatrixDouble matSkeleton
Definition: level_set.cpp:434
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: level_set.cpp:772
OpLhsSkeleton(boost::shared_ptr< SideData > side_data_ptr, boost::shared_ptr< FaceSideEle > side_fe_ptr)
Definition: level_set.cpp:577
boost::shared_ptr< VectorDouble > lPtr
Definition: level_set.cpp:394
MoFEMErrorCode iNtegrate(EntData &data)
Class dedicated to integrate operator.
Definition: level_set.cpp:583
boost::shared_ptr< MatrixDouble > velPtr
Definition: level_set.cpp:396
OpRhsDomain(const std::string field_name, boost::shared_ptr< VectorDouble > l_ptr, boost::shared_ptr< VectorDouble > l_dot_ptr, boost::shared_ptr< MatrixDouble > vel_ptr)
Definition: level_set.cpp:556
boost::shared_ptr< VectorDouble > lDotPtr
Definition: level_set.cpp:395
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: level_set.cpp:671
VectorDouble resSkelton
Definition: level_set.cpp:420
boost::shared_ptr< SideData > sideDataPtr
Definition: level_set.cpp:416
OpRhsSkeleton(boost::shared_ptr< SideData > side_data_ptr, boost::shared_ptr< FaceSideEle > side_fe_ptr)
Definition: level_set.cpp:571
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
Definition: level_set.cpp:418
data structure carrying information on skeleton on both sides.
Definition: level_set.cpp:70
std::array< VectorInt, 2 > indicesColSideMap
indices on columns for left hand-side
Definition: level_set.cpp:76
MatSideArray velMat
Definition: level_set.cpp:83
std::array< EntityHandle, 2 > feSideHandle
Definition: level_set.cpp:72
int currentFESide
current side counter
Definition: level_set.cpp:85
std::array< MatrixDouble, 2 > colBaseSideMap
Definition: level_set.cpp:78
std::array< MatrixDouble, 2 > rowBaseSideMap
Definition: level_set.cpp:77
std::array< VectorInt, 2 > indicesRowSideMap
indices on rows for left hand-side
Definition: level_set.cpp:74
VecSideArray lVec
Definition: level_set.cpp:82
std::array< int, 2 > senseMap
Definition: level_set.cpp:79
Use peculated errors on all levels while mesh projection.
Definition: level_set.cpp:286
WrapperClassErrorProjection(boost::shared_ptr< double > max_ptr)
Definition: level_set.cpp:287
MoFEMErrorCode runCalcs(LevelSet &level_set, int l)
Definition: level_set.cpp:291
MoFEMErrorCode setBits(LevelSet &level_set, int l)
Definition: level_set.cpp:290
double getThreshold(const double max)
Definition: level_set.cpp:306
MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l)
Definition: level_set.cpp:292
boost::shared_ptr< double > maxPtr
Definition: level_set.cpp:309
Wrapper executing stages while mesh refinement.
Definition: level_set.cpp:225
virtual double getThreshold(const double max)=0
virtual MoFEMErrorCode runCalcs(LevelSet &level_set, int l)=0
virtual MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l)=0
virtual MoFEMErrorCode setBits(LevelSet &level_set, int l)=0
Used to execute inital mesh approximation while mesh refinement.
Definition: level_set.cpp:237
double getThreshold(const double max)
Definition: level_set.cpp:273
MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l)
Definition: level_set.cpp:258
boost::shared_ptr< double > maxPtr
Definition: level_set.cpp:279
MoFEMErrorCode setBits(LevelSet &level_set, int l)
Definition: level_set.cpp:242
WrapperClassInitalSolution(boost::shared_ptr< double > max_ptr)
Definition: level_set.cpp:239
MoFEMErrorCode runCalcs(LevelSet &level_set, int l)
Definition: level_set.cpp:252
MoFEM::Interface & mField
integrate skeleton operators on khs
Definition: level_set.cpp:321
MoFEMErrorCode solveAdvection()
solve advection problem
Definition: level_set.cpp:1886
boost::shared_ptr< double > maxPtr
Definition: level_set.cpp:342
std::array< MatrixDouble, 2 > MatSideArray
Definition: level_set.cpp:64
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< 1, 1 > OpMassLL
Definition: level_set.cpp:326
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpSource< potential_velocity_field_dim, potential_velocity_field_dim > OpSourceV
Definition: level_set.cpp:332
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< potential_velocity_field_dim, potential_velocity_field_dim > OpMassVV
Definition: level_set.cpp:330
MoFEMErrorCode readMesh()
read mesh
Definition: level_set.cpp:475
MoFEMErrorCode initialiseFieldLevelSet(boost::function< double(double, double, double)> level_fun=get_level_set)
initialise field set
Definition: level_set.cpp:1647
MoFEMErrorCode setupProblem()
create fields, and set approximation order
Definition: level_set.cpp:533
boost::shared_ptr< FaceSideEle > getSideFE(boost::shared_ptr< SideData > side_data_ptr)
create side element to assemble data from sides
Definition: level_set.cpp:954
MoFEMErrorCode pushOpDomain()
push operators to integrate operators on domain
Definition: level_set.cpp:911
MoFEMErrorCode runProblem()
Definition: level_set.cpp:359
LevelSet(MoFEM::Interface &m_field)
Definition: level_set.cpp:58
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< 1 > OpScalarFieldL
Definition: level_set.cpp:334
static double get_velocity_potential(double x, double y, double z)
advection velocity field
std::array< VectorDouble, 2 > VecSideArray
Definition: level_set.cpp:63
MoFEMErrorCode pushOpSkeleton()
push operator to integrate on skeleton
Definition: level_set.cpp:1130
MoFEMErrorCode testOp()
test consistency between tangent matrix and the right hand side vectors
Definition: level_set.cpp:1558
MoFEMErrorCode dgProjection(const int prj_bit=projection_bit)
dg level set projection
Definition: level_set.cpp:2354
MoFEMErrorCode initialiseFieldVelocity(boost::function< double(double, double, double)> vel_fun=get_velocity_potential< SPACE_DIM >)
initialise potential velocity field
Definition: level_set.cpp:1766
static double get_level_set(const double x, const double y, const double z)
inital level set, i.e. advected filed
Definition: level_set.cpp:351
MoFEMErrorCode refineMesh(WrapperClass &&wp)
Definition: level_set.cpp:2096
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpSource< 1, 1 > OpSourceL
Definition: level_set.cpp:328
ForcesAndSourcesCore::UserDataOperator * getZeroLevelVelOp(boost::shared_ptr< MatrixDouble > vel_ptr)
Get operator calculating velocity on coarse mesh.
Definition: level_set.cpp:879
std::tuple< double, Tag > evaluateError()
evaluate error
Definition: level_set.cpp:1160
MoFEMErrorCode testSideFE()
test integration side elements
Definition: level_set.cpp:1343
Add operators pushing bases from local to physical configuration.
Managing BitRefLevels.
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
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.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure for User Loop Methods on finite elements
Basic algebra on fields.
Definition: FieldBlas.hpp:21
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
@ OPSPACE
operator do Work is execute on space data
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
Calculate curl of vector field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
Definition: Simple.hpp:313
intrusive_ptr for managing petsc objects
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:15
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
auto save_range