ElasticExample.hpp

Go to the documentation of this file.

/**
 * @file ElasticExample.hpp
 *
 * @brief Implementation of elastic example class
 *
 * @copyright Copyright (c) 2025
 *
 */
 
#include <ElasticSpring.hpp>
#include <FluidLevel.hpp>
#include <CalculateTraction.hpp>
#include <NaturalDomainBC.hpp>
#include <NaturalBoundaryBC.hpp>
#include <HookeOps.hpp>
#include <ElasticPostProc.hpp>
 
struct ElasticExample {
 
  ElasticExample(MoFEM::Interface &m_field) : mField(m_field) {}
  virtual ~ElasticExample() =
      default; // virtual destructor for proper cleanup in derived classes
 
  MoFEMErrorCode runProblem();
 
protected:
  using Pipeline = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>;
 
  MoFEM::Interface &mField;
  virtual MoFEMErrorCode readMesh();
  virtual MoFEMErrorCode setupProblem();
  virtual MoFEMErrorCode boundaryCondition();
  virtual MoFEMErrorCode assembleSystem();
  virtual MoFEMErrorCode solveSystem();
  virtual MoFEMErrorCode outputResults();
  virtual MoFEMErrorCode checkResults();
 
  /** @brief Add constitutive operators to the domain assembly pipelines. */
  virtual MoFEMErrorCode addDomainMaterialOps(Pipeline &lhs_pipeline,
                                              Pipeline &rhs_pipeline);
 
  // Auxiliary functions
  virtual MoFEMErrorCode setUpSolver(SmartPetscObj<KSP> solver);
  virtual MoFEMErrorCode kspSetUpAndSolve(SmartPetscObj<KSP> solver);
};
 
//! [Run problem]
MoFEMErrorCode ElasticExample::runProblem() {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR outputResults();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Read mesh]
MoFEMErrorCode ElasticExample::readMesh() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
  // Add meshsets if config file provided
  CHKERR mField.getInterface<MeshsetsManager>()->setMeshsetFromFile();
  MoFEMFunctionReturn(0);
}
//! [Read mesh]
 
//! [Set up problem]
MoFEMErrorCode ElasticExample::setupProblem() {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
  PetscInt choice_base_value = AINSWORTH;
  CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
                              LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
 
  FieldApproximationBase base;
  switch (choice_base_value) {
  case AINSWORTH:
    base = AINSWORTH_LEGENDRE_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set AINSWORTH_LEGENDRE_BASE for displacements";
    break;
  case DEMKOWICZ:
    base = DEMKOWICZ_JACOBI_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set DEMKOWICZ_JACOBI_BASE for displacements";
    break;
  default:
    base = LASTBASE;
    break;
  }
 
  // Add field
  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
  int order = 2;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
 
  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
 
  CHKERR simple->setFieldOrder("U", order);
  CHKERR simple->setFieldOrder("GEOMETRY", 2);
  CHKERR simple->setUp();
 
  auto project_ho_geometry = [&]() {
    Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
    return mField.loop_dofs("GEOMETRY", ent_method);
  };
  CHKERR project_ho_geometry();
 
  MoFEMFunctionReturn(0);
}
//! [Set up problem]
 
//! [Boundary condition]
MoFEMErrorCode ElasticExample::boundaryCondition() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
 
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
                                           "U", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
                                           "U", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
                                           "U", 2, 2);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                           "REMOVE_ALL", "U", 0, 3);
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "U");
 
  MoFEMFunctionReturn(0);
}
//! [Boundary condition]
 
//! [Push operators to pipeline]
MoFEMErrorCode ElasticExample::addDomainMaterialOps(Pipeline &lhs_pipeline,
                                                    Pipeline &rhs_pipeline) {
  MoFEMFunctionBegin;
  // Add LHS operator for elasticity (stiffness matrix)
  CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
      mField, lhs_pipeline, "U", "MAT_ELASTIC", Sev::verbose);
  // Add RHS operator for internal forces
  CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
      mField, rhs_pipeline, "U", "MAT_ELASTIC", Sev::verbose);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ElasticExample::assembleSystem() {
  MoFEMFunctionBegin;
  auto pip = mField.getInterface<PipelineManager>();
 
  //! [Integration rule]
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
 
  auto integration_rule_bc = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
 
  CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
  CHKERR pip->setDomainLhsIntegrationRule(integration_rule);
  CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
  CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
  //! [Integration rule]
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
 
  //! [Push domain stiffness matrix to pipeline]
  CHKERR addDomainMaterialOps(pip->getOpDomainLhsPipeline(),
                              pip->getOpDomainRhsPipeline());
  //! [Push domain stiffness matrix to pipeline]
 
  //! [Push Body forces]
  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
      pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
  //! [Push Body forces]
 
  //! [Push natural boundary conditions]
  // Add force boundary condition
  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
      pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
  // Add case for mix type of BCs
  CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
      pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::verbose);
  //! [Push natural boundary conditions]
  MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
 
MoFEMErrorCode ElasticExample::setUpSolver(SmartPetscObj<KSP>) { return 0; }
 
// ! [KSP set up]
MoFEMErrorCode ElasticExample::kspSetUpAndSolve(SmartPetscObj<KSP> solver) {
  MoFEMFunctionBegin;
 
  MOFEM_LOG_CHANNEL("TIMER");
  MOFEM_LOG_TAG("TIMER", "timer");
 
  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
  CHKERR KSPSetUp(solver);
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
 
  DM dm;
  CHKERR KSPGetDM(solver, &dm);
  auto D = createDMVector(dm);
  auto F = vectorDuplicate(D);
 
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
  CHKERR KSPSolve(solver, F, D);
  CHKERR mField.getInterface<TestingInterface>()->checkSolver(solver,
                                                              "elasticity");
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
 
  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
  MoFEMFunctionReturn(0);
};
// ! [KSP set up]
 
//! [Solve]
MoFEMErrorCode ElasticExample::solveSystem() {
  MoFEMFunctionBegin;
 
  auto pip = mField.getInterface<PipelineManager>();
  auto solver = pip->createKSP();
  CHKERR setUpSolver(solver);
  CHKERR KSPSetFromOptions(solver);
 
  auto set_essential_bc = [this]() {
    MoFEMFunctionBegin;
    // This is low level pushing finite elements (pipelines) to solver
    auto simple = mField.getInterface<Simple>();
    auto dm = simple->getDM();
    auto ksp_ctx_ptr = getDMKspCtx(dm);
 
    auto pre_proc_rhs = boost::make_shared<FEMethod>();
    auto post_proc_rhs = boost::make_shared<FEMethod>();
    auto post_proc_lhs = boost::make_shared<FEMethod>();
 
    auto get_pre_proc_hook = [this, pre_proc_rhs]() {
      return EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_rhs,
                                                       {});
    };
    pre_proc_rhs->preProcessHook = get_pre_proc_hook();
 
    auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
                                                           post_proc_rhs, 1.)();
      MoFEMFunctionReturnHot(0);
    };
 
    auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(mField,
                                                           post_proc_lhs, 1.)();
      MoFEMFunctionReturnHot(0);
    };
 
    post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
    post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
 
    ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
    ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
    ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_essential_bc();
  CHKERR kspSetUpAndSolve(solver);
 
  MoFEMFunctionReturn(0);
}
//! [Solve]
 
//! [Postprocess results]
MoFEMErrorCode ElasticExample::outputResults() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto evaluate_field_at_the_point = [&]() {
    MoFEMFunctionBegin;
 
    int coords_dim = 3;
    std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
    PetscBool do_eval_field = PETSC_FALSE;
    CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
                                    field_eval_coords.data(), &coords_dim,
                                    &do_eval_field);
 
    if (do_eval_field) {
 
      auto vector_field_ptr = boost::make_shared<MatrixDouble>();
      auto field_eval_data =
          mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
 
      CHKERR mField.getInterface<FieldEvaluatorInterface>()
          ->buildTree<SPACE_DIM>(field_eval_data, simple->getDomainFEName());
 
      field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
      auto no_rule = [](int, int, int) { return -1; };
      auto field_eval_fe_ptr = field_eval_data->feMethodPtr;
      field_eval_fe_ptr->getRuleHook = no_rule;
 
      field_eval_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", vector_field_ptr));
 
      CHKERR mField.getInterface<FieldEvaluatorInterface>()
          ->evalFEAtThePoint<SPACE_DIM>(
              field_eval_coords.data(), 1e-12, simple->getProblemName(),
              simple->getDomainFEName(), field_eval_data,
              mField.get_comm_rank(), mField.get_comm_rank(), nullptr, MF_EXIST,
              QUIET);
 
      if (vector_field_ptr->size1()) {
        auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vector_field_ptr);
        if constexpr (SPACE_DIM == 2)
          MOFEM_LOG("FieldEvaluator", Sev::inform)
              << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
        else
          MOFEM_LOG("FieldEvaluator", Sev::inform)
              << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
              << " U_Z: " << t_disp(2);
      }
 
      MOFEM_LOG_SYNCHRONISE(mField.get_comm());
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR evaluate_field_at_the_point();
  CHKERR ElasticPostProc::postProcessElasticResults<SPACE_DIM, DomainEle,
                                                    BoundaryEle, SideEle>(
      mField, simple->getDM(), simple->getDomainFEName(), "out_elastic.h5m");
  MoFEMFunctionReturn(0);
}
//! [Postprocess results]
 
//! [Check]
MoFEMErrorCode ElasticExample::checkResults() {
  MoFEMFunctionBegin;
  MOFEM_LOG_CHANNEL("WORLD");
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>();
  auto testing = mField.getInterface<TestingInterface>();
 
  // Core owns the interfaces retained by these registry hooks.
  auto check_residual_is_zero = [testing, simple, pip, &m_field = mField]() {
    MoFEMFunctionBegin;
    auto integration_rule = [](int, int, int p_data) { return 2 * p_data + 1; };
    CHKERR pip->setEvaluationIntegrationRule(integration_rule);
 
    auto &domain_evaluation_pipeline = pip->getOpEvaluationPipeline();
 
    if (domain_evaluation_pipeline.empty()) {
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          domain_evaluation_pipeline, {H1}, "GEOMETRY");
 
      auto op_loop_boundary = new OpLoopSide<BoundaryEle>(
          m_field, simple->getBoundaryFEName(), SPACE_DIM - 1, Sev::verbose);
      op_loop_boundary->getSideFEPtr()->getRuleHook = integration_rule;
      auto &boundary_evaluation_pipeline = op_loop_boundary->getOpPtrVector();
 
      CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
          boundary_evaluation_pipeline, {}, "GEOMETRY");
 
      // Add RHS operators for internal forces
      CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
          m_field, domain_evaluation_pipeline, "U", "MAT_ELASTIC",
          Sev::verbose);
 
      CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
          domain_evaluation_pipeline, m_field, "U", Sev::verbose);
 
      CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
          boundary_evaluation_pipeline, m_field, "U", -1, Sev::verbose);
      domain_evaluation_pipeline.push_back(op_loop_boundary);
    }
 
    auto dm = simple->getDM();
    auto res = createDMVector(dm);
    CHKERR VecSetDM(res, PETSC_NULLPTR);
 
    pip->getEvaluationFE()->f = res;
 
    CHKERR VecZeroEntries(res);
 
    CHKERR pip->loopFiniteElementsEvaluation();
    pip->getEvaluationFE()->f = PETSC_NULLPTR;
 
    CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(res);
    CHKERR VecAssemblyEnd(res);
 
    auto zero_residual_at_constrains = [&]() {
      MoFEMFunctionBegin;
      auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
      auto get_post_proc_hook_rhs = [&]() {
        MoFEMFunctionBegin;
        CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
            m_field, fe_post_proc_ptr, res)();
        CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
            m_field, fe_post_proc_ptr, 0, res)();
        MoFEMFunctionReturn(0);
      };
      fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
      CHKERR DMoFEMPostProcessFiniteElements(dm, fe_post_proc_ptr.get());
      MoFEMFunctionReturn(0);
    };
 
    CHKERR zero_residual_at_constrains();
 
    double nrm2;
    CHKERR VecNorm(res.get(), NORM_2, &nrm2);
    MOFEM_LOG_CHANNEL("WORLD");
    MOFEM_LOG_C("WORLD", Sev::inform, "residual = %3.4e\n", nrm2);
    CHKERR testing->checkValue<TestType::SCALAR>("residual norm", nrm2, 0.,
                                                 {1e-8, 0.});
    MoFEMFunctionReturn(0);
  };
 
  auto check_point_displacement = [testing, simple,
                                   comm = mField.get_comm()]() {
    MoFEMFunctionBegin;
    auto displacement_ptr = boost::make_shared<MatrixDouble>();
    bool point_found = false;
    CHKERR testing->addFieldEvaluator<SPACE_DIM>(
        simple->getProblemName(), simple->getDomainFEName(), {4.5, 2.5, 0.},
        [displacement_ptr](auto &pipeline) {
          MoFEMFunctionBegin;
          pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
              "U", displacement_ptr));
          MoFEMFunctionReturn(0);
        },
        1e-12, &point_found);
    const int local_found = point_found;
    int found;
    CHKERR MPI_Allreduce(&local_found, &found, 1, MPI_INT, MPI_SUM, comm);
    CHKERR testing->checkError<TestType::POINTWISE>(
        "point found on mesh", found ? 0. : 1., 0., {0., 0.});
    // Only owners contribute errors; every rank participates in the checks.
    double error_x = 0, error_y = 0;
    if (local_found) {
      using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
      auto t_disp =
          MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::get(
              *displacement_ptr, 1)();
      error_x = std::abs(t_disp(0) - 0.46);
      error_y = std::abs(t_disp(1) + 0.03);
    }
    CHKERR testing->checkError<TestType::POINTWISE>("U_X at point", error_x,
                                                    0.46, {1e-8, 0.});
    CHKERR testing->checkError<TestType::POINTWISE>("U_Y at point", error_y,
                                                    0.03, {1e-8, 0.});
    MoFEMFunctionReturn(0);
  };
 
  CHKERR testing->addTestCommandLine("", 1, "residual is zero",
                                     check_residual_is_zero);
  CHKERR testing->addTestCommandLine(
      "", 2, "point displacement",
      [check_residual_is_zero, check_point_displacement]() {
        MoFEMFunctionBegin;
        CHKERR check_residual_is_zero();
        CHKERR check_point_displacement();
        MoFEMFunctionReturn(0);
      });
 
  CHKERR testing->runTests();
 
  MoFEMFunctionReturn(0);
}
//! [Check]