CPMeshCut.cpp

Go to the documentation of this file.

/** \file CPSolvers.cpp
 */
 
/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
#ifdef WITH_TETGEN
 
#include <tetgen.h>
#ifdef REAL
#undef REAL
#endif
#endif
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include <BasicFiniteElements.hpp>
#include <Mortar.hpp>
#include <NeoHookean.hpp>
#include <Hooke.hpp>
 
#include <CrackFrontElement.hpp>
#include <ComplexConstArea.hpp>
#include <MWLS.hpp>
#include <GriffithForceElement.hpp>
 
#include <VolumeLengthQuality.hpp>
#include <CrackPropagation.hpp>
#include <CPMeshCut.hpp>
 
namespace FractureMechanics {
 
MoFEMErrorCode
CPMeshCut::query_interface(boost::typeindex::type_index type_index,
                           UnknownInterface **iface) const {
  *iface = const_cast<CPMeshCut *>(this);
  return 0;
}
 
CPMeshCut::CPMeshCut(CrackPropagation &cp) : cP(cp) {
  ierr = getInterfacesPtr();
  CHKERRABORT(PETSC_COMM_SELF, ierr);
 
  try {
    MoFEM::LogManager::getLog("CPMeshCutSelf");
  } catch (...) {
    auto core_log = logging::core::get();
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmSelf(), "CPMeshCutSelf"));
    LogManager::setLog("CPMeshCutSelf");
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmWorld(), "CPMeshCutWorld"));
    LogManager::setLog("CPMeshCutWorld");
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmSync(), "CPMeshCutSync"));
    LogManager::setLog("CPMeshCutSync");
    MOFEM_LOG_TAG("CPMeshCutSelf", "CPMeshCut");
    MOFEM_LOG_TAG("CPMeshCutWorld", "CPMeshCut");
    MOFEM_LOG_TAG("CPMeshCutSync", "CPMeshCut");
  }
 
  MOFEM_LOG_FUNCTION();
  MOFEM_LOG("CPMeshCutWorld", Sev::noisy) << "CPMeshCutSelf interface created";
}
 
MoFEMErrorCode CPMeshCut::getInterfacesPtr() {
  MoFEMFunctionBegin;
  CHKERR cP.mField.getInterface(cutMeshPtr);
  CHKERR cP.mField.getInterface(bitRefPtr);
  CHKERR cP.mField.getInterface(meshsetMngPtr);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::getOptions() {
  MoFEMFunctionBegin;
 
  PetscBool shift_flg;
  int nmax_shift = 3;
 
  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Mesh cut options", "none");
 
  edgesBlockSet = 100;
  CHKERR PetscOptionsInt("-edges_block_set", "edges side set", "",
                         edgesBlockSet, &edgesBlockSet, &edgesBlockSetFlg);
  vertexBlockSet = 101;
  CHKERR PetscOptionsInt("-vertex_block_set", "vertex block set", "",
                         vertexBlockSet, &vertexBlockSet, PETSC_NULL);
  //to resolve issues with newer cubit
  vertexNodeSet = 102;
  CHKERR PetscOptionsInt("-vertex_node_set", "vertex node set", "",
                         vertexNodeSet, &vertexNodeSet, PETSC_NULL);
  fractionLevel = 2;
  CHKERR PetscOptionsInt("-fraction_level", "fraction of merges merged", "",
                         fractionLevel, &fractionLevel, PETSC_NULL);
  tolCut = 1e-2;
  CHKERR PetscOptionsScalar("-tol_cut", "get tolerance for cut edges", "",
                            tolCut, &tolCut, PETSC_NULL);
  tolCutClose = 2e-1;
  CHKERR PetscOptionsScalar("-tol_cut_close", "get tolerance for cut edges", "",
                            tolCutClose, &tolCutClose, PETSC_NULL);
  tolTrimClose = 2e-1;
  CHKERR PetscOptionsScalar("-tol_trim_close", "get tolerance for trim edges",
                            "", tolTrimClose, &tolTrimClose, PETSC_NULL);
 
  debugCut = PETSC_FALSE;
  CHKERR PetscOptionsBool("-cut_debug_cut", "debug mesh cutting", "", debugCut,
                          &debugCut, NULL);
 
  removePathologicalFrontTris = PETSC_FALSE;
  CHKERR PetscOptionsBool(
      "-remove_pathological_front_tris",
      "remove triangles which have three nodes on crack front", "",
      removePathologicalFrontTris, &removePathologicalFrontTris, NULL);
 
  MOFEM_LOG_C("CPMeshCutWorld", Sev::inform,
              "### Input parameter: -fraction_level %d", fractionLevel);
  MOFEM_LOG_C("CPMeshCutWorld", Sev::inform,
              "### Input parameter: -tol_cut %6.4e", tolCut);
  MOFEM_LOG_C("CPMeshCutWorld", Sev::inform,
              "### Input parameter: -tol_cut_close %6.4e", tolCutClose);
  MOFEM_LOG_C("CPMeshCutWorld", Sev::inform,
              "### Input parameter: -tol_trim_close %6.4e", tolTrimClose);
 
  PetscBool crack_surf_mesh_flg;
  char crack_surf_mesh_name[255];
  CHKERR PetscOptionsString(
      "-my_cut_surface_file", "file with crack surface for initial cut", "",
      crack_surf_mesh_name, crack_surf_mesh_name, 255, &crack_surf_mesh_flg);
  if (crack_surf_mesh_flg)
    cutSurfMeshName = std::string(crack_surf_mesh_name);
 
  sHift[0] = sHift[1] = sHift[2] = 0;
  CHKERR PetscOptionsRealArray("-crack_shift", "shift surface by vector", "",
                               sHift, &nmax_shift, &shift_flg);
 
  cuttingSurfaceSidesetId = 400;
  CHKERR PetscOptionsInt(
      "-cut_surface_side_set", "sideset id with cutting surface", "",
      cuttingSurfaceSidesetId, &cuttingSurfaceSidesetId, PETSC_NULL);
 
  skinOfTheBodyId = 102;
  CHKERR PetscOptionsInt("-body_skin", "body skin sideset id", "",
                         skinOfTheBodyId, &skinOfTheBodyId, PETSC_NULL);
  crackSurfaceId = 200;
  CHKERR PetscOptionsInt("-crack_side_set", "sideset id with crack surface", "",
                         crackSurfaceId, &crackSurfaceId, PETSC_NULL);
  crackFrontId = 201;
  CHKERR PetscOptionsInt("-front_side_set", "sideset id with crack front", "",
                         crackFrontId, &crackFrontId, PETSC_NULL);
 
  nbRefBeforCut = 0;
  CHKERR PetscOptionsInt("-ref_before_cut", "number of refinements before cut",
                         "", nbRefBeforCut, &nbRefBeforCut, PETSC_NULL);
  MOFEM_LOG_C("CPMeshCutWorld", Sev::inform,
              "### Input parameter: -ref_before_cut %d", nbRefBeforCut);
 
  nbRefBeforTrim = 0;
  PetscBool flg = PETSC_FALSE;
  CHKERR PetscOptionsInt("-ref_before_tim", "number of refinements before trim",
                         "", nbRefBeforTrim, &nbRefBeforTrim, &flg);
  if (flg == PETSC_TRUE) {
    MOFEM_LOG("CPMeshCutWorld", Sev::inform)
        << "### Input parameter: -ref_before_tim is still valid but will be "
           "made obsolete in the next release, please use -ref_before_trim "
           "instead.";
  }
 
  CHKERR PetscOptionsInt("-ref_before_trim",
                         "number of refinements before trim", "",
                         nbRefBeforTrim, &nbRefBeforTrim, PETSC_NULL);
  MOFEM_LOG_C("CPMeshCutWorld", Sev::inform,
              "### Input parameter: -ref_before_trim %d", nbRefBeforTrim);
 
  crackedBodyBlockSetId = -1;
  CHKERR PetscOptionsInt(
      "-my_cracked_body_block_set", "blockset id of the cracked body", "",
      crackedBodyBlockSetId, &crackedBodyBlockSetId, PETSC_NULL);
 
  MOFEM_LOG_C("CPMeshCutWorld", Sev::inform,
              "### Input parameter: -my_cracked_body_block_set %d",
              crackedBodyBlockSetId);
 
  contactPrismsBlockSetId = 10000;
  CHKERR PetscOptionsInt(
      "-my_contact_prisms_block_set", "blockset id of contact prisms", "",
      contactPrismsBlockSetId, &contactPrismsBlockSetId, PETSC_NULL);
 
  ierr = PetscOptionsEnd();
  CHKERRG(ierr);
 
  if (shift_flg && nmax_shift != 3) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "three values expected");
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::setVolume(const BitRefLevel &bit) {
  MoFEMFunctionBegin;
  Range bit_tets;
  CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                 MBTET, bit_tets);
  CHKERR cutMeshPtr->setVolume(bit_tets);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::refineMesh(const Range *vol, const bool make_front,
                                     const int verb, const bool debug) {
  MoFEMFunctionBegin;
 
  if (vol) {
    CHKERR bitRefPtr->setBitRefLevel(
        *vol, BitRefLevel().set(BITREFLEVEL_SIZE - 1), true, verb);
    CHKERR cutMeshPtr->setVolume(*vol);
  } else
    // Set cutted front
    CHKERR setVolume(BitRefLevel().set(BITREFLEVEL_SIZE - 1));
 
  CHKERR setMeshOrgCoords();
 
  if (make_front)
    // make crack front from mesh skin
    CHKERR cutMeshPtr->makeFront(true);
 
  auto save_entities = [&](const std::string name, Range ents) {
    MoFEMFunctionBeginHot;
    if (!ents.empty()) {
      EntityHandle meshset;
      CHKERR cP.mField.get_moab().create_meshset(MESHSET_SET, meshset);
      CHKERR cP.mField.get_moab().add_entities(meshset, ents);
      CHKERR cP.mField.get_moab().write_file(name.c_str(), "VTK", "", &meshset,
                                             1);
      CHKERR cP.mField.get_moab().delete_entities(&meshset, 1);
    }
    MoFEMFunctionReturnHot(0);
  };
 
  // build tree
  CHKERR cutMeshPtr->buildTree();
  // refine and calculate level sets
  BitRefLevel bit = cP.mapBitLevel["mesh_cut"];
  std::size_t idx = 0;
  while (!bit.test(idx))
    ++idx;
  CHKERR cutMeshPtr->refineMesh(idx, nbRefBeforCut, nbRefBeforTrim, &fixedEdges,
                                verb, debug);
 
  auto set_last_bit_ref = [&]() {
    MoFEMFunctionBegin;
    if (nbRefBeforCut + nbRefBeforTrim) {
      BitRefLevel bits;
      for (int ll = 1; ll <= nbRefBeforCut + nbRefBeforTrim + 2; ++ll)
        bits.set(ll);
      BitRefLevel mask = bits;
      CHKERR cP.mField.delete_ents_by_bit_ref(BitRefLevel().set(20),
                                              BitRefLevel().set(20), true,
                                              QUIET, MF_NOT_THROW);
      Range ents;
      CHKERR bitRefPtr->getEntitiesByRefLevel(bits, mask, ents, verb);
      ents = subtract(ents, ents.subset_by_type(MBENTITYSET));
      CHKERR bitRefPtr->addBitRefLevel(ents, BitRefLevel().set(20));
    }
    MoFEMFunctionReturn(0);
  };
  CHKERR set_last_bit_ref();
 
  auto shift_after_ref = [&]() {
    MoFEMFunctionBegin;
    if (nbRefBeforCut + nbRefBeforTrim) {
      BitRefLevel mask;
      for (int ll = 0; ll <= nbRefBeforCut + nbRefBeforTrim + 2; ++ll)
        mask.set(ll);
      CHKERR bitRefPtr->shiftRightBitRef(nbRefBeforCut + nbRefBeforTrim, mask,
                                         verb, MF_NOT_THROW);
    }
    MoFEMFunctionReturn(0);
  };
  CHKERR shift_after_ref();
 
  // save coordinates to be later restored
  CHKERR getMeshOrgCoords();
 
  if (debug)
    CHKERR bitRefPtr->writeEntitiesAllBitLevelsByType(
        BitRefLevel().set(), MBTET, "all_after_ref_bits.vtk", "VTK", "");
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::setCutSurfaceFromFile() {
  MoFEMFunctionBegin;
  if (!cutSurfMeshName.empty()) {
    EntityHandle meshset_from_file;
    CHKERR meshsetMngPtr->getMeshset(cuttingSurfaceSidesetId, SIDESET,
                                     meshset_from_file);
    // FIXME:
    // if its cubit file (.cub) it will show the error :
    // [0]MOAB ERROR: -m-------------------- Error Message
    // ------------------------------------ [0]MOAB ERROR: Non-geometric entity
    // provided! [0]MOAB ERROR: set_sense() line 741 in src/GeomTopoTool.cpp
    // Failed to restore topology
    //
    // but the program will continue
    CHKERR cP.mField.get_moab().load_file(cutSurfMeshName.c_str(),
                                          &meshset_from_file, "");
    cutSurfMeshName.clear();
  }
  MoFEMFunctionReturn(0);
}
 
PetscErrorCode CPMeshCut::copySurface(const std::string save_mesh) {
  MoFEMFunctionBegin;
 
  sUrface.clear();
  fRont.clear();
 
  if (!meshsetMngPtr->checkMeshset(cuttingSurfaceSidesetId, SIDESET)) {
    PetscPrintf(PETSC_COMM_WORLD,
                "Sideset %d with cutting crack surface not found\n",
                cuttingSurfaceSidesetId);
    MoFEMFunctionReturnHot(0);
  }
  CHKERR meshsetMngPtr->getEntitiesByDimension(cuttingSurfaceSidesetId, SIDESET,
                                               2, sUrface, true);
 
  // Copy surface
  const_cast<Range &>(cutMeshPtr->getSurface()).clear();
  CHKERR cutMeshPtr->copySurface(sUrface, NULL, sHift, NULL, NULL, save_mesh);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::rebuildCrackSurface(const double factor,
                                              const std::string save_mesh,
                                              const int verb, bool debug) {
  Skinner skin(&cP.mField.get_moab());
  MoFEM::Interface &m_field = cP.mField;
  if (debugCut)
    debug = true;
  MoFEMFunctionBegin;
 
  BitRefLevel bit = cP.mapBitLevel["spatial_domain"];
 
  CHKERR setFixEdgesAndCorners(bit);
 
  Range bit_edges;
 
  Tag th_griffith_force;
  Tag th_griffith_force_projected;
  Tag th_w;
  Tag th_front_positions;
  Tag th_area_change;
  Tag th_material_force;
 
  boost::shared_ptr<OrientedBoxTreeTool> tree_body_skin_ptr;
  tree_body_skin_ptr = boost::shared_ptr<OrientedBoxTreeTool>(
      new OrientedBoxTreeTool(&cP.mField.get_moab(), "ROOTSETSURF", true));
  EntityHandle rootset_tree_body_skin;
 
  /// Maps nodes on old crack surface and new crack surface
  std::map<EntityHandle, EntityHandle> verts_map;
 
  // Seed pseudo-random generator with step. Algorithm is deterministic.
  srand(cP.startStep);
 
  struct SetFrontPositionsParameters {
 
    double tolRelSnapToFixedEdges; ///< Relative tolerance to snap edges to
                                   ///< fixed edges
    double tolAbsSnapToFixedEdges; ///< Relative tolerance to snap edges to
                                   ///< fixed edges
    double cornerFactor; ///< extension of the edge beyond corner when node on
                         ///< corner edge
    double skinFactor;   ///< extension of cutting surface when node on the skin
    double endNodeEdgeDeltaFactor; ///< End node edge extension
 
    SetFrontPositionsParameters() {
      tolRelSnapToFixedEdges = 0.1;
      tolAbsSnapToFixedEdges = 0.;
      cornerFactor = 0.2;
      skinFactor = 0.2;
      endNodeEdgeDeltaFactor = 0.2;
      ierr = getOptions();
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
    }
 
    MoFEMErrorCode getOptions() {
      MoFEMFunctionBegin;
 
      CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "cutting_",
                               "Set front positions", "none");
 
      CHKERR PetscOptionsScalar(
          "-surf_corner_factor",
          "extension of the edge beyond corner when node on corer edge", "",
          cornerFactor, &cornerFactor, PETSC_NULL);
      CHKERR PetscOptionsScalar(
          "-surf_skin_factor",
          "extension of cutting surface when node on the skin", "", skinFactor,
          &skinFactor, PETSC_NULL);
      CHKERR PetscOptionsScalar(
          "-end_node_edge_delta_factor",
          "extension of the front edge at the node on the surface", "",
          endNodeEdgeDeltaFactor, &endNodeEdgeDeltaFactor, PETSC_NULL);
      CHKERR PetscOptionsScalar("-snap_to_fixed_edge_rtol",
                                "snap to fixed edges (relative tolerances)", "",
                                tolRelSnapToFixedEdges, &tolRelSnapToFixedEdges,
                                PETSC_NULL);
      CHKERR PetscOptionsScalar("-snap_to_fixed_edge_atol",
                                "snap to fixed edges (absolute tolerances)", "",
                                tolAbsSnapToFixedEdges, &tolAbsSnapToFixedEdges,
                                PETSC_NULL);
 
      MOFEM_LOG_C("CPMeshCutWorld", Sev::inform,
                  "### Input parameter: -cutting_surf_corner_factor %6.4e",
                  cornerFactor);
      MOFEM_LOG_C("CPMeshCutWorld", Sev::inform,
                  "### Input parameter: -cutting_surf_skin_factor %6.4e",
                  skinFactor);
      MOFEM_LOG_C(
          "CPMeshCutWorld", Sev::inform,
          "### Input parameter: -cutting_end_node_edge_delta_factor %6.4e",
          endNodeEdgeDeltaFactor);
      MOFEM_LOG_C("CPMeshCutWorld", Sev::inform,
                  "### Input parameter: -cutting_snap_to_fixed_edge_rtol %6.4e",
                  tolRelSnapToFixedEdges);
      MOFEM_LOG_C("CPMeshCutWorld", Sev::inform,
                  "### Input parameter: -cutting_snap_to_fixed_edge_atol %6.4e",
                  tolAbsSnapToFixedEdges);
 
      ierr = PetscOptionsEnd();
      CHKERRG(ierr);
 
      MoFEMFunctionReturn(0);
    }
  };
 
  auto save_entities = [&](const std::string name, Range ents) {
    MoFEMFunctionBeginHot;
    if (!ents.empty()) {
      EntityHandle meshset;
      CHKERR cP.mField.get_moab().create_meshset(MESHSET_SET, meshset);
      CHKERR cP.mField.get_moab().add_entities(meshset, ents);
      CHKERR cP.mField.get_moab().write_file(name.c_str(), "VTK", "", &meshset,
                                             1);
      CHKERR cP.mField.get_moab().delete_entities(&meshset, 1);
    }
    MoFEMFunctionReturnHot(0);
  };
 
  auto get_tags = [&]() {
    MoFEMFunctionBegin;
    CHKERR cP.mField.get_moab().tag_get_handle("GRIFFITH_FORCE",
                                               th_griffith_force);
    CHKERR cP.mField.get_moab().tag_get_handle("GRIFFITH_FORCE_PROJECTED",
                                               th_griffith_force_projected);
    CHKERR cP.mField.get_moab().tag_get_handle("W", th_w);
    CHKERR cP.mField.get_moab().tag_get_handle("MATERIAL_FORCE",
                                               th_material_force);
 
    rval = cP.mField.get_moab().tag_get_handle("FRONT_POSITIONS",
                                               th_front_positions);
    if (rval == MB_SUCCESS)
      CHKERR cP.mField.get_moab().tag_delete(th_front_positions);
 
    double def_val[] = {0, 0, 0};
    CHKERR cP.mField.get_moab().tag_get_handle(
        "FRONT_POSITIONS", 3, MB_TYPE_DOUBLE, th_front_positions,
        MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
 
    rval = cP.mField.get_moab().tag_get_handle("AREA_CHANGE", th_area_change);
    if (rval == MB_SUCCESS)
      CHKERR cP.mField.get_moab().tag_delete(th_area_change);
 
    CHKERR cP.mField.get_moab().tag_get_handle(
        "AREA_CHANGE", 3, MB_TYPE_DOUBLE, th_area_change,
        MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
    MoFEMFunctionReturn(0);
  };
 
  auto get_surface = [this, &bit]() {
    Range surface;
    CHKERR meshsetMngPtr->getEntitiesByDimension(crackSurfaceId, SIDESET, 2,
                                                 surface, true);
    Range bit_faces;
    CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                   MBTRI, bit_faces);
    return intersect(surface, bit_faces);
  };
 
  auto get_crack_front = [&]() {
    Range crack_front;
    CHKERR meshsetMngPtr->getEntitiesByDimension(crackFrontId, SIDESET, 1,
                                                 crack_front, true);
    Range bit_edges;
    CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                   MBEDGE, bit_edges);
    return intersect(crack_front, bit_edges);
  };
 
  // blockset where crack propagation takes place
  auto get_cracked_body_blockset = [&]() {
    Range cracked_body_tets;
    CHKERR meshsetMngPtr->getEntitiesByDimension(
        crackedBodyBlockSetId, BLOCKSET, 3, cracked_body_tets, true);
    return cracked_body_tets;
  };
 
  // build tree to search distances to body surface
  auto build_body_skin_tree =
      [&](boost::shared_ptr<OrientedBoxTreeTool> &tree_body_skin_ptr,
          EntityHandle &rootset_tree_body_skin) {
        MoFEMFunctionBegin;
        Range last_bit_tets;
        CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(
            BitRefLevel().set(BITREFLEVEL_SIZE - 1), BitRefLevel().set(), MBTET,
            last_bit_tets);
        if (crackedBodyBlockSetId > 0)
          last_bit_tets = intersect(last_bit_tets, get_cracked_body_blockset());
        Range last_bit_tets_skin;
        CHKERR skin.find_skin(0, last_bit_tets, false, last_bit_tets_skin);
        CHKERR tree_body_skin_ptr->build(last_bit_tets_skin,
                                         rootset_tree_body_skin);
        MoFEMFunctionReturn(0);
      };
 
  // Set coords
  auto set_positions = [&](const Range nodes, Tag *th = nullptr) {
    MoFEMFunctionBegin;
    VectorDouble coords(3 * nodes.size());
    if (th)
      CHKERR cP.mField.get_moab().tag_get_data(*th, nodes, &*coords.begin());
    else
      CHKERR cP.mField.get_moab().get_coords(nodes, &*coords.begin());
    CHKERR cP.mField.get_moab().tag_set_data(th_front_positions, nodes,
                                             &*coords.begin());
    MoFEMFunctionReturn(0);
  };
 
  auto get_nodes = [this](Range edges) {
    Range verts;
    CHKERR cP.mField.get_moab().get_connectivity(edges, verts, true);
    return verts;
  };
 
  auto get_surface_skin = [this, &skin](auto &surface) {
    Range surface_skin;
    CHKERR skin.find_skin(0, surface, false, surface_skin);
    return surface_skin;
  };
 
  auto get_crack_front_ends = [this, &skin](Range crack_front) {
    Range end_nodes;
    CHKERR skin.find_skin(0, crack_front, false, end_nodes);
    return end_nodes;
  };
 
  auto get_edge_delta = [&](const EntityHandle edge, const EntityHandle node,
                            VectorDouble3 &delta) {
    MoFEMFunctionBegin;
    int num_nodes;
    const EntityHandle *conn;
    CHKERR cP.mField.get_moab().get_connectivity(edge, conn, num_nodes, true);
    VectorDouble3 n0(3), n1(3);
    if (conn[1] == node) {
      CHKERR cP.mField.get_moab().get_coords(&conn[0], 1, &n0[0]);
      CHKERR cP.mField.get_moab().get_coords(&conn[1], 1, &n1[0]);
    } else {
      CHKERR cP.mField.get_moab().get_coords(&conn[1], 1, &n0[0]);
      CHKERR cP.mField.get_moab().get_coords(&conn[0], 1, &n1[0]);
    }
    delta = n1 - n0;
    MoFEMFunctionReturn(0);
  };
 
  auto get_node_vec = [](std::vector<double> &v, const int nn) {
    return getVectorAdaptor(&v[3 * nn], 3);
  };
 
  auto check_if_node_is_in_range = [](const EntityHandle node,
                                      const Range &ents) {
    return (ents.find(node) != ents.end());
  };
 
  auto check_for_nan = [&](auto &v, std::string msg = "") {
    MoFEMFunctionBeginHot;
    for (auto d : v)
      if (std::isnan(d) || std::isinf(d))
        if (msg.size()) {
          SETERRQ2(cP.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                   "(%s) Not number %3.4e", msg.c_str(), d);
        } else {
          SETERRQ1(cP.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                   "Not number %3.4e", d);
        }
 
    MoFEMFunctionReturnHot(0);
  };
 
  auto set_tag = [&](const EntityHandle node, auto coords) {
    MoFEMFunctionBegin;
    CHKERR check_for_nan(coords, "set_tag");
    CHKERR cP.mField.get_moab().tag_set_data(th_front_positions, &node, 1,
                                             &coords[0]);
    MoFEMFunctionReturn(0);
  };
 
  auto get_str = [](auto v) { return boost::lexical_cast<std::string>(v); };
 
  auto get_tag_data = [this](Tag th, const Range &ents) {
    std::vector<double> v(3 * ents.size());
    CHKERR cP.mField.get_moab().tag_get_data(th, ents, &*v.begin());
    return v;
  };
 
  auto get_coords_vec = [this](const Range &ents) {
    std::vector<double> v(3 * ents.size());
    CHKERR cP.mField.get_moab().get_coords(ents, &*v.begin());
    return v;
  };
 
  auto get_prisms = [&]() {
    Range prisms;
    CHKERR cP.mField.get_moab().get_adjacencies(get_surface(), 3, false, prisms,
                                                moab::Interface::UNION);
    prisms = prisms.subset_by_type(MBPRISM);
    return prisms;
  };
 
  auto get_prims_surface = [&](Range prisms) {
    Range prisms_surface;
    for (auto prism : prisms) {
      EntityHandle face;
      CHKERR cP.mField.get_moab().side_element(prism, 2, 4, face);
      prisms_surface.insert(face);
    }
    return prisms_surface;
  };
 
  auto get_random_double = []() {
    return static_cast<double>(rand()) / static_cast<double>(RAND_MAX);
  };
 
  // Check if point is outside or inside of body
  auto is_point_is_outside = [&](const auto &coords) {
    VectorDouble3 random_ray(3);
    random_ray[0] =
        get_random_double() + std::numeric_limits<double>::epsilon();
    random_ray[1] = get_random_double();
    random_ray[2] = get_random_double();
    random_ray /= norm_2(random_ray);
    std::vector<double> distances_out;
    std::vector<EntityHandle> facets_out;
    CHKERR tree_body_skin_ptr->ray_intersect_triangles(
        distances_out, facets_out, rootset_tree_body_skin, 1e-12, &coords[0],
        &random_ray[0]);
    return ((facets_out.size() % 2) == 0);
  };
 
  // Check until results is consistent
  auto is_point_is_outside_with_check = [&](const auto &coords) {
    bool is_out = is_point_is_outside(coords);
    bool test_out;
    do {
      test_out = is_out;
      is_out = is_point_is_outside(coords);
    } while (test_out != is_out);
    return is_out;
  };
 
  auto get_proj_front_disp = [&](auto front_disp,
                                 auto griffith_force_projected) {
    const double mgn_force_projected = norm_2(griffith_force_projected);
    VectorDouble3 disp(3);
 
    if (mgn_force_projected > std::numeric_limits<double>::epsilon()) {
      griffith_force_projected /= mgn_force_projected;
      noalias(disp) = factor * griffith_force_projected *
                      fmax(0, inner_prod(front_disp, griffith_force_projected));
    } else {
      disp.clear();
    }
 
    return disp;
  };
 
  // get average area increment
  auto get_ave_a_change = [&](Range &crack_front_nodes,
                              std::vector<double> &area_change) {
    double ave_a_change = 0;
    int nn = 0;
    for (auto node : crack_front_nodes) {
      auto a_change = get_node_vec(area_change, nn);
      ave_a_change += norm_2(a_change);
    }
    ave_a_change /= crack_front_nodes.size();
    return ave_a_change;
  };
 
  auto get_ray_pout = [&](const auto coords, const auto ray) {
    std::vector<double> distances_out;
    std::vector<EntityHandle> facets_out;
    double ray_length = norm_2(ray);
    VectorDouble3 unit_ray = ray / ray_length;
    // send the ray in the direction of the edge
    CHKERR tree_body_skin_ptr->ray_intersect_triangles(
        distances_out, facets_out, rootset_tree_body_skin, 1e-12, &coords[0],
        &unit_ray[0], &ray_length);
    if (distances_out.empty())
      return -1.;
    else
      return *std::min_element(distances_out.begin(), distances_out.end());
  };
 
  CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(
      cP.mapBitLevel["spatial_domain"], BitRefLevel().set(), MBEDGE, bit_edges);
 
  SetFrontPositionsParameters front_params;
 
  CHKERR get_tags();
  sUrface = get_prims_surface(get_prisms());
 
  Range surface_skin = get_surface_skin(sUrface);
  Range crack_front = intersect(get_crack_front(), surface_skin);
  Range crack_front_nodes = get_nodes(crack_front);
  Range end_nodes = get_crack_front_ends(crack_front);
  Range surface_skin_nodes =
      subtract(subtract(get_nodes(surface_skin), crack_front_nodes),
               get_nodes(intersect(fixedEdges, surface_skin)));
  surface_skin_nodes.merge(intersect(cornerNodes, get_nodes(surface_skin)));
 
  CHKERR GriffithForceElement::SaveGriffithForceOnTags(
      cP.mField, th_area_change, sUrface, nullptr)();
  CHKERR set_positions(get_nodes(get_prisms()), nullptr);
  CHKERR build_body_skin_tree(tree_body_skin_ptr, rootset_tree_body_skin);
 
  if (debug)
    CHKERR save_entities("prisms_surface.vtk", get_prims_surface(get_prisms()));
  if (debug)
    CHKERR save_entities("surface_skin.vtk", surface_skin);
  if (debug)
    CHKERR save_entities("crack_front.vtk", crack_front);
  if (debug)
    CHKERR save_entities("crack_front_nodes.vtk", crack_front_nodes);
  if (debug)
    CHKERR save_entities("end_nodes.vtk", end_nodes);
 
  auto surface_snap = [&]() {
    MoFEMFunctionBegin;
    CHKERR cutMeshPtr->snapSurfaceToEdges(
        surface_skin, fixedEdges, front_params.tolRelSnapToFixedEdges,
        front_params.tolAbsSnapToFixedEdges, nullptr, debug);
    MoFEMFunctionReturn(0);
  };
 
  auto move_front_nodes = [&]() {
    MoFEMFunctionBegin;
 
    auto w = get_tag_data(th_w, crack_front_nodes);
    auto griffith_forces_projected =
        get_tag_data(th_griffith_force_projected, crack_front_nodes);
    auto area_change = get_tag_data(th_area_change, crack_front_nodes);
    auto mesh_node_positions = get_coords_vec(crack_front_nodes);
 
    const double ave_a_change =
        get_ave_a_change(crack_front_nodes, area_change);
 
    int nn = 0;
    for (auto node : crack_front_nodes) {
 
      auto front_disp = get_node_vec(w, nn);
      auto force_projected = get_node_vec(griffith_forces_projected, nn);
      auto a_change = get_node_vec(area_change, nn);
      auto position = get_node_vec(mesh_node_positions, nn);
 
      CHKERR check_for_nan(front_disp, "front_disp");
      CHKERR check_for_nan(force_projected, "force_projected");
      CHKERR check_for_nan(a_change, "a_change");
 
      VectorDouble3 disp = get_proj_front_disp(front_disp, force_projected);
 
      if (debug) {
        CHKERR check_for_nan(position);
        CHKERR check_for_nan(disp);
      }
 
      VectorDouble3 coords = position + disp;
 
      // check if nodes on fixed edges
      Range adj_crack_front_nodes_edges;
      CHKERR cP.mField.get_moab().get_adjacencies(&node, 1, 1, false,
                                                  adj_crack_front_nodes_edges,
                                                  moab::Interface::UNION);
      Range adj_crack_front_nodes_on_fixed_edges;
      adj_crack_front_nodes_on_fixed_edges =
          intersect(adj_crack_front_nodes_edges, fixedEdges);
      adj_crack_front_nodes_on_fixed_edges =
          intersect(adj_crack_front_nodes_on_fixed_edges, bit_edges);
 
      CHKERR check_for_nan(front_disp, "front_disp");
      CHKERR check_for_nan(coords, "coords");
      CHKERR check_for_nan(disp, "disp");
 
      enum Classifier {
        FREE = 0,
        FIXED_EDGE = 1 << 0,
        END_NODE = 1 << 1,
        FIXED_EDGE_KIND0 = 1 << 2,
        FIXED_EDGE_KIND1 = 1 << 3
      };
 
      unsigned int kind = FREE;
      if (check_if_node_is_in_range(node, get_nodes(fixedEdges)))
        kind |= FIXED_EDGE;
      if (check_if_node_is_in_range(node, end_nodes))
        kind |= END_NODE;
 
      if (kind == FREE) {
 
        // check pushing other nodes, which are not end nodes
        VectorDouble3 scaled_a_change =
            ave_a_change * a_change /
            (norm_2(a_change) + std::numeric_limits<double>::epsilon());
        if (is_point_is_outside_with_check(coords)) {
 
          coords += front_params.skinFactor * scaled_a_change;
          CHKERR check_for_nan(coords,
                               "Pushed through surface (out by factor)");
 
          MOFEM_LOG("CPMeshCutWorld", Sev::verbose)
              << "Pushed through surface (out by factor): ";
 
        } else {
 
          VectorDouble3 tmp_coords =
              coords + front_params.skinFactor * scaled_a_change;
          if (is_point_is_outside_with_check(tmp_coords)) {
            VectorDouble3 ray =
                scaled_a_change / (norm_2(scaled_a_change) +
                                   std::numeric_limits<double>::epsilon());
            const double dist = get_ray_pout(coords, scaled_a_change);
            if (dist > -0)
              coords += dist * ray;
 
            coords += front_params.skinFactor * scaled_a_change;
            CHKERR check_for_nan(
                coords, "Pushing through surface (out by skin factor)");
 
            MOFEM_LOG("CPMeshCutWorld", Sev::verbose)
                << "Pushing through surface (out by skin factor): ";
          }
        }
      }
 
      if (kind & END_NODE) {
 
        double min_dist = front_params.skinFactor * ave_a_change;
 
        if (kind & FIXED_EDGE) {
 
          VectorDouble3 pt_out(3);
          CHKERR cP.mField.getInterface<Tools>()->findMinDistanceFromTheEdges(
              &coords[0], 1, fixedEdges, &min_dist, &pt_out[0]);
          if (min_dist < front_params.skinFactor * ave_a_change)
            coords = pt_out;
 
          CHKERR check_for_nan(coords, "End node on the fixed edge");
 
          MOFEM_LOG("CPMeshCutWorld", Sev::verbose)
              << "End node on the fixed edge: ";
 
        } else {
          // Find corner edge and check if crack can go through it
 
          VectorDouble3 pt_out(3);
          CHKERR cP.mField.getInterface<Tools>()->findMinDistanceFromTheEdges(
              &coords[0], 1, fixedEdges, &min_dist, &pt_out[0]);
          if (min_dist < front_params.skinFactor * ave_a_change) {
            VectorDouble3 ray = pt_out - coords;
            if (inner_prod(a_change, ray) > 0) {
              VectorDouble3 scaled_ray =
                  front_params.skinFactor * ave_a_change * ray / norm_2(ray);
              coords += ray + scaled_ray;
              CHKERR check_for_nan(
                  coords, "Pushing through surface dist (close to edge)");
              MOFEM_LOG("CPMeshCutWorld", Sev::verbose)
                  << "Pushing through surface dist (close to edge): ";
            }
          }
        }
      }
 
      if (!adj_crack_front_nodes_on_fixed_edges.empty()) {
 
        const int nb_edges_on_fix_edges =
            adj_crack_front_nodes_on_fixed_edges.size() - 1;
        adj_crack_front_nodes_on_fixed_edges =
            intersect(adj_crack_front_nodes_on_fixed_edges, surface_skin);
 
        // adj node edge is on fixed edge and on the skin
        if (!adj_crack_front_nodes_on_fixed_edges.empty()) {
 
          const int nb_edges_on_fix_edges_and_surface_skin =
              adj_crack_front_nodes_on_fixed_edges.size();
          const int diff_nb_edges =
              nb_edges_on_fix_edges - nb_edges_on_fix_edges_and_surface_skin;
 
          // If is  >- in that case diff_nb_edges = 1
          // If is  >| in that case diff_nb_edges > 1
 
          if (diff_nb_edges > 1) {
 
            VectorDouble3 edge_delta;
            CHKERR get_edge_delta(adj_crack_front_nodes_on_fixed_edges[0], node,
                                  edge_delta);
            coords += front_params.cornerFactor * edge_delta;
            MOFEM_LOG("CPMeshCutWorld", Sev::verbose)
                << "Node on fixed edge (crack meeting fixed edge >|): ";
            CHKERR check_for_nan(
                coords, "Node on fixed edge (crack meeting fixed edge >|)");
 
            kind |= FIXED_EDGE_KIND0;
          }
 
        } else if ((kind & END_NODE) && (kind & FIXED_EDGE)) {
 
          // It s =|>
 
          // adj edge is on fixed edge but that edge is not on the skin
          Range adj_edges_on_skin =
              intersect(adj_crack_front_nodes_edges, surface_skin);
          if (!adj_edges_on_skin.empty()) {
 
            VectorDouble3 edge_delta;
            CHKERR get_edge_delta(adj_edges_on_skin[0], node, edge_delta);
            coords += front_params.cornerFactor * edge_delta;
            CHKERR check_for_nan(
                coords, "Node on fixed edge (crack crossing fixed edge =|>)");
 
            MOFEM_LOG("CPMeshCutWorld", Sev::verbose)
                << "Node on fixed edge (crack crossing fixed edge =|>): ";
 
            kind |= FIXED_EDGE_KIND1;
          }
        }
      }
 
      if (
 
          (kind & END_NODE) &&
 
          !(kind & FIXED_EDGE_KIND0) &&
 
          !(kind & FIXED_EDGE_KIND1)) {
 
        // is a end node on the skin, just go outside, in the direction
        // of crack front edged
        Range adj_edges = intersect(adj_crack_front_nodes_edges, crack_front);
        VectorDouble3 edge_delta;
        CHKERR get_edge_delta(adj_edges[0], node, edge_delta);
        coords += edge_delta * front_params.endNodeEdgeDeltaFactor;
        CHKERR check_for_nan(coords, "End node on the skin");
 
        // front_params.skinFactor;
        MOFEM_LOG("CPMeshCutWorld", Sev::verbose) << "End node on the skin: ";
      }
 
      MOFEM_LOG_C("CPMeshCutWorld", Sev::verbose, "Disp front %s at node %lu",
                  get_str(front_disp).c_str(), node);
 
      CHKERR set_tag(node, coords);
 
      ++nn;
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto smooth_front = [&]() {
    MoFEMFunctionBegin;
    Range smoothed_nodes = subtract(crack_front_nodes, end_nodes);
    std::vector<double> new_positions(3 * smoothed_nodes.size());
    int nn = 0;
    for (auto node : smoothed_nodes) {
      Range edges;
      CHKERR cP.mField.get_moab().get_adjacencies(&node, 1, 1, false, edges);
      edges = intersect(edges, crack_front);
      Range edges_nodes;
      CHKERR cP.mField.get_moab().get_connectivity(edges, edges_nodes, true);
      edges_nodes.erase(node);
      if (edges_nodes.size() != 2)
        SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                 "Expected two nodes but have %d", edges_nodes.size());
      VectorDouble3 node0(3), node1(3), node2(3);
      EntityHandle n0 = edges_nodes[0];
      EntityHandle n2 = edges_nodes[1];
      CHKERR cP.mField.get_moab().tag_get_data(th_front_positions, &n0, 1,
                                               &node0[0]);
      CHKERR cP.mField.get_moab().tag_get_data(th_front_positions, &node, 1,
                                               &node1[0]);
      CHKERR cP.mField.get_moab().tag_get_data(th_front_positions, &n2, 1,
                                               &node2[0]);
      VectorDouble3 delta0 = node0 - node1;
      VectorDouble3 delta2 = node2 - node1;
      const double l0 = norm_2(delta0);
      const double l2 = norm_2(delta2);
 
      if ((l0 / (l0 + l2)) < 0.5)
        node1 += 2 * (0.5 - (l0 / (l0 + l2))) * (delta2 / l2) * l0;
      else
        node1 += 2 * (0.5 - (l2 / (l0 + l2))) * (delta0 / l0) * l2;
 
      VectorAdaptor new_coords = get_node_vec(new_positions, nn);
      new_coords = node1;
      ++nn;
    }
    CHKERR cP.mField.get_moab().tag_set_data(th_front_positions, smoothed_nodes,
                                             &*new_positions.begin());
 
    MoFEMFunctionReturn(0);
  };
 
  auto delete_mofem_meshsets = [&]() {
    MoFEMFunctionBegin;
    CHKERR cP.mField.getInterface<MeshsetsManager>()->deleteMeshset(
        SIDESET, getSkinOfTheBodyId(), MF_ZERO);
    CHKERR cP.mField.getInterface<MeshsetsManager>()->deleteMeshset(
        SIDESET, getCrackSurfaceId(), MF_ZERO);
    CHKERR cP.mField.getInterface<MeshsetsManager>()->deleteMeshset(
        SIDESET, getCrackFrontId(), MF_ZERO);
    MoFEMFunctionReturn(0);
  };
 
  auto build_surface = [&]() {
    MoFEMFunctionBegin;
    // Create cutting surface from existing crack
    Tag th_interface_side;
    CHKERR cP.mField.get_moab().tag_get_handle("INTERFACE_SIDE",
                                               th_interface_side);
    Range new_surface;
    // Take nodes from top of the prims
    for (auto face : sUrface) {
 
      int side;
      CHKERR cP.mField.get_moab().tag_get_data(th_interface_side, &face, 1,
                                               &side);
      int num_nodes;
      const EntityHandle *conn;
      CHKERR cP.mField.get_moab().get_connectivity(face, conn, num_nodes, true);
 
      if (debug)
        if (!(conn[0] != conn[1] && conn[0] != conn[2] && conn[1] != conn[2]))
          SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                  "Imposible trinagle");
 
      MatrixDouble3by3 coords(3, 3);
      CHKERR cP.mField.get_moab().tag_get_data(th_front_positions, conn, 3,
                                               &coords(0, 0));
      if (debug)
        CHKERR check_for_nan(coords.data());
 
      EntityHandle new_verts[3] = {0, 0, 0};
      for (auto nn : {0, 1, 2}) {
        auto node_map = verts_map.find(conn[nn]);
        if (node_map == verts_map.end()) {
          EntityHandle new_node;
          CHKERR cP.mField.get_moab().create_vertex(&coords(nn, 0), new_node);
          verts_map[conn[nn]] = new_node;
        }
        new_verts[nn] = verts_map[conn[nn]];
      }
      for (auto nn : {0, 1, 2})
        if (new_verts[nn] == 0)
          SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                  "Imposible trinagle");
      if (side == 1) {
        EntityHandle n0 = new_verts[0];
        new_verts[0] = new_verts[1];
        new_verts[1] = n0;
      }
      EntityHandle ele;
      CHKERR cP.mField.get_moab().create_element(MBTRI, new_verts, 3, ele);
      new_surface.insert(ele);
    }
    sUrface.swap(new_surface);
    MoFEMFunctionReturn(0);
  };
 
  auto build_front = [&]() {
    MoFEMFunctionBegin;
    Range new_surface_skin;
    CHKERR skin.find_skin(0, sUrface, false, new_surface_skin);
    Range new_front;
    // Creating crack front from existing crack
    for (auto f : crack_front) {
      int num_nodes;
      const EntityHandle *conn;
      CHKERR cP.mField.get_moab().get_connectivity(f, conn, num_nodes, true);
      EntityHandle new_verts[num_nodes];
      for (int nn = 0; nn != 2; nn++) {
        auto node_map = verts_map.find(conn[nn]);
        if (node_map != verts_map.end())
          new_verts[nn] = node_map->second;
        else
          SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                  "Expected that vertex should be already created");
      }
      Range edges;
      CHKERR cP.mField.get_moab().get_adjacencies(new_verts, 2, 1, true, edges,
                                                  moab::Interface::INTERSECT);
      edges = intersect(edges, new_surface_skin);
      new_front.merge(edges);
    }
    fRont.swap(new_front);
    MoFEMFunctionReturn(0);
  };
 
  auto edges_snap = [&]() {
    MoFEMFunctionBegin;
    Range new_surface_skin;
    CHKERR skin.find_skin(0, sUrface, false, new_surface_skin);
    CHKERR cutMeshPtr->snapSurfaceToEdges(
        new_surface_skin, fixedEdges, front_params.tolRelSnapToFixedEdges,
        front_params.tolAbsSnapToFixedEdges, nullptr, debug);
    MoFEMFunctionReturn(0);
  };
 
  auto delete_meshsets_with_crack_surfaces = [&]() {
    MoFEMFunctionBegin;
    BitRefLevel mask;
    for (int ll = 0; ll != 19; ++ll)
      mask.set(ll);
    CHKERR bitRefPtr->shiftRightBitRef(19, mask, verb, MF_NOT_THROW);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR surface_snap();
  CHKERR move_front_nodes();
  CHKERR smooth_front();
  CHKERR delete_mofem_meshsets();
  CHKERR build_surface();
  if (debug)
    CHKERR save_entities("set_new_crack_surface.vtk", sUrface);
  CHKERR build_front();
  CHKERR edges_snap();
  if (debug)
    CHKERR save_entities("set_new_front.vtk", fRont);
  CHKERR delete_meshsets_with_crack_surfaces();
 
  // Set cutting inteface
  const_cast<Range &>(cutMeshPtr->getSurface()).clear();
  const_cast<Range &>(cutMeshPtr->getFront()).clear();
  CHKERR cutMeshPtr->setSurface(sUrface);
  CHKERR cutMeshPtr->setFront(fRont);
 
  if (!save_mesh.empty())
    CHKERR save_entities(save_mesh, unite(sUrface, fRont));
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::setFixEdgesAndCorners(const BitRefLevel bit) {
  MoFEMFunctionBegin;
  fixedEdges.clear();
  cornerNodes.clear();
 
  if (meshsetMngPtr->checkMeshset(edgesBlockSet, BLOCKSET))
    CHKERR meshsetMngPtr->getEntitiesByDimension(edgesBlockSet, BLOCKSET, 1,
                                                 fixedEdges, true);
 
  if (meshsetMngPtr->checkMeshset(vertexNodeSet, NODESET) &&
      meshsetMngPtr->checkMeshset(vertexBlockSet, BLOCKSET))
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "There are both BLOCKSET and NODESET with ID equal to input "
            "parameter -vertex_block_set");
 
  if (meshsetMngPtr->checkMeshset(vertexBlockSet, BLOCKSET))
    CHKERR meshsetMngPtr->getEntitiesByDimension(vertexBlockSet, BLOCKSET, 0,
                                                 cornerNodes, true);
 
  if (meshsetMngPtr->checkMeshset(vertexNodeSet, NODESET))
    CHKERR meshsetMngPtr->getEntitiesByDimension(vertexNodeSet, NODESET, 0,
                                                 cornerNodes, true);
 
  CHKERR bitRefPtr->filterEntitiesByRefLevel(bit, BitRefLevel().set(),
                                             fixedEdges);
  CHKERR bitRefPtr->filterEntitiesByRefLevel(bit, BitRefLevel().set(),
                                             cornerNodes);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::cutMesh(int &first_bit, const bool debug) {
  MoFEM::Interface &m_field = cP.mField;
  MoFEMFunctionBegin;
 
  auto save_mesh = [&](const std::string file, const Range ents) {
    MoFEMFunctionBegin;
    if (m_field.get_comm_rank() == 0) {
      if (!ents.empty()) {
        EntityHandle meshset;
        CHKERR m_field.get_moab().create_meshset(MESHSET_SET, meshset);
        CHKERR m_field.get_moab().add_entities(meshset, ents);
        CHKERR m_field.get_moab().write_file(file.c_str(), "VTK", "", &meshset,
                                             1);
        CHKERR m_field.get_moab().delete_entities(&meshset, 1);
      }
    }
    MoFEMFunctionReturn(0);
  };
 
  auto get_entities_from_contact_meshset = [&](BitRefLevel &bit) {
    MoFEMFunctionBegin;
    cP.contactMeshsetFaces.clear();
    Range tris;
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(cP.mField, BLOCKSET, cit)) {
      if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
        tris.clear();
        CHKERR meshsetMngPtr->getEntitiesByDimension(cit->getMeshsetId(),
                                                     BLOCKSET, 2, tris, true);
        cP.contactMeshsetFaces.merge(tris);
      }
    }
    CHKERR bitRefPtr->filterEntitiesByRefLevel(bit, BitRefLevel().set(),
                                               cP.contactMeshsetFaces);
    MoFEMFunctionReturn(0);
  };
 
  auto get_entities_from_constrained_interface = [&](BitRefLevel &bit) {
    MoFEMFunctionBegin;
    cP.constrainedInterface.clear();
    Range tris;
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(cP.mField, BLOCKSET, cit)) {
      if (cit->getName().compare(0, 15, "INT_CONSTRAINED") == 0) {
        tris.clear();
        CHKERR meshsetMngPtr->getEntitiesByDimension(cit->getMeshsetId(),
                                                     BLOCKSET, 2, tris, true);
        cP.constrainedInterface.merge(tris);
      }
    }
    CHKERR bitRefPtr->filterEntitiesByRefLevel(bit, BitRefLevel().set(),
                                               cP.constrainedInterface);
    MoFEMFunctionReturn(0);
  };
 
  Range all_constrained_faces;
  auto get_all_constrained_faces = [&](BitRefLevel &bit) {
    MoFEMFunctionBegin;
    all_constrained_faces.clear();
 
    if (!cP.ignoreContact) {
      CHKERR get_entities_from_contact_meshset(bit);
      all_constrained_faces.merge(cP.contactMeshsetFaces);
    }
 
    CHKERR get_entities_from_constrained_interface(bit);
    all_constrained_faces.merge(cP.constrainedInterface);
 
    MoFEMFunctionReturn(0);
  };
 
  if (debug)
    CHKERR save_mesh("crack_surface_before_merge.vtk",
                     cutMeshPtr->getSurface());
 
  if (debug) {
    CHKERR cP.mField.getInterface<BitRefManager>()->writeBitLevelByType(
        cP.mapBitLevel["mesh_cut"], BitRefLevel().set(), MBVERTEX,
        "verts_on_mesh_cut_level.vtk", "VTK", "", false);
    CHKERR cP.mField.getInterface<BitRefManager>()->writeBitLevelByType(
        cP.mapBitLevel["mesh_cut"], BitRefLevel().set(), MBEDGE,
        "edges_on_mesh_cut_level.vtk", "VTK", "", false);
    CHKERR cP.mField.getInterface<BitRefManager>()->writeBitLevelByType(
        cP.mapBitLevel["mesh_cut"], BitRefLevel().set(), MBTRI,
        "triangles_on_mesh_cut_level.vtk", "VTK", "", false);
    CHKERR cP.mField.getInterface<BitRefManager>()->writeBitLevelByType(
        cP.mapBitLevel["mesh_cut"], BitRefLevel().set(), MBTET,
        "tets_on_mesh_cut_level.vtk", "VTK", "", false);
  }
 
  CHKERR get_all_constrained_faces(cP.mapBitLevel["mesh_cut"]);
  if (!all_constrained_faces.empty()) {
    CHKERR cutMeshPtr->setConstrainSurface(all_constrained_faces);
    if (debug)
      CHKERR save_mesh("constrained_surface_before_merge.vtk",
                       all_constrained_faces);
  }
 
  // Cut mesh, trim surface and merge bad edges
  CHKERR setFixEdgesAndCorners(cP.mapBitLevel["mesh_cut"]);
 
  CHKERR cutMeshPtr->cutTrimAndMerge(first_bit, fractionLevel, nullptr, tolCut,
                                     tolCutClose, tolTrimClose, fixedEdges,
                                     cornerNodes, true, debug);
 
  BitRefLevel bit_after_merge = BitRefLevel().set(first_bit);
  if (removePathologicalFrontTris)
    CHKERR cutMeshPtr->removePathologicalFrontTris(
        bit_after_merge, const_cast<Range &>(cutMeshPtr->getMergedSurfaces()));
 
  if (debug)
    CHKERR save_mesh("crack_surface_after_merge.vtk",
                     cutMeshPtr->getMergedSurfaces());
 
  CHKERR get_all_constrained_faces(bit_after_merge);
  if (!all_constrained_faces.empty() && debug) {
    CHKERR save_mesh("constrained_surface_after_merge.vtk",
                     all_constrained_faces);
  }
 
  auto add_entities_to_crack_meshset = [&]() {
    MoFEMFunctionBegin;
    if (!meshsetMngPtr->checkMeshset(crackSurfaceId, SIDESET))
      CHKERR meshsetMngPtr->addMeshset(SIDESET, crackSurfaceId);
    Range crack_faces = cutMeshPtr->getMergedSurfaces();
 
    if (!all_constrained_faces.empty() && crackedBodyBlockSetId < 0) {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_FOUND,
              "Block set id of the cracked body is required for a simulation "
              "with contact and/or constrained interface, use "
              "'-my_cracked_body_block_set' parameter");
    }
    if (crackedBodyBlockSetId > 0) {
      Range cracked_body_tets, cracked_body_tris;
      CHKERR meshsetMngPtr->getEntitiesByDimension(
          crackedBodyBlockSetId, BLOCKSET, 3, cracked_body_tets, true);
      CHKERR m_field.get_moab().get_adjacencies(cracked_body_tets, 2, false,
                                                cracked_body_tris,
                                                moab::Interface::UNION);
      crack_faces = intersect(crack_faces, cracked_body_tris);
      if (debug)
        CHKERR save_mesh("crack_faces_intersected_with_body_block.vtk",
                         crack_faces);
    }
 
    CHKERR meshsetMngPtr->addEntitiesToMeshset(SIDESET, crackSurfaceId,
                                               crack_faces);
    MoFEMFunctionReturn(0);
  };
 
  auto clean_copied_cutting_surface_entities = [&]() {
    MoFEMFunctionBegin;
    Range surface_verts;
    CHKERR m_field.get_moab().get_connectivity(cutMeshPtr->getSurface(),
                                               surface_verts);
    Range surface_edges;
    CHKERR m_field.get_moab().get_adjacencies(cutMeshPtr->getSurface(), 1,
                                              false, surface_edges,
                                              moab::Interface::UNION);
    // Remove entities from meshsets
    Range meshsets;
    CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
                                                   false);
    for (Range::iterator mit = meshsets.begin(); mit != meshsets.end(); mit++) {
      CHKERR m_field.get_moab().remove_entities(*mit, cutMeshPtr->getSurface());
      CHKERR m_field.get_moab().remove_entities(*mit, surface_edges);
      CHKERR m_field.get_moab().remove_entities(*mit, surface_verts);
    }
 
    CHKERR m_field.get_moab().delete_entities(cutMeshPtr->getSurface());
    CHKERR m_field.get_moab().delete_entities(surface_edges);
    CHKERR m_field.get_moab().delete_entities(surface_verts);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_entities_to_crack_meshset();
  CHKERR clean_copied_cutting_surface_entities();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::findBodySkin(const BitRefLevel bit, const int verb,
                                       const bool debug,
                                       std::string file_name) {
  MoFEM::Interface &m_field = cP.mField;
  moab::Interface &moab = m_field.get_moab();
  MoFEMFunctionBegin;
 
  auto save_mesh = [&](const std::string file, const Range ents) {
    MoFEMFunctionBegin;
    if (m_field.get_comm_rank() == 0) {
      EntityHandle meshset;
      CHKERR moab.create_meshset(MESHSET_SET, meshset);
      CHKERR moab.add_entities(meshset, ents);
      CHKERR moab.write_file(file.c_str(), "VTK", "", &meshset, 1);
      CHKERR moab.delete_entities(&meshset, 1);
    }
    MoFEMFunctionReturn(0);
  };
 
  // Find body skin, take it from meshset if is found, if not create meshest,
  // and add skin entities
  cP.bodySkin.clear();
  if (!meshsetMngPtr->checkMeshset(skinOfTheBodyId, SIDESET)) {
    Range tets, prisms;
    CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                   MBTET, tets);
    CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                   MBPRISM, prisms);
    prisms = subtract(prisms, cP.contactElements);
    prisms = subtract(prisms, cP.mortarContactElements);
    Range prisms_faces;
    for (auto p : prisms) {
      EntityHandle f;
      CHKERR moab.side_element(p, 2, 3, f);
      prisms_faces.insert(f);
      CHKERR moab.side_element(p, 2, 4, f);
      prisms_faces.insert(f);
    }
    Range skin_faces;
    Skinner skin(&m_field.get_moab());
    CHKERR skin.find_skin(0, tets, false, skin_faces);
    skin_faces = subtract(skin_faces, prisms_faces);
    CHKERR meshsetMngPtr->addMeshset(SIDESET, skinOfTheBodyId);
    EntityHandle meshset;
    CHKERR meshsetMngPtr->getMeshset(skinOfTheBodyId, SIDESET, meshset);
    CHKERR moab.add_entities(meshset, skin_faces);
  }
 
  // take skin entities from meshset
  CHKERR meshsetMngPtr->getEntitiesByDimension(skinOfTheBodyId, SIDESET, 2,
                                               cP.bodySkin, true);
  CHKERR bitRefPtr->filterEntitiesByRefLevel(bit, BitRefLevel().set(),
                                             cP.bodySkin);
 
  if (verb >= VERY_VERBOSE)
    if (m_field.get_comm_rank() == 0)
      cerr << cP.bodySkin << endl;
 
  if (debug)
    save_mesh(file_name, cP.bodySkin);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::findCrack(const BitRefLevel bit, const int verb,
                                    const bool debug) {
  MoFEM::Interface &m_field = cP.mField;
  moab::Interface &moab = cP.mField.get_moab();
  MoFEMFunctionBegin;
 
  auto save_mesh = [&](const std::string file, const Range ents) {
    MoFEMFunctionBegin;
    if (m_field.get_comm_rank() == 0) {
      EntityHandle meshset;
      CHKERR moab.create_meshset(MESHSET_SET, meshset);
      CHKERR moab.add_entities(meshset, ents);
      CHKERR moab.write_file(file.c_str(), "VTK", "", &meshset, 1);
      CHKERR moab.delete_entities(&meshset, 1);
    }
    MoFEMFunctionReturn(0);
  };
 
  cP.crackFaces.clear();
  if (!meshsetMngPtr->checkMeshset(crackSurfaceId, SIDESET)) {
    PetscPrintf(PETSC_COMM_WORLD, "Sideset %d with crack surface not found\n",
                crackSurfaceId);
    MoFEMFunctionReturnHot(0);
  }
  CHKERR meshsetMngPtr->getEntitiesByDimension(crackSurfaceId, SIDESET, 2,
                                               cP.crackFaces, true);
  CHKERR bitRefPtr->filterEntitiesByRefLevel(bit, BitRefLevel().set(),
                                             cP.crackFaces);
 
  // Find crack front
  cP.crackFront.clear();
 
  if (!meshsetMngPtr->checkMeshset(crackFrontId, SIDESET)) {
    Range skin_faces_edges;
    Range constrained_faces =
        unite(cP.contactMeshsetFaces, cP.constrainedInterface);
    CHKERR moab.get_adjacencies(unite(cP.bodySkin, constrained_faces), 1, false,
                                skin_faces_edges, moab::Interface::UNION);
    Range crack_front;
    Skinner skin(&m_field.get_moab());
    CHKERR skin.find_skin(0, cP.crackFaces, false, crack_front);
    crack_front = subtract(crack_front, skin_faces_edges);
    CHKERR meshsetMngPtr->addMeshset(SIDESET, crackFrontId);
    EntityHandle meshset;
    CHKERR meshsetMngPtr->getMeshset(crackFrontId, SIDESET, meshset);
    CHKERR moab.add_entities(meshset, crack_front);
  }
 
  CHKERR meshsetMngPtr->getEntitiesByDimension(crackFrontId, SIDESET, 1,
                                               cP.crackFront, true);
 
  Range crack_faces_edges;
  CHKERR moab.get_adjacencies(cP.crackFaces, 1, false, crack_faces_edges,
                              moab::Interface::UNION);
  cP.crackFront = intersect(cP.crackFront, crack_faces_edges);
 
  if (verb >= VERBOSE)
    if (m_field.get_comm_rank() == 0)
      cerr << cP.crackFront << endl;
 
  if (debug)
    save_mesh("out_crack_and_front.vtk", unite(cP.crackFront, cP.crackFaces));
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::findCrackFromPrisms(const BitRefLevel bit,
                                              const int verb,
                                              const bool debug) {
  MoFEM::Interface &m_field = cP.mField;
  moab::Interface &moab = cP.mField.get_moab();
  MoFEMFunctionBegin;
 
  auto save_mesh = [&](const std::string file, const Range ents) {
    MoFEMFunctionBegin;
    if (m_field.get_comm_rank() == 0) {
      EntityHandle meshset;
      CHKERR moab.create_meshset(MESHSET_SET, meshset);
      CHKERR moab.add_entities(meshset, ents);
      CHKERR moab.write_file(file.c_str(), "VTK", "", &meshset, 1);
      CHKERR moab.delete_entities(&meshset, 1);
    }
    MoFEMFunctionReturn(0);
  };
 
  Range prisms_level;
  CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                 MBPRISM, prisms_level);
  prisms_level = subtract(prisms_level, cP.contactElements);
  prisms_level = subtract(prisms_level, cP.mortarContactElements);
 
  cP.oneSideCrackFaces.clear();
  cP.otherSideCrackFaces.clear();
  for (auto p : prisms_level) {
    EntityHandle face3, face4;
    CHKERR m_field.get_moab().side_element(p, 2, 3, face3);
    CHKERR m_field.get_moab().side_element(p, 2, 4, face4);
    cP.oneSideCrackFaces.insert(face3);
    cP.otherSideCrackFaces.insert(face4);
  }
  cP.crackFaces = unite(cP.otherSideCrackFaces, cP.oneSideCrackFaces);
  CHKERR meshsetMngPtr->deleteMeshset(SIDESET, crackSurfaceId);
  CHKERR meshsetMngPtr->addMeshset(SIDESET, crackSurfaceId);
  CHKERR meshsetMngPtr->addEntitiesToMeshset(SIDESET, crackSurfaceId,
                                             cP.crackFaces);
 
  if (cP.otherSideConstrains == PETSC_TRUE)
    cP.oneSideCrackFaces.swap(cP.otherSideCrackFaces);
 
  if (debug)
    if (m_field.get_comm_rank() == 0) {
      CHKERR save_mesh("out_one_side.vtk", cP.oneSideCrackFaces);
      if (!cP.otherSideCrackFaces.empty())
        CHKERR save_mesh("out_other_side.vtk", cP.otherSideCrackFaces);
    }
 
  if (!cP.otherSideCrackFaces.empty()) {
    if (cP.otherSideCrackFaces.size() != cP.oneSideCrackFaces.size())
      SETERRQ2(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
               "Wrong number of faces on both sides of crack %d != %d",
               cP.otherSideCrackFaces.size(), cP.oneSideCrackFaces.size());
  }
 
  // Find crack front
  cP.crackFront.clear();
 
  CHKERR meshsetMngPtr->deleteMeshset(SIDESET, crackFrontId, MF_ZERO);
  auto create_front_meshset = [&]() {
    MoFEMFunctionBegin;
    Range skin_faces_edges;
    CHKERR moab.get_adjacencies(unite(cP.bodySkin, cP.constrainedInterface), 1,
                                false, skin_faces_edges,
                                moab::Interface::UNION);
    Range crack_front;
    Skinner skin(&m_field.get_moab());
    CHKERR skin.find_skin(0, cP.oneSideCrackFaces, false, crack_front);
    crack_front = subtract(crack_front, skin_faces_edges);
    CHKERR meshsetMngPtr->addMeshset(SIDESET, crackFrontId);
    CHKERR meshsetMngPtr->addEntitiesToMeshset(SIDESET, crackFrontId,
                                               crack_front);
    MoFEMFunctionReturn(0);
  };
  CHKERR create_front_meshset();
  CHKERR meshsetMngPtr->getEntitiesByDimension(crackFrontId, SIDESET, 1,
                                               cP.crackFront, true);
 
  if (verb >= VERBOSE)
    if (m_field.get_comm_rank() == 0)
      cerr << cP.crackFront << endl;
 
  if (verb >= VERBOSE) {
    CHKERR PetscPrintf(m_field.get_comm(),
                       "number of Crack Surfaces Faces = %d\n",
                       cP.crackFaces.size());
    CHKERR PetscPrintf(m_field.get_comm(),
                       "number of Crack Surfaces Faces On One Side = %d\n",
                       cP.oneSideCrackFaces.size());
    CHKERR PetscPrintf(m_field.get_comm(),
                       "number of Crack Surfaces Faces On Other Side = %d\n",
                       cP.otherSideCrackFaces.size());
    CHKERR PetscPrintf(m_field.get_comm(), "number of Crack Front Edges = %d\n",
                       cP.crackFront.size());
    Range prisms_level;
    CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                   MBPRISM, prisms_level);
    CHKERR PetscPrintf(m_field.get_comm(), "number of Prisms = %d\n",
                       prisms_level.size());
  }
 
  if (debug)
    save_mesh("out_crack_and_front.vtk", unite(cP.crackFront, cP.crackFaces));
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::splitFaces(const int verb, const bool debug) {
  MoFEM::Interface &m_field = cP.mField;
  moab::Interface &moab = m_field.get_moab();
  PrismInterface *prism_interface;
  MoFEMFunctionBegin;
  CHKERR m_field.getInterface(prism_interface);
 
  BitRefLevel bit = cP.mapBitLevel["mesh_cut"];
  BitRefLevel bit_split = cP.mapBitLevel["spatial_domain"];
 
  // split sides
  EntityHandle meshset_interface;
  if (!meshsetMngPtr->checkMeshset(crackSurfaceId, SIDESET)) {
    PetscPrintf(PETSC_COMM_WORLD, "Sideset %d with crack surface not found\n",
                crackSurfaceId);
    Range ents;
    CHKERR bitRefPtr->getEntitiesByRefLevel(bit, BitRefLevel().set(), ents);
    CHKERR bitRefPtr->addBitRefLevel(ents, bit_split);
    MoFEMFunctionReturnHot(0);
  } else {
 
    // clear tags data created by previous use of prism interface
    Tag th_interface_side;
    rval =
        m_field.get_moab().tag_get_handle("INTERFACE_SIDE", th_interface_side);
    if (rval == MB_SUCCESS)
      rval = m_field.get_moab().tag_delete(th_interface_side);
 
    // create meshset and insert entities from bit level to it
    EntityHandle bit_meshset;
    CHKERR m_field.get_moab().create_meshset(MESHSET_SET, bit_meshset);
 
    CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                   MBTET, bit_meshset);
    CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                   MBPRISM, bit_meshset);
 
    auto save_mesh = [&](const std::string file, const EntityHandle meshset) {
      MoFEMFunctionBegin;
      if (m_field.get_comm_rank() == 0)
        CHKERR moab.write_file(file.c_str(), "VTK", "", &meshset, 1);
      MoFEMFunctionReturn(0);
    };
 
    CHKERR meshsetMngPtr->getMeshset(crackSurfaceId, SIDESET,
                                     meshset_interface);
    if (debug)
      save_mesh("crack_surface_before_split.vtk", meshset_interface);
 
    CHKERR prism_interface->getSides(meshset_interface, bit, true, QUIET);
    CHKERR prism_interface->splitSides(bit_meshset, bit_split,
                                       meshset_interface, true, true, QUIET);
    // add refined ent to cubit meshsets
    for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, cubit_it)) {
      EntityHandle update_meshset = cubit_it->meshset;
      CHKERR bitRefPtr->updateMeshsetByEntitiesChildren(
          update_meshset, bit, BitRefLevel().set(), bit_split, bit_split,
          update_meshset, MBMAXTYPE, true);
    }
    CHKERR moab.delete_entities(&bit_meshset, 1);
 
    if (cP.isPartitioned) {
      ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
      Tag part_tag = pcomm->part_tag();
      Range tagged_sets;
      CHKERR m_field.get_moab().get_entities_by_type_and_tag(
          0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
          moab::Interface::UNION);
      for (Range::iterator mit = tagged_sets.begin(); mit != tagged_sets.end();
           mit++) {
        int part;
        CHKERR moab.tag_get_data(part_tag, &*mit, 1, &part);
        CHKERR bitRefPtr->updateMeshsetByEntitiesChildren(
            *mit, bit, BitRefLevel().set(), bit_split, bit_split, *mit, MBTET,
            true);
        Range ents3d;
        CHKERR moab.get_entities_by_dimension(*mit, 3, ents3d, true);
        for (Range::iterator eit = ents3d.begin(); eit != ents3d.end(); eit++) {
          CHKERR moab.tag_set_data(part_tag, &*eit, 1, &part);
        }
        for (int dd = 0; dd != 3; dd++) {
          Range ents;
          CHKERR moab.get_adjacencies(ents3d, dd, false, ents,
                                      moab::Interface::UNION);
          for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
            CHKERR moab.tag_set_data(part_tag, &*eit, 1, &part);
          }
        }
      }
    }
 
    if (debug) {
      Range tets_level;
      CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(
          bit_split, BitRefLevel().set(), MBTET, tets_level);
      Range prisms_level;
      CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(
          bit_split, BitRefLevel().set(), MBPRISM, prisms_level);
      EntityHandle meshset_level;
      CHKERR moab.create_meshset(MESHSET_SET, meshset_level);
      CHKERR moab.add_entities(meshset_level, tets_level);
      CHKERR moab.write_file("out_tets_level.vtk", "VTK", "", &meshset_level,
                             1);
      CHKERR moab.delete_entities(&meshset_level, 1);
      EntityHandle meshset_prims_level;
      CHKERR moab.create_meshset(MESHSET_SET, meshset_prims_level);
      CHKERR moab.add_entities(meshset_prims_level, prisms_level);
      CHKERR moab.write_file("out_prisms_level.vtk", "VTK", "",
                             &meshset_prims_level, 1);
      CHKERR moab.delete_entities(&meshset_prims_level, 1);
    }
  }
  if (debug) {
    Range tets_level;
    CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit_split, bit_split, MBTET,
                                                   tets_level);
    EntityHandle meshset_level;
    CHKERR moab.create_meshset(MESHSET_SET, meshset_level);
    CHKERR moab.add_entities(meshset_level, tets_level);
    CHKERR moab.write_file("out_tets_level_only.vtk", "VTK", "", &meshset_level,
                           1);
    CHKERR moab.delete_entities(&meshset_level, 1);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::addMortarContactPrisms(const int verb,
                                                 const bool debug) {
  MoFEM::Interface &m_field = cP.mField;
  moab::Interface &moab = m_field.get_moab();
  MoFEMFunctionBegin;
 
  Range integration_triangles;
  cP.mortarContactMasterFaces.clear();
  cP.mortarContactSlaveFaces.clear();
  cP.mortarContactElements.clear();
 
  Range all_masters, all_slaves;
 
  std::vector<std::pair<Range, Range>> contact_surface_pairs; // <Slave, Master>
 
  CHKERR MortarContactStructures::findContactSurfacePairs(
      m_field, contact_surface_pairs);
 
  cP.contactSearchMultiIndexPtr =
      boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>(
          new ContactSearchKdTree::ContactCommonData_multiIndex());
 
  for (auto &contact_surface_pair : contact_surface_pairs) {
 
    Range slave_tris = contact_surface_pair.first;
    Range master_tris = contact_surface_pair.second;
 
    CHKERR bitRefPtr->filterEntitiesByRefLevel(cP.mapBitLevel["spatial_domain"],
                                               BitRefLevel().set(), slave_tris);
    CHKERR bitRefPtr->filterEntitiesByRefLevel(
        cP.mapBitLevel["spatial_domain"], BitRefLevel().set(), master_tris);
 
    all_slaves.merge(slave_tris);
    all_masters.merge(master_tris);
 
    if (slave_tris.empty() || master_tris.empty())
      continue;
 
    ContactSearchKdTree contact_search_kd_tree(m_field);
 
    // create kd_tree with master_surface only
    CHKERR contact_search_kd_tree.buildTree(master_tris);
 
    CHKERR contact_search_kd_tree.contactSearchAlgorithm(
        master_tris, slave_tris, cP.contactSearchMultiIndexPtr,
        cP.mortarContactElements);
  }
 
  if (all_slaves.empty() || all_masters.empty())
    MoFEMFunctionReturnHot(0);
 
  Range tris_from_prism;
  // find actual masters and slave used
  CHKERR m_field.get_moab().get_adjacencies(cP.mortarContactElements, 2, true,
                                            tris_from_prism,
                                            moab::Interface::UNION);
 
  tris_from_prism = tris_from_prism.subset_by_type(MBTRI);
  cP.mortarContactMasterFaces = intersect(tris_from_prism, all_masters);
  cP.mortarContactSlaveFaces = intersect(tris_from_prism, all_slaves);
 
  EntityHandle meshset_slave_master_prisms;
  CHKERR moab.create_meshset(MESHSET_SET, meshset_slave_master_prisms);
 
  CHKERR
  moab.add_entities(meshset_slave_master_prisms, cP.mortarContactElements);
 
  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
      meshset_slave_master_prisms, 3, cP.mapBitLevel["spatial_domain"]);
 
  typedef ContactSearchKdTree::ContactCommonData_multiIndex::index<
      ContactSearchKdTree::Prism_tag>::type::iterator ItMultIndexPrism;
 
  for (Range::iterator it_tri = cP.mortarContactElements.begin();
       it_tri != cP.mortarContactElements.end(); ++it_tri) {
 
    ItMultIndexPrism it_mult_index_prism =
        cP.contactSearchMultiIndexPtr->get<ContactSearchKdTree::Prism_tag>()
            .find(*it_tri);
 
    Range range_poly_tris;
    range_poly_tris.clear();
    range_poly_tris = it_mult_index_prism->get()->commonIntegratedTriangle;
    integration_triangles.insert(range_poly_tris.begin(),
                                 range_poly_tris.end());
  }
 
  EntityHandle ent_integration_vertex;
  CHKERR moab.create_meshset(MESHSET_SET, ent_integration_vertex);
 
  EntityHandle ent_integration_edge;
  CHKERR moab.create_meshset(MESHSET_SET, ent_integration_edge);
 
  Range edges, verts;
  CHKERR m_field.get_moab().get_adjacencies(integration_triangles, 1, true,
                                            edges, moab::Interface::UNION);
  CHKERR moab.add_entities(ent_integration_edge, edges);
 
  CHKERR m_field.get_moab().get_connectivity(integration_triangles, verts,
                                             true);
  CHKERR moab.add_entities(ent_integration_vertex, verts);
 
  EntityHandle ent_integration_tris;
  CHKERR moab.create_meshset(MESHSET_SET, ent_integration_tris);
  CHKERR moab.add_entities(ent_integration_tris, integration_triangles);
 
  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
      ent_integration_vertex, 0, cP.mapBitLevel["spatial_domain"]);
  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
      ent_integration_edge, 1, cP.mapBitLevel["spatial_domain"]);
  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
      ent_integration_tris, 2, cP.mapBitLevel["spatial_domain"]);
 
  CHKERR moab.delete_entities(&ent_integration_tris, 1);
 
  const std::string output_name = "slave_master_prisms.vtk";
  CHKERR moab.write_mesh(output_name.c_str(), &meshset_slave_master_prisms, 1);
 
  // get integration tris of each prism
  CHKERR moab.delete_entities(&meshset_slave_master_prisms, 1);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::insertContactInterface(const int verb,
                                                 const bool debug) {
  MoFEM::Interface &m_field = cP.mField;
  moab::Interface &moab = m_field.get_moab();
  PrismInterface *prism_interface;
  MoFEMFunctionBegin;
  CHKERR m_field.getInterface(prism_interface);
 
  std::vector<BitRefLevel> bit_levels;
  bit_levels.push_back(cP.mapBitLevel["mesh_cut"]);
 
  auto get_cracked_not_in_body_blockset = [&](auto ref_level_meshset) {
    Range ref_level_tet;
    CHKERR moab.get_entities_by_type(ref_level_meshset, MBTET, ref_level_tet,
                                     0);
    Range cracked_body_tets;
    CHKERR meshsetMngPtr->getEntitiesByDimension(
        crackedBodyBlockSetId, BLOCKSET, 3, cracked_body_tets, true);
    return subtract(ref_level_tet, cracked_body_tets);
  };
 
  std::size_t idx = 0;
  while (!bit_levels.back().test(idx))
    ++idx;
  int contact_lvl = idx + 1;
 
  // Remove parents, that should be no use anymore
  CHKERR reconstructMultiIndex(*m_field.get_ref_ents(),
                               RefEntity_change_parent(0));
 
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, cit)) {
    if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
      CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert %s (id: %d)\n",
                         cit->getName().c_str(), cit->getMeshsetId());
      EntityHandle cubit_meshset = cit->getMeshset();
 
      // get tet entities form back bit_level
      EntityHandle ref_level_meshset;
      CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
      CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(
          bit_levels.back(), BitRefLevel().set(), MBTET, ref_level_meshset);
      CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(
          bit_levels.back(), BitRefLevel().set(), MBPRISM, ref_level_meshset);
 
      // get faces and tets to split
      CHKERR prism_interface->getSides(
          cubit_meshset, bit_levels.back(),
          get_cracked_not_in_body_blockset(ref_level_meshset), true, verb);
      // set new bit level
      bit_levels.push_back(BitRefLevel().set(contact_lvl));
      // split faces and tets
      CHKERR prism_interface->splitSides(ref_level_meshset, bit_levels.back(),
                                         cubit_meshset, true, true, verb);
 
      // clean meshsets
      CHKERR moab.delete_entities(&ref_level_meshset, 1);
 
      // update cubit meshsets
      for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, cubit_it)) {
        EntityHandle updated_meshset = cubit_it->meshset;
        CHKERR bitRefPtr->updateMeshsetByEntitiesChildren(
            updated_meshset, bit_levels.front(), BitRefLevel().set(),
            bit_levels.back(), bit_levels.back(), updated_meshset, MBMAXTYPE,
            true);
      }
 
      if (cP.isPartitioned) {
        // FIXME check partitioned case
        ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
        Tag part_tag = pcomm->part_tag();
        Range tagged_sets;
        CHKERR m_field.get_moab().get_entities_by_type_and_tag(
            0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
            moab::Interface::UNION);
        for (Range::iterator mit = tagged_sets.begin();
             mit != tagged_sets.end(); mit++) {
          int part;
          CHKERR moab.tag_get_data(part_tag, &*mit, 1, &part);
          CHKERR bitRefPtr->updateMeshsetByEntitiesChildren(
              *mit, bit_levels.front(), BitRefLevel().set(), bit_levels.back(),
              bit_levels.back(), *mit, MBTET, true);
          Range ents3d;
          CHKERR moab.get_entities_by_dimension(*mit, 3, ents3d, true);
          for (Range::iterator eit = ents3d.begin(); eit != ents3d.end();
               eit++) {
            CHKERR moab.tag_set_data(part_tag, &*eit, 1, &part);
          }
          for (int dd = 0; dd != 3; dd++) {
            Range ents;
            CHKERR moab.get_adjacencies(ents3d, dd, false, ents,
                                        moab::Interface::UNION);
            for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
              CHKERR moab.tag_set_data(part_tag, &*eit, 1, &part);
            }
          }
        }
      }
 
      if (contact_lvl > 1) {
        CHKERR m_field.delete_ents_by_bit_ref(bit_levels.front(),
                                              bit_levels.front(), true);
      }
      BitRefLevel shift_mask;
      for (int ll = contact_lvl - 1; ll != 19; ++ll) {
        shift_mask.set(ll);
      }
      CHKERR m_field.getInterface<BitRefManager>()->shiftRightBitRef(
          1, shift_mask, verb);
      bit_levels.pop_back();
    }
  }
 
  if (meshsetMngPtr->checkMeshset(contactPrismsBlockSetId, BLOCKSET)) {
    CHKERR meshsetMngPtr->deleteMeshset(BLOCKSET, contactPrismsBlockSetId);
  }
  CHKERR meshsetMngPtr->addMeshset(BLOCKSET, contactPrismsBlockSetId);
 
  Range prisms_level;
  CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(
      bit_levels.back(), BitRefLevel().set(), MBPRISM, prisms_level);
 
  CHKERR meshsetMngPtr->addEntitiesToMeshset(BLOCKSET, contactPrismsBlockSetId,
                                             prisms_level);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::findContactFromPrisms(const BitRefLevel bit,
                                                const int verb,
                                                const bool debug) {
  MoFEM::Interface &m_field = cP.mField;
  moab::Interface &moab = cP.mField.get_moab();
  MoFEMFunctionBegin;
 
  cP.contactElements.clear();
  cP.contactSlaveFaces.clear();
  cP.contactMasterFaces.clear();
 
  CHKERR meshsetMngPtr->getEntitiesByDimension(
      contactPrismsBlockSetId, BLOCKSET, 3, cP.contactElements, true);
  CHKERR bitRefPtr->filterEntitiesByRefLevel(bit, BitRefLevel().set(),
                                             cP.contactElements);
 
  if (!cP.mortarContactElements.empty())
    cP.contactElements = subtract(cP.contactElements, cP.mortarContactElements);
 
  if (!cP.contactElements.empty()) {
    EntityHandle tri;
    for (auto p : cP.contactElements) {
      CHKERR moab.side_element(p, 2, 3, tri);
      cP.contactMasterFaces.insert(tri);
      CHKERR moab.side_element(p, 2, 4, tri);
      cP.contactSlaveFaces.insert(tri);
    }
  }
 
  Range slave_nodes, master_nodes;
  CHKERR moab.get_connectivity(cP.contactSlaveFaces, slave_nodes, false);
  CHKERR moab.get_connectivity(cP.contactMasterFaces, master_nodes, false);
 
  if (slave_nodes.size() < master_nodes.size()) {
    cP.contactBothSidesSlaveFaces = cP.contactSlaveFaces;
    cP.contactBothSidesMasterFaces = cP.contactMasterFaces;
  } else {
    cP.contactBothSidesSlaveFaces = cP.contactMasterFaces;
    cP.contactBothSidesMasterFaces = cP.contactSlaveFaces;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::refineCrackTip(const int front_id, const int verb,
                                         const bool debug) {
  MoFEM::Interface &m_field = cP.mField;
  moab::Interface &moab = m_field.get_moab();
  MeshRefinement *mesh_refiner_ptr;
  MoFEMFunctionBegin;
  CHKERR m_field.getInterface(mesh_refiner_ptr);
 
  BitRefLevel bit = cP.mapBitLevel["mesh_cut"];
 
  std::size_t idx = 0;
  while (!bit.test(idx))
    ++idx;
  BitRefLevel bit_ref = BitRefLevel().set(idx + 1);
 
  // loop to refine at crack tip
  for (int ll = 0; ll != cP.refAtCrackTip; ll++) {
 
    // get tets at bit level
    Range tets_level; // test at level
    CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                   MBTET, tets_level);
    Range edges_level; // edges at level
    CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                   MBEDGE, edges_level);
    if (!meshsetMngPtr->checkMeshset(front_id, SIDESET)) {
      PetscPrintf(PETSC_COMM_WORLD, "Sideset %d with crack front not found\n",
                  cuttingSurfaceSidesetId);
      MoFEMFunctionReturnHot(0);
    }
    Range crack_front; // crack front edges at level
    CHKERR meshsetMngPtr->getEntitiesByDimension(front_id, SIDESET, 1,
                                                 crack_front, true);
    crack_front = intersect(crack_front, edges_level);
    Range crack_front_nodes; // nodes at crack front
    CHKERR moab.get_connectivity(crack_front, crack_front_nodes, true);
 
    Range crack_front_nodes_tets; // tets adjacent to crack front nodes
    CHKERR m_field.get_moab().get_adjacencies(crack_front_nodes, 3, false,
                                              crack_front_nodes_tets,
                                              moab::Interface::UNION);
    crack_front_nodes_tets = intersect(crack_front_nodes_tets, tets_level);
    Range edges_to_refine; // edges which going to be refined
    CHKERR m_field.get_moab().get_adjacencies(crack_front_nodes_tets, 1, false,
                                              edges_to_refine,
                                              moab::Interface::UNION);
    edges_to_refine = intersect(edges_to_refine, edges_level);
    if (edges_to_refine.empty())
      continue;
 
    // create vertices in middle of refined edges
    CHKERR mesh_refiner_ptr->addVerticesInTheMiddleOfEdges(
        edges_to_refine, bit_ref, VERBOSE);
    // refine tets
    CHKERR mesh_refiner_ptr->refineTets(tets_level, bit_ref, false);
 
    // update meshsets by new entities
    // Meshsets are updated such that: if parent is in the meshset then its
    // chiled is in the meshset to.
    for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, cubit_it)) {
      EntityHandle cubit_meshset = cubit_it->meshset;
      CHKERR bitRefPtr->updateMeshsetByEntitiesChildren(
          cubit_meshset, bit, BitRefLevel().set(), bit_ref, bit_ref,
          cubit_meshset, MBMAXTYPE, true);
    }
 
    // update meshsets with partition mesh
    if (cP.isPartitioned) {
      ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
      Tag part_tag = pcomm->part_tag();
      Range tagged_sets;
      CHKERR m_field.get_moab().get_entities_by_type_and_tag(
          0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
          moab::Interface::UNION);
      for (Range::iterator mit = tagged_sets.begin(); mit != tagged_sets.end();
           mit++) {
        int part;
        CHKERR moab.tag_get_data(part_tag, &*mit, 1, &part);
        CHKERR bitRefPtr->updateMeshsetByEntitiesChildren(
            *mit, bit, BitRefLevel().set(), bit_ref, bit_ref, *mit, MBTET,
            true);
        Range ents3d;
        CHKERR moab.get_entities_by_dimension(*mit, 3, ents3d, true);
        for (Range::iterator eit = ents3d.begin(); eit != ents3d.end(); eit++) {
          CHKERR moab.tag_set_data(part_tag, &*eit, 1, &part);
        }
        for (int dd = 0; dd != 3; dd++) {
          Range ents;
          CHKERR moab.get_adjacencies(ents3d, dd, false, ents,
                                      moab::Interface::UNION);
          for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
            CHKERR moab.tag_set_data(part_tag, &*eit, 1, &part);
          }
        }
      }
    }
 
    // shift bits, if ents in on bit level 0, and only on that bit it is deleted
    BitRefLevel shift_mask;
    for (int ll = 0; ll != 19; ++ll)
      shift_mask.set(ll);
    CHKERR bitRefPtr->shiftRightBitRef(1, shift_mask);
 
    if (debug) {
      if (m_field.get_comm_rank() == 0) {
        std::string name;
        name = "out_crack_surface_level_" +
               boost::lexical_cast<std::string>(ll) + ".vtk";
        Range crack_faces;
        CHKERR meshsetMngPtr->getEntitiesByDimension(crackSurfaceId, SIDESET, 2,
                                                     crack_faces, true);
        EntityHandle meshset;
        CHKERR m_field.get_moab().create_meshset(MESHSET_SET, meshset);
        CHKERR m_field.get_moab().add_entities(meshset, crack_faces);
        CHKERR m_field.get_moab().write_file(name.c_str(), "VTK", "", &meshset,
                                             1);
        CHKERR m_field.get_moab().delete_entities(&meshset, 1);
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::refineAndSplit(const int verb, const bool debug) {
  MoFEM::Interface &m_field = cP.mField;
  MoFEMFunctionBegin;
 
  // clear tags data created by previous use of prism interface
  Tag th_interface_side;
  rval = m_field.get_moab().tag_get_handle("INTERFACE_SIDE", th_interface_side);
  if (rval == MB_SUCCESS)
    CHKERR m_field.get_moab().tag_delete(th_interface_side);
 
  auto delete_mofem_meshsets = [&]() {
    MoFEMFunctionBegin;
    CHKERR cP.mField.getInterface<MeshsetsManager>()->deleteMeshset(
        SIDESET, getSkinOfTheBodyId(), MF_ZERO);
    CHKERR cP.mField.getInterface<MeshsetsManager>()->deleteMeshset(
        SIDESET, getCrackFrontId(), MF_ZERO);
    MoFEMFunctionReturn(0);
  };
 
  BitRefLevel bit0 = cP.mapBitLevel["mesh_cut"];
  BitRefLevel bit1 = cP.mapBitLevel["spatial_domain"];
 
  CHKERR delete_mofem_meshsets();
  // find body skin
  CHKERR findBodySkin(bit0, verb, debug, "out_body_skin_bit0.vtk");
  // find/get crack surface
  CHKERR findCrack(bit0, verb, debug);
 
  // refine at crack trip
  if (cP.refAtCrackTip > 0) {
    if (cP.doCutMesh) {
      SETERRQ(
          cP.mField.get_comm(), MOFEM_IMPOSSIBLE_CASE,
          "Refinement at the crack front after the crack insertion ('my_ref') "
          "should not be used when the mesh cutting is on, use "
          "'ref_before_cut' or 'ref_before_trim' instead");
    }
    CHKERR refineCrackTip(crackFrontId, verb, debug);
    CHKERR delete_mofem_meshsets();
    CHKERR findBodySkin(bit0, verb, debug, "out_body_skin_bit0_ref.vtk");
    CHKERR findCrack(bit0, verb, debug);
  }
 
  if (debug)
    CHKERR bitRefPtr->writeBitLevelByType(bit0, BitRefLevel().set(), MBTET,
                                          "out_mesh_after_refine.vtk", "VTK",
                                          "");
 
  if (!cP.ignoreContact) {
 
    if (debug) {
      EntityHandle meshset_interface;
      CHKERR meshsetMngPtr->getMeshset(crackSurfaceId, SIDESET,
                                       meshset_interface);
      CHKERR cP.mField.get_moab().write_file(
          "crack_surface_before_contact_split.vtk", "VTK", "",
          &meshset_interface, 1);
    }
 
    // split faces to insert contact elements
    CHKERR insertContactInterface(verb, debug);
    if (debug) {
 
      EntityHandle meshset_interface;
      CHKERR meshsetMngPtr->getMeshset(crackSurfaceId, SIDESET,
                                       meshset_interface);
      CHKERR cP.mField.get_moab().write_file(
          "crack_surface_after_contact_split.vtk", "VTK", "",
          &meshset_interface, 1);
 
      CHKERR bitRefPtr->writeBitLevelByType(
          cP.mapBitLevel["mesh_cut"], BitRefLevel().set(), MBTET,
          "out_mesh_after_contact_split.vtk", "VTK", "");
      CHKERR bitRefPtr->writeBitLevelByType(
          cP.mapBitLevel["mesh_cut"], BitRefLevel().set(), MBPRISM,
          "out_mesh_contact_prisms.vtk", "VTK", "");
    }
  }
 
  // split faces to create crack
  CHKERR splitFaces(verb, debug);
 
  if (!cP.ignoreContact) {
    findContactFromPrisms(bit1, verb, debug);
  }
 
  if (!cP.ignoreContact) {
    CHKERR addMortarContactPrisms(verb, debug);
  }
 
  if (debug) {
    CHKERR bitRefPtr->writeBitLevelByType(bit1, BitRefLevel().set(), MBTET,
                                          "out_mesh_after_split.vtk", "VTK",
                                          "");
    CHKERR bitRefPtr->writeBitLevelByType(bit1, BitRefLevel().set(), MBPRISM,
                                          "out_mesh_after_split_prisms.vtk",
                                          "VTK", "");
  }
 
  // get information about topology
  CHKERR delete_mofem_meshsets();
 
  CHKERR setFixEdgesAndCorners(bit1);
  CHKERR findBodySkin(bit1, verb, debug, "out_body_skin_bit1.vtk");
  CHKERR findCrackFromPrisms(bit1, verb, debug);
  CHKERR getFrontEdgesAndElements(bit1, verb, debug);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::cutRefineAndSplit(const int verb, bool debug) {
  MoFEM::Interface &m_field = cP.mField;
  moab::Interface &moab = m_field.get_moab();
  if (debugCut)
    debug = true;
  MoFEMFunctionBegin;
 
  // FIXME: If will be needed in future one can cut, and split and then put
  // meshes on the right bit ref level.
 
  // Get entities not in database
  Range all_free_ents;
  CHKERR bitRefPtr->getAllEntitiesNotInDatabase(all_free_ents);
 
  // Cut mesh
  BitRefLevel bit = cP.mapBitLevel["mesh_cut"];
  std::size_t idx = 0;
  while (!bit.test(idx))
    ++idx;
  int first_bit = idx + 1;
 
  CHKERR cutMesh(first_bit, debug);
 
  // Squash bits
  BitRefLevel shift_mask;
  for (int ll = 0; ll != first_bit + 1; ++ll)
    shift_mask.set(ll);
  CHKERR bitRefPtr->shiftRightBitRef(first_bit, shift_mask, VERBOSE,
                                     MF_NOT_THROW);
 
  if (debug)
    CHKERR bitRefPtr->writeBitLevelByType(cP.mapBitLevel["mesh_cut"],
                                          BitRefLevel().set(), MBTET,
                                          "out_mesh_after_cut.vtk", "VTK", "");
 
  CHKERR refineAndSplit(verb, debug);
 
  if (debug)
    CHKERR bitRefPtr->writeBitLevelByType(
        cP.mapBitLevel["spatial_domain"], BitRefLevel().set(), MBTET,
        "out_mesh_after_split_and_refine.vtk", "VTK", "");
 
  // Delete left obsolete entities
 
  auto get_entities_to_delete = [&]() {
    Range ents_not_in_database;
    CHKERR moab.get_entities_by_handle(0, ents_not_in_database, false);
    ents_not_in_database = subtract(
        ents_not_in_database, ents_not_in_database.subset_by_type(MBENTITYSET));
    CHKERR bitRefPtr->filterEntitiesNotInDatabase(ents_not_in_database);
    return subtract(ents_not_in_database, all_free_ents);
  };
  Range ents_to_remove = get_entities_to_delete();
 
  auto remove_entities_from_meshsets = [&]() {
    MoFEMFunctionBegin;
    // Remove entities from meshsets
    Range meshsets;
    CHKERR moab.get_entities_by_type(0, MBENTITYSET, meshsets, false);
    for (auto m : meshsets)
      CHKERR moab.remove_entities(m, ents_to_remove);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR remove_entities_from_meshsets();
  CHKERR moab.delete_entities(ents_to_remove);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::getFrontEdgesAndElements(const BitRefLevel bit,
                                                   const int verb,
                                                   const bool debug) {
  MoFEM::Interface &m_field = cP.mField;
  moab::Interface &moab = m_field.get_moab();
  MoFEMFunctionBegin;
 
  Range level_edges;
  CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                 MBEDGE, level_edges);
 
  cP.crackFrontNodes.clear();
  CHKERR moab.get_connectivity(cP.crackFront, cP.crackFrontNodes, true);
  cP.crackFrontNodesEdges.clear();
  CHKERR moab.get_adjacencies(cP.crackFrontNodes, 1, false,
                              cP.crackFrontNodesEdges, moab::Interface::UNION);
  cP.crackFrontNodesEdges = subtract(cP.crackFrontNodesEdges, cP.crackFront);
  cP.crackFrontNodesEdges = intersect(cP.crackFrontNodesEdges, level_edges);
 
  Range level_tets;
  CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                 MBTET, level_tets);
 
  cP.crackFrontElements.clear();
  CHKERR moab.get_adjacencies(cP.crackFrontNodes, 3, false,
                              cP.crackFrontElements, moab::Interface::UNION);
  cP.crackFrontElements = intersect(cP.crackFrontElements, level_tets);
 
  if (debug) {
    EntityHandle out_meshset;
    CHKERR cP.mField.get_moab().create_meshset(MESHSET_SET, out_meshset);
    CHKERR cP.mField.get_moab().add_entities(out_meshset, cP.crackFront);
    CHKERR cP.mField.get_moab().add_entities(out_meshset, cP.crackFrontNodes);
    CHKERR cP.mField.get_moab().write_file("crack_front_nodes.vtk", "VTK", "",
                                           &out_meshset, 1);
    CHKERR cP.mField.get_moab().delete_entities(&out_meshset, 1);
  }
 
  if (verb >= VERY_VERBOSE) {
    cerr << "getFrontEdgesAndElements ";
    cerr << "crackFront\n" << cP.crackFront << endl;
    cerr << "crackFrontNodes\n" << cP.crackFrontNodes << endl;
    cerr << cP.crackFront.size() << " ";
    cerr << cP.crackFrontNodes.size() << " ";
    cerr << cP.crackFrontNodesEdges.size() << " ";
    cerr << cP.crackFrontElements.size() << endl;
  }
 
  if (verb >= NOISY) {
    cerr << cP.crackFrontNodes << endl;
    cerr << cP.crackFrontNodesEdges << endl;
    cerr << cP.crackFrontElements << endl;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::setMeshOrgCoords() {
  MoFEMFunctionBegin;
 
  BitRefLevel mask;
  for (int ll = 20; ll != BITREFLEVEL_SIZE; ++ll)
    mask.set(ll);
 
  // Get original mesh
  Range org_verts;
  CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(mask, BitRefLevel().set(),
                                                 MBVERTEX, org_verts);
 
  // Get original positions from tag
  Tag th;
  rval = cP.mField.get_moab().tag_get_handle("ORG_POSITION", th);
  if (rval == MB_SUCCESS) {
    originalCoords.resize(org_verts.size() * 3);
    CHKERR cP.mField.get_moab().tag_get_data(th, org_verts,
                                             &*originalCoords.begin());
    // Set cords
    CHKERR cP.mField.get_moab().set_coords(org_verts, &*originalCoords.begin());
  } else if (rval != MB_TAG_NOT_FOUND)
    SETERRQ(cP.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Can not get tag handle");
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::getMeshOrgCoords() {
  MoFEMFunctionBegin;
 
  BitRefLevel mask;
  for (int ll = 20; ll != BITREFLEVEL_SIZE; ++ll)
    mask.set(ll);
 
  // Get node coordinates on original mesh
  Range org_verts;
  CHKERR bitRefPtr->getEntitiesByTypeAndRefLevel(mask, BitRefLevel().set(),
                                                 MBVERTEX, org_verts);
  originalCoords.resize(org_verts.size() * 3);
  CHKERR cP.mField.get_moab().get_coords(org_verts, &*originalCoords.begin());
 
  // Save positions on tag
  double def_position[] = {0, 0, 0};
  Tag th;
  CHKERR cP.mField.get_moab().tag_get_handle("ORG_POSITION", 3, MB_TYPE_DOUBLE,
                                             th, MB_TAG_CREAT | MB_TAG_SPARSE,
                                             def_position);
  CHKERR cP.mField.get_moab().tag_set_data(th, org_verts,
                                           &*originalCoords.begin());
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPMeshCut::clearData() {
  MoFEMFunctionBegin;
 
  fixedEdges.clear();
  cornerNodes.clear();
  fRont.clear();
  sUrface.clear();
  originalCoords.clear();
 
  MoFEMFunctionReturn(0);
};
 
} // namespace FractureMechanics