VEC-8: Arc-Length

Note

Prerequisites of this tutorial include VEC-2: Nonlinear elastic

Note

Intended learning outcome:

• solving nonlinear vector-valued problem in MoFEM
• arc-length control

/**
 * @file nonlinear_elastic.cpp
 * @example nonlinear_elastic.cpp
 * @date 2025-10-5
 *
 * @copyright Anonymous authors (c) 2025 under the MIT license (see LICENSE for
 * details)
 */
 
#include <MoFEM.hpp>
#include <MatrixFunction.hpp>
 
using namespace MoFEM;
 
constexpr int SPACE_DIM =
    EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
 
using EntData = EntitiesFieldData::EntData;
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
using BoundaryEle =
    PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
 
using DomainNaturalBC =
    NaturalBC<DomainEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
using OpBodyForce =
    DomainNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
 
using BoundaryNaturalBC =
    NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
using OpForce = BoundaryNaturalBC::OpFlux<NaturalForceMeshsets, 1, SPACE_DIM>;
 
template <int DIM> struct PostProcEleByDim;
 
template <> struct PostProcEleByDim<2> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBase<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBase<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
};
 
template <> struct PostProcEleByDim<3> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBase<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBase<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
};
 
using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle;
using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy;
using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain;
using PostProcEdges =
    PostProcBrokenMeshInMoabBase<EdgeElementForcesAndSourcesCore>;
 
#include <HenckyOps.hpp>
using namespace HenckyOps;
#include <Monitor.hpp>
 
struct ArcTimeScale : public TimeScale {
  double getScale(const double time) override {
    if (snesPtr) {
      double lambda;
      CHK_THROW_MESSAGE(SNESNewtonALGetLoadParameter(snesPtr, &lambda),
                        "Cannot get load parameter from SNES");
      return lambda;
    } else {
      return time;
    }
  }
  static inline SmartPetscObj<SNES> snesPtr = nullptr;
};
 
using ExampleTimeScale = ArcTimeScale;
#include <ExampleArcLength.hpp>
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
  // Initialisation of MoFEM/PETSc and MOAB data structures
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  // Add logging channel for example
  auto core_log = logging::core::get();
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
  LogManager::setLog("EXAMPLE");
  MOFEM_LOG_TAG("EXAMPLE", "example");
 
  try {
 
    //! [Register MoFEM discrete manager in PETSc]
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    //! [Register MoFEM discrete manager in PETSc
 
    //! [Create MoAB]
    moab::Core mb_instance;              ///< mesh database
    moab::Interface &moab = mb_instance; ///< mesh database interface
    //! [Create MoAB]
 
    //! [Create MoFEM]
    MoFEM::Core core(moab);           ///< finite element database
    MoFEM::Interface &m_field = core; ///< finite element database interface
    //! [Create MoFEM]
 
    //! [Example]
    ExampleArcLength ex(m_field);
    CHKERR ex.runProblem();
    //! [Example]
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
}

/**
 * @file ExampleArcLength.hpp
 * @example ExampleArcLength.hpp
 * @date 2025-10-5
 *
 * @copyright Anonymous authors (c) 2025 under the MIT license (see LICENSE for
 * details)
 * 
 */
 
#include <ExampleNonlinearElastic.hpp>
 
struct ExampleArcLength : public ExampleNonlinearElastic {
 
  ExampleArcLength(MoFEM::Interface &m_field)
      : ExampleNonlinearElastic(m_field) {}
 
  MoFEMErrorCode runProblem();
 
protected:
  MoFEMErrorCode ArcLengthSolve();
  MoFEMErrorCode checkResults();
};
 
//! [Run problem]
MoFEMErrorCode ExampleArcLength::runProblem() {
  MoFEMFunctionBegin;
  
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR ArcLengthSolve();
  CHKERR gettingNorms();
  CHKERR outputResults();
  CHKERR checkResults();
 
  ExampleTimeScale::snesPtr.reset();
 
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Arc Length Solve]
MoFEMErrorCode ExampleArcLength::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",
        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 ExampleArcLength::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]