14static char help[] =
"...\n\n";
20 IntegrationType::GAUSS;
73 std::array<VectorInt, 2>
75 std::array<VectorInt, 2>
103 template <
int SPACE_DIM>
114 static double get_level_set(
const double x,
const double y,
const double z);
159 boost::shared_ptr<FaceSideEle>
160 getSideFE(boost::shared_ptr<SideData> side_data_ptr);
193 boost::function<
double(
double,
double,
double)> level_fun =
203 boost::function<
double(
double,
double,
double)> vel_fun =
204 get_velocity_potential<SPACE_DIM>);
245 simple->getBitRefLevel() =
266 ->synchroniseEntities(level);
275 return 0.05 * (*maxPtr);
300 ->synchroniseEntities(level);
332 I>::OpSource<potential_velocity_field_dim, potential_velocity_field_dim>;
334 I>::OpBaseTimesScalar<1>;
347double LevelSet::get_velocity_potential<2>(
double x,
double y,
double z) {
348 return (x * x - 0.25) * (y * y - 0.25);
352 constexpr double xc = 0.1;
353 constexpr double yc = 0.;
354 constexpr double zc = 0.;
355 constexpr double r = 0.2;
356 return std::sqrt(pow(x - xc, 2) + pow(y - yc, 2) + pow(z - zc, 2)) - r;
364 if constexpr (
debug) {
371 maxPtr = boost::make_shared<double>(0);
388 boost::shared_ptr<VectorDouble> l_ptr,
389 boost::shared_ptr<VectorDouble> l_dot_ptr,
390 boost::shared_ptr<MatrixDouble> vel_ptr);
394 boost::shared_ptr<VectorDouble>
lPtr;
401 boost::shared_ptr<MatrixDouble> vel_ptr);
411 boost::shared_ptr<FaceSideEle> side_fe_ptr);
417 boost::shared_ptr<FaceSideEle>
425 boost::shared_ptr<FaceSideEle> side_fe_ptr);
431 boost::shared_ptr<FaceSideEle>
437int main(
int argc,
char *argv[]) {
446 moab::Core moab_core;
447 moab::Interface &moab = moab_core;
454 DMType dm_name =
"DMMOFEM";
458 auto core_log = logging::core::get();
484 simple->getAddSkeletonFE() =
true;
485 simple->getAddBoundaryFE() =
true;
492 auto set_problem_bit = [&]() {
497 start_mask[s] =
true;
517 if constexpr (
debug) {
519 CHKERR bit_mng->writeBitLevelByDim(
521 (proc_str +
"level_base.vtk").c_str(),
"VTK",
"");
557 boost::shared_ptr<VectorDouble> l_ptr,
558 boost::shared_ptr<VectorDouble> l_dot_ptr,
559 boost::shared_ptr<MatrixDouble> vel_ptr)
561 lPtr(l_ptr), lDotPtr(l_dot_ptr), velPtr(vel_ptr) {}
564 boost::shared_ptr<MatrixDouble> vel_ptr)
572 boost::shared_ptr<SideData> side_data_ptr,
573 boost::shared_ptr<FaceSideEle> side_fe_ptr)
575 sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr) {}
578 boost::shared_ptr<SideData> side_data_ptr,
579 boost::shared_ptr<FaceSideEle> side_fe_ptr)
581 sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr) {}
586 const auto nb_int_points = getGaussPts().size2();
587 const auto nb_dofs = data.
getIndices().size();
588 const auto nb_base_func = data.
getN().size2();
592 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
597 auto t_w = getFTensor0IntegrationWeight();
598 for (
auto gg = 0; gg != nb_int_points; ++gg) {
599 const auto alpha = t_w * getMeasure();
600 auto res0 = alpha * t_l_dot;
602 t_res1(
i) = (alpha * t_l) * t_vel(
i);
608 auto &nf = this->locF;
611 for (; rr != nb_dofs; ++rr) {
612 nf(rr) += res0 * t_base - t_res1(
i) * t_diff_base(
i);
616 for (; rr < nb_base_func; ++rr) {
629 const auto nb_int_points = getGaussPts().size2();
630 const auto nb_base_func = row_data.
getN().size2();
631 const auto nb_row_dofs = row_data.
getIndices().size();
632 const auto nb_col_dofs = col_data.
getIndices().size();
634 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
639 auto t_w = getFTensor0IntegrationWeight();
640 for (
auto gg = 0; gg != nb_int_points; ++gg) {
641 const auto alpha = t_w * getMeasure();
642 const auto beta = alpha * getTSa();
645 auto &mat = this->locMat;
648 for (; rr != nb_row_dofs; ++rr) {
650 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
652 (beta * t_row_base - alpha * (t_row_diff_base(
i) * t_vel(
i))) *
659 for (; rr < nb_base_func; ++rr) {
676 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
677 const auto in_the_loop =
678 sideFEPtr->nInTheLoop;
680 auto not_side = [](
auto s) {
686 &*base_mat.data().begin());
689 if (in_the_loop > 0) {
692 auto t_normal = getFTensor1Normal();
693 const auto nb_gauss_pts = getGaussPts().size2();
698 const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
702 resSkelton.resize(nb_rows,
false);
706 const auto opposite_s0 = not_side(s0);
707 const auto sense_row = sideDataPtr->senseMap[s0];
709 const auto opposite_sense_row = sideDataPtr->senseMap[opposite_s0];
710 if (sense_row * opposite_sense_row > 0)
712 "Should be opposite sign");
716 const auto nb_row_base_functions =
717 sideDataPtr->rowBaseSideMap[s0].size2();
719 auto t_w = getFTensor0IntegrationWeight();
724 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
725 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
728 for (
auto &t_l : arr_t_l)
730 for (
auto &t_vel : arr_t_vel)
735 if (nb_gauss_pts != sideDataPtr->rowBaseSideMap[s0].size1())
737 "Inconsistent number of DOFs");
740 auto t_row_base =
get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
741 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
744 const auto dot = sense_row * (t_normal(
i) * t_vel(
i));
745 const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
746 const auto l_upwind = arr_t_l[l_upwind_side];
747 const auto res = t_w * dot * l_upwind;
751 for (; rr != nb_rows; ++rr) {
752 resSkelton[rr] += t_row_base * res;
755 for (; rr < nb_row_base_functions; ++rr) {
760 CHKERR ::VecSetValues(getTSf(),
761 sideDataPtr->indicesRowSideMap[s0].size(),
762 &*sideDataPtr->indicesRowSideMap[s0].begin(),
763 &*resSkelton.begin(), ADD_VALUES);
777 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
778 const auto in_the_loop =
779 sideFEPtr->nInTheLoop;
781 auto not_side = [](
auto s) {
787 &*base_mat.data().begin());
790 if (in_the_loop > 0) {
793 auto t_normal = getFTensor1Normal();
794 const auto nb_gauss_pts = getGaussPts().size2();
799 const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
804 const auto opposite_s0 = not_side(s0);
805 const auto sense_row = sideDataPtr->senseMap[s0];
808 const auto nb_row_base_functions =
809 sideDataPtr->rowBaseSideMap[s0].size2();
814 const auto nb_cols = sideDataPtr->indicesColSideMap[s1].size();
815 const auto sense_col = sideDataPtr->senseMap[s1];
818 matSkeleton.resize(nb_rows, nb_cols,
false);
821 auto t_w = getFTensor0IntegrationWeight();
823 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
824 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
827 for (
auto &t_vel : arr_t_vel)
831 auto t_row_base =
get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
832 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
836 const auto dot = sense_row * (t_normal(
i) * t_vel(
i));
837 const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
838 const auto sense_upwind = sideDataPtr->senseMap[l_upwind_side];
839 auto res = t_w * dot;
843 if (s1 == l_upwind_side) {
844 for (; rr != nb_rows; ++rr) {
845 auto get_ntensor = [](
auto &base_mat,
auto gg,
auto bb) {
846 double *ptr = &base_mat(gg, bb);
850 get_ntensor(sideDataPtr->colBaseSideMap[s1], gg, 0);
851 const auto res_row = res * t_row_base;
854 for (
size_t cc = 0; cc != nb_cols; ++cc) {
855 matSkeleton(rr, cc) += res_row * t_col_base;
860 for (; rr < nb_row_base_functions; ++rr) {
865 CHKERR ::MatSetValues(getTSB(),
866 sideDataPtr->indicesRowSideMap[s0].size(),
867 &*sideDataPtr->indicesRowSideMap[s0].begin(),
868 sideDataPtr->indicesColSideMap[s1].size(),
869 &*sideDataPtr->indicesColSideMap[s1].begin(),
870 &*matSkeleton.data().begin(), ADD_VALUES);
880 auto get_parent_vel_this = [&]() {
881 auto parent_fe_ptr = boost::make_shared<DomianParentEle>(
mField);
883 parent_fe_ptr->getOpPtrVector(), {potential_velocity_space});
884 parent_fe_ptr->getOpPtrVector().push_back(
887 return parent_fe_ptr;
890 auto get_parents_vel_fe_ptr = [&](
auto this_fe_ptr) {
891 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
893 parents_elems_ptr_vec.emplace_back(
894 boost::make_shared<DomianParentEle>(
mField));
896 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
901 return parents_elems_ptr_vec[0];
904 auto this_fe_ptr = get_parent_vel_this();
905 auto parent_fe_ptr = get_parents_vel_fe_ptr(this_fe_ptr);
917 pip->getOpDomainRhsPipeline().clear();
919 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
920 pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
922 pip->getDomainLhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
923 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
926 pip->getDomainRhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
927 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
931 auto l_ptr = boost::make_shared<VectorDouble>();
932 auto l_dot_ptr = boost::make_shared<VectorDouble>();
933 auto vel_ptr = boost::make_shared<MatrixDouble>();
936 pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
937 pip->getOpDomainRhsPipeline().push_back(
939 pip->getOpDomainRhsPipeline().push_back(
942 pip->getOpDomainRhsPipeline().push_back(
946 pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
948 pip->getOpDomainLhsPipeline().push_back(
new OpLhsDomain(
"L", vel_ptr));
953boost::shared_ptr<FaceSideEle>
958 auto l_ptr = boost::make_shared<VectorDouble>();
959 auto vel_ptr = boost::make_shared<MatrixDouble>();
962 OpSideData(boost::shared_ptr<SideData> side_data_ptr)
964 sideDataPtr(side_data_ptr) {
965 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE],
false);
966 for (
auto t = moab::CN::TypeDimensionMap[
SPACE_DIM].first;
967 t <= moab::CN::TypeDimensionMap[
SPACE_DIM].second; ++
t)
968 doEntities[
t] =
true;
976 if ((CN::Dimension(row_type) ==
SPACE_DIM) &&
977 (CN::Dimension(col_type) ==
SPACE_DIM)) {
979 auto reset = [&](
auto nb_in_loop) {
980 sideDataPtr->feSideHandle[nb_in_loop] = 0;
981 sideDataPtr->indicesRowSideMap[nb_in_loop].clear();
982 sideDataPtr->indicesColSideMap[nb_in_loop].clear();
983 sideDataPtr->rowBaseSideMap[nb_in_loop].clear();
984 sideDataPtr->colBaseSideMap[nb_in_loop].clear();
985 sideDataPtr->senseMap[nb_in_loop] = 0;
988 const auto nb_in_loop = getFEMethod()->nInTheLoop;
990 for (
auto s : {0, 1})
993 sideDataPtr->currentFESide = nb_in_loop;
994 sideDataPtr->senseMap[nb_in_loop] = getSkeletonSense();
1004 boost::shared_ptr<SideData> sideDataPtr;
1009 OpSideDataOnParent(boost::shared_ptr<SideData> side_data_ptr,
1010 boost::shared_ptr<VectorDouble> l_ptr,
1011 boost::shared_ptr<MatrixDouble> vel_ptr)
1013 sideDataPtr(side_data_ptr), lPtr(l_ptr), velPtr(vel_ptr) {
1014 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE],
false);
1015 for (
auto t = moab::CN::TypeDimensionMap[
SPACE_DIM].first;
1016 t <= moab::CN::TypeDimensionMap[
SPACE_DIM].second; ++
t)
1017 doEntities[
t] =
true;
1026 if ((CN::Dimension(row_type) ==
SPACE_DIM) &&
1027 (CN::Dimension(col_type) ==
SPACE_DIM)) {
1028 const auto nb_in_loop = sideDataPtr->currentFESide;
1029 sideDataPtr->feSideHandle[nb_in_loop] = getFEEntityHandle();
1030 sideDataPtr->indicesRowSideMap[nb_in_loop] = row_data.
getIndices();
1031 sideDataPtr->indicesColSideMap[nb_in_loop] = col_data.
getIndices();
1032 sideDataPtr->rowBaseSideMap[nb_in_loop] = row_data.
getN();
1033 sideDataPtr->colBaseSideMap[nb_in_loop] = col_data.
getN();
1034 (sideDataPtr->lVec)[nb_in_loop] = *lPtr;
1035 (sideDataPtr->velMat)[nb_in_loop] = *velPtr;
1038 if ((sideDataPtr->lVec)[nb_in_loop].size() !=
1039 (sideDataPtr->velMat)[nb_in_loop].size2())
1041 "Wrong number of integaration pts %d != %d",
1042 (sideDataPtr->lVec)[nb_in_loop].size(),
1043 (sideDataPtr->velMat)[nb_in_loop].size2());
1044 if ((sideDataPtr->velMat)[nb_in_loop].size1() !=
SPACE_DIM)
1046 "Wrong size of velocity vector size = %d",
1047 (sideDataPtr->velMat)[nb_in_loop].size1());
1051 (sideDataPtr->lVec)[1] = sideDataPtr->lVec[0];
1052 (sideDataPtr->velMat)[1] = (sideDataPtr->velMat)[0];
1055 if (sideDataPtr->rowBaseSideMap[0].size1() !=
1056 sideDataPtr->rowBaseSideMap[1].size1()) {
1058 "Wrong number of integration pt %d != %d",
1059 sideDataPtr->rowBaseSideMap[0].size1(),
1060 sideDataPtr->rowBaseSideMap[1].size1());
1062 if (sideDataPtr->colBaseSideMap[0].size1() !=
1063 sideDataPtr->colBaseSideMap[1].size1()) {
1065 "Wrong number of integration pt");
1078 boost::shared_ptr<SideData> sideDataPtr;
1079 boost::shared_ptr<VectorDouble> lPtr;
1080 boost::shared_ptr<MatrixDouble> velPtr;
1084 auto get_parent_this = [&]() {
1085 auto parent_fe_ptr = boost::make_shared<DomianParentEle>(
mField);
1087 parent_fe_ptr->getOpPtrVector(), {potential_velocity_space, L2});
1088 parent_fe_ptr->getOpPtrVector().push_back(
1090 parent_fe_ptr->getOpPtrVector().push_back(
1091 new OpSideDataOnParent(side_data_ptr, l_ptr, vel_ptr));
1092 return parent_fe_ptr;
1095 auto get_parents_fe_ptr = [&](
auto this_fe_ptr) {
1096 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
1098 parents_elems_ptr_vec.emplace_back(
1099 boost::make_shared<DomianParentEle>(
mField));
1101 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
1106 return parents_elems_ptr_vec[0];
1111 auto get_side_fe_ptr = [&]() {
1112 auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
1114 auto this_fe_ptr = get_parent_this();
1115 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
1117 side_fe_ptr->getOpPtrVector().push_back(
new OpSideData(side_data_ptr));
1119 side_fe_ptr->getOpPtrVector().push_back(
1127 return get_side_fe_ptr();
1135 pip->getOpSkeletonRhsPipeline().clear();
1137 pip->setSkeletonLhsIntegrationRule([](
int,
int,
int o) {
return 18; });
1138 pip->setSkeletonRhsIntegrationRule([](
int,
int,
int o) {
return 18; });
1140 pip->getSkeletonLhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
1141 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1144 pip->getSkeletonRhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
1145 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1149 auto side_data_ptr = boost::make_shared<SideData>();
1150 auto side_fe_ptr =
getSideFE(side_data_ptr);
1152 pip->getOpSkeletonRhsPipeline().push_back(
1154 pip->getOpSkeletonLhsPipeline().push_back(
1164 OpErrorSkel(boost::shared_ptr<FaceSideEle> side_fe_ptr,
1165 boost::shared_ptr<SideData> side_data_ptr,
1168 sideFEPtr(side_fe_ptr), sideDataPtr(side_data_ptr),
1169 errorSumPtr(error_sum_ptr), thError(th_error) {}
1175 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
1176 const auto in_the_loop =
1177 sideFEPtr->nInTheLoop;
1179 auto not_side = [](
auto s) {
1183 auto nb_gauss_pts = getGaussPts().size2();
1191 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
1192 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
1195 for (
auto &t_l : arr_t_l)
1197 for (
auto &t_vel : arr_t_vel)
1202 auto t_w = getFTensor0IntegrationWeight();
1203 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1204 e += t_w * getMeasure() *
1211 moab::Interface &moab =
1212 getNumeredEntFiniteElementPtr()->getBasicDataPtr()->moab;
1213 const void *tags_ptr[2];
1214 CHKERR moab.tag_get_by_ptr(thError, sideDataPtr->feSideHandle.data(), 2,
1216 for (
auto ff : {0, 1}) {
1217 *((
double *)tags_ptr[ff]) += e;
1219 CHKERR VecSetValue(errorSumPtr, 0, e, ADD_VALUES);
1226 boost::shared_ptr<FaceSideEle> sideFEPtr;
1227 boost::shared_ptr<SideData> sideDataPtr;
1238 MB_TAG_CREAT | MB_TAG_SPARSE,
1241 auto clear_tags = [&]() {
1250 auto evaluate_error = [&]() {
1252 auto skel_fe = boost::make_shared<BoundaryEle>(
mField);
1253 skel_fe->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
1254 auto side_data_ptr = boost::make_shared<SideData>();
1255 auto side_fe_ptr =
getSideFE(side_data_ptr);
1256 skel_fe->getOpPtrVector().push_back(
1257 new OpErrorSkel(side_fe_ptr, side_data_ptr, error_sum_ptr, th_error));
1260 skel_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1261 return fe_ptr->numeredEntFiniteElementPtr->
getBitRefLevel().test(
1266 simple->getSkeletonFEName(), skel_fe);
1271 auto assemble_and_sum = [](
auto vec) {
1279 auto propagate_error_to_parents = [&]() {
1283 auto fe_ptr = boost::make_shared<FEMethod>();
1284 fe_ptr->exeTestHook = [&](
FEMethod *fe_ptr) {
1285 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1289 fe_ptr->preProcessHook = []() {
return 0; };
1290 fe_ptr->postProcessHook = []() {
return 0; };
1291 fe_ptr->operatorHook = [&]() {
1294 auto fe_ent = fe_ptr->numeredEntFiniteElementPtr->getEnt();
1295 auto parent = fe_ptr->numeredEntFiniteElementPtr->getParentEnt();
1296 auto th_parent = fe_ptr->numeredEntFiniteElementPtr->getBasicDataPtr()
1297 ->th_RefParentHandle;
1300 CHKERR moab.tag_get_data(th_error, &fe_ent, 1, &error);
1303 [&](
auto fe_ent,
auto error) {
1306 CHKERR moab.tag_get_by_ptr(th_error, &fe_ent, 1,
1307 (
const void **)&e_ptr);
1311 CHKERR moab.tag_get_data(th_parent, &fe_ent, 1, &parent);
1312 if (parent != fe_ent && parent)
1313 CHKERR add_error(parent, *e_ptr);
1318 CHKERR add_error(parent, error);
1333 return std::make_tuple(assemble_and_sum(error_sum_ptr), th_error);
1351 DivergenceVol(boost::shared_ptr<VectorDouble> l_ptr,
1352 boost::shared_ptr<MatrixDouble> vel_ptr,
1359 const auto nb_dofs = data.
getIndices().size();
1361 const auto nb_gauss_pts = getGaussPts().size2();
1362 const auto t_w = getFTensor0IntegrationWeight();
1365 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
1367 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
1368 for (
int rr = 0; rr != nb_dofs; ++rr) {
1369 div += getMeasure() * t_w * t_l * (t_diff(
i) * t_vel(
i));
1376 CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
1382 boost::shared_ptr<VectorDouble> lPtr;
1383 boost::shared_ptr<MatrixDouble> velPtr;
1392 DivergenceSkeleton(boost::shared_ptr<SideData> side_data_ptr,
1393 boost::shared_ptr<FaceSideEle> side_fe_ptr,
1396 sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr), divVec(div_vec) {}
1403 &*base_mat.data().begin());
1406 auto not_side = [](
auto s) {
1411 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
1412 const auto in_the_loop =
1413 sideFEPtr->nInTheLoop;
1415 auto t_normal = getFTensor1Normal();
1416 const auto nb_gauss_pts = getGaussPts().size2();
1418 const auto nb_dofs = sideDataPtr->indicesRowSideMap[s0].size();
1420 auto t_base =
get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
1421 auto nb_row_base_functions = sideDataPtr->rowBaseSideMap[s0].size2();
1422 auto side_sense = sideDataPtr->senseMap[s0];
1423 auto opposite_s0 = not_side(s0);
1429 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
1430 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
1433 for (
auto &t_l : arr_t_l)
1435 for (
auto &t_vel : arr_t_vel)
1441 auto t_w = getFTensor0IntegrationWeight();
1442 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1446 const auto dot = (t_normal(
i) * t_vel(
i)) * side_sense;
1447 const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
1448 const auto l_upwind =
1449 arr_t_l[l_upwind_side];
1453 auto res = t_w * l_upwind * dot;
1457 for (; rr != nb_dofs; ++rr) {
1458 div += t_base * res;
1461 for (; rr < nb_row_base_functions; ++rr) {
1465 CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
1475 boost::shared_ptr<SideData> sideDataPtr;
1476 boost::shared_ptr<FaceSideEle> sideFEPtr;
1477 boost::shared_ptr<MatrixDouble> velPtr;
1481 auto vol_fe = boost::make_shared<DomainEle>(
mField);
1482 auto skel_fe = boost::make_shared<BoundaryEle>(
mField);
1484 vol_fe->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
1485 skel_fe->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
1490 auto l_ptr = boost::make_shared<VectorDouble>();
1491 auto vel_ptr = boost::make_shared<MatrixDouble>();
1492 auto side_data_ptr = boost::make_shared<SideData>();
1493 auto side_fe_ptr =
getSideFE(side_data_ptr);
1496 vol_fe->getOpPtrVector(), {potential_velocity_space, L2});
1497 vol_fe->getOpPtrVector().push_back(
1500 vol_fe->getOpPtrVector().push_back(
1501 new DivergenceVol(l_ptr, vel_ptr, div_vol_vec));
1503 skel_fe->getOpPtrVector().push_back(
1504 new DivergenceSkeleton(side_data_ptr, side_fe_ptr, div_skel_vec));
1507 auto dm =
simple->getDM();
1516 [](
double x,
double y,
double) {
return x - y; });
1518 [](
double x,
double y,
double) {
return x - y; });
1520 vol_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1521 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1524 skel_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1525 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1533 auto assemble_and_sum = [](
auto vec) {
1541 auto div_vol = assemble_and_sum(div_vol_vec);
1542 auto div_skel = assemble_and_sum(div_skel_vec);
1544 auto eps = std::abs((div_vol - div_skel) / (div_vol + div_skel));
1546 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Testing divergence volume: " << div_vol;
1547 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Testing divergence skeleton: " << div_skel
1548 <<
" relative difference: " <<
eps;
1550 constexpr double eps_err = 1e-6;
1553 "No consistency between skeleton integral and volume integral");
1571 auto post_proc = [&](
auto dm,
auto f_res,
auto out_name) {
1574 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
1578 auto l_vec = boost::make_shared<VectorDouble>();
1579 post_proc_fe->getOpPtrVector().push_back(
1582 post_proc_fe->getOpPtrVector().push_back(
1586 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1598 post_proc_fe->writeFile(out_name);
1602 constexpr double eps = 1e-4;
1605 opt->setRandomFields(
simple->getDM(), {{
"L", {-1, 1}}, {
"V", {-1, 1}}});
1606 auto dot_x = opt->setRandomFields(
simple->getDM(), {{
"L", {-1, 1}}});
1607 auto diff_x = opt->setRandomFields(
simple->getDM(), {{
"L", {-1, 1}}});
1609 auto test_domain_ops = [&](
auto fe_name,
auto lhs_pipeline,
1610 auto rhs_pipeline) {
1613 auto diff_res = opt->checkCentralFiniteDifference(
1614 simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
1617 if constexpr (
debug) {
1621 CHKERR post_proc(
simple->getDM(), diff_res,
"tangent_op_error.h5m");
1627 CHKERR VecNorm(diff_res, NORM_2, &fnorm);
1629 "Test consistency of tangent matrix %3.4e", fnorm);
1631 constexpr double err = 1e-9;
1634 "Norm of directional derivative too large err = %3.4e", fnorm);
1639 CHKERR test_domain_ops(
simple->getDomainFEName(), pip->getDomainLhsFE(),
1640 pip->getDomainRhsFE());
1641 CHKERR test_domain_ops(
simple->getSkeletonFEName(), pip->getSkeletonLhsFE(),
1642 pip->getSkeletonRhsFE());
1648 boost::function<
double(
double,
double,
double)> level_fun) {
1657 boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(
mField);
1658 boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(
mField);
1659 auto swap_fe = [&]() {
1660 lhs_fe.swap(pip->getDomainLhsFE());
1661 rhs_fe.swap(pip->getDomainRhsFE());
1665 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1666 pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1678 CHKERR prb_mng->removeDofsOnEntities(
"LEVELSET_POJECTION",
"L",
1680 auto test_bit_ele = [&](
FEMethod *fe_ptr) {
1681 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1684 pip->getDomainLhsFE()->exeTestHook = test_bit_ele;
1685 pip->getDomainRhsFE()->exeTestHook = test_bit_ele;
1688 pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
1690 pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
1691 pip->getOpDomainLhsPipeline().push_back(
new OpMassLL(
"L",
"L"));
1692 pip->getOpDomainRhsPipeline().push_back(
new OpSourceL(
"L", level_fun));
1696 auto ksp = pip->createKSP(sub_dm);
1697 CHKERR KSPSetDM(ksp, sub_dm);
1698 CHKERR KSPSetFromOptions(ksp);
1705 CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
1706 CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
1710 MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
1714 std::vector<Tag> tags{th_error};
1716 "PARALLEL=WRITE_PART", &fe_meshset, 1,
1717 &*tags.begin(), tags.size());
1720 auto post_proc = [&](
auto dm,
auto out_name,
auto th_error) {
1723 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
1724 post_proc_fe->setTagsToTransfer({th_error});
1725 post_proc_fe->exeTestHook = test_bit_ele;
1729 auto l_vec = boost::make_shared<VectorDouble>();
1730 auto l_grad_mat = boost::make_shared<MatrixDouble>();
1732 post_proc_fe->getOpPtrVector(), {potential_velocity_space, L2});
1733 post_proc_fe->getOpPtrVector().push_back(
1735 post_proc_fe->getOpPtrVector().push_back(
1738 post_proc_fe->getOpPtrVector().push_back(
1742 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1746 {{
"GradL", l_grad_mat}},
1754 post_proc_fe->writeFile(out_name);
1758 if constexpr (
debug)
1759 CHKERR post_proc(sub_dm,
"initial_level_set.h5m", th_error);
1767 boost::function<
double(
double,
double,
double)> vel_fun) {
1776 boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(
mField);
1777 boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(
mField);
1778 auto swap_fe = [&]() {
1779 lhs_fe.swap(pip->getDomainLhsFE());
1780 rhs_fe.swap(pip->getDomainRhsFE());
1784 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1785 pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1798 CHKERR prb_mng->removeDofsOnEntities(
"VELOCITY_PROJECTION",
"V",
1801 auto test_bit = [&](
FEMethod *fe_ptr) {
1802 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(0);
1804 pip->getDomainLhsFE()->exeTestHook = test_bit;
1805 pip->getDomainRhsFE()->exeTestHook = test_bit;
1808 pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
1810 pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
1812 pip->getOpDomainLhsPipeline().push_back(
new OpMassVV(
"V",
"V"));
1813 pip->getOpDomainRhsPipeline().push_back(
new OpSourceV(
"V", vel_fun));
1815 auto ksp = pip->createKSP(sub_dm);
1816 CHKERR KSPSetDM(ksp, sub_dm);
1817 CHKERR KSPSetFromOptions(ksp);
1824 CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
1825 CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
1828 auto post_proc = [&](
auto dm,
auto out_name) {
1831 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
1832 post_proc_fe->exeTestHook = test_bit;
1835 post_proc_fe->getOpPtrVector(), {potential_velocity_space});
1840 auto l_vec = boost::make_shared<VectorDouble>();
1841 auto potential_vec = boost::make_shared<VectorDouble>();
1842 auto velocity_mat = boost::make_shared<MatrixDouble>();
1844 post_proc_fe->getOpPtrVector().push_back(
1846 post_proc_fe->getOpPtrVector().push_back(
1850 post_proc_fe->getOpPtrVector().push_back(
1854 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1856 {{
"VelocityPotential", potential_vec}},
1858 {{
"Velocity", velocity_mat}},
1866 "3d case not implemented");
1871 post_proc_fe->writeFile(out_name);
1875 if constexpr (
debug)
1876 CHKERR post_proc(sub_dm,
"initial_velocity_potential.h5m");
1911 auto add_post_proc_fe = [&]() {
1912 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
1917 "Error", 1, MB_TYPE_DOUBLE, th_error, MB_TAG_CREAT | MB_TAG_SPARSE,
1919 post_proc_fe->setTagsToTransfer({th_error});
1921 post_proc_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1922 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1927 auto l_vec = boost::make_shared<VectorDouble>();
1928 auto vel_ptr = boost::make_shared<MatrixDouble>();
1931 post_proc_fe->getOpPtrVector(), {H1});
1932 post_proc_fe->getOpPtrVector().push_back(
1936 post_proc_fe->getOpPtrVector().push_back(
1940 post_proc_fe->getPostProcMesh(),
1942 post_proc_fe->getMapGaussPts(),
1951 return post_proc_fe;
1954 auto post_proc_fe = add_post_proc_fe();
1956 auto set_time_monitor = [&](
auto dm,
auto ts) {
1957 auto monitor_ptr = boost::make_shared<FEMethod>();
1959 monitor_ptr->preProcessHook = []() {
return 0; };
1960 monitor_ptr->operatorHook = []() {
return 0; };
1961 monitor_ptr->postProcessHook = [&]() {
1966 "Null pointer for post proc element");
1970 CHKERR post_proc_fe->writeFile(
1972 boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
".h5m");
1976 boost::shared_ptr<FEMethod> null;
1983 auto ts = pip->createTSIM(sub_dm);
1985 auto set_solution = [&](
auto ts) {
1994 auto monitor_pt = set_time_monitor(sub_dm, ts);
1995 CHKERR TSSetFromOptions(ts);
2003 auto ts_pre_step = [](TS ts) {
2010 MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
2012 auto get_norm = [&](
auto x) {
2014 CHKERR VecNorm(x, NORM_2, &nrm);
2018 auto set_solution = [&](
auto ts) {
2026 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
2027 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
2030 <<
"Problem " << prb_ptr->getName() <<
" solution vector norm "
2032 CHKERR TSSetSolution(ts, x);
2037 auto refine_and_project = [&](
auto ts) {
2055 ->removeDofsOnEntities(
"ADVECTION",
"L",
BitRefLevel().set(),
2061 auto ts_reset_theta = [&](
auto ts) {
2080 CHKERR refine_and_project(ts);
2081 CHKERR ts_reset_theta(ts);
2086 auto ts_post_step = [](TS ts) {
return 0; };
2088 CHKERR TSSetPreStep(ts, ts_pre_step);
2089 CHKERR TSSetPostStep(ts, ts_post_step);
2091 CHKERR TSSolve(ts, NULL);
2101 ParallelComm *pcomm =
2112 meshset_ptr->get_ptr(), 1);
2117 auto get_refined_elements_meshset = [&](
auto bit,
auto mask) {
2123 "get error handle");
2124 std::vector<double> errors(fe_ents.size());
2126 mField.
get_moab().tag_get_data(th_error, fe_ents, &*errors.begin()),
2128 auto it = std::max_element(errors.begin(), errors.end());
2130 MPI_Allreduce(&*it, &max, 1, MPI_DOUBLE, MPI_MAX,
mField.
get_comm());
2131 MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Max error: " << max;
2132 auto threshold = wp.getThreshold(max);
2134 std::vector<EntityHandle> fe_to_refine;
2135 fe_to_refine.reserve(fe_ents.size());
2137 auto fe_it = fe_ents.begin();
2138 auto error_it = errors.begin();
2139 for (
auto i = 0;
i != fe_ents.size(); ++
i) {
2140 if (*error_it > threshold) {
2141 fe_to_refine.push_back(*fe_it);
2148 ents.insert_list(fe_to_refine.begin(), fe_to_refine.end());
2150 ents,
nullptr,
NOISY);
2152 auto get_neighbours_by_bridge_vertices = [&](
auto &&ents) {
2156 moab::Interface::UNION);
2161 ents = get_neighbours_by_bridge_vertices(ents);
2167 "add entities to meshset");
2169 (proc_str +
"_fe_to_refine.vtk").c_str(),
"VTK",
"",
2170 meshset_ptr->get_ptr(), 1);
2178 auto refine_mesh = [&](
auto l,
auto &&fe_to_refine) {
2184 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2187 fe_to_refine = intersect(level_ents, fe_to_refine);
2189 level_ents = subtract(level_ents, fe_to_refine);
2192 Range fe_to_refine_children;
2193 bit_mng->updateRangeByChildren(fe_to_refine, fe_to_refine_children);
2195 fe_to_refine_children =
2196 fe_to_refine_children.subset_by_dimension(
SPACE_DIM);
2197 level_ents.merge(fe_to_refine_children);
2199 auto fix_neighbour_level = [&](
auto ll) {
2202 auto level_ll = level_ents;
2207 CHKERR skin.find_skin(0, level_ll,
false, skin_edges);
2210 for (
auto lll = 0; lll <= ll; ++lll) {
2211 CHKERR bit_mng->updateRangeByParent(skin_edges, skin_parents);
2215 for (
auto lll = 0; lll <= ll - 2; ++lll) {
2216 bad_bit[lll] =
true;
2219 Range skin_adj_ents;
2222 moab::Interface::UNION);
2225 skin_adj_ents = intersect(skin_adj_ents, level_ents);
2226 if (!skin_adj_ents.empty()) {
2227 level_ents = subtract(level_ents, skin_adj_ents);
2228 Range skin_adj_ents_children;
2229 bit_mng->updateRangeByChildren(skin_adj_ents, skin_adj_ents_children);
2230 level_ents.merge(skin_adj_ents_children);
2235 CHKERR fix_neighbour_level(
l);
2244 level_ents.subset_by_dimension(
SPACE_DIM), level_ents,
true);
2247 level_ents.subset_by_dimension(
SPACE_DIM), d,
false, level_ents,
2248 moab::Interface::UNION);
2260 CHKERR bit_mng->writeBitLevelByDim(
2262 (boost::lexical_cast<std::string>(
l) +
"_" + proc_str +
"_ref_mesh.vtk")
2271 auto set_skelton_bit = [&](
auto l) {
2281 Range level_edges_parents;
2282 CHKERR bit_mng->updateRangeByParent(level_edges, level_edges_parents);
2283 level_edges_parents =
2284 level_edges_parents.subset_by_dimension(
SPACE_DIM - 1);
2285 CHKERR bit_mng->filterEntitiesByRefLevel(
2289 auto parent_skeleton = intersect(level_edges, level_edges_parents);
2290 auto skeleton = subtract(level_edges, level_edges_parents);
2295 moab::Interface::UNION);
2303 CHKERR bit_mng->writeBitLevel(
2305 (boost::lexical_cast<std::string>(
l) +
"_" + proc_str +
"_skeleton.vtk")
2313 Range level0_current;
2317 Range level0_aggregate;
2323 start_mask[s] =
true;
2324 CHKERR bit_mng->lambdaBitRefLevel(
2339 CHKERR wp.setBits(*
this, 0);
2340 CHKERR wp.runCalcs(*
this, 0);
2342 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Process level: " <<
l;
2343 CHKERR refine_mesh(
l + 1, get_refined_elements_meshset(
2345 CHKERR set_skelton_bit(
l + 1);
2346 CHKERR wp.setAggregateBit(*
this,
l + 1);
2347 CHKERR wp.setBits(*
this,
l + 1);
2348 CHKERR wp.runCalcs(*
this,
l + 1);
2363 auto lhs_fe = boost::make_shared<DomainEle>(
mField);
2364 auto rhs_fe_prj = boost::make_shared<DomainEle>(
mField);
2365 auto rhs_fe_current = boost::make_shared<DomainEle>(
mField);
2367 lhs_fe->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
2368 rhs_fe_prj->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
2369 rhs_fe_current->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
2388 CHKERR bit_mng->updateRangeByParent(prj_ents, prj_ents);
2390 current_ents = subtract(current_ents, prj_ents);
2392 auto test_mesh_bit = [&](
FEMethod *fe_ptr) {
2393 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2396 auto test_prj_bit = [&](
FEMethod *fe_ptr) {
2397 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2400 auto test_current_bit = [&](
FEMethod *fe_ptr) {
2401 return current_ents.find(fe_ptr->getFEEntityHandle()) != current_ents.end();
2404 lhs_fe->exeTestHook = test_mesh_bit;
2405 rhs_fe_prj->exeTestHook = test_prj_bit;
2406 rhs_fe_current->exeTestHook = test_current_bit;
2410 CHKERR prb_mng->removeDofsOnEntities(
2411 "DG_PROJECTION",
"L",
BitRefLevel().set(), remove_mask,
nullptr, 0,
2415 lhs_fe->getOpPtrVector(), {potential_velocity_space, L2});
2416 lhs_fe->getOpPtrVector().push_back(
new OpMassLL(
"L",
"L"));
2418 auto l_vec = boost::make_shared<VectorDouble>();
2421 auto set_prj_from_child = [&](
auto rhs_fe_prj) {
2423 rhs_fe_prj->getOpPtrVector(), {potential_velocity_space, L2});
2426 rhs_fe_prj->getOpPtrVector().push_back(
2430 auto get_parent_this = [&]() {
2431 auto fe_parent_this = boost::make_shared<DomianParentEle>(
mField);
2432 fe_parent_this->getOpPtrVector().push_back(
2434 return fe_parent_this;
2439 auto get_parents_fe_ptr = [&](
auto this_fe_ptr) {
2440 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
2442 parents_elems_ptr_vec.emplace_back(
2443 boost::make_shared<DomianParentEle>(
mField));
2445 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
2451 return parents_elems_ptr_vec[0];
2454 auto this_fe_ptr = get_parent_this();
2455 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
2456 rhs_fe_prj->getOpPtrVector().push_back(
2463 auto set_prj_from_parent = [&](
auto rhs_fe_current) {
2466 auto get_parent_this = [&]() {
2467 auto fe_parent_this = boost::make_shared<DomianParentEle>(
mField);
2468 fe_parent_this->getOpPtrVector().push_back(
2470 return fe_parent_this;
2474 auto get_parents_fe_ptr = [&](
auto this_fe_ptr) {
2475 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
2477 parents_elems_ptr_vec.emplace_back(
2478 boost::make_shared<DomianParentEle>(
mField));
2480 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
2486 return parents_elems_ptr_vec[0];
2489 auto this_fe_ptr = get_parent_this();
2490 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
2500 rhs_fe_current->getOpPtrVector().push_back(reset_op_ptr);
2501 rhs_fe_current->getOpPtrVector().push_back(
2507 rhs_fe_current->getOpPtrVector().push_back(
new OpScalarFieldL(
"L", l_vec));
2510 set_prj_from_child(rhs_fe_prj);
2511 set_prj_from_parent(rhs_fe_current);
2513 boost::shared_ptr<FEMethod> null_fe;
2516 lhs_fe, null_fe, null_fe);
2520 rhs_fe_current, null_fe, null_fe);
2522 CHKERR KSPSetDM(ksp, sub_dm);
2524 CHKERR KSPSetDM(ksp, sub_dm);
2525 CHKERR KSPSetFromOptions(ksp);
2532 CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
2533 CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
2537 MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
2539 auto post_proc = [&](
auto dm,
auto out_name,
auto th_error) {
2542 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
2543 post_proc_fe->setTagsToTransfer({th_error});
2544 post_proc_fe->exeTestHook = test_mesh_bit;
2548 auto l_vec = boost::make_shared<VectorDouble>();
2549 auto l_grad_mat = boost::make_shared<MatrixDouble>();
2551 post_proc_fe->getOpPtrVector(), {potential_velocity_space, L2});
2552 post_proc_fe->getOpPtrVector().push_back(
2554 post_proc_fe->getOpPtrVector().push_back(
2557 post_proc_fe->getOpPtrVector().push_back(
2561 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2565 {{
"GradL", l_grad_mat}},
2573 post_proc_fe->writeFile(out_name);
2577 if constexpr (
debug)
2578 CHKERR post_proc(sub_dm,
"dg_projection.h5m", th_error);
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
BoundaryEle::UserDataOperator BoundaryEleOp
#define CATCH_ERRORS
Catch errors.
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
auto get_ntensor(T &base_mat)
ElementsAndOps< SPACE_DIM >::FaceSideEle FaceSideEle
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
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
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMSubDMSetUp_MoFEM(DM subdm)
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.
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
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
constexpr int skeleton_bit
skeleton elemets bit
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomianParentEle DomianParentEle
FTensor::Index< 'j', SPACE_DIM > j
constexpr int current_bit
dofs bit used to do calculations
constexpr FieldSpace potential_velocity_space
FTensor::Index< 'k', SPACE_DIM > k
FTensor::Index< 'i', SPACE_DIM > i
constexpr int SPACE_DIM
[Define dimension]
FaceSideEle::UserDataOperator FaceSideEleOp
constexpr int aggregate_projection_bit
all bits for projection problem
constexpr IntegrationType I
DomainEle::UserDataOperator DomainEleOp
LevelSet * level_set_raw_ptr
constexpr int aggregate_bit
all bits for advection problem
constexpr int projection_bit
constexpr size_t potential_velocity_field_dim
FTensor::Index< 'l', 3 > l
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
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.
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobina in TS solver.
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
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.
constexpr auto make_array(Arg &&...arg)
Create Array.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
auto getProblemPtr(DM dm)
get problem pointer from DM
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
constexpr double t
plate stiffness
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
OpLhsDomain(const std::string field_name, boost::shared_ptr< MatrixDouble > vel_ptr)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > velPtr
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
boost::shared_ptr< SideData > sideDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpLhsSkeleton(boost::shared_ptr< SideData > side_data_ptr, boost::shared_ptr< FaceSideEle > side_fe_ptr)
boost::shared_ptr< VectorDouble > lPtr
MoFEMErrorCode iNtegrate(EntData &data)
Class dedicated to integrate operator.
boost::shared_ptr< MatrixDouble > velPtr
OpRhsDomain(const std::string field_name, boost::shared_ptr< VectorDouble > l_ptr, boost::shared_ptr< VectorDouble > l_dot_ptr, boost::shared_ptr< MatrixDouble > vel_ptr)
boost::shared_ptr< VectorDouble > lDotPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< SideData > sideDataPtr
OpRhsSkeleton(boost::shared_ptr< SideData > side_data_ptr, boost::shared_ptr< FaceSideEle > side_fe_ptr)
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
data structure carrying information on skeleton on both sides.
std::array< VectorInt, 2 > indicesColSideMap
indices on columns for left hand-side
std::array< EntityHandle, 2 > feSideHandle
int currentFESide
current side counter
std::array< MatrixDouble, 2 > colBaseSideMap
std::array< MatrixDouble, 2 > rowBaseSideMap
std::array< VectorInt, 2 > indicesRowSideMap
indices on rows for left hand-side
std::array< int, 2 > senseMap
Use peculated errors on all levels while mesh projection.
WrapperClassErrorProjection(boost::shared_ptr< double > max_ptr)
MoFEMErrorCode runCalcs(LevelSet &level_set, int l)
MoFEMErrorCode setBits(LevelSet &level_set, int l)
double getThreshold(const double max)
MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l)
boost::shared_ptr< double > maxPtr
Wrapper executing stages while mesh refinement.
virtual double getThreshold(const double max)=0
virtual MoFEMErrorCode runCalcs(LevelSet &level_set, int l)=0
virtual MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l)=0
virtual MoFEMErrorCode setBits(LevelSet &level_set, int l)=0
Used to execute inital mesh approximation while mesh refinement.
double getThreshold(const double max)
MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l)
boost::shared_ptr< double > maxPtr
MoFEMErrorCode setBits(LevelSet &level_set, int l)
WrapperClassInitalSolution(boost::shared_ptr< double > max_ptr)
MoFEMErrorCode runCalcs(LevelSet &level_set, int l)
MoFEM::Interface & mField
integrate skeleton operators on khs
MoFEMErrorCode solveAdvection()
solve advection problem
boost::shared_ptr< double > maxPtr
std::array< MatrixDouble, 2 > MatSideArray
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< 1, 1 > OpMassLL
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpSource< potential_velocity_field_dim, potential_velocity_field_dim > OpSourceV
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< potential_velocity_field_dim, potential_velocity_field_dim > OpMassVV
MoFEMErrorCode readMesh()
read mesh
MoFEMErrorCode initialiseFieldLevelSet(boost::function< double(double, double, double)> level_fun=get_level_set)
initialise field set
MoFEMErrorCode setupProblem()
create fields, and set approximation order
boost::shared_ptr< FaceSideEle > getSideFE(boost::shared_ptr< SideData > side_data_ptr)
create side element to assemble data from sides
MoFEMErrorCode pushOpDomain()
push operators to integrate operators on domain
MoFEMErrorCode runProblem()
LevelSet(MoFEM::Interface &m_field)
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< 1 > OpScalarFieldL
static double get_velocity_potential(double x, double y, double z)
advection velocity field
std::array< VectorDouble, 2 > VecSideArray
MoFEMErrorCode pushOpSkeleton()
push operator to integrate on skeleton
MoFEMErrorCode testOp()
test consistency between tangent matrix and the right hand side vectors
MoFEMErrorCode dgProjection(const int prj_bit=projection_bit)
dg level set projection
MoFEMErrorCode initialiseFieldVelocity(boost::function< double(double, double, double)> vel_fun=get_velocity_potential< SPACE_DIM >)
initialise potential velocity field
static double get_level_set(const double x, const double y, const double z)
inital level set, i.e. advected filed
MoFEMErrorCode refineMesh(WrapperClass &&wp)
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpSource< 1, 1 > OpSourceL
ForcesAndSourcesCore::UserDataOperator * getZeroLevelVelOp(boost::shared_ptr< MatrixDouble > vel_ptr)
Get operator calculating velocity on coarse mesh.
std::tuple< double, Tag > evaluateError()
evaluate error
MoFEMErrorCode testSideFE()
test integration side elements
Add operators pushing bases from local to physical configuration.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
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
default operator for TRI element
@ 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.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Calculate curl of vector field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
intrusive_ptr for managing petsc objects
Interface for Time Stepping (TS) solver.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.