14static char help[] =
"...\n\n";
20 IntegrationType::GAUSS;
68 std::array<VectorInt, 2>
70 std::array<VectorInt, 2>
98 template <
int SPACE_DIM>
109 static double get_level_set(
const double x,
const double y,
const double z);
154 boost::shared_ptr<FaceSideEle>
155 getSideFE(boost::shared_ptr<SideData> side_data_ptr);
188 boost::function<
double(
double,
double,
double)> level_fun =
198 boost::function<
double(
double,
double,
double)> vel_fun =
199 get_velocity_potential<SPACE_DIM>);
240 simple->getBitRefLevel() =
261 ->synchroniseEntities(level);
270 return 0.05 * (*maxPtr);
295 ->synchroniseEntities(level);
327 I>::OpSource<potential_velocity_field_dim, potential_velocity_field_dim>;
342double LevelSet::get_velocity_potential<2>(
double x,
double y,
double z) {
343 return (x * x - 0.25) * (y * y - 0.25);
347 constexpr double xc = 0.1;
348 constexpr double yc = 0.;
349 constexpr double zc = 0.;
350 constexpr double r = 0.2;
351 return std::sqrt(pow(x - xc, 2) + pow(y - yc, 2) + pow(z - zc, 2)) - r;
359 if constexpr (
debug) {
366 maxPtr = boost::make_shared<double>(0);
383 boost::shared_ptr<VectorDouble> l_ptr,
384 boost::shared_ptr<VectorDouble> l_dot_ptr,
385 boost::shared_ptr<MatrixDouble> vel_ptr);
389 boost::shared_ptr<VectorDouble>
lPtr;
396 boost::shared_ptr<MatrixDouble> vel_ptr);
406 boost::shared_ptr<FaceSideEle> side_fe_ptr);
412 boost::shared_ptr<FaceSideEle>
420 boost::shared_ptr<FaceSideEle> side_fe_ptr);
426 boost::shared_ptr<FaceSideEle>
432int main(
int argc,
char *argv[]) {
441 moab::Core moab_core;
442 moab::Interface &moab = moab_core;
449 DMType dm_name =
"DMMOFEM";
453 auto core_log = logging::core::get();
479 simple->getAddSkeletonFE() =
true;
480 simple->getAddBoundaryFE() =
true;
486 auto set_problem_bit = [&]() {
491 start_mask[s] =
true;
509 if constexpr (
debug) {
511 CHKERR bit_mng->writeBitLevelByDim(
513 (proc_str +
"level_base.vtk").c_str(),
"VTK",
"");
549 boost::shared_ptr<VectorDouble> l_ptr,
550 boost::shared_ptr<VectorDouble> l_dot_ptr,
551 boost::shared_ptr<MatrixDouble> vel_ptr)
553 lPtr(l_ptr), lDotPtr(l_dot_ptr), velPtr(vel_ptr) {}
556 boost::shared_ptr<MatrixDouble> vel_ptr)
564 boost::shared_ptr<SideData> side_data_ptr,
565 boost::shared_ptr<FaceSideEle> side_fe_ptr)
567 sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr) {}
570 boost::shared_ptr<SideData> side_data_ptr,
571 boost::shared_ptr<FaceSideEle> side_fe_ptr)
573 sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr) {}
578 const auto nb_int_points = getGaussPts().size2();
579 const auto nb_dofs = data.
getIndices().size();
580 const auto nb_base_func = data.
getN().size2();
584 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
589 auto t_w = getFTensor0IntegrationWeight();
590 for (
auto gg = 0; gg != nb_int_points; ++gg) {
591 const auto alpha = t_w * getMeasure();
592 auto res0 = alpha * t_l_dot;
594 t_res1(
i) = (alpha * t_l) * t_vel(
i);
600 auto &nf = this->locF;
603 for (; rr != nb_dofs; ++rr) {
604 nf(rr) += res0 * t_base - t_res1(
i) * t_diff_base(
i);
608 for (; rr < nb_base_func; ++rr) {
621 const auto nb_int_points = getGaussPts().size2();
622 const auto nb_base_func = row_data.
getN().size2();
623 const auto nb_row_dofs = row_data.
getIndices().size();
624 const auto nb_col_dofs = col_data.
getIndices().size();
626 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
631 auto t_w = getFTensor0IntegrationWeight();
632 for (
auto gg = 0; gg != nb_int_points; ++gg) {
633 const auto alpha = t_w * getMeasure();
634 const auto beta = alpha * getTSa();
637 auto &mat = this->locMat;
640 for (; rr != nb_row_dofs; ++rr) {
642 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
644 (beta * t_row_base - alpha * (t_row_diff_base(
i) * t_vel(
i))) *
651 for (; rr < nb_base_func; ++rr) {
668 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
669 const auto in_the_loop =
670 sideFEPtr->nInTheLoop;
672 auto not_side = [](
auto s) {
678 &*base_mat.data().begin());
681 if (in_the_loop > 0) {
684 auto t_normal = getFTensor1Normal();
685 const auto nb_gauss_pts = getGaussPts().size2();
690 const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
694 resSkelton.resize(nb_rows,
false);
698 const auto opposite_s0 = not_side(s0);
699 const auto sense_row = sideDataPtr->senseMap[s0];
701 const auto opposite_sense_row = sideDataPtr->senseMap[opposite_s0];
702 if (sense_row * opposite_sense_row > 0)
704 "Should be opposite sign");
708 const auto nb_row_base_functions =
709 sideDataPtr->rowBaseSideMap[s0].size2();
711 auto t_w = getFTensor0IntegrationWeight();
716 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
717 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
720 for (
auto &t_l : arr_t_l)
722 for (
auto &t_vel : arr_t_vel)
727 if (nb_gauss_pts != sideDataPtr->rowBaseSideMap[s0].size1())
729 "Inconsistent number of DOFs");
732 auto t_row_base =
get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
733 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
736 const auto dot = sense_row * (t_normal(
i) * t_vel(
i));
737 const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
738 const auto l_upwind = arr_t_l[l_upwind_side];
739 const auto res = t_w * dot * l_upwind;
743 for (; rr != nb_rows; ++rr) {
744 resSkelton[rr] += t_row_base * res;
747 for (; rr < nb_row_base_functions; ++rr) {
752 CHKERR ::VecSetValues(getTSf(),
753 sideDataPtr->indicesRowSideMap[s0].size(),
754 &*sideDataPtr->indicesRowSideMap[s0].begin(),
755 &*resSkelton.begin(), ADD_VALUES);
769 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
770 const auto in_the_loop =
771 sideFEPtr->nInTheLoop;
773 auto not_side = [](
auto s) {
779 &*base_mat.data().begin());
782 if (in_the_loop > 0) {
785 auto t_normal = getFTensor1Normal();
786 const auto nb_gauss_pts = getGaussPts().size2();
791 const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
796 const auto opposite_s0 = not_side(s0);
797 const auto sense_row = sideDataPtr->senseMap[s0];
800 const auto nb_row_base_functions =
801 sideDataPtr->rowBaseSideMap[s0].size2();
806 const auto nb_cols = sideDataPtr->indicesColSideMap[s1].size();
807 const auto sense_col = sideDataPtr->senseMap[s1];
810 matSkeleton.resize(nb_rows, nb_cols,
false);
813 auto t_w = getFTensor0IntegrationWeight();
815 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
816 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
819 for (
auto &t_vel : arr_t_vel)
823 auto t_row_base =
get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
824 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
828 const auto dot = sense_row * (t_normal(
i) * t_vel(
i));
829 const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
830 const auto sense_upwind = sideDataPtr->senseMap[l_upwind_side];
831 auto res = t_w * dot;
835 if (s1 == l_upwind_side) {
836 for (; rr != nb_rows; ++rr) {
837 auto get_ntensor = [](
auto &base_mat,
auto gg,
auto bb) {
838 double *ptr = &base_mat(gg, bb);
842 get_ntensor(sideDataPtr->colBaseSideMap[s1], gg, 0);
843 const auto res_row = res * t_row_base;
846 for (
size_t cc = 0; cc != nb_cols; ++cc) {
847 matSkeleton(rr, cc) += res_row * t_col_base;
852 for (; rr < nb_row_base_functions; ++rr) {
857 CHKERR ::MatSetValues(getTSB(),
858 sideDataPtr->indicesRowSideMap[s0].size(),
859 &*sideDataPtr->indicesRowSideMap[s0].begin(),
860 sideDataPtr->indicesColSideMap[s1].size(),
861 &*sideDataPtr->indicesColSideMap[s1].begin(),
862 &*matSkeleton.data().begin(), ADD_VALUES);
872 auto get_parent_vel_this = [&]() {
873 auto parent_fe_ptr = boost::make_shared<DomianParentEle>(
mField);
875 parent_fe_ptr->getOpPtrVector(), {potential_velocity_space});
876 parent_fe_ptr->getOpPtrVector().push_back(
879 return parent_fe_ptr;
882 auto get_parents_vel_fe_ptr = [&](
auto this_fe_ptr) {
883 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
885 parents_elems_ptr_vec.emplace_back(
886 boost::make_shared<DomianParentEle>(
mField));
888 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
893 return parents_elems_ptr_vec[0];
896 auto this_fe_ptr = get_parent_vel_this();
897 auto parent_fe_ptr = get_parents_vel_fe_ptr(this_fe_ptr);
909 pip->getOpDomainRhsPipeline().clear();
911 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
912 pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
914 pip->getDomainLhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
915 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
918 pip->getDomainRhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
919 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
923 auto l_ptr = boost::make_shared<VectorDouble>();
924 auto l_dot_ptr = boost::make_shared<VectorDouble>();
925 auto vel_ptr = boost::make_shared<MatrixDouble>();
928 pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
929 pip->getOpDomainRhsPipeline().push_back(
931 pip->getOpDomainRhsPipeline().push_back(
934 pip->getOpDomainRhsPipeline().push_back(
938 pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
940 pip->getOpDomainLhsPipeline().push_back(
new OpLhsDomain(
"L", vel_ptr));
945boost::shared_ptr<FaceSideEle>
950 auto l_ptr = boost::make_shared<VectorDouble>();
951 auto vel_ptr = boost::make_shared<MatrixDouble>();
954 OpSideData(boost::shared_ptr<SideData> side_data_ptr)
956 sideDataPtr(side_data_ptr) {
957 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE],
false);
958 for (
auto t = moab::CN::TypeDimensionMap[
SPACE_DIM].first;
959 t <= moab::CN::TypeDimensionMap[
SPACE_DIM].second; ++
t)
960 doEntities[
t] =
true;
968 if ((CN::Dimension(row_type) ==
SPACE_DIM) &&
969 (CN::Dimension(col_type) ==
SPACE_DIM)) {
971 auto reset = [&](
auto nb_in_loop) {
972 sideDataPtr->feSideHandle[nb_in_loop] = 0;
973 sideDataPtr->indicesRowSideMap[nb_in_loop].clear();
974 sideDataPtr->indicesColSideMap[nb_in_loop].clear();
975 sideDataPtr->rowBaseSideMap[nb_in_loop].clear();
976 sideDataPtr->colBaseSideMap[nb_in_loop].clear();
977 sideDataPtr->senseMap[nb_in_loop] = 0;
980 const auto nb_in_loop = getFEMethod()->nInTheLoop;
982 for (
auto s : {0, 1})
985 sideDataPtr->currentFESide = nb_in_loop;
986 sideDataPtr->senseMap[nb_in_loop] = getSkeletonSense();
996 boost::shared_ptr<SideData> sideDataPtr;
1001 OpSideDataOnParent(boost::shared_ptr<SideData> side_data_ptr,
1002 boost::shared_ptr<VectorDouble> l_ptr,
1003 boost::shared_ptr<MatrixDouble> vel_ptr)
1005 sideDataPtr(side_data_ptr), lPtr(l_ptr), velPtr(vel_ptr) {
1006 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE],
false);
1007 for (
auto t = moab::CN::TypeDimensionMap[
SPACE_DIM].first;
1008 t <= moab::CN::TypeDimensionMap[
SPACE_DIM].second; ++
t)
1009 doEntities[
t] =
true;
1018 if ((CN::Dimension(row_type) ==
SPACE_DIM) &&
1019 (CN::Dimension(col_type) ==
SPACE_DIM)) {
1020 const auto nb_in_loop = sideDataPtr->currentFESide;
1021 sideDataPtr->feSideHandle[nb_in_loop] = getFEEntityHandle();
1022 sideDataPtr->indicesRowSideMap[nb_in_loop] = row_data.
getIndices();
1023 sideDataPtr->indicesColSideMap[nb_in_loop] = col_data.
getIndices();
1024 sideDataPtr->rowBaseSideMap[nb_in_loop] = row_data.
getN();
1025 sideDataPtr->colBaseSideMap[nb_in_loop] = col_data.
getN();
1026 (sideDataPtr->lVec)[nb_in_loop] = *lPtr;
1027 (sideDataPtr->velMat)[nb_in_loop] = *velPtr;
1030 if ((sideDataPtr->lVec)[nb_in_loop].size() !=
1031 (sideDataPtr->velMat)[nb_in_loop].size2())
1033 "Wrong number of integaration pts %d != %d",
1034 (sideDataPtr->lVec)[nb_in_loop].size(),
1035 (sideDataPtr->velMat)[nb_in_loop].size2());
1036 if ((sideDataPtr->velMat)[nb_in_loop].size1() !=
SPACE_DIM)
1038 "Wrong size of velocity vector size = %d",
1039 (sideDataPtr->velMat)[nb_in_loop].size1());
1043 (sideDataPtr->lVec)[1] = sideDataPtr->lVec[0];
1044 (sideDataPtr->velMat)[1] = (sideDataPtr->velMat)[0];
1047 if (sideDataPtr->rowBaseSideMap[0].size1() !=
1048 sideDataPtr->rowBaseSideMap[1].size1()) {
1050 "Wrong number of integration pt %d != %d",
1051 sideDataPtr->rowBaseSideMap[0].size1(),
1052 sideDataPtr->rowBaseSideMap[1].size1());
1054 if (sideDataPtr->colBaseSideMap[0].size1() !=
1055 sideDataPtr->colBaseSideMap[1].size1()) {
1057 "Wrong number of integration pt");
1070 boost::shared_ptr<SideData> sideDataPtr;
1071 boost::shared_ptr<VectorDouble> lPtr;
1072 boost::shared_ptr<MatrixDouble> velPtr;
1076 auto get_parent_this = [&]() {
1077 auto parent_fe_ptr = boost::make_shared<DomianParentEle>(
mField);
1079 parent_fe_ptr->getOpPtrVector(), {potential_velocity_space, L2});
1080 parent_fe_ptr->getOpPtrVector().push_back(
1082 parent_fe_ptr->getOpPtrVector().push_back(
1083 new OpSideDataOnParent(side_data_ptr, l_ptr, vel_ptr));
1084 return parent_fe_ptr;
1087 auto get_parents_fe_ptr = [&](
auto this_fe_ptr) {
1088 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
1090 parents_elems_ptr_vec.emplace_back(
1091 boost::make_shared<DomianParentEle>(
mField));
1093 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
1098 return parents_elems_ptr_vec[0];
1103 auto get_side_fe_ptr = [&]() {
1104 auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
1106 auto this_fe_ptr = get_parent_this();
1107 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
1109 side_fe_ptr->getOpPtrVector().push_back(
new OpSideData(side_data_ptr));
1111 side_fe_ptr->getOpPtrVector().push_back(
1119 return get_side_fe_ptr();
1127 pip->getOpSkeletonRhsPipeline().clear();
1129 pip->setSkeletonLhsIntegrationRule([](
int,
int,
int o) {
return 18; });
1130 pip->setSkeletonRhsIntegrationRule([](
int,
int,
int o) {
return 18; });
1132 pip->getSkeletonLhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
1133 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1136 pip->getSkeletonRhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
1137 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1141 auto side_data_ptr = boost::make_shared<SideData>();
1142 auto side_fe_ptr =
getSideFE(side_data_ptr);
1144 pip->getOpSkeletonRhsPipeline().push_back(
1146 pip->getOpSkeletonLhsPipeline().push_back(
1156 OpErrorSkel(boost::shared_ptr<FaceSideEle> side_fe_ptr,
1157 boost::shared_ptr<SideData> side_data_ptr,
1160 sideFEPtr(side_fe_ptr), sideDataPtr(side_data_ptr),
1161 errorSumPtr(error_sum_ptr), thError(th_error) {}
1167 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
1168 const auto in_the_loop =
1169 sideFEPtr->nInTheLoop;
1171 auto not_side = [](
auto s) {
1175 auto nb_gauss_pts = getGaussPts().size2();
1183 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
1184 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
1187 for (
auto &t_l : arr_t_l)
1189 for (
auto &t_vel : arr_t_vel)
1194 auto t_w = getFTensor0IntegrationWeight();
1195 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1196 e += t_w * getMeasure() *
1203 moab::Interface &moab =
1204 getNumeredEntFiniteElementPtr()->getBasicDataPtr()->moab;
1205 const void *tags_ptr[2];
1206 CHKERR moab.tag_get_by_ptr(thError, sideDataPtr->feSideHandle.data(), 2,
1208 for (
auto ff : {0, 1}) {
1209 *((
double *)tags_ptr[ff]) += e;
1211 CHKERR VecSetValue(errorSumPtr, 0, e, ADD_VALUES);
1218 boost::shared_ptr<FaceSideEle> sideFEPtr;
1219 boost::shared_ptr<SideData> sideDataPtr;
1230 MB_TAG_CREAT | MB_TAG_SPARSE,
1233 auto clear_tags = [&]() {
1242 auto evaluate_error = [&]() {
1244 auto skel_fe = boost::make_shared<BoundaryEle>(
mField);
1245 skel_fe->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
1246 auto side_data_ptr = boost::make_shared<SideData>();
1247 auto side_fe_ptr =
getSideFE(side_data_ptr);
1248 skel_fe->getOpPtrVector().push_back(
1249 new OpErrorSkel(side_fe_ptr, side_data_ptr, error_sum_ptr, th_error));
1252 skel_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1253 return fe_ptr->numeredEntFiniteElementPtr->
getBitRefLevel().test(
1258 simple->getSkeletonFEName(), skel_fe);
1263 auto assemble_and_sum = [](
auto vec) {
1271 auto propagate_error_to_parents = [&]() {
1275 auto fe_ptr = boost::make_shared<FEMethod>();
1276 fe_ptr->exeTestHook = [&](
FEMethod *fe_ptr) {
1277 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1281 fe_ptr->preProcessHook = []() {
return 0; };
1282 fe_ptr->postProcessHook = []() {
return 0; };
1283 fe_ptr->operatorHook = [&]() {
1286 auto fe_ent = fe_ptr->numeredEntFiniteElementPtr->getEnt();
1287 auto parent = fe_ptr->numeredEntFiniteElementPtr->getParentEnt();
1288 auto th_parent = fe_ptr->numeredEntFiniteElementPtr->getBasicDataPtr()
1289 ->th_RefParentHandle;
1292 CHKERR moab.tag_get_data(th_error, &fe_ent, 1, &error);
1295 [&](
auto fe_ent,
auto error) {
1298 CHKERR moab.tag_get_by_ptr(th_error, &fe_ent, 1,
1299 (
const void **)&e_ptr);
1303 CHKERR moab.tag_get_data(th_parent, &fe_ent, 1, &parent);
1304 if (parent != fe_ent && parent)
1305 CHKERR add_error(parent, *e_ptr);
1310 CHKERR add_error(parent, error);
1325 return std::make_tuple(assemble_and_sum(error_sum_ptr), th_error);
1343 DivergenceVol(boost::shared_ptr<VectorDouble> l_ptr,
1344 boost::shared_ptr<MatrixDouble> vel_ptr,
1351 const auto nb_dofs = data.
getIndices().size();
1353 const auto nb_gauss_pts = getGaussPts().size2();
1354 const auto t_w = getFTensor0IntegrationWeight();
1357 auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
1359 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
1360 for (
int rr = 0; rr != nb_dofs; ++rr) {
1361 div += getMeasure() * t_w * t_l * (t_diff(
i) * t_vel(
i));
1368 CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
1374 boost::shared_ptr<VectorDouble> lPtr;
1375 boost::shared_ptr<MatrixDouble> velPtr;
1384 DivergenceSkeleton(boost::shared_ptr<SideData> side_data_ptr,
1385 boost::shared_ptr<FaceSideEle> side_fe_ptr,
1388 sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr), divVec(div_vec) {}
1395 &*base_mat.data().begin());
1398 auto not_side = [](
auto s) {
1403 CHKERR loopSideFaces(
"dFE", *sideFEPtr);
1404 const auto in_the_loop =
1405 sideFEPtr->nInTheLoop;
1407 auto t_normal = getFTensor1Normal();
1408 const auto nb_gauss_pts = getGaussPts().size2();
1410 const auto nb_dofs = sideDataPtr->indicesRowSideMap[s0].size();
1412 auto t_base =
get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
1413 auto nb_row_base_functions = sideDataPtr->rowBaseSideMap[s0].size2();
1414 auto side_sense = sideDataPtr->senseMap[s0];
1415 auto opposite_s0 = not_side(s0);
1421 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
1422 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
1425 for (
auto &t_l : arr_t_l)
1427 for (
auto &t_vel : arr_t_vel)
1433 auto t_w = getFTensor0IntegrationWeight();
1434 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1438 const auto dot = (t_normal(
i) * t_vel(
i)) * side_sense;
1439 const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
1440 const auto l_upwind =
1441 arr_t_l[l_upwind_side];
1445 auto res = t_w * l_upwind * dot;
1449 for (; rr != nb_dofs; ++rr) {
1450 div += t_base * res;
1453 for (; rr < nb_row_base_functions; ++rr) {
1457 CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
1467 boost::shared_ptr<SideData> sideDataPtr;
1468 boost::shared_ptr<FaceSideEle> sideFEPtr;
1469 boost::shared_ptr<MatrixDouble> velPtr;
1473 auto vol_fe = boost::make_shared<DomainEle>(
mField);
1474 auto skel_fe = boost::make_shared<BoundaryEle>(
mField);
1476 vol_fe->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
1477 skel_fe->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
1482 auto l_ptr = boost::make_shared<VectorDouble>();
1483 auto vel_ptr = boost::make_shared<MatrixDouble>();
1484 auto side_data_ptr = boost::make_shared<SideData>();
1485 auto side_fe_ptr =
getSideFE(side_data_ptr);
1488 vol_fe->getOpPtrVector(), {potential_velocity_space, L2});
1489 vol_fe->getOpPtrVector().push_back(
1492 vol_fe->getOpPtrVector().push_back(
1493 new DivergenceVol(l_ptr, vel_ptr, div_vol_vec));
1495 skel_fe->getOpPtrVector().push_back(
1496 new DivergenceSkeleton(side_data_ptr, side_fe_ptr, div_skel_vec));
1499 auto dm =
simple->getDM();
1508 [](
double x,
double y,
double) {
return x - y; });
1510 [](
double x,
double y,
double) {
return x - y; });
1512 vol_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1513 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1516 skel_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1517 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1525 auto assemble_and_sum = [](
auto vec) {
1533 auto div_vol = assemble_and_sum(div_vol_vec);
1534 auto div_skel = assemble_and_sum(div_skel_vec);
1536 auto eps = std::abs((div_vol - div_skel) / (div_vol + div_skel));
1538 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Testing divergence volume: " << div_vol;
1539 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Testing divergence skeleton: " << div_skel
1540 <<
" relative difference: " <<
eps;
1542 constexpr double eps_err = 1e-6;
1545 "No consistency between skeleton integral and volume integral");
1563 auto post_proc = [&](
auto dm,
auto f_res,
auto out_name) {
1566 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
1570 auto l_vec = boost::make_shared<VectorDouble>();
1571 post_proc_fe->getOpPtrVector().push_back(
1574 post_proc_fe->getOpPtrVector().push_back(
1578 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1590 post_proc_fe->writeFile(out_name);
1594 constexpr double eps = 1e-4;
1597 opt->setRandomFields(
simple->getDM(), {{
"L", {-1, 1}}, {
"V", {-1, 1}}});
1598 auto dot_x = opt->setRandomFields(
simple->getDM(), {{
"L", {-1, 1}}});
1599 auto diff_x = opt->setRandomFields(
simple->getDM(), {{
"L", {-1, 1}}});
1601 auto test_domain_ops = [&](
auto fe_name,
auto lhs_pipeline,
1602 auto rhs_pipeline) {
1605 auto diff_res = opt->checkCentralFiniteDifference(
1606 simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
1609 if constexpr (
debug) {
1613 CHKERR post_proc(
simple->getDM(), diff_res,
"tangent_op_error.h5m");
1619 CHKERR VecNorm(diff_res, NORM_2, &fnorm);
1621 "Test consistency of tangent matrix %3.4e", fnorm);
1623 constexpr double err = 1e-9;
1626 "Norm of directional derivative too large err = %3.4e", fnorm);
1631 CHKERR test_domain_ops(
simple->getDomainFEName(), pip->getDomainLhsFE(),
1632 pip->getDomainRhsFE());
1633 CHKERR test_domain_ops(
simple->getSkeletonFEName(), pip->getSkeletonLhsFE(),
1634 pip->getSkeletonRhsFE());
1640 boost::function<
double(
double,
double,
double)> level_fun) {
1649 boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(
mField);
1650 boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(
mField);
1651 auto swap_fe = [&]() {
1652 lhs_fe.swap(pip->getDomainLhsFE());
1653 rhs_fe.swap(pip->getDomainRhsFE());
1657 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1658 pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1670 CHKERR prb_mng->removeDofsOnEntities(
"LEVELSET_POJECTION",
"L",
1672 auto test_bit_ele = [&](
FEMethod *fe_ptr) {
1673 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1676 pip->getDomainLhsFE()->exeTestHook = test_bit_ele;
1677 pip->getDomainRhsFE()->exeTestHook = test_bit_ele;
1680 pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
1682 pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
1683 pip->getOpDomainLhsPipeline().push_back(
new OpMassLL(
"L",
"L"));
1684 pip->getOpDomainRhsPipeline().push_back(
new OpSourceL(
"L", level_fun));
1688 auto ksp = pip->createKSP(sub_dm);
1689 CHKERR KSPSetDM(ksp, sub_dm);
1690 CHKERR KSPSetFromOptions(ksp);
1697 CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
1698 CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
1702 MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
1706 std::vector<Tag> tags{th_error};
1708 "PARALLEL=WRITE_PART", &fe_meshset, 1,
1709 &*tags.begin(), tags.size());
1712 auto post_proc = [&](
auto dm,
auto out_name,
auto th_error) {
1715 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
1716 post_proc_fe->setTagsToTransfer({th_error});
1717 post_proc_fe->exeTestHook = test_bit_ele;
1721 auto l_vec = boost::make_shared<VectorDouble>();
1722 auto l_grad_mat = boost::make_shared<MatrixDouble>();
1724 post_proc_fe->getOpPtrVector(), {potential_velocity_space, L2});
1725 post_proc_fe->getOpPtrVector().push_back(
1727 post_proc_fe->getOpPtrVector().push_back(
1730 post_proc_fe->getOpPtrVector().push_back(
1734 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1738 {{
"GradL", l_grad_mat}},
1746 post_proc_fe->writeFile(out_name);
1750 if constexpr (
debug)
1751 CHKERR post_proc(sub_dm,
"initial_level_set.h5m", th_error);
1759 boost::function<
double(
double,
double,
double)> vel_fun) {
1768 boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(
mField);
1769 boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(
mField);
1770 auto swap_fe = [&]() {
1771 lhs_fe.swap(pip->getDomainLhsFE());
1772 rhs_fe.swap(pip->getDomainRhsFE());
1776 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1777 pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
1790 CHKERR prb_mng->removeDofsOnEntities(
"VELOCITY_PROJECTION",
"V",
1793 auto test_bit = [&](
FEMethod *fe_ptr) {
1794 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(0);
1796 pip->getDomainLhsFE()->exeTestHook = test_bit;
1797 pip->getDomainRhsFE()->exeTestHook = test_bit;
1800 pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
1802 pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
1804 pip->getOpDomainLhsPipeline().push_back(
new OpMassVV(
"V",
"V"));
1805 pip->getOpDomainRhsPipeline().push_back(
new OpSourceV(
"V", vel_fun));
1807 auto ksp = pip->createKSP(sub_dm);
1808 CHKERR KSPSetDM(ksp, sub_dm);
1809 CHKERR KSPSetFromOptions(ksp);
1816 CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
1817 CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
1820 auto post_proc = [&](
auto dm,
auto out_name) {
1823 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
1824 post_proc_fe->exeTestHook = test_bit;
1827 post_proc_fe->getOpPtrVector(), {potential_velocity_space});
1832 auto l_vec = boost::make_shared<VectorDouble>();
1833 auto potential_vec = boost::make_shared<VectorDouble>();
1834 auto velocity_mat = boost::make_shared<MatrixDouble>();
1836 post_proc_fe->getOpPtrVector().push_back(
1838 post_proc_fe->getOpPtrVector().push_back(
1842 post_proc_fe->getOpPtrVector().push_back(
1846 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1848 {{
"VelocityPotential", potential_vec}},
1850 {{
"Velocity", velocity_mat}},
1858 "3d case not implemented");
1863 post_proc_fe->writeFile(out_name);
1867 if constexpr (
debug)
1868 CHKERR post_proc(sub_dm,
"initial_velocity_potential.h5m");
1903 auto add_post_proc_fe = [&]() {
1904 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
1909 "Error", 1, MB_TYPE_DOUBLE, th_error, MB_TAG_CREAT | MB_TAG_SPARSE,
1911 post_proc_fe->setTagsToTransfer({th_error});
1913 post_proc_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
1914 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1919 auto l_vec = boost::make_shared<VectorDouble>();
1920 auto vel_ptr = boost::make_shared<MatrixDouble>();
1923 post_proc_fe->getOpPtrVector(), {H1});
1924 post_proc_fe->getOpPtrVector().push_back(
1928 post_proc_fe->getOpPtrVector().push_back(
1932 post_proc_fe->getPostProcMesh(),
1934 post_proc_fe->getMapGaussPts(),
1943 return post_proc_fe;
1946 auto post_proc_fe = add_post_proc_fe();
1948 auto set_time_monitor = [&](
auto dm,
auto ts) {
1949 auto monitor_ptr = boost::make_shared<FEMethod>();
1951 monitor_ptr->preProcessHook = []() {
return 0; };
1952 monitor_ptr->operatorHook = []() {
return 0; };
1953 monitor_ptr->postProcessHook = [&]() {
1958 "Null pointer for post proc element");
1962 CHKERR post_proc_fe->writeFile(
1964 boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
".h5m");
1968 boost::shared_ptr<FEMethod> null;
1975 auto ts = pip->createTSIM(sub_dm);
1977 auto set_solution = [&](
auto ts) {
1986 auto monitor_pt = set_time_monitor(sub_dm, ts);
1987 CHKERR TSSetFromOptions(ts);
1995 auto ts_pre_step = [](TS ts) {
2002 MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
2004 auto get_norm = [&](
auto x) {
2006 CHKERR VecNorm(x, NORM_2, &nrm);
2010 auto set_solution = [&](
auto ts) {
2018 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
2019 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
2022 <<
"Problem " << prb_ptr->getName() <<
" solution vector norm "
2024 CHKERR TSSetSolution(ts, x);
2029 auto refine_and_project = [&](
auto ts) {
2047 ->removeDofsOnEntities(
"ADVECTION",
"L",
BitRefLevel().set(),
2053 auto ts_reset_theta = [&](
auto ts) {
2070 CHKERR refine_and_project(ts);
2071 CHKERR ts_reset_theta(ts);
2076 auto ts_post_step = [](TS ts) {
return 0; };
2078 CHKERR TSSetPreStep(ts, ts_pre_step);
2079 CHKERR TSSetPostStep(ts, ts_post_step);
2081 CHKERR TSSolve(ts, NULL);
2091 ParallelComm *pcomm =
2102 meshset_ptr->get_ptr(), 1);
2107 auto get_refined_elements_meshset = [&](
auto bit,
auto mask) {
2113 "get error handle");
2114 std::vector<double> errors(fe_ents.size());
2116 mField.
get_moab().tag_get_data(th_error, fe_ents, &*errors.begin()),
2118 auto it = std::max_element(errors.begin(), errors.end());
2120 MPI_Allreduce(&*it, &max, 1, MPI_DOUBLE, MPI_MAX,
mField.
get_comm());
2121 MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Max error: " << max;
2122 auto threshold = wp.getThreshold(max);
2124 std::vector<EntityHandle> fe_to_refine;
2125 fe_to_refine.reserve(fe_ents.size());
2127 auto fe_it = fe_ents.begin();
2128 auto error_it = errors.begin();
2129 for (
auto i = 0;
i != fe_ents.size(); ++
i) {
2130 if (*error_it > threshold) {
2131 fe_to_refine.push_back(*fe_it);
2138 ents.insert_list(fe_to_refine.begin(), fe_to_refine.end());
2142 auto get_neighbours_by_bridge_vertices = [&](
auto &&ents) {
2146 moab::Interface::UNION);
2151 ents = get_neighbours_by_bridge_vertices(ents);
2157 "add entities to meshset");
2159 (proc_str +
"_fe_to_refine.vtk").c_str(),
"VTK",
"",
2160 meshset_ptr->get_ptr(), 1);
2168 auto refine_mesh = [&](
auto l,
auto &&fe_to_refine) {
2174 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2177 fe_to_refine = intersect(level_ents, fe_to_refine);
2179 level_ents = subtract(level_ents, fe_to_refine);
2182 Range fe_to_refine_children;
2183 bit_mng->updateRangeByChildren(fe_to_refine, fe_to_refine_children);
2185 fe_to_refine_children =
2186 fe_to_refine_children.subset_by_dimension(
SPACE_DIM);
2187 level_ents.merge(fe_to_refine_children);
2189 auto fix_neighbour_level = [&](
auto ll) {
2191 auto level_ll = level_ents;
2195 CHKERR skin.find_skin(0, level_ll,
false, skin_edges);
2197 for (
auto lll = 0; lll <= ll; ++lll) {
2198 CHKERR bit_mng->updateRangeByParent(skin_edges, skin_parents);
2199 skin_edges = skin_parents;
2202 for (
auto lll = 0; lll <= ll - 2; ++lll) {
2203 bad_bit[lll] =
true;
2207 Range skin_adj_ents;
2210 moab::Interface::UNION);
2213 skin_adj_ents = intersect(skin_adj_ents, level_ents);
2214 if (!skin_adj_ents.empty()) {
2215 level_ents = subtract(level_ents, skin_adj_ents);
2216 Range skin_adj_ents_children;
2217 bit_mng->updateRangeByChildren(skin_adj_ents, skin_adj_ents_children);
2218 level_ents.merge(skin_adj_ents_children);
2223 CHKERR fix_neighbour_level(
l);
2232 level_ents.subset_by_dimension(
SPACE_DIM), level_ents,
true);
2235 level_ents.subset_by_dimension(
SPACE_DIM), d,
false, level_ents,
2236 moab::Interface::UNION);
2248 CHKERR bit_mng->writeBitLevelByDim(
2250 (boost::lexical_cast<std::string>(
l) +
"_" + proc_str +
"_ref_mesh.vtk")
2259 auto set_skelton_bit = [&](
auto l) {
2269 Range level_edges_parents;
2270 CHKERR bit_mng->updateRangeByParent(level_edges, level_edges_parents);
2271 level_edges_parents =
2272 level_edges_parents.subset_by_dimension(
SPACE_DIM - 1);
2273 CHKERR bit_mng->filterEntitiesByRefLevel(
2277 auto parent_skeleton = intersect(level_edges, level_edges_parents);
2278 auto skeleton = subtract(level_edges, level_edges_parents);
2283 moab::Interface::UNION);
2291 CHKERR bit_mng->writeBitLevel(
2293 (boost::lexical_cast<std::string>(
l) +
"_" + proc_str +
"_skeleton.vtk")
2301 Range level0_current;
2305 Range level0_aggregate;
2311 start_mask[s] =
true;
2312 CHKERR bit_mng->lambdaBitRefLevel(
2327 CHKERR wp.setBits(*
this, 0);
2328 CHKERR wp.runCalcs(*
this, 0);
2330 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Process level: " <<
l;
2331 CHKERR refine_mesh(
l + 1, get_refined_elements_meshset(
2333 CHKERR set_skelton_bit(
l + 1);
2334 CHKERR wp.setAggregateBit(*
this,
l + 1);
2335 CHKERR wp.setBits(*
this,
l + 1);
2336 CHKERR wp.runCalcs(*
this,
l + 1);
2351 auto lhs_fe = boost::make_shared<DomainEle>(
mField);
2352 auto rhs_fe_prj = boost::make_shared<DomainEle>(
mField);
2353 auto rhs_fe_current = boost::make_shared<DomainEle>(
mField);
2355 lhs_fe->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
2356 rhs_fe_prj->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
2357 rhs_fe_current->getRuleHook = [](int, int,
int o) {
return 3 *
o; };
2376 CHKERR bit_mng->updateRangeByParent(prj_ents, prj_ents);
2378 current_ents = subtract(current_ents, prj_ents);
2380 auto test_mesh_bit = [&](
FEMethod *fe_ptr) {
2381 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2384 auto test_prj_bit = [&](
FEMethod *fe_ptr) {
2385 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2388 auto test_current_bit = [&](
FEMethod *fe_ptr) {
2389 return current_ents.find(fe_ptr->getFEEntityHandle()) != current_ents.end();
2392 lhs_fe->exeTestHook = test_mesh_bit;
2393 rhs_fe_prj->exeTestHook = test_prj_bit;
2394 rhs_fe_current->exeTestHook = test_current_bit;
2398 CHKERR prb_mng->removeDofsOnEntities(
2399 "DG_PROJECTION",
"L",
BitRefLevel().set(), remove_mask,
nullptr, 0,
2403 lhs_fe->getOpPtrVector(), {potential_velocity_space, L2});
2404 lhs_fe->getOpPtrVector().push_back(
new OpMassLL(
"L",
"L"));
2406 auto l_vec = boost::make_shared<VectorDouble>();
2408 auto set_prj_from_child = [&](
auto rhs_fe_prj) {
2410 rhs_fe_prj->getOpPtrVector(), {potential_velocity_space, L2});
2411 rhs_fe_prj->getOpPtrVector().push_back(
2413 auto get_parent_this = [&]() {
2414 auto fe_parent_this = boost::make_shared<DomianParentEle>(
mField);
2415 fe_parent_this->getOpPtrVector().push_back(
2417 return fe_parent_this;
2420 auto get_parents_fe_ptr = [&](
auto this_fe_ptr) {
2421 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
2423 parents_elems_ptr_vec.emplace_back(
2424 boost::make_shared<DomianParentEle>(
mField));
2426 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
2432 return parents_elems_ptr_vec[0];
2435 auto this_fe_ptr = get_parent_this();
2436 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
2437 rhs_fe_prj->getOpPtrVector().push_back(
2443 auto set_prj_from_parent = [&](
auto rhs_fe_current) {
2444 auto get_parent_this = [&]() {
2445 auto fe_parent_this = boost::make_shared<DomianParentEle>(
mField);
2446 fe_parent_this->getOpPtrVector().push_back(
2448 return fe_parent_this;
2451 auto get_parents_fe_ptr = [&](
auto this_fe_ptr) {
2452 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
2454 parents_elems_ptr_vec.emplace_back(
2455 boost::make_shared<DomianParentEle>(
mField));
2457 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
2463 return parents_elems_ptr_vec[0];
2466 auto this_fe_ptr = get_parent_this();
2467 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
2477 rhs_fe_current->getOpPtrVector().push_back(reset_op_ptr);
2478 rhs_fe_current->getOpPtrVector().push_back(
2482 rhs_fe_current->getOpPtrVector().push_back(
new OpScalarFieldL(
"L", l_vec));
2485 set_prj_from_child(rhs_fe_prj);
2486 set_prj_from_parent(rhs_fe_current);
2488 boost::shared_ptr<FEMethod> null_fe;
2491 lhs_fe, null_fe, null_fe);
2495 rhs_fe_current, null_fe, null_fe);
2497 CHKERR KSPSetDM(ksp, sub_dm);
2499 CHKERR KSPSetDM(ksp, sub_dm);
2500 CHKERR KSPSetFromOptions(ksp);
2507 CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
2508 CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
2512 MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
2514 auto post_proc = [&](
auto dm,
auto out_name,
auto th_error) {
2517 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
2518 post_proc_fe->setTagsToTransfer({th_error});
2519 post_proc_fe->exeTestHook = test_mesh_bit;
2523 auto l_vec = boost::make_shared<VectorDouble>();
2524 auto l_grad_mat = boost::make_shared<MatrixDouble>();
2526 post_proc_fe->getOpPtrVector(), {potential_velocity_space, L2});
2527 post_proc_fe->getOpPtrVector().push_back(
2529 post_proc_fe->getOpPtrVector().push_back(
2532 post_proc_fe->getOpPtrVector().push_back(
2536 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2540 {{
"GradL", l_grad_mat}},
2548 post_proc_fe->writeFile(out_name);
2552 if constexpr (
debug)
2553 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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1, 1 > OpBaseTimesScalar
auto smartCreateDMMatrix(DM dm)
Get smart matrix from DM.
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.
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.
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
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
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomianParentEle DomianParentEle
FTensor::Index< 'j', SPACE_DIM > j
constexpr int current_bit
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
constexpr IntegrationType I
DomainEle::UserDataOperator DomainEleOp
LevelSet * level_set_raw_ptr
constexpr int aggregate_bit
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.
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
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 smartGetDMKspCtx(DM dm)
Get KSP context data structure used by DM.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
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 refinment.
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.