arc_length_interface.cpp

Go to the documentation of this file.

/**
 * \file arc_length_interface.cpp
 * \example mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp
 * \brief Example of arc-length with a cohesive interface element
 
 * \todo Make it work with multi-grid and distributed mesh
 
 * \todo Make more clever step adaptation
 */
 
#include <MoFEM.hpp>
#include <algorithm>
#include <limits>
using namespace MoFEM;
 
constexpr int SPACE_DIM = 3;
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEle =
    PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
using PostProcEle = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
 
struct BoundaryBCs {};
 
#include <ArcLengthTools.hpp>
#include <CohesiveInterfaceElement.hpp>
#include <HookeOps.hpp>
#include <InterfaceGapArcLengthControl.hpp>
#include <CohesiveNaturalBoundaryBC.hpp>
 
using namespace boost::numeric;
 
static char help[] = "\
 -my_file mesh file name\n\
 -my_sr reduction of step size\n\
 -my_its_d desired number of steps\n\
 -my_ms maximal number of steps\n\
 -gamma arc-length step adaptation exponent\n\
 -min_arc_length_step minimum arc-length step\n\
 -max_arc_length_step maximum arc-length step\n\
 -field_eval_coords x,y,z coordinates where displacement and stress are evaluated\n\n";
 
namespace CohesiveElement {
 
struct AssembleRhsVectors : public FEMethod {
 
  boost::shared_ptr<ArcLengthCtx> arcPtr;
 
  AssembleRhsVectors(boost::shared_ptr<ArcLengthCtx> &arc_ptr)
      : arcPtr(arc_ptr) {}
 
  MoFEMErrorCode preProcess() {
    MoFEMFunctionBegin;
 
    switch (snes_ctx) {
    case CTX_SNESNONE: {
    } break;
    case CTX_SNESSETFUNCTION: {
      CHKERR VecZeroEntries(snes_f);
      CHKERR VecGhostUpdateBegin(snes_f, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(snes_f, INSERT_VALUES, SCATTER_FORWARD);
    } break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode postProcess() {
    MoFEMFunctionBegin;
    switch (snes_ctx) {
    case CTX_SNESNONE: {
    } break;
    case CTX_SNESSETFUNCTION: {
      if (snes_f) {
        CHKERR VecGhostUpdateBegin(snes_f, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecGhostUpdateEnd(snes_f, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecAssemblyBegin(snes_f);
        CHKERR VecAssemblyEnd(snes_f);
        // Add the complete load tangent, including Neumann forces.
        CHKERR VecAXPY(snes_f, arcPtr->getFieldData(), arcPtr->F_lambda);
        MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "\tlambda = %6.4e",
                    arcPtr->getFieldData());
        // snes_f norm
        double fnorm;
        CHKERR VecNorm(snes_f, NORM_2, &fnorm);
        MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "\tpre-essential fnorm = %6.4e",
                    fnorm);
      }
    } break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
} // namespace CohesiveElement
 
using namespace CohesiveElement;
 
struct ArcLengthInterfaceExample {
 
  ArcLengthInterfaceExample(moab::Interface &moab) : mOab(moab) {}
 
  MoFEMErrorCode runProblem();
 
protected:
  moab::Interface &mOab;
  boost::shared_ptr<MoFEM::Core> corePtr;
  MoFEM::Interface *mField = nullptr;
  double *stepSizePtr = nullptr;
  int *stepPtr = nullptr;
 
  PetscScalar stepSizeReduction = 1.;
  PetscInt maxSteps = 5;
  PetscInt itsD = 6;
  PetscInt order = 2;
  const std::string naturalBcFeName = "NATURAL_BC";
 
  SmartPetscObj<Vec> f;
  SmartPetscObj<Vec> d;
  Vec dRaw = nullptr;
  SmartPetscObj<Mat> aij;
  SmartPetscObj<Mat> shellAij;
  SmartPetscObj<SNES> snes;
 
  boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>
      interfaceMaterials;
  boost::shared_ptr<ArcLengthCtx> arcCtx;
  boost::shared_ptr<ArcLengthIntElemFEMethod> arcMethodPtr;
  boost::shared_ptr<ArcLengthSnesCtx> snesCtxPtr;
  boost::shared_ptr<AssembleRhsVectors> prePostProcFePtr;
  boost::shared_ptr<CohesiveInterfaceElement> cohesiveElementsPtr;
  boost::shared_ptr<DomainEle> elasticRhsFe;
  boost::shared_ptr<DomainEle> elasticLhsFe;
  boost::shared_ptr<BoundaryEle> naturalBcFe;
  boost::shared_ptr<PostProcEle> postProcPtr;
  boost::scoped_ptr<ArcLengthMatShell> matCtx;
  boost::scoped_ptr<PCArcLengthCtx> pcCtx;
 
  std::array<double, SPACE_DIM> fieldEvalCoords{0., 0., 0.};
  PetscBool doEvalField = PETSC_FALSE;
  boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> fieldEvalData;
  boost::shared_ptr<MatrixDouble> fieldEvalDispPtr;
  boost::shared_ptr<MatrixDouble> fieldEvalStressPtr;
 
  double gamma = 0.5;
  double reduction = 1.;
 
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode boundaryCondition();
  MoFEMErrorCode assembleSystem();
  MoFEMErrorCode solveSystem();
  MoFEMErrorCode outputResults();
  MoFEMErrorCode checkResults();
  MoFEMErrorCode evaluateFieldAtPoint(const int load_step);
};
 
//! [Run problem]
MoFEMErrorCode ArcLengthInterfaceExample::runProblem() {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR outputResults();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Read mesh]
MoFEMErrorCode ArcLengthInterfaceExample::readMesh() {
  MoFEMFunctionBegin;
 
  CHKERR JsonConfigManager::setMeshFileFromJson();
 
  // Read the mesh file before constructing MoFEM::Core so that the MOAB
  // database can be populated first. JsonConfigManager::setMeshFileFromJson()
  // maps the primary JSON mesh to -file_name.
  PetscBool flg = PETSC_TRUE;
  char mesh_file_name[PETSC_MAX_PATH_LEN] = "";
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-my_file",
                               mesh_file_name, sizeof(mesh_file_name), &flg);
  if (flg != PETSC_TRUE) {
    CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-file_name",
                                 mesh_file_name, sizeof(mesh_file_name), &flg);
    if (flg != PETSC_TRUE)
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "*** ERROR -my_file or -file_name (MESH FILE NEEDED)");
  }
 
  // Read mesh to MOAB
  const char *option;
  option = "";
  CHKERR mOab.load_file(mesh_file_name, 0, option);
 
  // Data stored on mesh for restart
  Tag th_step_size, th_step;
  double def_step_size = 1;
  rval = mOab.tag_get_handle("_STEPSIZE", 1, MB_TYPE_DOUBLE, th_step_size,
                             MB_TAG_CREAT | MB_TAG_MESH, &def_step_size);
  if (rval == MB_ALREADY_ALLOCATED)
    rval = MB_SUCCESS;
  CHKERR rval;
  int def_step = 1;
  rval = mOab.tag_get_handle("_STEP", 1, MB_TYPE_INTEGER, th_step,
                             MB_TAG_CREAT | MB_TAG_MESH, &def_step);
  if (rval == MB_ALREADY_ALLOCATED)
    rval = MB_SUCCESS;
  CHKERR rval;
  const void *tag_data_step_size[1];
  EntityHandle root = 0;
  CHKERR mOab.tag_get_by_ptr(th_step_size, &root, 1, tag_data_step_size);
  double &step_size = *(double *)tag_data_step_size[0];
  const void *tag_data_step[1];
  CHKERR mOab.tag_get_by_ptr(th_step, &root, 1, tag_data_step);
  int &step = *(int *)tag_data_step[0];
  // end of data stored for restart
  MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "Start step %d and step_size = %6.4e",
              step, step_size);
 
  stepSizePtr = (double *)tag_data_step_size[0];
  stepPtr = (int *)tag_data_step[0];
 
  // Create MoFEM database after loading the mesh. This preserves loading of
  // JSON options by the Core constructor.
  corePtr = boost::make_shared<MoFEM::Core>(mOab);
  mField = corePtr.get();
 
  MoFEMFunctionReturn(0);
}
//! [Read mesh]
 
//! [Set up problem]
MoFEMErrorCode ArcLengthInterfaceExample::setupProblem() {
  MoFEMFunctionBegin;
  auto &m_field = *mField;
  auto &step = *stepPtr;
  PetscBool flg = PETSC_FALSE;
 
  // Constructing MoFEM::Core loads -json_config and inserts its "petsc" and
  // "mofem" sections into the PETSc options database. Read application
  // controls only after that point so they can be supplied by JSON.
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-my_sr",
                             &stepSizeReduction, &flg);
  if (flg != PETSC_TRUE) {
    stepSizeReduction = 1.;
  }
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_ms", &maxSteps,
                            &flg);
  if (flg != PETSC_TRUE) {
    maxSteps = 5;
  }
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_its_d", &itsD,
                            &flg);
  if (flg != PETSC_TRUE) {
    itsD = 6;
  }
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_order", &order,
                            &flg);
  if (flg != PETSC_TRUE) {
    order = 2;
  }
 
  MeshsetsManager *meshsets_manager_ptr;
  CHKERR m_field.getInterface(meshsets_manager_ptr);
  CHKERR meshsets_manager_ptr->setMeshsetFromFile();
 
  PrismInterface *interface_ptr;
  CHKERR m_field.getInterface(interface_ptr);
 
  Tag th_my_ref_level;
  BitRefLevel def_bit_level = 0;
  CHKERR m_field.get_moab().tag_get_handle(
      "_MY_REFINEMENT_LEVEL", sizeof(BitRefLevel), MB_TYPE_OPAQUE,
      th_my_ref_level, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_BYTES,
      &def_bit_level);
  const EntityHandle root_meshset = m_field.get_moab().get_root_set();
  BitRefLevel *ptr_bit_level0;
  CHKERR m_field.get_moab().tag_get_by_ptr(th_my_ref_level, &root_meshset, 1,
                                           (const void **)&ptr_bit_level0);
  BitRefLevel &bit_level0 = *ptr_bit_level0;
  BitRefLevel problem_bit_level = bit_level0;
 
  if (step == 1) {
 
    // ref meshset ref level 0
    CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
        0, 3, BitRefLevel().set(0));
 
    std::vector<BitRefLevel> bit_levels;
    bit_levels.push_back(BitRefLevel().set(0));
 
    int ll = 1;
 
    auto split_interface = [&](const EntityHandle cubit_meshset,
                               const int meshset_id) -> MoFEMErrorCode {
      MoFEMFunctionBegin;
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "Insert Interface %d", meshset_id);
      {
        // get tet entities form back bit_level
        EntityHandle ref_level_meshset = 0;
        CHKERR mOab.create_meshset(MESHSET_SET, ref_level_meshset);
        CHKERR m_field.getInterface<BitRefManager>()
            ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
                                           BitRefLevel().set(), MBTET,
                                           ref_level_meshset);
        CHKERR m_field.getInterface<BitRefManager>()
            ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
                                           BitRefLevel().set(), MBPRISM,
                                           ref_level_meshset);
        Range ref_level_tets;
        CHKERR mOab.get_entities_by_handle(ref_level_meshset, ref_level_tets,
                                           true);
        // get faces and test to split
        CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(), true,
                                       0);
        // set new bit level
        bit_levels.push_back(BitRefLevel().set(ll++));
        // split faces and
        CHKERR interface_ptr->splitSides(ref_level_meshset, bit_levels.back(),
                                         cubit_meshset, true, true, 0);
        // clean meshsets
        CHKERR mOab.delete_entities(&ref_level_meshset, 1);
      }
      // Update cubit meshsets
      for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
        EntityHandle cubit_meshset = ciit->meshset;
        CHKERR m_field.getInterface<BitRefManager>()
            ->updateMeshsetByEntitiesChildren(cubit_meshset, bit_levels.back(),
                                              cubit_meshset, MBMAXTYPE, true);
      }
      MoFEMFunctionReturn(0);
    };
 
    // Prefer Cubit INTERFACESET side sets. Use MAT_INTERF blocksets only as
    // a fallback for meshes configured through JsonConfigManager, since a
    // legacy mesh can store the same interface under both set types.
    std::set<int> split_interface_ids;
    bool split_from_sideset = false;
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             m_field, SIDESET | INTERFACESET, cit)) {
      if (split_interface_ids.insert(cit->getMeshsetId()).second) {
        CHKERR split_interface(cit->getMeshset(), cit->getMeshsetId());
        split_from_sideset = true;
      }
    }
    if (!split_from_sideset) {
      const std::string block_name = "MAT_INTERF";
      for (auto m : m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
               std::regex((boost::format("%s(.*)") % block_name).str()))) {
        if (split_interface_ids.insert(m->getMeshsetId()).second)
          CHKERR split_interface(m->getMeshset(), m->getMeshsetId());
      }
    }
 
    bit_level0 = bit_levels.back();
    problem_bit_level = bit_level0;
 
    /***/
    // Define problem
 
    // Fields
    CHKERR m_field.add_field("DISPLACEMENT", H1, AINSWORTH_LEGENDRE_BASE, 3);
    CHKERR
    m_field.add_field("GEOMETRY", H1, AINSWORTH_LEGENDRE_BASE, 3);
 
    CHKERR m_field.add_field("LAMBDA", NOFIELD, NOBASE, 1);
 
    // Field for ArcLength
    CHKERR
    m_field.add_field("X0_DISPLACEMENT", H1, AINSWORTH_LEGENDRE_BASE, 3);
 
    // FE
    CHKERR m_field.add_finite_element("ELASTIC");
 
    // Define rows/cols and element data
    CHKERR m_field.modify_finite_element_add_field_row("ELASTIC",
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_col("ELASTIC",
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_data("ELASTIC",
                                                        "DISPLACEMENT");
    CHKERR
    m_field.modify_finite_element_add_field_data("ELASTIC", "GEOMETRY");
    CHKERR m_field.modify_finite_element_add_field_row("ELASTIC", "LAMBDA");
    CHKERR m_field.modify_finite_element_add_field_col("ELASTIC", "LAMBDA");
    // this is for paremtis
    CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "LAMBDA");
 
    // FE Interface
    CHKERR m_field.add_finite_element("INTERFACE");
    CHKERR m_field.modify_finite_element_add_field_row("INTERFACE",
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_col("INTERFACE",
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_data("INTERFACE",
                                                        "DISPLACEMENT");
    CHKERR
    m_field.modify_finite_element_add_field_data("INTERFACE", "GEOMETRY");
 
    // FE ArcLength
    CHKERR m_field.add_finite_element("ARC_LENGTH");
 
    // Define rows/cols and element data
    CHKERR m_field.modify_finite_element_add_field_row("ARC_LENGTH", "LAMBDA");
    CHKERR m_field.modify_finite_element_add_field_col("ARC_LENGTH", "LAMBDA");
 
    // elem data
    CHKERR
    m_field.modify_finite_element_add_field_data("ARC_LENGTH", "LAMBDA");
 
    // define problems
    CHKERR m_field.add_problem("ELASTIC_MECHANICS");
 
    // set finite elements for problem
    CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
                                                     "ELASTIC");
    CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
                                                     "INTERFACE");
    CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
                                                     "ARC_LENGTH");
    // set refinement level for problem
    CHKERR m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",
                                                    problem_bit_level);
 
    /***/
    // Declare problem
 
    // add entities (by tets) to the field
    CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENT");
    CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "X0_DISPLACEMENT");
    CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "GEOMETRY");
 
    // add finite elements entities
    CHKERR m_field.add_ents_to_finite_element_by_bit_ref(
        problem_bit_level, BitRefLevel().set(), "ELASTIC", MBTET);
    CHKERR m_field.add_ents_to_finite_element_by_bit_ref(
        problem_bit_level, BitRefLevel().set(), "INTERFACE", MBPRISM);
 
    // Setting up LAMBDA field and ARC_LENGTH interface
    {
      // Add dummy no-field vertex
      EntityHandle no_field_vertex;
      {
        const double coords[] = {0, 0, 0};
        CHKERR m_field.get_moab().create_vertex(coords, no_field_vertex);
        Range range_no_field_vertex;
        range_no_field_vertex.insert(no_field_vertex);
        CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
            range_no_field_vertex, BitRefLevel().set());
 
        EntityHandle lambda_meshset = m_field.get_field_meshset("LAMBDA");
        CHKERR m_field.get_moab().add_entities(lambda_meshset,
                                               range_no_field_vertex);
      }
      // this entity will carray data for this finite element
      EntityHandle meshset_fe_arc_length;
      {
        CHKERR mOab.create_meshset(MESHSET_SET, meshset_fe_arc_length);
        CHKERR mOab.add_entities(meshset_fe_arc_length, &no_field_vertex, 1);
        CHKERR m_field.getInterface<BitRefManager>()->setBitLevelToMeshset(
            meshset_fe_arc_length, BitRefLevel().set());
      }
      // finally add created meshset to the ARC_LENGTH finite element
      CHKERR m_field.add_ents_to_finite_element_by_MESHSET(
          meshset_fe_arc_length, "ARC_LENGTH", false);
    }
 
    // set app. order
    // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
    // (Mark Ainsworth & Joe Coyle)
    CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
 
    CHKERR m_field.set_field_order(0, MBTET, "X0_DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBTRI, "X0_DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBEDGE, "X0_DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBVERTEX, "X0_DISPLACEMENT", 1);
 
    CHKERR m_field.set_field_order(0, MBTET, "GEOMETRY", 2);
    CHKERR m_field.set_field_order(0, MBTRI, "GEOMETRY", 2);
    CHKERR m_field.set_field_order(0, MBEDGE, "GEOMETRY", 2);
    CHKERR m_field.set_field_order(0, MBVERTEX, "GEOMETRY", 1);
  }
 
  if (!m_field.check_finite_element(naturalBcFeName)) {
    CHKERR m_field.add_finite_element(naturalBcFeName);
    CHKERR m_field.modify_finite_element_add_field_row(naturalBcFeName,
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_col(naturalBcFeName,
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_data(naturalBcFeName,
                                                        "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_data(naturalBcFeName,
                                                        "GEOMETRY");
    CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
                                                     naturalBcFeName);
 
    Range natural_bc_faces;
    auto meshsets_manager = m_field.getInterface<MeshsetsManager>();
    for (auto m : meshsets_manager->getCubitMeshsetPtr(NODESET | FORCESET)) {
      CHKERR mOab.get_entities_by_type(m->getMeshset(), MBTRI, natural_bc_faces,
                                       true);
    }
    for (auto m : meshsets_manager->getCubitMeshsetPtr(SIDESET | PRESSURESET)) {
      CHKERR mOab.get_entities_by_type(m->getMeshset(), MBTRI, natural_bc_faces,
                                       true);
    }
    for (auto m : meshsets_manager->getCubitMeshsetPtr(
             std::regex("(FORCE|PRESSURE)(.*)"))) {
      CHKERR mOab.get_entities_by_type(m->getMeshset(), MBTRI, natural_bc_faces,
                                       true);
    }
    CHKERR m_field.add_ents_to_finite_element_by_type(natural_bc_faces, MBTRI,
                                                      naturalBcFeName);
  }
 
  /****/
  // build database
 
  // build field
  CHKERR m_field.build_fields();
  Projection10NodeCoordsOnField ent_method_material(m_field, "GEOMETRY");
  CHKERR m_field.loop_dofs("GEOMETRY", ent_method_material);
 
  // build finite elemnts
  CHKERR m_field.build_finite_elements();
 
  // build adjacencies
  CHKERR m_field.build_adjacencies(problem_bit_level);
 
  /****/
  ProblemsManager *prb_mng_ptr;
  CHKERR m_field.getInterface(prb_mng_ptr);
  // build problem
  CHKERR prb_mng_ptr->buildProblem("ELASTIC_MECHANICS", true);
  // partition
  CHKERR prb_mng_ptr->partitionProblem("ELASTIC_MECHANICS");
  CHKERR prb_mng_ptr->partitionFiniteElements("ELASTIC_MECHANICS", false, 0,
                                              m_field.get_comm_size());
  // what are ghost nodes, see Petsc Manual
  CHKERR prb_mng_ptr->partitionGhostDofs("ELASTIC_MECHANICS");
 
  MoFEMFunctionReturn(0);
}
//! [Set up problem]
 
//! [Boundary condition]
MoFEMErrorCode ArcLengthInterfaceExample::boundaryCondition() {
  MoFEMFunctionBegin;
  auto &m_field = *mField;
 
  auto bc_mng = m_field.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      "ELASTIC_MECHANICS", "DISPLACEMENT");
 
  // print bcs
  MeshsetsManager *mmanager_ptr;
  CHKERR m_field.getInterface(mmanager_ptr);
  CHKERR mmanager_ptr->printDisplacementSet();
  CHKERR mmanager_ptr->printForceSet();
  // print block sets with materials
  CHKERR mmanager_ptr->printMaterialsSet();
 
  MoFEMFunctionReturn(0);
}
//! [Boundary condition]
 
//! [Assemble system]
MoFEMErrorCode ArcLengthInterfaceExample::assembleSystem() {
  MoFEMFunctionBegin;
  auto &m_field = *mField;
  auto &step = *stepPtr;
  auto &step_size = *stepSizePtr;
 
  // create matrices
  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("ELASTIC_MECHANICS",
                                                            COL, f);
  d = vectorDuplicate(f);
  dRaw = d;
  CHKERR m_field.getInterface<MatrixManager>()
      ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS", aij);
 
  auto meshsets_manager = m_field.getInterface<MeshsetsManager>();
  const std::string interface_block_name = "MAT_INTERF";
  for (auto m : meshsets_manager->getCubitMeshsetPtr(std::regex(
           (boost::format("%s(.*)") % interface_block_name).str()))) {
    MOFEM_LOG("ARC_LENGTH", Sev::inform) << *m;
    Mat_Interf mydata;
    CHKERR m->getAttributeDataStructure(mydata);
    MOFEM_LOG("ARC_LENGTH", Sev::inform) << mydata;
 
    interfaceMaterials.push_back(
        new CohesiveInterfaceElement::PhysicalEquation(m_field));
    interfaceMaterials.back().h = 1;
    interfaceMaterials.back().youngModulus = mydata.data.alpha;
    interfaceMaterials.back().beta = mydata.data.beta;
    interfaceMaterials.back().ft = mydata.data.ft;
    interfaceMaterials.back().Gf = mydata.data.Gf;
 
    EntityHandle meshset = m->getMeshset();
    Range tris;
    CHKERR mOab.get_entities_by_type(meshset, MBTRI, tris, true);
    Range ents3d;
    CHKERR mOab.get_adjacencies(tris, 3, false, ents3d, moab::Interface::UNION);
    interfaceMaterials.back().pRisms = ents3d.subset_by_type(MBPRISM);
  }
 
  arcCtx = boost::make_shared<ArcLengthCtx>(m_field, "ELASTIC_MECHANICS");
  arcMethodPtr =
      boost::make_shared<ArcLengthIntElemFEMethod>(m_field.get_moab(), arcCtx);
  snesCtxPtr = boost::make_shared<ArcLengthSnesCtx>(
      m_field, "ELASTIC_MECHANICS", arcCtx);
  prePostProcFePtr = boost::make_shared<AssembleRhsVectors>(arcCtx);
  auto &snes_ctx = *snesCtxPtr;
  auto &pre_post_proc_fe = *prePostProcFePtr;
 
  auto essential_pre_proc_zero_rhs = boost::make_shared<FEMethod>();
  auto essential_pre_proc_lhs = boost::make_shared<FEMethod>();
  auto essential_post_proc_rhs = boost::make_shared<FEMethod>();
  auto essential_post_proc_lhs = boost::make_shared<FEMethod>();
 
  struct ArcScale : public ScalingMethod {
    boost::shared_ptr<ArcLengthCtx> arcCtx;
 
    ArcScale(boost::shared_ptr<ArcLengthCtx> arc_ctx) : arcCtx(arc_ctx) {}
 
    double getScale(const double time) override {
      return arcCtx->getFieldData();
    }
  };
  auto arc_scale = boost::make_shared<ArcScale>(arcCtx);
 
  struct ArcZero : public ScalingMethod {
    double getScale(const double time) override {
      return 0.0;
    }
  };
  auto arc_zero = boost::make_shared<ArcZero>();
 
  auto get_essential_pre_proc_zero =
      [&m_field, essential_pre_proc_zero_rhs, arc_zero]() {
    return EssentialPreProc<DisplacementCubitBcData>(
        m_field, essential_pre_proc_zero_rhs, {arc_zero});
  };
  essential_pre_proc_zero_rhs->preProcessHook =
      get_essential_pre_proc_zero();
  auto get_essential_pre_proc_lhs =
      [&m_field, essential_pre_proc_lhs, arc_scale]() {
    return EssentialPreProc<DisplacementCubitBcData>(
        m_field, essential_pre_proc_lhs, {arc_scale});
  };
  essential_pre_proc_lhs->preProcessHook = get_essential_pre_proc_lhs();
 
  essential_post_proc_rhs->postProcessHook =
      [&m_field, essential_post_proc_rhs, arc_scale]() {
    MoFEMFunctionBegin;
    CHKERR EssentialPreProc<DisplacementCubitBcData>(
        m_field, essential_post_proc_rhs, {arc_scale})();
    CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
        m_field, essential_post_proc_rhs, 1.)();
    MoFEMFunctionReturn(0);
  };
 
  // ArcScale sets the reference prescribed field to lambda * F_lambda, so
  // the constrained residual is u_D - lambda * F_lambda,D.
  essential_post_proc_lhs->postProcessHook =
      EssentialPostProcLhs<DisplacementCubitBcData>(
          m_field, essential_post_proc_lhs, 1.);
 
  elasticRhsFe = boost::make_shared<DomainEle>(m_field);
  elasticLhsFe = boost::make_shared<DomainEle>(m_field);
  auto elastic_integration_rule = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
  elasticRhsFe->getRuleHook = elastic_integration_rule;
  elasticLhsFe->getRuleHook = elastic_integration_rule;
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      elasticRhsFe->getOpPtrVector(), {H1}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      elasticLhsFe->getOpPtrVector(), {H1}, "GEOMETRY");
  CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
      m_field, elasticRhsFe->getOpPtrVector(), "DISPLACEMENT", "MAT_ELASTIC",
      Sev::verbose, true);
  CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
      m_field, elasticLhsFe->getOpPtrVector(), "DISPLACEMENT", "MAT_ELASTIC",
      Sev::verbose);
 
  cohesiveElementsPtr = boost::make_shared<CohesiveInterfaceElement>(m_field);
  auto &cohesive_elements = *cohesiveElementsPtr;
  CHKERR cohesive_elements.addOps("DISPLACEMENT", interfaceMaterials);
 
  PetscInt M, N;
  CHKERR MatGetSize(aij, &M, &N);
  PetscInt m, n;
  CHKERR MatGetLocalSize(aij, &m, &n);
  matCtx.reset(new ArcLengthMatShell(aij, arcCtx, "ELASTIC_MECHANICS"));
  Mat shell_aij;
  CHKERR MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, (void *)matCtx.get(),
                        &shell_aij);
  shellAij = SmartPetscObj<Mat>(shell_aij);
  CHKERR MatShellSetOperation(shellAij, MATOP_MULT,
                              (void (*)(void))ArcLengthMatMultShellOp);
 
  using BoundaryRhsBCs = NaturalBC<BoundaryEle::UserDataOperator>::Assembly<
      PETSC>::LinearForm<GAUSS>;
  using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
 
  naturalBcFe = boost::make_shared<BoundaryEle>(m_field);
  naturalBcFe->getRuleHook = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
  CHKERR AddHOOps<2, 3, 3>::add(naturalBcFe->getOpPtrVector(), {NOSPACE},
                                "GEOMETRY");
  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
      naturalBcFe->getOpPtrVector(), m_field, "DISPLACEMENT", Sev::inform);
 
  snes = createSNES(PETSC_COMM_WORLD);
  CHKERR SNESSetApplicationContext(snes, &snes_ctx);
  CHKERR SNESSetFunction(snes, f, SnesRhs, &snes_ctx);
  CHKERR SNESSetJacobian(snes, shellAij, aij, SnesMat, &snes_ctx);
  CHKERR SNESSetFromOptions(snes);
 
  KSP ksp;
  CHKERR SNESGetKSP(snes, &ksp);
  PC pc;
  CHKERR KSPGetPC(ksp, &pc);
  pcCtx.reset(new PCArcLengthCtx(shellAij, aij, arcCtx));
  CHKERR PCSetType(pc, PCSHELL);
  CHKERR PCShellSetContext(pc, pcCtx.get());
  CHKERR PCShellSetApply(pc, PCApplyArcLength);
  CHKERR PCShellSetSetUp(pc, PCSetupArcLength);
 
  // Rhs
  SnesCtx::FEMethodsSequence &loops_to_do_Rhs = snes_ctx.getComputeRhs();
  snes_ctx.getPreProcComputeRhs().push_back(essential_pre_proc_zero_rhs);
  snes_ctx.getPreProcComputeRhs().push_back(&pre_post_proc_fe);
  loops_to_do_Rhs.push_back(
      SnesCtx::PairNameFEMethodPtr("INTERFACE", &cohesive_elements.getFeRhs()));
  loops_to_do_Rhs.push_back(
      SnesCtx::PairNameFEMethodPtr("ELASTIC", elasticRhsFe.get()));
  loops_to_do_Rhs.push_back(
      SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", arcMethodPtr.get()));
  // Add F_lambda before essential processing replaces the constrained rows
  // with u_D - lambda * F_lambda,D.
  snes_ctx.getPostProcComputeRhs().push_back(&pre_post_proc_fe);
  snes_ctx.getPostProcComputeRhs().push_back(essential_post_proc_rhs);
 
  // Mat
  SnesCtx::FEMethodsSequence &loops_to_do_Mat = snes_ctx.getSetOperators();
  snes_ctx.getPreProcSetOperators().push_back(essential_pre_proc_lhs);
  loops_to_do_Mat.push_back(
      SnesCtx::PairNameFEMethodPtr("INTERFACE", &cohesive_elements.getFeLhs()));
  loops_to_do_Mat.push_back(
      SnesCtx::PairNameFEMethodPtr("ELASTIC", elasticLhsFe.get()));
  loops_to_do_Mat.push_back(
      SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", arcMethodPtr.get()));
  snes_ctx.getPostProcSetOperators().push_back(essential_post_proc_lhs);
 
  // step = 1;
  if (step == 1) {
    step_size = stepSizeReduction;
  } else {
    reduction = stepSizeReduction;
    step++;
  }
 
  auto f_lambda_on_essential_bc = boost::make_shared<FEMethod>();
  f_lambda_on_essential_bc->preProcessHook =
      EssentialPreProc<DisplacementCubitBcData>(
          m_field, f_lambda_on_essential_bc, {}, false);
  // EssentialPreProc stores the unit prescribed field. Its derivative in
  // u_D - lambda * F_lambda,D has the opposite sign, hence vDiag = -1.
  f_lambda_on_essential_bc->postProcessHook =
      EssentialPostProcRhs<DisplacementCubitBcData>(
          m_field, f_lambda_on_essential_bc, -1., arcCtx->F_lambda);
 
  CHKERR VecZeroEntries(arcCtx->F_lambda);
  CHKERR VecGhostUpdateBegin(arcCtx->F_lambda, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(arcCtx->F_lambda, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR m_field.problem_basic_method_preProcess("ELASTIC_MECHANICS",
                                                 *f_lambda_on_essential_bc);
  elasticRhsFe->ksp_f = arcCtx->F_lambda;
  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
                                      *elasticRhsFe);
  naturalBcFe->ksp_f = arcCtx->F_lambda;
  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", naturalBcFeName,
                                      *naturalBcFe);
  CHKERR m_field.problem_basic_method_postProcess("ELASTIC_MECHANICS",
                                                  *f_lambda_on_essential_bc);
  CHKERR VecGhostUpdateBegin(arcCtx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateEnd(arcCtx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecAssemblyBegin(arcCtx->F_lambda);
  CHKERR VecAssemblyEnd(arcCtx->F_lambda);
 
  // F_lambda2
  CHKERR VecDot(arcCtx->F_lambda, arcCtx->F_lambda, &arcCtx->F_lambda2);
  MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "\tFlambda2 = %6.4e",
              arcCtx->F_lambda2);
 
  if (step > 1) {
    CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
        "ELASTIC_MECHANICS", COL, d, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
        "ELASTIC_MECHANICS", "DISPLACEMENT", "X0_DISPLACEMENT", COL, arcCtx->x0,
        INSERT_VALUES, SCATTER_FORWARD);
    double x0_nrm;
    CHKERR VecNorm(arcCtx->x0, NORM_2, &x0_nrm);
    MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                "\tRead x0_nrm = %6.4e dlambda = %6.4e", x0_nrm,
                arcCtx->dLambda);
    CHKERR arcCtx->setAlphaBeta(1, 0);
  } else {
    CHKERR arcCtx->setS(0);
    CHKERR arcCtx->setAlphaBeta(0, 1);
  }
  CHKERR SnesRhs(snes, d, f, &snes_ctx);
 
  postProcPtr = boost::make_shared<PostProcEle>(m_field);
  auto &post_proc = *postProcPtr;
  auto disp_ptr = boost::make_shared<MatrixDouble>();
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      post_proc.getOpPtrVector(), {H1}, "GEOMETRY");
  auto hooke_common_ptr =
      HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
          m_field, post_proc.getOpPtrVector(), "DISPLACEMENT", "MAT_ELASTIC",
          Sev::verbose);
  post_proc.getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>("DISPLACEMENT", disp_ptr));
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
  post_proc.getOpPtrVector().push_back(
      new OpPPMap(post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), {},
                  {{"DISPLACEMENT", disp_ptr}},
                  {{"DISPLACEMENT_GRAD", hooke_common_ptr->matGradPtr}},
                  {{"STRAIN", hooke_common_ptr->getMatStrain()},
                   {"STRESS", hooke_common_ptr->getMatCauchyStress()}}));
 
  int coords_dim = SPACE_DIM;
  CHKERR PetscOptionsGetRealArray(PETSC_NULLPTR, PETSC_NULLPTR,
                                  "-field_eval_coords", fieldEvalCoords.data(),
                                  &coords_dim, &doEvalField);
  if (doEvalField && coords_dim != SPACE_DIM)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_INVALID_DATA,
            "-field_eval_coords requires exactly three coordinates");
 
  if (doEvalField) {
    fieldEvalData =
        m_field.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
    CHKERR m_field.getInterface<FieldEvaluatorInterface>()
        ->buildTree<SPACE_DIM>(fieldEvalData, "ELASTIC");
    fieldEvalData->setEvalPoints(fieldEvalCoords.data(), 1);
 
    auto field_eval_fe = fieldEvalData->feMethodPtr;
    field_eval_fe->getRuleHook = [](int, int, int) { return -1; };
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        field_eval_fe->getOpPtrVector(), {H1}, "GEOMETRY");
 
    auto field_eval_hooke_common_ptr =
        HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
            m_field, field_eval_fe->getOpPtrVector(), "DISPLACEMENT",
            "MAT_ELASTIC", Sev::verbose);
    fieldEvalStressPtr = field_eval_hooke_common_ptr->getMatCauchyStress();
    fieldEvalDispPtr = boost::make_shared<MatrixDouble>();
    field_eval_fe->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("DISPLACEMENT",
                                                    fieldEvalDispPtr));
  }
 
  MoFEMFunctionReturn(0);
}
//! [Assemble system]
 
MoFEMErrorCode
ArcLengthInterfaceExample::evaluateFieldAtPoint(const int load_step) {
  MoFEMFunctionBegin;
  auto &m_field = *mField;
 
  if (!doEvalField)
    MoFEMFunctionReturnHot(0);
 
  CHKERR m_field.getInterface<FieldEvaluatorInterface>()
      ->evalFEAtThePoint<SPACE_DIM>(
          fieldEvalCoords.data(), 1e-12, "ELASTIC_MECHANICS", "ELASTIC",
          fieldEvalData, m_field.get_comm_rank(), m_field.get_comm_rank(),
          nullptr, MF_EXIST, QUIET);
 
  int point_found = fieldEvalDispPtr->size1() && fieldEvalStressPtr->size1();
  int global_point_found = 0;
  MPI_Allreduce(&point_found, &global_point_found, 1, MPI_INT, MPI_SUM,
                m_field.get_comm());
 
  if (point_found) {
    auto t_disp = getFTensor1FromMat<SPACE_DIM>(*fieldEvalDispPtr);
    auto t_stress = getFTensor2SymmetricFromMat<SPACE_DIM>(*fieldEvalStressPtr);
    MOFEM_LOG_C("SYNC", Sev::inform,
                "FieldEvaluator step %d lambda %6.4e point [%6.4e, %6.4e, "
                "%6.4e]",
                load_step, arcCtx->getFieldData(), fieldEvalCoords[0],
                fieldEvalCoords[1], fieldEvalCoords[2]);
    MOFEM_LOG_C("SYNC", Sev::inform,
                "FieldEvaluator displacement [Ux, Uy, Uz] = [%6.4e, "
                "%6.4e, %6.4e]",
                t_disp(0), t_disp(1), t_disp(2));
    MOFEM_LOG_C("SYNC", Sev::inform,
                "FieldEvaluator Cauchy stress [Sxx, Syy, Szz, Sxy, Syz, Sxz] "
                "= [%6.4e, %6.4e, %6.4e, %6.4e, %6.4e, %6.4e]",
                t_stress(0, 0), t_stress(1, 1), t_stress(2, 2), t_stress(0, 1),
                t_stress(1, 2), t_stress(0, 2));
  } else if (!global_point_found && !m_field.get_comm_rank()) {
    MOFEM_LOG_C("ARC_LENGTH", Sev::warning,
                "FieldEvaluator did not find point [%6.4e, %6.4e, %6.4e] in "
                "the ELASTIC domain",
                fieldEvalCoords[0], fieldEvalCoords[1], fieldEvalCoords[2]);
  }
  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
 
  MoFEMFunctionReturn(0);
}
 
//! [Solve]
MoFEMErrorCode ArcLengthInterfaceExample::solveSystem() {
  MoFEMFunctionBegin;
 
  double min_arc_length_step = std::numeric_limits<double>::epsilon();
  double max_arc_length_step = std::numeric_limits<double>::max();
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-gamma", &gamma,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-min_arc_length_step", &min_arc_length_step,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-max_arc_length_step", &max_arc_length_step,
                             PETSC_NULLPTR);
 
  if (min_arc_length_step <= 0.)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_INVALID_DATA,
            "-min_arc_length_step must be positive");
  if (max_arc_length_step <= 0.)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_INVALID_DATA,
            "-max_arc_length_step must be positive");
  if (min_arc_length_step > max_arc_length_step)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_INVALID_DATA,
            "-min_arc_length_step must not exceed -max_arc_length_step");
 
  auto &m_field = *mField;
  auto &step = *stepPtr;
  auto &step_size = *stepSizePtr;
 
  auto &cohesive_elements = *cohesiveElementsPtr;
  auto &post_proc = *postProcPtr;
 
  bool converged_state = false;
  for (; step < maxSteps; step++) {
 
    if (step == 1) {
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "Load Step %d step_size = %6.4e",
                  step, step_size);
      CHKERR arcCtx->setS(step_size);
      CHKERR arcCtx->setAlphaBeta(0, 1);
      CHKERR VecCopy(d, arcCtx->x0);
      double dlambda;
      CHKERR arcMethodPtr->calculate_init_dlambda(&dlambda);
      CHKERR arcMethodPtr->set_dlambda_to_x(dRaw, dlambda);
    } else if (step == 2) {
      CHKERR arcCtx->setAlphaBeta(1, 0);
      CHKERR arcMethodPtr->calculate_dx_and_dlambda(dRaw);
      CHKERR arcMethodPtr->calculate_lambda_int(step_size);
      CHKERR arcCtx->setS(step_size);
      double dlambda = arcCtx->dLambda;
      double dx_nrm;
      CHKERR VecNorm(arcCtx->dx, NORM_2, &dx_nrm);
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                  "Load Step %d step_size = %6.4e dlambda0 = %6.4e "
                  "dx_nrm = %6.4e dx2 = %6.4e",
                  step, step_size, dlambda, dx_nrm, arcCtx->dx2);
      CHKERR VecCopy(d, arcCtx->x0);
      CHKERR VecAXPY(d, 1., arcCtx->dx);
      CHKERR arcMethodPtr->set_dlambda_to_x(dRaw, dlambda);
    } else {
      CHKERR arcMethodPtr->calculate_dx_and_dlambda(dRaw);
      CHKERR arcMethodPtr->calculate_lambda_int(step_size);
      // step_size0_1/step_size0 = step_stize1/step_size
      // step_size0_1 = step_size0*(step_stize1/step_size)
      step_size *= reduction;
      CHKERR arcCtx->setS(step_size);
      double dlambda = reduction * arcCtx->dLambda;
      CHKERR VecScale(arcCtx->dx, reduction);
      double dx_nrm;
      CHKERR VecNorm(arcCtx->dx, NORM_2, &dx_nrm);
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                  "Load Step %d step_size = %6.4e dlambda0 = %6.4e "
                  "dx_nrm = %6.4e dx2 = %6.4e",
                  step, step_size, dlambda, dx_nrm, arcCtx->dx2);
      CHKERR VecCopy(d, arcCtx->x0);
      CHKERR VecAXPY(d, 1., arcCtx->dx);
      CHKERR arcMethodPtr->set_dlambda_to_x(dRaw, dlambda);
    }
 
    CHKERR SNESSolve(snes, PETSC_NULLPTR, d);
 
    // Distribute displacements on all processors
    CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
        "ELASTIC_MECHANICS", COL, d, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "INTERFACE",
                                        cohesive_elements.getFeHistory(), 0,
                                        m_field.get_comm_size());
    // Remove nodes of damaged prisms
    CHKERR arcMethodPtr->remove_damaged_prisms_nodes();
 
    int its;
    CHKERR SNESGetIterationNumber(snes, &its);
    MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "number of Newton iterations = %d",
                its);
 
    SNESConvergedReason reason;
    CHKERR SNESGetConvergedReason(snes, &reason);
 
    if (reason < 0) {
      CHKERR arcCtx->setAlphaBeta(1, 0);
      reduction = 0.1;
      converged_state = false;
      continue;
    } else {
      if (step > 1 && converged_state) {
        reduction = pow((double)itsD / (double)(its + 1), gamma);
        reduction =
            std::clamp(reduction, min_arc_length_step / std::abs(step_size),
                       max_arc_length_step / std::abs(step_size));
        MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "reduction step_size = %6.4e",
                    reduction);
      }
 
      // Save data on mesh
      CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
          "ELASTIC_MECHANICS", COL, d, INSERT_VALUES, SCATTER_REVERSE);
      CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
          "ELASTIC_MECHANICS", "DISPLACEMENT", "X0_DISPLACEMENT", COL,
          arcCtx->x0, INSERT_VALUES, SCATTER_REVERSE);
      converged_state = true;
    }
    //
    if (reason > 0) {
      CHKERR evaluateFieldAtPoint(step);
    }
 
    if (step % 1 == 0) {
 
      CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
                                          post_proc);
      std::ostringstream ss;
      ss << "out_values_" << step << ".h5m";
      CHKERR post_proc.writeFile(ss.str().c_str());
    }
  }
 
  MoFEMFunctionReturn(0);
}
//! [Solve]
 
//! [Postprocess results]
MoFEMErrorCode ArcLengthInterfaceExample::outputResults() {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}
//! [Postprocess results]
 
//! [Check]
MoFEMErrorCode ArcLengthInterfaceExample::checkResults() {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}
//! [Check]
 
int main(int argc, char *argv[]) {
 
  const string default_options = "-ksp_type fgmres \n"
                                 "-pc_type lu \n"
                                 "-pc_factor_mat_solver_type mumps\n"
                                 "-mat_mumps_icntl_20 0\n"
                                 "-ksp_monitor \n"
                                 "-ksp_atol 1e-10 \n"
                                 "-ksp_rtol 1e-10 \n"
                                 "-snes_monitor \n"
                                 "-snes_type newtonls \n"
                                 "-snes_linesearch_type l2 \n"
                                 "-snes_linesearch_monitor \n"
                                 "-snes_max_it 16 \n"
                                 "-snes_atol 1e-8 \n"
                                 "-snes_rtol 1e-8 \n"
                                 "-snes_converged_reason \n";
 
  string param_file = "param_file.petsc";
  if (!static_cast<bool>(ifstream(param_file))) {
    std::ofstream file(param_file.c_str(), std::ios::ate);
    if (file.is_open()) {
      file << default_options;
      file.close();
    }
  }
  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
 
  auto core_log = logging::core::get();
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "ARC_LENGTH"));
  LogManager::setLog("ARC_LENGTH");
  MOFEM_LOG_TAG("ARC_LENGTH", "ArcLength");
 
  try {
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
 
    ArcLengthInterfaceExample example(moab);
    CHKERR example.runProblem();
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}