|
| v0.14.0
|
- Examples
- level_set.cpp.
Definition at line 73 of file level_set.cpp.
◆ AssemblyBoundaryEleOp
◆ AssemblyDomainEleOp
◆ MatSideArray
◆ OpMassLL
◆ OpMassVV
◆ OpScalarFieldL
◆ OpSourceL
◆ OpSourceV
◆ ElementSide
◆ LevelSet()
◆ dgProjection()
dg level set projection
- Parameters
-
- Returns
- MoFEMErrorCode
- Examples
- level_set.cpp.
Definition at line 2401 of file level_set.cpp.
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);
2548 reset_op_ptr->doWorkRhsHook = [&](
DataOperator *op_ptr,
int, EntityType,
2550 l_vec->resize(
static_cast<DomainEleOp *
>(op_ptr)->getGaussPts().size2(),
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);
◆ evaluateError()
std::tuple< double, Tag > LevelSet::evaluateError |
( |
| ) |
|
|
private |
evaluate error
- Returns
- MoFEMErrorCode
- Examples
- level_set.cpp.
Definition at line 1203 of file level_set.cpp.
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);
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);
◆ get_level_set()
inital level set, i.e. advected filed
- Parameters
-
- Returns
- double
- Examples
- level_set.cpp.
Definition at line 378 of file level_set.cpp.
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;
◆ get_velocity_potential() [1/2]
advection velocity field
- Note
- in current implementation is assumed that advection field has zero normal component.
-
function define a vector velocity potential field, curl of potential field gives velocity, thus velocity is divergence free.
- Template Parameters
-
- Parameters
-
- Returns
- auto
◆ get_velocity_potential() [2/2]
Definition at line 374 of file level_set.cpp.
375 return (x * x - 0.25) * (y * y - 0.25);
◆ getSideFE()
boost::shared_ptr< FaceSideEle > LevelSet::getSideFE |
( |
boost::shared_ptr< SideData > |
side_data_ptr | ) |
|
|
private |
create side element to assemble data from sides
- Parameters
-
- Returns
- boost::shared_ptr<FaceSideEle>
- Examples
- level_set.cpp.
Definition at line 997 of file level_set.cpp.
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;
1015 MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
1016 EntityType col_type,
EntData &row_data,
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;
1064 MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
1065 EntityType col_type,
EntData &row_data,
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();
◆ getZeroLevelVelOp()
Get operator calculating velocity on coarse mesh.
- Parameters
-
- Returns
- DomainEleOp*
- Examples
- level_set.cpp.
Definition at line 922 of file level_set.cpp.
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);
◆ initialiseFieldLevelSet()
initialise field set
- Parameters
-
- Returns
- MoFEMErrorCode
- Examples
- level_set.cpp.
Definition at line 1693 of file level_set.cpp.
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);
◆ initialiseFieldVelocity()
initialise potential velocity field
- Parameters
-
- Returns
- MoFEMErrorCode
- Examples
- level_set.cpp.
Definition at line 1814 of file level_set.cpp.
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");
◆ pushOpDomain()
push operators to integrate operators on domain
- Returns
- MoFEMErrorCode
- Examples
- level_set.cpp.
Definition at line 954 of file level_set.cpp.
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(
986 new OpRhsDomain(
"L", l_ptr, l_dot_ptr, vel_ptr));
991 pip->getOpDomainLhsPipeline().push_back(
new OpLhsDomain(
"L", vel_ptr));
◆ pushOpSkeleton()
push operator to integrate on skeleton
- Returns
- MoFEMErrorCode
- Examples
- level_set.cpp.
Definition at line 1173 of file level_set.cpp.
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(
1196 new OpRhsSkeleton(side_data_ptr, side_fe_ptr));
1197 pip->getOpSkeletonLhsPipeline().push_back(
1198 new OpLhsSkeleton(side_data_ptr, side_fe_ptr));
◆ readMesh()
read mesh
- Returns
- MoFEMErrorCode
- Examples
- level_set.cpp.
Definition at line 502 of file level_set.cpp.
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",
"");
◆ refineMesh()
- Examples
- level_set.cpp.
Definition at line 2147 of file level_set.cpp.
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);
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);
◆ runProblem()
◆ setUpProblem()
◆ solveAdvection()
solve advection problem
- Returns
- * MoFEMErrorCode
- Examples
- level_set.cpp.
Definition at line 1934 of file level_set.cpp.
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);
◆ testOp()
test consistency between tangent matrix and the right hand side vectors
- Returns
- MoFEMErrorCode
- Examples
- level_set.cpp.
Definition at line 1602 of file level_set.cpp.
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());
◆ testSideFE()
test integration side elements
test side element
Check consistency between volume and skeleton integral.
- Returns
- MoFEMErrorCode
Check consistency between volume and skeleton integral
- Returns
- MoFEMErrorCode
calculate volume
calculate skeleton integral
Set up problem such that gradient of level set field is orthogonal to velocity field. Then volume and skeleton integral should yield the same value.
- Examples
- level_set.cpp.
Definition at line 1387 of file level_set.cpp.
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");
◆ maxPtr
boost::shared_ptr<double> LevelSet::maxPtr |
|
private |
◆ mField
The documentation for this struct was generated from the following file:
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
constexpr int aggregate_projection_bit
all bits for projection problem
Data on single entity (This is passed as argument to DataOperator::doWork)
#define MYPCOMM_INDEX
default communicator number PCOMM
MoFEMErrorCode readMesh()
read mesh
virtual EntityHandle get_finite_element_meshset(const std::string name) const =0
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
base operator to do operations at Gauss Pt. level
boost::shared_ptr< double > maxPtr
boost::shared_ptr< FaceSideEle > getSideFE(boost::shared_ptr< SideData > side_data_ptr)
create side element to assemble data from sides
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Problem manager is used to build and partition problems.
structure for User Loop Methods on finite elements
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< G >::OpMass< 1, DIM1 *DIM2 > OpMassLL
Operator to execute finite element instance on parent element. This operator is typically used to pro...
virtual MPI_Comm & get_comm() const =0
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
@ L2
field with C-1 continuity
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpBaseTimesVector< 1, DIM1 *DIM2, 1 > OpScalarFieldL
MoFEMErrorCode solveAdvection()
solve advection problem
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
virtual int get_comm_rank() const =0
constexpr int aggregate_bit
all bits for advection problem
FaceSideEle::UserDataOperator FaceSideEleOp
PipelineManager interface.
MoFEMErrorCode initialiseFieldLevelSet(boost::function< double(double, double, double)> level_fun=get_level_set)
initialise field set
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpSource< 1, DIM1 *DIM2 > OpSourceL
FTensor::Index< 'J', DIM1 > J
constexpr size_t potential_velocity_field_dim
auto createKSP(MPI_Comm comm)
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.
Simple interface for fast problem set-up.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
FTensor::Index< 'I', DIM1 > I
constexpr int FE_DIM
[Define dimension]
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
DeprecatedCoreInterface Interface
MoFEMErrorCode testOp()
test consistency between tangent matrix and the right hand side vectors
MoFEMErrorCode pushOpSkeleton()
push operator to integrate on skeleton
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpSource< potential_velocity_field_dim, potential_velocity_field_dim > OpSourceV
constexpr int current_bit
dofs bit used to do calculations
#define CHKERR
Inline error check.
auto createDMVector(DM dm)
Get smart vector from DM.
virtual moab::Interface & get_moab()=0
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
PetscErrorCode DMSubDMSetUp_MoFEM(DM subdm)
#define MOFEM_LOG_C(channel, severity, format,...)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
DomainEle::UserDataOperator DomainEleOp
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
MoFEMErrorCode setUpProblem()
create fields, and set approximation order
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Get value at integration points for scalar field.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto get_ntensor(T &base_mat)
LevelSet * level_set_raw_ptr
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.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< G >::OpMass< potential_velocity_field_dim, potential_velocity_field_dim > OpMassVV
constexpr double t
plate stiffness
Get time direvarive values at integration pts for tensor filed rank 2, i.e. matrix field.
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Add operators pushing bases from local to physical configuration.
boost::ptr_deque< UserDataOperator > & getOpSkeletonLhsPipeline()
Get the Op Skeleton Lhs Pipeline object.
FTensor::Index< 'i', SPACE_DIM > i
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
std::tuple< double, Tag > evaluateError()
evaluate error
MoFEM::Interface & mField
integrate skeleton operators on khs
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
constexpr int skeleton_bit
skeleton elements bit
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
MoFEMErrorCode initialiseFieldVelocity(boost::function< double(double, double, double)> vel_fun=get_velocity_potential< FE_DIM >)
initialise potential velocity field
MoFEMErrorCode refineMesh(WrapperClass &&wp)
MoFEMErrorCode dgProjection(const int prj_bit=projection_bit)
dg level set projection
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.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
auto getProblemPtr(DM dm)
get problem pointer from DM
@ MOFEM_DATA_INCONSISTENCY
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
ForcesAndSourcesCore::UserDataOperator * getZeroLevelVelOp(boost::shared_ptr< MatrixDouble > vel_ptr)
Get operator calculating velocity on coarse mesh.
const double D
diffusivity
MoFEMErrorCode pushOpDomain()
push operators to integrate operators on domain
constexpr int projection_bit
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
@ MOFEM_ATOM_TEST_INVALID
constexpr auto make_array(Arg &&...arg)
Create Array.
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
constexpr FieldSpace potential_velocity_space
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Calculate curl of vector field.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'l', 3 > l
Post post-proc data at points from hash maps.
MoFEMErrorCode testSideFE()
test integration side elements