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