LevelSet Struct Reference

Collaboration diagram for LevelSet:

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
using	OpMassLL
using	OpSourceL
using	OpMassVV
using	OpSourceV
using	OpScalarFieldL
using	AssemblyBoundaryEleOp

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

level_set.cpp.

Definition at line 73 of file level_set.cpp.

Member Typedef Documentation

AssemblyBoundaryEleOp

using LevelSet::AssemblyBoundaryEleOp

private

Initial value:

 
      FormsIntegrators<BoundaryEleOp>::Assembly<A>::OpBase

Definition at line 364 of file level_set.cpp.

AssemblyDomainEleOp

using LevelSet::AssemblyDomainEleOp

private

Initial value:

 
      FormsIntegrators<DomainEleOp>::Assembly<A>::OpBase

Definition at line 351 of file level_set.cpp.

MatSideArray

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

private

Examples

level_set.cpp.

Definition at line 80 of file level_set.cpp.

OpMassLL

using LevelSet::OpMassLL

private

Initial value:

 FormsIntegrators<DomainEleOp>::Assembly<A>::BiLinearForm<
      G>::OpMass<1, DIM1 * DIM2>

Examples

level_set.cpp.

Definition at line 353 of file level_set.cpp.

OpMassVV

using LevelSet::OpMassVV

private

Initial value:

 FormsIntegrators<DomainEleOp>::Assembly<A>::BiLinearForm<
      G>::OpMass<potential_velocity_field_dim, potential_velocity_field_dim>

Examples

level_set.cpp.

Definition at line 357 of file level_set.cpp.

OpScalarFieldL

using LevelSet::OpScalarFieldL

private

Initial value:

 FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm<
      G>::OpBaseTimesVector<1, DIM1 * DIM2, 1>

Examples

level_set.cpp.

Definition at line 361 of file level_set.cpp.

OpSourceL

using LevelSet::OpSourceL

private

Initial value:

 FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm<
      G>::OpSource<1, DIM1 * DIM2>

Examples

level_set.cpp.

Definition at line 355 of file level_set.cpp.

OpSourceV

using LevelSet::OpSourceV

private

Initial value:

 FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm<
      G>::OpSource<potential_velocity_field_dim, potential_velocity_field_dim>

Examples

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.

{ LEFT_SIDE = 0, RIGHT_SIDE = 1 };

Constructor & Destructor Documentation

LevelSet()

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

inline

Examples

level_set.cpp.

Definition at line 75 of file level_set.cpp.

: mField(m_field) {}

Member Function Documentation

dgProjection()

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

private

dg level set projection

Parameters

prj_bit
mesh_bit

Returns

MoFEMErrorCode

Examples

level_set.cpp.

Definition at line 2401 of file level_set.cpp.

                                                              {
  MoFEMFunctionBegin;
 
  // get operators tester
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>(); // get interface to
  auto bit_mng = mField.getInterface<BitRefManager>();
  auto prb_mng = mField.getInterface<ProblemsManager>();
 
  auto lhs_fe = boost::make_shared<DomainEle>(mField);
  auto rhs_fe_prj = boost::make_shared<DomainEle>(mField);
  auto rhs_fe_current = boost::make_shared<DomainEle>(mField);
 
  lhs_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
  rhs_fe_prj->getRuleHook = [](int, int, int o) { return 3 * o; };
  rhs_fe_current->getRuleHook = [](int, int, int o) { return 3 * o; };
 
  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "DG_PROJECTION");
  CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
  CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
  CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
  CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
  CHKERR DMSetUp(sub_dm);
 
  Range current_ents; // ents used to do calculations
  CHKERR bit_mng->getEntitiesByDimAndRefLevel(BitRefLevel().set(current_bit),
                                              BitRefLevel().set(), FE_DIM,
                                              current_ents);
  Range prj_ents; // ents from which data are projected
  CHKERR bit_mng->getEntitiesByDimAndRefLevel(
      BitRefLevel().set(projection_bit), BitRefLevel().set(), FE_DIM, prj_ents);
  for (auto l = 0; l != nb_levels; ++l) {
    CHKERR bit_mng->updateRangeByParent(prj_ents, prj_ents);
  }
  current_ents = subtract(
      current_ents, prj_ents); // only restric to entities needed projection
 
  auto test_mesh_bit = [&](FEMethod *fe_ptr) {
    return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
        current_bit);
  };
  auto test_prj_bit = [&](FEMethod *fe_ptr) {
    return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
        projection_bit);
  };
  auto test_current_bit = [&](FEMethod *fe_ptr) {
    return current_ents.find(fe_ptr->getFEEntityHandle()) != current_ents.end();
  };
 
  lhs_fe->exeTestHook =
      test_mesh_bit; // that element only is run when current bit is set
  rhs_fe_prj->exeTestHook =
      test_prj_bit; // that element is run only when projection bit is set
  rhs_fe_current->exeTestHook =
      test_current_bit; // that element is only run when current bit is set
 
  BitRefLevel remove_mask = BitRefLevel().set(current_bit);
  remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
  CHKERR prb_mng->removeDofsOnEntities(
      "DG_PROJECTION", "L", BitRefLevel().set(), remove_mask, nullptr, 0,
      MAX_DOFS_ON_ENTITY, 0, 100, NOISY,
      true); // remove all DOFs which are not
             // on current bit. This case works for L2 space
 
  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(lhs_fe->getOpPtrVector(),
                                                  {L2});
  lhs_fe->getOpPtrVector().push_back(
      new OpMassLL("L", "L")); // Assemble projection matrix
 
  // This assumes that projection mesh is finer, current mesh is coarsened.
  auto set_prj_from_child = [&](auto rhs_fe_prj) {
    CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
        rhs_fe_prj->getOpPtrVector(), {L2});
 
    // Evaluate field value on projection mesh
    auto l_vec = boost::make_shared<MatrixDouble>();
    rhs_fe_prj->getOpPtrVector().push_back(
        new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_vec));
 
    // This element is used to assemble
    auto get_parent_this = [&]() {
      auto fe_parent_this = boost::make_shared<DomianParentEle>(mField);
      fe_parent_this->getOpPtrVector().push_back(
          new OpScalarFieldL("L", l_vec));
      return fe_parent_this;
    };
 
    // Create levels of parent elements, until current element is reached, and
    // then assemble.
    auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
      std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
      for (int l = 0; l <= nb_levels; ++l)
        parents_elems_ptr_vec.emplace_back(
            boost::make_shared<DomianParentEle>(mField));
      for (auto l = 1; l <= nb_levels; ++l) {
        parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
            new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
                            BitRefLevel().set(current_bit).flip(), this_fe_ptr,
                            BitRefLevel().set(current_bit),
                            BitRefLevel().set()));
      }
      return parents_elems_ptr_vec[0];
    };
 
    auto this_fe_ptr = get_parent_this();
    auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
    rhs_fe_prj->getOpPtrVector().push_back(
        new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
                        BitRefLevel().set(current_bit).flip(), this_fe_ptr,
                        BitRefLevel().set(current_bit), BitRefLevel().set()));
  };
 
  // This assumed that current mesh is refined, and projection mesh is coarser
  auto set_prj_from_parent = [&](auto rhs_fe_current) {
 
    // Evaluate field value on projection mesh
    auto l_vec = boost::make_shared<MatrixDouble>();
 
    // Evaluate field on coarser element
    auto get_parent_this = [&]() {
      auto fe_parent_this = boost::make_shared<DomianParentEle>(mField);
      fe_parent_this->getOpPtrVector().push_back(
          new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_vec));
      return fe_parent_this;
    };
 
    // Create stack of evaluation on parent elements
    auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
      std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
      for (int l = 0; l <= nb_levels; ++l)
        parents_elems_ptr_vec.emplace_back(
            boost::make_shared<DomianParentEle>(mField));
      for (auto l = 1; l <= nb_levels; ++l) {
        parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
            new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
                            BitRefLevel().set(projection_bit).flip(),
                            this_fe_ptr, BitRefLevel().set(projection_bit),
                            BitRefLevel().set()));
      }
      return parents_elems_ptr_vec[0];
    };
 
    auto this_fe_ptr = get_parent_this();
    auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
 
    auto reset_op_ptr = new DomainEleOp(NOSPACE, DomainEleOp::OPSPACE);
    reset_op_ptr->doWorkRhsHook = [&](DataOperator *op_ptr, int, EntityType,
                                      EntData &) {
      l_vec->resize(static_cast<DomainEleOp *>(op_ptr)->getGaussPts().size2(),
                    false);
      l_vec->clear();
      return 0;
    };
    rhs_fe_current->getOpPtrVector().push_back(reset_op_ptr);
    rhs_fe_current->getOpPtrVector().push_back(
        new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
                        BitRefLevel().set(projection_bit).flip(), this_fe_ptr,
                        BitRefLevel(), BitRefLevel()));
 
    // At the end assemble of current finite element
    rhs_fe_current->getOpPtrVector().push_back(new OpScalarFieldL("L", l_vec));
  };
 
  set_prj_from_child(rhs_fe_prj);
  set_prj_from_parent(rhs_fe_current);
 
  boost::shared_ptr<FEMethod> null_fe;
  getDMKspCtx(sub_dm)->clearLoops();
  CHKERR DMMoFEMKSPSetComputeOperators(sub_dm, simple->getDomainFEName(),
                                       lhs_fe, null_fe, null_fe);
  CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(), rhs_fe_prj,
                                 null_fe, null_fe);
  CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(),
                                 rhs_fe_current, null_fe, null_fe);
  auto ksp = MoFEM::createKSP(mField.get_comm());
  CHKERR KSPSetDM(ksp, sub_dm);
 
  CHKERR KSPSetDM(ksp, sub_dm);
  CHKERR KSPSetFromOptions(ksp);
  CHKERR KSPSetUp(ksp);
 
  auto L = createDMVector(sub_dm);
  auto F = vectorDuplicate(L);
 
  CHKERR KSPSolve(ksp, F, L);
  CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
 
  auto [error, th_error] = evaluateError();
  MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
 
  auto post_proc = [&](auto dm, auto out_name, auto th_error) {
    MoFEMFunctionBegin;
    auto post_proc_fe =
        boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
    post_proc_fe->setTagsToTransfer({th_error});
    post_proc_fe->exeTestHook = test_mesh_bit;
 
    if constexpr (DIM1 == 1 && DIM2 == 1) {
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
      auto l_vec = boost::make_shared<VectorDouble>();
      auto l_grad_mat = boost::make_shared<MatrixDouble>();
      CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
          post_proc_fe->getOpPtrVector(), {L2});
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateScalarFieldValues("L", l_vec));
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateScalarFieldGradient<SPACE_DIM>("L", l_grad_mat));
 
      post_proc_fe->getOpPtrVector().push_back(
 
          new OpPPMap(
 
              post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
              {{"L", l_vec}},
 
              {{"GradL", l_grad_mat}},
 
              {}, {})
 
      );
    }
 
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                    post_proc_fe);
    post_proc_fe->writeFile(out_name);
    MoFEMFunctionReturn(0);
  };
 
  if constexpr (debug)
    CHKERR post_proc(sub_dm, "dg_projection.h5m", th_error);
 
  MoFEMFunctionReturn(0);
}

evaluateError()

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

private

evaluate error

Returns

MoFEMErrorCode

Examples

level_set.cpp.

Definition at line 1203 of file level_set.cpp.

                                              {
 
  struct OpErrorSkel : BoundaryEleOp {
 
    OpErrorSkel(boost::shared_ptr<FaceSideEle> side_fe_ptr,
                boost::shared_ptr<SideData> side_data_ptr,
                SmartPetscObj<Vec> error_sum_ptr, Tag th_error)
        : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
          sideFEPtr(side_fe_ptr), sideDataPtr(side_data_ptr),
          errorSumPtr(error_sum_ptr), thError(th_error) {}
 
    MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
      MoFEMFunctionBegin;
 
      // Collect data from side domain elements
      CHKERR loopSideFaces("dFE", *sideFEPtr);
      const auto in_the_loop =
          sideFEPtr->nInTheLoop; // return number of elements on the side
 
      auto not_side = [](auto s) {
        return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
      };
 
      auto nb_gauss_pts = getGaussPts().size2();
 
      for (auto s : {LEFT_SIDE, RIGHT_SIDE}) {
 
        auto arr_t_l = make_array(
            getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[LEFT_SIDE]),
            getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[RIGHT_SIDE]));
        auto arr_t_vel = make_array(
            getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[LEFT_SIDE]),
            getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[RIGHT_SIDE]));
 
        auto next = [&]() {
          for (auto &t_l : arr_t_l)
            ++t_l;
          for (auto &t_vel : arr_t_vel)
            ++t_vel;
        };
 
        double e = 0;
        auto t_w = getFTensor0IntegrationWeight();
        for (int gg = 0; gg != nb_gauss_pts; ++gg) {
          FTensor::Tensor2<double, DIM1, DIM2> t_diff;
          t_diff(I, J) = arr_t_l[LEFT_SIDE](I, J) - arr_t_l[RIGHT_SIDE](I, J);
          e += t_w * getMeasure() * t_diff(I, J) * t_diff(I, J);
          next();
          ++t_w;
        }
        e = std::sqrt(e);
 
        moab::Interface &moab =
            getNumeredEntFiniteElementPtr()->getBasicDataPtr()->moab;
        const void *tags_ptr[2];
        CHKERR moab.tag_get_by_ptr(thError, sideDataPtr->feSideHandle.data(), 2,
                                   tags_ptr);
        for (auto ff : {0, 1}) {
          *((double *)tags_ptr[ff]) += e;
        }
        CHKERR VecSetValue(errorSumPtr, 0, e, ADD_VALUES);
      };
 
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<FaceSideEle> sideFEPtr;
    boost::shared_ptr<SideData> sideDataPtr;
    SmartPetscObj<Vec> errorSumPtr;
    Tag thError;
  };
 
  auto simple = mField.getInterface<Simple>();
 
  auto error_sum_ptr = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
  Tag th_error;
  double def_val = 0;
  CHKERR mField.get_moab().tag_get_handle("Error", 1, MB_TYPE_DOUBLE, th_error,
                                          MB_TAG_CREAT | MB_TAG_SPARSE,
                                          &def_val);
 
  auto clear_tags = [&]() {
    MoFEMFunctionBegin;
    Range fe_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, FE_DIM, fe_ents);
    double zero;
    CHKERR mField.get_moab().tag_clear_data(th_error, fe_ents, &zero);
    MoFEMFunctionReturn(0);
  };
 
  auto evaluate_error = [&]() {
    MoFEMFunctionBegin;
    auto skel_fe = boost::make_shared<BoundaryEle>(mField);
    skel_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
    auto side_data_ptr = boost::make_shared<SideData>();
    auto side_fe_ptr = getSideFE(side_data_ptr);
    skel_fe->getOpPtrVector().push_back(
        new OpErrorSkel(side_fe_ptr, side_data_ptr, error_sum_ptr, th_error));
    auto simple = mField.getInterface<Simple>();
 
    skel_fe->exeTestHook = [&](FEMethod *fe_ptr) {
      return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
          skeleton_bit);
    };
 
    CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
                                    simple->getSkeletonFEName(), skel_fe);
 
    MoFEMFunctionReturn(0);
  };
 
  auto assemble_and_sum = [](auto vec) {
    CHK_THROW_MESSAGE(VecAssemblyBegin(vec), "assemble");
    CHK_THROW_MESSAGE(VecAssemblyEnd(vec), "assemble");
    double sum;
    CHK_THROW_MESSAGE(VecSum(vec, &sum), "assemble");
    return sum;
  };
 
  auto propagate_error_to_parents = [&]() {
    MoFEMFunctionBegin;
 
    auto &moab = mField.get_moab();
    auto fe_ptr = boost::make_shared<FEMethod>();
    fe_ptr->exeTestHook = [&](FEMethod *fe_ptr) {
      return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
          current_bit);
    };
 
    fe_ptr->preProcessHook = []() { return 0; };
    fe_ptr->postProcessHook = []() { return 0; };
    fe_ptr->operatorHook = [&]() {
      MoFEMFunctionBegin;
 
      auto fe_ent = fe_ptr->numeredEntFiniteElementPtr->getEnt();
      auto parent = fe_ptr->numeredEntFiniteElementPtr->getParentEnt();
      auto th_parent = fe_ptr->numeredEntFiniteElementPtr->getBasicDataPtr()
                           ->th_RefParentHandle;
 
      double error;
      CHKERR moab.tag_get_data(th_error, &fe_ent, 1, &error);
 
      boost::function<MoFEMErrorCode(EntityHandle, double)> add_error =
          [&](auto fe_ent, auto error) {
            MoFEMFunctionBegin;
            double *e_ptr;
            CHKERR moab.tag_get_by_ptr(th_error, &fe_ent, 1,
                                       (const void **)&e_ptr);
            (*e_ptr) += error;
 
            EntityHandle parent;
            CHKERR moab.tag_get_data(th_parent, &fe_ent, 1, &parent);
            if (parent != fe_ent && parent)
              CHKERR add_error(parent, *e_ptr);
 
            MoFEMFunctionReturn(0);
          };
 
      CHKERR add_error(parent, error);
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR DMoFEMLoopFiniteElements(simple->getDM(), simple->getDomainFEName(),
                                    fe_ptr);
 
    MoFEMFunctionReturn(0);
  };
 
  CHK_THROW_MESSAGE(clear_tags(), "clear error tags");
  CHK_THROW_MESSAGE(evaluate_error(), "evaluate error");
  CHK_THROW_MESSAGE(propagate_error_to_parents(), "propagate error");
 
  return std::make_tuple(assemble_and_sum(error_sum_ptr), th_error);
}

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

level_set.cpp.

Definition at line 378 of file level_set.cpp.

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

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

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.

                                                                       {
  return (x * x - 0.25) * (y * y - 0.25);
}

getSideFE()

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

private

create side element to assemble data from sides

Parameters

side_data_ptr

Returns

boost::shared_ptr<FaceSideEle>

Examples

level_set.cpp.

Definition at line 997 of file level_set.cpp.

                                                           {
 
  auto simple = mField.getInterface<Simple>();
 
  auto l_ptr = boost::make_shared<MatrixDouble>();
  auto vel_ptr = boost::make_shared<MatrixDouble>();
 
  struct OpSideData : public FaceSideEleOp {
    OpSideData(boost::shared_ptr<SideData> side_data_ptr)
        : FaceSideEleOp("L", "L", FaceSideEleOp::OPROWCOL),
          sideDataPtr(side_data_ptr) {
      std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
      for (auto t = moab::CN::TypeDimensionMap[FE_DIM].first;
           t <= moab::CN::TypeDimensionMap[FE_DIM].second; ++t)
        doEntities[t] = true;
      sYmm = false;
    }
 
    MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type, EntData &row_data,
                          EntData &col_data) {
      MoFEMFunctionBegin;
      if ((CN::Dimension(row_type) == FE_DIM) &&
          (CN::Dimension(col_type) == FE_DIM)) {
 
        auto reset = [&](auto nb_in_loop) {
          sideDataPtr->feSideHandle[nb_in_loop] = 0;
          sideDataPtr->indicesRowSideMap[nb_in_loop].clear();
          sideDataPtr->indicesColSideMap[nb_in_loop].clear();
          sideDataPtr->rowBaseSideMap[nb_in_loop].clear();
          sideDataPtr->colBaseSideMap[nb_in_loop].clear();
          sideDataPtr->senseMap[nb_in_loop] = 0;
        };
 
        const auto nb_in_loop = getFEMethod()->nInTheLoop;
        if (nb_in_loop == 0)
          for (auto s : {0, 1})
            reset(s);
 
        sideDataPtr->currentFESide = nb_in_loop;
        sideDataPtr->senseMap[nb_in_loop] = getSkeletonSense();
 
      } else {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Should not happen");
      }
 
      MoFEMFunctionReturn(0);
    };
 
  private:
    boost::shared_ptr<SideData> sideDataPtr;
  };
 
  struct OpSideDataOnParent : public DomainEleOp {
 
    OpSideDataOnParent(boost::shared_ptr<SideData> side_data_ptr,
                       boost::shared_ptr<MatrixDouble> l_ptr,
                       boost::shared_ptr<MatrixDouble> vel_ptr)
        : DomainEleOp("L", "L", DomainEleOp::OPROWCOL),
          sideDataPtr(side_data_ptr), lPtr(l_ptr), velPtr(vel_ptr) {
      std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
      for (auto t = moab::CN::TypeDimensionMap[FE_DIM].first;
           t <= moab::CN::TypeDimensionMap[FE_DIM].second; ++t)
        doEntities[t] = true;
      sYmm = false;
    }
 
    MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type, EntData &row_data,
                          EntData &col_data) {
      MoFEMFunctionBegin;
 
      if ((CN::Dimension(row_type) == FE_DIM) &&
          (CN::Dimension(col_type) == FE_DIM)) {
        const auto nb_in_loop = sideDataPtr->currentFESide;
        sideDataPtr->feSideHandle[nb_in_loop] = getFEEntityHandle();
        sideDataPtr->indicesRowSideMap[nb_in_loop] = row_data.getIndices();
        sideDataPtr->indicesColSideMap[nb_in_loop] = col_data.getIndices();
        sideDataPtr->rowBaseSideMap[nb_in_loop] = row_data.getN();
        sideDataPtr->colBaseSideMap[nb_in_loop] = col_data.getN();
        (sideDataPtr->lVec)[nb_in_loop] = *lPtr;
        (sideDataPtr->velMat)[nb_in_loop] = *velPtr;
 
#ifndef NDEBUG
        if ((sideDataPtr->lVec)[nb_in_loop].size2() !=
            (sideDataPtr->velMat)[nb_in_loop].size2())
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                   "Wrong number of integaration pts %d != %d",
                   (sideDataPtr->lVec)[nb_in_loop].size2(),
                   (sideDataPtr->velMat)[nb_in_loop].size2());
        if ((sideDataPtr->velMat)[nb_in_loop].size1() != SPACE_DIM)
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                   "Wrong size of velocity vector size = %d",
                   (sideDataPtr->velMat)[nb_in_loop].size1());
#endif
 
        if (!nb_in_loop) {
          (sideDataPtr->lVec)[1] = sideDataPtr->lVec[0];
          (sideDataPtr->velMat)[1] = (sideDataPtr->velMat)[0];
        } else {
#ifndef NDEBUG
          if (sideDataPtr->rowBaseSideMap[0].size1() !=
              sideDataPtr->rowBaseSideMap[1].size1()) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                     "Wrong number of integration pt %d != %d",
                     sideDataPtr->rowBaseSideMap[0].size1(),
                     sideDataPtr->rowBaseSideMap[1].size1());
          }
          if (sideDataPtr->colBaseSideMap[0].size1() !=
              sideDataPtr->colBaseSideMap[1].size1()) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                    "Wrong number of integration pt");
          }
#endif
        }
 
      } else {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Should not happen");
      }
 
      MoFEMFunctionReturn(0);
    };
 
  private:
    boost::shared_ptr<SideData> sideDataPtr;
    boost::shared_ptr<MatrixDouble> lPtr;
    boost::shared_ptr<MatrixDouble> velPtr;
  };
 
  // Calculate fields on param mesh bit element
  auto get_parent_this = [&]() {
    auto parent_fe_ptr = boost::make_shared<DomianParentEle>(mField);
    CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
        parent_fe_ptr->getOpPtrVector(), {L2});
    parent_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_ptr));
    parent_fe_ptr->getOpPtrVector().push_back(
        new OpSideDataOnParent(side_data_ptr, l_ptr, vel_ptr));
    return parent_fe_ptr;
  };
 
  auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
    std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
    for (int l = 0; l <= nb_levels; ++l)
      parents_elems_ptr_vec.emplace_back(
          boost::make_shared<DomianParentEle>(mField));
    for (auto l = 1; l <= nb_levels; ++l) {
      parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
          new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
                          BitRefLevel().set(current_bit).flip(), this_fe_ptr,
                          BitRefLevel().set(current_bit), BitRefLevel().set()));
    }
    return parents_elems_ptr_vec[0];
  };
 
  // Create aliased shared pointers, all elements are destroyed if side_fe_ptr
  // is destroyed
  auto get_side_fe_ptr = [&]() {
    auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
 
    auto this_fe_ptr = get_parent_this();
    auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
 
    side_fe_ptr->getOpPtrVector().push_back(new OpSideData(side_data_ptr));
    side_fe_ptr->getOpPtrVector().push_back(getZeroLevelVelOp(vel_ptr));
    side_fe_ptr->getOpPtrVector().push_back(
        new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
                        BitRefLevel().set(current_bit).flip(), this_fe_ptr,
                        BitRefLevel().set(current_bit), BitRefLevel().set()));
 
    return side_fe_ptr;
  };
 
  return get_side_fe_ptr();
};

getZeroLevelVelOp()

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

private

Get operator calculating velocity on coarse mesh.

Parameters

vel_ptr

Returns

DomainEleOp*

Examples

level_set.cpp.

Definition at line 922 of file level_set.cpp.

                                                                 {
  auto get_parent_vel_this = [&]() {
    auto parent_fe_ptr = boost::make_shared<DomianParentEle>(mField);
    CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
        parent_fe_ptr->getOpPtrVector(), {potential_velocity_space});
    parent_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateHcurlVectorCurl<potential_velocity_field_dim, SPACE_DIM>(
            "V", vel_ptr));
    return parent_fe_ptr;
  };
 
  auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr) {
    std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
    for (int l = 0; l <= nb_levels; ++l)
      parents_elems_ptr_vec.emplace_back(
          boost::make_shared<DomianParentEle>(mField));
    for (auto l = 1; l <= nb_levels; ++l) {
      parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
          new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
                          BitRefLevel().set(0).flip(), this_fe_ptr,
                          BitRefLevel().set(0), BitRefLevel().set()));
    }
    return parents_elems_ptr_vec[0];
  };
 
  auto this_fe_ptr = get_parent_vel_this();
  auto parent_fe_ptr = get_parents_vel_fe_ptr(this_fe_ptr);
  return new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
                         BitRefLevel().set(0).flip(), this_fe_ptr,
                         BitRefLevel().set(0), BitRefLevel().set());
}

initialiseFieldLevelSet()

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

private

initialise field set

Parameters

level_fun

Returns

MoFEMErrorCode

Examples

level_set.cpp.

Definition at line 1693 of file level_set.cpp.

                                                             {
  MoFEMFunctionBegin;
 
  // get operators tester
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>(); // get interface to
                                                     // pipeline manager
  auto prb_mng = mField.getInterface<ProblemsManager>();
 
  boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(mField);
  boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(mField);
  auto swap_fe = [&]() {
    lhs_fe.swap(pip->getDomainLhsFE());
    rhs_fe.swap(pip->getDomainRhsFE());
  };
  swap_fe();
 
  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 3 * o; });
  pip->setDomainRhsIntegrationRule([](int, int, int o) { return 3 * o; });
 
  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "LEVELSET_POJECTION");
  CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
  CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
  CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
  CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
  CHKERR DMSetUp(sub_dm);
 
  BitRefLevel remove_mask = BitRefLevel().set(current_bit);
  remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
  CHKERR prb_mng->removeDofsOnEntities("LEVELSET_POJECTION", "L",
                                       BitRefLevel().set(), remove_mask);
  auto test_bit_ele = [&](FEMethod *fe_ptr) {
    return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
        current_bit);
  };
  pip->getDomainLhsFE()->exeTestHook = test_bit_ele;
  pip->getDomainRhsFE()->exeTestHook = test_bit_ele;
 
  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainRhsPipeline(),
                                                  {L2});
  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainLhsPipeline(),
                                                  {L2});
  pip->getOpDomainLhsPipeline().push_back(new OpMassLL("L", "L"));
  pip->getOpDomainRhsPipeline().push_back(new OpSourceL("L", level_fun));
 
  CHKERR mField.getInterface<FieldBlas>()->setField(0, "L");
 
  auto ksp = pip->createKSP(sub_dm);
  CHKERR KSPSetDM(ksp, sub_dm);
  CHKERR KSPSetFromOptions(ksp);
  CHKERR KSPSetUp(ksp);
 
  auto L = createDMVector(sub_dm);
  auto F = vectorDuplicate(L);
 
  CHKERR KSPSolve(ksp, F, L);
  CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
 
  auto [error, th_error] = evaluateError();
  MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
#ifndef NDEBUG
  auto fe_meshset =
      mField.get_finite_element_meshset(simple->getDomainFEName());
  std::vector<Tag> tags{th_error};
  CHKERR mField.get_moab().write_file("error.h5m", "MOAB",
                                      "PARALLEL=WRITE_PART", &fe_meshset, 1,
                                      &*tags.begin(), tags.size());
#endif
 
  auto post_proc = [&](auto dm, auto out_name, auto th_error) {
    MoFEMFunctionBegin;
    auto post_proc_fe =
        boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
    post_proc_fe->setTagsToTransfer({th_error});
    post_proc_fe->exeTestHook = test_bit_ele;
 
    if constexpr (DIM1 == 1 && DIM2 == 1) {
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      auto l_vec = boost::make_shared<VectorDouble>();
      auto l_grad_mat = boost::make_shared<MatrixDouble>();
      CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
          post_proc_fe->getOpPtrVector(), {L2});
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateScalarFieldValues("L", l_vec));
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateScalarFieldGradient<SPACE_DIM>("L", l_grad_mat));
 
      post_proc_fe->getOpPtrVector().push_back(
 
          new OpPPMap(
 
              post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
              {{"L", l_vec}},
 
              {{"GradL", l_grad_mat}},
 
              {}, {})
 
      );
    }
 
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                    post_proc_fe);
    post_proc_fe->writeFile(out_name);
    MoFEMFunctionReturn(0);
  };
 
  if constexpr (debug)
    CHKERR post_proc(sub_dm, "initial_level_set.h5m", th_error);
 
  swap_fe();
 
  MoFEMFunctionReturn(0);
}

initialiseFieldVelocity()

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

private

initialise potential velocity field

Parameters

vel_fun

Returns

MoFEMErrorCode

Examples

level_set.cpp.

Definition at line 1814 of file level_set.cpp.

                                                           {
  MoFEMFunctionBegin;
 
  // get operators tester
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>(); // get interface to
                                                     // pipeline manager
  auto prb_mng = mField.getInterface<ProblemsManager>();
 
  boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(mField);
  boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(mField);
  auto swap_fe = [&]() {
    lhs_fe.swap(pip->getDomainLhsFE());
    rhs_fe.swap(pip->getDomainRhsFE());
  };
  swap_fe();
 
  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 3 * o; });
  pip->setDomainRhsIntegrationRule([](int, int, int o) { return 3 * o; });
 
  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "VELOCITY_PROJECTION");
  CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
  CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
  CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
  CHKERR DMMoFEMAddSubFieldRow(sub_dm, "V");
  CHKERR DMSetUp(sub_dm);
 
  // Velocities are calculated only on corse mesh
  BitRefLevel remove_mask = BitRefLevel().set(0);
  remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
  CHKERR prb_mng->removeDofsOnEntities("VELOCITY_PROJECTION", "V",
                                       BitRefLevel().set(), remove_mask);
 
  auto test_bit = [&](FEMethod *fe_ptr) {
    return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(0);
  };
  pip->getDomainLhsFE()->exeTestHook = test_bit;
  pip->getDomainRhsFE()->exeTestHook = test_bit;
 
  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainLhsPipeline(),
                                                  {potential_velocity_space});
  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainRhsPipeline(),
                                                  {potential_velocity_space});
 
  pip->getOpDomainLhsPipeline().push_back(new OpMassVV("V", "V"));
  pip->getOpDomainRhsPipeline().push_back(new OpSourceV("V", vel_fun));
 
  auto ksp = pip->createKSP(sub_dm);
  CHKERR KSPSetDM(ksp, sub_dm);
  CHKERR KSPSetFromOptions(ksp);
  CHKERR KSPSetUp(ksp);
 
  auto L = createDMVector(sub_dm);
  auto F = vectorDuplicate(L);
 
  CHKERR KSPSolve(ksp, F, L);
  CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
 
  auto post_proc = [&](auto dm, auto out_name) {
    MoFEMFunctionBegin;
    auto post_proc_fe =
        boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
    post_proc_fe->exeTestHook = test_bit;
 
    if constexpr (FE_DIM == 2) {
 
      CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
          post_proc_fe->getOpPtrVector(), {potential_velocity_space});
 
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      auto potential_vec = boost::make_shared<VectorDouble>();
      auto velocity_mat = boost::make_shared<MatrixDouble>();
 
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateScalarFieldValues("V", potential_vec));
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateHcurlVectorCurl<potential_velocity_field_dim,
                                         SPACE_DIM>("V", velocity_mat));
 
      post_proc_fe->getOpPtrVector().push_back(
 
          new OpPPMap(
 
              post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
              {{"VelocityPotential", potential_vec}},
 
              {{"Velocity", velocity_mat}},
 
              {}, {})
 
      );
 
    } else {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "3d case not implemented");
    }
 
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                    post_proc_fe);
    post_proc_fe->writeFile(out_name);
    MoFEMFunctionReturn(0);
  };
 
  if constexpr (debug)
    CHKERR post_proc(sub_dm, "initial_velocity_potential.h5m");
 
  swap_fe();
 
  MoFEMFunctionReturn(0);
}

pushOpDomain()

MoFEMErrorCode LevelSet::pushOpDomain ( )

private

push operators to integrate operators on domain

Returns

MoFEMErrorCode

Examples

level_set.cpp.

Definition at line 954 of file level_set.cpp.

                                      {
  MoFEMFunctionBegin;
  auto pip = mField.getInterface<PipelineManager>(); // get interface to
                                                     // pipeline manager
 
  pip->getOpDomainLhsPipeline().clear();
  pip->getOpDomainRhsPipeline().clear();
 
  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 3 * o; });
  pip->setDomainRhsIntegrationRule([](int, int, int o) { return 3 * o; });
 
  pip->getDomainLhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
    return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
        current_bit);
  };
  pip->getDomainRhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
    return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
        current_bit);
  };
 
  auto l_ptr = boost::make_shared<MatrixDouble>();
  auto l_dot_ptr = boost::make_shared<MatrixDouble>();
  auto vel_ptr = boost::make_shared<MatrixDouble>();
 
  pip->getOpDomainRhsPipeline().push_back(getZeroLevelVelOp(vel_ptr));
  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainRhsPipeline(),
                                                  {L2});
  pip->getOpDomainRhsPipeline().push_back(
      new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_ptr));
  pip->getOpDomainRhsPipeline().push_back(
      new OpCalculateTensor2FieldValuesDot<DIM1, DIM2>("L", l_dot_ptr));
  pip->getOpDomainRhsPipeline().push_back(
      new OpRhsDomain("L", l_ptr, l_dot_ptr, vel_ptr));
 
  pip->getOpDomainLhsPipeline().push_back(getZeroLevelVelOp(vel_ptr));
  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainLhsPipeline(),
                                                  {L2});
  pip->getOpDomainLhsPipeline().push_back(new OpLhsDomain("L", vel_ptr));
 
  MoFEMFunctionReturn(0);
}

pushOpSkeleton()

MoFEMErrorCode LevelSet::pushOpSkeleton ( )

private

push operator to integrate on skeleton

Returns

MoFEMErrorCode

Examples

level_set.cpp.

Definition at line 1173 of file level_set.cpp.

                                        {
  MoFEMFunctionBegin;
  auto pip = mField.getInterface<PipelineManager>(); // get interface to
 
  pip->getOpSkeletonLhsPipeline().clear();
  pip->getOpSkeletonRhsPipeline().clear();
 
  pip->setSkeletonLhsIntegrationRule([](int, int, int o) { return 18; });
  pip->setSkeletonRhsIntegrationRule([](int, int, int o) { return 18; });
 
  pip->getSkeletonLhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
    return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
        skeleton_bit);
  };
  pip->getSkeletonRhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
    return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
        skeleton_bit);
  };
 
  auto side_data_ptr = boost::make_shared<SideData>();
  auto side_fe_ptr = getSideFE(side_data_ptr);
 
  pip->getOpSkeletonRhsPipeline().push_back(
      new OpRhsSkeleton(side_data_ptr, side_fe_ptr));
  pip->getOpSkeletonLhsPipeline().push_back(
      new OpLhsSkeleton(side_data_ptr, side_fe_ptr));
 
  MoFEMFunctionReturn(0);
}

readMesh()

MoFEMErrorCode LevelSet::readMesh ( )

private

read mesh

Returns

MoFEMErrorCode

Examples

level_set.cpp.

Definition at line 502 of file level_set.cpp.

                                  {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  // get options from command line
  CHKERR simple->getOptions();
 
  // Only L2 field is set in this example. Two lines bellow forces simple
  // interface to creat lower dimension (edge) elements, despite that fact that
  // there is no field spanning on such elements. We need them for DG method.
  simple->getAddSkeletonFE() = true;
  simple->getAddBoundaryFE() = true;
 
  // load mesh file
  simple->getBitRefLevel() = BitRefLevel();
  CHKERR simple->loadFile();
 
  // Initialise bit ref levels
  auto set_problem_bit = [&]() {
    MoFEMFunctionBegin;
    auto bit0 = BitRefLevel().set(start_bit);
    BitRefLevel start_mask;
    for (auto s = 0; s != start_bit; ++s)
      start_mask[s] = true;
 
    auto bit_mng = mField.getInterface<BitRefManager>();
 
    Range level0;
    CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(0),
                                          BitRefLevel().set(), level0);
 
    CHKERR bit_mng->setNthBitRefLevel(level0, current_bit, true);
    CHKERR bit_mng->setNthBitRefLevel(level0, aggregate_bit, true);
    CHKERR bit_mng->setNthBitRefLevel(level0, skeleton_bit, true);
 
    // Set bits to build adjacencies between parents and children. That is used
    // by simple interface.
    simple->getBitAdjEnt() = BitRefLevel().set();
    simple->getBitAdjParent() = BitRefLevel().set();
    simple->getBitRefLevel() = BitRefLevel().set(current_bit);
    simple->getBitRefLevelMask() = BitRefLevel().set();
 
#ifndef NDEBUG
    if constexpr (debug) {
      auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
      CHKERR bit_mng->writeBitLevelByDim(
          BitRefLevel().set(0), BitRefLevel().set(), FE_DIM,
          (proc_str + "level_base.vtk").c_str(), "VTK", "");
    }
#endif
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_problem_bit();
 
  MoFEMFunctionReturn(0);
}

refineMesh()

MoFEMErrorCode LevelSet::refineMesh ( WrapperClass && wp )

private

Examples

level_set.cpp.

Definition at line 2147 of file level_set.cpp.

                                                     {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  auto bit_mng = mField.getInterface<BitRefManager>();
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
  auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
 
  auto set_bit = [](auto l) { return BitRefLevel().set(l); };
 
  auto save_range = [&](const std::string name, const Range &r) {
    MoFEMFunctionBegin;
    auto meshset_ptr = get_temp_meshset_ptr(mField.get_moab());
    CHKERR mField.get_moab().add_entities(*meshset_ptr, r);
    CHKERR mField.get_moab().write_file(name.c_str(), "VTK", "",
                                        meshset_ptr->get_ptr(), 1);
    MoFEMFunctionReturn(0);
  };
 
  // select domain elements to refine by threshold
  auto get_refined_elements_meshset = [&](auto bit, auto mask) {
    Range fe_ents;
    CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit, mask, FE_DIM, fe_ents);
 
    Tag th_error;
    CHK_MOAB_THROW(mField.get_moab().tag_get_handle("Error", th_error),
                   "get error handle");
    std::vector<double> errors(fe_ents.size());
    CHK_MOAB_THROW(
        mField.get_moab().tag_get_data(th_error, fe_ents, &*errors.begin()),
        "get tag data");
    auto it = std::max_element(errors.begin(), errors.end());
    double max;
    MPI_Allreduce(&*it, &max, 1, MPI_DOUBLE, MPI_MAX, mField.get_comm());
    MOFEM_LOG("LevelSet", Sev::inform) << "Max error: " << max;
    auto threshold = wp.getThreshold(max);
 
    std::vector<EntityHandle> fe_to_refine;
    fe_to_refine.reserve(fe_ents.size());
 
    auto fe_it = fe_ents.begin();
    auto error_it = errors.begin();
    for (auto i = 0; i != fe_ents.size(); ++i) {
      if (*error_it > threshold) {
        fe_to_refine.push_back(*fe_it);
      }
      ++fe_it;
      ++error_it;
    }
 
    Range ents;
    ents.insert_list(fe_to_refine.begin(), fe_to_refine.end());
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        ents, nullptr, NOISY);
 
    auto get_neighbours_by_bridge_vertices = [&](auto &&ents) {
      Range verts;
      CHKERR mField.get_moab().get_connectivity(ents, verts, true);
      CHKERR mField.get_moab().get_adjacencies(verts, FE_DIM, false, ents,
                                               moab::Interface::UNION);
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ents);
      return ents;
    };
 
    ents = get_neighbours_by_bridge_vertices(ents);
 
#ifndef NDEBUG
    if (debug) {
      auto meshset_ptr = get_temp_meshset_ptr(mField.get_moab());
      CHK_MOAB_THROW(mField.get_moab().add_entities(*meshset_ptr, ents),
                     "add entities to meshset");
      CHKERR mField.get_moab().write_file(
          (proc_str + "_fe_to_refine.vtk").c_str(), "VTK", "",
          meshset_ptr->get_ptr(), 1);
    }
#endif
 
    return ents;
  };
 
  // refine elements, and set bit ref level
  auto refine_mesh = [&](auto l, auto &&fe_to_refine) {
    Skinner skin(&mField.get_moab());
    MoFEMFunctionBegin;
 
    // get entities in "l-1" level
    Range level_ents;
    CHKERR bit_mng->getEntitiesByDimAndRefLevel(
        set_bit(start_bit + l - 1), BitRefLevel().set(), FE_DIM, level_ents);
    // select entities to refine
    fe_to_refine = intersect(level_ents, fe_to_refine);
    // select entities not to refine
    level_ents = subtract(level_ents, fe_to_refine);
 
    // for entities to refine get children, i.e. redlined entities
    Range fe_to_refine_children;
    bit_mng->updateRangeByChildren(fe_to_refine, fe_to_refine_children);
    // add entities to to level "l"
    fe_to_refine_children = fe_to_refine_children.subset_by_dimension(FE_DIM);
    level_ents.merge(fe_to_refine_children);
 
    auto fix_neighbour_level = [&](auto ll) {
      MoFEMFunctionBegin;
      // filter entities on level ll
      auto level_ll = level_ents;
      CHKERR bit_mng->filterEntitiesByRefLevel(set_bit(ll), BitRefLevel().set(),
                                               level_ll);
      // find skin of ll level
      Range skin_edges;
      CHKERR skin.find_skin(0, level_ll, false, skin_edges);
      // get parents of skin of level ll
      Range skin_parents;
      for (auto lll = 0; lll <= ll; ++lll) {
        CHKERR bit_mng->updateRangeByParent(skin_edges, skin_parents);
      }
      // filter parents on level ll - 1
      BitRefLevel bad_bit;
      for (auto lll = 0; lll <= ll - 2; ++lll) {
        bad_bit[lll] = true;
      }
      // get adjacents to parents
      Range skin_adj_ents;
      CHKERR mField.get_moab().get_adjacencies(
          skin_parents, FE_DIM, false, skin_adj_ents, moab::Interface::UNION);
      CHKERR bit_mng->filterEntitiesByRefLevel(bad_bit, BitRefLevel().set(),
                                               skin_adj_ents);
      skin_adj_ents = intersect(skin_adj_ents, level_ents);
      if (!skin_adj_ents.empty()) {
        level_ents = subtract(level_ents, skin_adj_ents);
        Range skin_adj_ents_children;
        bit_mng->updateRangeByChildren(skin_adj_ents, skin_adj_ents_children);
        level_ents.merge(skin_adj_ents_children);
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR fix_neighbour_level(l);
 
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        level_ents);
 
    // get lower dimension entities for level "l"
    for (auto d = 0; d != FE_DIM; ++d) {
      if (d == 0) {
        CHKERR mField.get_moab().get_connectivity(
            level_ents.subset_by_dimension(FE_DIM), level_ents, true);
      } else {
        CHKERR mField.get_moab().get_adjacencies(
            level_ents.subset_by_dimension(FE_DIM), d, false, level_ents,
            moab::Interface::UNION);
      }
    }
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        level_ents);
 
    // set bit ref level to level entities
    CHKERR bit_mng->setNthBitRefLevel(start_bit + l, false);
    CHKERR bit_mng->setNthBitRefLevel(level_ents, start_bit + l, true);
 
#ifndef NDEBUG
    auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
    CHKERR bit_mng->writeBitLevelByDim(
        set_bit(start_bit + l), BitRefLevel().set(), FE_DIM,
        (boost::lexical_cast<std::string>(l) + "_" + proc_str + "_ref_mesh.vtk")
            .c_str(),
        "VTK", "");
#endif
 
    MoFEMFunctionReturn(0);
  };
 
  // set skeleton
  auto set_skelton_bit = [&](auto l) {
    MoFEMFunctionBegin;
 
    // get entities of dim-1 on level "l"
    Range level_edges;
    CHKERR bit_mng->getEntitiesByDimAndRefLevel(
        set_bit(start_bit + l), BitRefLevel().set(), FE_DIM - 1, level_edges);
 
    // get parent of entities of level "l"
    Range level_edges_parents;
    CHKERR bit_mng->updateRangeByParent(level_edges, level_edges_parents);
    level_edges_parents = level_edges_parents.subset_by_dimension(FE_DIM - 1);
    CHKERR bit_mng->filterEntitiesByRefLevel(
        set_bit(start_bit + l), BitRefLevel().set(), level_edges_parents);
 
    // skeleton entities which do not have parents
    auto parent_skeleton = intersect(level_edges, level_edges_parents);
    auto skeleton = subtract(level_edges, level_edges_parents);
 
    // add adjacent domain entities
    CHKERR mField.get_moab().get_adjacencies(unite(parent_skeleton, skeleton),
                                             FE_DIM, false, skeleton,
                                             moab::Interface::UNION);
 
    // set levels
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(skeleton);
    CHKERR bit_mng->setNthBitRefLevel(skeleton_bit, false);
    CHKERR bit_mng->setNthBitRefLevel(skeleton, skeleton_bit, true);
 
#ifndef NDEBUG
    CHKERR bit_mng->writeBitLevel(
        set_bit(skeleton_bit), BitRefLevel().set(),
        (boost::lexical_cast<std::string>(l) + "_" + proc_str + "_skeleton.vtk")
            .c_str(),
        "VTK", "");
#endif
    MoFEMFunctionReturn(0);
  };
 
  // Reset bit sand set old current and aggregate bits as projection bits
  Range level0_current;
  CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(current_bit),
                                        BitRefLevel().set(), level0_current);
 
  Range level0_aggregate;
  CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(aggregate_bit),
                                        BitRefLevel().set(), level0_aggregate);
 
  BitRefLevel start_mask;
  for (auto s = 0; s != start_bit; ++s)
    start_mask[s] = true;
  CHKERR bit_mng->lambdaBitRefLevel(
      [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
  CHKERR bit_mng->setNthBitRefLevel(level0_current, projection_bit, true);
  CHKERR bit_mng->setNthBitRefLevel(level0_aggregate, aggregate_projection_bit,
                                    true);
 
  // Set zero bit ref level
  Range level0;
  CHKERR bit_mng->getEntitiesByRefLevel(set_bit(0), BitRefLevel().set(),
                                        level0);
  CHKERR bit_mng->setNthBitRefLevel(level0, start_bit, true);
  CHKERR bit_mng->setNthBitRefLevel(level0, current_bit, true);
  CHKERR bit_mng->setNthBitRefLevel(level0, aggregate_bit, true);
  CHKERR bit_mng->setNthBitRefLevel(level0, skeleton_bit, true);
 
  CHKERR wp.setBits(*this, 0);
  CHKERR wp.runCalcs(*this, 0);
  for (auto l = 0; l != nb_levels; ++l) {
    MOFEM_LOG("WORLD", Sev::inform) << "Process level: " << l;
    CHKERR refine_mesh(l + 1, get_refined_elements_meshset(
                                  set_bit(start_bit + l), BitRefLevel().set()));
    CHKERR set_skelton_bit(l + 1);
    CHKERR wp.setAggregateBit(*this, l + 1);
    CHKERR wp.setBits(*this, l + 1);
    CHKERR wp.runCalcs(*this, l + 1);
  }
 
  MoFEMFunctionReturn(0);
}

runProblem()

MoFEMErrorCode LevelSet::runProblem

(

)

Examples

level_set.cpp.

Definition at line 386 of file level_set.cpp.

                                    {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setUpProblem();
 
  if constexpr (debug) {
    CHKERR testSideFE();
    CHKERR testOp();
  }
 
  CHKERR initialiseFieldVelocity();
 
  maxPtr = boost::make_shared<double>(0);
  CHKERR refineMesh(WrapperClassInitalSolution(maxPtr));
 
  auto simple = mField.getInterface<Simple>();
  simple->getBitRefLevel() =
      BitRefLevel().set(skeleton_bit) | BitRefLevel().set(aggregate_bit);
  simple->getBitRefLevelMask() = BitRefLevel().set();
  simple->reSetUp(true);
 
  CHKERR solveAdvection();
 
  MoFEMFunctionReturn(0);
}

setUpProblem()

MoFEMErrorCode LevelSet::setUpProblem ( )

private

create fields, and set approximation order

Returns

MoFEMErrorCode

Examples

level_set.cpp.

Definition at line 560 of file level_set.cpp.

                                      {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  // Scalar fields and vector field is tested. Add more fields, i.e. vector
  // field if needed.
  CHKERR simple->addDomainField("L", L2, AINSWORTH_LEGENDRE_BASE, 1);
  CHKERR simple->addDomainField("V", potential_velocity_space,
                                AINSWORTH_LEGENDRE_BASE, 1);
 
  // set fields order, i.e. for most first cases order is sufficient.
  CHKERR simple->setFieldOrder("L", 4);
  CHKERR simple->setFieldOrder("V", 4);
 
  // setup problem
  CHKERR simple->setUp();
 
  MoFEMFunctionReturn(0);
}

solveAdvection()

MoFEMErrorCode LevelSet::solveAdvection ( )

private

solve advection problem

Returns

* MoFEMErrorCode

Examples

level_set.cpp.

Definition at line 1934 of file level_set.cpp.

                                        {
  MoFEMFunctionBegin;
 
  // get operators tester
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>(); // get interface to
  auto prb_mng = mField.getInterface<ProblemsManager>();
 
  CHKERR pushOpDomain();
  CHKERR pushOpSkeleton();
 
  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "ADVECTION");
  CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
  CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
  CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
  CHKERR DMMoFEMAddElement(sub_dm, simple->getSkeletonFEName());
  CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
  CHKERR DMSetUp(sub_dm);
 
  BitRefLevel remove_mask = BitRefLevel().set(current_bit);
  remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
  CHKERR prb_mng->removeDofsOnEntities("ADVECTION", "L", BitRefLevel().set(),
                                       remove_mask);
 
  auto add_post_proc_fe = [&]() {
    auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
 
    Tag th_error;
    double def_val = 0;
    CHKERR mField.get_moab().tag_get_handle(
        "Error", 1, MB_TYPE_DOUBLE, th_error, MB_TAG_CREAT | MB_TAG_SPARSE,
        &def_val);
    post_proc_fe->setTagsToTransfer({th_error});
 
    post_proc_fe->exeTestHook = [&](FEMethod *fe_ptr) {
      return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
          current_bit);
    };
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    auto vel_ptr = boost::make_shared<MatrixDouble>();
 
    CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
        post_proc_fe->getOpPtrVector(), {L2});
    post_proc_fe->getOpPtrVector().push_back(getZeroLevelVelOp(vel_ptr));
 
    if constexpr (DIM1 == 1 && DIM2 == 1) {
      auto l_vec = boost::make_shared<VectorDouble>();
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateScalarFieldValues("L", l_vec));
      post_proc_fe->getOpPtrVector().push_back(
 
          new OpPPMap(
 
              post_proc_fe->getPostProcMesh(),
 
              post_proc_fe->getMapGaussPts(),
 
              {{"L", l_vec}},
 
              {{"V", vel_ptr}},
 
              {}, {})
 
      );
    }
    return post_proc_fe;
  };
 
  auto post_proc_fe = add_post_proc_fe();
 
  auto set_time_monitor = [&](auto dm, auto ts) {
    auto monitor_ptr = boost::make_shared<FEMethod>();
 
    monitor_ptr->preProcessHook = []() { return 0; };
    monitor_ptr->operatorHook = []() { return 0; };
    monitor_ptr->postProcessHook = [&]() {
      MoFEMFunctionBegin;
 
      if (!post_proc_fe)
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "Null pointer for post proc element");
 
      CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                      post_proc_fe);
      CHKERR post_proc_fe->writeFile(
          "level_set_" +
          boost::lexical_cast<std::string>(monitor_ptr->ts_step) + ".h5m");
      MoFEMFunctionReturn(0);
    };
 
    boost::shared_ptr<FEMethod> null;
    DMMoFEMTSSetMonitor(sub_dm, ts, simple->getDomainFEName(), monitor_ptr,
                        null, null);
 
    return monitor_ptr;
  };
 
  auto ts = pip->createTSIM(sub_dm);
 
  auto set_solution = [&](auto ts) {
    MoFEMFunctionBegin;
    auto D = createDMVector(sub_dm);
    CHKERR DMoFEMMeshToLocalVector(sub_dm, D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR TSSetSolution(ts, D);
    MoFEMFunctionReturn(0);
  };
  CHKERR set_solution(ts);
 
  auto monitor_pt = set_time_monitor(sub_dm, ts);
  CHKERR TSSetFromOptions(ts);
 
  auto B = createDMMatrix(sub_dm);
  CHKERR TSSetIJacobian(ts, B, B, TsSetIJacobian, nullptr);
  level_set_raw_ptr = this;
 
  CHKERR TSSetUp(ts);
 
  auto ts_pre_step = [](TS ts) {
    auto &m_field = level_set_raw_ptr->mField;
    auto simple = m_field.getInterface<Simple>();
    auto bit_mng = m_field.getInterface<BitRefManager>();
    MoFEMFunctionBegin;
 
    auto [error, th_error] = level_set_raw_ptr->evaluateError();
    MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
 
    auto get_norm = [&](auto x) {
      double nrm;
      CHKERR VecNorm(x, NORM_2, &nrm);
      return nrm;
    };
 
    auto set_solution = [&](auto ts) {
      MoFEMFunctionBegin;
      DM dm;
      CHKERR TSGetDM(ts, &dm);
      auto prb_ptr = getProblemPtr(dm);
 
      auto x = createDMVector(dm);
      CHKERR DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
 
      MOFEM_LOG("LevelSet", Sev::inform)
          << "Problem " << prb_ptr->getName() << " solution vector norm "
          << get_norm(x);
      CHKERR TSSetSolution(ts, x);
 
      MoFEMFunctionReturn(0);
    };
 
    auto refine_and_project = [&](auto ts) {
      MoFEMFunctionBegin;
 
      CHKERR level_set_raw_ptr->refineMesh(
          WrapperClassErrorProjection(level_set_raw_ptr->maxPtr));
      simple->getBitRefLevel() = BitRefLevel().set(skeleton_bit) |
                                 BitRefLevel().set(aggregate_bit) |
                                 BitRefLevel().set(aggregate_projection_bit);
 
      simple->reSetUp(true);
      DM dm;
      CHKERR TSGetDM(ts, &dm);
      CHKERR DMSubDMSetUp_MoFEM(dm);
 
      BitRefLevel remove_mask = BitRefLevel().set(current_bit);
      remove_mask
          .flip(); // DOFs which are not on bit_domain_ele should be removed
      CHKERR level_set_raw_ptr->mField.getInterface<ProblemsManager>()
          ->removeDofsOnEntities("ADVECTION", "L", BitRefLevel().set(),
                                 remove_mask);
 
      MoFEMFunctionReturn(0);
    };
 
    auto ts_reset_theta = [&](auto ts) {
      MoFEMFunctionBegin;
      DM dm;
      CHKERR TSGetDM(ts, &dm);
 
      // FIXME: Look into vec-5 how to transfer internal theta method variables
 
      CHKERR TSReset(ts);
      CHKERR TSSetUp(ts);
 
      CHKERR level_set_raw_ptr->dgProjection(projection_bit);
      CHKERR set_solution(ts);
 
      auto B = createDMMatrix(dm);
      CHKERR TSSetIJacobian(ts, B, B, TsSetIJacobian, nullptr);
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR refine_and_project(ts);
    CHKERR ts_reset_theta(ts);
 
    MoFEMFunctionReturn(0);
  };
 
  auto ts_post_step = [](TS ts) { return 0; };
 
  CHKERR TSSetPreStep(ts, ts_pre_step);
  CHKERR TSSetPostStep(ts, ts_post_step);
 
  CHKERR TSSolve(ts, NULL);
 
  MoFEMFunctionReturn(0);
}

testOp()

MoFEMErrorCode LevelSet::testOp ( )

private

test consistency between tangent matrix and the right hand side vectors

Returns

MoFEMErrorCode

Examples

level_set.cpp.

Definition at line 1602 of file level_set.cpp.

                                {
  MoFEMFunctionBegin;
 
  // get operators tester
  auto simple = mField.getInterface<Simple>();
  auto opt = mField.getInterface<OperatorsTester>(); // get interface to
                                                     // OperatorsTester
  auto pip = mField.getInterface<PipelineManager>(); // get interface to
                                                     // pipeline manager
 
  CHKERR pushOpDomain();
  CHKERR pushOpSkeleton();
 
  auto post_proc = [&](auto dm, auto f_res, auto out_name) {
    MoFEMFunctionBegin;
    auto post_proc_fe =
        boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
 
    if constexpr (DIM1 == 1 && DIM2 == 1) {
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      auto l_vec = boost::make_shared<VectorDouble>();
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateScalarFieldValues("L", l_vec, f_res));
 
      post_proc_fe->getOpPtrVector().push_back(
 
          new OpPPMap(
 
              post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
              {{"L", l_vec}},
 
              {},
 
              {}, {})
 
      );
    }
 
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                    post_proc_fe);
    post_proc_fe->writeFile(out_name);
    MoFEMFunctionReturn(0);
  };
 
  constexpr double eps = 1e-4;
 
  auto x =
      opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}, {"V", {-1, 1}}});
  auto dot_x = opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}});
  auto diff_x = opt->setRandomFields(simple->getDM(), {{"L", {-1, 1}}});
 
  auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline,
                             auto rhs_pipeline) {
    MoFEMFunctionBegin;
 
    auto diff_res = opt->checkCentralFiniteDifference(
        simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
        SmartPetscObj<Vec>(), diff_x, 0, 1, eps);
 
    if constexpr (debug) {
      // Example how to plot direction in direction diff_x. If instead
      // directionalCentralFiniteDifference(...) diff_res is used, then error
      // on directive is plotted.
      CHKERR post_proc(simple->getDM(), diff_res, "tangent_op_error.h5m");
    }
 
    // Calculate norm of difference between directive calculated from finite
    // difference, and tangent matrix.
    double fnorm;
    CHKERR VecNorm(diff_res, NORM_2, &fnorm);
    MOFEM_LOG_C("LevelSet", Sev::inform,
                "Test consistency of tangent matrix %3.4e", fnorm);
 
    constexpr double err = 1e-9;
    if (fnorm > err)
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
               "Norm of directional derivative too large err = %3.4e", fnorm);
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
                         pip->getDomainRhsFE());
  CHKERR test_domain_ops(simple->getSkeletonFEName(), pip->getSkeletonLhsFE(),
                         pip->getSkeletonRhsFE());
 
  MoFEMFunctionReturn(0);
};

testSideFE()

MoFEMErrorCode LevelSet::testSideFE ( )

private

test integration side elements

test side element

Check consistency between volume and skeleton integral.

Returns

MoFEMErrorCode

Check consistency between volume and skeleton integral

Returns

MoFEMErrorCode

calculate volume

calculate skeleton integral

Set up problem such that gradient of level set field is orthogonal to velocity field. Then volume and skeleton integral should yield the same value.

Examples

level_set.cpp.

Definition at line 1387 of file level_set.cpp.

                                    {
  MoFEMFunctionBegin;
 
  /**
   * @brief calculate volume
   *
   */
  struct DivergenceVol : public DomainEleOp {
    DivergenceVol(boost::shared_ptr<MatrixDouble> l_ptr,
                  boost::shared_ptr<MatrixDouble> vel_ptr,
                  SmartPetscObj<Vec> div_vec)
        : DomainEleOp("L", DomainEleOp::OPROW), lPtr(l_ptr), velPtr(vel_ptr),
          divVec(div_vec) {}
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
      const auto nb_dofs = data.getIndices().size();
      if (nb_dofs) {
        const auto nb_gauss_pts = getGaussPts().size2();
        const auto t_w = getFTensor0IntegrationWeight();
        auto t_diff = data.getFTensor1DiffN<SPACE_DIM>();
        auto t_l = getFTensor2FromMat<DIM1, DIM2>(*lPtr);
        auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
        double div = 0;
        for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
          for (int rr = 0; rr != nb_dofs; ++rr) {
            div += getMeasure() * t_w * t_l(I, I) * (t_diff(i) * t_vel(i));
            ++t_diff;
          }
          ++t_w;
          ++t_l;
          ++t_vel;
        }
        CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
      }
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<MatrixDouble> lPtr;
    boost::shared_ptr<MatrixDouble> velPtr;
    SmartPetscObj<Vec> divVec;
  };
 
  /**
   * @brief calculate skeleton integral
   *
   */
  struct DivergenceSkeleton : public BoundaryEleOp {
    DivergenceSkeleton(boost::shared_ptr<SideData> side_data_ptr,
                       boost::shared_ptr<FaceSideEle> side_fe_ptr,
                       SmartPetscObj<Vec> div_vec)
        : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
          sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr), divVec(div_vec) {}
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
 
      auto get_ntensor = [](auto &base_mat) {
        return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(
            &*base_mat.data().begin());
      };
 
      auto not_side = [](auto s) {
        return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
      };
 
      // Collect data from side domain elements
      CHKERR loopSideFaces("dFE", *sideFEPtr);
      const auto in_the_loop =
          sideFEPtr->nInTheLoop; // return number of elements on the side
 
      auto t_normal = getFTensor1Normal();
      const auto nb_gauss_pts = getGaussPts().size2();
      for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
        const auto nb_dofs = sideDataPtr->indicesRowSideMap[s0].size();
        if (nb_dofs) {
          auto t_base = get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
          auto nb_row_base_functions = sideDataPtr->rowBaseSideMap[s0].size2();
          auto side_sense = sideDataPtr->senseMap[s0];
          auto opposite_s0 = not_side(s0);
 
          auto arr_t_l = make_array(
              getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[LEFT_SIDE]),
              getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[RIGHT_SIDE]));
          auto arr_t_vel = make_array(
              getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[LEFT_SIDE]),
              getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[RIGHT_SIDE]));
 
          auto next = [&]() {
            for (auto &t_l : arr_t_l)
              ++t_l;
            for (auto &t_vel : arr_t_vel)
              ++t_vel;
          };
 
          double div = 0;
 
          auto t_w = getFTensor0IntegrationWeight();
          for (int gg = 0; gg != nb_gauss_pts; ++gg) {
            FTensor::Tensor1<double, SPACE_DIM> t_vel;
            t_vel(i) =
                (arr_t_vel[LEFT_SIDE](i) + arr_t_vel[RIGHT_SIDE](i)) / 2.;
            const auto dot = (t_normal(i) * t_vel(i)) * side_sense;
            const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
            const auto l_upwind =
                arr_t_l[l_upwind_side]; //< assume that field is continues,
                                        // initialisation field has to be smooth
                                        // and exactly approximated by approx
                                        // base
            auto res = t_w * l_upwind(I, I) * dot;
            ++t_w;
            next();
            int rr = 0;
            for (; rr != nb_dofs; ++rr) {
              div += t_base * res;
              ++t_base;
            }
            for (; rr < nb_row_base_functions; ++rr) {
              ++t_base;
            }
          }
          CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
        }
        if (!in_the_loop)
          break;
      }
 
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<SideData> sideDataPtr;
    boost::shared_ptr<FaceSideEle> sideFEPtr;
    boost::shared_ptr<MatrixDouble> velPtr;
    SmartPetscObj<Vec> divVec;
  };
 
  auto vol_fe = boost::make_shared<DomainEle>(mField);
  auto skel_fe = boost::make_shared<BoundaryEle>(mField);
 
  vol_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
  skel_fe->getRuleHook = [](int, int, int o) { return 3 * o; };
 
  auto div_vol_vec = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
  auto div_skel_vec = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
 
  auto l_ptr = boost::make_shared<MatrixDouble>();
  auto vel_ptr = boost::make_shared<MatrixDouble>();
  auto side_data_ptr = boost::make_shared<SideData>();
  auto side_fe_ptr = getSideFE(side_data_ptr);
 
  CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(vol_fe->getOpPtrVector(),
                                                  {L2});
  vol_fe->getOpPtrVector().push_back(
      new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_ptr));
  vol_fe->getOpPtrVector().push_back(getZeroLevelVelOp(vel_ptr));
  vol_fe->getOpPtrVector().push_back(
      new DivergenceVol(l_ptr, vel_ptr, div_vol_vec));
 
  skel_fe->getOpPtrVector().push_back(
      new DivergenceSkeleton(side_data_ptr, side_fe_ptr, div_skel_vec));
 
  auto simple = mField.getInterface<Simple>();
  auto dm = simple->getDM();
 
  /**
   * 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.
   */
 
  CHKERR initialiseFieldVelocity(
      [](double x, double y, double) { return x - y; });
  CHKERR initialiseFieldLevelSet(
      [](double x, double y, double) { return x - y; });
 
  vol_fe->exeTestHook = [&](FEMethod *fe_ptr) {
    return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
        current_bit);
  };
  skel_fe->exeTestHook = [&](FEMethod *fe_ptr) {
    return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
        skeleton_bit);
  };
 
  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), vol_fe);
  CHKERR DMoFEMLoopFiniteElements(dm, simple->getSkeletonFEName(), skel_fe);
  CHKERR DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(), skel_fe);
 
  auto assemble_and_sum = [](auto vec) {
    CHK_THROW_MESSAGE(VecAssemblyBegin(vec), "assemble");
    CHK_THROW_MESSAGE(VecAssemblyEnd(vec), "assemble");
    double sum;
    CHK_THROW_MESSAGE(VecSum(vec, &sum), "assemble");
    return sum;
  };
 
  auto div_vol = assemble_and_sum(div_vol_vec);
  auto div_skel = assemble_and_sum(div_skel_vec);
 
  auto eps = std::abs((div_vol - div_skel) / (div_vol + div_skel));
 
  MOFEM_LOG("WORLD", Sev::inform) << "Testing divergence volume: " << div_vol;
  MOFEM_LOG("WORLD", Sev::inform) << "Testing divergence skeleton: " << div_skel
                                  << " relative difference: " << eps;
 
  constexpr double eps_err = 1e-6;
  if (eps > eps_err)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "No consistency between skeleton integral and volume integral");
 
  MoFEMFunctionReturn(0);
};

Member Data Documentation

maxPtr

boost::shared_ptr<double> LevelSet::maxPtr

private

Examples

level_set.cpp.

Definition at line 370 of file level_set.cpp.

mField

MoFEM::Interface& LevelSet::mField

private

integrate skeleton operators on khs

Examples

level_set.cpp.

Definition at line 349 of file level_set.cpp.

---

The documentation for this struct was generated from the following file:

• tutorials/adv-3/level_set.cpp