ArcLengthExample.hpp

Go to the documentation of this file.

/**
 * @file ArcLengthExample.hpp
 * @example mofem/tutorials/vec-9_arc_length/src/ArcLengthExample.hpp
 * @date 2025-10-5
 *
 * @copyright Anonymous authors (c) 2025 under the MIT license (see LICENSE for
 * details)
 * 
 */
 
#include <NonlinearElasticExample.hpp>
 
struct ArcLengthExample : public NonlinearElasticExample {
 
  ArcLengthExample(MoFEM::Interface &m_field)
      : NonlinearElasticExample(m_field) {}
 
  MoFEMErrorCode runProblem();
 
protected:
  MoFEMErrorCode ArcLengthSolve();
  MoFEMErrorCode checkResults();
};
 
//! [Run problem]
MoFEMErrorCode ArcLengthExample::runProblem() {
  MoFEMFunctionBegin;
  
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
 
  PetscBool test_op = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_op", &test_op,
                             PETSC_NULLPTR);
 
  if (useAdolcMaterial == PETSC_TRUE) {
    CHKERR assembleSystemAdolc();
    if (test_op == PETSC_TRUE) {
      CHKERR opTest();
    }
  } else {
    CHKERR assembleSystemHencky();
    if (test_op == PETSC_TRUE) {
      CHKERR opTest();
    }
  }
 
  CHKERR ArcLengthSolve();
  CHKERR gettingNorms();
  CHKERR outputResults();
  CHKERR checkResults();
 
  ExampleTimeScale::snesPtr.reset();
 
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Arc Length Solve]
MoFEMErrorCode ArcLengthExample::ArcLengthSolve() {
  MoFEMFunctionBegin;
  auto *simple = mField.getInterface<Simple>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto dm = simple->getDM();
  auto snes = pipeline_mng->createSNES();
  ExampleTimeScale::snesPtr = snes;
 
  PetscBool post_proc_vol;
  PetscBool post_proc_skin;
 
  if constexpr (SPACE_DIM == 2) {
    post_proc_vol = PETSC_TRUE;
    post_proc_skin = PETSC_FALSE;
  } else {
    post_proc_vol = PETSC_FALSE;
    post_proc_skin = PETSC_TRUE;
  }
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_vol",
                             &post_proc_vol, PETSC_NULLPTR);
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin",
                             &post_proc_skin, PETSC_NULLPTR);
 
  // Setup postprocessing
  auto create_post_proc_fe = [dm, this, simple, post_proc_vol,
                              post_proc_skin]() {
    auto post_proc_ele_domain = [dm, this](auto &pip_domain) {
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip_domain, {H1},
                                                            "GEOMETRY");
      auto common_ptr = commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
          mField, pip_domain, "U", "MAT_ELASTIC", Sev::inform);
      return common_ptr;
    };
 
    auto post_proc_map = [this](auto &pip, auto u_ptr, auto common_ptr) {
      MoFEMFunctionBegin;
 
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      auto X_ptr = boost::make_shared<MatrixDouble>();
 
      pip->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
 
      pip->getOpPtrVector().push_back(
 
          new OpPPMap(pip->getPostProcMesh(), pip->getMapGaussPts(), {},
                      {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
                      {{"GRAD", common_ptr->matGradPtr},
                       {"FIRST_PIOLA", common_ptr->getMatFirstPiolaStress()}},
                      {}));
      MoFEMFunctionReturn(0);
    };
 
    auto push_post_proc_bdy = [dm, this, simple, post_proc_skin,
                               &post_proc_ele_domain,
                               &post_proc_map](auto &pip_bdy) {
      if (post_proc_skin == PETSC_FALSE)
        return boost::shared_ptr<PostProcEleBdy>();
 
      auto u_ptr = boost::make_shared<MatrixDouble>();
      pip_bdy->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
      auto op_loop_side =
          new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
      auto common_ptr = post_proc_ele_domain(op_loop_side->getOpPtrVector());
      pip_bdy->getOpPtrVector().push_back(op_loop_side);
 
      CHKERR post_proc_map(pip_bdy, u_ptr, common_ptr);
 
      return pip_bdy;
    };
 
    auto push_post_proc_domain = [dm, this, simple, post_proc_vol,
                                  &post_proc_ele_domain,
                                  &post_proc_map](auto &pip_domain) {
      if (post_proc_vol == PETSC_FALSE)
        return boost::shared_ptr<PostProcEleDomain>();
 
      auto u_ptr = boost::make_shared<MatrixDouble>();
      pip_domain->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
      auto common_ptr = post_proc_ele_domain(pip_domain->getOpPtrVector());
 
      CHKERR post_proc_map(pip_domain, u_ptr, common_ptr);
 
      return pip_domain;
    };
 
    auto post_proc_fe_domain = boost::make_shared<PostProcEleDomain>(mField);
    auto post_proc_fe_bdy = boost::make_shared<PostProcEleBdy>(mField);
 
    return std::make_pair(push_post_proc_domain(post_proc_fe_domain),
                          push_post_proc_bdy(post_proc_fe_bdy));
  };
 
  auto add_arc_length_load_tangent_operator = [&](auto snes) {
    MoFEMFunctionBegin;
    auto bdy_pipeline = boost::make_shared<BoundaryEle>(mField);
    auto domain_pipeline = boost::make_shared<DomainEle>(mField);
 
    auto integration_rule = [](int, int, int approx_order) {
      return 2 * (approx_order - 1);
    };
    bdy_pipeline->getRuleHook = integration_rule;
    domain_pipeline->getRuleHook = integration_rule;
 
    struct TimeOne : TimeScale {
      double getScale(const double) override { return -1.; }
    };
    auto time_one = boost::make_shared<TimeOne>();
 
    // The trick is here, that is assumed that load is applied linearly
    // proportional to the load factor (time in this case). So derivative of the
    // load with respect to the load factor is the same as the load at time
    // = 1.0
 
    CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
        bdy_pipeline->getOpPtrVector(), mField, "U", {time_one}, "FORCE",
        "PRESSURE", Sev::inform);
    //! [Define gravity vector]
    CHKERR DomainNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
        domain_pipeline->getOpPtrVector(), mField, "U", {time_one},
        "BODY_FORCE", Sev::inform);
 
    auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
    snes_ctx_ptr->getLoadTangent().push_back(
        {simple->getBoundaryFEName(), bdy_pipeline});
    snes_ctx_ptr->getLoadTangent().push_back(
        {simple->getDomainFEName(), domain_pipeline});
 
    auto time_scale = boost::make_shared<ExampleTimeScale>();
 
    auto pre_proc_ptr = boost::make_shared<FEMethod>();
    auto get_bc_hook_rhs = [this, pre_proc_ptr, time_scale]() {
      MoFEMFunctionBeginHot;
      CHKERR EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_ptr,
                                                       {time_scale}, false)();
      MoFEMFunctionReturnHot(0);
    };
    pre_proc_ptr->preProcessHook = get_bc_hook_rhs;
    snes_ctx_ptr->getPreProcLoadTangent().push_back(pre_proc_ptr);
    auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
    auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
      MoFEMFunctionBegin;
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
          mField, post_proc_rhs_ptr, 0.)();
      MoFEMFunctionReturn(0);
    };
    post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
    snes_ctx_ptr->getPostProcLoadTangent().push_back(post_proc_rhs_ptr);
 
    CHKERR SNESNewtonALSetFunction(snes, SnesLoadTangent, snes_ctx_ptr.get());
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_extra_finite_elements_to_solver_pipelines = [&]() {
    MoFEMFunctionBegin;
 
    auto pre_proc_ptr = boost::make_shared<FEMethod>();
    auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
    auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
 
    auto time_scale = boost::make_shared<ExampleTimeScale>();
 
    auto get_bc_hook_rhs = [this, pre_proc_ptr, time_scale]() {
      MoFEMFunctionBeginHot;
      CHKERR EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_ptr,
                                                       {time_scale}, false)();
      MoFEMFunctionReturnHot(0);
    };
 
    pre_proc_ptr->preProcessHook = get_bc_hook_rhs;
 
    auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
      MoFEMFunctionBegin;
      CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
          mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
          mField, post_proc_rhs_ptr, 1.)();
      MoFEMFunctionReturn(0);
    };
    auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
      MoFEMFunctionBegin;
      CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
          mField, post_proc_lhs_ptr, 1.)();
      MoFEMFunctionReturn(0);
    };
    post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
    post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs;
 
    // This is low level pushing finite elements (pipelines) to solver
    auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
    snes_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_ptr);
    snes_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_ptr);
    snes_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs_ptr);
    snes_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs_ptr);
    MoFEMFunctionReturn(0);
  };
 
  auto set_log_and_post_hook = [this](auto post_proc_fe, auto post_proc_bdy) {
    MoFEMFunctionBegin;
 
    auto make_vtk = [this, post_proc_fe, post_proc_bdy](auto iter,
                                                        auto cache_ptr) {
      MoFEMFunctionBegin;
      auto simple = mField.getInterface<Simple>();
      auto dm = simple->getDM();
      if (post_proc_fe) {
        CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe, cache_ptr);
        CHKERR post_proc_fe->writeFile(
            "out_step_" + boost::lexical_cast<std::string>(iter) + ".h5m");
      }
      if (post_proc_bdy) {
        CHKERR DMoFEMLoopFiniteElements(dm, "bFE", post_proc_bdy, cache_ptr);
        CHKERR post_proc_bdy->writeFile(
            "out_skin_" + boost::lexical_cast<std::string>(iter) + ".h5m");
      }
      MoFEMFunctionReturn(0);
    };
 
    auto hook = [this, make_vtk](SNES snes, Vec x, Vec F,
                                 boost::shared_ptr<CacheTuple> cache_ptr,
                                 void *ctx) -> MoFEMErrorCode {
      MoFEMFunctionBegin;
 
      double atol, rtol, stol;
      int maxit, maxf;
      CHKERR SNESGetTolerances(snes, &atol, &rtol, &stol, &maxit, &maxf);
 
      double lambda;
      CHKERR SNESNewtonALGetLoadParameter(snes, &lambda);
 
      double nrm;
      CHKERR VecNorm(F, NORM_2, &nrm);
 
      int iter;
      CHKERR SNESGetIterationNumber(snes, &iter);
      SNESConvergedReason reason;
      CHKERR SNESGetConvergedReason(snes, &reason);
      MOFEM_LOG_C("EXAMPLE", Sev::inform,
                  "SNES iteration %d, reason %d nrm2 %g lambda %g", iter,
                  reason, nrm, lambda);
 
      CHKERR make_vtk(iter, cache_ptr);
 
      MoFEMFunctionReturn(0);
    };
 
    auto simple = mField.getInterface<Simple>();
    auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
    snes_ctx_ptr->getRhsHook() = hook;
 
    MoFEMFunctionReturn(0);
  };
 
  // Set arc length solver
  auto D = createDMVector(simple->getDM());
  CHKERR SNESSetType(snes, SNESNEWTONLS);
  CHKERR SNESSetFromOptions(snes);
 
  auto B = createDMMatrix(dm);
  CHKERR SNESSetJacobian(snes, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
  // Add extra finite elements to SNES solver pipelines to resolve essential
  // boundary conditions
  CHKERR add_extra_finite_elements_to_solver_pipelines();
  CHKERR add_arc_length_load_tangent_operator(snes);
 
  auto [post_proc_fe, post_proc_bdy] = create_post_proc_fe();
  CHKERR set_log_and_post_hook(post_proc_fe, post_proc_bdy);
 
  CHKERR SNESSolve(snes, PETSC_NULLPTR, D);
 
  MoFEMFunctionReturn(0);
}
//! [Arc Length Solve]
 
//! [Check]
MoFEMErrorCode ArcLengthExample::checkResults() {
  MoFEMFunctionBegin;
  PetscInt test_nb = 0;
  PetscBool test_flg = PETSC_FALSE;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-test", &test_nb, &test_flg);
 
  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("EXAMPLE", Sev::verbose) << "Regression norm " << nrm2;
    double regression_value = 0;
    switch (test_nb) {
    case 1:
      regression_value = 1.215134e-04;
      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);
}
//! [Check]