v0.14.0
Classes | Public Member Functions | Private Types | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
LevelSet Struct Reference
Collaboration diagram for LevelSet:
[legend]

Classes

struct  OpLhsDomain
 
struct  OpLhsSkeleton
 
struct  OpRhsDomain
 
struct  OpRhsSkeleton
 
struct  SideData
 data structure carrying information on skeleton on both sides. More...
 
struct  WrapperClass
 Wrapper executing stages while mesh refinement. More...
 
struct  WrapperClassErrorProjection
 Use peculated errors on all levels while mesh projection. More...
 
struct  WrapperClassInitalSolution
 Used to execute inital mesh approximation while mesh refinement. More...
 

Public Member Functions

 LevelSet (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 

Private Types

enum  ElementSide { LEFT_SIDE = 0, RIGHT_SIDE = 1 }
 
using MatSideArray = std::array< MatrixDouble, 2 >
 
using AssemblyDomainEleOp = FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase
 
using OpMassLL = FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< G >::OpMass< 1, DIM1 *DIM2 >
 
using OpSourceL = FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpSource< 1, DIM1 *DIM2 >
 
using OpMassVV = FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< G >::OpMass< potential_velocity_field_dim, potential_velocity_field_dim >
 
using OpSourceV = FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpSource< potential_velocity_field_dim, potential_velocity_field_dim >
 
using OpScalarFieldL = FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpBaseTimesVector< 1, DIM1 *DIM2, 1 >
 
using AssemblyBoundaryEleOp = FormsIntegrators< BoundaryEleOp >::Assembly< A >::OpBase
 

Private Member Functions

MoFEMErrorCode readMesh ()
 read mesh More...
 
MoFEMErrorCode setUpProblem ()
 create fields, and set approximation order More...
 
MoFEMErrorCode pushOpDomain ()
 push operators to integrate operators on domain More...
 
std::tuple< double, Tag > evaluateError ()
 evaluate error More...
 
ForcesAndSourcesCore::UserDataOperatorgetZeroLevelVelOp (boost::shared_ptr< MatrixDouble > vel_ptr)
 Get operator calculating velocity on coarse mesh. More...
 
boost::shared_ptr< FaceSideElegetSideFE (boost::shared_ptr< SideData > side_data_ptr)
 create side element to assemble data from sides More...
 
MoFEMErrorCode pushOpSkeleton ()
 push operator to integrate on skeleton More...
 
MoFEMErrorCode testSideFE ()
 test integration side elements More...
 
MoFEMErrorCode testOp ()
 test consistency between tangent matrix and the right hand side vectors More...
 
MoFEMErrorCode initialiseFieldLevelSet (boost::function< double(double, double, double)> level_fun=get_level_set)
 initialise field set More...
 
MoFEMErrorCode initialiseFieldVelocity (boost::function< double(double, double, double)> vel_fun=get_velocity_potential< FE_DIM >)
 initialise potential velocity field More...
 
MoFEMErrorCode dgProjection (const int prj_bit=projection_bit)
 dg level set projection More...
 
MoFEMErrorCode solveAdvection ()
 solve advection problem More...
 
MoFEMErrorCode refineMesh (WrapperClass &&wp)
 

Static Private Member Functions

template<int FE_DIM>
static double get_velocity_potential (double x, double y, double z)
 advection velocity field More...
 
static double get_level_set (const double x, const double y, const double z)
 inital level set, i.e. advected filed More...
 
template<>
double get_velocity_potential (double x, double y, double z)
 

Private Attributes

MoFEM::InterfacemField
 integrate skeleton operators on khs More...
 
boost::shared_ptr< doublemaxPtr
 

Detailed Description

Examples
level_set.cpp.

Definition at line 73 of file level_set.cpp.

Member Typedef Documentation

◆ AssemblyBoundaryEleOp

Definition at line 365 of file level_set.cpp.

◆ AssemblyDomainEleOp

Definition at line 352 of file level_set.cpp.

◆ MatSideArray

using LevelSet::MatSideArray = std::array<MatrixDouble, 2>
private

Definition at line 80 of file level_set.cpp.

◆ OpMassLL

Examples
level_set.cpp.

Definition at line 354 of file level_set.cpp.

◆ OpMassVV

Examples
level_set.cpp.

Definition at line 358 of file level_set.cpp.

◆ OpScalarFieldL

using LevelSet::OpScalarFieldL = FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm< G>::OpBaseTimesVector<1, DIM1 * DIM2, 1>
private
Examples
level_set.cpp.

Definition at line 362 of file level_set.cpp.

◆ OpSourceL

using LevelSet::OpSourceL = FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm< G>::OpSource<1, DIM1 * DIM2>
private
Examples
level_set.cpp.

Definition at line 356 of file level_set.cpp.

◆ OpSourceV

Examples
level_set.cpp.

Definition at line 360 of file level_set.cpp.

Member Enumeration Documentation

◆ ElementSide

enum LevelSet::ElementSide
private
Enumerator
LEFT_SIDE 
RIGHT_SIDE 

Definition at line 367 of file level_set.cpp.

367 { LEFT_SIDE = 0, RIGHT_SIDE = 1 };

Constructor & Destructor Documentation

◆ LevelSet()

LevelSet::LevelSet ( MoFEM::Interface m_field)
inline

Definition at line 75 of file level_set.cpp.

75 : mField(m_field) {}

Member Function Documentation

◆ dgProjection()

MoFEMErrorCode LevelSet::dgProjection ( const int  prj_bit = projection_bit)
private

dg level set projection

Parameters
prj_bit
mesh_bit
Returns
MoFEMErrorCode
Examples
level_set.cpp.

Definition at line 2401 of file level_set.cpp.

2401  {
2403 
2404  // get operators tester
2405  auto simple = mField.getInterface<Simple>();
2406  auto pip = mField.getInterface<PipelineManager>(); // get interface to
2407  auto bit_mng = mField.getInterface<BitRefManager>();
2408  auto prb_mng = mField.getInterface<ProblemsManager>();
2409 
2410  auto lhs_fe = boost::make_shared<DomainEle>(mField);
2411  auto rhs_fe_prj = boost::make_shared<DomainEle>(mField);
2412  auto rhs_fe_current = boost::make_shared<DomainEle>(mField);
2413 
2414  lhs_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
2415  rhs_fe_prj->getRuleHook = [](int, int, int o) { return 3 * o; };
2416  rhs_fe_current->getRuleHook = [](int, int, int o) { return 3 * o; };
2417 
2418  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
2419  CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "DG_PROJECTION");
2420  CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
2421  CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
2422  CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
2423  CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
2424  CHKERR DMSetUp(sub_dm);
2425 
2426  Range current_ents; // ents used to do calculations
2427  CHKERR bit_mng->getEntitiesByDimAndRefLevel(BitRefLevel().set(current_bit),
2428  BitRefLevel().set(), FE_DIM,
2429  current_ents);
2430  Range prj_ents; // ents from which data are projected
2431  CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2432  BitRefLevel().set(projection_bit), BitRefLevel().set(), FE_DIM, prj_ents);
2433  for (auto l = 0; l != nb_levels; ++l) {
2434  CHKERR bit_mng->updateRangeByParent(prj_ents, prj_ents);
2435  }
2436  current_ents = subtract(
2437  current_ents, prj_ents); // only restric to entities needed projection
2438 
2439  auto test_mesh_bit = [&](FEMethod *fe_ptr) {
2440  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2441  current_bit);
2442  };
2443  auto test_prj_bit = [&](FEMethod *fe_ptr) {
2444  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2445  projection_bit);
2446  };
2447  auto test_current_bit = [&](FEMethod *fe_ptr) {
2448  return current_ents.find(fe_ptr->getFEEntityHandle()) != current_ents.end();
2449  };
2450 
2451  lhs_fe->exeTestHook =
2452  test_mesh_bit; // that element only is run when current bit is set
2453  rhs_fe_prj->exeTestHook =
2454  test_prj_bit; // that element is run only when projection bit is set
2455  rhs_fe_current->exeTestHook =
2456  test_current_bit; // that element is only run when current bit is set
2457 
2458  BitRefLevel remove_mask = BitRefLevel().set(current_bit);
2459  remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
2460  CHKERR prb_mng->removeDofsOnEntities(
2461  "DG_PROJECTION", "L", BitRefLevel().set(), remove_mask, nullptr, 0,
2462  MAX_DOFS_ON_ENTITY, 0, 100, NOISY,
2463  true); // remove all DOFs which are not
2464  // on current bit. This case works for L2 space
2465 
2466  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(lhs_fe->getOpPtrVector(),
2467  {L2});
2468  lhs_fe->getOpPtrVector().push_back(
2469  new OpMassLL("L", "L")); // Assemble projection matrix
2470 
2471  // This assumes that projection mesh is finer, current mesh is coarsened.
2472  auto set_prj_from_child = [&](auto rhs_fe_prj) {
2474  rhs_fe_prj->getOpPtrVector(), {L2});
2475 
2476  // Evaluate field value on projection mesh
2477  auto l_vec = boost::make_shared<MatrixDouble>();
2478  rhs_fe_prj->getOpPtrVector().push_back(
2480 
2481  // This element is used to assemble
2482  auto get_parent_this = [&]() {
2483  auto fe_parent_this = boost::make_shared<DomianParentEle>(mField);
2484  fe_parent_this->getOpPtrVector().push_back(
2485  new OpScalarFieldL("L", l_vec));
2486  return fe_parent_this;
2487  };
2488 
2489  // Create levels of parent elements, until current element is reached, and
2490  // then assemble.
2491  auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
2492  std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
2493  for (int l = 0; l <= nb_levels; ++l)
2494  parents_elems_ptr_vec.emplace_back(
2495  boost::make_shared<DomianParentEle>(mField));
2496  for (auto l = 1; l <= nb_levels; ++l) {
2497  parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
2498  new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
2499  BitRefLevel().set(current_bit).flip(), this_fe_ptr,
2500  BitRefLevel().set(current_bit),
2501  BitRefLevel().set()));
2502  }
2503  return parents_elems_ptr_vec[0];
2504  };
2505 
2506  auto this_fe_ptr = get_parent_this();
2507  auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
2508  rhs_fe_prj->getOpPtrVector().push_back(
2509  new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
2510  BitRefLevel().set(current_bit).flip(), this_fe_ptr,
2511  BitRefLevel().set(current_bit), BitRefLevel().set()));
2512  };
2513 
2514  // This assumed that current mesh is refined, and projection mesh is coarser
2515  auto set_prj_from_parent = [&](auto rhs_fe_current) {
2516 
2517  // Evaluate field value on projection mesh
2518  auto l_vec = boost::make_shared<MatrixDouble>();
2519 
2520  // Evaluate field on coarser element
2521  auto get_parent_this = [&]() {
2522  auto fe_parent_this = boost::make_shared<DomianParentEle>(mField);
2523  fe_parent_this->getOpPtrVector().push_back(
2525  return fe_parent_this;
2526  };
2527 
2528  // Create stack of evaluation on parent elements
2529  auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
2530  std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
2531  for (int l = 0; l <= nb_levels; ++l)
2532  parents_elems_ptr_vec.emplace_back(
2533  boost::make_shared<DomianParentEle>(mField));
2534  for (auto l = 1; l <= nb_levels; ++l) {
2535  parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
2536  new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
2537  BitRefLevel().set(projection_bit).flip(),
2538  this_fe_ptr, BitRefLevel().set(projection_bit),
2539  BitRefLevel().set()));
2540  }
2541  return parents_elems_ptr_vec[0];
2542  };
2543 
2544  auto this_fe_ptr = get_parent_this();
2545  auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
2546 
2547  auto reset_op_ptr = new DomainEleOp(NOSPACE, DomainEleOp::OPSPACE);
2548  reset_op_ptr->doWorkRhsHook = [&](DataOperator *op_ptr, int, EntityType,
2549  EntData &) {
2550  l_vec->resize(static_cast<DomainEleOp *>(op_ptr)->getGaussPts().size2(),
2551  false);
2552  l_vec->clear();
2553  return 0;
2554  };
2555  rhs_fe_current->getOpPtrVector().push_back(reset_op_ptr);
2556  rhs_fe_current->getOpPtrVector().push_back(
2557  new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
2558  BitRefLevel().set(projection_bit).flip(), this_fe_ptr,
2559  BitRefLevel(), BitRefLevel()));
2560 
2561  // At the end assemble of current finite element
2562  rhs_fe_current->getOpPtrVector().push_back(new OpScalarFieldL("L", l_vec));
2563  };
2564 
2565  set_prj_from_child(rhs_fe_prj);
2566  set_prj_from_parent(rhs_fe_current);
2567 
2568  boost::shared_ptr<FEMethod> null_fe;
2569  getDMKspCtx(sub_dm)->clearLoops();
2570  CHKERR DMMoFEMKSPSetComputeOperators(sub_dm, simple->getDomainFEName(),
2571  lhs_fe, null_fe, null_fe);
2572  CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(), rhs_fe_prj,
2573  null_fe, null_fe);
2574  CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(),
2575  rhs_fe_current, null_fe, null_fe);
2576  auto ksp = MoFEM::createKSP(mField.get_comm());
2577  CHKERR KSPSetDM(ksp, sub_dm);
2578 
2579  CHKERR KSPSetDM(ksp, sub_dm);
2580  CHKERR KSPSetFromOptions(ksp);
2581  CHKERR KSPSetUp(ksp);
2582 
2583  auto L = createDMVector(sub_dm);
2584  auto F = vectorDuplicate(L);
2585 
2586  CHKERR KSPSolve(ksp, F, L);
2587  CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
2588  CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
2589  CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
2590 
2591  auto [error, th_error] = evaluateError();
2592  MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
2593 
2594  auto post_proc = [&](auto dm, auto out_name, auto th_error) {
2596  auto post_proc_fe =
2597  boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
2598  post_proc_fe->setTagsToTransfer({th_error});
2599  post_proc_fe->exeTestHook = test_mesh_bit;
2600 
2601  if constexpr (DIM1 == 1 && DIM2 == 1) {
2603  auto l_vec = boost::make_shared<VectorDouble>();
2604  auto l_grad_mat = boost::make_shared<MatrixDouble>();
2606  post_proc_fe->getOpPtrVector(), {L2});
2607  post_proc_fe->getOpPtrVector().push_back(
2608  new OpCalculateScalarFieldValues("L", l_vec));
2609  post_proc_fe->getOpPtrVector().push_back(
2610  new OpCalculateScalarFieldGradient<SPACE_DIM>("L", l_grad_mat));
2611 
2612  post_proc_fe->getOpPtrVector().push_back(
2613 
2614  new OpPPMap(
2615 
2616  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2617 
2618  {{"L", l_vec}},
2619 
2620  {{"GradL", l_grad_mat}},
2621 
2622  {}, {})
2623 
2624  );
2625  }
2626 
2627  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
2628  post_proc_fe);
2629  post_proc_fe->writeFile(out_name);
2631  };
2632 
2633  if constexpr (debug)
2634  CHKERR post_proc(sub_dm, "dg_projection.h5m", th_error);
2635 
2637 }

◆ 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.

1203  {
1204 
1205  struct OpErrorSkel : BoundaryEleOp {
1206 
1207  OpErrorSkel(boost::shared_ptr<FaceSideEle> side_fe_ptr,
1208  boost::shared_ptr<SideData> side_data_ptr,
1209  SmartPetscObj<Vec> error_sum_ptr, Tag th_error)
1210  : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
1211  sideFEPtr(side_fe_ptr), sideDataPtr(side_data_ptr),
1212  errorSumPtr(error_sum_ptr), thError(th_error) {}
1213 
1214  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1216 
1217  // Collect data from side domain elements
1218  CHKERR loopSideFaces("dFE", *sideFEPtr);
1219  const auto in_the_loop =
1220  sideFEPtr->nInTheLoop; // return number of elements on the side
1221 
1222  auto not_side = [](auto s) {
1223  return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
1224  };
1225 
1226  auto nb_gauss_pts = getGaussPts().size2();
1227 
1228  for (auto s : {LEFT_SIDE, RIGHT_SIDE}) {
1229 
1230  auto arr_t_l = make_array(
1231  getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[LEFT_SIDE]),
1232  getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[RIGHT_SIDE]));
1233  auto arr_t_vel = make_array(
1234  getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[LEFT_SIDE]),
1235  getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[RIGHT_SIDE]));
1236 
1237  auto next = [&]() {
1238  for (auto &t_l : arr_t_l)
1239  ++t_l;
1240  for (auto &t_vel : arr_t_vel)
1241  ++t_vel;
1242  };
1243 
1244  double e = 0;
1245  auto t_w = getFTensor0IntegrationWeight();
1246  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1248  t_diff(I, J) = arr_t_l[LEFT_SIDE](I, J) - arr_t_l[RIGHT_SIDE](I, J);
1249  e += t_w * getMeasure() * t_diff(I, J) * t_diff(I, J);
1250  next();
1251  ++t_w;
1252  }
1253  e = std::sqrt(e);
1254 
1255  moab::Interface &moab =
1256  getNumeredEntFiniteElementPtr()->getBasicDataPtr()->moab;
1257  const void *tags_ptr[2];
1258  CHKERR moab.tag_get_by_ptr(thError, sideDataPtr->feSideHandle.data(), 2,
1259  tags_ptr);
1260  for (auto ff : {0, 1}) {
1261  *((double *)tags_ptr[ff]) += e;
1262  }
1263  CHKERR VecSetValue(errorSumPtr, 0, e, ADD_VALUES);
1264  };
1265 
1267  }
1268 
1269  private:
1270  boost::shared_ptr<FaceSideEle> sideFEPtr;
1271  boost::shared_ptr<SideData> sideDataPtr;
1272  SmartPetscObj<Vec> errorSumPtr;
1273  Tag thError;
1274  };
1275 
1276  auto simple = mField.getInterface<Simple>();
1277 
1278  auto error_sum_ptr = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
1279  Tag th_error;
1280  double def_val = 0;
1281  CHKERR mField.get_moab().tag_get_handle("Error", 1, MB_TYPE_DOUBLE, th_error,
1282  MB_TAG_CREAT | MB_TAG_SPARSE,
1283  &def_val);
1284 
1285  auto clear_tags = [&]() {
1287  Range fe_ents;
1288  CHKERR mField.get_moab().get_entities_by_dimension(0, FE_DIM, fe_ents);
1289  double zero;
1290  CHKERR mField.get_moab().tag_clear_data(th_error, fe_ents, &zero);
1292  };
1293 
1294  auto evaluate_error = [&]() {
1296  auto skel_fe = boost::make_shared<BoundaryEle>(mField);
1297  skel_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
1298  auto side_data_ptr = boost::make_shared<SideData>();
1299  auto side_fe_ptr = getSideFE(side_data_ptr);
1300  skel_fe->getOpPtrVector().push_back(
1301  new OpErrorSkel(side_fe_ptr, side_data_ptr, error_sum_ptr, th_error));
1302  auto simple = mField.getInterface<Simple>();
1303 
1304  skel_fe->exeTestHook = [&](FEMethod *fe_ptr) {
1305  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1306  skeleton_bit);
1307  };
1308 
1310  simple->getSkeletonFEName(), skel_fe);
1311 
1313  };
1314 
1315  auto assemble_and_sum = [](auto vec) {
1316  CHK_THROW_MESSAGE(VecAssemblyBegin(vec), "assemble");
1317  CHK_THROW_MESSAGE(VecAssemblyEnd(vec), "assemble");
1318  double sum;
1319  CHK_THROW_MESSAGE(VecSum(vec, &sum), "assemble");
1320  return sum;
1321  };
1322 
1323  auto propagate_error_to_parents = [&]() {
1325 
1326  auto &moab = mField.get_moab();
1327  auto fe_ptr = boost::make_shared<FEMethod>();
1328  fe_ptr->exeTestHook = [&](FEMethod *fe_ptr) {
1329  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1330  current_bit);
1331  };
1332 
1333  fe_ptr->preProcessHook = []() { return 0; };
1334  fe_ptr->postProcessHook = []() { return 0; };
1335  fe_ptr->operatorHook = [&]() {
1337 
1338  auto fe_ent = fe_ptr->numeredEntFiniteElementPtr->getEnt();
1339  auto parent = fe_ptr->numeredEntFiniteElementPtr->getParentEnt();
1340  auto th_parent = fe_ptr->numeredEntFiniteElementPtr->getBasicDataPtr()
1341  ->th_RefParentHandle;
1342 
1343  double error;
1344  CHKERR moab.tag_get_data(th_error, &fe_ent, 1, &error);
1345 
1346  boost::function<MoFEMErrorCode(EntityHandle, double)> add_error =
1347  [&](auto fe_ent, auto error) {
1349  double *e_ptr;
1350  CHKERR moab.tag_get_by_ptr(th_error, &fe_ent, 1,
1351  (const void **)&e_ptr);
1352  (*e_ptr) += error;
1353 
1354  EntityHandle parent;
1355  CHKERR moab.tag_get_data(th_parent, &fe_ent, 1, &parent);
1356  if (parent != fe_ent && parent)
1357  CHKERR add_error(parent, *e_ptr);
1358 
1360  };
1361 
1362  CHKERR add_error(parent, error);
1363 
1365  };
1366 
1367  CHKERR DMoFEMLoopFiniteElements(simple->getDM(), simple->getDomainFEName(),
1368  fe_ptr);
1369 
1371  };
1372 
1373  CHK_THROW_MESSAGE(clear_tags(), "clear error tags");
1374  CHK_THROW_MESSAGE(evaluate_error(), "evaluate error");
1375  CHK_THROW_MESSAGE(propagate_error_to_parents(), "propagate error");
1376 
1377  return std::make_tuple(assemble_and_sum(error_sum_ptr), th_error);
1378 }

◆ get_level_set()

double LevelSet::get_level_set ( const double  x,
const double  y,
const double  z 
)
staticprivate

inital level set, i.e. advected filed

Parameters
x
y
z
Returns
double
Examples
level_set.cpp.

Definition at line 378 of file level_set.cpp.

378  {
379  constexpr double xc = 0.1;
380  constexpr double yc = 0.;
381  constexpr double zc = 0.;
382  constexpr double r = 0.2;
383  return std::sqrt(pow(x - xc, 2) + pow(y - yc, 2) + pow(z - zc, 2)) - r;
384 }

◆ get_velocity_potential() [1/2]

template<int FE_DIM>
static double LevelSet::get_velocity_potential ( double  x,
double  y,
double  z 
)
staticprivate

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
FE_DIM
Parameters
x
y
z
Returns
auto

◆ get_velocity_potential() [2/2]

template<>
double LevelSet::get_velocity_potential ( double  x,
double  y,
double  z 
)
staticprivate

Definition at line 374 of file level_set.cpp.

374  {
375  return (x * x - 0.25) * (y * y - 0.25);
376 }

◆ getSideFE()

boost::shared_ptr< FaceSideEle > LevelSet::getSideFE ( boost::shared_ptr< SideData side_data_ptr)
private

create side element to assemble data from sides

Parameters
side_data_ptr
Returns
boost::shared_ptr<FaceSideEle>
Examples
level_set.cpp.

Definition at line 997 of file level_set.cpp.

997  {
998 
999  auto simple = mField.getInterface<Simple>();
1000 
1001  auto l_ptr = boost::make_shared<MatrixDouble>();
1002  auto vel_ptr = boost::make_shared<MatrixDouble>();
1003 
1004  struct OpSideData : public FaceSideEleOp {
1005  OpSideData(boost::shared_ptr<SideData> side_data_ptr)
1006  : FaceSideEleOp("L", "L", FaceSideEleOp::OPROWCOL),
1007  sideDataPtr(side_data_ptr) {
1008  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
1009  for (auto t = moab::CN::TypeDimensionMap[FE_DIM].first;
1010  t <= moab::CN::TypeDimensionMap[FE_DIM].second; ++t)
1011  doEntities[t] = true;
1012  sYmm = false;
1013  }
1014 
1015  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
1016  EntityType col_type, EntData &row_data,
1017  EntData &col_data) {
1019  if ((CN::Dimension(row_type) == FE_DIM) &&
1020  (CN::Dimension(col_type) == FE_DIM)) {
1021 
1022  auto reset = [&](auto nb_in_loop) {
1023  sideDataPtr->feSideHandle[nb_in_loop] = 0;
1024  sideDataPtr->indicesRowSideMap[nb_in_loop].clear();
1025  sideDataPtr->indicesColSideMap[nb_in_loop].clear();
1026  sideDataPtr->rowBaseSideMap[nb_in_loop].clear();
1027  sideDataPtr->colBaseSideMap[nb_in_loop].clear();
1028  sideDataPtr->senseMap[nb_in_loop] = 0;
1029  };
1030 
1031  const auto nb_in_loop = getFEMethod()->nInTheLoop;
1032  if (nb_in_loop == 0)
1033  for (auto s : {0, 1})
1034  reset(s);
1035 
1036  sideDataPtr->currentFESide = nb_in_loop;
1037  sideDataPtr->senseMap[nb_in_loop] = getSkeletonSense();
1038 
1039  } else {
1040  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Should not happen");
1041  }
1042 
1044  };
1045 
1046  private:
1047  boost::shared_ptr<SideData> sideDataPtr;
1048  };
1049 
1050  struct OpSideDataOnParent : public DomainEleOp {
1051 
1052  OpSideDataOnParent(boost::shared_ptr<SideData> side_data_ptr,
1053  boost::shared_ptr<MatrixDouble> l_ptr,
1054  boost::shared_ptr<MatrixDouble> vel_ptr)
1055  : DomainEleOp("L", "L", DomainEleOp::OPROWCOL),
1056  sideDataPtr(side_data_ptr), lPtr(l_ptr), velPtr(vel_ptr) {
1057  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
1058  for (auto t = moab::CN::TypeDimensionMap[FE_DIM].first;
1059  t <= moab::CN::TypeDimensionMap[FE_DIM].second; ++t)
1060  doEntities[t] = true;
1061  sYmm = false;
1062  }
1063 
1064  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
1065  EntityType col_type, EntData &row_data,
1066  EntData &col_data) {
1068 
1069  if ((CN::Dimension(row_type) == FE_DIM) &&
1070  (CN::Dimension(col_type) == FE_DIM)) {
1071  const auto nb_in_loop = sideDataPtr->currentFESide;
1072  sideDataPtr->feSideHandle[nb_in_loop] = getFEEntityHandle();
1073  sideDataPtr->indicesRowSideMap[nb_in_loop] = row_data.getIndices();
1074  sideDataPtr->indicesColSideMap[nb_in_loop] = col_data.getIndices();
1075  sideDataPtr->rowBaseSideMap[nb_in_loop] = row_data.getN();
1076  sideDataPtr->colBaseSideMap[nb_in_loop] = col_data.getN();
1077  (sideDataPtr->lVec)[nb_in_loop] = *lPtr;
1078  (sideDataPtr->velMat)[nb_in_loop] = *velPtr;
1079 
1080 #ifndef NDEBUG
1081  if ((sideDataPtr->lVec)[nb_in_loop].size2() !=
1082  (sideDataPtr->velMat)[nb_in_loop].size2())
1083  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1084  "Wrong number of integaration pts %d != %d",
1085  (sideDataPtr->lVec)[nb_in_loop].size2(),
1086  (sideDataPtr->velMat)[nb_in_loop].size2());
1087  if ((sideDataPtr->velMat)[nb_in_loop].size1() != SPACE_DIM)
1088  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1089  "Wrong size of velocity vector size = %d",
1090  (sideDataPtr->velMat)[nb_in_loop].size1());
1091 #endif
1092 
1093  if (!nb_in_loop) {
1094  (sideDataPtr->lVec)[1] = sideDataPtr->lVec[0];
1095  (sideDataPtr->velMat)[1] = (sideDataPtr->velMat)[0];
1096  } else {
1097 #ifndef NDEBUG
1098  if (sideDataPtr->rowBaseSideMap[0].size1() !=
1099  sideDataPtr->rowBaseSideMap[1].size1()) {
1100  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1101  "Wrong number of integration pt %d != %d",
1102  sideDataPtr->rowBaseSideMap[0].size1(),
1103  sideDataPtr->rowBaseSideMap[1].size1());
1104  }
1105  if (sideDataPtr->colBaseSideMap[0].size1() !=
1106  sideDataPtr->colBaseSideMap[1].size1()) {
1107  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1108  "Wrong number of integration pt");
1109  }
1110 #endif
1111  }
1112 
1113  } else {
1114  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Should not happen");
1115  }
1116 
1118  };
1119 
1120  private:
1121  boost::shared_ptr<SideData> sideDataPtr;
1122  boost::shared_ptr<MatrixDouble> lPtr;
1123  boost::shared_ptr<MatrixDouble> velPtr;
1124  };
1125 
1126  // Calculate fields on param mesh bit element
1127  auto get_parent_this = [&]() {
1128  auto parent_fe_ptr = boost::make_shared<DomianParentEle>(mField);
1130  parent_fe_ptr->getOpPtrVector(), {L2});
1131  parent_fe_ptr->getOpPtrVector().push_back(
1133  parent_fe_ptr->getOpPtrVector().push_back(
1134  new OpSideDataOnParent(side_data_ptr, l_ptr, vel_ptr));
1135  return parent_fe_ptr;
1136  };
1137 
1138  auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
1139  std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
1140  for (int l = 0; l <= nb_levels; ++l)
1141  parents_elems_ptr_vec.emplace_back(
1142  boost::make_shared<DomianParentEle>(mField));
1143  for (auto l = 1; l <= nb_levels; ++l) {
1144  parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
1145  new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
1146  BitRefLevel().set(current_bit).flip(), this_fe_ptr,
1147  BitRefLevel().set(current_bit), BitRefLevel().set()));
1148  }
1149  return parents_elems_ptr_vec[0];
1150  };
1151 
1152  // Create aliased shared pointers, all elements are destroyed if side_fe_ptr
1153  // is destroyed
1154  auto get_side_fe_ptr = [&]() {
1155  auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
1156 
1157  auto this_fe_ptr = get_parent_this();
1158  auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
1159 
1160  side_fe_ptr->getOpPtrVector().push_back(new OpSideData(side_data_ptr));
1161  side_fe_ptr->getOpPtrVector().push_back(getZeroLevelVelOp(vel_ptr));
1162  side_fe_ptr->getOpPtrVector().push_back(
1163  new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
1164  BitRefLevel().set(current_bit).flip(), this_fe_ptr,
1165  BitRefLevel().set(current_bit), BitRefLevel().set()));
1166 
1167  return side_fe_ptr;
1168  };
1169 
1170  return get_side_fe_ptr();
1171 };

◆ getZeroLevelVelOp()

ForcesAndSourcesCore::UserDataOperator * LevelSet::getZeroLevelVelOp ( boost::shared_ptr< MatrixDouble >  vel_ptr)
private

Get operator calculating velocity on coarse mesh.

Parameters
vel_ptr
Returns
DomainEleOp*
Examples
level_set.cpp.

Definition at line 922 of file level_set.cpp.

922  {
923  auto get_parent_vel_this = [&]() {
924  auto parent_fe_ptr = boost::make_shared<DomianParentEle>(mField);
926  parent_fe_ptr->getOpPtrVector(), {potential_velocity_space});
927  parent_fe_ptr->getOpPtrVector().push_back(
929  "V", vel_ptr));
930  return parent_fe_ptr;
931  };
932 
933  auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr) {
934  std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
935  for (int l = 0; l <= nb_levels; ++l)
936  parents_elems_ptr_vec.emplace_back(
937  boost::make_shared<DomianParentEle>(mField));
938  for (auto l = 1; l <= nb_levels; ++l) {
939  parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
940  new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
941  BitRefLevel().set(0).flip(), this_fe_ptr,
942  BitRefLevel().set(0), BitRefLevel().set()));
943  }
944  return parents_elems_ptr_vec[0];
945  };
946 
947  auto this_fe_ptr = get_parent_vel_this();
948  auto parent_fe_ptr = get_parents_vel_fe_ptr(this_fe_ptr);
949  return new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
950  BitRefLevel().set(0).flip(), this_fe_ptr,
951  BitRefLevel().set(0), BitRefLevel().set());
952 }

◆ initialiseFieldLevelSet()

MoFEMErrorCode LevelSet::initialiseFieldLevelSet ( boost::function< double(double, double, double)>  level_fun = get_level_set)
private

initialise field set

Parameters
level_fun
Returns
MoFEMErrorCode
Examples
level_set.cpp.

Definition at line 1693 of file level_set.cpp.

1694  {
1696 
1697  // get operators tester
1698  auto simple = mField.getInterface<Simple>();
1699  auto pip = mField.getInterface<PipelineManager>(); // get interface to
1700  // pipeline manager
1701  auto prb_mng = mField.getInterface<ProblemsManager>();
1702 
1703  boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(mField);
1704  boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(mField);
1705  auto swap_fe = [&]() {
1706  lhs_fe.swap(pip->getDomainLhsFE());
1707  rhs_fe.swap(pip->getDomainRhsFE());
1708  };
1709  swap_fe();
1710 
1711  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 3 * o; });
1712  pip->setDomainRhsIntegrationRule([](int, int, int o) { return 3 * o; });
1713 
1714  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
1715  CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "LEVELSET_POJECTION");
1716  CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
1717  CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
1718  CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
1719  CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
1720  CHKERR DMSetUp(sub_dm);
1721 
1722  BitRefLevel remove_mask = BitRefLevel().set(current_bit);
1723  remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
1724  CHKERR prb_mng->removeDofsOnEntities("LEVELSET_POJECTION", "L",
1725  BitRefLevel().set(), remove_mask);
1726  auto test_bit_ele = [&](FEMethod *fe_ptr) {
1727  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1728  current_bit);
1729  };
1730  pip->getDomainLhsFE()->exeTestHook = test_bit_ele;
1731  pip->getDomainRhsFE()->exeTestHook = test_bit_ele;
1732 
1733  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainRhsPipeline(),
1734  {L2});
1735  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainLhsPipeline(),
1736  {L2});
1737  pip->getOpDomainLhsPipeline().push_back(new OpMassLL("L", "L"));
1738  pip->getOpDomainRhsPipeline().push_back(new OpSourceL("L", level_fun));
1739 
1740  CHKERR mField.getInterface<FieldBlas>()->setField(0, "L");
1741 
1742  auto ksp = pip->createKSP(sub_dm);
1743  CHKERR KSPSetDM(ksp, sub_dm);
1744  CHKERR KSPSetFromOptions(ksp);
1745  CHKERR KSPSetUp(ksp);
1746 
1747  auto L = createDMVector(sub_dm);
1748  auto F = vectorDuplicate(L);
1749 
1750  CHKERR KSPSolve(ksp, F, L);
1751  CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
1752  CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
1753  CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
1754 
1755  auto [error, th_error] = evaluateError();
1756  MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
1757 #ifndef NDEBUG
1758  auto fe_meshset =
1759  mField.get_finite_element_meshset(simple->getDomainFEName());
1760  std::vector<Tag> tags{th_error};
1761  CHKERR mField.get_moab().write_file("error.h5m", "MOAB",
1762  "PARALLEL=WRITE_PART", &fe_meshset, 1,
1763  &*tags.begin(), tags.size());
1764 #endif
1765 
1766  auto post_proc = [&](auto dm, auto out_name, auto th_error) {
1768  auto post_proc_fe =
1769  boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
1770  post_proc_fe->setTagsToTransfer({th_error});
1771  post_proc_fe->exeTestHook = test_bit_ele;
1772 
1773  if constexpr (DIM1 == 1 && DIM2 == 1) {
1775 
1776  auto l_vec = boost::make_shared<VectorDouble>();
1777  auto l_grad_mat = boost::make_shared<MatrixDouble>();
1779  post_proc_fe->getOpPtrVector(), {L2});
1780  post_proc_fe->getOpPtrVector().push_back(
1781  new OpCalculateScalarFieldValues("L", l_vec));
1782  post_proc_fe->getOpPtrVector().push_back(
1783  new OpCalculateScalarFieldGradient<SPACE_DIM>("L", l_grad_mat));
1784 
1785  post_proc_fe->getOpPtrVector().push_back(
1786 
1787  new OpPPMap(
1788 
1789  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1790 
1791  {{"L", l_vec}},
1792 
1793  {{"GradL", l_grad_mat}},
1794 
1795  {}, {})
1796 
1797  );
1798  }
1799 
1800  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
1801  post_proc_fe);
1802  post_proc_fe->writeFile(out_name);
1804  };
1805 
1806  if constexpr (debug)
1807  CHKERR post_proc(sub_dm, "initial_level_set.h5m", th_error);
1808 
1809  swap_fe();
1810 
1812 }

◆ initialiseFieldVelocity()

MoFEMErrorCode LevelSet::initialiseFieldVelocity ( boost::function< double(double, double, double)>  vel_fun = get_velocity_potential<FE_DIM>)
private

initialise potential velocity field

Parameters
vel_fun
Returns
MoFEMErrorCode
Examples
level_set.cpp.

Definition at line 1814 of file level_set.cpp.

1815  {
1817 
1818  // get operators tester
1819  auto simple = mField.getInterface<Simple>();
1820  auto pip = mField.getInterface<PipelineManager>(); // get interface to
1821  // pipeline manager
1822  auto prb_mng = mField.getInterface<ProblemsManager>();
1823 
1824  boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(mField);
1825  boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(mField);
1826  auto swap_fe = [&]() {
1827  lhs_fe.swap(pip->getDomainLhsFE());
1828  rhs_fe.swap(pip->getDomainRhsFE());
1829  };
1830  swap_fe();
1831 
1832  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 3 * o; });
1833  pip->setDomainRhsIntegrationRule([](int, int, int o) { return 3 * o; });
1834 
1835  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
1836  CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "VELOCITY_PROJECTION");
1837  CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
1838  CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
1839  CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
1840  CHKERR DMMoFEMAddSubFieldRow(sub_dm, "V");
1841  CHKERR DMSetUp(sub_dm);
1842 
1843  // Velocities are calculated only on corse mesh
1844  BitRefLevel remove_mask = BitRefLevel().set(0);
1845  remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
1846  CHKERR prb_mng->removeDofsOnEntities("VELOCITY_PROJECTION", "V",
1847  BitRefLevel().set(), remove_mask);
1848 
1849  auto test_bit = [&](FEMethod *fe_ptr) {
1850  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(0);
1851  };
1852  pip->getDomainLhsFE()->exeTestHook = test_bit;
1853  pip->getDomainRhsFE()->exeTestHook = test_bit;
1854 
1855  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainLhsPipeline(),
1856  {potential_velocity_space});
1857  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainRhsPipeline(),
1858  {potential_velocity_space});
1859 
1860  pip->getOpDomainLhsPipeline().push_back(new OpMassVV("V", "V"));
1861  pip->getOpDomainRhsPipeline().push_back(new OpSourceV("V", vel_fun));
1862 
1863  auto ksp = pip->createKSP(sub_dm);
1864  CHKERR KSPSetDM(ksp, sub_dm);
1865  CHKERR KSPSetFromOptions(ksp);
1866  CHKERR KSPSetUp(ksp);
1867 
1868  auto L = createDMVector(sub_dm);
1869  auto F = vectorDuplicate(L);
1870 
1871  CHKERR KSPSolve(ksp, F, L);
1872  CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
1873  CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
1874  CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
1875 
1876  auto post_proc = [&](auto dm, auto out_name) {
1878  auto post_proc_fe =
1879  boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
1880  post_proc_fe->exeTestHook = test_bit;
1881 
1882  if constexpr (FE_DIM == 2) {
1883 
1885  post_proc_fe->getOpPtrVector(), {potential_velocity_space});
1886 
1888 
1889  auto potential_vec = boost::make_shared<VectorDouble>();
1890  auto velocity_mat = boost::make_shared<MatrixDouble>();
1891 
1892  post_proc_fe->getOpPtrVector().push_back(
1893  new OpCalculateScalarFieldValues("V", potential_vec));
1894  post_proc_fe->getOpPtrVector().push_back(
1896  SPACE_DIM>("V", velocity_mat));
1897 
1898  post_proc_fe->getOpPtrVector().push_back(
1899 
1900  new OpPPMap(
1901 
1902  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1903 
1904  {{"VelocityPotential", potential_vec}},
1905 
1906  {{"Velocity", velocity_mat}},
1907 
1908  {}, {})
1909 
1910  );
1911 
1912  } else {
1913  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1914  "3d case not implemented");
1915  }
1916 
1917  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
1918  post_proc_fe);
1919  post_proc_fe->writeFile(out_name);
1921  };
1922 
1923  if constexpr (debug)
1924  CHKERR post_proc(sub_dm, "initial_velocity_potential.h5m");
1925 
1926  swap_fe();
1927 
1929 }

◆ pushOpDomain()

MoFEMErrorCode LevelSet::pushOpDomain ( )
private

push operators to integrate operators on domain

Returns
MoFEMErrorCode
Examples
level_set.cpp.

Definition at line 954 of file level_set.cpp.

954  {
956  auto pip = mField.getInterface<PipelineManager>(); // get interface to
957  // pipeline manager
958 
959  pip->getOpDomainLhsPipeline().clear();
960  pip->getOpDomainRhsPipeline().clear();
961 
962  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 3 * o; });
963  pip->setDomainRhsIntegrationRule([](int, int, int o) { return 3 * o; });
964 
965  pip->getDomainLhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
966  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
967  current_bit);
968  };
969  pip->getDomainRhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
970  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
971  current_bit);
972  };
973 
974  auto l_ptr = boost::make_shared<MatrixDouble>();
975  auto l_dot_ptr = boost::make_shared<MatrixDouble>();
976  auto vel_ptr = boost::make_shared<MatrixDouble>();
977 
978  pip->getOpDomainRhsPipeline().push_back(getZeroLevelVelOp(vel_ptr));
979  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainRhsPipeline(),
980  {L2});
981  pip->getOpDomainRhsPipeline().push_back(
983  pip->getOpDomainRhsPipeline().push_back(
985  pip->getOpDomainRhsPipeline().push_back(
986  new OpRhsDomain("L", l_ptr, l_dot_ptr, vel_ptr));
987 
988  pip->getOpDomainLhsPipeline().push_back(getZeroLevelVelOp(vel_ptr));
989  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainLhsPipeline(),
990  {L2});
991  pip->getOpDomainLhsPipeline().push_back(new OpLhsDomain("L", vel_ptr));
992 
994 }

◆ pushOpSkeleton()

MoFEMErrorCode LevelSet::pushOpSkeleton ( )
private

push operator to integrate on skeleton

Returns
MoFEMErrorCode
Examples
level_set.cpp.

Definition at line 1173 of file level_set.cpp.

1173  {
1175  auto pip = mField.getInterface<PipelineManager>(); // get interface to
1176 
1177  pip->getOpSkeletonLhsPipeline().clear();
1178  pip->getOpSkeletonRhsPipeline().clear();
1179 
1180  pip->setSkeletonLhsIntegrationRule([](int, int, int o) { return 18; });
1181  pip->setSkeletonRhsIntegrationRule([](int, int, int o) { return 18; });
1182 
1183  pip->getSkeletonLhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
1184  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1185  skeleton_bit);
1186  };
1187  pip->getSkeletonRhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
1188  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1189  skeleton_bit);
1190  };
1191 
1192  auto side_data_ptr = boost::make_shared<SideData>();
1193  auto side_fe_ptr = getSideFE(side_data_ptr);
1194 
1195  pip->getOpSkeletonRhsPipeline().push_back(
1196  new OpRhsSkeleton(side_data_ptr, side_fe_ptr));
1197  pip->getOpSkeletonLhsPipeline().push_back(
1198  new OpLhsSkeleton(side_data_ptr, side_fe_ptr));
1199 
1201 }

◆ readMesh()

MoFEMErrorCode LevelSet::readMesh ( )
private

read mesh

Returns
MoFEMErrorCode
Examples
level_set.cpp.

Definition at line 502 of file level_set.cpp.

502  {
504  auto simple = mField.getInterface<Simple>();
505  // get options from command line
506  CHKERR simple->getOptions();
507 
508  // Only L2 field is set in this example. Two lines bellow forces simple
509  // interface to creat lower dimension (edge) elements, despite that fact that
510  // there is no field spanning on such elements. We need them for DG method.
511  simple->getAddSkeletonFE() = true;
512  simple->getAddBoundaryFE() = true;
513 
514  // load mesh file
515  simple->getBitRefLevel() = BitRefLevel();
516  CHKERR simple->loadFile();
517 
518  // Initialise bit ref levels
519  auto set_problem_bit = [&]() {
521  auto bit0 = BitRefLevel().set(start_bit);
522  BitRefLevel start_mask;
523  for (auto s = 0; s != start_bit; ++s)
524  start_mask[s] = true;
525 
526  auto bit_mng = mField.getInterface<BitRefManager>();
527 
528  Range level0;
529  CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(0),
530  BitRefLevel().set(), level0);
531 
532  CHKERR bit_mng->setNthBitRefLevel(level0, current_bit, true);
533  CHKERR bit_mng->setNthBitRefLevel(level0, aggregate_bit, true);
534  CHKERR bit_mng->setNthBitRefLevel(level0, skeleton_bit, true);
535 
536  // Set bits to build adjacencies between parents and children. That is used
537  // by simple interface.
538  simple->getBitAdjEnt() = BitRefLevel().set();
539  simple->getBitAdjParent() = BitRefLevel().set();
540  simple->getBitRefLevel() = BitRefLevel().set(current_bit);
541  simple->getBitRefLevelMask() = BitRefLevel().set();
542 
543 #ifndef NDEBUG
544  if constexpr (debug) {
545  auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
546  CHKERR bit_mng->writeBitLevelByDim(
547  BitRefLevel().set(0), BitRefLevel().set(), FE_DIM,
548  (proc_str + "level_base.vtk").c_str(), "VTK", "");
549  }
550 #endif
551 
553  };
554 
555  CHKERR set_problem_bit();
556 
558 }

◆ refineMesh()

MoFEMErrorCode LevelSet::refineMesh ( WrapperClass &&  wp)
private
Examples
level_set.cpp.

Definition at line 2147 of file level_set.cpp.

2147  {
2149 
2150  auto simple = mField.getInterface<Simple>();
2151  auto bit_mng = mField.getInterface<BitRefManager>();
2152  ParallelComm *pcomm =
2153  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2154  auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
2155 
2156  auto set_bit = [](auto l) { return BitRefLevel().set(l); };
2157 
2158  auto save_range = [&](const std::string name, const Range &r) {
2160  auto meshset_ptr = get_temp_meshset_ptr(mField.get_moab());
2161  CHKERR mField.get_moab().add_entities(*meshset_ptr, r);
2162  CHKERR mField.get_moab().write_file(name.c_str(), "VTK", "",
2163  meshset_ptr->get_ptr(), 1);
2165  };
2166 
2167  // select domain elements to refine by threshold
2168  auto get_refined_elements_meshset = [&](auto bit, auto mask) {
2169  Range fe_ents;
2170  CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit, mask, FE_DIM, fe_ents);
2171 
2172  Tag th_error;
2173  CHK_MOAB_THROW(mField.get_moab().tag_get_handle("Error", th_error),
2174  "get error handle");
2175  std::vector<double> errors(fe_ents.size());
2177  mField.get_moab().tag_get_data(th_error, fe_ents, &*errors.begin()),
2178  "get tag data");
2179  auto it = std::max_element(errors.begin(), errors.end());
2180  double max;
2181  MPI_Allreduce(&*it, &max, 1, MPI_DOUBLE, MPI_MAX, mField.get_comm());
2182  MOFEM_LOG("LevelSet", Sev::inform) << "Max error: " << max;
2183  auto threshold = wp.getThreshold(max);
2184 
2185  std::vector<EntityHandle> fe_to_refine;
2186  fe_to_refine.reserve(fe_ents.size());
2187 
2188  auto fe_it = fe_ents.begin();
2189  auto error_it = errors.begin();
2190  for (auto i = 0; i != fe_ents.size(); ++i) {
2191  if (*error_it > threshold) {
2192  fe_to_refine.push_back(*fe_it);
2193  }
2194  ++fe_it;
2195  ++error_it;
2196  }
2197 
2198  Range ents;
2199  ents.insert_list(fe_to_refine.begin(), fe_to_refine.end());
2200  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2201  ents, nullptr, NOISY);
2202 
2203  auto get_neighbours_by_bridge_vertices = [&](auto &&ents) {
2204  Range verts;
2205  CHKERR mField.get_moab().get_connectivity(ents, verts, true);
2206  CHKERR mField.get_moab().get_adjacencies(verts, FE_DIM, false, ents,
2207  moab::Interface::UNION);
2208  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ents);
2209  return ents;
2210  };
2211 
2212  ents = get_neighbours_by_bridge_vertices(ents);
2213 
2214 #ifndef NDEBUG
2215  if (debug) {
2216  auto meshset_ptr = get_temp_meshset_ptr(mField.get_moab());
2217  CHK_MOAB_THROW(mField.get_moab().add_entities(*meshset_ptr, ents),
2218  "add entities to meshset");
2219  CHKERR mField.get_moab().write_file(
2220  (proc_str + "_fe_to_refine.vtk").c_str(), "VTK", "",
2221  meshset_ptr->get_ptr(), 1);
2222  }
2223 #endif
2224 
2225  return ents;
2226  };
2227 
2228  // refine elements, and set bit ref level
2229  auto refine_mesh = [&](auto l, auto &&fe_to_refine) {
2230  Skinner skin(&mField.get_moab());
2232 
2233  // get entities in "l-1" level
2234  Range level_ents;
2235  CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2236  set_bit(start_bit + l - 1), BitRefLevel().set(), FE_DIM, level_ents);
2237  // select entities to refine
2238  fe_to_refine = intersect(level_ents, fe_to_refine);
2239  // select entities not to refine
2240  level_ents = subtract(level_ents, fe_to_refine);
2241 
2242  // for entities to refine get children, i.e. redlined entities
2243  Range fe_to_refine_children;
2244  bit_mng->updateRangeByChildren(fe_to_refine, fe_to_refine_children);
2245  // add entities to to level "l"
2246  fe_to_refine_children = fe_to_refine_children.subset_by_dimension(FE_DIM);
2247  level_ents.merge(fe_to_refine_children);
2248 
2249  auto fix_neighbour_level = [&](auto ll) {
2251  // filter entities on level ll
2252  auto level_ll = level_ents;
2253  CHKERR bit_mng->filterEntitiesByRefLevel(set_bit(ll), BitRefLevel().set(),
2254  level_ll);
2255  // find skin of ll level
2256  Range skin_edges;
2257  CHKERR skin.find_skin(0, level_ll, false, skin_edges);
2258  // get parents of skin of level ll
2259  Range skin_parents;
2260  for (auto lll = 0; lll <= ll; ++lll) {
2261  CHKERR bit_mng->updateRangeByParent(skin_edges, skin_parents);
2262  }
2263  // filter parents on level ll - 1
2264  BitRefLevel bad_bit;
2265  for (auto lll = 0; lll <= ll - 2; ++lll) {
2266  bad_bit[lll] = true;
2267  }
2268  // get adjacents to parents
2269  Range skin_adj_ents;
2270  CHKERR mField.get_moab().get_adjacencies(
2271  skin_parents, FE_DIM, false, skin_adj_ents, moab::Interface::UNION);
2272  CHKERR bit_mng->filterEntitiesByRefLevel(bad_bit, BitRefLevel().set(),
2273  skin_adj_ents);
2274  skin_adj_ents = intersect(skin_adj_ents, level_ents);
2275  if (!skin_adj_ents.empty()) {
2276  level_ents = subtract(level_ents, skin_adj_ents);
2277  Range skin_adj_ents_children;
2278  bit_mng->updateRangeByChildren(skin_adj_ents, skin_adj_ents_children);
2279  level_ents.merge(skin_adj_ents_children);
2280  }
2282  };
2283 
2284  CHKERR fix_neighbour_level(l);
2285 
2286  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2287  level_ents);
2288 
2289  // get lower dimension entities for level "l"
2290  for (auto d = 0; d != FE_DIM; ++d) {
2291  if (d == 0) {
2292  CHKERR mField.get_moab().get_connectivity(
2293  level_ents.subset_by_dimension(FE_DIM), level_ents, true);
2294  } else {
2295  CHKERR mField.get_moab().get_adjacencies(
2296  level_ents.subset_by_dimension(FE_DIM), d, false, level_ents,
2297  moab::Interface::UNION);
2298  }
2299  }
2300  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2301  level_ents);
2302 
2303  // set bit ref level to level entities
2304  CHKERR bit_mng->setNthBitRefLevel(start_bit + l, false);
2305  CHKERR bit_mng->setNthBitRefLevel(level_ents, start_bit + l, true);
2306 
2307 #ifndef NDEBUG
2308  auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
2309  CHKERR bit_mng->writeBitLevelByDim(
2310  set_bit(start_bit + l), BitRefLevel().set(), FE_DIM,
2311  (boost::lexical_cast<std::string>(l) + "_" + proc_str + "_ref_mesh.vtk")
2312  .c_str(),
2313  "VTK", "");
2314 #endif
2315 
2317  };
2318 
2319  // set skeleton
2320  auto set_skelton_bit = [&](auto l) {
2322 
2323  // get entities of dim-1 on level "l"
2324  Range level_edges;
2325  CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2326  set_bit(start_bit + l), BitRefLevel().set(), FE_DIM - 1, level_edges);
2327 
2328  // get parent of entities of level "l"
2329  Range level_edges_parents;
2330  CHKERR bit_mng->updateRangeByParent(level_edges, level_edges_parents);
2331  level_edges_parents = level_edges_parents.subset_by_dimension(FE_DIM - 1);
2332  CHKERR bit_mng->filterEntitiesByRefLevel(
2333  set_bit(start_bit + l), BitRefLevel().set(), level_edges_parents);
2334 
2335  // skeleton entities which do not have parents
2336  auto parent_skeleton = intersect(level_edges, level_edges_parents);
2337  auto skeleton = subtract(level_edges, level_edges_parents);
2338 
2339  // add adjacent domain entities
2340  CHKERR mField.get_moab().get_adjacencies(unite(parent_skeleton, skeleton),
2341  FE_DIM, false, skeleton,
2342  moab::Interface::UNION);
2343 
2344  // set levels
2345  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(skeleton);
2346  CHKERR bit_mng->setNthBitRefLevel(skeleton_bit, false);
2347  CHKERR bit_mng->setNthBitRefLevel(skeleton, skeleton_bit, true);
2348 
2349 #ifndef NDEBUG
2350  CHKERR bit_mng->writeBitLevel(
2351  set_bit(skeleton_bit), BitRefLevel().set(),
2352  (boost::lexical_cast<std::string>(l) + "_" + proc_str + "_skeleton.vtk")
2353  .c_str(),
2354  "VTK", "");
2355 #endif
2357  };
2358 
2359  // Reset bit sand set old current and aggregate bits as projection bits
2360  Range level0_current;
2361  CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(current_bit),
2362  BitRefLevel().set(), level0_current);
2363 
2364  Range level0_aggregate;
2365  CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(aggregate_bit),
2366  BitRefLevel().set(), level0_aggregate);
2367 
2368  BitRefLevel start_mask;
2369  for (auto s = 0; s != start_bit; ++s)
2370  start_mask[s] = true;
2371  CHKERR bit_mng->lambdaBitRefLevel(
2372  [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
2373  CHKERR bit_mng->setNthBitRefLevel(level0_current, projection_bit, true);
2374  CHKERR bit_mng->setNthBitRefLevel(level0_aggregate, aggregate_projection_bit,
2375  true);
2376 
2377  // Set zero bit ref level
2378  Range level0;
2379  CHKERR bit_mng->getEntitiesByRefLevel(set_bit(0), BitRefLevel().set(),
2380  level0);
2381  CHKERR bit_mng->setNthBitRefLevel(level0, start_bit, true);
2382  CHKERR bit_mng->setNthBitRefLevel(level0, current_bit, true);
2383  CHKERR bit_mng->setNthBitRefLevel(level0, aggregate_bit, true);
2384  CHKERR bit_mng->setNthBitRefLevel(level0, skeleton_bit, true);
2385 
2386  CHKERR wp.setBits(*this, 0);
2387  CHKERR wp.runCalcs(*this, 0);
2388  for (auto l = 0; l != nb_levels; ++l) {
2389  MOFEM_LOG("WORLD", Sev::inform) << "Process level: " << l;
2390  CHKERR refine_mesh(l + 1, get_refined_elements_meshset(
2391  set_bit(start_bit + l), BitRefLevel().set()));
2392  CHKERR set_skelton_bit(l + 1);
2393  CHKERR wp.setAggregateBit(*this, l + 1);
2394  CHKERR wp.setBits(*this, l + 1);
2395  CHKERR wp.runCalcs(*this, l + 1);
2396  }
2397 
2399 }

◆ runProblem()

MoFEMErrorCode LevelSet::runProblem ( )
Examples
level_set.cpp.

Definition at line 386 of file level_set.cpp.

386  {
388  CHKERR readMesh();
390 
391  if constexpr (debug) {
392  CHKERR testSideFE();
393  CHKERR testOp();
394  }
395 
397 
398  maxPtr = boost::make_shared<double>(0);
399  CHKERR refineMesh(WrapperClassInitalSolution(maxPtr));
400 
401  auto simple = mField.getInterface<Simple>();
402  simple->getBitRefLevel() =
404  simple->getBitRefLevelMask() = BitRefLevel().set();
405  simple->reSetUp(true);
406 
408 
410 }

◆ setUpProblem()

MoFEMErrorCode LevelSet::setUpProblem ( )
private

create fields, and set approximation order

Returns
MoFEMErrorCode
Examples
level_set.cpp.

Definition at line 560 of file level_set.cpp.

560  {
562  auto simple = mField.getInterface<Simple>();
563  // Scalar fields and vector field is tested. Add more fields, i.e. vector
564  // field if needed.
565  CHKERR simple->addDomainField("L", L2, AINSWORTH_LEGENDRE_BASE, 1);
566  CHKERR simple->addDomainField("V", potential_velocity_space,
568 
569  // set fields order, i.e. for most first cases order is sufficient.
570  CHKERR simple->setFieldOrder("L", 4);
571  CHKERR simple->setFieldOrder("V", 4);
572 
573  // setup problem
574  CHKERR simple->setUp();
575 
577 }

◆ solveAdvection()

MoFEMErrorCode LevelSet::solveAdvection ( )
private

solve advection problem

Returns
* MoFEMErrorCode
Examples
level_set.cpp.

Definition at line 1934 of file level_set.cpp.

1934  {
1936 
1937  // get operators tester
1938  auto simple = mField.getInterface<Simple>();
1939  auto pip = mField.getInterface<PipelineManager>(); // get interface to
1940  auto prb_mng = mField.getInterface<ProblemsManager>();
1941 
1942  CHKERR pushOpDomain();
1944 
1945  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
1946  CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "ADVECTION");
1947  CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
1948  CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
1949  CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
1950  CHKERR DMMoFEMAddElement(sub_dm, simple->getSkeletonFEName());
1951  CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
1952  CHKERR DMSetUp(sub_dm);
1953 
1954  BitRefLevel remove_mask = BitRefLevel().set(current_bit);
1955  remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
1956  CHKERR prb_mng->removeDofsOnEntities("ADVECTION", "L", BitRefLevel().set(),
1957  remove_mask);
1958 
1959  auto add_post_proc_fe = [&]() {
1960  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
1961 
1962  Tag th_error;
1963  double def_val = 0;
1964  CHKERR mField.get_moab().tag_get_handle(
1965  "Error", 1, MB_TYPE_DOUBLE, th_error, MB_TAG_CREAT | MB_TAG_SPARSE,
1966  &def_val);
1967  post_proc_fe->setTagsToTransfer({th_error});
1968 
1969  post_proc_fe->exeTestHook = [&](FEMethod *fe_ptr) {
1970  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1971  current_bit);
1972  };
1973 
1975 
1976  auto vel_ptr = boost::make_shared<MatrixDouble>();
1977 
1979  post_proc_fe->getOpPtrVector(), {L2});
1980  post_proc_fe->getOpPtrVector().push_back(getZeroLevelVelOp(vel_ptr));
1981 
1982  if constexpr (DIM1 == 1 && DIM2 == 1) {
1983  auto l_vec = boost::make_shared<VectorDouble>();
1984  post_proc_fe->getOpPtrVector().push_back(
1985  new OpCalculateScalarFieldValues("L", l_vec));
1986  post_proc_fe->getOpPtrVector().push_back(
1987 
1988  new OpPPMap(
1989 
1990  post_proc_fe->getPostProcMesh(),
1991 
1992  post_proc_fe->getMapGaussPts(),
1993 
1994  {{"L", l_vec}},
1995 
1996  {{"V", vel_ptr}},
1997 
1998  {}, {})
1999 
2000  );
2001  }
2002  return post_proc_fe;
2003  };
2004 
2005  auto post_proc_fe = add_post_proc_fe();
2006 
2007  auto set_time_monitor = [&](auto dm, auto ts) {
2008  auto monitor_ptr = boost::make_shared<FEMethod>();
2009 
2010  monitor_ptr->preProcessHook = []() { return 0; };
2011  monitor_ptr->operatorHook = []() { return 0; };
2012  monitor_ptr->postProcessHook = [&]() {
2014 
2015  if (!post_proc_fe)
2016  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
2017  "Null pointer for post proc element");
2018 
2019  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
2020  post_proc_fe);
2021  CHKERR post_proc_fe->writeFile(
2022  "level_set_" +
2023  boost::lexical_cast<std::string>(monitor_ptr->ts_step) + ".h5m");
2025  };
2026 
2027  boost::shared_ptr<FEMethod> null;
2028  DMMoFEMTSSetMonitor(sub_dm, ts, simple->getDomainFEName(), monitor_ptr,
2029  null, null);
2030 
2031  return monitor_ptr;
2032  };
2033 
2034  auto ts = pip->createTSIM(sub_dm);
2035 
2036  auto set_solution = [&](auto ts) {
2038  auto D = createDMVector(sub_dm);
2039  CHKERR DMoFEMMeshToLocalVector(sub_dm, D, INSERT_VALUES, SCATTER_FORWARD);
2040  CHKERR TSSetSolution(ts, D);
2042  };
2043  CHKERR set_solution(ts);
2044 
2045  auto monitor_pt = set_time_monitor(sub_dm, ts);
2046  CHKERR TSSetFromOptions(ts);
2047 
2048  auto B = createDMMatrix(sub_dm);
2049  CHKERR TSSetIJacobian(ts, B, B, TsSetIJacobian, nullptr);
2050  level_set_raw_ptr = this;
2051 
2052  CHKERR TSSetUp(ts);
2053 
2054  auto ts_pre_step = [](TS ts) {
2055  auto &m_field = level_set_raw_ptr->mField;
2056  auto simple = m_field.getInterface<Simple>();
2057  auto bit_mng = m_field.getInterface<BitRefManager>();
2059 
2060  auto [error, th_error] = level_set_raw_ptr->evaluateError();
2061  MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
2062 
2063  auto get_norm = [&](auto x) {
2064  double nrm;
2065  CHKERR VecNorm(x, NORM_2, &nrm);
2066  return nrm;
2067  };
2068 
2069  auto set_solution = [&](auto ts) {
2071  DM dm;
2072  CHKERR TSGetDM(ts, &dm);
2073  auto prb_ptr = getProblemPtr(dm);
2074 
2075  auto x = createDMVector(dm);
2076  CHKERR DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_FORWARD);
2077  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
2078  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
2079 
2080  MOFEM_LOG("LevelSet", Sev::inform)
2081  << "Problem " << prb_ptr->getName() << " solution vector norm "
2082  << get_norm(x);
2083  CHKERR TSSetSolution(ts, x);
2084 
2086  };
2087 
2088  auto refine_and_project = [&](auto ts) {
2090 
2092  WrapperClassErrorProjection(level_set_raw_ptr->maxPtr));
2093  simple->getBitRefLevel() = BitRefLevel().set(skeleton_bit) |
2094  BitRefLevel().set(aggregate_bit) |
2096 
2097  simple->reSetUp(true);
2098  DM dm;
2099  CHKERR TSGetDM(ts, &dm);
2101 
2102  BitRefLevel remove_mask = BitRefLevel().set(current_bit);
2103  remove_mask
2104  .flip(); // DOFs which are not on bit_domain_ele should be removed
2106  ->removeDofsOnEntities("ADVECTION", "L", BitRefLevel().set(),
2107  remove_mask);
2108 
2110  };
2111 
2112  auto ts_reset_theta = [&](auto ts) {
2114  DM dm;
2115  CHKERR TSGetDM(ts, &dm);
2116 
2117  // FIXME: Look into vec-5 how to transfer internal theta method variables
2118 
2119  CHKERR TSReset(ts);
2120  CHKERR TSSetUp(ts);
2121 
2123  CHKERR set_solution(ts);
2124 
2125  auto B = createDMMatrix(dm);
2126  CHKERR TSSetIJacobian(ts, B, B, TsSetIJacobian, nullptr);
2127 
2129  };
2130 
2131  CHKERR refine_and_project(ts);
2132  CHKERR ts_reset_theta(ts);
2133 
2135  };
2136 
2137  auto ts_post_step = [](TS ts) { return 0; };
2138 
2139  CHKERR TSSetPreStep(ts, ts_pre_step);
2140  CHKERR TSSetPostStep(ts, ts_post_step);
2141 
2142  CHKERR TSSolve(ts, NULL);
2143 
2145 }

◆ testOp()

MoFEMErrorCode LevelSet::testOp ( )
private

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.

1602  {
1604 
1605  // get operators tester
1606  auto simple = mField.getInterface<Simple>();
1607  auto opt = mField.getInterface<OperatorsTester>(); // get interface to
1608  // OperatorsTester
1609  auto pip = mField.getInterface<PipelineManager>(); // get interface to
1610  // pipeline manager
1611 
1612  CHKERR pushOpDomain();
1614 
1615  auto post_proc = [&](auto dm, auto f_res, auto out_name) {
1617  auto post_proc_fe =
1618  boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
1619 
1620  if constexpr (DIM1 == 1 && DIM2 == 1) {
1622 
1623  auto l_vec = boost::make_shared<VectorDouble>();
1624  post_proc_fe->getOpPtrVector().push_back(
1625  new OpCalculateScalarFieldValues("L", l_vec, f_res));
1626 
1627  post_proc_fe->getOpPtrVector().push_back(
1628 
1629  new OpPPMap(
1630 
1631  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1632 
1633  {{"L", l_vec}},
1634 
1635  {},
1636 
1637  {}, {})
1638 
1639  );
1640  }
1641 
1642  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
1643  post_proc_fe);
1644  post_proc_fe->writeFile(out_name);
1646  };
1647 
1648  constexpr double eps = 1e-4;
1649 
1650  auto x =
1651  opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}, {"V", {-1, 1}}});
1652  auto dot_x = opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}});
1653  auto diff_x = opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}});
1654 
1655  auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline,
1656  auto rhs_pipeline) {
1658 
1659  auto diff_res = opt->checkCentralFiniteDifference(
1660  simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
1661  SmartPetscObj<Vec>(), diff_x, 0, 1, eps);
1662 
1663  if constexpr (debug) {
1664  // Example how to plot direction in direction diff_x. If instead
1665  // directionalCentralFiniteDifference(...) diff_res is used, then error
1666  // on directive is plotted.
1667  CHKERR post_proc(simple->getDM(), diff_res, "tangent_op_error.h5m");
1668  }
1669 
1670  // Calculate norm of difference between directive calculated from finite
1671  // difference, and tangent matrix.
1672  double fnorm;
1673  CHKERR VecNorm(diff_res, NORM_2, &fnorm);
1674  MOFEM_LOG_C("LevelSet", Sev::inform,
1675  "Test consistency of tangent matrix %3.4e", fnorm);
1676 
1677  constexpr double err = 1e-9;
1678  if (fnorm > err)
1679  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1680  "Norm of directional derivative too large err = %3.4e", fnorm);
1681 
1683  };
1684 
1685  CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
1686  pip->getDomainRhsFE());
1687  CHKERR test_domain_ops(simple->getSkeletonFEName(), pip->getSkeletonLhsFE(),
1688  pip->getSkeletonRhsFE());
1689 
1691 };

◆ testSideFE()

MoFEMErrorCode LevelSet::testSideFE ( )
private

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.

1387  {
1389 
1390  /**
1391  * @brief calculate volume
1392  *
1393  */
1394  struct DivergenceVol : public DomainEleOp {
1395  DivergenceVol(boost::shared_ptr<MatrixDouble> l_ptr,
1396  boost::shared_ptr<MatrixDouble> vel_ptr,
1397  SmartPetscObj<Vec> div_vec)
1398  : DomainEleOp("L", DomainEleOp::OPROW), lPtr(l_ptr), velPtr(vel_ptr),
1399  divVec(div_vec) {}
1400  MoFEMErrorCode doWork(int side, EntityType type,
1403  const auto nb_dofs = data.getIndices().size();
1404  if (nb_dofs) {
1405  const auto nb_gauss_pts = getGaussPts().size2();
1406  const auto t_w = getFTensor0IntegrationWeight();
1407  auto t_diff = data.getFTensor1DiffN<SPACE_DIM>();
1408  auto t_l = getFTensor2FromMat<DIM1, DIM2>(*lPtr);
1409  auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
1410  double div = 0;
1411  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1412  for (int rr = 0; rr != nb_dofs; ++rr) {
1413  div += getMeasure() * t_w * t_l(I, I) * (t_diff(i) * t_vel(i));
1414  ++t_diff;
1415  }
1416  ++t_w;
1417  ++t_l;
1418  ++t_vel;
1419  }
1420  CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
1421  }
1423  }
1424 
1425  private:
1426  boost::shared_ptr<MatrixDouble> lPtr;
1427  boost::shared_ptr<MatrixDouble> velPtr;
1428  SmartPetscObj<Vec> divVec;
1429  };
1430 
1431  /**
1432  * @brief calculate skeleton integral
1433  *
1434  */
1435  struct DivergenceSkeleton : public BoundaryEleOp {
1436  DivergenceSkeleton(boost::shared_ptr<SideData> side_data_ptr,
1437  boost::shared_ptr<FaceSideEle> side_fe_ptr,
1438  SmartPetscObj<Vec> div_vec)
1439  : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
1440  sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr), divVec(div_vec) {}
1441  MoFEMErrorCode doWork(int side, EntityType type,
1444 
1445  auto get_ntensor = [](auto &base_mat) {
1447  &*base_mat.data().begin());
1448  };
1449 
1450  auto not_side = [](auto s) {
1451  return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
1452  };
1453 
1454  // Collect data from side domain elements
1455  CHKERR loopSideFaces("dFE", *sideFEPtr);
1456  const auto in_the_loop =
1457  sideFEPtr->nInTheLoop; // return number of elements on the side
1458 
1459  auto t_normal = getFTensor1Normal();
1460  const auto nb_gauss_pts = getGaussPts().size2();
1461  for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
1462  const auto nb_dofs = sideDataPtr->indicesRowSideMap[s0].size();
1463  if (nb_dofs) {
1464  auto t_base = get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
1465  auto nb_row_base_functions = sideDataPtr->rowBaseSideMap[s0].size2();
1466  auto side_sense = sideDataPtr->senseMap[s0];
1467  auto opposite_s0 = not_side(s0);
1468 
1469  auto arr_t_l = make_array(
1470  getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[LEFT_SIDE]),
1471  getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[RIGHT_SIDE]));
1472  auto arr_t_vel = make_array(
1473  getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[LEFT_SIDE]),
1474  getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[RIGHT_SIDE]));
1475 
1476  auto next = [&]() {
1477  for (auto &t_l : arr_t_l)
1478  ++t_l;
1479  for (auto &t_vel : arr_t_vel)
1480  ++t_vel;
1481  };
1482 
1483  double div = 0;
1484 
1485  auto t_w = getFTensor0IntegrationWeight();
1486  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1488  t_vel(i) =
1489  (arr_t_vel[LEFT_SIDE](i) + arr_t_vel[RIGHT_SIDE](i)) / 2.;
1490  const auto dot = (t_normal(i) * t_vel(i)) * side_sense;
1491  const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
1492  const auto l_upwind =
1493  arr_t_l[l_upwind_side]; //< assume that field is continues,
1494  // initialisation field has to be smooth
1495  // and exactly approximated by approx
1496  // base
1497  auto res = t_w * l_upwind(I, I) * dot;
1498  ++t_w;
1499  next();
1500  int rr = 0;
1501  for (; rr != nb_dofs; ++rr) {
1502  div += t_base * res;
1503  ++t_base;
1504  }
1505  for (; rr < nb_row_base_functions; ++rr) {
1506  ++t_base;
1507  }
1508  }
1509  CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
1510  }
1511  if (!in_the_loop)
1512  break;
1513  }
1514 
1516  }
1517 
1518  private:
1519  boost::shared_ptr<SideData> sideDataPtr;
1520  boost::shared_ptr<FaceSideEle> sideFEPtr;
1521  boost::shared_ptr<MatrixDouble> velPtr;
1522  SmartPetscObj<Vec> divVec;
1523  };
1524 
1525  auto vol_fe = boost::make_shared<DomainEle>(mField);
1526  auto skel_fe = boost::make_shared<BoundaryEle>(mField);
1527 
1528  vol_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
1529  skel_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
1530 
1531  auto div_vol_vec = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
1532  auto div_skel_vec = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
1533 
1534  auto l_ptr = boost::make_shared<MatrixDouble>();
1535  auto vel_ptr = boost::make_shared<MatrixDouble>();
1536  auto side_data_ptr = boost::make_shared<SideData>();
1537  auto side_fe_ptr = getSideFE(side_data_ptr);
1538 
1539  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(vol_fe->getOpPtrVector(),
1540  {L2});
1541  vol_fe->getOpPtrVector().push_back(
1543  vol_fe->getOpPtrVector().push_back(getZeroLevelVelOp(vel_ptr));
1544  vol_fe->getOpPtrVector().push_back(
1545  new DivergenceVol(l_ptr, vel_ptr, div_vol_vec));
1546 
1547  skel_fe->getOpPtrVector().push_back(
1548  new DivergenceSkeleton(side_data_ptr, side_fe_ptr, div_skel_vec));
1549 
1550  auto simple = mField.getInterface<Simple>();
1551  auto dm = simple->getDM();
1552 
1553  /**
1554  * Set up problem such that gradient of level set field is orthogonal to
1555  * velocity field. Then volume and skeleton integral should yield the same
1556  * value.
1557  */
1558 
1560  [](double x, double y, double) { return x - y; });
1562  [](double x, double y, double) { return x - y; });
1563 
1564  vol_fe->exeTestHook = [&](FEMethod *fe_ptr) {
1565  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1566  current_bit);
1567  };
1568  skel_fe->exeTestHook = [&](FEMethod *fe_ptr) {
1569  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1570  skeleton_bit);
1571  };
1572 
1573  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), vol_fe);
1574  CHKERR DMoFEMLoopFiniteElements(dm, simple->getSkeletonFEName(), skel_fe);
1575  CHKERR DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(), skel_fe);
1576 
1577  auto assemble_and_sum = [](auto vec) {
1578  CHK_THROW_MESSAGE(VecAssemblyBegin(vec), "assemble");
1579  CHK_THROW_MESSAGE(VecAssemblyEnd(vec), "assemble");
1580  double sum;
1581  CHK_THROW_MESSAGE(VecSum(vec, &sum), "assemble");
1582  return sum;
1583  };
1584 
1585  auto div_vol = assemble_and_sum(div_vol_vec);
1586  auto div_skel = assemble_and_sum(div_skel_vec);
1587 
1588  auto eps = std::abs((div_vol - div_skel) / (div_vol + div_skel));
1589 
1590  MOFEM_LOG("WORLD", Sev::inform) << "Testing divergence volume: " << div_vol;
1591  MOFEM_LOG("WORLD", Sev::inform) << "Testing divergence skeleton: " << div_skel
1592  << " relative difference: " << eps;
1593 
1594  constexpr double eps_err = 1e-6;
1595  if (eps > eps_err)
1596  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1597  "No consistency between skeleton integral and volume integral");
1598 
1600 };

Member Data Documentation

◆ maxPtr

boost::shared_ptr<double> LevelSet::maxPtr
private
Examples
level_set.cpp.

Definition at line 370 of file level_set.cpp.

◆ mField

MoFEM::Interface& LevelSet::mField
private

integrate skeleton operators on khs

Examples
level_set.cpp.

Definition at line 346 of file level_set.cpp.


The documentation for this struct was generated from the following file:
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
aggregate_projection_bit
constexpr int aggregate_projection_bit
all bits for projection problem
Definition: level_set.cpp:70
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
LevelSet::readMesh
MoFEMErrorCode readMesh()
read mesh
Definition: level_set.cpp:502
MoFEM::CoreInterface::get_finite_element_meshset
virtual EntityHandle get_finite_element_meshset(const std::string name) const =0
MoFEM::DMMoFEMKSPSetComputeRHS
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
Definition: DMMoFEM.cpp:637
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition: DataOperators.hpp:24
LevelSet::maxPtr
boost::shared_ptr< double > maxPtr
Definition: level_set.cpp:370
LevelSet::getSideFE
boost::shared_ptr< FaceSideEle > getSideFE(boost::shared_ptr< SideData > side_data_ptr)
create side element to assemble data from sides
Definition: level_set.cpp:997
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::OpCalculateTensor2FieldValues
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition: UserDataOperators.hpp:887
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
NOISY
@ NOISY
Definition: definitions.h:224
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
LevelSet::OpMassLL
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< G >::OpMass< 1, DIM1 *DIM2 > OpMassLL
Definition: level_set.cpp:354
MoFEM::OpRunParent
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Definition: MeshProjectionDataOperators.hpp:17
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
sdf_hertz_2d_axisymm_plane.zc
int zc
Definition: sdf_hertz_2d_axisymm_plane.py:8
MoFEM::getDMKspCtx
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
Definition: DMMoFEM.hpp:1113
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
LevelSet::OpScalarFieldL
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpBaseTimesVector< 1, DIM1 *DIM2, 1 > OpScalarFieldL
Definition: level_set.cpp:362
LevelSet::solveAdvection
MoFEMErrorCode solveAdvection()
solve advection problem
Definition: level_set.cpp:1934
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
nb_levels
constexpr int nb_levels
Definition: level_set.cpp:58
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
aggregate_bit
constexpr int aggregate_bit
all bits for advection problem
Definition: level_set.cpp:66
FaceSideEleOp
FaceSideEle::UserDataOperator FaceSideEleOp
Definition: level_set.cpp:45
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
LevelSet::initialiseFieldLevelSet
MoFEMErrorCode initialiseFieldLevelSet(boost::function< double(double, double, double)> level_fun=get_level_set)
initialise field set
Definition: level_set.cpp:1693
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
LevelSet::OpSourceL
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpSource< 1, DIM1 *DIM2 > OpSourceL
Definition: level_set.cpp:356
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
potential_velocity_field_dim
constexpr size_t potential_velocity_field_dim
Definition: level_set.cpp:50
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:261
MoFEM::DMMoFEMKSPSetComputeOperators
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
Definition: DMMoFEM.cpp:678
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
sdf.r
int r
Definition: sdf.py:8
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1056
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
I
FTensor::Index< 'I', DIM1 > I
Definition: level_set.cpp:29
FE_DIM
constexpr int FE_DIM
[Define dimension]
Definition: level_set.cpp:19
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
LevelSet::testOp
MoFEMErrorCode testOp()
test consistency between tangent matrix and the right hand side vectors
Definition: level_set.cpp:1602
LevelSet::pushOpSkeleton
MoFEMErrorCode pushOpSkeleton()
push operator to integrate on skeleton
Definition: level_set.cpp:1173
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:438
SPACE_DIM
constexpr int SPACE_DIM
Definition: level_set.cpp:20
LevelSet::OpSourceV
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpSource< potential_velocity_field_dim, potential_velocity_field_dim > OpSourceV
Definition: level_set.cpp:360
current_bit
constexpr int current_bit
dofs bit used to do calculations
Definition: level_set.cpp:63
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
BoundaryEleOp
MoFEM::DMSubDMSetUp_MoFEM
PetscErrorCode DMSubDMSetUp_MoFEM(DM subdm)
Definition: DMMoFEM.cpp:1350
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1886
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
DIM1
constexpr int DIM1
Definition: level_set.cpp:21
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Definition: level_set.cpp:43
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:215
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
LevelSet::setUpProblem
MoFEMErrorCode setUpProblem()
create fields, and set approximation order
Definition: level_set.cpp:560
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
start_bit
constexpr int start_bit
Definition: level_set.cpp:60
Poisson2DiscontGalerkinOperators::get_ntensor
auto get_ntensor(T &base_mat)
Definition: PoissonDiscontinousGalerkin.hpp:90
level_set_raw_ptr
LevelSet * level_set_raw_ptr
Definition: level_set.cpp:1931
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
MoFEM::TsSetIJacobian
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition: TsCtx.cpp:165
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: PipelineManager.hpp:749
LevelSet::OpMassVV
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< G >::OpMass< potential_velocity_field_dim, potential_velocity_field_dim > OpMassVV
Definition: level_set.cpp:358
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
MoFEM::OpCalculateTensor2FieldValuesDot
Get time direvarive values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition: UserDataOperators.hpp:902
MoFEM::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition: OperatorsTester.hpp:21
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
MoFEM::PipelineManager::getOpSkeletonLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpSkeletonLhsPipeline()
Get the Op Skeleton Lhs Pipeline object.
Definition: PipelineManager.hpp:845
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: level_set.cpp:579
MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:249
LevelSet::evaluateError
std::tuple< double, Tag > evaluateError()
evaluate error
Definition: level_set.cpp:1203
sdf_hertz_2d_axisymm_plane.yc
int yc
Definition: sdf_hertz_2d_axisymm_plane.py:7
LevelSet::mField
MoFEM::Interface & mField
integrate skeleton operators on khs
Definition: level_set.cpp:346
Range
DomainEleOp
save_range
auto save_range
Definition: thermo_elastic.cpp:144
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
FTensor::Tensor0
Definition: Tensor0.hpp:16
skeleton_bit
constexpr int skeleton_bit
skeleton elements bit
Definition: level_set.cpp:65
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1318
LevelSet::initialiseFieldVelocity
MoFEMErrorCode initialiseFieldVelocity(boost::function< double(double, double, double)> vel_fun=get_velocity_potential< FE_DIM >)
initialise potential velocity field
Definition: level_set.cpp:1814
LevelSet::refineMesh
MoFEMErrorCode refineMesh(WrapperClass &&wp)
Definition: level_set.cpp:2147
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
LevelSet::dgProjection
MoFEMErrorCode dgProjection(const int prj_bit=projection_bit)
dg level set projection
Definition: level_set.cpp:2401
debug
constexpr bool debug
Definition: level_set.cpp:53
MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::getProblemPtr
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition: DMMoFEM.hpp:1044
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Simple::getBitRefLevel
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
Definition: Simple.hpp:354
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
DIM2
constexpr int DIM2
Definition: level_set.cpp:22
LevelSet::LEFT_SIDE
@ LEFT_SIDE
Definition: level_set.cpp:367
LevelSet::getZeroLevelVelOp
ForcesAndSourcesCore::UserDataOperator * getZeroLevelVelOp(boost::shared_ptr< MatrixDouble > vel_ptr)
Get operator calculating velocity on coarse mesh.
Definition: level_set.cpp:922
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
LevelSet::pushOpDomain
MoFEMErrorCode pushOpDomain()
push operators to integrate operators on domain
Definition: level_set.cpp:954
sdf_hertz_2d_axisymm_plane.xc
int xc
Definition: sdf_hertz_2d_axisymm_plane.py:6
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
projection_bit
constexpr int projection_bit
Definition: level_set.cpp:68
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::make_array
constexpr auto make_array(Arg &&...arg)
Create Array.
Definition: Templates.hpp:1990
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
MoFEM::SmartPetscObj< Vec >
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
potential_velocity_space
constexpr FieldSpace potential_velocity_space
Definition: level_set.cpp:49
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
LevelSet::RIGHT_SIDE
@ RIGHT_SIDE
Definition: level_set.cpp:367
MoFEM::OpCalculateHcurlVectorCurl
Calculate curl of vector field.
Definition: UserDataOperators.hpp:2488
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
F
@ F
Definition: free_surface.cpp:394
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
LevelSet::testSideFE
MoFEMErrorCode testSideFE()
test integration side elements
Definition: level_set.cpp:1387