v0.16.3
Loading...
Searching...
No Matches
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
 
MoFEMErrorCode setUpProblem ()
 create fields, and set approximation order
 
MoFEMErrorCode pushOpDomain ()
 push operators to integrate operators on domain
 
std::tuple< double, Tag > evaluateError ()
 evaluate error
 
ForcesAndSourcesCore::UserDataOperator * getZeroLevelVelOp (boost::shared_ptr< MatrixDouble > vel_ptr)
 Get operator calculating velocity on coarse mesh.
 
boost::shared_ptr< FaceSideEle > getSideFE (boost::shared_ptr< SideData > side_data_ptr)
 create side element to assemble data from sides
 
MoFEMErrorCode pushOpSkeleton ()
 push operator to integrate on skeleton
 
MoFEMErrorCode testSideFE ()
 test integration side elements
 
MoFEMErrorCode testOp ()
 test consistency between tangent matrix and the right hand side vectors
 
MoFEMErrorCode initialiseFieldLevelSet (boost::function< double(double, double, double)> level_fun=get_level_set)
 initialise field set
 
MoFEMErrorCode initialiseFieldVelocity (boost::function< double(double, double, double)> vel_fun=get_velocity_potential< FE_DIM >)
 initialise potential velocity field
 
MoFEMErrorCode dgProjection (const int prj_bit=projection_bit)
 dg level set projection
 
MoFEMErrorCode solveAdvection ()
 solve advection problem
 
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
 
static double get_level_set (const double x, const double y, const double z)
 inital level set, i.e. advected field
 
template<>
double get_velocity_potential (double x, double y, double z)
 

Private Attributes

MoFEM::Interface & mField
 integrate skeleton operators on khs
 
boost::shared_ptr< double > maxPtr
 

Detailed Description

Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 73 of file level_set.cpp.

Member Typedef Documentation

◆ AssemblyBoundaryEleOp

using LevelSet::AssemblyBoundaryEleOp = FormsIntegrators<BoundaryEleOp>::Assembly<A>::OpBase
private

Definition at line 364 of file level_set.cpp.

◆ AssemblyDomainEleOp

using LevelSet::AssemblyDomainEleOp = FormsIntegrators<DomainEleOp>::Assembly<A>::OpBase
private

Definition at line 351 of file level_set.cpp.

◆ MatSideArray

using LevelSet::MatSideArray = std::array<MatrixDouble, 2>
private
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 80 of file level_set.cpp.

◆ OpMassLL

using LevelSet::OpMassLL = FormsIntegrators<DomainEleOp>::Assembly<A>::BiLinearForm< G>::OpMass<1, DIM1 * DIM2>
private
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 353 of file level_set.cpp.

◆ OpMassVV

using LevelSet::OpMassVV = FormsIntegrators<DomainEleOp>::Assembly<A>::BiLinearForm< G>::OpMass<potential_velocity_field_dim, potential_velocity_field_dim>
private
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 357 of file level_set.cpp.

◆ OpScalarFieldL

using LevelSet::OpScalarFieldL = FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm< G>::OpBaseTimesVector<1, DIM1 * DIM2, 1>
private
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 361 of file level_set.cpp.

◆ OpSourceL

using LevelSet::OpSourceL = FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm< G>::OpSource<1, DIM1 * DIM2>
private
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 355 of file level_set.cpp.

◆ OpSourceV

using LevelSet::OpSourceV = FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm< G>::OpSource<potential_velocity_field_dim, potential_velocity_field_dim>
private
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 359 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) {}
LevelSet::mField
MoFEM::Interface & mField
integrate skeleton operators on khs
Definition level_set.cpp:349

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
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 2388 of file level_set.cpp.

2388 {
2389 MoFEMFunctionBegin;
2390
2391 // get operators tester
2392 auto simple = mField.getInterface<Simple>();
2393 auto bit_mng = mField.getInterface<BitRefManager>();
2394 auto prb_mng = mField.getInterface<ProblemsManager>();
2395
2396 auto lhs_fe = boost::make_shared<DomainEle>(mField);
2397 auto rhs_fe_prj = boost::make_shared<DomainEle>(mField);
2398 auto rhs_fe_current = boost::make_shared<DomainEle>(mField);
2399
2400 lhs_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
2401 rhs_fe_prj->getRuleHook = [](int, int, int o) { return 3 * o; };
2402 rhs_fe_current->getRuleHook = [](int, int, int o) { return 3 * o; };
2403
2404 auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
2405 CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "DG_PROJECTION");
2406 CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
2407 CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
2408 CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
2409 CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
2410 CHKERR DMSetUp(sub_dm);
2411
2412 Range current_ents; // ents used to do calculations
2413 CHKERR bit_mng->getEntitiesByDimAndRefLevel(BitRefLevel().set(current_bit),
2414 BitRefLevel().set(), FE_DIM,
2415 current_ents);
2416 Range prj_ents; // ents from which data are projected
2417 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2418 BitRefLevel().set(projection_bit), BitRefLevel().set(), FE_DIM, prj_ents);
2419 for (auto l = 0; l != nb_levels; ++l) {
2420 CHKERR bit_mng->updateRangeByParent(prj_ents, prj_ents);
2421 }
2422 current_ents = subtract(
2423 current_ents, prj_ents); // only restric to entities needed projection
2424
2425 auto test_mesh_bit = [&](FEMethod *fe_ptr) {
2426 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2427 current_bit);
2428 };
2429 auto test_prj_bit = [&](FEMethod *fe_ptr) {
2430 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2431 projection_bit);
2432 };
2433 auto test_current_bit = [&](FEMethod *fe_ptr) {
2434 return current_ents.find(fe_ptr->getFEEntityHandle()) != current_ents.end();
2435 };
2436
2437 lhs_fe->exeTestHook =
2438 test_mesh_bit; // that element only is run when current bit is set
2439 rhs_fe_prj->exeTestHook =
2440 test_prj_bit; // that element is run only when projection bit is set
2441 rhs_fe_current->exeTestHook =
2442 test_current_bit; // that element is only run when current bit is set
2443
2444 BitRefLevel remove_mask = BitRefLevel().set(current_bit);
2445 remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
2446 CHKERR prb_mng->removeDofsOnEntities(
2447 "DG_PROJECTION", "L", BitRefLevel().set(), remove_mask, nullptr, 0,
2448 MAX_DOFS_ON_ENTITY, 0, 100, NOISY,
2449 true); // remove all DOFs which are not
2450 // on current bit. This case works for L2 space
2451
2452 CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(lhs_fe->getOpPtrVector(),
2453 {L2});
2454 lhs_fe->getOpPtrVector().push_back(
2455 new OpMassLL("L", "L")); // Assemble projection matrix
2456
2457 auto make_tensor_mat_to_vec_mat_op = [](auto src_mat_ptr, auto dst_mat_ptr) {
2458 auto op_ptr = new DomainEleOp(NOSPACE, DomainEleOp::OPSPACE);
2459 op_ptr->doWorkRhsHook = [src_mat_ptr, dst_mat_ptr](DataOperator *op_ptr,
2460 int, EntityType,
2461 EntData &) {
2462 MoFEMFunctionBegin;
2463 const auto nb_gauss_pts =
2464 static_cast<DomainEleOp *>(op_ptr)->getGaussPts().size2();
2465 auto &src_mat = *src_mat_ptr;
2466 auto &dst_mat = *dst_mat_ptr;
2467#ifndef NDEBUG
2468 if (src_mat.size1() != nb_gauss_pts)
2469 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2470 "Wrong number of integration pts %zu != %zu", src_mat.size1(),
2471 nb_gauss_pts);
2472 if (src_mat.size2() != DIM1 * DIM2)
2473 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2474 "Wrong tensor components layout %zu != %d", src_mat.size2(),
2475 DIM1 * DIM2);
2476#endif
2477 dst_mat.resize(nb_gauss_pts, DIM1 * DIM2, false);
2478 for (size_t gg = 0; gg != nb_gauss_pts; ++gg)
2479 for (size_t dd = 0; dd != DIM1 * DIM2; ++dd)
2480 dst_mat(gg, dd) = src_mat(gg, dd);
2481 MoFEMFunctionReturn(0);
2482 };
2483 return op_ptr;
2484 };
2485
2486 // This assumes that projection mesh is finer, current mesh is coarsened.
2487 auto set_prj_from_child = [&](auto rhs_fe_prj) {
2488 CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
2489 rhs_fe_prj->getOpPtrVector(), {L2});
2490
2491 // Evaluate field value on projection mesh
2492 auto l_tensor_mat = boost::make_shared<MatrixDouble>();
2493 auto l_vec_mat = boost::make_shared<MatrixDouble>();
2494 rhs_fe_prj->getOpPtrVector().push_back(
2495 new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_tensor_mat));
2496 rhs_fe_prj->getOpPtrVector().push_back(
2497 make_tensor_mat_to_vec_mat_op(l_tensor_mat, l_vec_mat));
2498
2499 // This element is used to assemble
2500 auto get_parent_this = [&]() {
2501 auto fe_parent_this = boost::make_shared<DomainParentEle>(mField);
2502 fe_parent_this->getOpPtrVector().push_back(
2503 new OpScalarFieldL("L", l_vec_mat));
2504 return fe_parent_this;
2505 };
2506
2507 // Create levels of parent elements, until current element is reached, and
2508 // then assemble.
2509 auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
2510 std::vector<boost::shared_ptr<DomainParentEle>> parents_elems_ptr_vec;
2511 for (int l = 0; l <= nb_levels; ++l)
2512 parents_elems_ptr_vec.emplace_back(
2513 boost::make_shared<DomainParentEle>(mField));
2514 for (auto l = 1; l <= nb_levels; ++l) {
2515 parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
2516 new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
2517 BitRefLevel().set(current_bit).flip(), this_fe_ptr,
2518 BitRefLevel().set(current_bit),
2519 BitRefLevel().set()));
2520 }
2521 return parents_elems_ptr_vec[0];
2522 };
2523
2524 auto this_fe_ptr = get_parent_this();
2525 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
2526 rhs_fe_prj->getOpPtrVector().push_back(
2527 new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
2528 BitRefLevel().set(current_bit).flip(), this_fe_ptr,
2529 BitRefLevel().set(current_bit), BitRefLevel().set()));
2530 };
2531
2532 // This assumed that current mesh is refined, and projection mesh is coarser
2533 auto set_prj_from_parent = [&](auto rhs_fe_current) {
2534
2535 // Evaluate field value on projection mesh
2536 auto l_tensor_mat = boost::make_shared<MatrixDouble>();
2537 auto l_vec_mat = boost::make_shared<MatrixDouble>();
2538
2539 // Evaluate field on coarser element
2540 auto get_parent_this = [&]() {
2541 auto fe_parent_this = boost::make_shared<DomainParentEle>(mField);
2542 fe_parent_this->getOpPtrVector().push_back(
2543 new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_tensor_mat));
2544 fe_parent_this->getOpPtrVector().push_back(
2545 make_tensor_mat_to_vec_mat_op(l_tensor_mat, l_vec_mat));
2546 return fe_parent_this;
2547 };
2548
2549 // Create stack of evaluation on parent elements
2550 auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
2551 std::vector<boost::shared_ptr<DomainParentEle>> parents_elems_ptr_vec;
2552 for (int l = 0; l <= nb_levels; ++l)
2553 parents_elems_ptr_vec.emplace_back(
2554 boost::make_shared<DomainParentEle>(mField));
2555 for (auto l = 1; l <= nb_levels; ++l) {
2556 parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
2557 new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
2558 BitRefLevel().set(projection_bit).flip(),
2559 this_fe_ptr, BitRefLevel().set(projection_bit),
2560 BitRefLevel().set()));
2561 }
2562 return parents_elems_ptr_vec[0];
2563 };
2564
2565 auto this_fe_ptr = get_parent_this();
2566 auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
2567
2568 auto reset_op_ptr = new DomainEleOp(NOSPACE, DomainEleOp::OPSPACE);
2569 reset_op_ptr->doWorkRhsHook = [&](DataOperator *op_ptr, int, EntityType,
2570 EntData &) {
2571 l_vec_mat->resize(
2572 static_cast<DomainEleOp *>(op_ptr)->getGaussPts().size2(),
2573 DIM1 * DIM2, false);
2574 l_vec_mat->clear();
2575 return 0;
2576 };
2577 rhs_fe_current->getOpPtrVector().push_back(reset_op_ptr);
2578 rhs_fe_current->getOpPtrVector().push_back(
2579 new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
2580 BitRefLevel().set(projection_bit).flip(), this_fe_ptr,
2581 BitRefLevel(), BitRefLevel()));
2582
2583 // At the end assemble of current finite element
2584 rhs_fe_current->getOpPtrVector().push_back(
2585 new OpScalarFieldL("L", l_vec_mat));
2586 };
2587
2588 set_prj_from_child(rhs_fe_prj);
2589 set_prj_from_parent(rhs_fe_current);
2590
2591 boost::shared_ptr<FEMethod> null_fe;
2592 getDMKspCtx(sub_dm)->clearLoops();
2593 CHKERR DMMoFEMKSPSetComputeOperators(sub_dm, simple->getDomainFEName(),
2594 lhs_fe, null_fe, null_fe);
2595 CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(), rhs_fe_prj,
2596 null_fe, null_fe);
2597 CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(),
2598 rhs_fe_current, null_fe, null_fe);
2599 auto ksp = MoFEM::createKSP(mField.get_comm());
2600 CHKERR KSPSetDM(ksp, sub_dm);
2601
2602 CHKERR KSPSetDM(ksp, sub_dm);
2603 CHKERR KSPSetFromOptions(ksp);
2604 CHKERR KSPSetUp(ksp);
2605
2606 auto L = createDMVector(sub_dm);
2607 auto F = vectorDuplicate(L);
2608
2609 CHKERR KSPSolve(ksp, F, L);
2610 CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
2611 CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
2612 CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
2613
2614 auto [error, th_error] = evaluateError();
2615 MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
2616
2617 auto post_proc = [&](auto dm, auto out_name, auto th_error) {
2618 MoFEMFunctionBegin;
2619 auto post_proc_fe =
2620 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
2621 post_proc_fe->setTagsToTransfer({th_error});
2622 post_proc_fe->exeTestHook = test_mesh_bit;
2623
2624 if constexpr (DIM1 == 1 && DIM2 == 1) {
2625 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
2626 auto l_vec = boost::make_shared<VectorDouble>();
2627 auto l_grad_mat = boost::make_shared<MatrixDouble>();
2628 CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
2629 post_proc_fe->getOpPtrVector(), {L2});
2630 post_proc_fe->getOpPtrVector().push_back(
2631 new OpCalculateScalarFieldValues("L", l_vec));
2632 post_proc_fe->getOpPtrVector().push_back(
2633 new OpCalculateScalarFieldGradient<SPACE_DIM>("L", l_grad_mat));
2634
2635 post_proc_fe->getOpPtrVector().push_back(
2636
2637 new OpPPMap(
2638
2639 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2640
2641 {{"L", l_vec}},
2642
2643 {{"GradL", l_grad_mat}},
2644
2645 {}, {})
2646
2647 );
2648 }
2649
2650 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
2651 post_proc_fe);
2652 post_proc_fe->writeFile(out_name);
2653 MoFEMFunctionReturn(0);
2654 };
2655
2656 if constexpr (debug)
2657 CHKERR post_proc(sub_dm, "dg_projection.h5m", th_error);
2658
2659 MoFEMFunctionReturn(0);
2660}
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
NOISY
@ NOISY
Definition definitions.h:224
MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition definitions.h:249
NOSPACE
@ NOSPACE
Definition definitions.h:83
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
F
@ F
Definition free_surface.cpp:625
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
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
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:627
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:576
MoFEM::createDMVector
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
Definition DMMoFEM.hpp:1237
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:668
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
FE_DIM
constexpr int FE_DIM
[Define dimension]
Definition level_set.cpp:19
current_bit
constexpr int current_bit
dofs bit used to do calculations
Definition level_set.cpp:63
nb_levels
constexpr int nb_levels
Definition level_set.cpp:58
DIM2
constexpr int DIM2
Definition level_set.cpp:22
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Definition level_set.cpp:43
projection_bit
constexpr int projection_bit
Definition level_set.cpp:68
DIM1
constexpr int DIM1
Definition level_set.cpp:21
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:23
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition ddTensor0.hpp:33
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::createKSP
auto createKSP(MPI_Comm comm)
Definition PetscSmartObj.hpp:269
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition DMMoFEM.cpp:434
MoFEM::debug
static const bool debug
Definition ConstrainMatrixCtx.cpp:9
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition PetscSmartObj.hpp:223
MoFEM::getDMKspCtx
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
Definition DMMoFEM.hpp:1251
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition PetscSmartObj.hpp:143
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition photon_diffusion.cpp:40
LevelSet::OpMassLL
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< G >::OpMass< 1, DIM1 *DIM2 > OpMassLL
Definition level_set.cpp:354
LevelSet::OpScalarFieldL
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpBaseTimesVector< 1, DIM1 *DIM2, 1 > OpScalarFieldL
Definition level_set.cpp:362
LevelSet::evaluateError
std::tuple< double, Tag > evaluateError()
evaluate error
Definition level_set.cpp:1200
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition HODataOperators.hpp:642
MoFEM::BitRefManager
Managing BitRefLevels.
Definition BitRefManager.hpp:21
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition DataOperators.hpp:24
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition EntitiesFieldData.hpp:150
MoFEM::FEMethod
Structure for user loop methods on finite elements.
Definition LoopMethods.hpp:657
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Definition UserDataOperators.hpp:1445
MoFEM::OpCalculateScalarFieldValues
Specialization for double precision scalar field values calculation.
Definition UserDataOperators.hpp:157
MoFEM::OpCalculateTensor2FieldValues
Get values at integration pts for tensor field rank 2, i.e. matrix field.
Definition UserDataOperators.hpp:1034
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition PostProcBrokenMeshInMoabBase.hpp:711
MoFEM::OpRunParent
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Definition MeshProjectionDataOperators.hpp:17
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition ProblemsManager.hpp:133
MoFEM::Simple
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93

◆ evaluateError()

std::tuple< double, Tag > LevelSet::evaluateError ( )
private

evaluate error

Returns
MoFEMErrorCode
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 1200 of file level_set.cpp.

1200 {
1201
1202 struct OpErrorSkel : BoundaryEleOp {
1203
1204 OpErrorSkel(boost::shared_ptr<FaceSideEle> side_fe_ptr,
1205 boost::shared_ptr<SideData> side_data_ptr,
1206 SmartPetscObj<Vec> error_sum_ptr, Tag th_error)
1207 : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
1208 sideFEPtr(side_fe_ptr), sideDataPtr(side_data_ptr),
1209 errorSumPtr(error_sum_ptr), thError(th_error) {}
1210
1211 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1212 MoFEMFunctionBegin;
1213
1214 // Collect data from side domain elements
1215 CHKERR loopSideFaces("dFE", *sideFEPtr);
1216
1217 auto nb_gauss_pts = getGaussPts().size2();
1218
1219 for ([[maybe_unused]] auto s : {LEFT_SIDE, RIGHT_SIDE}) {
1220
1221 auto arr_t_l = make_array(
1222 getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[LEFT_SIDE]),
1223 getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[RIGHT_SIDE]));
1224 auto arr_t_vel = make_array(
1225 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[LEFT_SIDE]),
1226 getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[RIGHT_SIDE]));
1227
1228 auto next = [&]() {
1229 for (auto &t_l : arr_t_l)
1230 ++t_l;
1231 for (auto &t_vel : arr_t_vel)
1232 ++t_vel;
1233 };
1234
1235 double e = 0;
1236 auto t_w = getFTensor0IntegrationWeight();
1237 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1238 FTensor::Tensor2<double, DIM1, DIM2> t_diff;
1239 t_diff(I, J) = arr_t_l[LEFT_SIDE](I, J) - arr_t_l[RIGHT_SIDE](I, J);
1240 e += t_w * getMeasure() * t_diff(I, J) * t_diff(I, J);
1241 next();
1242 ++t_w;
1243 }
1244 e = std::sqrt(e);
1245
1246 moab::Interface &moab =
1247 getNumeredEntFiniteElementPtr()->getBasicDataPtr()->moab;
1248 const void *tags_ptr[2];
1249 CHKERR moab.tag_get_by_ptr(thError, sideDataPtr->feSideHandle.data(), 2,
1250 tags_ptr);
1251 for (auto ff : {0, 1}) {
1252 *((double *)tags_ptr[ff]) += e;
1253 }
1254 CHKERR VecSetValue(errorSumPtr, 0, e, ADD_VALUES);
1255 };
1256
1257 MoFEMFunctionReturn(0);
1258 }
1259
1260 private:
1261 boost::shared_ptr<FaceSideEle> sideFEPtr;
1262 boost::shared_ptr<SideData> sideDataPtr;
1263 SmartPetscObj<Vec> errorSumPtr;
1264 Tag thError;
1265 };
1266
1267 auto simple = mField.getInterface<Simple>();
1268
1269 auto error_sum_ptr = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
1270 Tag th_error;
1271 double def_val = 0;
1272 CHKERR mField.get_moab().tag_get_handle("Error", 1, MB_TYPE_DOUBLE, th_error,
1273 MB_TAG_CREAT | MB_TAG_SPARSE,
1274 &def_val);
1275
1276 auto clear_tags = [&]() {
1277 MoFEMFunctionBegin;
1278 Range fe_ents;
1279 CHKERR mField.get_moab().get_entities_by_dimension(0, FE_DIM, fe_ents);
1280 double zero = 0;
1281 CHKERR mField.get_moab().tag_clear_data(th_error, fe_ents, &zero);
1282 MoFEMFunctionReturn(0);
1283 };
1284
1285 auto evaluate_error = [&]() {
1286 MoFEMFunctionBegin;
1287 auto skel_fe = boost::make_shared<BoundaryEle>(mField);
1288 skel_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
1289 auto side_data_ptr = boost::make_shared<SideData>();
1290 auto side_fe_ptr = getSideFE(side_data_ptr);
1291 skel_fe->getOpPtrVector().push_back(
1292 new OpErrorSkel(side_fe_ptr, side_data_ptr, error_sum_ptr, th_error));
1293 auto simple = mField.getInterface<Simple>();
1294
1295 skel_fe->exeTestHook = [&](FEMethod *fe_ptr) {
1296 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1297 skeleton_bit);
1298 };
1299
1300 CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
1301 simple->getSkeletonFEName(), skel_fe);
1302
1303 MoFEMFunctionReturn(0);
1304 };
1305
1306 auto assemble_and_sum = [](auto vec) {
1307 CHK_THROW_MESSAGE(VecAssemblyBegin(vec), "assemble");
1308 CHK_THROW_MESSAGE(VecAssemblyEnd(vec), "assemble");
1309 double sum;
1310 CHK_THROW_MESSAGE(VecSum(vec, &sum), "assemble");
1311 return sum;
1312 };
1313
1314 auto propagate_error_to_parents = [&]() {
1315 MoFEMFunctionBegin;
1316
1317 auto &moab = mField.get_moab();
1318 auto fe_ptr = boost::make_shared<FEMethod>();
1319 fe_ptr->exeTestHook = [&](FEMethod *fe_ptr) {
1320 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1321 current_bit);
1322 };
1323
1324 fe_ptr->preProcessHook = []() { return 0; };
1325 fe_ptr->postProcessHook = []() { return 0; };
1326 fe_ptr->operatorHook = [&]() {
1327 MoFEMFunctionBegin;
1328
1329 auto fe_ent = fe_ptr->numeredEntFiniteElementPtr->getEnt();
1330 auto parent = fe_ptr->numeredEntFiniteElementPtr->getParentEnt();
1331 auto th_parent = fe_ptr->numeredEntFiniteElementPtr->getBasicDataPtr()
1332 ->th_RefParentHandle;
1333
1334 double error;
1335 CHKERR moab.tag_get_data(th_error, &fe_ent, 1, &error);
1336
1337 boost::function<MoFEMErrorCode(EntityHandle, double)> add_error =
1338 [&](auto fe_ent, auto error) {
1339 MoFEMFunctionBegin;
1340 double *e_ptr;
1341 CHKERR moab.tag_get_by_ptr(th_error, &fe_ent, 1,
1342 (const void **)&e_ptr);
1343 (*e_ptr) += error;
1344
1345 EntityHandle parent;
1346 CHKERR moab.tag_get_data(th_parent, &fe_ent, 1, &parent);
1347 if (parent != fe_ent && parent)
1348 CHKERR add_error(parent, *e_ptr);
1349
1350 MoFEMFunctionReturn(0);
1351 };
1352
1353 CHKERR add_error(parent, error);
1354
1355 MoFEMFunctionReturn(0);
1356 };
1357
1358 CHKERR DMoFEMLoopFiniteElements(simple->getDM(), simple->getDomainFEName(),
1359 fe_ptr);
1360
1361 MoFEMFunctionReturn(0);
1362 };
1363
1364 CHK_THROW_MESSAGE(clear_tags(), "clear error tags");
1365 CHK_THROW_MESSAGE(evaluate_error(), "evaluate error");
1366 CHK_THROW_MESSAGE(propagate_error_to_parents(), "propagate error");
1367
1368 return std::make_tuple(assemble_and_sum(error_sum_ptr), th_error);
1369}
type
std::string type
Definition JsonConfigManager.cpp:28
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition definitions.h:612
skeleton_bit
constexpr int skeleton_bit
skeleton elements bit
Definition level_set.cpp:65
J
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition PetscSmartObj.hpp:204
MoFEM::make_array
constexpr auto make_array(Arg &&...arg)
Create Array.
Definition Templates.hpp:2528
I
constexpr IntegrationType I
Definition operators_tests.cpp:31
LEFT_SIDE
@ LEFT_SIDE
Definition plate.cpp:92
RIGHT_SIDE
@ RIGHT_SIDE
Definition plate.cpp:92
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:995
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::Simple::getBitRefLevel
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
Definition Simple.hpp:415
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85

◆ get_level_set()

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

inital level set, i.e. advected field

Parameters
x
y
z
Returns
double
Examples
mofem/tutorials/adv-3_level_set/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}
sdf_hertz_2d_axisymm_plane.xc
int xc
Definition sdf_hertz_2d_axisymm_plane.py:6
sdf_hertz_2d_axisymm_plane.zc
int zc
Definition sdf_hertz_2d_axisymm_plane.py:8
sdf_hertz_2d_axisymm_plane.yc
int yc
Definition sdf_hertz_2d_axisymm_plane.py:7
sdf.r
int r
Definition sdf.py:205

◆ 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
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

◆ 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
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 995 of file level_set.cpp.

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

◆ 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
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 920 of file level_set.cpp.

920 {
921 auto get_parent_vel_this = [&]() {
922 auto parent_fe_ptr = boost::make_shared<DomainParentEle>(mField);
923 CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
924 parent_fe_ptr->getOpPtrVector(), {potential_velocity_space});
925 parent_fe_ptr->getOpPtrVector().push_back(
926 new OpCalculateHcurlVectorCurl<potential_velocity_field_dim, SPACE_DIM>(
927 "V", vel_ptr));
928 return parent_fe_ptr;
929 };
930
931 auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr) {
932 std::vector<boost::shared_ptr<DomainParentEle>> parents_elems_ptr_vec;
933 for (int l = 0; l <= nb_levels; ++l)
934 parents_elems_ptr_vec.emplace_back(
935 boost::make_shared<DomainParentEle>(mField));
936 for (auto l = 1; l <= nb_levels; ++l) {
937 parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
938 new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
939 BitRefLevel().set(0).flip(), this_fe_ptr,
940 BitRefLevel().set(0), BitRefLevel().set()));
941 }
942 return parents_elems_ptr_vec[0];
943 };
944
945 auto this_fe_ptr = get_parent_vel_this();
946 auto parent_fe_ptr = get_parents_vel_fe_ptr(this_fe_ptr);
947 return new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
948 BitRefLevel().set(0).flip(), this_fe_ptr,
949 BitRefLevel().set(0), BitRefLevel().set());
950}
MoFEM::OpCalculateHcurlVectorCurl
Calculate curl of vector field.
Definition UserDataOperators.hpp:2929

◆ 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
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 1684 of file level_set.cpp.

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

◆ 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
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 1805 of file level_set.cpp.

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

◆ pushOpDomain()

MoFEMErrorCode LevelSet::pushOpDomain ( )
private

push operators to integrate operators on domain

Returns
MoFEMErrorCode
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 952 of file level_set.cpp.

952 {
953 MoFEMFunctionBegin;
954 auto pip = mField.getInterface<PipelineManager>(); // get interface to
955 // pipeline manager
956
957 pip->getOpDomainLhsPipeline().clear();
958 pip->getOpDomainRhsPipeline().clear();
959
960 pip->setDomainLhsIntegrationRule([](int, int, int o) { return 3 * o; });
961 pip->setDomainRhsIntegrationRule([](int, int, int o) { return 3 * o; });
962
963 pip->getDomainLhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
964 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
965 current_bit);
966 };
967 pip->getDomainRhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
968 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
969 current_bit);
970 };
971
972 auto l_ptr = boost::make_shared<MatrixDouble>();
973 auto l_dot_ptr = boost::make_shared<MatrixDouble>();
974 auto vel_ptr = boost::make_shared<MatrixDouble>();
975
976 pip->getOpDomainRhsPipeline().push_back(getZeroLevelVelOp(vel_ptr));
977 CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainRhsPipeline(),
978 {L2});
979 pip->getOpDomainRhsPipeline().push_back(
980 new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_ptr));
981 pip->getOpDomainRhsPipeline().push_back(
982 new OpCalculateTensor2FieldValuesDot<DIM1, DIM2>("L", l_dot_ptr));
983 pip->getOpDomainRhsPipeline().push_back(
984 new OpRhsDomain("L", l_ptr, l_dot_ptr, vel_ptr));
985
986 pip->getOpDomainLhsPipeline().push_back(getZeroLevelVelOp(vel_ptr));
987 CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainLhsPipeline(),
988 {L2});
989 pip->getOpDomainLhsPipeline().push_back(new OpLhsDomain("L", vel_ptr));
990
991 MoFEMFunctionReturn(0);
992}
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition PipelineManager.hpp:1265
MoFEM::OpCalculateTensor2FieldValuesDot
Get time direvarive values at integration pts for tensor field rank 2, i.e. matrix field.
Definition UserDataOperators.hpp:1047

◆ pushOpSkeleton()

MoFEMErrorCode LevelSet::pushOpSkeleton ( )
private

push operator to integrate on skeleton

Returns
MoFEMErrorCode
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 1170 of file level_set.cpp.

1170 {
1171 MoFEMFunctionBegin;
1172 auto pip = mField.getInterface<PipelineManager>(); // get interface to
1173
1174 pip->getOpSkeletonLhsPipeline().clear();
1175 pip->getOpSkeletonRhsPipeline().clear();
1176
1177 pip->setSkeletonLhsIntegrationRule([](int, int, int o) { return 18; });
1178 pip->setSkeletonRhsIntegrationRule([](int, int, int o) { return 18; });
1179
1180 pip->getSkeletonLhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
1181 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1182 skeleton_bit);
1183 };
1184 pip->getSkeletonRhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
1185 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1186 skeleton_bit);
1187 };
1188
1189 auto side_data_ptr = boost::make_shared<SideData>();
1190 auto side_fe_ptr = getSideFE(side_data_ptr);
1191
1192 pip->getOpSkeletonRhsPipeline().push_back(
1193 new OpRhsSkeleton(side_data_ptr, side_fe_ptr));
1194 pip->getOpSkeletonLhsPipeline().push_back(
1195 new OpLhsSkeleton(side_data_ptr, side_fe_ptr));
1196
1197 MoFEMFunctionReturn(0);
1198}
MoFEM::PipelineManager::getOpSkeletonLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpSkeletonLhsPipeline()
Get the Op Skeleton Lhs Pipeline object.
Definition PipelineManager.hpp:1433

◆ readMesh()

MoFEMErrorCode LevelSet::readMesh ( )
private

read mesh

Returns
MoFEMErrorCode
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 502 of file level_set.cpp.

502 {
503 MoFEMFunctionBegin;
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 = [&]() {
520 MoFEMFunctionBegin;
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
552 MoFEMFunctionReturn(0);
553 };
554
555 CHKERR set_problem_bit();
556
557 MoFEMFunctionReturn(0);
558}
start_bit
constexpr int start_bit
Definition level_set.cpp:60
aggregate_bit
constexpr int aggregate_bit
all bits for advection problem
Definition level_set.cpp:66
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180

◆ refineMesh()

MoFEMErrorCode LevelSet::refineMesh ( WrapperClass &&  wp)
private
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 2137 of file level_set.cpp.

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

◆ runProblem()

MoFEMErrorCode LevelSet::runProblem ( )
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 386 of file level_set.cpp.

386 {
387 MoFEMFunctionBegin;
388 CHKERR readMesh();
389 CHKERR setUpProblem();
390
391 if constexpr (debug) {
392 CHKERR testSideFE();
393 CHKERR testOp();
394 }
395
396 CHKERR initialiseFieldVelocity();
397
398 maxPtr = boost::make_shared<double>(0);
399 CHKERR refineMesh(WrapperClassInitalSolution(maxPtr));
400
401 auto simple = mField.getInterface<Simple>();
402 simple->getBitRefLevel() =
403 BitRefLevel().set(skeleton_bit) | BitRefLevel().set(aggregate_bit);
404 simple->getBitRefLevelMask() = BitRefLevel().set();
405 simple->reSetUp(true);
406
407 CHKERR solveAdvection();
408
409 MoFEMFunctionReturn(0);
410}
LevelSet::solveAdvection
MoFEMErrorCode solveAdvection()
solve advection problem
Definition level_set.cpp:1925
LevelSet::maxPtr
boost::shared_ptr< double > maxPtr
Definition level_set.cpp:370
LevelSet::readMesh
MoFEMErrorCode readMesh()
read mesh
Definition level_set.cpp:502
LevelSet::setUpProblem
MoFEMErrorCode setUpProblem()
create fields, and set approximation order
Definition level_set.cpp:560
LevelSet::testOp
MoFEMErrorCode testOp()
test consistency between tangent matrix and the right hand side vectors
Definition level_set.cpp:1593
LevelSet::refineMesh
MoFEMErrorCode refineMesh(WrapperClass &&wp)
Definition level_set.cpp:2137
LevelSet::testSideFE
MoFEMErrorCode testSideFE()
test integration side elements
Definition level_set.cpp:1378
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:1805

◆ setUpProblem()

MoFEMErrorCode LevelSet::setUpProblem ( )
private

create fields, and set approximation order

Returns
MoFEMErrorCode
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 560 of file level_set.cpp.

560 {
561 MoFEMFunctionBegin;
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,
567 AINSWORTH_LEGENDRE_BASE, 1);
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
576 MoFEMFunctionReturn(0);
577}
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
L2
@ L2
field with C-1 continuity
Definition definitions.h:88
potential_velocity_space
constexpr FieldSpace potential_velocity_space
Definition level_set.cpp:49
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition Simple.cpp:261

◆ solveAdvection()

MoFEMErrorCode LevelSet::solveAdvection ( )
private

solve advection problem

Returns
* MoFEMErrorCode
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 1925 of file level_set.cpp.

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

◆ testOp()

MoFEMErrorCode LevelSet::testOp ( )
private

test consistency between tangent matrix and the right hand side vectors

Returns
MoFEMErrorCode
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 1593 of file level_set.cpp.

1593 {
1594 MoFEMFunctionBegin;
1595
1596 // get operators tester
1597 auto simple = mField.getInterface<Simple>();
1598 auto opt = mField.getInterface<OperatorsTester>(); // get interface to
1599 // OperatorsTester
1600 auto pip = mField.getInterface<PipelineManager>(); // get interface to
1601 // pipeline manager
1602
1603 CHKERR pushOpDomain();
1604 CHKERR pushOpSkeleton();
1605
1606 auto post_proc = [&](auto dm, auto f_res, auto out_name) {
1607 MoFEMFunctionBegin;
1608 auto post_proc_fe =
1609 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
1610
1611 if constexpr (DIM1 == 1 && DIM2 == 1) {
1612 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
1613
1614 auto l_vec = boost::make_shared<VectorDouble>();
1615 post_proc_fe->getOpPtrVector().push_back(
1616 new OpCalculateScalarFieldValues("L", l_vec, f_res));
1617
1618 post_proc_fe->getOpPtrVector().push_back(
1619
1620 new OpPPMap(
1621
1622 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1623
1624 {{"L", l_vec}},
1625
1626 {},
1627
1628 {}, {})
1629
1630 );
1631 }
1632
1633 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
1634 post_proc_fe);
1635 post_proc_fe->writeFile(out_name);
1636 MoFEMFunctionReturn(0);
1637 };
1638
1639 constexpr double eps = 1e-4;
1640
1641 auto x =
1642 opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}, {"V", {-1, 1}}});
1643 auto dot_x = opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}});
1644 auto diff_x = opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}});
1645
1646 auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline,
1647 auto rhs_pipeline) {
1648 MoFEMFunctionBegin;
1649
1650 auto diff_res = opt->checkCentralFiniteDifference(
1651 simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
1652 SmartPetscObj<Vec>(), diff_x, 0, 1, eps);
1653
1654 if constexpr (debug) {
1655 // Example how to plot direction in direction diff_x. If instead
1656 // directionalCentralFiniteDifference(...) diff_res is used, then error
1657 // on directive is plotted.
1658 CHKERR post_proc(simple->getDM(), diff_res, "tangent_op_error.h5m");
1659 }
1660
1661 // Calculate norm of difference between directive calculated from finite
1662 // difference, and tangent matrix.
1663 double fnorm;
1664 CHKERR VecNorm(diff_res, NORM_2, &fnorm);
1665 MOFEM_LOG_C("LevelSet", Sev::inform,
1666 "Test consistency of tangent matrix %3.4e", fnorm);
1667
1668 constexpr double err = 1e-9;
1669 if (fnorm > err)
1670 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1671 "Norm of directional derivative too large err = %3.4e", fnorm);
1672
1673 MoFEMFunctionReturn(0);
1674 };
1675
1676 CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
1677 pip->getDomainRhsFE());
1678 CHKERR test_domain_ops(simple->getSkeletonFEName(), pip->getSkeletonLhsFE(),
1679 pip->getSkeletonRhsFE());
1680
1681 MoFEMFunctionReturn(0);
1682};
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition LogManager.hpp:319
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
MoFEM::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition OperatorsTester.hpp:21

◆ 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
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 1378 of file level_set.cpp.

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

Member Data Documentation

◆ maxPtr

boost::shared_ptr<double> LevelSet::maxPtr
private
Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 370 of file level_set.cpp.

◆ mField

MoFEM::Interface& LevelSet::mField
private

integrate skeleton operators on khs

Examples
mofem/tutorials/adv-3_level_set/level_set.cpp.

Definition at line 349 of file level_set.cpp.

The documentation for this struct was generated from the following file:
Generated by Doxygen 1.9.8 and hosted at University of Glasgow