COR-12: Cohesive Interface with Arc-Length Control

Note

Prerequisites of this tutorial include VEC-0: Linear elasticity and VEC-9: Arc-Length.

Intended learning outcomes:

• defining cohesive interface elements on prism entities
• evolving an interface damage law
• solving a softening problem with arc-length control
• assembling essential and natural boundary conditions

/**
 * \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;
}

/** \file CohesiveInterfaceElement.hpp
  \brief Implementation of linear interface element
 
*/
 
namespace CohesiveElement {
 
/** \brief Cohesive element implementation
 
  \bug Interface element not working with HO geometry.
*/
struct CohesiveInterfaceElement {
 
  struct CommonData {
    MatrixDouble gapGlob;
    MatrixDouble gapLoc;
    ublas::vector<MatrixDouble> R;
  };
  CommonData commonData;
 
  struct MyPrism : public MoFEM::FlatPrismElementForcesAndSourcesCore {
    MyPrism(MoFEM::Interface &m_field)
        : MoFEM::FlatPrismElementForcesAndSourcesCore(m_field) {
      meshPositionsFieldName = "GEOMETRY";
    }
    int getRule(int order) { return 2 * order; };
  };
  MyPrism feRhs;
  MyPrism feLhs;
  MyPrism feHistory;
 
  CohesiveInterfaceElement(MoFEM::Interface &m_field)
      : feRhs(m_field), feLhs(m_field), feHistory(m_field){};
 
  virtual ~CohesiveInterfaceElement() {}
 
  MyPrism &getFeRhs() { return feRhs; }
  MyPrism &getFeLhs() { return feLhs; }
  MyPrism &getFeHistory() { return feHistory; }
 
  /** \brief Constitutive (physical) equation for interface
 
    This is linear degradation model. Material parameters are: strength
    \f$f_t\f$, interface fracture energy \f$G_f\f$, elastic material stiffness
    \f$E\f$. Parameter \f$\beta\f$ controls how interface opening is calculated.
 
    Model parameter is interface penalty thickness \f$h\f$.
 
    */
  struct PhysicalEquation {
 
    MoFEM::Interface &mField;
    bool isInitialised;
    PhysicalEquation(MoFEM::Interface &m_field)
        : mField(m_field), isInitialised(false){};
 
    virtual ~PhysicalEquation() {}
 
    double h, youngModulus, beta, ft, Gf;
    Range pRisms;
    Tag thKappa, thDamagedPrism;
 
    double E0, g0, kappa1;
 
    /** \brief Initialize history variable data
 
    Create tag on the prism/interface to store damage history variable
 
    */
    MoFEMErrorCode iNitailise(const FEMethod *fe_method) {
      MoFEMFunctionBeginHot;
 
      double def_damaged = 0;
      CHKERR mField.get_moab().tag_get_handle(
          "DAMAGED_PRISM", 1, MB_TYPE_INTEGER, thDamagedPrism,
          MB_TAG_CREAT | MB_TAG_SPARSE, &def_damaged);
      const int def_len = 0;
      CHKERR mField.get_moab().tag_get_handle(
          "_KAPPA", def_len, MB_TYPE_DOUBLE, thKappa,
          MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN, nullptr);
      E0 = youngModulus / h;
      g0 = ft / E0;
      kappa1 = 2 * Gf / ft;
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Calculate gap opening
 
    \f[
    g = \sqrt{ g_n^2 + \beta(g_{s1}^2 + g_{s2}^2)}
    \f]
 
    */
    double calcG(int gg, MatrixDouble gap_loc) {
      return sqrt(pow(gap_loc(gg, 0), 2) +
                  beta * (pow(gap_loc(gg, 1), 2) + pow(gap_loc(gg, 2), 2)));
    }
 
    double *kappaPtr;
    int kappaSize;
 
    /** \brief Get pointer from the mesh to histoy variables \f$\kappa\f$
    */
    MoFEMErrorCode getKappa(int nb_gauss_pts, const FEMethod *fe_method) {
      MoFEMFunctionBeginHot;
      EntityHandle ent = fe_method->numeredEntFiniteElementPtr->getEnt();
 
      rval = mField.get_moab().tag_get_by_ptr(
          thKappa, &ent, 1, (const void **)&kappaPtr, &kappaSize);
      if (rval != MB_SUCCESS || kappaSize != nb_gauss_pts) {
        VectorDouble kappa;
        kappa.resize(nb_gauss_pts);
        kappa.clear();
        int tag_size[1];
        tag_size[0] = nb_gauss_pts;
        void const *tag_data[] = {&kappa[0]};
        CHKERR mField.get_moab().tag_set_by_ptr(thKappa, &ent, 1, tag_data,
                                                tag_size);
        CHKERR mField.get_moab().tag_get_by_ptr(
            thKappa, &ent, 1, (const void **)&kappaPtr, &kappaSize);
      }
      MoFEMFunctionReturnHot(0);
    }
 
    MatrixDouble Dglob, Dloc;
 
    /** \brief Calculate stiffness material matrix
 
    \f[
    \mathbf{D}_\textrm{loc} = (1-\Omega) \mathbf{I} E_0
    \f]
    where \f$E_0\f$ is initial interface penalty stiffness
 
    \f[
    \mathbf{D}_\textrm{glob} = \mathbf{R}^\textrm{T} \mathbf{D}_\textrm{loc}\mathbf{R}
    \f]
 
    */
    MoFEMErrorCode calcDglob(const double omega, MatrixDouble &R) {
      MoFEMFunctionBeginHot;
      Dglob.resize(3, 3);
      Dloc.resize(3, 3);
      Dloc.clear();
      double E = (1 - omega) * E0;
      Dloc(0, 0) = E;
      Dloc(1, 1) = E;
      Dloc(2, 2) = E;
      Dglob = prod(Dloc, R);
      Dglob = prod(trans(R), Dglob);
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Calculate damage
 
    \f[
    \Omega = \frac{1}{2} \frac{(2 G_f E_0+f_t^2)\kappa}{(ft+E_0 \kappa)G_f}
    \f]
 
    */
    MoFEMErrorCode calcOmega(const double kappa, double &omega) {
      MoFEMFunctionBeginHot;
      omega = 0;
      if (kappa >= kappa1) {
        omega = 1;
        MoFEMFunctionReturnHot(0);
      } else if (kappa > 0) {
        double a = (2.0 * Gf * E0 + ft * ft) * kappa;
        double b = (ft + E0 * kappa) * Gf;
        omega = 0.5 * a / b;
      }
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Calculate tangent material stiffness
    */
    MoFEMErrorCode calcTangetDglob(const double omega, double g,
                                   const VectorDouble &gap_loc,
                                   MatrixDouble &R) {
      MoFEMFunctionBeginHot;
      Dglob.resize(3, 3);
      Dloc.resize(3, 3);
      double domega =
          0.5 * (2 * Gf * E0 + ft * ft) / ((ft + (g - ft / E0) * E0) * Gf) -
          0.5 * ((g - ft / E0) * (2 * Gf * E0 + ft * ft) * E0) /
              (pow(ft + (g - ft / E0) * E0, 2) * Gf);
      Dloc.resize(3, 3);
      //r0
      Dloc(0, 0) = (1 - omega) * E0 - domega * E0 * gap_loc[0] * gap_loc[0] / g;
      Dloc(0, 1) = -domega * E0 * gap_loc[0] * beta * gap_loc[1] / g;
      Dloc(0, 2) = -domega * E0 * gap_loc[0] * beta * gap_loc[2] / g;
      //r1
      Dloc(1, 0) = -domega * E0 * gap_loc[1] * gap_loc[0] / g;
      Dloc(1, 1) =
          (1 - omega) * E0 - domega * E0 * gap_loc[1] * beta * gap_loc[1] / g;
      Dloc(1, 2) = -domega * E0 * gap_loc[1] * beta * gap_loc[2] / g;
      //r2
      Dloc(2, 0) = -domega * E0 * gap_loc[2] * gap_loc[0] / g;
      Dloc(2, 1) = -domega * E0 * gap_loc[2] * beta * gap_loc[1] / g;
      Dloc(2, 2) =
          (1 - omega) * E0 - domega * E0 * gap_loc[2] * beta * gap_loc[2] / g;
      Dglob = prod(Dloc, R);
      Dglob = prod(trans(R), Dglob);
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Calculate tractions
 
    \f[
    \mathbf{t} = \mathbf{D}_\textrm{glob}\mathbf{g}
    \f]
 
    */
    virtual MoFEMErrorCode calculateTraction(VectorDouble &traction, int gg,
                                             CommonData &common_data,
                                             const FEMethod *fe_method) {
      MoFEMFunctionBeginHot;
 
      if (!isInitialised) {
        CHKERR iNitailise(fe_method);
        isInitialised = true;
      }
      if (gg == 0) {
        CHKERR getKappa(common_data.gapGlob.size1(), fe_method);
      }
      double g = calcG(gg, common_data.gapLoc);
      double kappa = fmax(g - g0, kappaPtr[gg]);
      double omega = 0;
      CHKERR calcOmega(kappa, omega);
      CHKERR calcDglob(omega, common_data.R[gg]);
      traction.resize(3);
      ublas::matrix_row<MatrixDouble> gap_glob(common_data.gapGlob, gg);
      noalias(traction) = prod(Dglob, gap_glob);
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Calculate tangent stiffness
    */
    virtual MoFEMErrorCode
    calculateTangentStiffeness(MatrixDouble &tangent_matrix, int gg,
                               CommonData &common_data,
                               const FEMethod *fe_method) {
      MoFEMFunctionBeginHot;
 
      try {
        if (!isInitialised) {
          CHKERR iNitailise(fe_method);
          isInitialised = true;
        }
        if (gg == 0) {
          CHKERR getKappa(common_data.gapGlob.size1(), fe_method);
        }
        double g = calcG(gg, common_data.gapLoc);
        double kappa = fmax(g - g0, kappaPtr[gg]);
        double omega = 0;
        CHKERR calcOmega(kappa, omega);
        int iter;
        CHKERR SNESGetIterationNumber(fe_method->snes, &iter);
        if ((kappa <= kappaPtr[gg]) || (kappa >= kappa1) || (iter <= 1)) {
          CHKERR calcDglob(omega, common_data.R[gg]);
        } else {
          ublas::matrix_row<MatrixDouble> g_loc(common_data.gapLoc, gg);
          CHKERR calcTangetDglob(omega, g, g_loc, common_data.R[gg]);
        }
        tangent_matrix.resize(3, 3);
        noalias(tangent_matrix) = Dglob;
      } catch (const std::exception &ex) {
        std::ostringstream ss;
        ss << "throw in method: " << ex.what() << std::endl;
        SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
      }
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Update history variables when converged
    */
    virtual MoFEMErrorCode updateHistory(CommonData &common_data,
                                         const FEMethod *fe_method) {
      MoFEMFunctionBeginHot;
 
      if (!isInitialised) {
        CHKERR iNitailise(fe_method);
        isInitialised = true;
      }
      CHKERR getKappa(common_data.gapGlob.size1(), fe_method);
      bool all_gauss_pts_damaged = true;
      for (unsigned int gg = 0; gg < common_data.gapGlob.size1(); gg++) {
        double omega = 0;
        double g = calcG(gg, common_data.gapLoc);
        double kappa = fmax(g - g0, kappaPtr[gg]);
        kappaPtr[gg] = kappa;
        CHKERR calcOmega(kappa, omega);
        //if(omega < 1.) {
        all_gauss_pts_damaged = false;
        //}
      }
      if (all_gauss_pts_damaged) {
        EntityHandle ent = fe_method->numeredEntFiniteElementPtr->getEnt();
        int set_prism_as_demaged = 1;
        CHKERR mField.get_moab().tag_set_data(thDamagedPrism, &ent, 1,
                                              &set_prism_as_demaged);
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Set negative sign to shape functions on face 4
    */
  struct OpSetSignToShapeFunctions
      : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    OpSetSignToShapeFunctions(const std::string field_name)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROW) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBeginHot;
      if (data.getN().size1() == 0)
        MoFEMFunctionReturnHot(0);
      if (data.getN().size2() == 0)
        MoFEMFunctionReturnHot(0);
      switch (type) {
      case MBVERTEX:
        for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
          for (int nn = 3; nn < 6; nn++) {
            data.getN()(gg, nn) *= -1;
          }
        }
        break;
      case MBEDGE:
        if (side < 3)
          MoFEMFunctionReturnHot(0);
        data.getN() *= -1;
        break;
      case MBTRI:
        if (side == 3)
          MoFEMFunctionReturnHot(0);
        data.getN() *= -1;
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, 1, "data inconsitency");
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Operator calculate gap, normal vector and rotation matrix
  */
  struct OpCalculateGapGlobal
      : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    OpCalculateGapGlobal(const std::string field_name, CommonData &common_data)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
          commonData(common_data) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBeginHot;
      try {
        int nb_dofs = data.getIndices().size();
        if (nb_dofs == 0)
          MoFEMFunctionReturnHot(0);
        int nb_gauss_pts = data.getN().size1();
        if (type == MBVERTEX) {
          commonData.R.resize(nb_gauss_pts);
          for (int gg = 0; gg < nb_gauss_pts; gg++) {
            commonData.R[gg].resize(3, 3);
            double nrm2_normal = 0;
            double nrm2_tangent1 = 0;
            double nrm2_tangent2 = 0;
            for (int dd = 0; dd < 3; dd++) {
              nrm2_normal += pow(getNormalsAtGaussPtsF3()(gg, dd), 2);
              nrm2_tangent1 += pow(getTangent1AtGaussPtF3()(gg, dd), 2);
              nrm2_tangent2 += pow(getTangent2AtGaussPtF3()(gg, dd), 2);
            }
            nrm2_normal = sqrt(nrm2_normal);
            nrm2_tangent1 = sqrt(nrm2_tangent1);
            nrm2_tangent2 = sqrt(nrm2_tangent2);
            for (int dd = 0; dd < 3; dd++) {
              commonData.R[gg](0, dd) =
                  getNormalsAtGaussPtsF3()(gg, dd) / nrm2_normal;
              commonData.R[gg](1, dd) =
                  getTangent1AtGaussPtF3()(gg, dd) / nrm2_tangent1;
              commonData.R[gg](2, dd) =
                  getTangent2AtGaussPtF3()(gg, dd) / nrm2_tangent2;
            }
          }
        }
        if (type == MBVERTEX) {
          commonData.gapGlob.resize(nb_gauss_pts, 3);
          commonData.gapGlob.clear();
        }
        for (int gg = 0; gg < nb_gauss_pts; gg++) {
          for (int dd = 0; dd < 3; dd++) {
            commonData.gapGlob(gg, dd) += cblas_ddot(
                nb_dofs / 3, &data.getN(gg)[0], 1, &data.getFieldData()[dd], 3);
          }
        }
      } catch (const std::exception &ex) {
        std::ostringstream ss;
        ss << "throw in method: " << ex.what() << std::endl;
        SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Operator calculate gap in local coordinate system
  */
  struct OpCalculateGapLocal
      : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    OpCalculateGapLocal(const std::string field_name, CommonData &common_data)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
          commonData(common_data) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBeginHot;
      try {
        if (type == MBVERTEX) {
          int nb_gauss_pts = data.getN().size1();
          commonData.gapLoc.resize(nb_gauss_pts, 3);
          for (int gg = 0; gg < nb_gauss_pts; gg++) {
            ublas::matrix_row<MatrixDouble> gap_glob(commonData.gapGlob, gg);
            ublas::matrix_row<MatrixDouble> gap_loc(commonData.gapLoc, gg);
            gap_loc = prod(commonData.R[gg], gap_glob);
          }
        }
      } catch (const std::exception &ex) {
        std::ostringstream ss;
        ss << "throw in method: " << ex.what() << std::endl;
        SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Operator calculate right hand side vector
  */
  struct OpRhs : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    PhysicalEquation &physicalEqations;
    OpRhs(const std::string field_name, CommonData &common_data,
          PhysicalEquation &physical_eqations)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
          commonData(common_data), physicalEqations(physical_eqations) {}
 
    VectorDouble traction, Nf;
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBeginHot;
 
      try {
        int nb_dofs = data.getIndices().size();
        if (nb_dofs == 0)
          MoFEMFunctionReturnHot(0);
        if (physicalEqations.pRisms.find(
                getNumeredEntFiniteElementPtr()->getEnt()) ==
            physicalEqations.pRisms.end()) {
          MoFEMFunctionReturnHot(0);
        }
        Nf.resize(nb_dofs);
        Nf.clear();
        int nb_gauss_pts = data.getN().size1();
        for (int gg = 0; gg < nb_gauss_pts; gg++) {
          CHKERR physicalEqations.calculateTraction(traction, gg, commonData,
                                                    getFEMethod());
          double w = getGaussPts()(2, gg) *
                     cblas_dnrm2(3, &getNormalsAtGaussPtsF3()(gg, 0), 1) * 0.5;
          for (int nn = 0; nn < nb_dofs / 3; nn++) {
            for (int dd = 0; dd < 3; dd++) {
              Nf[3 * nn + dd] += w * data.getN(gg)[nn] * traction[dd];
            }
          }
        }
        CHKERR VecSetValues(getFEMethod()->snes_f, data.getIndices().size(),
                            &data.getIndices()[0], &Nf[0], ADD_VALUES);
      } catch (const std::exception &ex) {
        std::ostringstream ss;
        ss << "throw in method: " << ex.what() << std::endl;
        SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Operator calculate element stiffens matrix
  */
  struct OpLhs : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    PhysicalEquation &physicalEqations;
    OpLhs(const std::string field_name, CommonData &common_data,
          PhysicalEquation &physical_eqations)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROWCOL),
          commonData(common_data), physicalEqations(physical_eqations) {
      sYmm = false;
    }
 
    MatrixDouble K, D, ND;
    MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type,
                          EntitiesFieldData::EntData &row_data,
                          EntitiesFieldData::EntData &col_data) {
      MoFEMFunctionBeginHot;
 
      try {
        int nb_row = row_data.getIndices().size();
        if (nb_row == 0)
          MoFEMFunctionReturnHot(0);
        int nb_col = col_data.getIndices().size();
        if (nb_col == 0)
          MoFEMFunctionReturnHot(0);
        if (physicalEqations.pRisms.find(
                getNumeredEntFiniteElementPtr()->getEnt()) ==
            physicalEqations.pRisms.end()) {
          MoFEMFunctionReturnHot(0);
        }
        ND.resize(nb_row, 3);
        K.resize(nb_row, nb_col);
        K.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for (int gg = 0; gg < nb_gauss_pts; gg++) {
          CHKERR physicalEqations.calculateTangentStiffeness(D, gg, commonData,
                                                             getFEMethod());
          double w = getGaussPts()(2, gg) *
                     cblas_dnrm2(3, &getNormalsAtGaussPtsF3()(gg, 0), 1) * 0.5;
          ND.clear();
          for (int nn = 0; nn < nb_row / 3; nn++) {
            for (int dd = 0; dd < 3; dd++) {
              for (int DD = 0; DD < 3; DD++) {
                ND(3 * nn + dd, DD) += row_data.getN(gg)[nn] * D(dd, DD);
              }
            }
          }
          for (int nn = 0; nn < nb_row / 3; nn++) {
            for (int dd = 0; dd < 3; dd++) {
              for (int NN = 0; NN < nb_col / 3; NN++) {
                for (int DD = 0; DD < 3; DD++) {
                  K(3 * nn + dd, 3 * NN + DD) +=
                      w * ND(3 * nn + dd, DD) * col_data.getN(gg)[NN];
                }
              }
            }
          }
        }
        CHKERR MatSetValues(getFEMethod()->snes_B, nb_row,
                            &row_data.getIndices()[0], nb_col,
                            &col_data.getIndices()[0], &K(0, 0), ADD_VALUES);
      } catch (const std::exception &ex) {
        std::ostringstream ss;
        ss << "throw in method: " << ex.what() << std::endl;
        SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Operator update history variables
  */
  struct OpHistory
      : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    PhysicalEquation &physicalEqations;
    OpHistory(const std::string field_name, CommonData &common_data,
              PhysicalEquation &physical_eqations)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
          commonData(common_data), physicalEqations(physical_eqations) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBeginHot;
 
      if (type != MBVERTEX)
        MoFEMFunctionReturnHot(0);
      if (physicalEqations.pRisms.find(
              getNumeredEntFiniteElementPtr()->getEnt()) ==
          physicalEqations.pRisms.end()) {
        MoFEMFunctionReturnHot(0);
      }
      CHKERR physicalEqations.updateHistory(commonData, getFEMethod());
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Driver function settting all operators needed for interface element
  */
  MoFEMErrorCode
  addOps(const std::string field_name,
         boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>
             &interfaces) {
    MoFEMFunctionBeginHot;
 
    //Rhs
    feRhs.getOpPtrVector().push_back(new OpSetSignToShapeFunctions(field_name));
    feRhs.getOpPtrVector().push_back(
        new OpCalculateGapGlobal(field_name, commonData));
    feRhs.getOpPtrVector().push_back(
        new OpCalculateGapLocal(field_name, commonData));
    //Lhs
    feLhs.getOpPtrVector().push_back(new OpSetSignToShapeFunctions(field_name));
    feLhs.getOpPtrVector().push_back(
        new OpCalculateGapGlobal(field_name, commonData));
    feLhs.getOpPtrVector().push_back(
        new OpCalculateGapLocal(field_name, commonData));
    //History
    feHistory.getOpPtrVector().push_back(
        new OpSetSignToShapeFunctions(field_name));
    feHistory.getOpPtrVector().push_back(
        new OpCalculateGapGlobal(field_name, commonData));
    feHistory.getOpPtrVector().push_back(
        new OpCalculateGapLocal(field_name, commonData));
 
    //add equations/data for physical interfaces
    boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator pit;
    for (pit = interfaces.begin(); pit != interfaces.end(); pit++) {
      feRhs.getOpPtrVector().push_back(new OpRhs(field_name, commonData, *pit));
      feLhs.getOpPtrVector().push_back(new OpLhs(field_name, commonData, *pit));
      feHistory.getOpPtrVector().push_back(
          new OpHistory(field_name, commonData, *pit));
    }
 
    MoFEMFunctionReturnHot(0);
  }
};
 
} // namespace CohesiveElement

/** \file InterfaceGapArcLengthControl.hpp
  \brief Implementation of arc-length control for cohesive elements
 
  Arc-length in that version controls gap opening
 
*/
 
namespace CohesiveElement {
 
struct ArcLengthIntElemFEMethod : public FEMethod {
  moab::Interface &mOab;
 
  boost::shared_ptr<ArcLengthCtx> arcPtr;
 
  Vec GhostLambdaInt;
  Range Faces3, Faces4;
  Range Edges3, Edges4;
  Range Nodes3, Nodes4;
 
  Tag thDamagedPrism;
 
  ArcLengthIntElemFEMethod(moab::Interface &moab,
                           boost::shared_ptr<ArcLengthCtx> arcptr)
      : FEMethod(), mOab(moab), arcPtr(arcptr) {
    PetscInt ghosts[1] = {0};
    ParallelComm *pcomm = ParallelComm::get_pcomm(&mOab, MYPCOMM_INDEX);
    if (pcomm->rank() == 0) {
      CHKERRABORT(PETSC_COMM_WORLD, VecCreateGhost(PETSC_COMM_WORLD, 1, 1, 0,
                                                   ghosts, &GhostLambdaInt));
    } else {
      CHKERRABORT(PETSC_COMM_WORLD, VecCreateGhost(PETSC_COMM_WORLD, 0, 1, 1,
                                                   ghosts, &GhostLambdaInt));
    }
    Range prisms;
    CHKERRABORT(PETSC_COMM_WORLD,
                mOab.get_entities_by_type(0, MBPRISM, prisms, false));
    for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
      EntityHandle f3, f4;
      CHKERRABORT(PETSC_COMM_WORLD, mOab.side_element(*pit, 2, 3, f3));
      CHKERRABORT(PETSC_COMM_WORLD, mOab.side_element(*pit, 2, 4, f4));
      Faces3.insert(f3);
      Faces4.insert(f4);
    }
 
    CHKERRABORT(PETSC_COMM_WORLD,
                mOab.get_adjacencies(Faces3, 1, false, Edges3));
    CHKERRABORT(PETSC_COMM_WORLD,
                mOab.get_adjacencies(Faces4, 1, false, Edges4));
    CHKERRABORT(PETSC_COMM_WORLD, mOab.get_connectivity(Faces3, Nodes3, true));
    CHKERRABORT(PETSC_COMM_WORLD, mOab.get_connectivity(Faces4, Nodes4, true));
    //Faces3.insert(Edges3.begin(),Edges3.end());
    Faces3.insert(Nodes3.begin(), Nodes3.end());
    //Faces4.insert(Edges4.begin(),Edges4.end());
    Faces4.insert(Nodes4.begin(), Nodes4.end());
 
    double def_damaged = 0;
    CHKERRABORT(
        PETSC_COMM_WORLD,
        mOab.tag_get_handle("DAMAGED_PRISM", 1, MB_TYPE_INTEGER, thDamagedPrism,
                            MB_TAG_CREAT | MB_TAG_SPARSE, &def_damaged));
  }
  ~ArcLengthIntElemFEMethod() {
    CHKERRABORT(PETSC_COMM_WORLD, VecDestroy(&GhostLambdaInt));
  }
 
  /** \brief remove nodes of prims which are fully damaged
    *
    */
  MoFEMErrorCode remove_damaged_prisms_nodes() {
    MoFEMFunctionBeginHot;
    Range prisms;
    CHKERR mOab.get_entities_by_type(0, MBPRISM, prisms, false);
    std::vector<int> is_prism_damaged(prisms.size());
    CHKERR mOab.tag_get_data(thDamagedPrism, prisms,
                             &*is_prism_damaged.begin());
    Range::iterator pit = prisms.begin();
    std::vector<int>::iterator vit = is_prism_damaged.begin();
    for (; pit != prisms.end(); pit++, vit++) {
      if (*vit > 0) {
        Range nodes;
        CHKERR mOab.get_connectivity(&*pit, 1, nodes, true);
        for (Range::iterator nit = nodes.begin(); nit != nodes.end(); nit++) {
          Faces3.erase(*nit);
          Faces4.erase(*nit);
        }
      }
    }
    MoFEMFunctionReturnHot(0);
  }
 
  double lambda_int;
  MoFEMErrorCode preProcess() {
    MoFEMFunctionBeginHot;
    switch (snes_ctx) {
    case CTX_SNESSETFUNCTION: {
      CHKERR calculate_dx_and_dlambda(snes_x);
      CHKERR calculate_db();
      CHKERR calculate_lambda_int(lambda_int);
    } break;
    default:
      break;
    }
    MoFEMFunctionReturnHot(0);
  }
 
  MoFEMErrorCode calculate_lambda_int(double &_lambda_int_) {
    MoFEMFunctionBeginHot;
    ParallelComm *pcomm = ParallelComm::get_pcomm(&mOab, MYPCOMM_INDEX);
    NumeredDofEntityByLocalIdx::iterator dit, hi_dit;
    dit = problemPtr->getNumeredRowDofsPtr()
              ->get<PetscLocalIdx_mi_tag>()
              .lower_bound(0);
    hi_dit = problemPtr->getNumeredRowDofsPtr()
                 ->get<PetscLocalIdx_mi_tag>()
                 .upper_bound(problemPtr->getNbLocalDofsRow());
    double *array;
    double *array_int_lambda;
    CHKERR VecZeroEntries(GhostLambdaInt);
    CHKERR VecGhostUpdateBegin(GhostLambdaInt, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(GhostLambdaInt, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGetArray(arcPtr->dx, &array);
    CHKERR VecGetArray(GhostLambdaInt, &array_int_lambda);
    array_int_lambda[0] = 0;
    for (; dit != hi_dit; dit++) {
      if (dit->get()->getEntType() != MBVERTEX)
        continue;
      if (pcomm->rank() != dit->get()->getPart())
        continue;
      if (Nodes3.find(dit->get()->getEnt()) != Nodes3.end()) {
        array_int_lambda[0] += array[dit->get()->getPetscLocalDofIdx()];
      }
      if (Nodes4.find(dit->get()->getEnt()) != Nodes4.end()) {
        array_int_lambda[0] -= array[dit->get()->getPetscLocalDofIdx()];
      }
    }
    CHKERR VecRestoreArray(arcPtr->dx, &array);
    CHKERR VecRestoreArray(GhostLambdaInt, &array_int_lambda);
    CHKERR VecGhostUpdateBegin(GhostLambdaInt, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(GhostLambdaInt, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateBegin(GhostLambdaInt, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(GhostLambdaInt, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGetArray(GhostLambdaInt, &array_int_lambda);
    _lambda_int_ = arcPtr->alpha * array_int_lambda[0] +
                   arcPtr->dLambda * arcPtr->beta * sqrt(arcPtr->F_lambda2);
    CHKERR VecRestoreArray(GhostLambdaInt, &array_int_lambda);
    MoFEMFunctionReturnHot(0);
  }
 
  virtual MoFEMErrorCode calculate_db() {
    MoFEMFunctionBeginHot;
    CHKERR VecZeroEntries(arcPtr->db);
    CHKERR VecGhostUpdateBegin(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(arcPtr->db, INSERT_VALUES, SCATTER_FORWARD);
    NumeredDofEntityByLocalIdx::iterator dit, hi_dit;
    dit = problemPtr->getNumeredRowDofsPtr()
              ->get<PetscLocalIdx_mi_tag>()
              .lower_bound(0);
    hi_dit = problemPtr->getNumeredRowDofsPtr()
                 ->get<PetscLocalIdx_mi_tag>()
                 .upper_bound(problemPtr->getNbLocalDofsRow() +
                              problemPtr->getNbGhostDofsRow());
    double *array;
    CHKERR VecGetArray(arcPtr->db, &array);
    for (; dit != hi_dit; dit++) {
      if (dit->get()->getEntType() != MBVERTEX) {
        array[dit->get()->getPetscLocalDofIdx()] = 0;
        continue;
      }
      if (Nodes3.find(dit->get()->getEnt()) != Nodes3.end()) {
        array[dit->get()->getPetscLocalDofIdx()] = +arcPtr->alpha;
      } else if (Nodes4.find(dit->get()->getEnt()) != Nodes4.end()) {
        array[dit->get()->getPetscLocalDofIdx()] = -arcPtr->alpha;
      } else
        array[dit->get()->getPetscLocalDofIdx()] = 0;
    }
    CHKERR VecRestoreArray(arcPtr->db, &array);
    MoFEMFunctionReturnHot(0);
  }
 
  MoFEMErrorCode operator()() {
    MoFEMFunctionBeginHot;
 
    switch (snes_ctx) {
    case CTX_SNESSETFUNCTION: {
      //calculate residual for arc length row
      arcPtr->res_lambda = lambda_int - arcPtr->s;
      CHKERR VecSetValue(snes_f, arcPtr->getPetscGlobalDofIdx(),
                         arcPtr->res_lambda, ADD_VALUES);
      MOFEM_LOG_C(
          "ARC_LENGTH", Sev::inform,
          "\tres_lambda = %6.4e lambda_int = %6.4e s = %6.4e",
          arcPtr->res_lambda, lambda_int, arcPtr->s);
    } break;
    case CTX_SNESSETJACOBIAN: {
      //calculate diagonal therm
      arcPtr->dIag = arcPtr->beta * sqrt(arcPtr->F_lambda2);
      CHKERR MatSetValue(snes_B, arcPtr->getPetscGlobalDofIdx(),
                         arcPtr->getPetscGlobalDofIdx(), 1, ADD_VALUES);
    } break;
    default:
      break;
    }
 
    MoFEMFunctionReturnHot(0);
  }
 
  MoFEMErrorCode postProcess() {
    MoFEMFunctionBeginHot;
    switch (snes_ctx) {
    case CTX_SNESSETJACOBIAN: {
      CHKERR VecGhostUpdateBegin(arcPtr->ghostDiag, INSERT_VALUES,
                                 SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(arcPtr->ghostDiag, INSERT_VALUES,
                               SCATTER_FORWARD);
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "\tdiag = %6.4e", arcPtr->dIag);
      CHKERR MatAssemblyBegin(snes_B, MAT_FLUSH_ASSEMBLY);
      CHKERR MatAssemblyEnd(snes_B, MAT_FLUSH_ASSEMBLY);
    } break;
    default:
      break;
    }
    MoFEMFunctionReturnHot(0);
  }
 
  MoFEMErrorCode calculate_dx_and_dlambda(Vec &x) {
    MoFEMFunctionBeginHot;
    //dx
    CHKERR VecCopy(x, arcPtr->dx);
    CHKERR VecAXPY(arcPtr->dx, -1, arcPtr->x0);
    //if LAMBDA dof is on this partition
    if (arcPtr->getPetscLocalDofIdx() != -1) {
      double *array;
      CHKERR VecGetArray(arcPtr->dx, &array);
      arcPtr->dLambda = array[arcPtr->getPetscLocalDofIdx()];
      array[arcPtr->getPetscLocalDofIdx()] = 0;
      CHKERR VecRestoreArray(arcPtr->dx, &array);
    }
    //brodcast dlambda
    CHKERR VecGhostUpdateBegin(arcPtr->ghosTdLambda, INSERT_VALUES,
                               SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(arcPtr->ghosTdLambda, INSERT_VALUES,
                             SCATTER_FORWARD);
    //calculate dx2 (dot product)
    CHKERR VecDot(arcPtr->dx, arcPtr->dx, &arcPtr->dx2);
    MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                "\tdlambda = %6.4e dx2 = %6.4e", arcPtr->dLambda,
                arcPtr->dx2);
    MoFEMFunctionReturnHot(0);
  }
 
  MoFEMErrorCode calculate_init_dlambda(double *dlambda) {
    MoFEMFunctionBeginHot;
 
    *dlambda = arcPtr->s / (arcPtr->beta * sqrt(arcPtr->F_lambda2));
    MOFEM_LOG_C(
        "ARC_LENGTH", Sev::inform,
        "\tInit dlambda = %6.4e s = %6.4e beta = %6.4e F_lambda2 = %6.4e",
        *dlambda, arcPtr->s, arcPtr->beta, arcPtr->F_lambda2);
    double a = *dlambda;
    if (a - a != 0) {
      std::ostringstream sss;
      sss << "s " << arcPtr->s << " " << arcPtr->beta << " "
          << arcPtr->F_lambda2;
      SETERRQ(PETSC_COMM_SELF, 1, sss.str().c_str());
    }
 
    MoFEMFunctionReturnHot(0);
  }
 
  MoFEMErrorCode set_dlambda_to_x(Vec &x, double dlambda) {
    MoFEMFunctionBeginHot;
 
    if (arcPtr->getPetscLocalDofIdx() != -1) {
      double *array;
      CHKERR VecGetArray(x, &array);
      double lambda_old = array[arcPtr->getPetscLocalDofIdx()];
      if (!(dlambda == dlambda)) {
        std::ostringstream sss;
        sss << "s " << arcPtr->s << " " << arcPtr->beta << " "
            << arcPtr->F_lambda2;
        SETERRQ(PETSC_COMM_SELF, 1, sss.str().c_str());
      }
      array[arcPtr->getPetscLocalDofIdx()] = lambda_old + dlambda;
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                  "\tlambda = %6.4e, %6.4e (%6.4e)", lambda_old,
                  array[arcPtr->getPetscLocalDofIdx()], dlambda);
      CHKERR VecRestoreArray(x, &array);
    }
 
    MoFEMFunctionReturnHot(0);
  }
};
 
} // namespace CohesiveElement