FreeSurface Struct Reference

Collaboration diagram for FreeSurface:

Public Member Functions
	FreeSurface (MoFEM::Interface &m_field)
	Constructor.
MoFEMErrorCode	runProblem ()
	Main function to run the complete free surface simulation.
MoFEMErrorCode	makeRefProblem ()
	Create refined problem for mesh adaptation.

Public Attributes
MoFEM::Interface &	mField

Private Member Functions
MoFEMErrorCode	readMesh ()
	Read mesh from input file.
MoFEMErrorCode	setupProblem ()
	Setup problem fields and parameters.
MoFEMErrorCode	boundaryCondition ()
	Apply boundary conditions and initialize fields.
MoFEMErrorCode	projectData ()
	Project solution data between mesh levels.
MoFEMErrorCode	assembleSystem ()
	Assemble system operators and matrices.
MoFEMErrorCode	solveSystem ()
	Solve the time-dependent free surface problem.
std::vector< Range >	findEntitiesCrossedByPhaseInterface (size_t overlap)
	Find entities on refinement levels.
Range	findParentsToRefine (Range ents, BitRefLevel level, BitRefLevel mask)
	Find parent entities that need refinement.
MoFEMErrorCode	refineMesh (size_t overlap)
	Perform adaptive mesh refinement.
MoFEMErrorCode	setParentDofs (boost::shared_ptr< FEMethod > fe_top, std::string field_name, ForcesAndSourcesCore::UserDataOperator::OpType op, BitRefLevel child_ent_bit, boost::function< boost::shared_ptr< ForcesAndSourcesCore >()> get_elem, int verbosity, LogManager::SeverityLevel sev)
	Create hierarchy of elements run on parent levels.
MoFEMErrorCode	checkResults ()
	Check results for correctness.

Friends
struct	TSPrePostProc

Detailed Description

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 797 of file free_surface.cpp.

Constructor & Destructor Documentation

FreeSurface()

FreeSurface::FreeSurface ( MoFEM::Interface &  m_field )

inline

Constructor.

Parameters

m_field

Reference to MoFEM interface for finite element operations

Definition at line 804 of file free_surface.cpp.

: mField(m_field) {}

Member Function Documentation

assembleSystem()

MoFEMErrorCode FreeSurface::assembleSystem ( )

private

Assemble system operators and matrices.

[Data projection]

Returns

MoFEMErrorCode Success/failure code

Sets up finite element operators for:

• Domain integration (momentum, phase field, mass conservation)
• Boundary integration (surface tension, wetting angle, Lagrange multipliers)
• Parent-child mesh hierarchies for adaptive refinement

[Push operators to pip]

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 2030 of file free_surface.cpp.

                                           {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  using UDO = ForcesAndSourcesCore::UserDataOperator;
 
  auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
    return setParentDofs(
        fe, field, op, bit(get_skin_parent_bit()),
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new DomainParentEle(mField));
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
  };
 
  auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
    return setParentDofs(
        fe, field, op, bit(get_skin_parent_bit()),
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new BoundaryParentEle(mField));
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
  };
 
  auto test_bit_child = [](FEMethod *fe_ptr) {
    return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
        get_start_bit() + nb_levels - 1);
  };
 
  auto dot_u_ptr = boost::make_shared<MatrixDouble>(); // time derivative of velocity u ie du/dt
  auto u_ptr = boost::make_shared<MatrixDouble>(); // velocity u
  auto grad_u_ptr = boost::make_shared<MatrixDouble>(); // velocity gradient tensor ie grad(u)
  auto dot_h_ptr = boost::make_shared<VectorDouble>(); // time derivative of phase field ie dh/dt
  auto h_ptr = boost::make_shared<VectorDouble>(); // phase field variable h
  auto grad_h_ptr = boost::make_shared<MatrixDouble>(); // phase field gradient ie grad(h)
  auto g_ptr = boost::make_shared<VectorDouble>(); // chemical potential g
  auto grad_g_ptr = boost::make_shared<MatrixDouble>(); // chemical potential gradient ie grad(g)
  auto lambda_ptr = boost::make_shared<VectorDouble>(); // Lagrange multiplier lambda
  auto p_ptr = boost::make_shared<VectorDouble>(); // pressure p
  auto div_u_ptr = boost::make_shared<VectorDouble>(); // divergence of velocity ie div(u)
 
  // Push element from reference configuration to current configuration in 3d
  // space
  auto set_domain_general = [&](auto fe) {
    MoFEMFunctionBegin;
    auto &pip = fe->getOpPtrVector();
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
 
    CHKERR setParentDofs(
        fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new DomainParentEle(mField));
          AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
              fe_parent->getOpPtrVector(), {H1});
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
 
    CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
    CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
    pip.push_back(
        new OpCalculateVectorFieldValuesDot<U_FIELD_DIM>("U", dot_u_ptr));
    pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
    pip.push_back(new OpCalculateVectorFieldGradient<U_FIELD_DIM, SPACE_DIM>(
        "U", grad_u_ptr));
 
    switch (coord_type) {
    case CARTESIAN:
      pip.push_back(
          new OpCalculateDivergenceVectorFieldValues<SPACE_DIM, CARTESIAN>(
              "U", div_u_ptr));
      break;
    case CYLINDRICAL:
      pip.push_back(
          new OpCalculateDivergenceVectorFieldValues<SPACE_DIM, CYLINDRICAL>(
              "U", div_u_ptr));
      break;
    default:
      SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
              "Coordinate system not implemented");
    }
 
    CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
    pip.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
    pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
    pip.push_back(
        new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
 
    CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
    pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
    pip.push_back(
        new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
 
    CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P");
    pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
    MoFEMFunctionReturn(0);
  };
 
  auto set_domain_rhs = [&](auto fe) {
    MoFEMFunctionBegin;
    auto &pip = fe->getOpPtrVector();
 
    CHKERR set_domain_general(fe);
 
    CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
    pip.push_back(new OpRhsU("U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
                             grad_h_ptr, g_ptr, p_ptr));
 
    CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
    pip.push_back(new OpRhsH<false>("H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
                                    grad_g_ptr));
 
    CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
    pip.push_back(new OpRhsG<false>("G", h_ptr, grad_h_ptr, g_ptr));
 
    CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
    pip.push_back(new OpDomainAssembleScalar(
        "P", div_u_ptr, [](const double r, const double, const double) {
          return cylindrical(r);
        }));
    pip.push_back(new OpDomainAssembleScalar(
        "P", p_ptr, [](const double r, const double, const double) {
          return eps * cylindrical(r);
        }));
    MoFEMFunctionReturn(0);
  };
 
  auto set_domain_lhs = [&](auto fe) {
    MoFEMFunctionBegin;
    auto &pip = fe->getOpPtrVector();
 
    CHKERR set_domain_general(fe);
 
    CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
    {
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
      pip.push_back(new OpLhsU_dU("U", u_ptr, grad_u_ptr, h_ptr));
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
      pip.push_back(
          new OpLhsU_dH("U", "H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
 
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
      pip.push_back(new OpLhsU_dG("U", "G", grad_h_ptr));
    }
 
    CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
    {
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
      pip.push_back(new OpLhsH_dU("H", "U", grad_h_ptr));
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
      pip.push_back(new OpLhsH_dH<false>("H", u_ptr, h_ptr, grad_g_ptr));
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
      pip.push_back(new OpLhsH_dG<false>("H", "G", h_ptr));
    }
 
    CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
    {
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
      pip.push_back(new OpLhsG_dH<false>("G", "H", h_ptr));
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
      pip.push_back(new OpLhsG_dG("G"));
    }
 
    CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
    {
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
 
      switch (coord_type) {
      case CARTESIAN:
        pip.push_back(new OpMixScalarTimesDiv<CARTESIAN>(
            "P", "U",
            [](const double r, const double, const double) {
              return cylindrical(r);
            },
            true, false));
        break;
      case CYLINDRICAL:
        pip.push_back(new OpMixScalarTimesDiv<CYLINDRICAL>(
            "P", "U",
            [](const double r, const double, const double) {
              return cylindrical(r);
            },
            true, false));
        break;
      default:
        SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
                "Coordinate system not implemented");
      }
 
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P");
      pip.push_back(new OpDomainMassP("P", "P", [](double r, double, double) {
        return eps * cylindrical(r);
      }));
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto get_block_name = [](auto name) {
    return boost::format("%s(.*)") % "WETTING_ANGLE";
  };
 
  auto get_blocks = [&](auto &&name) {
    return mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
        std::regex(name.str()));
  };
 
  auto set_boundary_rhs = [&](auto fe) {
    MoFEMFunctionBegin;
    auto &pip = fe->getOpPtrVector();
 
    CHKERR setParentDofs(
        fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new BoundaryParentEle(mField));
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
 
    CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "U");
    pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
 
    CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "L");
    pip.push_back(new OpCalculateScalarFieldValues("L", lambda_ptr));
    pip.push_back(new OpNormalConstrainRhs("L", u_ptr));
 
    CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpFluidFlux>::add(
        pip, mField, "L", {}, "INFLUX",
        [](double r, double, double) { return cylindrical(r); }, Sev::inform);
 
    CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
    pip.push_back(new OpNormalForceRhs("U", lambda_ptr));
 
    auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
    if (wetting_block.size()) {
      // push operators to the side element which is called from op_bdy_side
      auto op_bdy_side =
          new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
      op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
 
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_bdy_side->getOpPtrVector(), {H1});
 
      CHKERR setParentDofs(
          op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
          bit(get_skin_parent_bit()),
 
          [&]() {
            boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
                new DomainParentEle(mField));
            AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
                fe_parent->getOpPtrVector(), {H1});
            return fe_parent;
          },
 
          QUIET, Sev::noisy);
 
      CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
                                     "H");
      op_bdy_side->getOpPtrVector().push_back(
          new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
      // push bdy side op
      pip.push_back(op_bdy_side);
 
      // push operators for rhs wetting angle
      for (auto &b : wetting_block) {
        Range force_edges;
        std::vector<double> attr_vec;
        CHKERR b->getMeshsetIdEntitiesByDimension(
            mField.get_moab(), SPACE_DIM - 1, force_edges, true);
        b->getAttributes(attr_vec);
        if (attr_vec.size() != 1)
          SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                  "Should be one attribute");
        MOFEM_LOG("FS", Sev::inform) << "Wetting angle: " << attr_vec.front();
        // need to find the attributes and pass to operator
        CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
        CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "G");
        pip.push_back(new OpWettingAngleRhs(
            "G", grad_h_ptr, boost::make_shared<Range>(force_edges),
            attr_vec.front()));
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto set_boundary_lhs = [&](auto fe) {
    MoFEMFunctionBegin;
    auto &pip = fe->getOpPtrVector();
 
    CHKERR setParentDofs(
        fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new BoundaryParentEle(mField));
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
 
    CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
    CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "U");
    pip.push_back(new OpNormalConstrainLhs("L", "U"));
 
    auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
    if (wetting_block.size()) {
      auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
      auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
 
      // push operators to the side element which is called from op_bdy_side
      auto op_bdy_side =
          new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
      op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
 
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_bdy_side->getOpPtrVector(), {H1});
 
      CHKERR setParentDofs(
          op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
          bit(get_skin_parent_bit()),
 
          [&]() {
            boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
                new DomainParentEle(mField));
            AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
                fe_parent->getOpPtrVector(), {H1});
            return fe_parent;
          },
 
          QUIET, Sev::noisy);
 
      CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
                                     "H");
      CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
                                     "H");
      op_bdy_side->getOpPtrVector().push_back(
          new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
      op_bdy_side->getOpPtrVector().push_back(
          new OpLoopSideGetDataForSideEle("H", col_ind_ptr, col_diff_base_ptr));
 
      // push bdy side op
      pip.push_back(op_bdy_side);
 
      // push operators for lhs wetting angle
      for (auto &b : wetting_block) {
        Range force_edges;
        std::vector<double> attr_vec;
        CHKERR b->getMeshsetIdEntitiesByDimension(
            mField.get_moab(), SPACE_DIM - 1, force_edges, true);
        b->getAttributes(attr_vec);
        if (attr_vec.size() != 1)
          SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                  "Should be one attribute");
        MOFEM_LOG("FS", Sev::inform)
            << "wetting angle edges size " << force_edges.size();
 
        CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
        CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "G");
        pip.push_back(new OpWettingAngleLhs(
            "G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
            boost::make_shared<Range>(force_edges), attr_vec.front()));
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto *pip_mng = mField.getInterface<PipelineManager>();
 
  CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
  CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
  CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
  CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
 
  CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
  CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
 
  pip_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
  pip_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
  pip_mng->getBoundaryLhsFE()->exeTestHook = test_bit_child;
  pip_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
 
  MoFEMFunctionReturn(0);
}

boundaryCondition()

MoFEMErrorCode FreeSurface::boundaryCondition ( )

private

Apply boundary conditions and initialize fields.

[Set up problem]

Returns

MoFEMErrorCode Success/failure code

• Initializes phase field using analytical or Python functions
• Solves initialization problem for consistent initial conditions
• Performs mesh refinement based on interface location
• Removes DOFs for boundary conditions (symmetry, fixed, etc.)

[Boundary condition]

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 1109 of file free_surface.cpp.

                                              {
  MoFEMFunctionBegin;
 
#ifdef PYTHON_INIT_SURFACE
  auto get_py_surface_init = []() {
    auto py_surf_init = boost::make_shared<SurfacePython>();
    CHKERR py_surf_init->surfaceInit("surface.py");
    surfacePythonWeakPtr = py_surf_init;
    return py_surf_init;
  };
  auto py_surf_init = get_py_surface_init();
#endif
 
  auto simple = mField.getInterface<Simple>();
  auto pip_mng = mField.getInterface<PipelineManager>();
  auto bc_mng = mField.getInterface<BcManager>();
  auto bit_mng = mField.getInterface<BitRefManager>();
  auto dm = simple->getDM();
 
  using UDO = ForcesAndSourcesCore::UserDataOperator;
 
  auto reset_bits = [&]() {
    MoFEMFunctionBegin;
    BitRefLevel start_mask;
    for (auto s = 0; s != get_start_bit(); ++s)
      start_mask[s] = true;
    // reset bit ref levels
    CHKERR bit_mng->lambdaBitRefLevel(
        [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
    Range level0;
    CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
    CHKERR bit_mng->setNthBitRefLevel(level0, get_current_bit(), true);
    CHKERR bit_mng->setNthBitRefLevel(level0, get_projection_bit(), true);
    MoFEMFunctionReturn(0);
  };
 
  auto add_parent_field = [&](auto fe, auto op, auto field) {
    return setParentDofs(
        fe, field, op, bit(get_skin_parent_bit()),
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new DomainParentEle(mField));
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
  };
 
  auto h_ptr = boost::make_shared<VectorDouble>();
  auto grad_h_ptr = boost::make_shared<MatrixDouble>();
  auto g_ptr = boost::make_shared<VectorDouble>();
  auto grad_g_ptr = boost::make_shared<MatrixDouble>();
 
  auto set_generic = [&](auto fe) {
    MoFEMFunctionBegin;
    auto &pip = fe->getOpPtrVector();
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
 
    CHKERR setParentDofs(
        fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new DomainParentEle(mField));
          AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
              fe_parent->getOpPtrVector(), {H1});
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
 
    CHKERR add_parent_field(fe, UDO::OPCOL, "H");
    pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
    pip.push_back(
        new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
 
    CHKERR add_parent_field(fe, UDO::OPCOL, "G");
    pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
    pip.push_back(
        new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
 
    MoFEMFunctionReturn(0);
  };
 
  auto post_proc = [&](auto exe_test) {
    MoFEMFunctionBegin;
    auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
    post_proc_fe->exeTestHook = exe_test;
 
    CHKERR set_generic(post_proc_fe);
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(post_proc_fe->getPostProcMesh(),
                    post_proc_fe->getMapGaussPts(),
 
                    {{"H", h_ptr}, {"G", g_ptr}},
 
                    {{"GRAD_H", grad_h_ptr}, {"GRAD_G", grad_g_ptr}},
 
                    {},
 
                    {}
 
                    )
 
    );
 
    CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
    CHKERR post_proc_fe->writeFile("out_init.h5m");
 
    MoFEMFunctionReturn(0);
  };
 
  auto solve_init = [&](auto exe_test) {
    MoFEMFunctionBegin;
 
    pip_mng->getOpDomainRhsPipeline().clear();
    pip_mng->getOpDomainLhsPipeline().clear();
 
    auto set_domain_rhs = [&](auto fe) {
      MoFEMFunctionBegin;
      CHKERR set_generic(fe);
      auto &pip = fe->getOpPtrVector();
 
      CHKERR add_parent_field(fe, UDO::OPROW, "H");
      pip.push_back(new OpRhsH<true>("H", nullptr, nullptr, h_ptr, grad_h_ptr,
                                     grad_g_ptr));
      CHKERR add_parent_field(fe, UDO::OPROW, "G");
      pip.push_back(new OpRhsG<true>("G", h_ptr, grad_h_ptr, g_ptr));
      MoFEMFunctionReturn(0);
    };
 
    auto set_domain_lhs = [&](auto fe) {
      MoFEMFunctionBegin;
 
      CHKERR set_generic(fe);
      auto &pip = fe->getOpPtrVector();
 
      CHKERR add_parent_field(fe, UDO::OPROW, "H");
      CHKERR add_parent_field(fe, UDO::OPCOL, "H");
      pip.push_back(new OpLhsH_dH<true>("H", nullptr, h_ptr, grad_g_ptr));
 
      CHKERR add_parent_field(fe, UDO::OPCOL, "G");
      pip.push_back(new OpLhsH_dG<true>("H", "G", h_ptr));
 
      CHKERR add_parent_field(fe, UDO::OPROW, "G");
      pip.push_back(new OpLhsG_dG("G"));
 
      CHKERR add_parent_field(fe, UDO::OPCOL, "H");
      pip.push_back(new OpLhsG_dH<true>("G", "H", h_ptr));
      MoFEMFunctionReturn(0);
    };
 
    auto create_subdm = [&]() {
      auto level_ents_ptr = boost::make_shared<Range>();
      CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
          bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
 
      DM subdm;
      CHKERR DMCreate(mField.get_comm(), &subdm);
      CHKERR DMSetType(subdm, "DMMOFEM");
      CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_INIT");
      CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
      CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
      CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_TRUE);
 
      for (auto f : {"H", "G"}) {
        CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
        CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
      }
      CHKERR DMSetUp(subdm);
 
      if constexpr (debug) {
        if (mField.get_comm_size() == 1) {
          auto dm_ents = get_dofs_ents_all(SmartPetscObj<DM>(subdm, true));
          CHKERR save_range(mField.get_moab(), "sub_dm_init_ents_verts.h5m",
                            dm_ents.subset_by_type(MBVERTEX));
          dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
          CHKERR save_range(mField.get_moab(), "sub_dm_init_ents.h5m", dm_ents);
        }
      }
 
      return SmartPetscObj<DM>(subdm);
    };
 
    auto subdm = create_subdm();
 
    pip_mng->getDomainRhsFE().reset();
    pip_mng->getDomainLhsFE().reset();
    CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
    CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
    pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
    pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
 
    CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
    CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
 
    auto D = createDMVector(subdm);
    auto snes = pip_mng->createSNES(subdm);
    auto snes_ctx_ptr = getDMSnesCtx(subdm);
 
    auto set_section_monitor = [&](auto solver) {
      MoFEMFunctionBegin;
      CHKERR SNESMonitorSet(solver,
                            (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
                                               void *))MoFEMSNESMonitorFields,
                            (void *)(snes_ctx_ptr.get()), nullptr);
 
      MoFEMFunctionReturn(0);
    };
 
    auto print_section_field = [&]() {
      MoFEMFunctionBegin;
 
      auto section =
          mField.getInterface<ISManager>()->sectionCreate("SUB_INIT");
      PetscInt num_fields;
      CHKERR PetscSectionGetNumFields(section, &num_fields);
      for (int f = 0; f < num_fields; ++f) {
        const char *field_name;
        CHKERR PetscSectionGetFieldName(section, f, &field_name);
        MOFEM_LOG("FS", Sev::inform)
            << "Field " << f << " " << std::string(field_name);
      }
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_section_monitor(snes);
    CHKERR print_section_field();
 
    for (auto f : {"U", "P", "H", "G", "L"}) {
      CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
    }
 
    CHKERR SNESSetFromOptions(snes);
    auto B = createDMMatrix(subdm);
    CHKERR SNESSetJacobian(snes, B, B,
                                 PETSC_NULLPTR, PETSC_NULLPTR);
    CHKERR SNESSolve(snes, PETSC_NULLPTR, D);
 
    CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR reset_bits();
  CHKERR solve_init(
      [](FEMethod *fe_ptr) { return get_fe_bit(fe_ptr).test(0); });
  CHKERR refineMesh(refine_overlap);
 
  for (auto f : {"U", "P", "H", "G", "L"}) {
    CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
  }
  CHKERR solve_init([](FEMethod *fe_ptr) {
    return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
  });
 
  CHKERR post_proc([](FEMethod *fe_ptr) {
    return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
  });
 
  pip_mng->getOpDomainRhsPipeline().clear();
  pip_mng->getOpDomainLhsPipeline().clear();
 
  // Remove DOFs where boundary conditions are set
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
                                           "U", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
                                           "L", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
                                           "U", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
                                           "L", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "U",
                                           0, SPACE_DIM);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "L",
                                           0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "ZERO",
                                           "L", 0, 0);
 
  MoFEMFunctionReturn(0);
}

checkResults()

MoFEMErrorCode FreeSurface::checkResults ( )

private

Check results for correctness.

[Check]

Returns

MoFEMErrorCode Success/failure code

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 3650 of file free_surface.cpp.

                                         {
  MoFEMFunctionBegin;
  PetscInt test_nb = 0;
  PetscBool test_flg = PETSC_FALSE;
  double regression_value = 0.0;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-test", &test_nb, &test_flg);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-regression_value", &regression_value, PETSC_NULLPTR);
 
  if (test_flg) {
    auto simple = mField.getInterface<Simple>();
    auto T = createDMVector(simple->getDM());
    CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
                                   SCATTER_FORWARD);
    double nrm2;
    CHKERR VecNorm(T, NORM_2, &nrm2);
    MOFEM_LOG("FS", Sev::verbose) << "Regression norm " << nrm2;
    switch (test_nb) {
    case 1:
      break;
    default:
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong test number.");
      break;
    }
    if (fabs(nrm2 - regression_value) > 1e-2)
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
              "Regression test field; wrong norm value. %6.4e != %6.4e", nrm2,
              regression_value);
  }
  MoFEMFunctionReturn(0);
}

findEntitiesCrossedByPhaseInterface()

std::vector< Range > FreeSurface::findEntitiesCrossedByPhaseInterface ( size_t  overlap )

private

Find entities on refinement levels.

Parameters

overlap

Level of overlap around phase interface

Returns

Vector of entity ranges for each refinement level

Identifies mesh entities that cross the phase interface by analyzing phase field values at vertices. Returns entities that need refinement at each hierarchical level with specified overlap.

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 2924 of file free_surface.cpp.

                                                               {
 
  auto &moab = mField.get_moab();
  auto bit_mng = mField.getInterface<BitRefManager>();
  auto comm_mng = mField.getInterface<CommInterface>();
 
  Range vertices;
  CHK_THROW_MESSAGE(bit_mng->getEntitiesByTypeAndRefLevel(
                        bit(0), BitRefLevel().set(), MBVERTEX, vertices),
                    "can not get vertices on bit 0");
 
  auto &dofs_mi = mField.get_dofs()->get<Unique_mi_tag>();
  auto field_bit_number = mField.get_field_bit_number("H");
 
  Range plus_range, minus_range;
  std::vector<EntityHandle> plus, minus;
 
  // get vertices on level 0 on plus and minus side
  for (auto p = vertices.pair_begin(); p != vertices.pair_end(); ++p) {
 
    const auto f = p->first;
    const auto s = p->second;
 
    // Lowest Dof UId for given field (field bit number) on entity f
    const auto lo_uid = DofEntity::getLoFieldEntityUId(field_bit_number, f);
    const auto hi_uid = DofEntity::getHiFieldEntityUId(field_bit_number, s);
    auto it = dofs_mi.lower_bound(lo_uid);
    const auto hi_it = dofs_mi.upper_bound(hi_uid);
 
    plus.clear();
    minus.clear();
    plus.reserve(std::distance(it, hi_it));
    minus.reserve(std::distance(it, hi_it));
 
    for (; it != hi_it; ++it) {
      const auto v = (*it)->getFieldData();
      if (v > 0)
        plus.push_back((*it)->getEnt());
      else
        minus.push_back((*it)->getEnt());
    }
 
    plus_range.insert_list(plus.begin(), plus.end());
    minus_range.insert_list(minus.begin(), minus.end());
  }
 
  MOFEM_LOG_CHANNEL("SYNC");
  MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
      << "Plus range " << plus_range << endl;
  MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
      << "Minus range " << minus_range << endl;
  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
 
  auto get_elems = [&](auto &ents, auto bit, auto mask) {
    Range adj;
    CHK_MOAB_THROW(moab.get_adjacencies(ents, SPACE_DIM, false, adj,
                                        moab::Interface::UNION),
                   "can not get adjacencies");
    CHK_THROW_MESSAGE(bit_mng->filterEntitiesByRefLevel(bit, mask, adj),
                      "can not filter elements with bit 0");
    return adj;
  };
 
  CHKERR comm_mng->synchroniseEntities(plus_range);
  CHKERR comm_mng->synchroniseEntities(minus_range);
 
  std::vector<Range> ele_plus(nb_levels), ele_minus(nb_levels);
  ele_plus[0] = get_elems(plus_range, bit(0), BitRefLevel().set());
  ele_minus[0] = get_elems(minus_range, bit(0), BitRefLevel().set());
  auto common = intersect(ele_plus[0], ele_minus[0]);
  ele_plus[0] = subtract(ele_plus[0], common);
  ele_minus[0] = subtract(ele_minus[0], common);
 
  auto get_children = [&](auto &p, auto &c) {
    MoFEMFunctionBegin;
    CHKERR bit_mng->updateRangeByChildren(p, c);
    c = c.subset_by_dimension(SPACE_DIM);
    MoFEMFunctionReturn(0);
  };
 
  for (auto l = 1; l != nb_levels; ++l) {
    CHK_THROW_MESSAGE(get_children(ele_plus[l - 1], ele_plus[l]),
                      "get children");
    CHK_THROW_MESSAGE(get_children(ele_minus[l - 1], ele_minus[l]),
                      "get children");
  }
 
  auto get_level = [&](auto &p, auto &m, auto z, auto bit, auto mask) {
    Range l;
    CHK_THROW_MESSAGE(
        bit_mng->getEntitiesByDimAndRefLevel(bit, mask, SPACE_DIM, l),
        "can not get vertices on bit");
    l = subtract(l, p);
    l = subtract(l, m);
    for (auto f = 0; f != z; ++f) {
      Range conn;
      CHK_MOAB_THROW(moab.get_connectivity(l, conn, true), "");
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(conn);
      l = get_elems(conn, bit, mask);
    }
    return l;
  };
 
  std::vector<Range> vec_levels(nb_levels);
  for (auto l = nb_levels - 1; l >= 0; --l) {
    vec_levels[l] = get_level(ele_plus[l], ele_minus[l], 2 * overlap, bit(l),
                              BitRefLevel().set());
  }
 
  if constexpr (debug) {
    if (mField.get_comm_size() == 1) {
      for (auto l = 0; l != nb_levels; ++l) {
        std::string name = (boost::format("out_r%d.h5m") % l).str();
        CHK_THROW_MESSAGE(save_range(mField.get_moab(), name, vec_levels[l]),
                          "save mesh");
      }
    }
  }
 
  return vec_levels;
}

findParentsToRefine()

Range FreeSurface::findParentsToRefine ( Range  ents, BitRefLevel  level, BitRefLevel  mask  )

private

Find parent entities that need refinement.

Parameters

ents	Child entities requiring refinement
level	Bit level for entity marking
mask	Bit mask for filtering

Returns

Range of parent entities to refine

Traverses mesh hierarchy to find parent entities that should be refined to accommodate interface tracking requirements.

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

makeRefProblem()

MoFEMErrorCode FreeSurface::makeRefProblem

(

)

Create refined problem for mesh adaptation.

Returns

MoFEMErrorCode Success/failure code

Creates mesh refinement data structures needed for adaptive meshing

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

projectData()

MoFEMErrorCode FreeSurface::projectData ( )

private

Project solution data between mesh levels.

[Boundary condition]

Returns

MoFEMErrorCode Success/failure code

Projects field data from coarse to fine mesh levels during refinement. Handles both time stepping vectors (for theta method) and regular fields. Uses L2 projection with mass matrix assembly.

[Data projection]

Zero DOFs, used by FieldBlas

get TSTheta data operators

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 1402 of file free_surface.cpp.

                                        {
  MoFEMFunctionBegin;
 
  // FIXME: Functionality in this method should be improved (projection should
  // be done field by field), generalized, and move to become core
  // functionality.
 
  auto simple = mField.getInterface<Simple>();
  auto pip_mng = mField.getInterface<PipelineManager>();
  auto bit_mng = mField.getInterface<BitRefManager>();
  auto field_blas = mField.getInterface<FieldBlas>();
 
  // Store all existing elements pipelines, replace them by data projection
  // pipelines, to put them back when projection is done.
  auto fe_domain_lhs = pip_mng->getDomainLhsFE();
  auto fe_domain_rhs = pip_mng->getDomainRhsFE();
  auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
  auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
 
  pip_mng->getDomainLhsFE().reset();
  pip_mng->getDomainRhsFE().reset();
  pip_mng->getBoundaryLhsFE().reset();
  pip_mng->getBoundaryRhsFE().reset();
 
  using UDO = ForcesAndSourcesCore::UserDataOperator;
 
  // extract field data for domain parent element
  auto add_parent_field_domain = [&](auto fe, auto op, auto field, auto bit) {
    return setParentDofs(
        fe, field, op, bit,
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new DomainParentEle(mField));
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
  };
 
  // extract field data for boundary parent element
  auto add_parent_field_bdy = [&](auto fe, auto op, auto field, auto bit) {
    return setParentDofs(
        fe, field, op, bit,
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new BoundaryParentEle(mField));
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
  };
 
  // run this element on element with given bit, otherwise run other nested
  // element
  auto create_run_parent_op = [&](auto parent_fe_ptr, auto this_fe_ptr,
                                  auto fe_bit) {
    auto parent_mask = fe_bit;
    parent_mask.flip();
    return new OpRunParent(parent_fe_ptr, BitRefLevel().set(), parent_mask,
                           this_fe_ptr, fe_bit, BitRefLevel().set(), QUIET,
                           Sev::inform);
  };
 
  // create hierarchy of nested elements
  auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr, auto fe_bit) {
    std::vector<boost::shared_ptr<DomainParentEle>> parents_elems_ptr_vec;
    for (int l = 0; l < nb_levels; ++l)
      parents_elems_ptr_vec.emplace_back(
          boost::make_shared<DomainParentEle>(mField));
    for (auto l = 1; l < nb_levels; ++l) {
      parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
          create_run_parent_op(parents_elems_ptr_vec[l], this_fe_ptr, fe_bit));
    }
    return parents_elems_ptr_vec[0];
  };
 
  // solve projection problem
  auto solve_projection = [&](auto exe_test) {
    MoFEMFunctionBegin;
 
    auto set_domain_rhs = [&](auto fe) {
      MoFEMFunctionBegin;
      auto &pip = fe->getOpPtrVector();
 
      auto u_ptr = boost::make_shared<MatrixDouble>();
      auto p_ptr = boost::make_shared<VectorDouble>();
      auto h_ptr = boost::make_shared<VectorDouble>();
      auto g_ptr = boost::make_shared<VectorDouble>();
 
      auto eval_fe_ptr = boost::make_shared<DomainParentEle>(mField);
      {
        auto &pip = eval_fe_ptr->getOpPtrVector();
        CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
        CHKERR setParentDofs(
            eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
 
            [&]() {
              boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
                  new DomainParentEle(mField));
              return fe_parent;
            },
 
            QUIET, Sev::noisy);
        // That can be done much smarter, by block, field by field. For
        // simplicity is like that.
        CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPCOL, "U",
                                       bit(get_skin_projection_bit()));
        pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
        CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPCOL, "P",
                                       bit(get_skin_projection_bit()));
        pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
        CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPCOL, "H",
                                       bit(get_skin_projection_bit()));
        pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
        CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPCOL, "G",
                                       bit(get_skin_projection_bit()));
        pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
      }
      auto parent_eval_fe_ptr =
          get_parents_vel_fe_ptr(eval_fe_ptr, bit(get_projection_bit()));
      pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
                                         bit(get_projection_bit())));
 
      auto assemble_fe_ptr = boost::make_shared<DomainParentEle>(mField);
      {
        auto &pip = assemble_fe_ptr->getOpPtrVector();
        CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
        CHKERR setParentDofs(
            assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
            [&]() {
              boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
                  new DomainParentEle(mField));
              return fe_parent;
            },
 
            QUIET, Sev::noisy);
        CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "U",
                                       bit(get_skin_parent_bit()));
        pip.push_back(new OpDomainAssembleVector("U", u_ptr));
        CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "P",
                                       bit(get_skin_parent_bit()));
        pip.push_back(new OpDomainAssembleScalar("P", p_ptr));
        CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "H",
                                       bit(get_skin_parent_bit()));
        pip.push_back(new OpDomainAssembleScalar("H", h_ptr));
        CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "G",
                                       bit(get_skin_parent_bit()));
        pip.push_back(new OpDomainAssembleScalar("G", g_ptr));
      }
      auto parent_assemble_fe_ptr =
          get_parents_vel_fe_ptr(assemble_fe_ptr, bit(get_current_bit()));
      pip.push_back(create_run_parent_op(
          parent_assemble_fe_ptr, assemble_fe_ptr, bit(get_current_bit())));
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_domain_lhs = [&](auto fe) {
      MoFEMFunctionBegin;
 
      auto &pip = fe->getOpPtrVector();
 
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
 
      CHKERR setParentDofs(
          fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
          [&]() {
            boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
                new DomainParentEle(mField));
            return fe_parent;
          },
 
          QUIET, Sev::noisy);
 
      // That can be done much smarter, by block, field by field. For simplicity
      // is like that.
      CHKERR add_parent_field_domain(fe, UDO::OPROW, "U",
                                     bit(get_skin_parent_bit()));
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U",
                                     bit(get_skin_parent_bit()));
      pip.push_back(new OpDomainMassU("U", "U"));
      CHKERR add_parent_field_domain(fe, UDO::OPROW, "P",
                                     bit(get_skin_parent_bit()));
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P",
                                     bit(get_skin_parent_bit()));
      pip.push_back(new OpDomainMassP("P", "P"));
      CHKERR add_parent_field_domain(fe, UDO::OPROW, "H",
                                     bit(get_skin_parent_bit()));
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H",
                                     bit(get_skin_parent_bit()));
      pip.push_back(new OpDomainMassH("H", "H"));
      CHKERR add_parent_field_domain(fe, UDO::OPROW, "G",
                                     bit(get_skin_parent_bit()));
      CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G",
                                     bit(get_skin_parent_bit()));
      pip.push_back(new OpDomainMassG("G", "G"));
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_bdy_rhs = [&](auto fe) {
      MoFEMFunctionBegin;
      auto &pip = fe->getOpPtrVector();
 
      auto l_ptr = boost::make_shared<VectorDouble>();
 
      auto eval_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
      {
        auto &pip = eval_fe_ptr->getOpPtrVector();
        CHKERR setParentDofs(
            eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
 
            [&]() {
              boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
                  new BoundaryParentEle(mField));
              return fe_parent;
            },
 
            QUIET, Sev::noisy);
        // That can be done much smarter, by block, field by field. For
        // simplicity is like that.
        CHKERR add_parent_field_bdy(eval_fe_ptr, UDO::OPCOL, "L",
                                    bit(get_skin_projection_bit()));
        pip.push_back(new OpCalculateScalarFieldValues("L", l_ptr));
      }
      auto parent_eval_fe_ptr =
          get_parents_vel_fe_ptr(eval_fe_ptr, bit(get_projection_bit()));
      pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
                                         bit(get_projection_bit())));
 
      auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
      {
        auto &pip = assemble_fe_ptr->getOpPtrVector();
        CHKERR setParentDofs(
            assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
            [&]() {
              boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
                  new BoundaryParentEle(mField));
              return fe_parent;
            },
 
            QUIET, Sev::noisy);
 
        struct OpLSize : public BoundaryEleOp {
          OpLSize(boost::shared_ptr<VectorDouble> l_ptr)
              : BoundaryEleOp(NOSPACE, DomainEleOp::OPSPACE), lPtr(l_ptr) {}
          MoFEMErrorCode doWork(int, EntityType, EntData &) {
            MoFEMFunctionBegin;
            if (lPtr->size() != getGaussPts().size2()) {
              lPtr->resize(getGaussPts().size2());
              lPtr->clear();
            }
            MoFEMFunctionReturn(0);
          }
 
        private:
          boost::shared_ptr<VectorDouble> lPtr;
        };
 
        pip.push_back(new OpLSize(l_ptr));
 
        CHKERR add_parent_field_bdy(assemble_fe_ptr, UDO::OPROW, "L",
                                    bit(get_skin_parent_bit()));
        pip.push_back(new OpBoundaryAssembleScalar("L", l_ptr));
      }
      auto parent_assemble_fe_ptr =
          get_parents_vel_fe_ptr(assemble_fe_ptr, bit(get_current_bit()));
      pip.push_back(create_run_parent_op(
          parent_assemble_fe_ptr, assemble_fe_ptr, bit(get_current_bit())));
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_bdy_lhs = [&](auto fe) {
      MoFEMFunctionBegin;
 
      auto &pip = fe->getOpPtrVector();
 
      CHKERR setParentDofs(
          fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
          [&]() {
            boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
                new BoundaryParentEle(mField));
            return fe_parent;
          },
 
          QUIET, Sev::noisy);
 
      // That can be done much smarter, by block, field by field. For simplicity
      // is like that.
      CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L",
                                  bit(get_skin_parent_bit()));
      CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "L",
                                  bit(get_skin_parent_bit()));
      pip.push_back(new OpBoundaryMassL("L", "L"));
 
      MoFEMFunctionReturn(0);
    };
 
    auto create_subdm = [&]() -> SmartPetscObj<DM> {
      auto level_ents_ptr = boost::make_shared<Range>();
      CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
          bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
 
      auto get_prj_ents = [&]() {
        Range prj_mesh;
        CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_projection_bit()),
                                                    BitRefLevel().set(),
                                                    SPACE_DIM, prj_mesh);
        auto common_ents = intersect(prj_mesh, *level_ents_ptr);
        prj_mesh = subtract(unite(*level_ents_ptr, prj_mesh), common_ents)
                       .subset_by_dimension(SPACE_DIM);
 
        return prj_mesh;
      };
 
      auto prj_ents = get_prj_ents();
 
      if (get_global_size(prj_ents.size())) {
 
        auto rebuild = [&]() {
          auto prb_mng = mField.getInterface<ProblemsManager>();
          MoFEMFunctionBegin;
 
          std::vector<std::string> fields{"U", "P", "H", "G", "L"};
          std::map<std::string, boost::shared_ptr<Range>> range_maps{
 
              {"U", level_ents_ptr},
              {"P", level_ents_ptr},
              {"H", level_ents_ptr},
              {"G", level_ents_ptr},
              {"L", level_ents_ptr}
 
          };
 
          CHKERR prb_mng->buildSubProblem("SUB_SOLVER", fields, fields,
                                          simple->getProblemName(), PETSC_TRUE,
                                          &range_maps, &range_maps);
 
          // partition problem
          CHKERR prb_mng->partitionFiniteElements("SUB_SOLVER", true, 0,
                                                  mField.get_comm_size());
          // set ghost nodes
          CHKERR prb_mng->partitionGhostDofsOnDistributedMesh("SUB_SOLVER");
 
          MoFEMFunctionReturn(0);
        };
 
        MOFEM_LOG("FS", Sev::verbose) << "Create projection problem";
 
        CHKERR rebuild();
 
        auto dm = simple->getDM();
        DM subdm;
        CHKERR DMCreate(mField.get_comm(), &subdm);
        CHKERR DMSetType(subdm, "DMMOFEM");
        CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
        return SmartPetscObj<DM>(subdm);
      }
 
      MOFEM_LOG("FS", Sev::inform) << "Nothing to project";
 
      return SmartPetscObj<DM>();
    };
 
    auto create_dummy_dm = [&]() {
      auto dummy_dm = createDM(mField.get_comm(), "DMMOFEM");
      CHK_THROW_MESSAGE(DMMoFEMCreateMoFEM(dummy_dm, &mField,
                                           simple->getProblemName().c_str(),
                                           BitRefLevel(), BitRefLevel()),
                        "create dummy dm");
      return dummy_dm;
    };
 
    auto subdm = create_subdm();
    if (subdm) {
 
      pip_mng->getDomainRhsFE().reset();
      pip_mng->getDomainLhsFE().reset();
      pip_mng->getBoundaryRhsFE().reset();
      pip_mng->getBoundaryLhsFE().reset();
      CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
      CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
      CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
      CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
      pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
      pip_mng->getDomainRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
        return get_fe_bit(fe_ptr).test(nb_levels - 1);
      };
      pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
      pip_mng->getBoundaryRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
        return get_fe_bit(fe_ptr).test(nb_levels - 1);
      };
 
      CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
      CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
      CHKERR set_bdy_rhs(pip_mng->getCastBoundaryRhsFE());
      CHKERR set_bdy_lhs(pip_mng->getCastBoundaryLhsFE());
 
      auto D = createDMVector(subdm);
      auto F = vectorDuplicate(D);
 
      auto ksp = pip_mng->createKSP(subdm);
      CHKERR KSPSetFromOptions(ksp);
      CHKERR KSPSetUp(ksp);
 
      // get vector norm
      auto get_norm = [&](auto x) {
        double nrm;
        CHKERR VecNorm(x, NORM_2, &nrm);
        return nrm;
      };
 
      /**
       * @brief Zero DOFs, used by FieldBlas
       */
      auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
        MoFEMFunctionBegin;
        for (auto &v : ent_ptr->getEntFieldData()) {
          v = 0;
        }
        MoFEMFunctionReturn(0);
      };
 
      auto solve = [&](auto S) {
        MoFEMFunctionBegin;
        CHKERR VecZeroEntries(S);
        CHKERR VecZeroEntries(F);
        CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR KSPSolve(ksp, F, S);
        CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
 
 
 
        MoFEMFunctionReturn(0);
      };
 
      MOFEM_LOG("FS", Sev::inform) << "Solve projection";
      CHKERR solve(D);
 
      auto glob_x = createDMVector(simple->getDM());
      auto sub_x = createDMVector(subdm);
      auto dummy_dm = create_dummy_dm();
 
      /**
       * @brief get TSTheta data operators
       */
      auto apply_restrict = [&]() {
        auto get_is = [](auto v) {
          IS iy;
          auto create = [&]() {
            MoFEMFunctionBegin;
            int n, ystart;
            CHKERR VecGetLocalSize(v, &n);
            CHKERR VecGetOwnershipRange(v, &ystart, NULL);
            CHKERR ISCreateStride(PETSC_COMM_SELF, n, ystart, 1, &iy);
            MoFEMFunctionReturn(0);
          };
          CHK_THROW_MESSAGE(create(), "create is");
          return SmartPetscObj<IS>(iy);
        };
 
        auto iy = get_is(glob_x);
        auto s = createVecScatter(glob_x, PETSC_NULLPTR, glob_x, iy);
 
        CHK_THROW_MESSAGE(
            DMSubDomainRestrict(simple->getDM(), s, PETSC_NULLPTR, dummy_dm),
            "restrict");
        Vec X0, Xdot;
        CHK_THROW_MESSAGE(DMGetNamedGlobalVector(dummy_dm, "TSTheta_X0", &X0),
                          "get X0");
        CHK_THROW_MESSAGE(
            DMGetNamedGlobalVector(dummy_dm, "TSTheta_Xdot", &Xdot),
            "get Xdot");
 
        auto forward_ghost = [](auto g) {
          MoFEMFunctionBegin;
          CHKERR VecGhostUpdateBegin(g, INSERT_VALUES, SCATTER_FORWARD);
          CHKERR VecGhostUpdateEnd(g, INSERT_VALUES, SCATTER_FORWARD);
          MoFEMFunctionReturn(0);
        };
 
        CHK_THROW_MESSAGE(forward_ghost(X0), "");
        CHK_THROW_MESSAGE(forward_ghost(Xdot), "");
 
        if constexpr (debug) {
          MOFEM_LOG("FS", Sev::inform)
              << "Reverse restrict: X0 " << get_norm(X0) << " Xdot "
              << get_norm(Xdot);
        }
 
        return std::vector<Vec>{X0, Xdot};
      };
 
      auto ts_solver_vecs = apply_restrict();
 
      if (ts_solver_vecs.size()) {
 
        for (auto v : ts_solver_vecs) {
          MOFEM_LOG("FS", Sev::inform) << "Solve projection vector";
 
          CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
                                         SCATTER_REVERSE);
          CHKERR solve(sub_x);
 
          for (auto f : {"U", "P", "H", "G", "L"}) {
            MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
            CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
          }
 
          CHKERR DMoFEMMeshToLocalVector(subdm, sub_x, INSERT_VALUES,
                                         SCATTER_REVERSE);
          CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
                                         SCATTER_FORWARD);
 
          MOFEM_LOG("FS", Sev::inform) << "Norm V " << get_norm(v);
        }
 
        CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_X0",
                                          &ts_solver_vecs[0]);
        CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_Xdot",
                                          &ts_solver_vecs[1]);
      }
 
      for (auto f : {"U", "P", "H", "G", "L"}) {
        MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
        CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
      }
      CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  // postprocessing (only for debugging purposes)
  auto post_proc = [&](auto exe_test) {
    MoFEMFunctionBegin;
    auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
    auto &pip = post_proc_fe->getOpPtrVector();
    post_proc_fe->exeTestHook = exe_test;
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
 
    CHKERR setParentDofs(
        post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new DomainParentEle(mField));
          AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
              fe_parent->getOpPtrVector(), {H1});
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    auto u_ptr = boost::make_shared<MatrixDouble>();
    auto p_ptr = boost::make_shared<VectorDouble>();
    auto h_ptr = boost::make_shared<VectorDouble>();
    auto g_ptr = boost::make_shared<VectorDouble>();
 
    CHKERR add_parent_field_domain(post_proc_fe, UDO::OPCOL, "U",
                                   bit(get_skin_parent_bit()));
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
    CHKERR add_parent_field_domain(post_proc_fe, UDO::OPCOL, "P",
                                   bit(get_skin_parent_bit()));
    pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
    CHKERR add_parent_field_domain(post_proc_fe, UDO::OPCOL, "H",
                                   bit(get_skin_parent_bit()));
    pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
    CHKERR add_parent_field_domain(post_proc_fe, UDO::OPCOL, "G",
                                   bit(get_skin_parent_bit()));
    pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(post_proc_fe->getPostProcMesh(),
                    post_proc_fe->getMapGaussPts(),
 
                    {{"P", p_ptr}, {"H", h_ptr}, {"G", g_ptr}},
 
                    {{"U", u_ptr}},
 
                    {},
 
                    {}
 
                    )
 
    );
 
    auto dm = simple->getDM();
    CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
    CHKERR post_proc_fe->writeFile("out_projection.h5m");
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR solve_projection([](FEMethod *fe_ptr) {
    return get_fe_bit(fe_ptr).test(get_current_bit());
  });
 
  if constexpr (debug) {
    CHKERR post_proc([](FEMethod *fe_ptr) {
      return get_fe_bit(fe_ptr).test(get_current_bit());
    });
  }
 
  fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
  fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
  fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
  fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
 
  MoFEMFunctionReturn(0);
}

readMesh()

MoFEMErrorCode FreeSurface::readMesh ( )

private

Read mesh from input file.

[Run programme]

Returns

MoFEMErrorCode Success/failure code

Loads mesh using Simple interface and sets up parent adjacencies for hierarchical mesh refinement

[Read mesh]

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 988 of file free_surface.cpp.

                                     {
  MoFEMFunctionBegin;
  MOFEM_LOG("FS", Sev::inform) << "Read mesh for problem";
  auto simple = mField.getInterface<Simple>();
 
  simple->getParentAdjacencies() = true;
  simple->getBitRefLevel() = BitRefLevel();
 
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
 
  MoFEMFunctionReturn(0);
}

refineMesh()

MoFEMErrorCode FreeSurface::refineMesh ( size_t  overlap )

private

Perform adaptive mesh refinement.

Parameters

overlap

Number of element layers around interface to refine

Returns

MoFEMErrorCode Success/failure code

Performs hierarchical mesh refinement around the phase interface:

• Identifies entities crossing interface
• Creates refinement hierarchy
  
• Sets up skin elements between levels
• Updates bit level markings for computation

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 3046 of file free_surface.cpp.

                                                     {
  MoFEMFunctionBegin;
 
  auto bit_mng = mField.getInterface<BitRefManager>();
 
  BitRefLevel start_mask;
  for (auto s = 0; s != get_start_bit(); ++s)
    start_mask[s] = true;
 
  // store prev_level
  Range prev_level;
  CHKERR bit_mng->getEntitiesByRefLevel(bit(get_current_bit()),
                                        BitRefLevel().set(), prev_level);
  Range prev_level_skin;
  CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_parent_bit()),
                                        BitRefLevel().set(), prev_level_skin);
  // reset bit ref levels
  CHKERR bit_mng->lambdaBitRefLevel(
      [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
  CHKERR bit_mng->setNthBitRefLevel(prev_level, get_projection_bit(), true);
  CHKERR bit_mng->setNthBitRefLevel(prev_level_skin, get_skin_projection_bit(),
                                    true);
 
  // set refinement levels
  auto set_levels = [&](auto &&
                            vec_levels /*entities are refined on each level*/) {
    MoFEMFunctionBegin;
 
    // start with zero level, which is the coarsest mesh
    Range level0;
    CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
    CHKERR bit_mng->setNthBitRefLevel(level0, get_start_bit(), true);
 
    // get lower dimension entities
    auto get_adj = [&](auto ents) {
      Range conn;
      CHK_MOAB_THROW(mField.get_moab().get_connectivity(ents, conn, true),
                     "get conn");
      for (auto d = 1; d != SPACE_DIM; ++d) {
        CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
                           ents.subset_by_dimension(SPACE_DIM), d, false, conn,
                           moab::Interface::UNION),
                       "get adj");
      }
      ents.merge(conn);
      return ents;
    };
 
    // set bit levels
    for (auto l = 1; l != nb_levels; ++l) {
      Range level_prev;
      CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_start_bit() + l - 1),
                                                  BitRefLevel().set(),
                                                  SPACE_DIM, level_prev);
      Range parents;
      CHKERR bit_mng->updateRangeByParent(vec_levels[l], parents);
      // subtract entities from previous level, which are refined, so should be
      // not there
      level_prev = subtract(level_prev, parents);
      // and instead add their children
      level_prev.merge(vec_levels[l]);
      // set bit to each level
      CHKERR bit_mng->setNthBitRefLevel(level_prev, get_start_bit() + l, true);
    }
 
    // set bit levels to lower dimension entities
    for (auto l = 1; l != nb_levels; ++l) {
      Range level;
      CHKERR bit_mng->getEntitiesByDimAndRefLevel(
          bit(get_start_bit() + l), BitRefLevel().set(), SPACE_DIM, level);
      level = get_adj(level);
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(level);
      CHKERR bit_mng->setNthBitRefLevel(level, get_start_bit() + l, true);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  // resolve skin between refined levels
  auto set_skins = [&]() {
    MoFEMFunctionBegin;
 
    moab::Skinner skinner(&mField.get_moab());
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
 
    // get skin of bit level
    auto get_bit_skin = [&](BitRefLevel bit, BitRefLevel mask) {
      Range bit_ents;
      CHK_THROW_MESSAGE(
          mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
              bit, mask, SPACE_DIM, bit_ents),
          "can't get bit level");
      Range bit_skin;
      CHK_MOAB_THROW(skinner.find_skin(0, bit_ents, false, bit_skin),
                     "can't get skin");
      return bit_skin;
    };
 
    auto get_level_skin = [&]() {
      Range skin;
      BitRefLevel bit_prev;
      for (auto l = 1; l != nb_levels; ++l) {
        auto skin_level_mesh = get_bit_skin(bit(l), BitRefLevel().set());
        // filter (remove) all entities which are on partition borders
        CHKERR pcomm->filter_pstatus(skin_level_mesh,
                                     PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                     PSTATUS_NOT, -1, nullptr);
        auto skin_level =
            get_bit_skin(bit(get_start_bit() + l), BitRefLevel().set());
        skin_level = subtract(skin_level,
                              skin_level_mesh); // get only internal skins, not
                                                // on the body boundary
        // get lower dimension adjacencies (FIXME: add edges if 3D)
        Range skin_level_verts;
        CHKERR mField.get_moab().get_connectivity(skin_level, skin_level_verts,
                                                  true);
        skin_level.merge(skin_level_verts);
 
        // remove previous level
        bit_prev.set(l - 1);
        Range level_prev;
        CHKERR bit_mng->getEntitiesByRefLevel(bit_prev, BitRefLevel().set(),
                                              level_prev);
        skin.merge(subtract(skin_level, level_prev));
      }
 
      return skin;
    };
 
    auto resolve_shared = [&](auto &&skin) {
      Range tmp_skin = skin;
 
      map<int, Range> map_procs;
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          tmp_skin, &map_procs);
 
      Range from_other_procs; // entities which also exist on other processors
      for (auto &m : map_procs) {
        if (m.first != mField.get_comm_rank()) {
          from_other_procs.merge(m.second);
        }
      }
 
      auto common = intersect(
          skin, from_other_procs); // entities which are on internal skin
      skin.merge(from_other_procs);
 
      // entities which are on processor borders, and several processors are not
      // true skin.
      if (!common.empty()) {
        // skin is internal exist on other procs
        skin = subtract(skin, common);
      }
 
      return skin;
    };
 
    // get parents of entities
    auto get_parent_level_skin = [&](auto skin) {
      Range skin_parents;
      CHKERR bit_mng->updateRangeByParent(
          skin.subset_by_dimension(SPACE_DIM - 1), skin_parents);
      Range skin_parent_verts;
      CHKERR mField.get_moab().get_connectivity(skin_parents, skin_parent_verts,
                                                true);
      skin_parents.merge(skin_parent_verts);
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          skin_parents);
      return skin_parents;
    };
 
    auto child_skin = resolve_shared(get_level_skin());
    auto parent_skin = get_parent_level_skin(child_skin);
 
    child_skin = subtract(child_skin, parent_skin);
    CHKERR bit_mng->setNthBitRefLevel(child_skin, get_skin_child_bit(), true);
    CHKERR bit_mng->setNthBitRefLevel(parent_skin, get_skin_parent_bit(), true);
 
    MoFEMFunctionReturn(0);
  };
 
  // take last level, remove childs on boarder, and set bit
  auto set_current = [&]() {
    MoFEMFunctionBegin;
    Range last_level;
    CHKERR bit_mng->getEntitiesByRefLevel(bit(get_start_bit() + nb_levels - 1),
                                          BitRefLevel().set(), last_level);
    Range skin_child;
    CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_child_bit()),
                                          BitRefLevel().set(), skin_child);
 
    last_level = subtract(last_level, skin_child);
    CHKERR bit_mng->setNthBitRefLevel(last_level, get_current_bit(), true);
    MoFEMFunctionReturn(0);
  };
 
  // set bits to levels
  CHKERR set_levels(findEntitiesCrossedByPhaseInterface(overlap));
  // set bits to skin
  CHKERR set_skins();
  // set current level bit
  CHKERR set_current();
 
  if constexpr (debug) {
    if (mField.get_comm_size() == 1) {
      for (auto l = 0; l != nb_levels; ++l) {
        std::string name = (boost::format("out_level%d.h5m") % l).str();
        CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_start_bit() + l),
                                      BitRefLevel().set(), name.c_str(), "MOAB",
                                      "PARALLEL=WRITE_PART");
      }
      CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_current_bit()),
                                    BitRefLevel().set(), "current_bit.h5m",
                                    "MOAB", "PARALLEL=WRITE_PART");
      CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_projection_bit()),
                                    BitRefLevel().set(), "projection_bit.h5m",
                                    "MOAB", "PARALLEL=WRITE_PART");
 
      CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_child_bit()),
                                    BitRefLevel().set(), "skin_child_bit.h5m",
                                    "MOAB", "PARALLEL=WRITE_PART");
      CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_parent_bit()),
                                    BitRefLevel().set(), "skin_parent_bit.h5m",
                                    "MOAB", "PARALLEL=WRITE_PART");
    }
  }
 
  MoFEMFunctionReturn(0);
};

runProblem()

MoFEMErrorCode FreeSurface::runProblem

(

)

Main function to run the complete free surface simulation.

[Run programme]

Returns

MoFEMErrorCode Success/failure code

Executes the complete simulation workflow:

1. Read mesh from file
2. Setup problem (fields, approximation orders, parameters)
3. Apply boundary conditions and initialize fields
4. Assemble system matrices and operators
5. Solve time-dependent problem

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 975 of file free_surface.cpp.

                                       {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}

setParentDofs()

MoFEMErrorCode FreeSurface::setParentDofs ( boost::shared_ptr< FEMethod >  fe_top, std::string  field_name, ForcesAndSourcesCore::UserDataOperator::OpType  op, BitRefLevel  child_ent_bit, boost::function< boost::shared_ptr< ForcesAndSourcesCore >()>  get_elem, int  verbosity, LogManager::SeverityLevel  sev  )

private

Create hierarchy of elements run on parent levels.

Parameters

fe_top	Pipeline element to which hierarchy is attached
field_name	Name of field for DOF extraction ("" for all fields)
op	Type of operator OPSPACE, OPROW, OPCOL or OPROWCOL
child_ent_bit	Bit level marking child entities
get_elem	Lambda function to create element on hierarchy
verbosity	Verbosity level for debugging output
sev	Severity level for logging

Returns

MoFEMErrorCode Success/failure code

Sets up parent-child relationships for hierarchical mesh refinement. Allows access to field data from parent mesh levels during computation on refined child meshes. Essential for projection and interpolation.

Collect data from parent elements to child

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 3277 of file free_surface.cpp.

                                                {
  MoFEMFunctionBegin;
 
  /**
   * @brief Collect data from parent elements to child
   */
  boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
      add_parent_level =
          [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
            // Evaluate if not last parent element
            if (level > 0) {
 
              // Create domain parent FE
              auto fe_ptr_current = get_elem();
 
              // Call next level
              add_parent_level(
                  boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
                      fe_ptr_current),
                  level - 1);
 
              // Add data to curent fe level
              if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
 
                // Only base
                parent_fe_pt->getOpPtrVector().push_back(
 
                    new OpAddParentEntData(
 
                        H1, op, fe_ptr_current,
 
                        BitRefLevel().set(), BitRefLevel().set(),
 
                        child_ent_bit, BitRefLevel().set(),
 
                        verbosity, sev));
 
              } else {
 
                // Filed data
                parent_fe_pt->getOpPtrVector().push_back(
 
                    new OpAddParentEntData(
 
                        field_name, op, fe_ptr_current,
 
                        BitRefLevel().set(), BitRefLevel().set(),
 
                        child_ent_bit, BitRefLevel().set(),
 
                        verbosity, sev));
              }
            }
          };
 
  add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
                   nb_levels);
 
  MoFEMFunctionReturn(0);
}

setupProblem()

MoFEMErrorCode FreeSurface::setupProblem ( )

private

Setup problem fields and parameters.

[Read mesh]

Returns

MoFEMErrorCode Success/failure code

Creates finite element fields:

• U: Velocity field (H1, vector)
• P: Pressure field (H1, scalar)
• H: Phase/order field (H1, scalar)
• G: Chemical potential (H1, scalar)
• L: Lagrange multiplier for boundary conditions (H1, scalar)

Sets approximation orders and reads physical parameters from command line

[Set up problem]

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 1004 of file free_surface.cpp.

                                         {
  MoFEMFunctionBegin;
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-nb_levels", &nb_levels,
                            PETSC_NULLPTR);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-refine_overlap", &refine_overlap,
                            PETSC_NULLPTR);
 
  const char *coord_type_names[] = {"cartesian", "polar", "cylindrical",
                                    "spherical"};
  CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-coords", coord_type_names,
                              LAST_COORDINATE_SYSTEM, &coord_type, PETSC_NULLPTR);
 
  MOFEM_LOG("FS", Sev::inform) << "Approximation order = " << order;
  MOFEM_LOG("FS", Sev::inform)
      << "Number of refinement levels nb_levels = " << nb_levels;
  nb_levels += 1;
 
  auto simple = mField.getInterface<Simple>();
 
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-a0", &a0, PETSC_NULLPTR); // Acceleration
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-rho_m", &rho_m, PETSC_NULLPTR); // Density minus phase
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-mu_m", &mu_m, PETSC_NULLPTR); // Viscosity minus phase
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-rho_p", &rho_p, PETSC_NULLPTR); // Density plus phase
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-mu_p", &mu_p, PETSC_NULLPTR); // Viscosity plus phase
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-lambda", &lambda, PETSC_NULLPTR); // Surface tension
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-W", &W, PETSC_NULLPTR); // Height of the well in energy functional
 
  rho_ave = (rho_p + rho_m) / 2;
  rho_diff = (rho_p - rho_m) / 2;
  mu_ave = (mu_p + mu_m) / 2;
  mu_diff = (mu_p - mu_m) / 2;
 
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-h", &h, PETSC_NULLPTR);
  eta = h;
  eta2 = eta * eta;
  kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
 
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-md", &md, PETSC_NULLPTR);
 
  MOFEM_LOG("FS", Sev::inform) << "Acceleration a0 = " << a0;
  MOFEM_LOG("FS", Sev::inform) << "\"Minus\" phase density rho_m = " << rho_m;
  MOFEM_LOG("FS", Sev::inform) << "\"Minus\" phase viscosity mu_m = " << mu_m;
  MOFEM_LOG("FS", Sev::inform) << "\"Plus\" phase density rho_p = " << rho_p;
  MOFEM_LOG("FS", Sev::inform) << "\"Plus\" phase viscosity mu_p = " << mu_p;
  MOFEM_LOG("FS", Sev::inform) << "Surface tension lambda = " << lambda;
  MOFEM_LOG("FS", Sev::inform)
      << "Height of the well in energy functional W = " << W;
  MOFEM_LOG("FS", Sev::inform) << "Characteristic mesh size h = " << h;
  MOFEM_LOG("FS", Sev::inform) << "Mobility md = " << md;
 
  MOFEM_LOG("FS", Sev::inform) << "Average density rho_ave = " << rho_ave;
  MOFEM_LOG("FS", Sev::inform) << "Difference density rho_diff = " << rho_diff;
  MOFEM_LOG("FS", Sev::inform) << "Average viscosity mu_ave = " << mu_ave;
  MOFEM_LOG("FS", Sev::inform) << "Difference viscosity mu_diff = " << mu_diff;
  MOFEM_LOG("FS", Sev::inform) << "kappa = " << kappa;
 
 
  // Fields on domain
 
  // Velocity field
  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, U_FIELD_DIM);
  // Pressure field
  CHKERR simple->addDomainField("P", H1, AINSWORTH_LEGENDRE_BASE, 1);
  // Order/phase fild
  CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
  // Chemical potential
  CHKERR simple->addDomainField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
 
  // Field on boundary
  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE,
                                  U_FIELD_DIM);
  CHKERR simple->addBoundaryField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
  CHKERR simple->addBoundaryField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
  // Lagrange multiplier which constrains slip conditions
  CHKERR simple->addBoundaryField("L", H1, AINSWORTH_LEGENDRE_BASE, 1);
 
  CHKERR simple->setFieldOrder("U", order);
  CHKERR simple->setFieldOrder("P", order - 1);
  CHKERR simple->setFieldOrder("H", order);
  CHKERR simple->setFieldOrder("G", order);
  CHKERR simple->setFieldOrder("L", order);
 
  // Initialise bit ref levels
  auto set_problem_bit = [&]() {
    MoFEMFunctionBegin;
    // 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();
    simple->getBitRefLevelMask() = BitRefLevel().set();
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_problem_bit();
 
  CHKERR simple->setUp();
 
  MoFEMFunctionReturn(0);
}

solveSystem()

MoFEMErrorCode FreeSurface::solveSystem ( )

private

Solve the time-dependent free surface problem.

[Solve]

Returns

MoFEMErrorCode Success/failure code

Creates and configures time stepping solver with:

• Sub-problem DM for refined mesh
• Post-processing for output visualization
  
• Custom preconditioner and monitors
• TSTheta implicit time integration

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 2498 of file free_surface.cpp.

                                        {
  MoFEMFunctionBegin;
 
  using UDO = ForcesAndSourcesCore::UserDataOperator;
 
  auto simple = mField.getInterface<Simple>();
  auto pip_mng = mField.getInterface<PipelineManager>();
 
  auto create_solver_dm = [&](auto dm) -> SmartPetscObj<DM> {
    DM subdm;
 
    auto setup_subdm = [&](auto dm) {
      MoFEMFunctionBegin;
      auto simple = mField.getInterface<Simple>();
      auto bit_mng = mField.getInterface<BitRefManager>();
      auto dm = simple->getDM();
      auto level_ents_ptr = boost::make_shared<Range>();
      CHKERR bit_mng->getEntitiesByRefLevel(
          bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
      CHKERR DMCreate(mField.get_comm(), &subdm);
      CHKERR DMSetType(subdm, "DMMOFEM");
      CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
      CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
      CHKERR DMMoFEMAddElement(subdm, simple->getBoundaryFEName());
      CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
      CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_FALSE);
      for (auto f : {"U", "P", "H", "G", "L"}) {
        CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
        CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
      }
      CHKERR DMSetUp(subdm);
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(setup_subdm(dm), "create subdm");
 
    return SmartPetscObj<DM>(subdm);
  };
 
  auto get_fe_post_proc = [&](auto post_proc_mesh) {
    auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
      return setParentDofs(
          fe, field, op, bit(get_skin_parent_bit()),
 
          [&]() {
            boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
                new DomainParentEle(mField));
            return fe_parent;
          },
 
          QUIET, Sev::noisy);
    };
 
    auto post_proc_fe =
        boost::make_shared<PostProcEleDomainCont>(mField, post_proc_mesh);
    post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
      return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
          get_start_bit() + nb_levels - 1);
    };
 
    auto u_ptr = boost::make_shared<MatrixDouble>();
    auto grad_u_ptr = boost::make_shared<MatrixDouble>();
    auto h_ptr = boost::make_shared<VectorDouble>();
    auto grad_h_ptr = boost::make_shared<MatrixDouble>();
    auto p_ptr = boost::make_shared<VectorDouble>();
    auto g_ptr = boost::make_shared<VectorDouble>();
    auto grad_g_ptr = boost::make_shared<MatrixDouble>();
 
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        post_proc_fe->getOpPtrVector(), {H1});
 
    CHKERR setParentDofs(
        post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new DomainParentEle(mField));
          AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
              fe_parent->getOpPtrVector(), {H1});
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
 
    CHKERR add_parent_field_domain(post_proc_fe, UDO::OPCOL, "U");
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateVectorFieldGradient<U_FIELD_DIM, SPACE_DIM>("U",
                                                                   grad_u_ptr));
 
    CHKERR add_parent_field_domain(post_proc_fe, UDO::OPCOL, "H");
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("H", h_ptr));
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
 
    CHKERR add_parent_field_domain(post_proc_fe, UDO::OPCOL, "P");
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("P", p_ptr));
 
    CHKERR add_parent_field_domain(post_proc_fe, UDO::OPCOL, "G");
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("G", g_ptr));
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
 
    using OpPPMap = OpPostProcMapInMoab<2, 2>;
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(
            post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
            {{"H", h_ptr}, {"P", p_ptr}, {"G", g_ptr}},
 
            {{"U", u_ptr}, {"H_GRAD", grad_h_ptr}, {"G_GRAD", grad_g_ptr}},
 
            {{"GRAD_U", grad_u_ptr}},
 
            {}
 
            )
 
    );
 
    return post_proc_fe;
  };
 
  auto get_bdy_post_proc_fe = [&](auto post_proc_mesh) {
    auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
      return setParentDofs(
          fe, field, op, bit(get_skin_parent_bit()),
 
          [&]() {
            boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
                new BoundaryParentEle(mField));
            return fe_parent;
          },
 
          QUIET, Sev::noisy);
    };
 
    auto post_proc_fe =
        boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
    post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
      return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
          get_start_bit() + nb_levels - 1);
    };
 
    CHKERR setParentDofs(
        post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new BoundaryParentEle(mField));
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
 
    struct OpGetNormal : public BoundaryEleOp {
 
      OpGetNormal(boost::shared_ptr<VectorDouble> l_ptr,
                  boost::shared_ptr<MatrixDouble> n_ptr)
          : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), ptrL(l_ptr),
            ptrNormal(n_ptr) {}
 
      MoFEMErrorCode doWork(int side, EntityType type,
                            EntitiesFieldData::EntData &data) {
        MoFEMFunctionBegin;
        auto t_l = getFTensor0FromVec(*ptrL);
        auto t_n_fe = getFTensor1NormalsAtGaussPts();
        ptrNormal->resize(SPACE_DIM, getGaussPts().size2(), false);
        auto t_n = getFTensor1FromMat<SPACE_DIM>(*ptrNormal);
        for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
          t_n(i) = t_n_fe(i) * t_l / std::sqrt(t_n_fe(i) * t_n_fe(i));
          ++t_n_fe;
          ++t_l;
          ++t_n;
        }
        MoFEMFunctionReturn(0);
      };
 
    protected:
      boost::shared_ptr<VectorDouble> ptrL;
      boost::shared_ptr<MatrixDouble> ptrNormal;
    };
 
    auto u_ptr = boost::make_shared<MatrixDouble>();
    auto p_ptr = boost::make_shared<VectorDouble>();
    auto lambda_ptr = boost::make_shared<VectorDouble>();
    auto normal_l_ptr = boost::make_shared<MatrixDouble>();
 
    CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPCOL, "U");
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
 
    CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPCOL, "L");
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("L", lambda_ptr));
    post_proc_fe->getOpPtrVector().push_back(
        new OpGetNormal(lambda_ptr, normal_l_ptr));
 
    CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPCOL, "P");
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("P", p_ptr));
 
    auto op_ptr = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
    op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
                                EntityType type,
                                EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
      MoFEMFunctionReturn(0);
    };
 
    using OpPPMap = OpPostProcMapInMoab<2, 2>;
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(post_proc_fe->getPostProcMesh(),
                    post_proc_fe->getMapGaussPts(),
 
                    OpPPMap::DataMapVec{{"P", p_ptr}},
 
                    OpPPMap::DataMapMat{{"U", u_ptr}, {"L", normal_l_ptr}},
 
                    OpPPMap::DataMapMat(),
 
                    OpPPMap::DataMapMat()
 
                        )
 
    );
 
    return post_proc_fe;
  };
 
  auto get_lift_fe = [&]() {
    auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
      return setParentDofs(
          fe, field, op, bit(get_skin_parent_bit()),
 
          [&]() {
            boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
                new BoundaryParentEle(mField));
            return fe_parent;
          },
 
          QUIET, Sev::noisy);
    };
 
    auto fe = boost::make_shared<BoundaryEle>(mField);
    fe->exeTestHook = [](FEMethod *fe_ptr) {
      return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
          get_start_bit() + nb_levels - 1);
    };
 
    auto lift_ptr = boost::make_shared<VectorDouble>();
    auto p_ptr = boost::make_shared<VectorDouble>();
    auto ents_ptr = boost::make_shared<Range>();
 
    CHKERR setParentDofs(
        fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
 
        [&]() {
          boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
              new BoundaryParentEle(mField));
          return fe_parent;
        },
 
        QUIET, Sev::noisy);
 
    std::vector<const CubitMeshSets *> vec_ptr;
    CHKERR mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
        std::regex("LIFT"), vec_ptr);
    for (auto m_ptr : vec_ptr) {
      auto meshset = m_ptr->getMeshset();
      Range ents;
      CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM - 1,
                                                         ents, true);
      ents_ptr->merge(ents);
    }
 
    MOFEM_LOG("FS", Sev::noisy) << "Lift ents " << (*ents_ptr);
    CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
    CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "L");
    fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("L", p_ptr));
    fe->getOpPtrVector().push_back(
        new OpCalculateLift("L", p_ptr, lift_ptr, ents_ptr));
 
    return std::make_pair(fe, lift_ptr);
  };
 
  auto set_post_proc_monitor = [&](auto ts) {
    MoFEMFunctionBegin;
    DM dm;
    CHKERR TSGetDM(ts, &dm);
    boost::shared_ptr<FEMethod> null_fe;
    auto post_proc_mesh = boost::make_shared<moab::Core>();
    auto monitor_ptr = boost::make_shared<Monitor>(
        SmartPetscObj<DM>(dm, true), post_proc_mesh,
        get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
        get_lift_fe());
    CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
                               null_fe, monitor_ptr);
    MoFEMFunctionReturn(0);
  };
 
  auto dm = simple->getDM();
  auto ts = createTS(mField.get_comm());
  CHKERR TSSetDM(ts, dm);
 
  auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
  tsPrePostProc = ts_pre_post_proc;
 
  if (auto ptr = tsPrePostProc.lock()) {
 
    ptr->fsRawPtr = this;
    ptr->solverSubDM = create_solver_dm(simple->getDM());
    ptr->globSol = createDMVector(dm);
    CHKERR DMoFEMMeshToLocalVector(dm, ptr->globSol, INSERT_VALUES,
                                   SCATTER_FORWARD);
    CHKERR VecAssemblyBegin(ptr->globSol);
    CHKERR VecAssemblyEnd(ptr->globSol);
 
    auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
 
    CHKERR set_post_proc_monitor(sub_ts);
 
    // Add monitor to time solver
    CHKERR TSSetFromOptions(ts);
    CHKERR ptr->tsSetUp(ts);
    CHKERR TSSetUp(ts);
 
    auto print_fields_in_section = [&]() {
      MoFEMFunctionBegin;
      auto section = mField.getInterface<ISManager>()->sectionCreate(
          simple->getProblemName());
      PetscInt num_fields;
      CHKERR PetscSectionGetNumFields(section, &num_fields);
      for (int f = 0; f < num_fields; ++f) {
        const char *field_name;
        CHKERR PetscSectionGetFieldName(section, f, &field_name);
        MOFEM_LOG("FS", Sev::inform)
            << "Field " << f << " " << std::string(field_name);
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR print_fields_in_section();
 
    CHKERR TSSolve(ts, ptr->globSol);
  }
 
  MoFEMFunctionReturn(0);
}

Friends And Related Symbol Documentation

TSPrePostProc

friend struct TSPrePostProc

friend

Definition at line 964 of file free_surface.cpp.

Member Data Documentation

mField

MoFEM::Interface& FreeSurface::mField

Examples

mofem/tutorials/vec-5_free_surface/free_surface.cpp.

Definition at line 829 of file free_surface.cpp.

---

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

• tutorials/vec-5_free_surface/free_surface.cpp