mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp

Example of arc-length with a cohesive interface element.

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

/**
 * \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>
using namespace MoFEM;
 
constexpr int SPACE_DIM = 3;
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
 
#include <ArcLengthTools.hpp>
#include <CohesiveInterfaceElement.hpp>
#include <HookeOps.hpp>
#include <InterfaceGapArcLengthControl.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\n";
 
#define DATAFILENAME "load_disp.txt"
 
namespace CohesiveElement {
 
struct ArcLengthElement : public ArcLengthIntElemFEMethod {
  MoFEM::Interface &mField;
  Range postProcNodes;
  ArcLengthElement(MoFEM::Interface &m_field,
                   boost::shared_ptr<ArcLengthCtx> &arc_ptr)
      : ArcLengthIntElemFEMethod(m_field.get_moab(), arc_ptr), mField(m_field) {
 
    for (_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(mField, "LoadPath", cit)) {
      EntityHandle meshset = cit->getMeshset();
      Range nodes;
      CHKERRABORT(PETSC_COMM_WORLD,
                  mOab.get_entities_by_type(meshset, MBVERTEX, nodes, true));
      postProcNodes.merge(nodes);
    }
 
    MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "Nb. PostProcNodes %lu",
                postProcNodes.size());
  };
 
  MoFEMErrorCode postProcessLoadPath() {
    MoFEMFunctionBegin;
    FILE *datafile;
    CHKERR PetscFOpen(PETSC_COMM_SELF, DATAFILENAME, "a+", &datafile);
    const auto bit_number_lambda = mField.get_field_bit_number("LAMBDA");
 
    boost::shared_ptr<NumeredDofEntity_multiIndex> numered_dofs_rows =
        problemPtr->getNumeredRowDofsPtr();
    auto lit = numered_dofs_rows->lower_bound(
        FieldEntity::getLoBitNumberUId(bit_number_lambda));
    auto hi_lit = numered_dofs_rows->upper_bound(
        FieldEntity::getHiBitNumberUId(bit_number_lambda));
 
    if (lit == numered_dofs_rows->end()) {
      fclose(datafile);
      MoFEMFunctionReturnHot(0);
    } else if (std::distance(lit, hi_lit) != 1) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Only one DOF is expected");
    }
 
    Range::iterator nit = postProcNodes.begin();
    for (; nit != postProcNodes.end(); nit++) {
      NumeredDofEntityByEnt::iterator dit, hi_dit;
      dit = numered_dofs_rows->get<Ent_mi_tag>().lower_bound(*nit);
      hi_dit = numered_dofs_rows->get<Ent_mi_tag>().upper_bound(*nit);
      double coords[3];
      CHKERR mOab.get_coords(&*nit, 1, coords);
      for (; dit != hi_dit; dit++) {
        MOFEM_LOG_C(
            "ARC_LENGTH", Sev::inform,
            "%s [ %d ] %6.4e -> %s [ %d ] %6.4e -> %3.4f %3.4f %3.4f",
            lit->get()->getName().c_str(), lit->get()->getDofCoeffIdx(),
            lit->get()->getFieldData(), dit->get()->getName().c_str(),
            dit->get()->getDofCoeffIdx(), dit->get()->getFieldData(), coords[0],
            coords[1], coords[2]);
        CHKERR PetscFPrintf(PETSC_COMM_WORLD, datafile, "%6.4e %6.4e ",
                            dit->get()->getFieldData(),
                            lit->get()->getFieldData());
      }
    }
    CHKERR PetscFPrintf(PETSC_COMM_WORLD, datafile, "\n");
    fclose(datafile);
    MoFEMFunctionReturn(0);
  }
};
 
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 F_lambda
        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, "\tfnorm = %6.4e", fnorm);
      }
    } break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
} // namespace CohesiveElement
 
using namespace CohesiveElement;
 
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 {
 
    CHKERR JsonConfigManager::setMeshFileFromJson();
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    int rank;
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 
    // Reade parameters from line command
    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)");
    }
 
    PetscScalar step_size_reduction;
    CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-my_sr",
                               &step_size_reduction, &flg);
    if (flg != PETSC_TRUE) {
      step_size_reduction = 1.;
    }
 
    PetscInt max_steps;
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_ms",
                              &max_steps, &flg);
    if (flg != PETSC_TRUE) {
      max_steps = 5;
    }
 
    int its_d;
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-my_its_d", &its_d, &flg);
    if (flg != PETSC_TRUE) {
      its_d = 6;
    }
 
    PetscInt order;
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_order", &order,
                              &flg);
    if (flg != PETSC_TRUE) {
      order = 2;
    }
 
    // Check if new start or restart. If new start, delete previous
    // load_disp.txt
    if (std::string(mesh_file_name).find("restart") == std::string::npos) {
      remove(DATAFILENAME);
    }
 
    // 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);
 
    // Create MoFEM 2database
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    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;
 
      for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
               m_field, SIDESET | INTERFACESET, cit)) {
        MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "Insert Interface %d",
                    cit->getMeshsetId());
        EntityHandle cubit_meshset = cit->getMeshset();
        {
          // 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);
        }
      }
 
      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("MESH_NODE_POSITIONS", 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", "MESH_NODE_POSITIONS");
      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", "MESH_NODE_POSITIONS");
 
      // 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, "MESH_NODE_POSITIONS");
 
      // 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, "MESH_NODE_POSITIONS", 2);
      CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
      CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
      CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
    }
 
    const std::string natural_bc_fe_name = "NATURAL_BC";
    if (!m_field.check_finite_element(natural_bc_fe_name)) {
      CHKERR m_field.add_finite_element(natural_bc_fe_name);
      CHKERR m_field.modify_finite_element_add_field_row(natural_bc_fe_name,
                                                         "DISPLACEMENT");
      CHKERR m_field.modify_finite_element_add_field_col(natural_bc_fe_name,
                                                         "DISPLACEMENT");
      CHKERR m_field.modify_finite_element_add_field_data(natural_bc_fe_name,
                                                          "DISPLACEMENT");
      CHKERR m_field.modify_finite_element_add_field_data(
          natural_bc_fe_name, "MESH_NODE_POSITIONS");
      CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
                                                       natural_bc_fe_name);
 
      Range natural_bc_faces;
      for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,
                                                      NODESET | FORCESET, it)) {
        CHKERR moab.get_entities_by_type(it->meshset, MBTRI, natural_bc_faces,
                                         true);
      }
      for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
               m_field, SIDESET | PRESSURESET, it)) {
        CHKERR moab.get_entities_by_type(it->meshset, MBTRI, natural_bc_faces,
                                         true);
      }
      CHKERR m_field.add_ents_to_finite_element_by_type(natural_bc_faces, MBTRI,
                                                        natural_bc_fe_name);
    }
 
    /****/
    // build database
 
    // build field
    CHKERR m_field.build_fields();
    Projection10NodeCoordsOnField ent_method_material(m_field,
                                                      "MESH_NODE_POSITIONS");
    CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", 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");
 
    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();
 
    // create matrices
    SmartPetscObj<Vec> F;
    CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(
        "ELASTIC_MECHANICS", COL, F);
    auto D = vectorDuplicate(F);
    Vec D_raw = D;
    SmartPetscObj<Mat> Aij;
    CHKERR m_field.getInterface<MatrixManager>()
        ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS",
                                                        Aij);
 
    // Assemble F and Aij
    double young_modulus = 1;
    boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>
        interface_materials;
 
    // FIXME this in fact allow only for one type of interface,
    // problem is Young Modulus in interface mayoung_modulusterial
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
      MOFEM_LOG("ARC_LENGTH", Sev::inform) << *it;
 
      // Get block name
      string name = it->getName();
 
      if (name.compare(0, 11, "MAT_ELASTIC") == 0) {
        Mat_Elastic mydata;
        CHKERR it->getAttributeDataStructure(mydata);
        MOFEM_LOG("ARC_LENGTH", Sev::inform) << mydata;
        young_modulus = mydata.data.Young;
      } else if (name.compare(0, 10, "MAT_INTERF") == 0) {
        Mat_Interf mydata;
        CHKERR it->getAttributeDataStructure(mydata);
        MOFEM_LOG("ARC_LENGTH", Sev::inform) << mydata;
 
        interface_materials.push_back(
            new CohesiveInterfaceElement::PhysicalEquation(m_field));
        interface_materials.back().h = 1;
        interface_materials.back().youngModulus = mydata.data.alpha;
        interface_materials.back().beta = mydata.data.beta;
        interface_materials.back().ft = mydata.data.ft;
        interface_materials.back().Gf = mydata.data.Gf;
 
        EntityHandle meshset = it->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);
        interface_materials.back().pRisms = ents3d.subset_by_type(MBPRISM);
      }
    }
 
    { // FIXME
      boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator
          pit = interface_materials.begin();
      for (; pit != interface_materials.end(); pit++) {
        pit->youngModulus = young_modulus;
      }
    }
 
    boost::shared_ptr<ArcLengthCtx> arc_ctx = boost::shared_ptr<ArcLengthCtx>(
        new ArcLengthCtx(m_field, "ELASTIC_MECHANICS"));
    boost::scoped_ptr<ArcLengthElement> my_arc_method_ptr(
        new ArcLengthElement(m_field, arc_ctx));
    ArcLengthSnesCtx snes_ctx(m_field, "ELASTIC_MECHANICS", arc_ctx);
    AssembleRhsVectors pre_post_proc_fe(arc_ctx);
 
    auto essential_pre_proc = boost::make_shared<FEMethod>();
    auto essential_post_proc_rhs = boost::make_shared<FEMethod>();
    auto essential_post_proc_lhs = boost::make_shared<FEMethod>();
 
    auto get_essential_pre_proc = [&m_field, essential_pre_proc]() {
      return EssentialPreProc<DisplacementCubitBcData>(m_field,
                                                       essential_pre_proc, {});
    };
    essential_pre_proc->preProcessHook = get_essential_pre_proc();
 
    essential_post_proc_rhs->postProcessHook =
        EssentialPostProcRhs<DisplacementCubitBcData>(
            m_field, essential_post_proc_rhs, 1.);
    essential_post_proc_lhs->postProcessHook =
        EssentialPostProcLhs<DisplacementCubitBcData>(
            m_field, essential_post_proc_lhs, 1.);
 
    auto elastic_rhs_fe = boost::make_shared<DomainEle>(m_field);
    auto elastic_lhs_fe = boost::make_shared<DomainEle>(m_field);
    auto elastic_integration_rule = [](int, int, int approx_order) {
      return 2 * approx_order + 1;
    };
    elastic_rhs_fe->getRuleHook = elastic_integration_rule;
    elastic_lhs_fe->getRuleHook = elastic_integration_rule;
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        elastic_rhs_fe->getOpPtrVector(), {H1}, "MESH_NODE_POSITIONS");
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        elastic_lhs_fe->getOpPtrVector(), {H1}, "MESH_NODE_POSITIONS");
    CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
        m_field, elastic_rhs_fe->getOpPtrVector(), "DISPLACEMENT",
        "MAT_ELASTIC", Sev::verbose, true);
    CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
        m_field, elastic_lhs_fe->getOpPtrVector(), "DISPLACEMENT",
        "MAT_ELASTIC", Sev::verbose);
 
    CohesiveInterfaceElement cohesive_elements(m_field);
    CHKERR cohesive_elements.addOps("DISPLACEMENT", interface_materials);
 
    PetscInt M, N;
    CHKERR MatGetSize(Aij, &M, &N);
    PetscInt m, n;
    CHKERR MatGetLocalSize(Aij, &m, &n);
    boost::scoped_ptr<ArcLengthMatShell> mat_ctx(
        new ArcLengthMatShell(Aij, arc_ctx, "ELASTIC_MECHANICS"));
    Mat shell_aij;
    CHKERR MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, (void *)mat_ctx.get(),
                          &shell_aij);
    SmartPetscObj<Mat> ShellAij(shell_aij);
    CHKERR MatShellSetOperation(ShellAij, MATOP_MULT,
                                (void (*)(void))ArcLengthMatMultShellOp);
 
    using BoundaryEle = PipelineManager::ElementsAndOpsByDim<3>::BoundaryEle;
    using BoundaryRhsBCs = NaturalBC<BoundaryEle::UserDataOperator>::Assembly<
        PETSC>::LinearForm<GAUSS>;
    using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux<NaturalForceMeshsets, 1, 3>;
 
    auto natural_bc_fe = boost::make_shared<BoundaryEle>(m_field);
    natural_bc_fe->getRuleHook = [](int, int, int approx_order) {
      return 2 * approx_order + 1;
    };
    CHKERR AddHOOps<2, 3, 3>::add(natural_bc_fe->getOpPtrVector(), {NOSPACE},
                                  "MESH_NODE_POSITIONS");
    CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
        natural_bc_fe->getOpPtrVector(), m_field, "DISPLACEMENT", {}, "FORCE",
        "PRESSURE", Sev::inform);
 
    boost::scoped_ptr<PCArcLengthCtx> pc_ctx;
    auto 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);
    pc_ctx.reset(new PCArcLengthCtx(ShellAij, Aij, arc_ctx));
    CHKERR PCSetType(pc, PCSHELL);
    CHKERR PCShellSetContext(pc, pc_ctx.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);
    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", elastic_rhs_fe.get()));
    loops_to_do_Rhs.push_back(
        SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", my_arc_method_ptr.get()));
    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);
    loops_to_do_Mat.push_back(SnesCtx::PairNameFEMethodPtr(
        "INTERFACE", &cohesive_elements.getFeLhs()));
    loops_to_do_Mat.push_back(
        SnesCtx::PairNameFEMethodPtr("ELASTIC", elastic_lhs_fe.get()));
    loops_to_do_Mat.push_back(
        SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", my_arc_method_ptr.get()));
    snes_ctx.getPostProcSetOperators().push_back(essential_post_proc_lhs);
 
    double gamma = 0.5, reduction = 1;
    // step = 1;
    if (step == 1) {
      step_size = step_size_reduction;
    } else {
      reduction = step_size_reduction;
      step++;
    }
 
    CHKERR VecZeroEntries(arc_ctx->F_lambda);
    CHKERR VecGhostUpdateBegin(arc_ctx->F_lambda, INSERT_VALUES,
                               SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(arc_ctx->F_lambda, INSERT_VALUES, SCATTER_FORWARD);
    natural_bc_fe->ksp_f = arc_ctx->F_lambda;
    CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", natural_bc_fe_name,
                                        *natural_bc_fe);
    CHKERR VecGhostUpdateBegin(arc_ctx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(arc_ctx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(arc_ctx->F_lambda);
    CHKERR VecAssemblyEnd(arc_ctx->F_lambda);
 
    auto zero_f_lambda_on_essential_bc = boost::make_shared<FEMethod>();
    zero_f_lambda_on_essential_bc->postProcessHook =
        EssentialPostProcRhs<DisplacementCubitBcData>(
            m_field, zero_f_lambda_on_essential_bc, 0., arc_ctx->F_lambda);
    CHKERR m_field.problem_basic_method_postProcess(
        "ELASTIC_MECHANICS", *zero_f_lambda_on_essential_bc);
 
    // F_lambda2
    CHKERR VecDot(arc_ctx->F_lambda, arc_ctx->F_lambda, &arc_ctx->F_lambda2);
    MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "\tFlambda2 = %6.4e",
                arc_ctx->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,
          arc_ctx->x0, INSERT_VALUES, SCATTER_FORWARD);
      double x0_nrm;
      CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                  "\tRead x0_nrm = %6.4e dlambda = %6.4e", x0_nrm,
                  arc_ctx->dLambda);
      CHKERR arc_ctx->setAlphaBeta(1, 0);
    } else {
      CHKERR arc_ctx->setS(0);
      CHKERR arc_ctx->setAlphaBeta(0, 1);
    }
    CHKERR SnesRhs(snes, D, F, &snes_ctx);
 
    PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(
        m_field);
    auto disp_ptr = boost::make_shared<MatrixDouble>();
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        post_proc.getOpPtrVector(), {H1}, "MESH_NODE_POSITIONS");
    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()}}));
 
    bool converged_state = false;
    for (; step < max_steps; step++) {
 
      if (step == 1) {
        MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                    "Load Step %d step_size = %6.4e", step, step_size);
        CHKERR arc_ctx->setS(step_size);
        CHKERR arc_ctx->setAlphaBeta(0, 1);
        CHKERR VecCopy(D, arc_ctx->x0);
        double dlambda;
        CHKERR my_arc_method_ptr->calculate_init_dlambda(&dlambda);
        CHKERR my_arc_method_ptr->set_dlambda_to_x(D_raw, dlambda);
      } else if (step == 2) {
        CHKERR arc_ctx->setAlphaBeta(1, 0);
        CHKERR my_arc_method_ptr->calculate_dx_and_dlambda(D_raw);
        CHKERR my_arc_method_ptr->calculate_lambda_int(step_size);
        CHKERR arc_ctx->setS(step_size);
        double dlambda = arc_ctx->dLambda;
        double dx_nrm;
        CHKERR VecNorm(arc_ctx->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, arc_ctx->dx2);
        CHKERR VecCopy(D, arc_ctx->x0);
        CHKERR VecAXPY(D, 1., arc_ctx->dx);
        CHKERR my_arc_method_ptr->set_dlambda_to_x(D_raw, dlambda);
      } else {
        CHKERR my_arc_method_ptr->calculate_dx_and_dlambda(D_raw);
        CHKERR my_arc_method_ptr->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 arc_ctx->setS(step_size);
        double dlambda = reduction * arc_ctx->dLambda;
        CHKERR VecScale(arc_ctx->dx, reduction);
        double dx_nrm;
        CHKERR VecNorm(arc_ctx->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, arc_ctx->dx2);
        CHKERR VecCopy(D, arc_ctx->x0);
        CHKERR VecAXPY(D, 1., arc_ctx->dx);
        CHKERR my_arc_method_ptr->set_dlambda_to_x(D_raw, 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 my_arc_method_ptr->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 arc_ctx->setAlphaBeta(1, 0);
        reduction = 0.1;
        converged_state = false;
        continue;
      } else {
        if (step > 1 && converged_state) {
          reduction = pow((double)its_d / (double)(its + 1), gamma);
          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,
            arc_ctx->x0, INSERT_VALUES, SCATTER_REVERSE);
        converged_state = true;
      }
      //
      if (reason > 0) {
        FILE *datafile;
        CHKERR PetscFOpen(PETSC_COMM_SELF, DATAFILENAME, "a+", &datafile);
        CHKERR PetscFPrintf(PETSC_COMM_WORLD, datafile, "%d %d ", reason, its);
        fclose(datafile);
        CHKERR my_arc_method_ptr->postProcessLoadPath();
      }
 
      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());
      }
    }
 
    // Save data on mesh
    CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
        "ELASTIC_MECHANICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}