14static char help[] =
"...\n\n";
34 IntegrationType::GAUSS;
89 std::array<VectorInt, 2>
91 std::array<VectorInt, 2>
119 template <
int FE_DIM>
130 static double get_level_set(
const double x,
const double y,
const double z);
166 ForcesAndSourcesCore::UserDataOperator *
175 boost::shared_ptr<FaceSideEle>
176 getSideFE(boost::shared_ptr<SideData> side_data_ptr);
209 boost::function<
double(
double,
double,
double)> level_fun =
219 boost::function<
double(
double,
double,
double)> vel_fun =
220 get_velocity_potential<FE_DIM>);
274 simple->getBitRefLevel() =
295 ->synchroniseEntities(level);
304 return 0.05 * (*maxPtr);
329 ->synchroniseEntities(level);
356 G>::OpSource<1, DIM1 * DIM2>;
360 G>::OpSource<potential_velocity_field_dim, potential_velocity_field_dim>;
362 G>::OpBaseTimesVector<1, DIM1 * DIM2, 1>;
374double LevelSet::get_velocity_potential<2>(
double x,
double y,
double z) {
375 return (x * x - 0.25) * (y * y - 0.25);
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;
391 if constexpr (
debug) {
398 maxPtr = boost::make_shared<double>(0);
402 simple->getBitRefLevel() =
415 boost::shared_ptr<MatrixDouble> l_ptr,
416 boost::shared_ptr<MatrixDouble> l_dot_ptr,
417 boost::shared_ptr<MatrixDouble> vel_ptr);
421 boost::shared_ptr<MatrixDouble>
lPtr;
428 boost::shared_ptr<MatrixDouble> vel_ptr);
438 boost::shared_ptr<FaceSideEle> side_fe_ptr);
444 boost::shared_ptr<FaceSideEle>
452 boost::shared_ptr<FaceSideEle> side_fe_ptr);
458 boost::shared_ptr<FaceSideEle>
464int main(
int argc,
char *argv[]) {
473 moab::Core moab_core;
474 moab::Interface &moab = moab_core;
481 DMType dm_name =
"DMMOFEM";
485 auto core_log = logging::core::get();
511 simple->getAddSkeletonFE() =
true;
512 simple->getAddBoundaryFE() =
true;
519 auto set_problem_bit = [&]() {
524 start_mask[s] =
true;
544 if constexpr (
debug) {
546 CHKERR bit_mng->writeBitLevelByDim(
548 (proc_str +
"level_base.vtk").c_str(),
"VTK",
"");
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) {}
591 boost::shared_ptr<MatrixDouble> vel_ptr)
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) {}
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) {}
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();
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);
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);
636 auto &nf = this->locF;
637 auto t_nf = getFTensor2FromPtr<DIM1, DIM2>(&*nf.begin());
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);
646 for (; rr < nb_base_func; ++rr) {
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();
664 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
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();
675 auto &mat = this->locMat;
678 for (; rr != nb_row_dofs; ++rr) {
680 auto t_mat = getFTensor2FromPtr<DIM1, DIM2>(&mat(rr *
DIM1, 0));
681 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
683 (beta * t_row_base - alpha * (t_row_diff_base(
i) * t_vel(
i))) *
691 for (; rr < nb_base_func; ++rr) {
708 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
709 const auto in_the_loop =
710 sideFEPtr->nInTheLoop;
712 auto not_side = [](
auto s) {
718 &*base_mat.data().begin());
721 if (in_the_loop > 0) {
724 auto t_normal = getFTensor1Normal();
725 const auto nb_gauss_pts = getGaussPts().size2();
730 const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
734 resSkelton.resize(nb_rows,
false);
738 const auto opposite_s0 = not_side(s0);
739 const auto sense_row = sideDataPtr->senseMap[s0];
741 const auto opposite_sense_row = sideDataPtr->senseMap[opposite_s0];
742 if (sense_row * opposite_sense_row > 0)
744 "Should be opposite sign");
748 const auto nb_row_base_functions =
749 sideDataPtr->rowBaseSideMap[s0].size2();
751 auto t_w = getFTensor0IntegrationWeight();
753 getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[
LEFT_SIDE]),
754 getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[
RIGHT_SIDE]));
756 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
757 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
760 for (
auto &t_l : arr_t_l)
762 for (
auto &t_vel : arr_t_vel)
767 if (nb_gauss_pts != sideDataPtr->rowBaseSideMap[s0].size1())
769 "Inconsistent number of DOFs");
772 auto t_row_base =
get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
773 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
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);
784 auto t_res_skeleton =
785 getFTensor2FromPtr<DIM1, DIM2>(&*resSkelton.data().begin());
787 for (; rr != nb_rows; ++rr) {
788 t_res_skeleton(
I,
J) += t_row_base * t_res(
I,
J);
792 for (; rr < nb_row_base_functions; ++rr) {
797 CHKERR ::VecSetValues(getTSf(),
798 sideDataPtr->indicesRowSideMap[s0].size(),
799 &*sideDataPtr->indicesRowSideMap[s0].begin(),
800 &*resSkelton.begin(), ADD_VALUES);
814 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
815 const auto in_the_loop =
816 sideFEPtr->nInTheLoop;
818 auto not_side = [](
auto s) {
824 &*base_mat.data().begin());
827 if (in_the_loop > 0) {
830 auto t_normal = getFTensor1Normal();
831 const auto nb_gauss_pts = getGaussPts().size2();
836 const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
841 const auto opposite_s0 = not_side(s0);
842 const auto sense_row = sideDataPtr->senseMap[s0];
845 const auto nb_row_base_functions =
846 sideDataPtr->rowBaseSideMap[s0].size2();
851 const auto nb_cols = sideDataPtr->indicesColSideMap[s1].size();
852 const auto sense_col = sideDataPtr->senseMap[s1];
855 matSkeleton.resize(nb_rows, nb_cols,
false);
858 auto t_w = getFTensor0IntegrationWeight();
860 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
861 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
864 for (
auto &t_vel : arr_t_vel)
868 auto t_row_base =
get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
869 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
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;
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);
888 get_ntensor(sideDataPtr->colBaseSideMap[s1], gg, 0);
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;
896 for (
size_t cc = 0; cc != nb_cols; ++cc) {
897 t_mat_skeleton(
I,
J) += t_res_row(
I,
J) * t_col_base;
903 for (; rr < nb_row_base_functions; ++rr) {
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);
921ForcesAndSourcesCore::UserDataOperator *
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(
930 return parent_fe_ptr;
933 auto get_parents_vel_fe_ptr = [&](
auto this_fe_ptr) {
934 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
936 parents_elems_ptr_vec.emplace_back(
937 boost::make_shared<DomianParentEle>(
mField));
939 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
944 return parents_elems_ptr_vec[0];
947 auto this_fe_ptr = get_parent_vel_this();
948 auto parent_fe_ptr = get_parents_vel_fe_ptr(this_fe_ptr);
960 pip->getOpDomainRhsPipeline().clear();
962 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
963 pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
965 pip->getDomainLhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
966 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
969 pip->getDomainRhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
970 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
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>();
981 pip->getOpDomainRhsPipeline().push_back(
983 pip->getOpDomainRhsPipeline().push_back(
985 pip->getOpDomainRhsPipeline().push_back(
991 pip->getOpDomainLhsPipeline().push_back(
new OpLhsDomain(
"L", vel_ptr));
996boost::shared_ptr<FaceSideEle>
1001 auto l_ptr = boost::make_shared<MatrixDouble>();
1002 auto vel_ptr = boost::make_shared<MatrixDouble>();
1005 OpSideData(boost::shared_ptr<SideData> side_data_ptr)
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;
1019 if ((CN::Dimension(row_type) ==
FE_DIM) &&
1020 (CN::Dimension(col_type) ==
FE_DIM)) {
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;
1031 const auto nb_in_loop = getFEMethod()->nInTheLoop;
1032 if (nb_in_loop == 0)
1033 for (
auto s : {0, 1})
1036 sideDataPtr->currentFESide = nb_in_loop;
1037 sideDataPtr->senseMap[nb_in_loop] = getSkeletonSense();
1047 boost::shared_ptr<SideData> sideDataPtr;
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;
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;
1081 if ((sideDataPtr->lVec)[nb_in_loop].size2() !=
1082 (sideDataPtr->velMat)[nb_in_loop].size2())
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)
1089 "Wrong size of velocity vector size = %d",
1090 (sideDataPtr->velMat)[nb_in_loop].size1());
1094 (sideDataPtr->lVec)[1] = sideDataPtr->lVec[0];
1095 (sideDataPtr->velMat)[1] = (sideDataPtr->velMat)[0];
1098 if (sideDataPtr->rowBaseSideMap[0].size1() !=
1099 sideDataPtr->rowBaseSideMap[1].size1()) {
1101 "Wrong number of integration pt %d != %d",
1102 sideDataPtr->rowBaseSideMap[0].size1(),
1103 sideDataPtr->rowBaseSideMap[1].size1());
1105 if (sideDataPtr->colBaseSideMap[0].size1() !=
1106 sideDataPtr->colBaseSideMap[1].size1()) {
1108 "Wrong number of integration pt");
1121 boost::shared_ptr<SideData> sideDataPtr;
1122 boost::shared_ptr<MatrixDouble> lPtr;
1123 boost::shared_ptr<MatrixDouble> velPtr;
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;
1138 auto get_parents_fe_ptr = [&](
auto this_fe_ptr) {
1139 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
1141 parents_elems_ptr_vec.emplace_back(
1142 boost::make_shared<DomianParentEle>(
mField));
1144 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
1149 return parents_elems_ptr_vec[0];
1154 auto get_side_fe_ptr = [&]() {
1155 auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
1157 auto this_fe_ptr = get_parent_this();
1158 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
1160 side_fe_ptr->getOpPtrVector().push_back(
new OpSideData(side_data_ptr));
1162 side_fe_ptr->getOpPtrVector().push_back(
1170 return get_side_fe_ptr();
1178 pip->getOpSkeletonRhsPipeline().clear();
1180 pip->setSkeletonLhsIntegrationRule([](
int,
int,
int o) {
return 18; });
1181 pip->setSkeletonRhsIntegrationRule([](
int,
int,
int o) {
return 18; });
1183 pip->getSkeletonLhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
1184 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1187 pip->getSkeletonRhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
1188 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1192 auto side_data_ptr = boost::make_shared<SideData>();
1193 auto side_fe_ptr =
getSideFE(side_data_ptr);
1195 pip->getOpSkeletonRhsPipeline().push_back(
1197 pip->getOpSkeletonLhsPipeline().push_back(
1207 OpErrorSkel(boost::shared_ptr<FaceSideEle> side_fe_ptr,
1208 boost::shared_ptr<SideData> side_data_ptr,
1211 sideFEPtr(side_fe_ptr), sideDataPtr(side_data_ptr),
1212 errorSumPtr(error_sum_ptr), thError(th_error) {}
1218 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
1219 const auto in_the_loop =
1220 sideFEPtr->nInTheLoop;
1222 auto not_side = [](
auto s) {
1226 auto nb_gauss_pts = getGaussPts().size2();
1231 getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[
LEFT_SIDE]),
1232 getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[
RIGHT_SIDE]));
1234 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
1235 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
1238 for (
auto &t_l : arr_t_l)
1240 for (
auto &t_vel : arr_t_vel)
1245 auto t_w = getFTensor0IntegrationWeight();
1246 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1249 e += t_w * getMeasure() * t_diff(
I,
J) * t_diff(
I,
J);
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,
1260 for (
auto ff : {0, 1}) {
1261 *((
double *)tags_ptr[ff]) += e;
1263 CHKERR VecSetValue(errorSumPtr, 0, e, ADD_VALUES);
1270 boost::shared_ptr<FaceSideEle> sideFEPtr;
1271 boost::shared_ptr<SideData> sideDataPtr;
1282 MB_TAG_CREAT | MB_TAG_SPARSE,
1285 auto clear_tags = [&]() {
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));
1304 skel_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1305 return fe_ptr->numeredEntFiniteElementPtr->
getBitRefLevel().test(
1310 simple->getSkeletonFEName(), skel_fe);
1315 auto assemble_and_sum = [](
auto vec) {
1323 auto propagate_error_to_parents = [&]() {
1327 auto fe_ptr = boost::make_shared<FEMethod>();
1328 fe_ptr->exeTestHook = [&](
FEMethod *fe_ptr) {
1329 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1333 fe_ptr->preProcessHook = []() {
return 0; };
1334 fe_ptr->postProcessHook = []() {
return 0; };
1335 fe_ptr->operatorHook = [&]() {
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;
1344 CHKERR moab.tag_get_data(th_error, &fe_ent, 1, &error);
1347 [&](
auto fe_ent,
auto error) {
1350 CHKERR moab.tag_get_by_ptr(th_error, &fe_ent, 1,
1351 (
const void **)&e_ptr);
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);
1362 CHKERR add_error(parent, error);
1377 return std::make_tuple(assemble_and_sum(error_sum_ptr), th_error);
1395 DivergenceVol(boost::shared_ptr<MatrixDouble> l_ptr,
1396 boost::shared_ptr<MatrixDouble> vel_ptr,
1403 const auto nb_dofs = data.
getIndices().size();
1405 const auto nb_gauss_pts = getGaussPts().size2();
1406 const auto t_w = getFTensor0IntegrationWeight();
1408 auto t_l = getFTensor2FromMat<DIM1, DIM2>(*lPtr);
1409 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
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));
1420 CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
1426 boost::shared_ptr<MatrixDouble> lPtr;
1427 boost::shared_ptr<MatrixDouble> velPtr;
1436 DivergenceSkeleton(boost::shared_ptr<SideData> side_data_ptr,
1437 boost::shared_ptr<FaceSideEle> side_fe_ptr,
1440 sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr), divVec(div_vec) {}
1447 &*base_mat.data().begin());
1450 auto not_side = [](
auto s) {
1455 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
1456 const auto in_the_loop =
1457 sideFEPtr->nInTheLoop;
1459 auto t_normal = getFTensor1Normal();
1460 const auto nb_gauss_pts = getGaussPts().size2();
1462 const auto nb_dofs = sideDataPtr->indicesRowSideMap[s0].size();
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);
1470 getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[
LEFT_SIDE]),
1471 getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[
RIGHT_SIDE]));
1473 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
1474 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
1477 for (
auto &t_l : arr_t_l)
1479 for (
auto &t_vel : arr_t_vel)
1485 auto t_w = getFTensor0IntegrationWeight();
1486 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
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];
1497 auto res = t_w * l_upwind(
I,
I) * dot;
1501 for (; rr != nb_dofs; ++rr) {
1502 div += t_base * res;
1505 for (; rr < nb_row_base_functions; ++rr) {
1509 CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
1519 boost::shared_ptr<SideData> sideDataPtr;
1520 boost::shared_ptr<FaceSideEle> sideFEPtr;
1521 boost::shared_ptr<MatrixDouble> velPtr;
1525 auto vol_fe = boost::make_shared<DomainEle>(
mField);
1526 auto skel_fe = boost::make_shared<BoundaryEle>(
mField);
1528 vol_fe->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
1529 skel_fe->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
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);
1541 vol_fe->getOpPtrVector().push_back(
1544 vol_fe->getOpPtrVector().push_back(
1545 new DivergenceVol(l_ptr, vel_ptr, div_vol_vec));
1547 skel_fe->getOpPtrVector().push_back(
1548 new DivergenceSkeleton(side_data_ptr, side_fe_ptr, div_skel_vec));
1551 auto dm =
simple->getDM();
1560 [](
double x,
double y,
double) {
return x - y; });
1562 [](
double x,
double y,
double) {
return x - y; });
1564 vol_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1565 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1568 skel_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1569 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1577 auto assemble_and_sum = [](
auto vec) {
1585 auto div_vol = assemble_and_sum(div_vol_vec);
1586 auto div_skel = assemble_and_sum(div_skel_vec);
1588 auto eps = std::abs((div_vol - div_skel) / (div_vol + div_skel));
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;
1594 constexpr double eps_err = 1e-6;
1597 "No consistency between skeleton integral and volume integral");
1615 auto post_proc = [&](
auto dm,
auto f_res,
auto out_name) {
1618 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
1620 if constexpr (
DIM1 == 1 &&
DIM2 == 1) {
1623 auto l_vec = boost::make_shared<VectorDouble>();
1624 post_proc_fe->getOpPtrVector().push_back(
1627 post_proc_fe->getOpPtrVector().push_back(
1631 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1644 post_proc_fe->writeFile(out_name);
1648 constexpr double eps = 1e-4;
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}}});
1655 auto test_domain_ops = [&](
auto fe_name,
auto lhs_pipeline,
1656 auto rhs_pipeline) {
1659 auto diff_res = opt->checkCentralFiniteDifference(
1660 simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
1663 if constexpr (
debug) {
1667 CHKERR post_proc(
simple->getDM(), diff_res,
"tangent_op_error.h5m");
1673 CHKERR VecNorm(diff_res, NORM_2, &fnorm);
1675 "Test consistency of tangent matrix %3.4e", fnorm);
1677 constexpr double err = 1e-9;
1680 "Norm of directional derivative too large err = %3.4e", fnorm);
1685 CHKERR test_domain_ops(
simple->getDomainFEName(), pip->getDomainLhsFE(),
1686 pip->getDomainRhsFE());
1687 CHKERR test_domain_ops(
simple->getSkeletonFEName(), pip->getSkeletonLhsFE(),
1688 pip->getSkeletonRhsFE());
1694 boost::function<
double(
double,
double,
double)> level_fun) {
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());
1711 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1712 pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1724 CHKERR prb_mng->removeDofsOnEntities(
"LEVELSET_POJECTION",
"L",
1726 auto test_bit_ele = [&](
FEMethod *fe_ptr) {
1727 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1730 pip->getDomainLhsFE()->exeTestHook = test_bit_ele;
1731 pip->getDomainRhsFE()->exeTestHook = test_bit_ele;
1737 pip->getOpDomainLhsPipeline().push_back(
new OpMassLL(
"L",
"L"));
1738 pip->getOpDomainRhsPipeline().push_back(
new OpSourceL(
"L", level_fun));
1742 auto ksp = pip->createKSP(sub_dm);
1743 CHKERR KSPSetDM(ksp, sub_dm);
1744 CHKERR KSPSetFromOptions(ksp);
1751 CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
1752 CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
1756 MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
1760 std::vector<Tag> tags{th_error};
1762 "PARALLEL=WRITE_PART", &fe_meshset, 1,
1763 &*tags.begin(), tags.size());
1766 auto post_proc = [&](
auto dm,
auto out_name,
auto th_error) {
1769 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
1770 post_proc_fe->setTagsToTransfer({th_error});
1771 post_proc_fe->exeTestHook = test_bit_ele;
1773 if constexpr (
DIM1 == 1 &&
DIM2 == 1) {
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(
1782 post_proc_fe->getOpPtrVector().push_back(
1785 post_proc_fe->getOpPtrVector().push_back(
1789 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1793 {{
"GradL", l_grad_mat}},
1802 post_proc_fe->writeFile(out_name);
1806 if constexpr (
debug)
1807 CHKERR post_proc(sub_dm,
"initial_level_set.h5m", th_error);
1815 boost::function<
double(
double,
double,
double)> vel_fun) {
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());
1832 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1833 pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1846 CHKERR prb_mng->removeDofsOnEntities(
"VELOCITY_PROJECTION",
"V",
1849 auto test_bit = [&](
FEMethod *fe_ptr) {
1850 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(0);
1852 pip->getDomainLhsFE()->exeTestHook = test_bit;
1853 pip->getDomainRhsFE()->exeTestHook = test_bit;
1856 {potential_velocity_space});
1858 {potential_velocity_space});
1860 pip->getOpDomainLhsPipeline().push_back(
new OpMassVV(
"V",
"V"));
1861 pip->getOpDomainRhsPipeline().push_back(
new OpSourceV(
"V", vel_fun));
1863 auto ksp = pip->createKSP(sub_dm);
1864 CHKERR KSPSetDM(ksp, sub_dm);
1865 CHKERR KSPSetFromOptions(ksp);
1872 CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
1873 CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
1876 auto post_proc = [&](
auto dm,
auto out_name) {
1879 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
1880 post_proc_fe->exeTestHook = test_bit;
1882 if constexpr (
FE_DIM == 2) {
1885 post_proc_fe->getOpPtrVector(), {potential_velocity_space});
1889 auto potential_vec = boost::make_shared<VectorDouble>();
1890 auto velocity_mat = boost::make_shared<MatrixDouble>();
1892 post_proc_fe->getOpPtrVector().push_back(
1894 post_proc_fe->getOpPtrVector().push_back(
1898 post_proc_fe->getOpPtrVector().push_back(
1902 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1904 {{
"VelocityPotential", potential_vec}},
1906 {{
"Velocity", velocity_mat}},
1914 "3d case not implemented");
1919 post_proc_fe->writeFile(out_name);
1923 if constexpr (
debug)
1924 CHKERR post_proc(sub_dm,
"initial_velocity_potential.h5m");
1959 auto add_post_proc_fe = [&]() {
1960 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
1965 "Error", 1, MB_TYPE_DOUBLE, th_error, MB_TAG_CREAT | MB_TAG_SPARSE,
1967 post_proc_fe->setTagsToTransfer({th_error});
1969 post_proc_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1970 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1976 auto vel_ptr = boost::make_shared<MatrixDouble>();
1979 post_proc_fe->getOpPtrVector(), {L2});
1982 if constexpr (
DIM1 == 1 &&
DIM2 == 1) {
1983 auto l_vec = boost::make_shared<VectorDouble>();
1984 post_proc_fe->getOpPtrVector().push_back(
1986 post_proc_fe->getOpPtrVector().push_back(
1990 post_proc_fe->getPostProcMesh(),
1992 post_proc_fe->getMapGaussPts(),
2002 return post_proc_fe;
2005 auto post_proc_fe = add_post_proc_fe();
2007 auto set_time_monitor = [&](
auto dm,
auto ts) {
2008 auto monitor_ptr = boost::make_shared<FEMethod>();
2010 monitor_ptr->preProcessHook = []() {
return 0; };
2011 monitor_ptr->operatorHook = []() {
return 0; };
2012 monitor_ptr->postProcessHook = [&]() {
2017 "Null pointer for post proc element");
2021 CHKERR post_proc_fe->writeFile(
2023 boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
".h5m");
2027 boost::shared_ptr<FEMethod> null;
2034 auto ts = pip->createTSIM(sub_dm);
2036 auto set_solution = [&](
auto ts) {
2045 auto monitor_pt = set_time_monitor(sub_dm, ts);
2046 CHKERR TSSetFromOptions(ts);
2054 auto ts_pre_step = [](TS ts) {
2061 MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
2063 auto get_norm = [&](
auto x) {
2065 CHKERR VecNorm(x, NORM_2, &nrm);
2069 auto set_solution = [&](
auto ts) {
2077 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
2078 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
2081 <<
"Problem " << prb_ptr->getName() <<
" solution vector norm "
2083 CHKERR TSSetSolution(ts, x);
2088 auto refine_and_project = [&](
auto ts) {
2106 ->removeDofsOnEntities(
"ADVECTION",
"L",
BitRefLevel().set(),
2112 auto ts_reset_theta = [&](
auto ts) {
2131 CHKERR refine_and_project(ts);
2132 CHKERR ts_reset_theta(ts);
2137 auto ts_post_step = [](TS ts) {
return 0; };
2139 CHKERR TSSetPreStep(ts, ts_pre_step);
2140 CHKERR TSSetPostStep(ts, ts_post_step);
2142 CHKERR TSSolve(ts, NULL);
2152 ParallelComm *pcomm =
2163 meshset_ptr->get_ptr(), 1);
2168 auto get_refined_elements_meshset = [&](
auto bit,
auto mask) {
2170 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
bit, mask,
FE_DIM, fe_ents);
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()),
2179 auto it = std::max_element(errors.begin(), errors.end());
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);
2185 std::vector<EntityHandle> fe_to_refine;
2186 fe_to_refine.reserve(fe_ents.size());
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);
2199 ents.insert_list(fe_to_refine.begin(), fe_to_refine.end());
2201 ents,
nullptr,
NOISY);
2203 auto get_neighbours_by_bridge_vertices = [&](
auto &&ents) {
2207 moab::Interface::UNION);
2212 ents = get_neighbours_by_bridge_vertices(ents);
2218 "add entities to meshset");
2220 (proc_str +
"_fe_to_refine.vtk").c_str(),
"VTK",
"",
2221 meshset_ptr->get_ptr(), 1);
2229 auto refine_mesh = [&](
auto l,
auto &&fe_to_refine) {
2235 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2238 fe_to_refine = intersect(level_ents, fe_to_refine);
2240 level_ents = subtract(level_ents, fe_to_refine);
2243 Range fe_to_refine_children;
2244 bit_mng->updateRangeByChildren(fe_to_refine, fe_to_refine_children);
2246 fe_to_refine_children = fe_to_refine_children.subset_by_dimension(
FE_DIM);
2247 level_ents.merge(fe_to_refine_children);
2249 auto fix_neighbour_level = [&](
auto ll) {
2252 auto level_ll = level_ents;
2257 CHKERR skin.find_skin(0, level_ll,
false, skin_edges);
2260 for (
auto lll = 0; lll <= ll; ++lll) {
2261 CHKERR bit_mng->updateRangeByParent(skin_edges, skin_parents);
2265 for (
auto lll = 0; lll <= ll - 2; ++lll) {
2266 bad_bit[lll] =
true;
2269 Range skin_adj_ents;
2271 skin_parents,
FE_DIM,
false, skin_adj_ents, moab::Interface::UNION);
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);
2284 CHKERR fix_neighbour_level(
l);
2290 for (
auto d = 0; d !=
FE_DIM; ++d) {
2293 level_ents.subset_by_dimension(
FE_DIM), level_ents,
true);
2296 level_ents.subset_by_dimension(
FE_DIM), d,
false, level_ents,
2297 moab::Interface::UNION);
2309 CHKERR bit_mng->writeBitLevelByDim(
2311 (boost::lexical_cast<std::string>(
l) +
"_" + proc_str +
"_ref_mesh.vtk")
2320 auto set_skelton_bit = [&](
auto l) {
2325 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
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(
2336 auto parent_skeleton = intersect(level_edges, level_edges_parents);
2337 auto skeleton = subtract(level_edges, level_edges_parents);
2342 moab::Interface::UNION);
2350 CHKERR bit_mng->writeBitLevel(
2352 (boost::lexical_cast<std::string>(
l) +
"_" + proc_str +
"_skeleton.vtk")
2360 Range level0_current;
2364 Range level0_aggregate;
2370 start_mask[s] =
true;
2371 CHKERR bit_mng->lambdaBitRefLevel(
2386 CHKERR wp.setBits(*
this, 0);
2387 CHKERR wp.runCalcs(*
this, 0);
2389 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Process level: " <<
l;
2390 CHKERR refine_mesh(
l + 1, get_refined_elements_meshset(
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);
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);
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; };
2431 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2434 CHKERR bit_mng->updateRangeByParent(prj_ents, prj_ents);
2436 current_ents = subtract(
2437 current_ents, prj_ents);
2439 auto test_mesh_bit = [&](
FEMethod *fe_ptr) {
2440 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2443 auto test_prj_bit = [&](
FEMethod *fe_ptr) {
2444 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2447 auto test_current_bit = [&](
FEMethod *fe_ptr) {
2448 return current_ents.find(fe_ptr->getFEEntityHandle()) != current_ents.end();
2451 lhs_fe->exeTestHook =
2453 rhs_fe_prj->exeTestHook =
2455 rhs_fe_current->exeTestHook =
2460 CHKERR prb_mng->removeDofsOnEntities(
2461 "DG_PROJECTION",
"L",
BitRefLevel().set(), remove_mask,
nullptr, 0,
2468 lhs_fe->getOpPtrVector().push_back(
2472 auto set_prj_from_child = [&](
auto rhs_fe_prj) {
2474 rhs_fe_prj->getOpPtrVector(), {L2});
2477 auto l_vec = boost::make_shared<MatrixDouble>();
2478 rhs_fe_prj->getOpPtrVector().push_back(
2482 auto get_parent_this = [&]() {
2483 auto fe_parent_this = boost::make_shared<DomianParentEle>(
mField);
2484 fe_parent_this->getOpPtrVector().push_back(
2486 return fe_parent_this;
2491 auto get_parents_fe_ptr = [&](
auto this_fe_ptr) {
2492 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
2494 parents_elems_ptr_vec.emplace_back(
2495 boost::make_shared<DomianParentEle>(
mField));
2497 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
2503 return parents_elems_ptr_vec[0];
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(
2515 auto set_prj_from_parent = [&](
auto rhs_fe_current) {
2518 auto l_vec = boost::make_shared<MatrixDouble>();
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;
2529 auto get_parents_fe_ptr = [&](
auto this_fe_ptr) {
2530 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
2532 parents_elems_ptr_vec.emplace_back(
2533 boost::make_shared<DomianParentEle>(
mField));
2535 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
2541 return parents_elems_ptr_vec[0];
2544 auto this_fe_ptr = get_parent_this();
2545 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
2555 rhs_fe_current->getOpPtrVector().push_back(reset_op_ptr);
2556 rhs_fe_current->getOpPtrVector().push_back(
2562 rhs_fe_current->getOpPtrVector().push_back(
new OpScalarFieldL(
"L", l_vec));
2565 set_prj_from_child(rhs_fe_prj);
2566 set_prj_from_parent(rhs_fe_current);
2568 boost::shared_ptr<FEMethod> null_fe;
2571 lhs_fe, null_fe, null_fe);
2575 rhs_fe_current, null_fe, null_fe);
2577 CHKERR KSPSetDM(ksp, sub_dm);
2579 CHKERR KSPSetDM(ksp, sub_dm);
2580 CHKERR KSPSetFromOptions(ksp);
2587 CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
2588 CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
2592 MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
2594 auto post_proc = [&](
auto dm,
auto out_name,
auto th_error) {
2597 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
2598 post_proc_fe->setTagsToTransfer({th_error});
2599 post_proc_fe->exeTestHook = test_mesh_bit;
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(
2609 post_proc_fe->getOpPtrVector().push_back(
2612 post_proc_fe->getOpPtrVector().push_back(
2616 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2620 {{
"GradL", l_grad_mat}},
2629 post_proc_fe->writeFile(out_name);
2633 if constexpr (
debug)
2634 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.
auto get_ntensor(T &base_mat)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
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 FE_DIM
[Define dimension]
constexpr int skeleton_bit
skeleton elements bit
FTensor::Index< 'I', DIM1 > I
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
FTensor::Index< 'J', DIM1 > J
FaceSideEle::UserDataOperator FaceSideEleOp
constexpr int aggregate_projection_bit
all bits for projection problem
LevelSet * level_set_raw_ptr
constexpr int aggregate_bit
all bits for advection problem
PipelineManager::ElementsAndOpsByDim< FE_DIM >::DomianParentEle DomianParentEle
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 jacobian 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.
DomainEle::UserDataOperator DomainEleOp
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
constexpr IntegrationType I
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< MatrixDouble > lPtr
MoFEMErrorCode iNtegrate(EntData &data)
Class dedicated to integrate operator.
boost::shared_ptr< MatrixDouble > velPtr
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)
boost::shared_ptr< MatrixDouble > 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)
Run calculations.
MoFEMErrorCode setBits(LevelSet &level_set, int l)
Set bit ref level to problem.
double getThreshold(const double max)
MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l)
Add bit to current element, so it aggregate all previious current elements.
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
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.
double getThreshold(const double max)
MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l)
Add bit to current element, so it aggregate all previious current elements.
boost::shared_ptr< double > maxPtr
MoFEMErrorCode setBits(LevelSet &level_set, int l)
Set bit ref level to problem.
WrapperClassInitalSolution(boost::shared_ptr< double > max_ptr)
MoFEMErrorCode runCalcs(LevelSet &level_set, int l)
Run calculations.
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< G >::OpMass< 1, DIM1 *DIM2 > OpMassLL
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 >::LinearForm< G >::OpSource< potential_velocity_field_dim, potential_velocity_field_dim > OpSourceV
MoFEMErrorCode readMesh()
read mesh
MoFEMErrorCode setUpProblem()
create fields, and set approximation order
MoFEMErrorCode initialiseFieldLevelSet(boost::function< double(double, double, double)> level_fun=get_level_set)
initialise field set
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)
static double get_velocity_potential(double x, double y, double z)
advection velocity field
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< G >::OpMass< potential_velocity_field_dim, potential_velocity_field_dim > OpMassVV
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
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpSource< 1, DIM1 *DIM2 > OpSourceL
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)
ForcesAndSourcesCore::UserDataOperator * getZeroLevelVelOp(boost::shared_ptr< MatrixDouble > vel_ptr)
Get operator calculating velocity on coarse mesh.
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpBaseTimesVector< 1, DIM1 *DIM2, 1 > OpScalarFieldL
std::tuple< double, Tag > evaluateError()
evaluate error
MoFEMErrorCode testSideFE()
test integration side elements
MoFEMErrorCode initialiseFieldVelocity(boost::function< double(double, double, double)> vel_fun=get_velocity_potential< FE_DIM >)
initialise potential velocity field
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
friend class UserDataOperator
@ 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.
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.
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.
friend class UserDataOperator