EshelbianFracture.cpp

Go to the documentation of this file.

/**
 * @file EshelbianFracture.cpp
 * @author your name (you@domain.com)
 * @brief
 * @version 0.1
 * @date 2025-12-21
 *
 * @copyright Copyright (c) 2025
 *
 */
 
namespace EshelbianPlasticity {
 
MoFEMErrorCode EshelbianCore::calculateFaceMaterialForce(
    const int tag, TS ts, SmartPetscObj<Vec> *adjoint_gradient_vector) {
  MoFEMFunctionBegin;
 
  constexpr bool debug = false;
 
  auto get_tags_vec = [&](std::vector<std::pair<std::string, int>> names) {
    std::vector<Tag> tags;
    tags.reserve(names.size());
    auto create_and_clean = [&]() {
      MoFEMFunctionBegin;
      for (auto n : names) {
        tags.push_back(Tag());
        auto &tag = tags.back();
        auto &moab = mField.get_moab();
        auto rval = moab.tag_get_handle(n.first.c_str(), tag);
        if (rval == MB_SUCCESS) {
          moab.tag_delete(tag);
        }
        double def_val[] = {0., 0., 0.};
        CHKERR moab.tag_get_handle(n.first.c_str(), n.second, MB_TYPE_DOUBLE,
                                   tag, MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
      }
      MoFEMFunctionReturn(0);
    };
    CHK_THROW_MESSAGE(create_and_clean(), "create_and_clean");
    return tags;
  };
 
  enum ExhangeTags {
    MATERIALFORCE,
    ADJOINT_MATERIALFORCE,
    AREAGROWTH,
    GRIFFITHFORCE,
    ADJOINT_GRIFFITHFORCE,
    FACEPRESSURE
  };
 
  auto tags = get_tags_vec({{"MaterialForce", 3},
                            {"AdjointMaterialForce", 3},
                            {"AreaGrowth", 3},
                            {"GriffithForce", 1},
                            {"AdjointGriffithForce", 1},
                            {"FacePressure", 1}});
 
  auto calculate_material_forces = [&]() {
    MoFEMFunctionBegin;
 
    /**
     * @brief Create element to integration faces energies
     */
    auto get_face_material_force_fe = [&]() {
      using FaceEle = FaceElementForcesAndSourcesCore;
      auto fe_ptr = boost::make_shared<FaceEle>(mField);
      fe_ptr->getRuleHook = [](int, int, int) { return -1; };
      fe_ptr->setRuleHook =
          SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
      if (ts != PETSC_NULLPTR) {
        fe_ptr->data_ctx |= PetscData::CTX_SET_TIME;
        CHKERR TSGetTime(ts, &(fe_ptr->ts_t));
        CHKERR TSGetTimeStep(ts, &(fe_ptr->ts_dt));
      }
      // hybrid disp, evaluated on face first
      EshelbianPlasticity::AddHOOps<2, 2, 3>::add(
          fe_ptr->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
      fe_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          hybridSpatialDisp, dataAtPts->getHybridDispAtPts()));
      fe_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
              hybridSpatialDisp, dataAtPts->getGradHybridDispAtPts()));
      auto op_loop_domain_side =
          new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
              mField, elementVolumeName, SPACE_DIM, Sev::noisy);
      fe_ptr->getOpPtrVector().push_back(op_loop_domain_side);
      fe_ptr->getOpPtrVector().push_back(new OpFaceMaterialForce(dataAtPts));
 
      // evaluated in side domain, that is op_loop_domain_side
      op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
          boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
 
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
          materialH1Positions, frontAdjEdges, nullptr, nullptr, nullptr);
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
              piolaStress, dataAtPts->getApproxPAtPts()));
 
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(
              rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
      CHKERR physicalEquations->pushMaterialForceFields(
          *this, op_loop_domain_side->getOpPtrVector(), dataAtPts);
 
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpFaceSideMaterialForce(dataAtPts));
 
      return fe_ptr;
    };
 
    auto integrate_face_material_force_fe = [&](auto &&face_energy_fe) {
      MoFEMFunctionBegin;
      CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
          dM, skeletonElement, face_energy_fe, 0, mField.get_comm_size());
 
      auto face_exchange = CommInterface::createEntitiesPetscVector(
          mField.get_comm(), mField.get_moab(), 2, 3, Sev::inform);
 
      auto print_loc_size = [this](auto v, auto str, auto sev) {
        MoFEMFunctionBegin;
        int size;
        CHKERR VecGetLocalSize(v.second, &size);
        int low, high;
        CHKERR VecGetOwnershipRange(v.second, &low, &high);
        MOFEM_LOG("EPSYNC", sev) << str << " local size " << size << " ( "
                                 << low << " " << high << " ) ";
        MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), sev);
        MoFEMFunctionReturn(0);
      };
      CHKERR print_loc_size(face_exchange, "material face_exchange",
                            Sev::verbose);
 
      CHKERR CommInterface::updateEntitiesPetscVector(
          mField.get_moab(), face_exchange, tags[ExhangeTags::MATERIALFORCE]);
      CHKERR CommInterface::updateEntitiesPetscVector(
          mField.get_moab(), faceExchange, tags[ExhangeTags::FACEPRESSURE]);
 
      #ifndef NDEBUG
      if (debug) {
        CHKERR save_range(mField.get_moab(),
                          "front_skin_faces_material_force_" +
                              std::to_string(mField.get_comm_rank()) + ".vtk",
                          *skeletonFaces);
      }
      #endif
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR integrate_face_material_force_fe(get_face_material_force_fe());
 
    MoFEMFunctionReturn(0);
  };
 
  auto get_conn = [&](auto e) {
    Range conn;
    CHK_MOAB_THROW(mField.get_moab().get_connectivity(&e, 1, conn, true),
                   "get connectivity");
    return conn;
  };
 
  auto get_conn_range = [&](auto e) {
    Range conn;
    CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
                   "get connectivity");
    return conn;
  };
 
  auto get_adj = [&](auto e, auto dim) {
    Range adj;
    CHK_MOAB_THROW(mField.get_moab().get_adjacencies(&e, 1, dim, true, adj),
                   "get adj");
    return adj;
  };
 
  auto get_adj_range = [&](auto e, auto dim) {
    Range adj;
    CHK_MOAB_THROW(mField.get_moab().get_adjacencies(e, dim, true, adj,
                                                     moab::Interface::UNION),
                   "get adj");
    return adj;
  };
 
  auto get_vector_tag_data = [&](auto r, auto th) {
    MatrixDouble tag_data(r.size(), 3, false);
    CHK_MOAB_THROW(
        mField.get_moab().tag_get_data(th, r, tag_data.data().data()),
        "get vector tag data");
    return tag_data;
  };
 
  auto calculate_edge_direction = [&](auto e) {
    const EntityHandle *conn;
    int num_nodes;
    CHK_MOAB_THROW(
        mField.get_moab().get_connectivity(e, conn, num_nodes, true),
        "get connectivity");
    std::array<double, 6> coords;
    CHK_MOAB_THROW(mField.get_moab().get_coords(conn, num_nodes, coords.data()),
                   "get coords");
    FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p0{
        &coords[0], &coords[1], &coords[2]};
    FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p1{
        &coords[3], &coords[4], &coords[5]};
    FTensor::Tensor1<double, SPACE_DIM> t_dir;
    FTENSOR_INDEX(SPACE_DIM, i);
    t_dir(i) = t_p1(i) - t_p0(i);
    return t_dir;
  };
 
  auto average_vector_tag_at_edge = [&](auto th) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEXES(SPACE_DIM, I, J, K);
 
    for (auto e : *frontEdges) {
      auto conn = get_conn(e);
      auto data = get_vector_tag_data(conn, th);
      auto t_node = getFTensor1FromPtr<SPACE_DIM>(data.data().data());
      FTensor::Tensor1<double, SPACE_DIM> t_edge_material_force{0., 0., 0.};
      for (auto n : conn) {
        NOT_USED(n);
        t_edge_material_force(I) += t_node(I);
        ++t_node;
      }
      t_edge_material_force(I) /= conn.size();
 
      FTensor::Tensor1<double, SPACE_DIM> t_edge_direction =
          calculate_edge_direction(e);
      t_edge_direction.normalize();
 
      // Project the averaged vector to the plane normal to the front edge.
      FTensor::Tensor1<double, SPACE_DIM> t_cross;
      t_cross(K) = FTensor::levi_civita(I, J, K) * t_edge_direction(I) *
                   t_edge_material_force(J);
      t_edge_material_force(K) =
          FTensor::levi_civita(I, J, K) * t_edge_direction(J) * t_cross(I);
 
      CHKERR mField.get_moab().tag_set_data(th, &e, 1,
                                            &t_edge_material_force(0));
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto average_material_force_at_edge = [&](auto th) {
    MoFEMFunctionBegin;
 
    if (mField.get_comm_rank() == 0) {
      CHKERR average_vector_tag_at_edge(th);
 
// #ifndef NDEBUG
//       if (debug) {
        int ts_step;
        CHKERR TSGetStepNumber(ts, &ts_step);
        CHKERR save_range(mField.get_moab(),
                          "front_edges_material_force_" +
                              std::to_string(ts_step) + ".vtk",
                          *frontEdges);
//       }
// #endif
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_force_through_node = [&](auto nb_J_integral_contours) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEX(SPACE_DIM, I);
 
    if (mField.get_comm_rank() == 0) {
      auto front_nodes = get_conn_range(*frontEdges);
      Range all_skin_faces;
 
      for (auto n : front_nodes) {
        auto adj_tets = get_adj(n, SPACE_DIM);
        for (int ll = 0; ll < nb_J_integral_contours; ++ll) {
          auto conn = get_conn_range(adj_tets);
          adj_tets = get_adj_range(conn, SPACE_DIM);
        }
 
        auto skin_faces = get_skin(mField, adj_tets);
        auto material_forces =
            get_vector_tag_data(skin_faces, tags[ExhangeTags::MATERIALFORCE]);
 
#ifndef NDEBUG
        if (debug) {
          all_skin_faces.merge(skin_faces);
        }
#endif
 
        auto t_face_T =
            getFTensor1FromPtr<SPACE_DIM>(material_forces.data().data());
        FTensor::Tensor1<double, SPACE_DIM> t_node_force{0., 0., 0.};
        for (auto face : skin_faces) {
 
          FTensor::Tensor1<double, SPACE_DIM> t_face_force_tmp{0., 0., 0.};
          t_face_force_tmp(I) = t_face_T(I);
          ++t_face_T;
 
          auto face_tets = intersect(get_adj(face, SPACE_DIM), adj_tets);
 
          if (face_tets.empty()) {
            continue;
          }
 
          if (face_tets.size() != 1) {
            CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                              "face_tets.size() != 1");
          }
 
          int side_number, sense, offset;
          CHK_MOAB_THROW(mField.get_moab().side_number(face_tets[0], face,
                                                       side_number, sense,
                                                       offset),
                         "moab side number");
          t_face_force_tmp(I) *= sense;
          t_node_force(I) += t_face_force_tmp(I);
        }
 
        t_node_force(I) /= griffithEnergy; // scale all by griffith energy
        CHK_MOAB_THROW(
            mField.get_moab().tag_set_data(tags[ExhangeTags::MATERIALFORCE],
                                           &n, 1, &t_node_force(0)),
            "set data");
      }
 
#ifndef NDEBUG
      if (debug) {
        int ts_step;
        CHKERR TSGetStepNumber(ts, &ts_step);
        CHKERR save_range(mField.get_moab(),
                          "front_skin_faces_material_force_" +
                              std::to_string(ts_step) + ".vtk",
                          all_skin_faces);
      }
#endif
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto get_adj_tets_for_contour = [&](auto n, auto nb_J_integral_contours) {
    auto adj_tets = get_adj(n, SPACE_DIM);
    for (int ll = 0; ll < nb_J_integral_contours; ++ll) {
      auto conn = get_conn_range(adj_tets);
      adj_tets = get_adj_range(conn, SPACE_DIM);
    }
    return adj_tets;
  };
 
  auto get_front_node_adj_crack_faces = [&](auto n) {
    return intersect(get_adj(n, SPACE_DIM - 1), *crackFaces);
  };
 
  auto calculate_crack_area_growth_face = [&](auto nb_J_integral_contours) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEXES(SPACE_DIM, I, J, K, L);
 
    if (mField.get_comm_rank() == 0) {
      auto front_nodes = get_conn_range(*frontEdges);
      auto body_edges = get_range_from_block(mField, "EDGES", 1);
      Range body_ents;
      CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                         body_ents);
      auto body_skin = get_skin(mField, body_ents);
      auto body_skin_conn = get_conn_range(body_skin);
 
      auto calculate_seed_area_growth = [&](auto n, auto &adj_faces) {
        // if skin is on body surface, project the direction on it
        FTensor::Tensor1<double, SPACE_DIM> t_project{0., 0., 0.};
        auto boundary_node = intersect(Range(n, n), body_skin_conn);
        if (boundary_node.size()) {
          auto faces = intersect(get_adj(n, SPACE_DIM - 1), body_skin);
          for (auto f : faces) {
            FTensor::Tensor1<double, 3> t_normal_face;
            CHKERR mField.getInterface<Tools>()->getTriNormal(
                f, &t_normal_face(0));
            t_project(I) += t_normal_face(I);
          }
          t_project.normalize();
        }
 
        // calculate surface projection matrix
        FTensor::Tensor2<double, 3, 3> t_Q;
        auto t_kd = FTensor::Kronecker_Delta<double>();
        t_Q(I, J) = t_kd(I, J);
        if (boundary_node.size()) {
          t_Q(I, J) -= t_project(I) * t_project(J);
        }
 
        FTensor::Tensor1<double, 3> t_area_dir{0., 0., 0.};
        for (auto f : adj_faces) {
          int num_nodes;
          const EntityHandle *conn;
          CHKERR mField.get_moab().get_connectivity(f, conn, num_nodes, true);
          std::array<double, 9> coords;
          CHKERR mField.get_moab().get_coords(conn, num_nodes, coords.data());
          FTensor::Tensor1<double, 3> t_face_normal;
          FTensor::Tensor3<double, 3, 3, 3> t_d_normal;
          CHKERR mField.getInterface<Tools>()->getTriNormal(
              coords.data(), &t_face_normal(0), &t_d_normal(0, 0, 0));
          auto n_it = std::find(conn, conn + num_nodes, n);
          auto n_index = std::distance(conn, n_it);
 
          FTensor::Tensor2<double, 3, 3> t_face_hessian{
              t_d_normal(0, n_index, 0), t_d_normal(0, n_index, 1),
              t_d_normal(0, n_index, 2),
 
              t_d_normal(1, n_index, 0), t_d_normal(1, n_index, 1),
              t_d_normal(1, n_index, 2),
 
              t_d_normal(2, n_index, 0), t_d_normal(2, n_index, 1),
              t_d_normal(2, n_index, 2)};
 
          FTensor::Tensor2<double, 3, 3> t_projected_hessian;
          t_projected_hessian(I, J) =
              t_Q(I, K) * (t_face_hessian(K, L) * t_Q(L, J));
          t_face_normal.normalize();
          t_area_dir(K) += t_face_normal(I) * t_projected_hessian(I, K) / 2.;
        }
 
        return t_area_dir;
      };
 
      auto get_crack_area_growth_seed_nodes = [&](auto &adj_tets) {
        // This gets all 1D edges adjacent to the current tetrahedral patch
        // adj_tets, then keeps only edges that are either crack-front edges or
        // special body edges from the "EDGES" block. So adj_edges is the local
        // edge stencil relevant to crack growth near the current front node.
        auto adj_edges = intersect(get_adj_range(adj_tets, 1),
                                   unite(*frontEdges, body_edges));
 
        // This collects all vertices/nodes connected to those edges. These are
        // candidate seed nodes whose local crack-area growth contribution will
        // be accumulated.
        auto seed_n = get_conn_range(adj_edges);
        auto skin_adj_edges = get_skin(mField, adj_edges);
        skin_adj_edges = subtract(skin_adj_edges, body_skin_conn);
        seed_n = subtract(seed_n, skin_adj_edges);
 
        return std::make_pair(seed_n, skin_adj_edges);
      };
 
      auto calculate_front_node_area_growth = [&](auto &adj_tets) {
        auto [seed_n, skin_adj_edges] =
            get_crack_area_growth_seed_nodes(adj_tets);
 
        FTensor::Tensor1<double, SPACE_DIM> t_area_dir{0., 0., 0.};
        auto add_area_growth_direction = [&](auto sn, double weight) {
          auto adj_faces = intersect(get_adj(sn, SPACE_DIM - 1), *crackFaces);
          if (adj_faces.empty()) {
            return;
          }
 
          auto t_area_dir_sn = calculate_seed_area_growth(sn, adj_faces);
          t_area_dir(I) += weight * t_area_dir_sn(I);
        };
 
        for (auto sn : seed_n) {
          add_area_growth_direction(sn, 1.);
        }
        for (auto sn : skin_adj_edges) {
          add_area_growth_direction(sn, 0.5);
        }
 
        return t_area_dir;
      };
 
      for (auto n : front_nodes) {
        auto front_node_adj_faces = get_front_node_adj_crack_faces(n);
        if (front_node_adj_faces.empty()) {
          continue;
        }
 
        auto adj_tets = get_adj_tets_for_contour(n, nb_J_integral_contours);
        auto t_area_dir = calculate_front_node_area_growth(adj_tets);
 
        CHK_MOAB_THROW(
            mField.get_moab().tag_set_data(tags[ExhangeTags::AREAGROWTH], &n, 1,
                                           &t_area_dir(0)),
            "set data");
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_crack_area_growth_no_face = [&](auto nb_J_integral_contours,
                                                 auto material_force_tag) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEX(SPACE_DIM, I);
 
    if (mField.get_comm_rank() == 0) {
      auto front_nodes = get_conn_range(*frontEdges);
      auto body_edges = get_range_from_block(mField, "EDGES", 1);
      Range body_ents;
      CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                         body_ents);
      auto body_skin = get_skin(mField, body_ents);
      auto body_skin_conn = get_conn_range(body_skin);
 
      auto calculate_seed_area_growth = [&](auto n, auto &t_node_force) {
        auto adj_edges =
            intersect(get_adj(n, 1), unite(*frontEdges, body_edges));
        double l = 0;
        for (auto e : adj_edges) {
          auto t_dir = calculate_edge_direction(e);
          l += t_dir.l2();
        }
        l /= 2;
 
        FTensor::Tensor1<double, SPACE_DIM> t_area_dir{0., 0., 0.};
        FTensor::Tensor1<double, SPACE_DIM> t_node_force_tmp;
        t_node_force_tmp(I) = t_node_force(I);
        t_node_force_tmp.normalize();
        t_area_dir(I) = -t_node_force_tmp(I);
        t_area_dir(I) *= l / 2;
        return t_area_dir;
      };
 
      auto get_crack_area_growth_seed_nodes = [&](auto &adj_tets) {
        auto adj_edges = intersect(get_adj_range(adj_tets, 1),
                                   unite(*frontEdges, body_edges));
        auto seed_n = get_conn_range(adj_edges);
        auto skin_adj_edges = get_skin(mField, adj_edges);
        skin_adj_edges = subtract(skin_adj_edges, body_skin_conn);
        seed_n = subtract(seed_n, skin_adj_edges);
 
        return std::make_pair(seed_n, skin_adj_edges);
      };
 
      auto calculate_front_node_area_growth = [&](auto &adj_tets,
                                                  auto &t_node_force) {
        auto [seed_n, skin_adj_edges] =
            get_crack_area_growth_seed_nodes(adj_tets);
 
        FTensor::Tensor1<double, SPACE_DIM> t_area_dir{0., 0., 0.};
        auto add_area_growth_direction = [&](auto sn, double weight) {
          auto t_area_dir_sn = calculate_seed_area_growth(sn, t_node_force);
          t_area_dir(I) += weight * t_area_dir_sn(I);
        };
 
        for (auto sn : seed_n) {
          add_area_growth_direction(sn, 1.);
        }
        for (auto sn : skin_adj_edges) {
          add_area_growth_direction(sn, 0.5);
        }
 
        return t_area_dir;
      };
 
      for (auto n : front_nodes) {
        auto front_node_adj_faces = get_front_node_adj_crack_faces(n);
        if (front_node_adj_faces.empty()) {
          FTensor::Tensor1<double, SPACE_DIM> t_node_force;
          CHKERR mField.get_moab().tag_get_data(tags[material_force_tag], &n, 1,
                                                &t_node_force(0));
 
          auto adj_tets = get_adj_tets_for_contour(n, nb_J_integral_contours);
          auto t_area_dir =
              calculate_front_node_area_growth(adj_tets, t_node_force);
 
          CHK_MOAB_THROW(
              mField.get_moab().tag_set_data(tags[ExhangeTags::AREAGROWTH], &n,
                                             1, &t_area_dir(0)),
              "set data");
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto update_crack_area_growth_edges = [&]() {
    MoFEMFunctionBegin;
 
    if (mField.get_comm_rank() == 0) {
      CHKERR average_vector_tag_at_edge(tags[ExhangeTags::AREAGROWTH]);
    }
 
    auto area_growth_edge_exchange = CommInterface::createEntitiesPetscVector(
        mField.get_comm(), mField.get_moab(), 1, 3, Sev::inform);
    CHKERR CommInterface::updateEntitiesPetscVector(
        mField.get_moab(), area_growth_edge_exchange,
        tags[ExhangeTags::AREAGROWTH]);
 
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_griffith_force = [&](ExhangeTags material_force_tag,
                                      ExhangeTags griffith_force_tag) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEXES(SPACE_DIM, I, J, K);
 
    if (mField.get_comm_rank() == 0) {
      auto front_nodes = get_conn_range(*frontEdges);
      Range all_front_faces;
 
      for (auto n : front_nodes) {
        FTensor::Tensor1<double, SPACE_DIM> t_node_force;
        CHKERR mField.get_moab().tag_get_data(tags[material_force_tag], &n, 1,
                                              &t_node_force(0));
        FTensor::Tensor1<double, SPACE_DIM> t_area_dir;
        CHKERR mField.get_moab().tag_get_data(tags[ExhangeTags::AREAGROWTH], &n,
                                              1, &t_area_dir(0));
 
        auto griffith =
            -t_node_force(I) * t_area_dir(I) / (t_area_dir(K) * t_area_dir(K));
        CHK_MOAB_THROW(mField.get_moab().tag_set_data(
                           tags[griffith_force_tag], &n, 1, &griffith),
                       "set data");
      }
 
      for (auto e : *frontEdges) {
        FTensor::Tensor1<double, SPACE_DIM> t_edge_force;
        CHKERR mField.get_moab().tag_get_data(tags[material_force_tag], &e, 1,
                                              &t_edge_force(0));
        FTensor::Tensor1<double, SPACE_DIM> t_edge_area_dir;
        CHKERR mField.get_moab().tag_get_data(tags[ExhangeTags::AREAGROWTH], &e,
                                              1, &t_edge_area_dir(0));
        double griffith_energy =
            -t_edge_force(I) * t_edge_area_dir(I) /
            (t_edge_area_dir(K) * t_edge_area_dir(K));
        CHKERR mField.get_moab().tag_set_data(tags[griffith_force_tag], &e, 1,
                                              &griffith_energy);
      }
 
      for (auto e : *frontEdges) {
        auto adj_faces = get_adj(e, SPACE_DIM - 1);
 
        if (debug) {
          all_front_faces.merge(adj_faces);
        }
 
        FTensor::Tensor1<double, SPACE_DIM> t_edge_force;
        CHKERR mField.get_moab().tag_get_data(tags[material_force_tag], &e, 1,
                                              &t_edge_force(0));
        FTensor::Tensor1<double, SPACE_DIM> t_edge_direction =
            calculate_edge_direction(e);
        t_edge_direction.normalize();
 
        FTensor::Tensor1<double, SPACE_DIM> t_cross;
        t_cross(K) = FTensor::levi_civita(I, J, K) * t_edge_direction(I) *
                     t_edge_force(J);
 
        for (auto f : adj_faces) {
          FTensor::Tensor1<double, SPACE_DIM> t_normal;
          CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
          t_normal.normalize();
          int side_number, sense, offset;
          CHKERR mField.get_moab().side_number(f, e, side_number, sense, offset);
          auto dot = -sense * t_cross(I) * t_normal(I);
          CHK_MOAB_THROW(mField.get_moab().tag_set_data(
                             tags[griffith_force_tag], &f, 1, &dot),
                         "set data");
        }
      }
 
#ifndef NDEBUG
      if (debug) {
        int ts_step;
        CHKERR TSGetStepNumber(ts, &ts_step);
        CHKERR save_range(mField.get_moab(),
                          "front_faces_material_force_" +
                              std::to_string(ts_step) + ".vtk",
                          all_front_faces);
      }
#endif
    }
 
    auto vector_edge_exchange = CommInterface::createEntitiesPetscVector(
        mField.get_comm(), mField.get_moab(), 1, 3, Sev::inform);
    CHKERR CommInterface::updateEntitiesPetscVector(
        mField.get_moab(), vector_edge_exchange, tags[material_force_tag]);
    auto &scalar_edge_exchange = edgeExchange;
    CHKERR CommInterface::updateEntitiesPetscVector(
        mField.get_moab(), scalar_edge_exchange, tags[griffith_force_tag]);
 
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_griffith_force_simplified = [&](auto material_force_tag,
                                                 auto griffith_force_tag) {
    MoFEMFunctionBegin;
 
    if (mField.get_comm_rank() == 0) {
      auto front_nodes = get_conn_range(*frontEdges);
 
      for (auto n : front_nodes) {
        FTensor::Tensor1<double, SPACE_DIM> t_node_force;
        CHKERR mField.get_moab().tag_get_data(tags[material_force_tag], &n, 1,
                                              &t_node_force(0));
 
        auto adj_edges = intersect(get_adj(n, 1), *frontEdges);
        double adj_edges_length = 0.;
        for (auto e : adj_edges) {
          auto t_edge_dir = calculate_edge_direction(e);
          adj_edges_length += t_edge_dir.l2();
        }
 
        const double nodal_front_length = adj_edges_length / 2.;
        if (nodal_front_length <= 0.) {
          CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                            "Front node has zero adjacent front edge length");
        }
 
        double griffith_energy = t_node_force.l2() / nodal_front_length;
        CHK_MOAB_THROW(mField.get_moab().tag_set_data(tags[griffith_force_tag],
                                                      &n, 1, &griffith_energy),
                       "set data");
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_adjoint_material_force = [&]() {
    MoFEMFunctionBegin;
 
    if (ts == PETSC_NULLPTR) {
      SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "TS is required to calculate adjoint material force");
    }
 
    auto g = createDMVector(dmMaterial, RowColData::ROW);
    CHKERR VecZeroEntries(g);
 
    // Disabled until the topological objective provides a material-compatible
    // tangent for Neo-Hookean models. Keep g zero so the remaining exchange
    // path cannot consume uninitialised vector data.
    // auto topological_tao_ctx = createTopologicalTAOCtx(
    //     this, SetIntegrationAtFrontVolume(frontVertices, frontAdjEdges),
    //     SetIntegrationAtFrontFace(frontVertices, frontAdjEdges),
    //     SmartPetscObj<TS>(ts, true));
    // double obj_value = 0;
    // CHKERR evaluateGradient(topological_tao_ctx.get(), &obj_value, g,
    //                         ObjectiveModelType::HENCKY_MODEL);
 
    auto set_vertex_exchange_from_gradient = [&]() {
      MoFEMFunctionBegin;
 
      CHKERR VecZeroEntries(vertexExchange.second);
      CHKERR VecGhostUpdateBegin(vertexExchange.second, INSERT_VALUES,
                                 SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(vertexExchange.second, INSERT_VALUES,
                               SCATTER_FORWARD);
 
      auto *problem_ptr = getProblemPtr(dmMaterial);
      auto &dofs =
          problem_ptr->getNumeredRowDofsPtr()->get<Unique_mi_tag>();
      const auto field_bit = mField.get_field_bit_number(materialH1Positions);
 
      double *g_array;
      double *exchange_array;
      CHKERR VecGetArray(g, &g_array);
      CHKERR VecGetArray(vertexExchange.second, &exchange_array);
 
      auto ptr = exchange_array; // vector values are arranged as entries, That
                                 // is key idea behind vertexExchange vector.
      for (auto v : vertexExchange.first.first) {
        std::array<double, SPACE_DIM> values = {0., 0., 0.};
        auto lo =
            dofs.lower_bound(DofEntity::getLoFieldEntityUId(field_bit, v));
        auto hi =
            dofs.upper_bound(DofEntity::getHiFieldEntityUId(field_bit, v));
        for (; lo != hi; ++lo) {
          if (!(*lo)->getHasLocalIndex())
            continue;
          const auto coeff = (*lo)->getDofCoeffIdx();
          if (coeff < SPACE_DIM)
            values[coeff] = g_array[(*lo)->getPetscLocalDofIdx()];
        }
        for (int d = 0; d != SPACE_DIM; ++d, ++ptr) {
          *ptr = values[d];
        }
      }
 
      CHKERR VecRestoreArray(vertexExchange.second, &exchange_array);
      CHKERR VecRestoreArray(g, &g_array);
 
      if (adjoint_gradient_vector != nullptr) {
        (*adjoint_gradient_vector) = g;
      }
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_vertex_exchange_from_gradient();
 
    CHKERR CommInterface::setTagFromVector(
        mField.get_moab(), vertexExchange,
        tags[ExhangeTags::ADJOINT_MATERIALFORCE]);
    CHKERR CommInterface::updateEntitiesPetscVector(
        mField.get_moab(), vertexExchange,
        tags[ExhangeTags::ADJOINT_MATERIALFORCE]);
 
    MoFEMFunctionReturn(0);
  };
 
  auto print_results = [&](auto nb_J_integral_conturs, bool print_material,
                           bool print_adjoint) {
    MoFEMFunctionBegin;
 
    if (!print_material && !print_adjoint) {
      MoFEMFunctionReturnHot(0);
    }
 
    auto get_conn_range = [&](auto e) {
      Range conn;
      CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
                     "get connectivity");
      return conn;
    };
 
    auto get_tag_data = [&](auto &ents, auto tag, auto dim) {
      std::vector<double> data(ents.size() * dim);
      CHK_MOAB_THROW(mField.get_moab().tag_get_data(tag, ents, data.data()),
                     "get data");
      return data;
    };
 
    if (mField.get_comm_rank() == 0) {
      auto at_nodes = [&]() {
        MoFEMFunctionBegin;
        auto conn = get_conn_range(*frontEdges);
        std::vector<double> material_force;
        std::vector<double> adjoint_material_force;
        auto area_growth = get_tag_data(conn, tags[ExhangeTags::AREAGROWTH], 3);
        std::vector<double> griffith_force;
        std::vector<double> adjoint_griffith_force;
        if (print_material) {
          material_force =
              get_tag_data(conn, tags[ExhangeTags::MATERIALFORCE], 3);
          griffith_force =
              get_tag_data(conn, tags[ExhangeTags::GRIFFITHFORCE], 1);
        }
        if (print_adjoint) {
          adjoint_material_force =
              get_tag_data(conn, tags[ExhangeTags::ADJOINT_MATERIALFORCE], 3);
          adjoint_griffith_force =
              get_tag_data(conn, tags[ExhangeTags::ADJOINT_GRIFFITHFORCE], 1);
        }
        std::vector<double> coords(conn.size() * 3);
        CHK_MOAB_THROW(mField.get_moab().get_coords(conn, coords.data()),
                       "get coords");
        MOFEM_LOG("EPSELF", Sev::inform) << "Force results at nodes";
        MOFEM_LOG("EPSELF", Sev::inform)
              << std::left << std::setw(10) << "kind" << std::right
              << std::setw(9) << "node" << std::setw(18) << "coord_x"
            << std::setw(18) << "coord_y" << std::setw(18) << "coord_z"
            << std::setw(18) << "force_x" << std::setw(18) << "force_y"
            << std::setw(18) << "force_z" << std::setw(18) << "area_x"
            << std::setw(18) << "area_y" << std::setw(18) << "area_z"
            << std::setw(18) << "griffith" << std::setw(10) << "contour";
 
        auto print_row = [&](const char *kind, const auto &force,
                             const auto &griffith, const size_t i) {
          MOFEM_LOG("EPSELF", Sev::inform)
              << std::left << std::setw(10) << kind << std::right
              << std::setw(9) << conn[i] << std::scientific
              << std::setprecision(10) << std::setw(18) << coords[i * 3 + 0]
              << std::setw(18) << coords[i * 3 + 1] << std::setw(18)
              << coords[i * 3 + 2] << std::setw(18) << force[i * 3 + 0]
              << std::setw(18) << force[i * 3 + 1] << std::setw(18)
              << force[i * 3 + 2] << std::setw(18) << area_growth[i * 3 + 0]
              << std::setw(18) << area_growth[i * 3 + 1] << std::setw(18)
              << area_growth[i * 3 + 2] << std::setw(18) << griffith[i]
              << std::defaultfloat << std::setprecision(6) << std::setw(10)
              << nb_J_integral_conturs;
        };
 
        for (size_t i = 0; i < conn.size(); ++i) {
          if (print_material) {
            print_row("material", material_force, griffith_force, i);
          }
          if (print_adjoint) {
            print_row("adjoint", adjoint_material_force, adjoint_griffith_force,
                      i);
          }
        }
 
        MoFEMFunctionReturn(0);
      };
 
      at_nodes();
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR calculate_material_forces();
 
  PetscBool all_contours = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "",
                             "-calculate_J_integral_all_levels", &all_contours,
                             PETSC_NULLPTR); // for backward compatibility
  CHKERR PetscOptionsGetBool(
      PETSC_NULLPTR, "", "-calculate_J_integral_all_contours", &all_contours,
      PETSC_NULLPTR); // new name
 
  if (all_contours == PETSC_TRUE) {
    for (int l = 0; l < nbJIntegralContours; ++l) {
      CHKERR calculate_force_through_node(l);
      CHKERR average_material_force_at_edge(tags[ExhangeTags::MATERIALFORCE]);
      CHKERR calculate_crack_area_growth_face(l);
      CHKERR calculate_crack_area_growth_no_face(l, ExhangeTags::MATERIALFORCE);
      CHKERR update_crack_area_growth_edges();
      CHKERR calculate_griffith_force(ExhangeTags::MATERIALFORCE,
                                      ExhangeTags::GRIFFITHFORCE);
      CHKERR print_results(l, true, false);
    }
  }
 
  PetscBool has_nonzero_ts_solution = PETSC_FALSE;
 
  if (ts != PETSC_NULLPTR) {
    Vec ts_solution = PETSC_NULLPTR;
    CHKERR TSGetSolution(ts, &ts_solution);
    if (ts_solution != PETSC_NULLPTR) {
      PetscReal ts_solution_norm = 0.0;
      CHKERR VecNorm(ts_solution, NORM_2, &ts_solution_norm);
      has_nonzero_ts_solution =
          (ts_solution_norm > PETSC_MACHINE_EPSILON) ? PETSC_TRUE : PETSC_FALSE;
    }
 
    if (has_nonzero_ts_solution == PETSC_TRUE) {
      CHKERR calculate_adjoint_material_force();
    }
  }
 
  CHKERR calculate_force_through_node(nbJIntegralContours);
  CHKERR average_material_force_at_edge(tags[ExhangeTags::MATERIALFORCE]);
  CHKERR calculate_crack_area_growth_face(nbJIntegralContours);
  CHKERR calculate_crack_area_growth_no_face(nbJIntegralContours,
                                             ExhangeTags::MATERIALFORCE);
  CHKERR update_crack_area_growth_edges();
  CHKERR calculate_griffith_force(ExhangeTags::MATERIALFORCE,
                                  ExhangeTags::GRIFFITHFORCE);
  if (has_nonzero_ts_solution == PETSC_TRUE) {
    CHKERR calculate_griffith_force(ExhangeTags::ADJOINT_MATERIALFORCE,
                                    ExhangeTags::ADJOINT_GRIFFITHFORCE);
    CHKERR calculate_griffith_force_simplified(
        ExhangeTags::ADJOINT_MATERIALFORCE, ExhangeTags::ADJOINT_GRIFFITHFORCE);
  }
  CHKERR print_results(nbJIntegralContours, true, true);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::calculateOrientation(const int tag,
                                                   bool set_orientation) {
  MoFEMFunctionBegin;
 
  constexpr bool debug = false;
  (void)debug;
  constexpr auto sev = Sev::verbose;
 
  Range body_ents;
  CHKERR mField.get_moab().get_entities_by_dimension(0, 3, body_ents);
  auto body_skin = get_skin(mField, body_ents);
  Range body_skin_edges;
  CHKERR mField.get_moab().get_adjacencies(body_skin, 1, false, body_skin_edges,
                                           moab::Interface::UNION);
  Range boundary_skin_verts;
  CHKERR mField.get_moab().get_connectivity(body_skin_edges,
                                            boundary_skin_verts, true);
 
  auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
  Range geometry_edges_verts;
  CHKERR mField.get_moab().get_connectivity(geometry_edges,
                                            geometry_edges_verts, true);
  Range crack_faces_verts;
  CHKERR mField.get_moab().get_connectivity(*crackFaces, crack_faces_verts,
                                            true);
  Range crack_faces_edges;
  CHKERR mField.get_moab().get_adjacencies(
      *crackFaces, 1, true, crack_faces_edges, moab::Interface::UNION);
  Range crack_faces_tets;
  CHKERR mField.get_moab().get_adjacencies(
      *crackFaces, 3, true, crack_faces_tets, moab::Interface::UNION);
 
  Range front_verts;
  CHKERR mField.get_moab().get_connectivity(*frontEdges, front_verts, true);
  Range front_faces;
  CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2, true, front_faces,
                                           moab::Interface::UNION);
  Range front_verts_edges;
  CHKERR mField.get_moab().get_adjacencies(
      front_verts, 1, true, front_verts_edges, moab::Interface::UNION);
 
  auto get_tags_vec = [&](auto tag_name, int dim) {
    std::vector<Tag> tags(1);
 
    if (dim > 3)
      CHK_THROW_MESSAGE(MOFEM_ATOM_TEST_INVALID, "dim>3");
 
    auto create_and_clean = [&]() {
      MoFEMFunctionBegin;
      auto &moab = mField.get_moab();
      auto rval = moab.tag_get_handle(tag_name, tags[0]);
      if (rval == MB_SUCCESS) {
        moab.tag_delete(tags[0]);
      }
      double def_val[] = {0., 0., 0.};
      CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
                                 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(create_and_clean(), "create_and_clean");
 
    return tags;
  };
 
  auto get_adj_front = [&](bool subtract_crack) {
    Range adj_front;
    CHKERR mField.get_moab().get_adjacencies(*frontEdges, SPACE_DIM - 1, true,
                                             adj_front, moab::Interface::UNION);
    if (subtract_crack)
      adj_front = subtract(adj_front, *crackFaces);
    return adj_front;
  };
 
  MOFEM_LOG_CHANNEL("SELF");
 
  auto th_front_position = get_tags_vec("FrontPosition", 3);
  auto th_max_face_energy = get_tags_vec("MaxFaceEnergy", 1);
 
  if (mField.get_comm_rank() == 0) {
 
    auto get_layers_for_sides = [&](auto &side) {
      std::vector<Range> layers;
      auto get = [&]() {
        MoFEMFunctionBegin;
 
        auto get_adj = [&](auto &r, int dim) {
          Range adj;
          CHKERR mField.get_moab().get_adjacencies(r, dim, true, adj,
                                                   moab::Interface::UNION);
          return adj;
        };
 
        auto get_tets = [&](auto r) { return get_adj(r, SPACE_DIM); };
 
        Range front_nodes;
        CHKERR mField.get_moab().get_connectivity(*frontEdges, front_nodes,
                                                  true);
        Range front_faces = get_adj(front_nodes, 2);
        front_faces = subtract(front_faces, *crackFaces);
        auto front_tets = get_tets(front_nodes);
        auto front_side = intersect(side, front_tets);
        layers.push_back(front_side);
        for (;;) {
          auto adj_faces = get_skin(mField, layers.back());
          adj_faces = intersect(adj_faces, front_faces);
          auto adj_faces_tets = get_tets(adj_faces);
          adj_faces_tets = intersect(adj_faces_tets, front_tets);
          layers.push_back(unite(layers.back(), adj_faces_tets));
          if (layers.back().size() == layers[layers.size() - 2].size()) {
            break;
          }
        }
        MoFEMFunctionReturn(0);
      };
      CHK_THROW_MESSAGE(get(), "get_layers_for_sides");
      return layers;
    };
 
    auto sides_pair = get_two_sides_of_crack_surface(mField, *crackFaces);
    auto layers_top = get_layers_for_sides(sides_pair.first);
    auto layers_bottom = get_layers_for_sides(sides_pair.second);
 
#ifndef NDEBUG
    if (debug) {
      auto get_crack_adj_tets = [&](auto r) {
        Range crack_faces_conn;
        CHKERR mField.get_moab().get_connectivity(r, crack_faces_conn);
        Range crack_faces_conn_tets;
        CHKERR mField.get_moab().get_adjacencies(
            crack_faces_conn, SPACE_DIM, true, crack_faces_conn_tets,
            moab::Interface::UNION);
        return crack_faces_conn_tets;
      };
      CHKERR save_range(
          mField.get_moab(),
          "crack_tets_" +
              boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
          get_crack_adj_tets(*crackFaces));
      CHKERR save_range(mField.get_moab(), "sides_first.vtk", sides_pair.first);
      CHKERR save_range(mField.get_moab(), "sides_second.vtk",
                        sides_pair.second);
      MOFEM_LOG("EP", sev) << "Nb. layers " << layers_top.size();
      int l = 0;
      for (auto &r : layers_top) {
        MOFEM_LOG("EP", sev) << "Layer " << l << " size " << r.size();
        CHKERR save_range(
            mField.get_moab(),
            "layers_top_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
        ++l;
      }
 
      l = 0;
      for (auto &r : layers_bottom) {
        MOFEM_LOG("EP", sev) << "Layer " << l << " size " << r.size();
        CHKERR save_range(
            mField.get_moab(),
            "layers_bottom_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
        ++l;
      }
    }
#endif
 
    auto get_cross = [&](auto t_dir, auto f) {
      FTensor::Tensor1<double, SPACE_DIM> t_normal;
      CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
      t_normal.normalize();
      FTensor::Tensor1<double, SPACE_DIM> t_cross;
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
      FTENSOR_INDEX(SPACE_DIM, k);
      t_cross(i) = FTensor::levi_civita(i, j, k) * t_normal(j) * t_dir(k);
      return t_cross;
    };
 
    auto get_sense = [&](auto f, auto e) {
      int side, sense, offset;
      CHK_MOAB_THROW(mField.get_moab().side_number(f, e, side, sense, offset),
                     "get sense");
      return std::make_tuple(side, sense, offset);
    };
 
    auto calculate_edge_direction = [&](auto e, auto normalize = true) {
      const EntityHandle *conn;
      int num_nodes;
      CHKERR mField.get_moab().get_connectivity(e, conn, num_nodes, true);
      std::array<double, 6> coords;
      CHKERR mField.get_moab().get_coords(conn, num_nodes, coords.data());
      FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p0{
          &coords[0], &coords[1], &coords[2]};
      FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p1{
          &coords[3], &coords[4], &coords[5]};
      FTensor::Tensor1<double, SPACE_DIM> t_dir;
      FTENSOR_INDEX(SPACE_DIM, i);
      t_dir(i) = t_p1(i) - t_p0(i);
      if (normalize)
        t_dir.normalize();
      return t_dir;
    };
 
    auto evaluate_face_energy_and_set_orientation = [&](auto front_edges,
                                                        auto front_faces,
                                                        auto &sides_pair,
                                                        auto th_position) {
      MoFEMFunctionBegin;
 
      Tag th_face_energy;
      Tag th_material_force;
      switch (energyReleaseSelector) {
      case GRIFFITH_FORCE:
      case GRIFFITH_SKELETON:
        CHKERR mField.get_moab().tag_get_handle("GriffithForce",
                                                th_face_energy);
        // CHKERR mField.get_moab().tag_get_handle("MaterialForce",
        //                                         th_material_force);
        CHKERR mField.get_moab().tag_get_handle("MaterialForce",
                                                th_material_force);
 
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
                "Unknown energy release selector");
      };
 
      /**
       * Iterate over front edges, get adjacent faces, find maximal face energy.
       * Maximal face energy is stored in the edge. Maximal face energy is
       * magnitude of edge Griffith force.
       */
      auto find_maximal_face_energy = [&](auto front_edges, auto front_faces,
                                          auto &edge_face_max_energy_map) {
        MoFEMFunctionBegin;
 
        Range body_ents;
        CHKERR mField.get_moab().get_entities_by_dimension(0, 3, body_ents);
        auto body_skin = get_skin(mField, body_ents);
 
        Range max_faces;
 
        for (auto e : front_edges) {
 
          double griffith_force;
          CHKERR mField.get_moab().tag_get_data(th_face_energy, &e, 1,
                                                &griffith_force);
 
          Range faces;
          CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
          faces = subtract(intersect(faces, front_faces), body_skin);
          std::vector<double> face_energy(faces.size());
          CHKERR mField.get_moab().tag_get_data(th_face_energy, faces,
                                                face_energy.data());
          auto max_energy_it =
              std::max_element(face_energy.begin(), face_energy.end());
          double max_energy =
              max_energy_it != face_energy.end() ? *max_energy_it : 0;
 
          edge_face_max_energy_map[e] =
              std::make_tuple(faces[max_energy_it - face_energy.begin()],
                              griffith_force, static_cast<double>(0));
          MOFEM_LOG("EP", Sev::inform)
              << "Edge " << e << " griffith force " << griffith_force
              << " max face energy " << max_energy << " factor "
              << max_energy / griffith_force;
 
          max_faces.insert(faces[max_energy_it - face_energy.begin()]);
        }
 
#ifndef NDEBUG
        if (debug) {
          CHKERR save_range(
              mField.get_moab(),
              "max_faces_" +
                  boost::lexical_cast<std::string>(mField.get_comm_rank()) +
                  ".vtk",
              max_faces);
        }
#endif
 
        MoFEMFunctionReturn(0);
      };
 
      /**
       * For each front edge, find maximal face energy and orientation. This is
       * by finding angle between edge material force and maximal face normal
       *
       */
      auto calculate_face_orientation = [&](auto &edge_face_max_energy_map) {
        MoFEMFunctionBegin;
 
        auto up_down_face = [&](
 
                                auto &face_angle_map_up,
                                auto &face_angle_map_down
 
                            ) {
          MoFEMFunctionBegin;
 
          for (auto &m : edge_face_max_energy_map) {
            auto e = m.first;
            auto [max_face, energy, opt_angle] = m.second;
 
            Range faces;
            CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
            faces = intersect(faces, front_faces);
            Range adj_tets; // tetrahedrons adjacent to the face
            CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
                                                     false, adj_tets,
                                                     moab::Interface::UNION);
            if (adj_tets.size()) {
 
              Range adj_tets; // tetrahedrons adjacent to the face
              CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
                                                       false, adj_tets,
                                                       moab::Interface::UNION);
              if (adj_tets.size()) {
 
                Range adj_tets_faces;
                // get faces
                CHKERR mField.get_moab().get_adjacencies(
                    adj_tets, SPACE_DIM - 1, false, adj_tets_faces,
                    moab::Interface::UNION);
                adj_tets_faces = intersect(adj_tets_faces, faces);
                FTENSOR_INDEX(SPACE_DIM, i);
 
                // cross product of face normal and edge direction
                auto t_cross_max =
                    get_cross(calculate_edge_direction(e, true), max_face);
                auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
                t_cross_max(i) *= sense_max;
 
                for (auto t : adj_tets) {
                  Range adj_tets_faces;
                  CHKERR mField.get_moab().get_adjacencies(
                      &t, 1, SPACE_DIM - 1, false, adj_tets_faces);
                  adj_tets_faces = intersect(adj_tets_faces, faces);
                  adj_tets_faces =
                      subtract(adj_tets_faces, Range(max_face, max_face));
 
                  if (adj_tets_faces.size() == 1) {
 
                    // cross product of adjacent face normal and edge
                    // direction
                    auto t_cross = get_cross(calculate_edge_direction(e, true),
                                             adj_tets_faces[0]);
                    auto [side, sense, offset] =
                        get_sense(adj_tets_faces[0], e);
                    t_cross(i) *= sense;
                    double dot = t_cross(i) * t_cross_max(i);
                    auto angle = std::acos(dot);
 
                    double face_energy;
                    CHKERR mField.get_moab().tag_get_data(
                        th_face_energy, adj_tets_faces, &face_energy);
 
                    auto [side_face, sense_face, offset_face] =
                        get_sense(t, max_face);
 
                    if (sense_face > 0) {
                      face_angle_map_up[e] = std::make_tuple(face_energy, angle,
                                                             adj_tets_faces[0]);
 
                    } else {
                      face_angle_map_down[e] = std::make_tuple(
                          face_energy, -angle, adj_tets_faces[0]);
                    }
                  }
                }
              }
            }
          }
 
          MoFEMFunctionReturn(0);
        };
 
        auto calc_optimal_angle = [&](
 
                                      auto &face_angle_map_up,
                                      auto &face_angle_map_down
 
                                  ) {
          MoFEMFunctionBegin;
 
          for (auto &m : edge_face_max_energy_map) {
            auto e = m.first;
            auto &[max_face, e0, a0] = m.second;
 
            if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
 
              if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
                  face_angle_map_down.find(e) == face_angle_map_down.end()) {
                // Do nothing
              } else {
 
                switch (energyReleaseSelector) {
                case GRIFFITH_FORCE:
                case GRIFFITH_SKELETON: {
 
                  Tag th_material_force;
                  CHKERR mField.get_moab().tag_get_handle("MaterialForce",
                                                          th_material_force);
                  FTensor::Tensor1<double, SPACE_DIM> t_material_force;
                  CHKERR mField.get_moab().tag_get_data(
                      th_material_force, &e, 1, &t_material_force(0));
                  auto material_force_magnitude = t_material_force.l2();
                  if (material_force_magnitude <
                      std::numeric_limits<double>::epsilon()) {
                    a0 = 0;
 
                  } else {
 
                    auto t_edge_dir = calculate_edge_direction(e, true);
                    auto t_cross_max = get_cross(t_edge_dir, max_face);
                    auto [side, sense, offset] = get_sense(max_face, e);
                    t_cross_max(sense) *= sense;
 
                    FTENSOR_INDEX(SPACE_DIM, I);
                    FTENSOR_INDEX(SPACE_DIM, J);
                    FTENSOR_INDEX(SPACE_DIM, K);
 
                    t_material_force.normalize();
                    t_cross_max.normalize();
                    FTensor::Tensor1<double, SPACE_DIM> t_cross;
                    t_cross(I) = FTensor::levi_civita(I, J, K) *
                                 t_material_force(J) * t_cross_max(K);
                    a0 = -std::asin(t_cross(I) * t_edge_dir(I));
 
                    MOFEM_LOG("EP", sev)
                        << "Optimal angle " << a0 << " energy " << e0;
                  }
                  break;
                }
                default: {
 
                  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
                          "Unknown energy release selector");
                }
                }
              }
            }
          }
 
          MoFEMFunctionReturn(0);
        };
 
        std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
            face_angle_map_up;
        std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
            face_angle_map_down;
        CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
        CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
 
#ifndef NDEBUG
        if (debug) {
          auto th_angle = get_tags_vec("Angle", 1);
          Range up;
          for (auto &m : face_angle_map_up) {
            auto [e, a, face] = m.second;
            up.insert(face);
            CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
          }
          Range down;
          for (auto &m : face_angle_map_down) {
            auto [e, a, face] = m.second;
            down.insert(face);
            CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
          }
 
          Range max_energy_faces;
          for (auto &m : edge_face_max_energy_map) {
            auto [face, e, angle] = m.second;
            max_energy_faces.insert(face);
            CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1,
                                                  &angle);
          }
          if (mField.get_comm_rank() == 0) {
            CHKERR save_range(mField.get_moab(), "up_faces.vtk", up);
            CHKERR save_range(mField.get_moab(), "down_faces.vtk", down);
            CHKERR save_range(mField.get_moab(), "max_energy_faces.vtk",
                              max_energy_faces);
          }
        }
#endif // NDEBUG
 
        MoFEMFunctionReturn(0);
      };
 
      auto get_conn = [&](auto e) {
        Range conn;
        CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
                       "get conn");
        return conn;
      };
 
      auto get_adj = [&](auto e, auto dim) {
        Range adj;
        CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
                           e, dim, false, adj, moab::Interface::UNION),
                       "get adj");
        return adj;
      };
 
      auto get_coords = [&](auto v) {
        FTensor::Tensor1<double, 3> t_coords;
        CHK_MOAB_THROW(mField.get_moab().get_coords(v, &t_coords(0)),
                       "get coords");
        return t_coords;
      };
 
      // calculate normal of the max energy face
      auto get_rotated_normal = [&](auto e, auto f, auto angle) {
        FTENSOR_INDEX(SPACE_DIM, i);
        FTENSOR_INDEX(SPACE_DIM, j);
        auto t_edge_dir = calculate_edge_direction(e, true);
        auto [side, sense, offset] = get_sense(f, e);
        t_edge_dir(i) *= sense;
        t_edge_dir.normalize();
        t_edge_dir(i) *= angle;
        auto t_R = LieGroups::SO3::exp(t_edge_dir, angle);
        FTensor::Tensor1<double, SPACE_DIM> t_normal;
        mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
        FTensor::Tensor1<double, SPACE_DIM> t_rotated_normal;
        t_rotated_normal(i) = t_R(i, j) * t_normal(j);
        return std::make_tuple(t_normal, t_rotated_normal);
      };
 
      auto set_coord = [&](auto v, auto &adj_vertex_tets_verts, auto &coords,
                           auto &t_move, auto gamma) {
        auto index = adj_vertex_tets_verts.index(v);
        if (index >= 0) {
          for (auto ii : {0, 1, 2}) {
            coords[3 * index + ii] += gamma * t_move(ii);
          }
          return true;
        }
        return false;
      };
 
      auto tets_quality = [&](auto quality, auto &adj_vertex_tets_verts,
                              auto &adj_vertex_tets, auto &coords) {
        for (auto t : adj_vertex_tets) {
          const EntityHandle *conn;
          int num_nodes;
          CHKERR mField.get_moab().get_connectivity(t, conn, num_nodes, true);
          std::array<double, 12> tet_coords;
          for (auto n = 0; n != 4; ++n) {
            auto index = adj_vertex_tets_verts.index(conn[n]);
            if (index < 0) {
              CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong index");
            }
            for (auto ii = 0; ii != 3; ++ii) {
              tet_coords[3 * n + ii] = coords[3 * index + ii];
            }
          }
          double q = Tools::volumeLengthQuality(tet_coords.data());
          if (!std::isnormal(q))
            q = -2;
          quality = std::min(quality, q);
        };
 
        return quality;
      };
 
      auto calculate_free_face_node_displacement =
          [&](auto &edge_face_max_energy_map) {
            // get edges adjacent to vertex along which nodes are moving
            auto get_vertex_edges = [&](auto vertex) {
              Range vertex_edges; // edges adjacent to vertex
 
              auto impl = [&]() {
                MoFEMFunctionBegin;
                CHKERR mField.get_moab().get_adjacencies(vertex, 1, false,
                                                         vertex_edges);
                vertex_edges = subtract(vertex_edges, front_verts_edges);
 
                if (boundary_skin_verts.size() &&
                    boundary_skin_verts.find(vertex[0]) !=
                        boundary_skin_verts.end()) {
                  MOFEM_LOG("EP", sev) << "Boundary vertex";
                  vertex_edges = intersect(vertex_edges, body_skin_edges);
                }
                if (geometry_edges_verts.size() &&
                    geometry_edges_verts.find(vertex[0]) !=
                        geometry_edges_verts.end()) {
                  MOFEM_LOG("EP", sev) << "Geometry edge vertex";
                  vertex_edges = intersect(vertex_edges, geometry_edges);
                }
                if (crack_faces_verts.size() &&
                    crack_faces_verts.find(vertex[0]) !=
                        crack_faces_verts.end()) {
                  MOFEM_LOG("EP", sev) << "Crack face vertex";
                  vertex_edges = intersect(vertex_edges, crack_faces_edges);
                }
                MoFEMFunctionReturn(0);
              };
 
              CHK_THROW_MESSAGE(impl(), "get_vertex_edges");
 
              return vertex_edges;
            };
 
            // vector of rotated faces, edge along node is moved, moved edge,
            // moved displacement, quality, cardinality, gamma
            using Bundle = std::vector<
 
                std::tuple<EntityHandle, EntityHandle, EntityHandle,
                           FTensor::Tensor1<double, 3>, double, double, double>
 
                >;
            std::map<EntityHandle, Bundle> edge_bundle_map;
 
            for (auto &m : edge_face_max_energy_map) {
 
              auto edge = m.first;
              auto &[max_face, energy, opt_angle] = m.second;
 
              // calculate rotation of max energy face
              auto [t_normal, t_rotated_normal] =
                  get_rotated_normal(edge, max_face, opt_angle);
 
              auto front_vertex = get_conn(Range(m.first, m.first));
              auto adj_tets = get_adj(Range(max_face, max_face), 3);
              auto adj_tets_faces = get_adj(adj_tets, 2);
              auto adj_front_faces = subtract(
                  intersect(get_adj(Range(edge, edge), 2), adj_tets_faces),
                  *crackFaces);
              if (adj_front_faces.size() > 3)
                CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                                  "adj_front_faces.size()>3");
 
              FTensor::Tensor1<double, SPACE_DIM> t_material_force;
              CHKERR mField.get_moab().tag_get_data(th_material_force, &edge, 1,
                                                    &t_material_force(0));
              std::vector<double> griffith_energy(adj_front_faces.size());
              CHKERR mField.get_moab().tag_get_data(
                  th_face_energy, adj_front_faces, griffith_energy.data());
 
              auto set_edge_bundle = [&](auto min_gamma) {
                for (auto rotated_f : adj_front_faces) {
 
                  double rotated_face_energy =
                      griffith_energy[adj_front_faces.index(rotated_f)];
 
                  auto vertex = subtract(get_conn(Range(rotated_f, rotated_f)),
                                         front_vertex);
                  if (vertex.size() != 1) {
                    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                                      "Wrong number of vertex to move");
                  }
                  auto front_vertex_edges_vertex = get_conn(
                      intersect(get_adj(front_vertex, 1), crack_faces_edges));
                  vertex = subtract(
                      vertex, front_vertex_edges_vertex); // vertex free to move
                  if (vertex.empty()) {
                    continue;
                  }
 
                  auto face_cardinality = [&](auto f, auto &seen_front_edges) {
                    auto whole_front =
                        unite(*frontEdges,
                              subtract(body_skin_edges, crack_faces_edges));
                    auto faces = Range(f, f);
                    int c = 0;
                    for (; c < 10; ++c) {
                      auto front_edges =
                          subtract(get_adj(faces, 1), seen_front_edges);
                      if (front_edges.size() == 0) {
                        return 0;
                      }
                      auto front_connected_edges =
                          intersect(front_edges, whole_front);
                      if (front_connected_edges.size()) {
                        seen_front_edges.merge(front_connected_edges);
                        return c;
                      }
                      faces.merge(get_adj(front_edges, 2));
                      ++c;
                    }
                    return c;
                  };
 
                  Range seen_edges = Range(edge, edge);
                  double rotated_face_cardinality = face_cardinality(
                      rotated_f,
                      seen_edges); // add cardinality of max energy
                                   // face to rotated face cardinality
                  // rotated_face_cardinality +=
                  //     face_cardinality(max_face, seen_edges);
                  rotated_face_cardinality = std::max(rotated_face_cardinality,
                                                      1.); // at least one edge
 
                  auto t_vertex_coords = get_coords(vertex);
                  auto vertex_edges = get_vertex_edges(vertex);
 
                  EntityHandle f0 = front_vertex[0];
                  EntityHandle f1 = front_vertex[1];
                  FTensor::Tensor1<double, 3> t_v_e0, t_v_e1;
                  CHKERR mField.get_moab().get_coords(&f0, 1, &t_v_e0(0));
                  CHKERR mField.get_moab().get_coords(&f1, 1, &t_v_e1(0));
 
                  FTENSOR_INDEX(SPACE_DIM, i);
                  for (auto e_used_to_move_detection : vertex_edges) {
                    auto edge_conn = get_conn(Range(e_used_to_move_detection,
                                                    e_used_to_move_detection));
                    edge_conn = subtract(edge_conn, vertex);
                    // Find displacement of the edge such that dot porduct with
                    // normal is zero.
                    //
                    // { (t_v0 - t_vertex_coords) + gamma * (t_v3 -
                    // t_vertex_coords) } * n = 0
                    // where t_v0 is the edge vertex, t_v3 is the edge end
                    // point, n is the rotated normal of the face gamma is the
                    // factor by which the edge is moved
                    FTensor::Tensor1<double, 3> t_v0;
                    t_v0(i) = (t_v_e0(i) + t_v_e1(i)) / 2;
                    FTensor::Tensor1<double, 3> t_v3;
                    CHKERR mField.get_moab().get_coords(edge_conn, &t_v3(0));
                    auto a =
                        (t_v0(i) - t_vertex_coords(i)) * t_rotated_normal(i);
                    auto b =
                        (t_v3(i) - t_vertex_coords(i)) * t_rotated_normal(i);
                    auto gamma = a / b;
 
                    constexpr double eps =
                        std::numeric_limits<double>::epsilon();
                    if (std::isnormal(gamma) && gamma < 1.0 - eps &&
                        gamma > -0.1) {
                      FTensor::Tensor1<double, 3> t_move;
                      t_move(i) = gamma * (t_v3(i) - t_vertex_coords(i));
 
                      auto check_rotated_face_directoon = [&]() {
                        FTensor::Tensor1<double, SPACE_DIM> t_delta;
                        t_delta(i) = t_vertex_coords(i) + t_move(i) - t_v0(i);
                        t_delta.normalize();
                        auto dot =
                            (t_material_force(i) / t_material_force.l2()) *
                            t_delta(i);
                        return -dot > 0 ? true : false;
                      };
 
                      if (check_rotated_face_directoon()) {
 
                        MOFEM_LOG("EP", Sev::verbose)
                            << "Crack edge " << edge << " moved face "
                            << rotated_f
                            << " edge: " << e_used_to_move_detection
                            << " face direction/energy " << rotated_face_energy
                            << " face cardinality " << rotated_face_cardinality
                            << " gamma: " << gamma;
 
                        auto &bundle = edge_bundle_map[edge];
                        bundle.emplace_back(rotated_f, e_used_to_move_detection,
                                            vertex[0], t_move, 1,
                                            rotated_face_cardinality, gamma);
                      }
                    }
                  }
                }
              };
 
              set_edge_bundle(std::numeric_limits<double>::epsilon());
              if (edge_bundle_map[edge].empty()) {
                set_edge_bundle(-1.);
              }
            }
 
            return edge_bundle_map;
          };
 
      auto get_sort_by_energy = [&](auto &edge_face_max_energy_map) {
        std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
            sort_by_energy;
 
        for (auto &m : edge_face_max_energy_map) {
          auto e = m.first;
          auto &[max_face, energy, opt_angle] = m.second;
          auto abs_energy = std::abs(energy);
          sort_by_energy[abs_energy] = std::make_tuple(e, max_face, opt_angle);
        }
 
        return sort_by_energy;
      };
 
      auto set_tag = [&](auto &&adj_edges_map, auto &&sort_by_energy) {
        MoFEMFunctionBegin;
 
        Tag th_face_pressure;
        CHK_MOAB_THROW(
            mField.get_moab().tag_get_handle("FacePressure", th_face_pressure),
            "get tag");
        auto get_face_pressure = [&](auto face) {
          double pressure;
          CHK_MOAB_THROW(mField.get_moab().tag_get_data(th_face_pressure, &face,
                                                        1, &pressure),
                         "get rag data");
          return pressure;
        };
 
        MOFEM_LOG("EPSELF", Sev::inform)
            << "Number of edges to check " << sort_by_energy.size();
 
        enum face_energy { POSITIVE, NEGATIVE };
        constexpr bool skip_negative = true;
 
        for (auto fe : {face_energy::POSITIVE, face_energy::NEGATIVE}) {
 
          std::vector<double> energies;
          double max_pressure = -1;
          // check max energies and average all energies along the crack front
          // extract max pressure along the crack front
          for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
               ++it) {
            auto energy = it->first;
            auto [max_edge, max_face, opt_angle] = it->second;
 
            auto face_pressure = get_face_pressure(max_face);
            MOFEM_LOG("EPSELF", Sev::inform)
                << "Faces to check: " << max_face << " energy " << energy
                << " face pressure " << face_pressure;
 
            const bool pressure_check =
                propagateUnderCompression || face_pressure > crackingAtol;
            if (energy > 0 && pressure_check) {
              energies.push_back(energy);
            }
            max_pressure = std::max(max_pressure, face_pressure);
          }
 
          double average_energy = 0;
          if (!energies.empty()) {
            average_energy =
                std::accumulate(energies.begin(), energies.end(), 0.) /
                energies.size();
          }
 
          MOFEM_LOG("EPSELF", Sev::inform)
              << "Average energy Griffiths energy of crack front "
              << average_energy;
 
          bool positive_pressure_face_found = false;
 
          // iterate edges wih maximal energy, and make them seed. Such edges,
          // will most likely will have also smallest node displacement
          for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
               ++it) {
 
            auto energy = it->first;
            auto [max_edge, max_face, opt_angle] = it->second;
 
            auto face_pressure = get_face_pressure(max_face);
            if (skip_negative) {
              if (fe == face_energy::POSITIVE) {
                if (face_pressure <
                    -(crackingAtol + crackingRtol * std::abs(max_pressure))) {
                  MOFEM_LOG("EPSELF", Sev::inform)
                      << "Skip negative face " << max_face << " with energy "
                      << energy << " and pressure " << face_pressure;
                  continue;
                }
              }
            }
 
            if (fe == face_energy::POSITIVE)
              positive_pressure_face_found = true;
 
            MOFEM_LOG("EPSELF", Sev::inform)
                << "Check face " << max_face << " edge " << max_edge
                << " energy " << energy << " optimal angle " << opt_angle
                << " face pressure " << face_pressure;
 
            // store energy of max face
            if (!average_energy) {
              MOFEM_LOG("EPSELF", Sev::warning)
                  << "Average energy is zero, setting max Griffiths energy to "
                     "current energy "
                  << energy;
              average_energy = energy;
            }
            avgGriffithsEnergy = average_energy;
            auto jt = adj_edges_map.find(max_edge);
            if (jt == adj_edges_map.end()) {
              MOFEM_LOG("EPSELF", Sev::warning)
                  << "Edge " << max_edge << " not found in adj_edges_map";
              continue;
            }
            auto &bundle = jt->second;
 
            auto find_max_in_bundle_impl = [&](auto edge, auto &bundle,
                                               auto gamma) {
              FTENSOR_INDEX(SPACE_DIM, i);
 
              EntityHandle vertex_max = 0;
              EntityHandle face_max = 0;
              EntityHandle move_edge_max = 0;
              double max_quality = -2;
              double max_quality_evaluated = -2;
              double min_cardinality = std::numeric_limits<double>::max();
 
              FTensor::Tensor1<double, SPACE_DIM> t_move_last{0., 0., 0.};
 
              for (auto &b : bundle) {
                auto &[face, move_edge, vertex, t_move, quality, cardinality,
                       edge_gamma] = b;
 
                auto adj_vertex_tets = get_adj(Range(vertex, vertex), 3);
                auto adj_vertex_tets_verts = get_conn(adj_vertex_tets);
                std::vector<double> coords(3 * adj_vertex_tets_verts.size());
                CHK_MOAB_THROW(mField.get_moab().get_coords(
                                   adj_vertex_tets_verts, coords.data()),
                               "get coords");
 
                set_coord(vertex, adj_vertex_tets_verts, coords, t_move, gamma);
                quality = tets_quality(quality, adj_vertex_tets_verts,
                                       adj_vertex_tets, coords);
 
                auto eval_quality = [](auto q, auto c, auto edge_gamma) {
                  if (q < 0) {
                    return q;
                  } else {
                    return ((edge_gamma < 0) ? (q / 2) : q) / pow(c, 2);
                  }
                };
 
                if (eval_quality(quality, cardinality, edge_gamma) >=
                    max_quality_evaluated) {
                  max_quality = quality;
                  min_cardinality = cardinality;
                  vertex_max = vertex;
                  face_max = face;
                  move_edge_max = move_edge;
                  t_move_last(i) = t_move(i);
                  max_quality_evaluated =
                      eval_quality(max_quality, min_cardinality, edge_gamma);
                }
              }
 
              return std::make_tuple(vertex_max, face_max, t_move_last,
                                     max_quality, min_cardinality);
            };
 
            auto find_max_in_bundle = [&](auto edge, auto &bundle) {
              auto b_org_bundle = bundle;
              auto r = find_max_in_bundle_impl(edge, bundle, 1.);
              auto &[vertex_max, face_max, t_move_last, max_quality,
                     cardinality] = r;
              if (max_quality < 0) {
                for (double gamma = 0.95; gamma >= 0.45; gamma -= 0.05) {
                  bundle = b_org_bundle;
                  r = find_max_in_bundle_impl(edge, bundle, gamma);
                  auto &[vertex_max, face_max, t_move_last, max_quality,
                         cardinality] = r;
                  MOFEM_LOG("EPSELF", Sev::warning)
                      << "Back tracking: gamma " << gamma << " edge " << edge
                      << " quality " << max_quality << " cardinality "
                      << cardinality;
                  if (max_quality > 0.01) {
                    FTENSOR_INDEX(SPACE_DIM, I);
                    t_move_last(I) *= gamma;
                    return r;
                  }
                }
                FTENSOR_INDEX(SPACE_DIM, I);
                t_move_last(I) = 0;
              }
              return r;
            };
 
            // set tags with displacement of node and face energy
            auto set_tag_to_vertex_and_face = [&](auto &&r, auto &quality) {
              MoFEMFunctionBegin;
              auto &[v, f, t_move, q, cardinality] = r;
 
              if ((q > 0 && std::isnormal(q)) && energy > 0) {
 
                MOFEM_LOG("EPSELF", Sev::inform)
                    << "Set tag: vertex " << v << " face " << f << " "
                    << max_edge << " move " << t_move << " energy " << energy
                    << " quality " << q << " cardinality " << cardinality;
                CHKERR mField.get_moab().tag_set_data(th_position[0], &v, 1,
                                                      &t_move(0));
                CHKERR mField.get_moab().tag_set_data(th_max_face_energy[0], &f,
                                                      1, &energy);
              }
 
              quality = q;
              MoFEMFunctionReturn(0);
            };
 
            double quality = -2;
            CHKERR set_tag_to_vertex_and_face(
 
                find_max_in_bundle(max_edge, bundle),
 
                quality
 
            );
 
            if (quality > 0 && std::isnormal(quality) && energy > 0) {
              MOFEM_LOG("EPSELF", Sev::inform)
                  << "Crack face set with quality: " << quality;
              MoFEMFunctionReturnHot(0);
            }
          }
 
          if (fe == face_energy::POSITIVE && !positive_pressure_face_found) {
            if (!propagateUnderCompression) {
              potentialCrackArrest = true;
              MOFEM_LOG("EPSELF", Sev::warning)
                  << "POTENTIAL ARREST: No suitable face found with positive "
                     "face pressure to propagate crack";
            } else {
              MOFEM_LOG("EPSELF", Sev::warning)
                  << "POTENTIAL ARREST: No suitable face found with positive "
                     "face pressure to propagate crack; continuing because "
                     "propagation under compression is enabled";
            }
          }
 
          if (!skip_negative)
            break;
        }
 
        MoFEMFunctionReturn(0);
      };
 
      // map: {edge, {face, energy, optimal_angle}}
      MOFEM_LOG("EP", sev) << "Calculate orientation";
      std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
          edge_face_max_energy_map;
      CHKERR find_maximal_face_energy(front_edges, front_faces,
                                      edge_face_max_energy_map);
      CHKERR calculate_face_orientation(edge_face_max_energy_map);
 
      MOFEM_LOG("EP", sev) << "Calculate node positions";
      CHKERR set_tag(
 
          calculate_free_face_node_displacement(edge_face_max_energy_map),
          get_sort_by_energy(edge_face_max_energy_map)
 
      );
 
      MoFEMFunctionReturn(0);
    };
 
    auto get_max_griffith_force = [&](auto r) {
      auto &moab = mField.get_moab();
      std::vector<double> gc(r.size());
      Tag th_gc;
      CHKERR moab.tag_get_handle("GriffithForce", th_gc);
      CHKERR moab.tag_get_data(th_gc, r, gc.data());
      double max_griffith_force = 0;
      for (size_t i = 0; i < r.size(); ++i) {
        max_griffith_force = std::max(max_griffith_force, std::abs(gc[i]));
      }
      return max_griffith_force;
    };
 
    MOFEM_LOG("EP", sev) << "Front edges " << frontEdges->size();
    if (std::abs(get_max_griffith_force(get_adj_front(true))) >
        std::numeric_limits<double>::epsilon()) {
      CHKERR evaluate_face_energy_and_set_orientation(
          *frontEdges, get_adj_front(true), sides_pair, th_front_position);
    } else {
      auto adj_front = get_adj_front(true);
      double zero[] = {0., 0., 0.};
      CHKERR mField.get_moab().tag_clear_data(th_front_position[0], adj_front,
                                              zero);
    }
  }
 
  // exchange positions and energies from processor zero to all other
  CHKERR VecZeroEntries(vertexExchange.second);
  CHKERR VecGhostUpdateBegin(vertexExchange.second, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(vertexExchange.second, INSERT_VALUES,
                           SCATTER_FORWARD);
  CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
      mField.get_moab(), vertexExchange, th_front_position[0]);
  CHKERR VecZeroEntries(faceExchange.second);
  CHKERR VecGhostUpdateBegin(faceExchange.second, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(faceExchange.second, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
      mField.get_moab(), faceExchange, th_max_face_energy[0]);
 
  auto get_max_moved_faces = [&]() {
    Range max_moved_faces;
    auto adj_front = get_adj_front(false);
    std::vector<double> face_energy(adj_front.size());
    CHKERR mField.get_moab().tag_get_data(th_max_face_energy[0], adj_front,
                                          face_energy.data());
    for (int i = 0; i != adj_front.size(); ++i) {
      if (face_energy[i] > std::numeric_limits<double>::epsilon()) {
        max_moved_faces.insert(adj_front[i]);
      }
    }
 
    return boost::make_shared<Range>(max_moved_faces);
  };
 
  // move all faces with energy larger than 0
  maxMovedFaces = get_max_moved_faces();
  MOFEM_LOG("EP", sev) << "Number of of moved faces: " << maxMovedFaces->size();
 
#ifndef NDEBUG
  if (debug) {
    CHKERR save_range(
        mField.get_moab(),
        "max_moved_faces_" +
            boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
        *maxMovedFaces);
  }
#endif
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setNewFrontCoordinates() {
  MoFEMFunctionBegin;
 
  if (!maxMovedFaces)
    MoFEMFunctionReturnHot(0);
 
  Tag th_front_position;
  auto rval =
      mField.get_moab().tag_get_handle("FrontPosition", th_front_position);
  if (rval == MB_SUCCESS && maxMovedFaces) {
    Range verts;
    CHKERR mField.get_moab().get_connectivity(*maxMovedFaces, verts, true);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(verts);
    std::vector<double> coords(3 * verts.size());
    CHKERR mField.get_moab().get_coords(verts, coords.data());
    std::vector<double> pos(3 * verts.size());
    CHKERR mField.get_moab().tag_get_data(th_front_position, verts, pos.data());
    for (int i = 0; i != 3 * verts.size(); ++i) {
      coords[i] += pos[i];
    }
    CHKERR mField.get_moab().set_coords(verts, coords.data());
    double zero[] = {0., 0., 0.};
    CHKERR mField.get_moab().tag_clear_data(th_front_position, verts, zero);
  }
 
#ifndef NDEBUG
  constexpr bool debug = false;
  if (debug) {
 
    CHKERR save_range(
        mField.get_moab(),
        "set_coords_faces_" +
            boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
        *maxMovedFaces);
  }
#endif
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::addCrackSurfaces(const bool debug) {
  MoFEMFunctionBegin;
 
  constexpr bool potential_crack_debug = false;
  if constexpr (potential_crack_debug) {
 
    auto add_ents = get_range_from_block(mField, "POTENTIAL", SPACE_DIM - 1);
    Range crack_front_verts;
    CHKERR mField.get_moab().get_connectivity(*frontEdges, crack_front_verts,
                                              true);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        crack_front_verts);
    Range crack_front_faces;
    CHKERR mField.get_moab().get_adjacencies(crack_front_verts, SPACE_DIM - 1,
                                             true, crack_front_faces,
                                             moab::Interface::UNION);
    crack_front_faces = intersect(crack_front_faces, add_ents);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        crack_front_faces);
    CHKERR mField.getInterface<MeshsetsManager>()->addEntitiesToMeshset(
        BLOCKSET, addCrackMeshsetId, crack_front_faces);
  }
 
  auto get_crack_faces = [&]() {
    if (maxMovedFaces) {
      return unite(*crackFaces, *maxMovedFaces);
    } else {
      return *crackFaces;
    }
  };
 
  auto get_extended_crack_faces = [&]() {
    auto get_faces_of_crack_front_verts = [&](auto crack_faces_org) {
      ParallelComm *pcomm =
          ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
 
      Range crack_faces;
 
      if (!pcomm->rank()) {
 
        auto get_nodes = [&](auto &&e) {
          Range nodes;
          CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, nodes, true),
                         "get connectivity");
          return nodes;
        };
 
        auto get_adj = [&](auto &&e, auto dim,
                           auto t = moab::Interface::UNION) {
          Range adj;
          CHK_MOAB_THROW(
              mField.get_moab().get_adjacencies(e, dim, true, adj, t),
              "get adj");
          return adj;
        };
 
        Range body_ents;
        CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                           body_ents);
        auto body_skin = get_skin(mField, body_ents);
        auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
        auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
        auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
        auto front_block_nodes = get_nodes(front_block_edges);
 
        size_t s;
        do {
          s = crack_faces.size();
 
          auto crack_face_nodes = get_nodes(crack_faces_org);
          auto crack_faces_edges =
              get_adj(crack_faces_org, 1, moab::Interface::UNION);
 
          auto crack_skin = get_skin(mField, crack_faces_org);
          front_block_edges = subtract(front_block_edges, crack_skin);
          auto crack_skin_nodes = get_nodes(crack_skin);
          crack_skin_nodes.merge(front_block_nodes);
 
          auto crack_skin_faces =
              get_adj(crack_skin, 2, moab::Interface::UNION);
          crack_skin_faces =
              subtract(subtract(crack_skin_faces, crack_faces_org), body_skin);
 
          crack_faces = crack_faces_org;
          for (auto f : crack_skin_faces) {
            auto edges = intersect(
                get_adj(Range(f, f), 1, moab::Interface::UNION), crack_skin);
 
            // if other edge is part of body skin, e.g. crack punching through
            // body surface
            if (edges.size() == 2) {
              edges.merge(
                  intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
                            body_skin_edges));
            }
 
            if (edges.size() == 2) {
              auto edge_conn = get_nodes(Range(edges));
              auto faces = intersect(get_adj(edges, 2, moab::Interface::UNION),
                                     crack_faces_org);
              if (faces.size() == 2) {
                auto edge0_conn = get_nodes(Range(edges[0], edges[0]));
                auto edge1_conn = get_nodes(Range(edges[1], edges[1]));
                auto edges_conn = intersect(intersect(edge0_conn, edge1_conn),
                                            crack_skin_nodes); // node at apex
                if (edges_conn.size() == 1) {
 
                  auto node_edges =
                      subtract(intersect(get_adj(edges_conn, 1,
                                                 moab::Interface::INTERSECT),
                                         crack_faces_edges),
                               crack_skin); // nodes on crack surface, but not
                                            // at the skin
 
                  if (node_edges.size()) {
                    FTENSOR_INDEX(SPACE_DIM, i);
                    FTensor::Tensor1<double, SPACE_DIM> t_v0;
                    CHKERR mField.get_moab().get_coords(edges_conn, &t_v0(0));
 
                    auto get_t_dir = [&](auto e_conn) {
                      auto other_node = subtract(e_conn, edges_conn);
                      FTensor::Tensor1<double, SPACE_DIM> t_dir;
                      CHKERR mField.get_moab().get_coords(other_node,
                                                          &t_dir(0));
                      t_dir(i) -= t_v0(i);
                      return t_dir;
                    };
 
                    FTensor::Tensor1<double, SPACE_DIM> t_ave_dir;
                    t_ave_dir(i) =
                        get_t_dir(edge0_conn)(i) + get_t_dir(edge1_conn)(i);
 
                    FTensor::Tensor1<double, SPACE_DIM> t_crack_surface_ave_dir;
                    t_crack_surface_ave_dir(i) = 0;
                    for (auto e : node_edges) {
                      auto e_conn = get_nodes(Range(e, e));
                      auto t_dir = get_t_dir(e_conn);
                      t_crack_surface_ave_dir(i) += t_dir(i);
                    }
 
                    auto dot = t_ave_dir(i) * t_crack_surface_ave_dir(i);
                    // ave edges is in opposite direction to crack surface, so
                    // thus crack is not turning back
                    if (dot < -std::numeric_limits<double>::epsilon()) {
                      crack_faces.insert(f);
                    }
                  } else {
                    crack_faces.insert(f);
                  }
                }
              }
            } else if (edges.size() == 3) {
              crack_faces.insert(f);
            }
 
            // if other edge is part of geometry edge, e.g. keyway
            if (edges.size() == 1) {
              edges.merge(
                  intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
                            geometry_edges));
              edges.merge(
                  intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
                            front_block_edges));
              if (edges.size() == 2) {
                crack_faces.insert(f);
                continue;
              }
            }
          }
 
          crack_faces_org = crack_faces;
 
        } while (s != crack_faces.size());
      };
 
      return crack_faces; // send_type(mField, crack_faces, MBTRI);
    };
 
    return get_faces_of_crack_front_verts(get_crack_faces());
  };
 
  if (debug) {
    CHKERR save_range(mField.get_moab(), "new_crack_surface_debug.vtk",
                      get_extended_crack_faces());
  }
 
  auto reconstruct_crack_faces = [&](auto crack_faces) {
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
 
    auto impl = [&]() {
      MoFEMFunctionBegin;
 
      Range new_crack_faces;
      if (!pcomm->rank()) {
 
        auto get_nodes = [&](auto &&e) {
          Range nodes;
          CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, nodes, true),
                         "get connectivity");
          return nodes;
        };
 
        auto get_adj = [&](auto &&e, auto dim,
                           auto t = moab::Interface::UNION) {
          Range adj;
          CHK_MOAB_THROW(
              mField.get_moab().get_adjacencies(e, dim, true, adj, t),
              "get adj");
          return adj;
        };
 
        auto get_test_on_crack_surface = [&]() {
          auto crack_faces_nodes =
              get_nodes(crack_faces); // nodes on crac faces
          auto crack_faces_tets =
              get_adj(crack_faces_nodes, 3,
                      moab::Interface::UNION); // adjacent
                                               // tets to
                                               // crack
                                               // faces throug nodes
          auto crack_faces_tets_nodes =
              get_nodes(crack_faces_tets); // nodes on crack faces tets
          crack_faces_tets_nodes =
              subtract(crack_faces_tets_nodes, crack_faces_nodes);
          crack_faces_tets =
              subtract(crack_faces_tets, get_adj(crack_faces_tets_nodes, 3,
                                                 moab::Interface::UNION));
          new_crack_faces =
              get_adj(crack_faces_tets, 2,
                      moab::Interface::UNION); // adjacency faces to crack
                                               // faces through tets
          new_crack_faces.merge(crack_faces);  // add original crack faces
 
          return std::make_tuple(new_crack_faces, crack_faces_tets);
        };
 
        auto carck_faces_test_edges = [&](auto faces, auto tets) {
          auto adj_tets_faces = get_adj(tets, 2, moab::Interface::UNION);
          auto adj_faces_edges = get_adj(subtract(faces, adj_tets_faces), 1,
                                         moab::Interface::UNION);
          auto adj_tets_edges = get_adj(tets, 1, moab::Interface::UNION);
          auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
          auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
          adj_faces_edges.merge(geometry_edges);    // geometry edges
          adj_faces_edges.merge(front_block_edges); // front block edges
 
          auto boundary_tets_edges = intersect(adj_tets_edges, adj_faces_edges);
          auto boundary_test_nodes = get_nodes(boundary_tets_edges);
          auto boundary_test_nodes_edges =
              get_adj(boundary_test_nodes, 1, moab::Interface::UNION);
          auto boundary_test_nodes_edges_nodes = subtract(
              get_nodes(boundary_test_nodes_edges), boundary_test_nodes);
 
          boundary_tets_edges =
              subtract(boundary_test_nodes_edges,
                       get_adj(boundary_test_nodes_edges_nodes, 1,
                               moab::Interface::UNION));
 
          Range body_ents;
          CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                             body_ents);
          auto body_skin = get_skin(mField, body_ents);
 
          auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
          body_skin_edges = intersect(get_adj(tets, 1, moab::Interface::UNION),
                                      body_skin_edges);
          body_skin = intersect(body_skin, adj_tets_faces);
          body_skin_edges = subtract(
              body_skin_edges, get_adj(body_skin, 1, moab::Interface::UNION));
 
          save_range(mField.get_moab(), "body_skin_edges.vtk", body_skin_edges);
          for (auto e : body_skin_edges) {
            auto adj_tet = intersect(
                get_adj(Range(e, e), 3, moab::Interface::INTERSECT), tets);
            if (adj_tet.size() == 1) {
              boundary_tets_edges.insert(e);
            }
          }
 
          return boundary_tets_edges;
        };
 
        auto p = get_test_on_crack_surface();
        auto &[new_crack_faces, crack_faces_tets] = p;
 
        if (debug) {
          CHKERR save_range(mField.get_moab(), "hole_crack_faces_debug.vtk",
                            crack_faces);
          CHKERR save_range(mField.get_moab(), "new_crack_faces_debug.vtk",
                            new_crack_faces);
          CHKERR save_range(mField.get_moab(), "new_crack_tets_debug.vtk",
                            crack_faces_tets);
        }
 
        auto boundary_tets_edges =
            carck_faces_test_edges(new_crack_faces, crack_faces_tets);
        CHKERR save_range(mField.get_moab(), "boundary_tets_edges.vtk",
                          boundary_tets_edges);
 
        auto resolve_surface = [&](auto boundary_tets_edges,
                                   auto crack_faces_tets) {
          auto boundary_tets_edges_nodes = get_nodes(boundary_tets_edges);
          auto crack_faces_tets_faces =
              get_adj(crack_faces_tets, 2, moab::Interface::UNION);
 
          Range all_removed_faces;
          Range all_removed_tets;
          int counter = 0;
 
          int size = 0;
          while (size != crack_faces_tets.size()) {
            auto tets_faces =
                get_adj(crack_faces_tets, 2, moab::Interface::UNION);
            auto skin_tets = get_skin(mField, crack_faces_tets);
            auto skin_skin =
                get_skin(mField, subtract(crack_faces_tets_faces, tets_faces));
            auto skin_skin_nodes = get_nodes(skin_skin);
 
            size = crack_faces_tets.size();
            MOFEM_LOG("SELF", Sev::inform)
                << "Crack faces tets size " << crack_faces_tets.size()
                << " crack faces size " << crack_faces_tets_faces.size();
            auto skin_tets_nodes = subtract(
                get_nodes(skin_tets),
                boundary_tets_edges_nodes); // not remove tets which are
                                            // adjagasent to crack faces nodes
            skin_tets_nodes = subtract(skin_tets_nodes, skin_skin_nodes);
 
            Range removed_nodes;
            Range tets_to_remove;
            Range faces_to_remove;
            for (auto n : skin_tets_nodes) {
              auto tets =
                  intersect(get_adj(Range(n, n), 3, moab::Interface::INTERSECT),
                            crack_faces_tets);
              if (tets.size() == 0) {
                continue;
              }
 
              auto hole_detetction = [&]() {
                auto adj_tets =
                    get_adj(Range(n, n), 3, moab::Interface::INTERSECT);
                adj_tets =
                    subtract(adj_tets,
                             crack_faces_tets); // tetst adjacent to the node
                                                // but not part of crack surface
                if (adj_tets.size() == 0) {
                  return std::make_pair(
                      intersect(
                          get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
                          tets_faces),
                      tets);
                }
 
                std::vector<Range> tets_groups;
                auto test_adj_tets = adj_tets;
                while (test_adj_tets.size()) {
                  auto seed_size = 0;
                  Range seed = Range(test_adj_tets[0], test_adj_tets[0]);
                  while (seed.size() != seed_size) {
                    auto adj_faces =
                        subtract(get_adj(seed, 2, moab::Interface::UNION),
                                 tets_faces); // edges which are not
                                              // part of the node
                    seed_size = seed.size();
                    seed.merge(
                        intersect(get_adj(adj_faces, 3, moab::Interface::UNION),
                                  test_adj_tets));
                  }
                  tets_groups.push_back(seed);
                  test_adj_tets = subtract(test_adj_tets, seed);
                }
                if (tets_groups.size() == 1) {
 
                  return std::make_pair(
                      intersect(
                          get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
                          tets_faces),
                      tets);
                }
 
                Range tets_to_remove;
                Range faces_to_remove;
                for (auto &r : tets_groups) {
                  auto f = get_adj(r, 2, moab::Interface::UNION);
                  auto t = intersect(get_adj(f, 3, moab::Interface::UNION),
                                     crack_faces_tets); // tets
 
                  if (f.size() > faces_to_remove.size() ||
                      faces_to_remove.size() == 0) {
                    faces_to_remove = f;
                    tets_to_remove = t; // largest group of tets
                  }
                }
                MOFEM_LOG("EPSELF", Sev::inform)
                    << "Hole detection: faces to remove "
                    << faces_to_remove.size() << " tets to remove "
                    << tets_to_remove.size();
                return std::make_pair(faces_to_remove, tets_to_remove);
              };
 
              if (tets.size() < tets_to_remove.size() ||
                  tets_to_remove.size() == 0) {
                removed_nodes = Range(n, n);
                auto [h_faces_to_remove, h_tets_to_remove] =
                    hole_detetction(); // find faces and tets to remove
                faces_to_remove = h_faces_to_remove;
                tets_to_remove = h_tets_to_remove;
 
                // intersect(
                //     get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
                //     tets_faces);
 
              } // find tets which is largest adjacencty size, so that it is
                // removed first, and then faces are removed
              all_removed_faces.merge(faces_to_remove);
              all_removed_tets.merge(tets_to_remove);
            }
 
            crack_faces_tets = subtract(crack_faces_tets, tets_to_remove);
            crack_faces_tets_faces =
                subtract(crack_faces_tets_faces, faces_to_remove);
 
            if (debug) {
              save_range(mField.get_moab(),
                         "removed_nodes_" +
                             boost::lexical_cast<std::string>(counter) + ".vtk",
                         removed_nodes);
              save_range(mField.get_moab(),
                         "faces_to_remove_" +
                             boost::lexical_cast<std::string>(counter) + ".vtk",
                         faces_to_remove);
              save_range(mField.get_moab(),
                         "tets_to_remove_" +
                             boost::lexical_cast<std::string>(counter) + ".vtk",
                         tets_to_remove);
              save_range(mField.get_moab(),
                         "crack_faces_tets_faces_" +
                             boost::lexical_cast<std::string>(counter) + ".vtk",
                         crack_faces_tets_faces);
              save_range(mField.get_moab(),
                         "crack_faces_tets_" +
                             boost::lexical_cast<std::string>(counter) + ".vtk",
                         crack_faces_tets);
            }
            counter++;
          }
 
          auto cese_internal_faces = [&]() {
            MoFEMFunctionBegin;
            auto skin_tets = get_skin(mField, crack_faces_tets);
            auto adj_faces = get_adj(skin_tets, 2, moab::Interface::UNION);
            adj_faces =
                subtract(adj_faces, skin_tets); // remove skin tets faces
            auto adj_tets = get_adj(adj_faces, 3,
                                    moab::Interface::UNION); // tets which are
                                                             // adjacent to skin
            crack_faces_tets =
                subtract(crack_faces_tets,
                         adj_tets); // remove tets which are adjacent to
                                    // skin, so that they are not removed
            crack_faces_tets_faces =
                subtract(crack_faces_tets_faces, adj_faces);
 
            all_removed_faces.merge(adj_faces);
            all_removed_tets.merge(adj_tets);
 
            MOFEM_LOG("EPSELF", Sev::inform)
                << "Remove internal faces size " << adj_faces.size()
                << " tets size " << adj_tets.size();
            MoFEMFunctionReturn(0);
          };
 
          auto case_only_one_free_edge = [&]() {
            MoFEMFunctionBegin;
 
            for (auto t : Range(crack_faces_tets)) {
 
              auto adj_faces = get_adj(
                  Range(t, t), 2,
                  moab::Interface::UNION); // faces of tet which can be removed
              auto crack_surface_edges =
                  get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
                                   adj_faces),
                          1,
                          moab::Interface::UNION); // edges not on the tet but
                                                   // on crack surface
              auto adj_edges =
                  subtract(get_adj(Range(t, t), 1, moab::Interface::INTERSECT),
                           crack_surface_edges); // free edges
              adj_edges = subtract(
                  adj_edges,
                  boundary_tets_edges); // edges which are not part of gemetry
 
              if (adj_edges.size() == 1) {
                crack_faces_tets =
                    subtract(crack_faces_tets,
                             Range(t, t)); // remove tets which are adjacent to
                                           // skin, so that they are not removed
 
                auto faces_to_remove =
                    get_adj(adj_edges, 2, moab::Interface::UNION); // faces
                                                                   // which can
                                                                   // be removed
                crack_faces_tets_faces =
                    subtract(crack_faces_tets_faces, faces_to_remove);
 
                all_removed_faces.merge(faces_to_remove);
                all_removed_tets.merge(Range(t, t));
 
                MOFEM_LOG("EPSELF", Sev::inform) << "Remove free one edges ";
              }
            }
 
            crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
            crack_faces_tets_faces =
                subtract(crack_faces_tets_faces, all_removed_faces);
 
            MoFEMFunctionReturn(0);
          };
 
          auto cese_flat_tet = [&](auto max_adj_edges) {
            MoFEMFunctionBegin;
 
            Range body_ents;
            CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                               body_ents);
            auto body_skin = get_skin(mField, body_ents);
            auto body_skin_edges =
                get_adj(body_skin, 1, moab::Interface::UNION);
 
            for (auto t : Range(crack_faces_tets)) {
 
              auto adj_faces = get_adj(
                  Range(t, t), 2,
                  moab::Interface::UNION); // faces of tet which can be removed
              auto crack_surface_edges =
                  get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
                                   adj_faces),
                          1,
                          moab::Interface::UNION); // edges not on the tet but
                                                   // on crack surface
              auto adj_edges =
                  subtract(get_adj(Range(t, t), 1, moab::Interface::INTERSECT),
                           crack_surface_edges); // free edges
              adj_edges = subtract(adj_edges, body_skin_edges);
 
              auto tet_edges = get_adj(Range(t, t), 1,
                                       moab::Interface::UNION); // edges of
                                                                // tet
              tet_edges = subtract(tet_edges, adj_edges);
 
              for (auto e : tet_edges) {
                constexpr int opposite_edge[] = {5, 3, 4, 1, 2, 0};
                auto get_side = [&](auto e) {
                  int side, sense, offset;
                  CHK_MOAB_THROW(
                      mField.get_moab().side_number(t, e, side, sense, offset),
                      "get side number failed");
                  return side;
                };
                auto get_side_ent = [&](auto side) {
                  EntityHandle side_edge;
                  CHK_MOAB_THROW(
                      mField.get_moab().side_element(t, 1, side, side_edge),
                      "get side failed");
                  return side_edge;
                };
                adj_edges.erase(get_side_ent(opposite_edge[get_side(e)]));
              }
 
              if (adj_edges.size() <= max_adj_edges) {
 
                double dot = 1;
                Range faces_to_remove;
                for (auto e : adj_edges) {
                  auto edge_adj_faces =
                      get_adj(Range(e, e), 2, moab::Interface::UNION);
                  edge_adj_faces = intersect(edge_adj_faces, adj_faces);
                  if (edge_adj_faces.size() != 2) {
                    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                                      "Adj faces size is not 2 for edge " +
                                          boost::lexical_cast<std::string>(e));
                  }
 
                  auto get_normal = [&](auto f) {
                    FTensor::Tensor1<double, SPACE_DIM> t_n;
                    CHK_THROW_MESSAGE(
                        mField.getInterface<Tools>()->getTriNormal(f, &t_n(0)),
                        "get tri normal failed");
                    return t_n;
                  };
                  auto t_n0 = get_normal(edge_adj_faces[0]);
                  auto t_n1 = get_normal(edge_adj_faces[1]);
                  auto get_sense = [&](auto f) {
                    int side, sense, offset;
                    CHK_MOAB_THROW(mField.get_moab().side_number(t, f, side,
                                                                 sense, offset),
                                   "get side number failed");
                    return sense;
                  };
                  auto sense0 = get_sense(edge_adj_faces[0]);
                  auto sense1 = get_sense(edge_adj_faces[1]);
                  t_n0.normalize();
                  t_n1.normalize();
 
                  FTENSOR_INDEX(SPACE_DIM, i);
                  auto dot_e = (sense0 * sense1) * t_n0(i) * t_n1(i);
                  if (dot_e < dot || e == adj_edges[0]) {
                    dot = dot_e;
                    faces_to_remove = edge_adj_faces;
                  }
                }
 
                all_removed_faces.merge(faces_to_remove);
                all_removed_tets.merge(Range(t, t));
 
                MOFEM_LOG("EPSELF", Sev::inform)
                    << "Remove free edges on flat tet, with considered nb. of "
                       "edges "
                    << adj_edges.size();
              }
            }
 
            crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
            crack_faces_tets_faces =
                subtract(crack_faces_tets_faces, all_removed_faces);
 
            MoFEMFunctionReturn(0);
          };
 
          CHK_THROW_MESSAGE(case_only_one_free_edge(),
                            "Case only one free edge failed");
          for (auto max_adj_edges : {0, 1, 2, 3}) {
            CHK_THROW_MESSAGE(cese_flat_tet(max_adj_edges),
                              "Case only one free edge failed");
          }
          CHK_THROW_MESSAGE(cese_internal_faces(),
                            "Case internal faces failed");
 
          if (debug) {
            save_range(mField.get_moab(),
                       "crack_faces_tets_faces_" +
                           boost::lexical_cast<std::string>(counter) + ".vtk",
                       crack_faces_tets_faces);
            save_range(mField.get_moab(),
                       "crack_faces_tets_" +
                           boost::lexical_cast<std::string>(counter) + ".vtk",
                       crack_faces_tets);
          }
 
          return std::make_tuple(crack_faces_tets_faces, crack_faces_tets,
                                 all_removed_faces, all_removed_tets);
        };
 
        auto [resolved_faces, resolved_tets, all_removed_faces,
              all_removed_tets] =
            resolve_surface(boundary_tets_edges, crack_faces_tets);
        resolved_faces.merge(subtract(crack_faces, all_removed_faces));
        if (debug) {
          CHKERR save_range(mField.get_moab(), "resolved_faces.vtk",
                            resolved_faces);
          CHKERR save_range(mField.get_moab(), "resolved_tets.vtk",
                            resolved_tets);
        }
 
        crack_faces = resolved_faces;
      }
 
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(impl(), "resolve new crack surfaces");
 
    return crack_faces; // send_type(mField, crack_faces, MBTRI);
  };
 
  auto resolve_consisten_crack_extension = [&]() {
    MoFEMFunctionBegin;
    auto crack_meshset =
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
            addCrackMeshsetId, BLOCKSET);
    auto meshset = crack_meshset->getMeshset();
 
    if (!mField.get_comm_rank() && !noCrackExtension) {
      Range old_crack_faces;
      CHKERR mField.get_moab().get_entities_by_type(meshset, MBTRI,
                                                    old_crack_faces);
      auto extendeded_crack_faces = get_extended_crack_faces();
      auto reconstructed_crack_faces =
          subtract(reconstruct_crack_faces(extendeded_crack_faces),
                   subtract(*crackFaces, old_crack_faces));
      if (nbCrackFaces >= reconstructed_crack_faces.size()) {
        MOFEM_LOG("EPSELF", Sev::warning)
            << "No new crack faces to add, skipping adding to meshset";
        extendeded_crack_faces = subtract(
            extendeded_crack_faces, subtract(*crackFaces, old_crack_faces));
        MOFEM_LOG("EPSELF", Sev::inform)
            << "Number crack faces size (extended) "
            << extendeded_crack_faces.size();
        CHKERR mField.get_moab().clear_meshset(&meshset, 1);
        CHKERR mField.get_moab().add_entities(meshset, extendeded_crack_faces);
      } else {
        CHKERR mField.get_moab().clear_meshset(&meshset, 1);
        CHKERR mField.get_moab().add_entities(meshset,
                                              reconstructed_crack_faces);
        MOFEM_LOG("EPSELF", Sev::inform)
            << "Number crack faces size (reconstructed) "
            << reconstructed_crack_faces.size();
        nbCrackFaces = reconstructed_crack_faces.size();
      }
    }
 
    Range crack_faces;
    if (!mField.get_comm_rank()) {
      CHKERR mField.get_moab().get_entities_by_type(meshset, MBTRI,
                                                    crack_faces);
    }
    crack_faces = send_type(mField, crack_faces, MBTRI);
    if (mField.get_comm_rank()) {
      CHKERR mField.get_moab().clear_meshset(&meshset, 1);
      CHKERR mField.get_moab().add_entities(meshset, crack_faces);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR resolve_consisten_crack_extension();
 
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode EshelbianCore::createCrackSurfaceMeshset() {
  MoFEMFunctionBegin;
  auto meshset_mng = mField.getInterface<MeshsetsManager>();
  while (meshset_mng->checkMeshset(addCrackMeshsetId, BLOCKSET))
    ++addCrackMeshsetId;
  MOFEM_LOG("EP", Sev::inform)
      << "Crack added surface meshset " << addCrackMeshsetId;
  CHKERR meshset_mng->addMeshset(BLOCKSET, addCrackMeshsetId, "CRACK_COMPUTED");
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode
EshelbianCore::calculateCrackArea(boost::shared_ptr<double> &area_ptr) {
  MoFEMFunctionBegin;
 
  if (!area_ptr) {
    // initialize area
    area_ptr = boost::shared_ptr<double>(new double(0.0));
  }
 
  int success;
  *area_ptr = 0;
  if (mField.get_comm_rank() == 0) {
    MOFEM_LOG("EP", Sev::inform) << "Calculate crack area";
    auto crack_faces = get_range_from_block(mField, "CRACK", SPACE_DIM - 1);
    for (auto f : crack_faces) {
      *area_ptr += mField.getInterface<Tools>()->getTriArea(f);
    }
    success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0, mField.get_comm());
  } else {
    success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0, mField.get_comm());
  }
  if (success != MPI_SUCCESS) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "bcast failed");
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity