EshelbianPlasticity.cpp

Go to the documentation of this file.

/**
 * \file EshelbianPlasticity.cpp
 * \example mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp
 *
 * \brief Eshelbian plasticity implementation
 *
 * \copyright 2024. Various authors, some of them anonymous contributors under
 * MiT core contributors license agreement.
 */
 
#define SINGULARITY
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <CGGTonsorialBubbleBase.hpp>
 
#include <EshelbianPlasticity.hpp>
#include <boost/math/constants/constants.hpp>
 
#include <cholesky.hpp>
#ifdef ENABLE_PYTHON_BINDING
  #include <boost/python.hpp>
  #include <boost/python/def.hpp>
  #include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
 
  #pragma message "With ENABLE_PYTHON_BINDING"
 
#else
 
  #pragma message "Without ENABLE_PYTHON_BINDING"
 
#endif
 
#include <EshelbianAux.hpp>
#include <EshelbianContact.hpp>
 
extern "C" {
#include <phg-quadrule/quad.h>
}
 
#include <queue>
 
namespace EshelbianPlasticity {
struct VolUserDataOperatorStabAssembly
    : public VolumeElementForcesAndSourcesCore::UserDataOperator {
  using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
struct FaceUserDataOperatorStabAssembly
    : public FaceElementForcesAndSourcesCore::UserDataOperator {
  using FaceElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
 
} // namespace EshelbianPlasticity
 
static auto send_type(MoFEM::Interface &m_field, Range r,
                      const EntityType type) {
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
 
  auto dim = CN::Dimension(type);
 
  std::vector<int> sendcounts(pcomm->size());
  std::vector<int> displs(pcomm->size());
  std::vector<int> sendbuf(r.size());
  if (pcomm->rank() == 0) {
    for (auto p = 1; p != pcomm->size(); p++) {
      auto part_ents = m_field.getInterface<CommInterface>()
                           ->getPartEntities(m_field.get_moab(), p)
                           .subset_by_dimension(SPACE_DIM);
      Range faces;
      CHKERR m_field.get_moab().get_adjacencies(part_ents, dim, true, faces,
                                                moab::Interface::UNION);
      faces = intersect(faces, r);
      sendcounts[p] = faces.size();
      displs[p] = sendbuf.size();
      for (auto f : faces) {
        auto id = id_from_handle(f);
        sendbuf.push_back(id);
      }
    }
  }
 
  int recv_data;
  MPI_Scatter(sendcounts.data(), 1, MPI_INT, &recv_data, 1, MPI_INT, 0,
              pcomm->comm());
  std::vector<int> recvbuf(recv_data);
  MPI_Scatterv(sendbuf.data(), sendcounts.data(), displs.data(), MPI_INT,
               recvbuf.data(), recv_data, MPI_INT, 0, pcomm->comm());
 
  if (pcomm->rank() > 0) {
    Range r;
    for (auto &f : recvbuf) {
      r.insert(ent_form_type_and_id(type, f));
    }
    return r;
  }
 
  return r;
}
 
static auto get_range_from_block(MoFEM::Interface &m_field,
                                 const std::string block_name, int dim) {
  Range r;
 
  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
  auto bcs = mesh_mng->getCubitMeshsetPtr(
 
      std::regex((boost::format("%s(.*)") % block_name).str())
 
  );
 
  for (auto bc : bcs) {
    Range faces;
    CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
                                                       faces, true),
                   "get meshset ents");
    r.merge(faces);
  }
 
  return r;
};
 
static auto get_range_from_block_map(MoFEM::Interface &m_field,
                                     const std::string block_name, int dim) {
  std::map<std::string, Range> r;
 
  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
  auto bcs = mesh_mng->getCubitMeshsetPtr(
 
      std::regex((boost::format("%s(.*)") % block_name).str())
 
  );
 
  for (auto bc : bcs) {
    Range faces;
    CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
                                                       faces, true),
                   "get meshset ents");
    r[bc->getName()] = faces;
  }
 
  return r;
}
 
static auto get_block_meshset(MoFEM::Interface &m_field, const int ms_id,
                              const unsigned int cubit_bc_type) {
  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
  EntityHandle meshset;
  CHKERR mesh_mng->getMeshset(ms_id, cubit_bc_type, meshset);
  return meshset;
};
 
static auto save_range(moab::Interface &moab, const std::string name,
                       const Range r, std::vector<Tag> tags = {}) {
  MoFEMFunctionBegin;
  auto out_meshset = get_temp_meshset_ptr(moab);
  CHKERR moab.add_entities(*out_meshset, r);
  if (r.size()) {
    CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1,
                           tags.data(), tags.size());
  } else {
    MOFEM_LOG("SELF", Sev::warning) << "Empty range for " << name;
  }
  MoFEMFunctionReturn(0);
};
 
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin) {
  Range boundary_ents;
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
  CHK_MOAB_THROW(pcomm->filter_pstatus(skin,
                                       PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                       PSTATUS_NOT, -1, &boundary_ents),
                 "filter_pstatus");
  return boundary_ents;
};
 
static auto filter_owners(MoFEM::Interface &m_field, Range skin) {
  Range owner_ents;
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
  CHK_MOAB_THROW(pcomm->filter_pstatus(skin, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1,
                                       &owner_ents),
                 "filter_pstatus");
  return owner_ents;
};
 
static auto get_skin(MoFEM::Interface &m_field, Range body_ents) {
  Skinner skin(&m_field.get_moab());
  Range skin_ents;
  CHK_MOAB_THROW(skin.find_skin(0, body_ents, false, skin_ents), "find_skin");
  return skin_ents;
};
 
static auto get_crack_front_edges(MoFEM::Interface &m_field,
                                  Range crack_faces) {
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
  auto &moab = m_field.get_moab();
  Range crack_skin_without_bdy;
  if (pcomm->rank() == 0) {
    Range crack_edges;
    CHKERR moab.get_adjacencies(crack_faces, 1, true, crack_edges,
                                moab::Interface::UNION);
    auto crack_skin = get_skin(m_field, crack_faces);
    Range body_ents;
    CHK_MOAB_THROW(
        m_field.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents),
        "get_entities_by_dimension");
    auto body_skin = get_skin(m_field, body_ents);
    Range body_skin_edges;
    CHK_MOAB_THROW(moab.get_adjacencies(body_skin, 1, true, body_skin_edges,
                                        moab::Interface::UNION),
                   "get_adjacencies");
    crack_skin_without_bdy = subtract(crack_skin, body_skin_edges);
    auto front_edges_map = get_range_from_block_map(m_field, "FRONT", 1);
    for (auto &m : front_edges_map) {
      auto add_front = subtract(m.second, crack_edges);
      auto i = intersect(m.second, crack_edges);
      if (i.empty()) {
        crack_skin_without_bdy.merge(add_front);
      } else {
        auto i_skin = get_skin(m_field, i);
        Range adj_i_skin;
        CHKERR moab.get_adjacencies(i_skin, 1, true, adj_i_skin,
                                    moab::Interface::UNION);
        adj_i_skin = subtract(intersect(adj_i_skin, m.second), crack_edges);
        crack_skin_without_bdy.merge(adj_i_skin);
      }
    }
  }
  return send_type(m_field, crack_skin_without_bdy, MBEDGE);
}
 
static auto get_two_sides_of_crack_surface(MoFEM::Interface &m_field,
                                           Range crack_faces) {
 
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
 
  MOFEM_LOG("EP", Sev::noisy) << "get_two_sides_of_crack_surface";
 
  if (!pcomm->rank()) {
 
    auto impl = [&](auto &saids) {
      MoFEMFunctionBegin;
 
      auto &moab = m_field.get_moab();
 
      auto get_adj = [&](auto e, auto dim) {
        Range adj;
        CHK_MOAB_THROW(m_field.get_moab().get_adjacencies(
                           e, dim, true, adj, moab::Interface::UNION),
                       "get adj");
        return adj;
      };
 
      auto get_conn = [&](auto e) {
        Range conn;
        CHK_MOAB_THROW(m_field.get_moab().get_connectivity(e, conn, true),
                       "get connectivity");
        return conn;
      };
 
      constexpr bool debug = false;
      Range body_ents;
      CHKERR m_field.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                          body_ents);
      auto body_skin = get_skin(m_field, body_ents);
      auto body_skin_edges = get_adj(body_skin, 1);
 
      auto crack_skin =
          subtract(get_skin(m_field, crack_faces), body_skin_edges);
      auto crack_skin_conn = get_conn(crack_skin);
      auto crack_skin_conn_edges = get_adj(crack_skin_conn, 1);
      auto crack_edges = get_adj(crack_faces, 1);
      crack_edges = subtract(crack_edges, crack_skin);
      auto all_tets = get_adj(crack_edges, 3);
      crack_edges = subtract(crack_edges, crack_skin_conn_edges);
      auto crack_conn = get_conn(crack_edges);
      all_tets.merge(get_adj(crack_conn, 3));
 
      if (debug) {
        CHKERR save_range(m_field.get_moab(), "crack_faces.vtk", crack_faces);
        CHKERR save_range(m_field.get_moab(), "all_crack_tets.vtk", all_tets);
        CHKERR save_range(m_field.get_moab(), "crack_edges_all.vtk",
                          crack_edges);
      }
 
      if (crack_faces.size()) {
        auto grow = [&](auto r) {
          auto crack_faces_conn = get_conn(crack_faces);
          Range v;
          auto size_r = 0;
          while (size_r != r.size() && r.size() > 0) {
            size_r = r.size();
            CHKERR moab.get_connectivity(r, v, true);
            v = subtract(v, crack_faces_conn);
            if (v.size()) {
              CHKERR moab.get_adjacencies(v, SPACE_DIM, true, r,
                                          moab::Interface::UNION);
              r = intersect(r, all_tets);
            }
            if (r.empty()) {
              break;
            }
          }
          return r;
        };
 
        Range all_tets_ord = all_tets;
        while (all_tets.size()) {
          Range faces = get_adj(unite(saids.first, saids.second), 2);
          faces = subtract(crack_faces, faces);
          if (faces.size()) {
            Range tets;
            auto fit = faces.begin();
            for (; fit != faces.end(); ++fit) {
              tets = intersect(get_adj(Range(*fit, *fit), 3), all_tets);
              if (tets.size() == 2) {
                break;
              }
            }
            if (tets.empty()) {
              break;
            } else {
              saids.first.insert(tets[0]);
              saids.first = grow(saids.first);
              all_tets = subtract(all_tets, saids.first);
              if (tets.size() == 2) {
                saids.second.insert(tets[1]);
                saids.second = grow(saids.second);
                all_tets = subtract(all_tets, saids.second);
              }
            }
          } else {
            break;
          }
        }
 
        saids.first = subtract(all_tets_ord, saids.second);
        saids.second = subtract(all_tets_ord, saids.first);
      }
 
      MoFEMFunctionReturn(0);
    };
 
    std::pair<Range, Range> saids;
    if (crack_faces.size())
      CHK_THROW_MESSAGE(impl(saids), "get crack both sides");
    return saids;
  }
 
  MOFEM_LOG("EP", Sev::noisy) << "get_two_sides_of_crack_surface <- done";
 
  return std::pair<Range, Range>();
}
 
namespace EshelbianPlasticity {
 
struct SetIntegrationAtFrontVolume {
 
  SetIntegrationAtFrontVolume(boost::shared_ptr<Range> front_nodes,
                              boost::shared_ptr<Range> front_edges)
      : frontNodes(front_nodes), frontEdges(front_edges){};
 
  MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                            int order_col, int order_data) {
    MoFEMFunctionBegin;
 
    constexpr bool debug = false;
 
    constexpr int numNodes = 4;
    constexpr int numEdges = 6;
    constexpr int refinementLevels = 6;
 
    auto &m_field = fe_raw_ptr->mField;
    auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
    auto fe_handle = fe_ptr->getFEEntityHandle();
 
    auto set_base_quadrature = [&]() {
      MoFEMFunctionBegin;
      int rule = 2 * order_data + 1;
      if (rule < QUAD_3D_TABLE_SIZE) {
        if (QUAD_3D_TABLE[rule]->dim != 3) {
          SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
                  "wrong dimension");
        }
        if (QUAD_3D_TABLE[rule]->order < rule) {
          SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
                   "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
        }
        const size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
        auto &gauss_pts = fe_ptr->gaussPts;
        gauss_pts.resize(4, nb_gauss_pts, false);
        cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
                    &gauss_pts(0, 0), 1);
        cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
                    &gauss_pts(1, 0), 1);
        cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
                    &gauss_pts(2, 0), 1);
        cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
                    &gauss_pts(3, 0), 1);
        auto &data = fe_ptr->dataOnElement[H1];
        data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
                                                              false);
        double *shape_ptr =
            &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
        cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
                    1);
      } else {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_base_quadrature();
 
    if (frontNodes->size() && EshelbianCore::setSingularity) {
 
      auto get_singular_nodes = [&]() {
        int num_nodes;
        const EntityHandle *conn;
        CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
                                                   true);
        std::bitset<numNodes> singular_nodes;
        for (auto nn = 0; nn != numNodes; ++nn) {
          if (frontNodes->find(conn[nn]) != frontNodes->end()) {
            singular_nodes.set(nn);
          } else {
            singular_nodes.reset(nn);
          }
        }
        return singular_nodes;
      };
 
      auto get_singular_edges = [&]() {
        std::bitset<numEdges> singular_edges;
        for (int ee = 0; ee != numEdges; ee++) {
          EntityHandle edge;
          CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
          if (frontEdges->find(edge) != frontEdges->end()) {
            singular_edges.set(ee);
          } else {
            singular_edges.reset(ee);
          }
        }
        return singular_edges;
      };
 
      auto set_gauss_pts = [&](auto &ref_gauss_pts) {
        MoFEMFunctionBegin;
        fe_ptr->gaussPts.swap(ref_gauss_pts);
        const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
        auto &data = fe_ptr->dataOnElement[H1];
        data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
        double *shape_ptr =
            &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
        CHKERR ShapeMBTET(shape_ptr, &fe_ptr->gaussPts(0, 0),
                          &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
                          nb_gauss_pts);
        MoFEMFunctionReturn(0);
      };
 
      auto singular_nodes = get_singular_nodes();
      if (singular_nodes.count()) {
        auto it_map_ref_coords = mapRefCoords.find(singular_nodes.to_ulong());
        if (it_map_ref_coords != mapRefCoords.end()) {
          CHKERR set_gauss_pts(it_map_ref_coords->second);
          MoFEMFunctionReturnHot(0);
        } else {
 
          auto refine_quadrature = [&]() {
            MoFEMFunctionBegin;
 
            const int max_level = refinementLevels;
            EntityHandle tet;
 
            moab::Core moab_ref;
            double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
            EntityHandle nodes[4];
            for (int nn = 0; nn != 4; nn++)
              CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
            CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
            MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
            MoFEM::Interface &m_field_ref = mofem_ref_core;
            {
              Range tets(tet, tet);
              Range edges;
              CHKERR m_field_ref.get_moab().get_adjacencies(
                  tets, 1, true, edges, moab::Interface::UNION);
              CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
                  tets, BitRefLevel().set(0), false, VERBOSE);
            }
 
            Range nodes_at_front;
            for (int nn = 0; nn != numNodes; nn++) {
              if (singular_nodes[nn]) {
                EntityHandle ent;
                CHKERR moab_ref.side_element(tet, 0, nn, ent);
                nodes_at_front.insert(ent);
              }
            }
 
            auto singular_edges = get_singular_edges();
 
            EntityHandle meshset;
            CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
            for (int ee = 0; ee != numEdges; ee++) {
              if (singular_edges[ee]) {
                EntityHandle ent;
                CHKERR moab_ref.side_element(tet, 1, ee, ent);
                CHKERR moab_ref.add_entities(meshset, &ent, 1);
              }
            }
 
            // refine mesh
            auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
            for (int ll = 0; ll != max_level; ll++) {
              Range edges;
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                                 BitRefLevel().set(), MBEDGE,
                                                 edges);
              Range ref_edges;
              CHKERR moab_ref.get_adjacencies(
                  nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
              ref_edges = intersect(ref_edges, edges);
              Range ents;
              CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
              ref_edges = intersect(ref_edges, ents);
              Range tets;
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->getEntitiesByTypeAndRefLevel(
                      BitRefLevel().set(ll), BitRefLevel().set(), MBTET, tets);
              CHKERR m_ref->addVerticesInTheMiddleOfEdges(
                  ref_edges, BitRefLevel().set(ll + 1));
              CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->updateMeshsetByEntitiesChildren(meshset,
                                                    BitRefLevel().set(ll + 1),
                                                    meshset, MBEDGE, true);
            }
 
            // get ref coords
            Range tets;
            CHKERR m_field_ref.getInterface<BitRefManager>()
                ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
                                               BitRefLevel().set(), MBTET,
                                               tets);
 
            if (debug) {
              CHKERR save_range(moab_ref, "ref_tets.vtk", tets);
            }
 
            MatrixDouble ref_coords(tets.size(), 12, false);
            int tt = 0;
            for (Range::iterator tit = tets.begin(); tit != tets.end();
                 tit++, tt++) {
              int num_nodes;
              const EntityHandle *conn;
              CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
              CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
            }
 
            auto &data = fe_ptr->dataOnElement[H1];
            const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
            MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
            MatrixDouble &shape_n =
                data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
            int gg = 0;
            for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
              double *tet_coords = &ref_coords(tt, 0);
              double det = Tools::tetVolume(tet_coords);
              det *= 6;
              for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
                for (int dd = 0; dd != 3; dd++) {
                  ref_gauss_pts(dd, gg) =
                      shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
                      shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
                      shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
                      shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
                }
                ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
              }
            }
 
            mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
            CHKERR set_gauss_pts(mapRefCoords[singular_nodes.to_ulong()]);
 
            MoFEMFunctionReturn(0);
          };
 
          CHKERR refine_quadrature();
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  struct Fe : public ForcesAndSourcesCore {
    using ForcesAndSourcesCore::dataOnElement;
 
  private:
    using ForcesAndSourcesCore::ForcesAndSourcesCore;
  };
 
  boost::shared_ptr<Range> frontNodes;
  boost::shared_ptr<Range> frontEdges;
 
  static inline std::map<long int, MatrixDouble> mapRefCoords;
};
 
struct SetIntegrationAtFrontFace {
 
  using FunRule = boost::function<int(int)>;
  FunRule funRule = [](int p) { return 2 * (p + 1); };
 
  SetIntegrationAtFrontFace(boost::shared_ptr<Range> front_nodes,
                            boost::shared_ptr<Range> front_edges)
      : frontNodes(front_nodes), frontEdges(front_edges){};
 
  SetIntegrationAtFrontFace(boost::shared_ptr<Range> front_nodes,
                            boost::shared_ptr<Range> front_edges,
                            FunRule fun_rule)
      : frontNodes(front_nodes), frontEdges(front_edges), funRule(fun_rule){};
 
  MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                            int order_col, int order_data) {
    MoFEMFunctionBegin;
 
    constexpr bool debug = false;
 
    constexpr int numNodes = 3;
    constexpr int numEdges = 3;
    constexpr int refinementLevels = 6;
 
    auto &m_field = fe_raw_ptr->mField;
    auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
    auto fe_handle = fe_ptr->getFEEntityHandle();
 
    auto set_base_quadrature = [&]() {
      MoFEMFunctionBegin;
      int rule = funRule(order_data);
      if (rule < QUAD_2D_TABLE_SIZE) {
        if (QUAD_2D_TABLE[rule]->dim != 2) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
        }
        if (QUAD_2D_TABLE[rule]->order < rule) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                   "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
        }
        const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
        fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
        cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
                    &fe_ptr->gaussPts(0, 0), 1);
        cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
                    &fe_ptr->gaussPts(1, 0), 1);
        cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
                    &fe_ptr->gaussPts(2, 0), 1);
        auto &data = fe_ptr->dataOnElement[H1];
        data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
                                                              false);
        double *shape_ptr =
            &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
        cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
                    1);
        data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
        std::copy(
            Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
            data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
 
      } else {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_base_quadrature();
 
    if (frontNodes->size() && EshelbianCore::setSingularity) {
 
      auto get_singular_nodes = [&]() {
        int num_nodes;
        const EntityHandle *conn;
        CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
                                                   true);
        std::bitset<numNodes> singular_nodes;
        for (auto nn = 0; nn != numNodes; ++nn) {
          if (frontNodes->find(conn[nn]) != frontNodes->end()) {
            singular_nodes.set(nn);
          } else {
            singular_nodes.reset(nn);
          }
        }
        return singular_nodes;
      };
 
      auto get_singular_edges = [&]() {
        std::bitset<numEdges> singular_edges;
        for (int ee = 0; ee != numEdges; ee++) {
          EntityHandle edge;
          CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
          if (frontEdges->find(edge) != frontEdges->end()) {
            singular_edges.set(ee);
          } else {
            singular_edges.reset(ee);
          }
        }
        return singular_edges;
      };
 
      auto set_gauss_pts = [&](auto &ref_gauss_pts) {
        MoFEMFunctionBegin;
        fe_ptr->gaussPts.swap(ref_gauss_pts);
        const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
        auto &data = fe_ptr->dataOnElement[H1];
        data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
        double *shape_ptr =
            &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
        CHKERR ShapeMBTRI(shape_ptr, &fe_ptr->gaussPts(0, 0),
                          &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
        MoFEMFunctionReturn(0);
      };
 
      auto singular_nodes = get_singular_nodes();
      if (singular_nodes.count()) {
        auto it_map_ref_coords = mapRefCoords.find(singular_nodes.to_ulong());
        if (it_map_ref_coords != mapRefCoords.end()) {
          CHKERR set_gauss_pts(it_map_ref_coords->second);
          MoFEMFunctionReturnHot(0);
        } else {
 
          auto refine_quadrature = [&]() {
            MoFEMFunctionBegin;
 
            const int max_level = refinementLevels;
 
            moab::Core moab_ref;
            double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
            EntityHandle nodes[numNodes];
            for (int nn = 0; nn != numNodes; nn++)
              CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
            EntityHandle tri;
            CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
            MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
            MoFEM::Interface &m_field_ref = mofem_ref_core;
            {
              Range tris(tri, tri);
              Range edges;
              CHKERR m_field_ref.get_moab().get_adjacencies(
                  tris, 1, true, edges, moab::Interface::UNION);
              CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
                  tris, BitRefLevel().set(0), false, VERBOSE);
            }
 
            Range nodes_at_front;
            for (int nn = 0; nn != numNodes; nn++) {
              if (singular_nodes[nn]) {
                EntityHandle ent;
                CHKERR moab_ref.side_element(tri, 0, nn, ent);
                nodes_at_front.insert(ent);
              }
            }
 
            auto singular_edges = get_singular_edges();
 
            EntityHandle meshset;
            CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
            for (int ee = 0; ee != numEdges; ee++) {
              if (singular_edges[ee]) {
                EntityHandle ent;
                CHKERR moab_ref.side_element(tri, 1, ee, ent);
                CHKERR moab_ref.add_entities(meshset, &ent, 1);
              }
            }
 
            // refine mesh
            auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
            for (int ll = 0; ll != max_level; ll++) {
              Range edges;
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                                 BitRefLevel().set(), MBEDGE,
                                                 edges);
              Range ref_edges;
              CHKERR moab_ref.get_adjacencies(
                  nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
              ref_edges = intersect(ref_edges, edges);
              Range ents;
              CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
              ref_edges = intersect(ref_edges, ents);
              Range tris;
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->getEntitiesByTypeAndRefLevel(
                      BitRefLevel().set(ll), BitRefLevel().set(), MBTRI, tris);
              CHKERR m_ref->addVerticesInTheMiddleOfEdges(
                  ref_edges, BitRefLevel().set(ll + 1));
              CHKERR m_ref->refineTris(tris, BitRefLevel().set(ll + 1));
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->updateMeshsetByEntitiesChildren(meshset,
                                                    BitRefLevel().set(ll + 1),
                                                    meshset, MBEDGE, true);
            }
 
            // get ref coords
            Range tris;
            CHKERR m_field_ref.getInterface<BitRefManager>()
                ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
                                               BitRefLevel().set(), MBTRI,
                                               tris);
 
            if (debug) {
              CHKERR save_range(moab_ref, "ref_tris.vtk", tris);
              SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "debug");
            }
 
            MatrixDouble ref_coords(tris.size(), 9, false);
            int tt = 0;
            for (Range::iterator tit = tris.begin(); tit != tris.end();
                 tit++, tt++) {
              int num_nodes;
              const EntityHandle *conn;
              CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
              CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
            }
 
            auto &data = fe_ptr->dataOnElement[H1];
            const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
            MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
            MatrixDouble &shape_n =
                data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
            int gg = 0;
            for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
              double *tri_coords = &ref_coords(tt, 0);
              FTensor::Tensor1<double, 3> t_normal;
              CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
              auto det = t_normal.l2();
              for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
                for (int dd = 0; dd != 2; dd++) {
                  ref_gauss_pts(dd, gg) =
                      shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
                      shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
                      shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
                }
                ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
              }
            }
 
            mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
            CHKERR set_gauss_pts(mapRefCoords[singular_nodes.to_ulong()]);
 
            MoFEMFunctionReturn(0);
          };
 
          CHKERR refine_quadrature();
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  struct Fe : public ForcesAndSourcesCore {
    using ForcesAndSourcesCore::dataOnElement;
 
  private:
    using ForcesAndSourcesCore::ForcesAndSourcesCore;
  };
 
  boost::shared_ptr<Range> frontNodes;
  boost::shared_ptr<Range> frontEdges;
 
  static inline std::map<long int, MatrixDouble> mapRefCoords;
};
 
double EshelbianCore::exponentBase = exp(1);
boost::function<double(const double)> EshelbianCore::f = EshelbianCore::f_log_e;
boost::function<double(const double)> EshelbianCore::d_f =
    EshelbianCore::d_f_log_e;
boost::function<double(const double)> EshelbianCore::dd_f =
    EshelbianCore::dd_f_log_e;
boost::function<double(const double)> EshelbianCore::inv_f =
    EshelbianCore::inv_f_log_e;
boost::function<double(const double)> EshelbianCore::inv_d_f =
    EshelbianCore::inv_d_f_log_e;
boost::function<double(const double)> EshelbianCore::inv_dd_f =
    EshelbianCore::inv_dd_f_log_e;
 
MoFEMErrorCode
EshelbianCore::query_interface(boost::typeindex::type_index type_index,
                               UnknownInterface **iface) const {
  *iface = const_cast<EshelbianCore *>(this);
  return 0;
}
 
MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
 
  if (evalRhs)
    CHKERR evaluateRhs(data);
 
  if (evalLhs)
    CHKERR evaluateLhs(data);
 
  MoFEMFunctionReturn(0);
}
 
EshelbianCore::EshelbianCore(MoFEM::Interface &m_field) : mField(m_field) {
  CHK_THROW_MESSAGE(getOptions(), "getOptions failed");
}
 
EshelbianCore::~EshelbianCore() = default;
 
MoFEMErrorCode EshelbianCore::getOptions() {
  MoFEMFunctionBegin;
  const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
  const char *list_symm[] = {"symm", "not_symm"};
  const char *list_release[] = {"griffith_force", "griffith_skeleton"};
  const char *list_stretches[] = {"linear", "log", "log_quadratic"};
  PetscInt choice_rot = EshelbianCore::rotSelector;
  PetscInt choice_grad = EshelbianCore::gradApproximator;
  PetscInt choice_symm = EshelbianCore::symmetrySelector;
  PetscInt choice_release = EshelbianCore::energyReleaseSelector;
  PetscInt choice_stretch = StretchSelector::LOG;
  char analytical_expr_file_name[255] = "analytical_expr.py";
 
  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
                           "none");
  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
                         spaceOrder, &spaceOrder, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
                         spaceH1Order, &spaceH1Order, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
                         "", materialOrder, &materialOrder, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
                            &alphaU, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
                            &alphaW, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_omega", "rot viscosity", "",
                            alphaOmega, &alphaOmega, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
                            &alphaRho, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
                           LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
                           PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-grad", "gradient of defamation approximate", "",
                           list_rots, NO_H1_CONFIGURATION + 1,
                           list_rots[choice_grad], &choice_grad, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-symm", "symmetric variant", "", list_symm, 2,
                           list_symm[choice_symm], &choice_symm, PETSC_NULLPTR);
 
  CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
                            &EshelbianCore::exponentBase, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-stretches", "stretches", "", list_stretches,
                           StretchSelector::STRETCH_SELECTOR_LAST,
                           list_stretches[choice_stretch], &choice_stretch,
                           PETSC_NULLPTR);
 
  CHKERR PetscOptionsBool("-no_stretch", "do not solve for stretch", "",
                          noStretch, &noStretch, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-set_singularity", "set singularity", "",
                          setSingularity, &setSingularity, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-l2_user_base_scale", "streach scale", "",
                          l2UserBaseScale, &l2UserBaseScale, PETSC_NULLPTR);
  
  // dynamic relaxation
  CHKERR PetscOptionsBool("-dynamic_relaxation", "dynamic time relaxation", "",
                          dynamicRelaxation, &dynamicRelaxation, PETSC_NULLPTR);
 
  // contact parameters
  CHKERR PetscOptionsInt("-contact_max_post_proc_ref_level", "refinement level",
                         "", contactRefinementLevels, &contactRefinementLevels,
                         PETSC_NULLPTR);
 
  // cracking parameters
  CHKERR PetscOptionsBool("-cracking_on", "cracking ON", "", crackingOn,
                          &crackingOn, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-cracking_start_time", "cracking start time", "",
                            crackingStartTime, &crackingStartTime, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-griffith_energy", "Griffith energy", "",
                            griffithEnergy, &griffithEnergy, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-energy_release_variant", "energy release variant",
                           "", list_release, 2, list_release[choice_release],
                           &choice_release, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-nb_J_integral_levels", "Number of J integarl levels",
                         "", nbJIntegralContours, &nbJIntegralContours,
                         PETSC_NULLPTR); // backward compatibility
  CHKERR PetscOptionsInt(
      "-nb_J_integral_contours", "Number of J integral contours", "",
      nbJIntegralContours, &nbJIntegralContours, PETSC_NULLPTR);
 
  // internal stress
  char tag_name[255] = "";
  CHKERR PetscOptionsString("-internal_stress_tag_name",
                            "internal stress tag name", "", "", tag_name, 255,
                            PETSC_NULLPTR);
  internalStressTagName = string(tag_name);
  CHKERR PetscOptionsInt(
      "-internal_stress_interp_order", "internal stress interpolation order",
      "", internalStressInterpOrder, &internalStressInterpOrder, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-internal_stress_voigt", "Voigt index notation", "",
                          internalStressVoigt, &internalStressVoigt,
                          PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-analytical_expr_file",
                               analytical_expr_file_name, 255, PETSC_NULLPTR);
 
  PetscOptionsEnd();
 
  if (internalStressInterpOrder != 0 && internalStressInterpOrder != 1) {
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
             "Unsupported internal stress interpolation order %d",
             internalStressInterpOrder);
  }
 
  if (setSingularity) {
    l2UserBaseScale = PETSC_TRUE;
  } 
 
  EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
  EshelbianCore::gradApproximator = static_cast<RotSelector>(choice_grad);
  EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
  EshelbianCore::symmetrySelector = static_cast<SymmetrySelector>(choice_symm);
  EshelbianCore::energyReleaseSelector =
      static_cast<EnergyReleaseSelector>(choice_release);
 
  switch (EshelbianCore::stretchSelector) {
  case StretchSelector::LINEAR:
    EshelbianCore::f = f_linear;
    EshelbianCore::d_f = d_f_linear;
    EshelbianCore::dd_f = dd_f_linear;
    EshelbianCore::inv_f = inv_f_linear;
    EshelbianCore::inv_d_f = inv_d_f_linear;
    EshelbianCore::inv_dd_f = inv_dd_f_linear;
    break;
  case StretchSelector::LOG:
    if (std::fabs(EshelbianCore::exponentBase - exp(1)) >
        std::numeric_limits<float>::epsilon()) {
      EshelbianCore::f = f_log;
      EshelbianCore::d_f = d_f_log;
      EshelbianCore::dd_f = dd_f_log;
      EshelbianCore::inv_f = inv_f_log;
      EshelbianCore::inv_d_f = inv_d_f_log;
      EshelbianCore::inv_dd_f = inv_dd_f_log;
    } else {
      EshelbianCore::f = f_log_e;
      EshelbianCore::d_f = d_f_log_e;
      EshelbianCore::dd_f = dd_f_log_e;
      EshelbianCore::inv_f = inv_f_log;
      EshelbianCore::inv_d_f = inv_d_f_log;
      EshelbianCore::inv_dd_f = inv_dd_f_log;
    }
    break;
  case StretchSelector::LOG_QUADRATIC:
    EshelbianCore::f = f_log_e_quadratic;
    EshelbianCore::d_f = d_f_log_e_quadratic;
    EshelbianCore::dd_f = dd_f_log_e_quadratic;
    EshelbianCore::inv_f = [](const double x) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "No logarithmic quadratic stretch for this case");
      return 0;
    };
    EshelbianCore::inv_d_f = EshelbianCore::inv_f;
    EshelbianCore::inv_dd_f = EshelbianCore::inv_f;
    break; // no stretch, do not use stretch functions
  default:
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
    break;
  };
 
  MOFEM_LOG("EP", Sev::inform) << "spaceOrder: -space_order " << spaceOrder;
  MOFEM_LOG("EP", Sev::inform)
      << "spaceH1Order: -space_h1_order " << spaceH1Order;
  MOFEM_LOG("EP", Sev::inform)
      << "materialOrder: -material_order " << materialOrder;
  MOFEM_LOG("EP", Sev::inform) << "alphaU: -viscosity_alpha_u " << alphaU;
  MOFEM_LOG("EP", Sev::inform) << "alphaW: -viscosity_alpha_w " << alphaW;
  MOFEM_LOG("EP", Sev::inform)
      << "alphaOmega: -viscosity_alpha_omega " << alphaOmega;
  MOFEM_LOG("EP", Sev::inform) << "alphaRho: -density_alpha_rho " << alphaRho;
  MOFEM_LOG("EP", Sev::inform)
      << "Rotations: -rotations " << list_rots[EshelbianCore::rotSelector];
  MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
                               << list_rots[EshelbianCore::gradApproximator];
  MOFEM_LOG("EP", Sev::inform)
      << "Symmetry: -symm " << list_symm[EshelbianCore::symmetrySelector];
  if (exponentBase != exp(1))
    MOFEM_LOG("EP", Sev::inform)
        << "Base exponent: -exponent_base " << EshelbianCore::exponentBase;
  else
    MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
  MOFEM_LOG("EP", Sev::inform)
      << "Stretch: -stretches " << list_stretches[choice_stretch];
  MOFEM_LOG("EP", Sev::inform) << "No stretch: -no_stretch " << (noStretch)
      ? "yes"
      : "no";
 
  MOFEM_LOG("EP", Sev::inform)
          << "Dynamic relaxation: -dynamic_relaxation " << (dynamicRelaxation)
      ? "yes"
      : "no";
  MOFEM_LOG("EP", Sev::inform)
          << "Singularity: -set_singularity " << (setSingularity)
      ? "yes"
      : "no";
  MOFEM_LOG("EP", Sev::inform)
          << "L2 user base scale: -l2_user_base_scale " << (l2UserBaseScale)
      ? "yes"
      : "no";
 
  MOFEM_LOG("EP", Sev::inform) << "Cracking on: -cracking_on " << (crackingOn)
      ? "yes"
      : "no";
  MOFEM_LOG("EP", Sev::inform)
      << "Cracking start time: -cracking_start_time " << crackingStartTime;
  MOFEM_LOG("EP", Sev::inform)
      << "Griffith energy: -griffith_energy " << griffithEnergy;
  MOFEM_LOG("EP", Sev::inform)
      << "Energy release variant: -energy_release_variant "
      << list_release[EshelbianCore::energyReleaseSelector];
  MOFEM_LOG("EP", Sev::inform)
      << "Number of J integral contours: -nb_J_integral_contours "
      << nbJIntegralContours;
 
#ifdef ENABLE_PYTHON_BINDING
  auto file_exists = [](std::string myfile) {
    std::ifstream file(myfile.c_str());
    if (file) {
      return true;
    }
    return false;
  };
 
  if (file_exists(analytical_expr_file_name)) {
    MOFEM_LOG("EP", Sev::inform) << analytical_expr_file_name << " file found";
 
    AnalyticalExprPythonPtr = boost::make_shared<AnalyticalExprPython>();
    CHKERR AnalyticalExprPythonPtr->analyticalExprInit(
        analytical_expr_file_name);
    AnalyticalExprPythonWeakPtr = AnalyticalExprPythonPtr;
  } else {
    MOFEM_LOG("EP", Sev::warning)
        << analytical_expr_file_name << " file NOT found";
  }
#endif
 
  if (spaceH1Order == -1)
    spaceH1Order = spaceOrder;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::addFields(const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  auto get_tets = [&]() {
    Range tets;
    CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
    return tets;
  };
 
  auto get_tets_skin = [&]() {
    Range tets_skin_part;
    Skinner skin(&mField.get_moab());
    CHKERR skin.find_skin(0, get_tets(), false, tets_skin_part);
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    Range tets_skin;
    CHKERR pcomm->filter_pstatus(tets_skin_part,
                                 PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &tets_skin);
    return tets_skin;
  };
 
  auto subtract_boundary_conditions = [&](auto &&tets_skin) {
    // That mean, that hybrid field on all faces on which traction is applied,
    // on other faces, or enforcing displacements as
    // natural boundary condition.
    if (bcSpatialTractionVecPtr)
      for (auto &v : *bcSpatialTractionVecPtr) {
        tets_skin = subtract(tets_skin, v.faces);
      }
    if (bcSpatialAnalyticalTractionVecPtr)
      for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
        tets_skin = subtract(tets_skin, v.faces);
      }      
 
    if (bcSpatialPressureVecPtr)
      for (auto &v : *bcSpatialPressureVecPtr) {
        tets_skin = subtract(tets_skin, v.faces);
      }
 
    return tets_skin;
  };
 
  auto add_blockset = [&](auto block_name, auto &&tets_skin) {
    auto crack_faces =
        get_range_from_block(mField, "block_name", SPACE_DIM - 1);
    tets_skin.merge(crack_faces);
    return tets_skin;
  };
 
  auto subtract_blockset = [&](auto block_name, auto &&tets_skin) {
    auto contact_range =
        get_range_from_block(mField, block_name, SPACE_DIM - 1);
    tets_skin = subtract(tets_skin, contact_range);
    return tets_skin;
  };
 
  auto get_stress_trace_faces = [&](auto &&tets_skin) {
    Range faces;
    CHKERR mField.get_moab().get_adjacencies(get_tets(), SPACE_DIM - 1, true,
                                             faces, moab::Interface::UNION);
    Range trace_faces = subtract(faces, tets_skin);
    return trace_faces;
  };
 
  auto tets = get_tets();
 
  // remove also contact faces, i.e. that is also kind of hybrid field but
  // named but used to enforce contact conditions
  auto trace_faces = get_stress_trace_faces(
 
      subtract_blockset("CONTACT",
                        subtract_boundary_conditions(get_tets_skin()))
 
  );
 
  contactFaces = boost::make_shared<Range>(intersect(
      trace_faces, get_range_from_block(mField, "CONTACT", SPACE_DIM - 1)));
  skeletonFaces =
      boost::make_shared<Range>(subtract(trace_faces, *contactFaces));
  // materialSkeletonFaces =
  //     boost::make_shared<Range>(get_stress_trace_faces(Range()));
 
#ifndef NDEBUG
  if (contactFaces->size())
    CHKERR save_range(mField.get_moab(),
                      "contact_faces_" +
                          std::to_string(mField.get_comm_rank()) + ".vtk",
                      *contactFaces);
  if (skeletonFaces->size())
    CHKERR save_range(mField.get_moab(),
                      "skeleton_faces_" +
                          std::to_string(mField.get_comm_rank()) + ".vtk",
                      *skeletonFaces);
    // if (skeletonFaces->size())
    //   CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
    //                     *materialSkeletonFaces);
#endif
 
  auto add_broken_hdiv_field = [this, meshset](const std::string field_name,
                                               const int order) {
    MoFEMFunctionBegin;
 
    constexpr FieldApproximationBase base = DEMKOWICZ_JACOBI_BASE;
 
    auto get_side_map_hdiv = [&]() {
      return std::vector<
 
          std::pair<EntityType,
                    std::function<MoFEMErrorCode(BaseFunction::DofsSideMap &)>
 
                    >>{
 
          {MBTET,
           [&](BaseFunction::DofsSideMap &dofs_side_map) -> MoFEMErrorCode {
             return TetPolynomialBase::setDofsSideMap(HDIV, DISCONTINUOUS, base,
                                                      dofs_side_map);
           }}
 
      };
    };
 
    CHKERR mField.add_broken_field(field_name, HDIV, base, SPACE_DIM,
                                   get_side_map_hdiv(), MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_l2_field = [this, meshset](const std::string field_name,
                                      const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_h1_field = [this, meshset](const std::string field_name,
                                      const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
    CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_l2_field_by_range = [this](const std::string field_name,
                                      const int order, const int dim,
                                      const int field_dim, Range &&r) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, field_dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(r);
    CHKERR mField.add_ents_to_field_by_dim(r, dim, field_name);
    CHKERR mField.set_field_order(r, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_bubble_field = [this, meshset](const std::string field_name,
                                          const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
                            MF_ZERO);
    // Modify field
    auto field_ptr = mField.get_field_structure(field_name);
    auto field_order_table =
        const_cast<Field *>(field_ptr)->getFieldOrderTable();
    auto get_cgg_bubble_order_zero = [](int p) { return 0; };
    auto get_cgg_bubble_order_tet = [](int p) {
      return NBVOLUMETET_CCG_BUBBLE(p);
    };
    field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
    field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
    field_order_table[MBTRI] = get_cgg_bubble_order_zero;
    field_order_table[MBTET] = get_cgg_bubble_order_tet;
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_user_l2_field = [this, meshset](const std::string field_name,
                                           const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, L2, USER_BASE, dim, MB_TAG_DENSE,
                            MF_ZERO);
    // Modify field
    auto field_ptr = mField.get_field_structure(field_name);
    auto field_order_table =
        const_cast<Field *>(field_ptr)->getFieldOrderTable();
    auto zero_dofs = [](int p) { return 0; };
    auto dof_l2_tet = [](int p) { return NBVOLUMETET_L2(p); };
    field_order_table[MBVERTEX] = zero_dofs;
    field_order_table[MBEDGE] = zero_dofs;
    field_order_table[MBTRI] = zero_dofs;
    field_order_table[MBTET] = dof_l2_tet;
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  // spatial fields
  CHKERR add_broken_hdiv_field(piolaStress, spaceOrder);
  CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
  CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
  CHKERR add_user_l2_field(rotAxis, spaceOrder - 1, 3);
  CHKERR add_user_l2_field(stretchTensor, noStretch ? -1 : spaceOrder, 6);
 
  if (!skeletonFaces)
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
  if (!contactFaces)
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No contact faces");
 
  auto get_hybridised_disp = [&]() {
    auto faces = *skeletonFaces;
    auto skin = subtract_boundary_conditions(get_tets_skin());
    for (auto &bc : *bcSpatialNormalDisplacementVecPtr) {
      faces.merge(intersect(bc.faces, skin));
    }
    return faces;
  };
 
  CHKERR add_l2_field_by_range(hybridSpatialDisp, spaceOrder - 1, 2, 3,
                               get_hybridised_disp());
  CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
                               Range(*contactFaces));
 
  // spatial displacement
  CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
  // material positions
  CHKERR add_h1_field(materialH1Positions, 2, 3);
 
  // Eshelby stress
  // CHKERR add_broken_hdiv_field(eshelbyStress, spaceOrder);
  // CHKERR add_l2_field(materialL2Disp, spaceOrder - 1, 3);
  // CHKERR add_l2_field_by_range(hybridMaterialDisp, spaceOrder - 1, 2, 3,
  //                              Range(*materialSkeletonFaces));
 
  CHKERR mField.build_fields();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::projectGeometry(const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  Range meshset_ents;
  CHKERR mField.get_moab().get_entities_by_handle(meshset, meshset_ents);
 
  auto project_ho_geometry = [&](auto field) {
    Projection10NodeCoordsOnField ent_method(mField, materialH1Positions);
    return mField.loop_dofs(field, ent_method);
  };
  CHKERR project_ho_geometry(materialH1Positions);
 
  auto get_adj_front_edges = [&](auto &front_edges) {
    Range front_crack_nodes;
    Range crack_front_edges_with_both_nodes_not_at_front;
 
    if (mField.get_comm_rank() == 0) {
      auto &moab = mField.get_moab();
      CHK_MOAB_THROW(
          moab.get_connectivity(front_edges, front_crack_nodes, true),
          "get_connectivity failed");
      Range crack_front_edges;
      CHK_MOAB_THROW(moab.get_adjacencies(front_crack_nodes, SPACE_DIM - 2,
                                          false, crack_front_edges,
                                          moab::Interface::UNION),
                     "get_adjacencies failed");
      Range crack_front_edges_nodes;
      CHK_MOAB_THROW(moab.get_connectivity(crack_front_edges,
                                           crack_front_edges_nodes, true),
                     "get_connectivity failed");
      // those nodes are hannging nodes
      crack_front_edges_nodes =
          subtract(crack_front_edges_nodes, front_crack_nodes);
      Range crack_front_edges_with_both_nodes_not_at_front;
      CHK_MOAB_THROW(
          moab.get_adjacencies(crack_front_edges_nodes, 1, false,
                               crack_front_edges_with_both_nodes_not_at_front,
                               moab::Interface::UNION),
          "get_adjacencies failed");
      // those edges are have one node not at the crack front
      crack_front_edges_with_both_nodes_not_at_front =
          intersect(crack_front_edges,
                    crack_front_edges_with_both_nodes_not_at_front);
    }
 
    front_crack_nodes = send_type(mField, front_crack_nodes, MBVERTEX);
    crack_front_edges_with_both_nodes_not_at_front = send_type(
        mField, crack_front_edges_with_both_nodes_not_at_front, MBEDGE);
 
    return std::make_pair(boost::make_shared<Range>(front_crack_nodes),
                          boost::make_shared<Range>(
                              crack_front_edges_with_both_nodes_not_at_front));
  };
 
  crackFaces = boost::make_shared<Range>(
      get_range_from_block(mField, "CRACK", SPACE_DIM - 1));
  frontEdges =
      boost::make_shared<Range>(get_crack_front_edges(mField, *crackFaces));
  auto [front_vertices, front_adj_edges] = get_adj_front_edges(*frontEdges);
  frontVertices = front_vertices;
  frontAdjEdges = front_adj_edges;
 
  MOFEM_LOG("EP", Sev::inform)
      << "Number of crack faces: " << crackFaces->size();
  MOFEM_LOG("EP", Sev::inform)
      << "Number of front edges: " << frontEdges->size();
  MOFEM_LOG("EP", Sev::inform)
      << "Number of front vertices: " << frontVertices->size();
  MOFEM_LOG("EP", Sev::inform)
      << "Number of front adjacent edges: " << frontAdjEdges->size();
 
#ifndef NDEBUG
  if (crackingOn) {
    auto rank = mField.get_comm_rank();
    // CHKERR save_range(mField.get_moab(),
    //                   (boost::format("meshset_ents_%d.vtk") % rank).str(),
    //                   meshset_ents);
    CHKERR save_range(mField.get_moab(),
                      (boost::format("crack_faces_%d.vtk") % rank).str(),
                      *crackFaces);
    CHKERR save_range(mField.get_moab(),
                      (boost::format("front_edges_%d.vtk") % rank).str(),
                      *frontEdges);
    // CHKERR save_range(mField.get_moab(),
    //                   (boost::format("front_vertices_%d.vtk") % rank).str(),
    //                   *frontVertices);
    // CHKERR save_range(mField.get_moab(),
    //                   (boost::format("front_adj_edges_%d.vtk") % rank).str(),
    //                   *frontAdjEdges);
  }
#endif // NDEBUG
 
  auto set_singular_dofs = [&](auto &front_adj_edges, auto &front_vertices) {
    MoFEMFunctionBegin;
    auto &moab = mField.get_moab();
 
    double eps = 1;
    double beta = 0;
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-singularity_eps", &beta,
                                 PETSC_NULLPTR);
    MOFEM_LOG("EP", Sev::inform) << "Singularity eps " << beta;
    eps -= beta;
 
    auto field_blas = mField.getInterface<FieldBlas>();
    auto lambda =
        [&](boost::shared_ptr<FieldEntity> field_entity_ptr) -> MoFEMErrorCode {
      MoFEMFunctionBegin;
      FTENSOR_INDEX(3, i);
      FTENSOR_INDEX(3, j);
 
      auto nb_dofs = field_entity_ptr->getEntFieldData().size();
      if (nb_dofs == 0) {
        MoFEMFunctionReturnHot(0);
      }
 
#ifndef NDEBUG
      if (field_entity_ptr->getNbOfCoeffs() != 3)
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "Expected 3 coefficients per edge");
      if (nb_dofs % 3 != 0)
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "Expected multiple of 3 coefficients per edge");
#endif // NDEBUG
 
      auto get_conn = [&]() {
        int num_nodes;
        const EntityHandle *conn;
        CHKERR moab.get_connectivity(field_entity_ptr->getEnt(), conn,
                                     num_nodes, false);
        return std::make_pair(conn, num_nodes);
      };
 
      auto get_dir = [&](auto &&conn_p) {
        auto [conn, num_nodes] = conn_p;
        double coords[6];
        CHKERR moab.get_coords(conn, num_nodes, coords);
        FTensor::Tensor1<double, 3> t_edge_dir{coords[3] - coords[0],
                                               coords[4] - coords[1],
                                               coords[5] - coords[2]};
        return t_edge_dir;
      };
 
      auto get_singularity_dof = [&](auto &&conn_p, auto &&t_edge_dir) {
        auto [conn, num_nodes] = conn_p;
        FTensor::Tensor1<double, 3> t_singularity_dof{0., 0., 0.};
        if (front_vertices.find(conn[0]) != front_vertices.end()) {
          t_singularity_dof(i) = t_edge_dir(i) * (-eps);
        } else if (front_vertices.find(conn[1]) != front_vertices.end()) {
          t_singularity_dof(i) = t_edge_dir(i) * eps;
        }
        return t_singularity_dof;
      };
 
      auto t_singularity_dof =
          get_singularity_dof(get_conn(), get_dir(get_conn()));
 
      auto field_data = field_entity_ptr->getEntFieldData();
      FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_dof{
          &field_data[0], &field_data[1], &field_data[2]};
 
      t_dof(i) =  t_singularity_dof(i);
      ++t_dof;
      for (auto n = 1; n < field_data.size() / 3; ++n) {    
        t_dof(i) = 0;
        ++t_dof;
      }
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR field_blas->fieldLambdaOnEntities(lambda, materialH1Positions,
                                             &front_adj_edges);
 
    MoFEMFunctionReturn(0);
  };
 
  if (setSingularity)
    CHKERR set_singular_dofs(*frontAdjEdges, *frontVertices);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EshelbianCore::addVolumeFiniteElement(const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  // set finite element fields
  auto add_field_to_fe = [this](const std::string fe,
                                const std::string field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
    MoFEMFunctionReturn(0);
  };
 
  if (!mField.check_finite_element(elementVolumeName)) {
    CHKERR mField.add_finite_element(elementVolumeName, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(meshset, MBTET,
                                                     elementVolumeName);
 
    CHKERR add_field_to_fe(elementVolumeName, piolaStress);
    CHKERR add_field_to_fe(elementVolumeName, bubbleField);
    if (!noStretch)
      CHKERR add_field_to_fe(elementVolumeName, stretchTensor);
    CHKERR add_field_to_fe(elementVolumeName, rotAxis);
    CHKERR add_field_to_fe(elementVolumeName, spatialL2Disp);
    CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
    CHKERR add_field_to_fe(elementVolumeName, contactDisp);
    CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
                                                       materialH1Positions);
 
    // build finite elements data structures
    CHKERR mField.build_finite_elements(elementVolumeName);
  }
 
  // if (!mField.check_finite_element(materialVolumeElement)) {
 
  //   Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM - 2);
 
  //   Range front_elements;
  //   for (auto l = 0; l != frontLayers; ++l) {
  //     Range front_elements_layer;
  //     CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM, true,
  //                                              front_elements_layer,
  //                                              moab::Interface::UNION);
  //     front_elements.merge(front_elements_layer);
  //     front_edges.clear();
  //     CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
  //                                              SPACE_DIM - 2, true,
  //                                              front_edges,
  //                                              moab::Interface::UNION);
  //   }
 
  //   CHKERR mField.add_finite_element(materialVolumeElement, MF_ZERO);
  //   CHKERR mField.add_ents_to_finite_element_by_type(front_elements, MBTET,
  //                                                    materialVolumeElement);
  //   // CHKERR add_field_to_fe(materialVolumeElement, eshelbyStress);
  //   // CHKERR add_field_to_fe(materialVolumeElement, materialL2Disp);
  //   CHKERR mField.build_finite_elements(materialVolumeElement);
  // }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EshelbianCore::addBoundaryFiniteElement(const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  Range meshset_ents;
  CHKERR mField.get_moab().get_entities_by_handle(meshset, meshset_ents);
 
  auto set_fe_adjacency = [&](auto fe_name) {
    MoFEMFunctionBegin;
    parentAdjSkeletonFunctionDim2 =
        boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
            bitAdjParent, bitAdjParentMask, bitAdjEnt, bitAdjEntMask);
    CHKERR mField.modify_finite_element_adjacency_table(
        fe_name, MBTRI, *parentAdjSkeletonFunctionDim2);
    MoFEMFunctionReturn(0);
  };
 
  // set finite element fields
  auto add_field_to_fe = [this](const std::string fe,
                                const std::string field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, materialH1Positions);
    MoFEMFunctionReturn(0);
  };
 
  if (!mField.check_finite_element(naturalBcElement)) {
 
    Range natural_bc_elements;
    if (bcSpatialDispVecPtr) {
      for (auto &v : *bcSpatialDispVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialRotationVecPtr) {
      for (auto &v : *bcSpatialRotationVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialNormalDisplacementVecPtr) {
      for (auto &v : *bcSpatialNormalDisplacementVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialAnalyticalDisplacementVecPtr) {
      for (auto &v : *bcSpatialAnalyticalDisplacementVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialTractionVecPtr) {
      for (auto &v : *bcSpatialTractionVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialAnalyticalTractionVecPtr) {
      for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialPressureVecPtr) {
      for (auto &v : *bcSpatialPressureVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    natural_bc_elements = intersect(natural_bc_elements, meshset_ents);
 
    CHKERR mField.add_finite_element(naturalBcElement, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
                                                     naturalBcElement);
    CHKERR add_field_to_fe(naturalBcElement, piolaStress);
    CHKERR add_field_to_fe(naturalBcElement, hybridSpatialDisp);
    CHKERR set_fe_adjacency(naturalBcElement);
    CHKERR mField.build_finite_elements(naturalBcElement);
  }
 
  auto get_skin = [&](auto &body_ents) {
    Skinner skin(&mField.get_moab());
    Range skin_ents;
    CHKERR skin.find_skin(0, body_ents, false, skin_ents);
    return skin_ents;
  };
 
  auto filter_true_skin = [&](auto &&skin) {
    Range boundary_ents;
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &boundary_ents);
    return boundary_ents;
  };
 
  if (!mField.check_finite_element(skinElement)) {
 
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM,
                                                       body_ents);
    auto skin = filter_true_skin(get_skin(body_ents));
 
    CHKERR mField.add_finite_element(skinElement, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(skin, MBTRI, skinElement);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       spatialH1Disp);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       materialH1Positions);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       contactDisp);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       hybridSpatialDisp);
    // CHKERR add_field_to_fe(skinElement, hybridSpatialDisp);
    // CHKERR add_field_to_fe(skinElement, hybridMaterialDisp);
 
    CHKERR mField.build_finite_elements(skinElement);
  }
 
  if (!mField.check_finite_element(contactElement)) {
    if (contactFaces) {
      MOFEM_LOG("EP", Sev::inform)
          << "Contact elements " << contactFaces->size();
      CHKERR mField.add_finite_element(contactElement, MF_ZERO);
      CHKERR mField.add_ents_to_finite_element_by_type(*contactFaces, MBTRI,
                                                       contactElement);
      CHKERR add_field_to_fe(contactElement, piolaStress);
      CHKERR add_field_to_fe(contactElement, hybridSpatialDisp);
      CHKERR add_field_to_fe(contactElement, contactDisp);
      CHKERR add_field_to_fe(contactElement, spatialH1Disp);
      CHKERR set_fe_adjacency(contactElement);
      CHKERR mField.build_finite_elements(contactElement);
    }
  }
 
  if (!mField.check_finite_element(skeletonElement)) {
    if (!skeletonFaces)
      SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
    MOFEM_LOG("EP", Sev::inform)
        << "Skeleton elements " << skeletonFaces->size();
    CHKERR mField.add_finite_element(skeletonElement, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(*skeletonFaces, MBTRI,
                                                     skeletonElement);
    CHKERR add_field_to_fe(skeletonElement, piolaStress);
    CHKERR add_field_to_fe(skeletonElement, hybridSpatialDisp);
    CHKERR add_field_to_fe(skeletonElement, contactDisp);
    CHKERR add_field_to_fe(skeletonElement, spatialH1Disp);
    CHKERR set_fe_adjacency(skeletonElement);
    CHKERR mField.build_finite_elements(skeletonElement);
  }
 
  //   if (!mField.check_finite_element(materialSkeletonElement)) {
 
  //     Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM -
  //     2);
 
  //     Range front_elements;
  //     for (auto l = 0; l != frontLayers; ++l) {
  //       Range front_elements_layer;
  //       CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM,
  //       true,
  //                                                front_elements_layer,
  //                                                moab::Interface::UNION);
  //       front_elements.merge(front_elements_layer);
  //       front_edges.clear();
  //       CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
  //                                                SPACE_DIM - 2, true,
  //                                                front_edges,
  //                                                moab::Interface::UNION);
  //     }
  //     Range body_ents;
  //     CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
  //     body_ents); Range front_elements_faces; CHKERR
  //     mField.get_moab().get_adjacencies(front_elements, SPACE_DIM - 1,
  //                                              true, front_elements_faces,
  //                                              moab::Interface::UNION);
  //     auto body_skin = filter_true_skin(get_skin(body_ents));
  //     auto front_elements_skin = filter_true_skin(get_skin(front_elements));
  //     Range material_skeleton_faces =
  //         subtract(front_elements_faces, front_elements_skin);
  //     material_skeleton_faces.merge(intersect(front_elements_skin,
  //     body_skin));
 
  // #ifndef NDEBUG
  //     CHKERR save_range(mField.get_moab(), "front_elements.vtk",
  //     front_elements); CHKERR save_range(mField.get_moab(),
  //     "front_elements_skin.vtk",
  //                       front_elements_skin);
  //     CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
  //                       material_skeleton_faces);
  // #endif
 
  //     CHKERR mField.add_finite_element(materialSkeletonElement, MF_ZERO);
  //     CHKERR mField.add_ents_to_finite_element_by_type(
  //         material_skeleton_faces, MBTRI, materialSkeletonElement);
  //     // CHKERR add_field_to_fe(materialSkeletonElement, eshelbyStress);
  //     // CHKERR add_field_to_fe(materialSkeletonElement, hybridMaterialDisp);
  //     CHKERR set_fe_adjacency(materialSkeletonElement);
  //     CHKERR mField.build_finite_elements(materialSkeletonElement);
  //   }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::addDMs(const BitRefLevel bit,
                                     const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  // find adjacencies between finite elements and dofs
  CHKERR mField.build_adjacencies(bit, QUIET);
 
  // Create coupled problem
  dM = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
                            BitRefLevel().set());
  CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
  CHKERR DMMoFEMAddElement(dM, elementVolumeName);
  CHKERR DMMoFEMAddElement(dM, naturalBcElement);
  CHKERR DMMoFEMAddElement(dM, skeletonElement);
  CHKERR DMMoFEMAddElement(dM, contactElement);
  CHKERR DMMoFEMAddElement(dM, skinElement);
 
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
  CHKERR DMSetUp(dM);
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
 
  auto remove_dofs_on_broken_skin = [&](const std::string prb_name) {
    MoFEMFunctionBegin;
    for (int d : {0, 1, 2}) {
      std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
      CHKERR mField.getInterface<ProblemsManager>()
          ->getSideDofsOnBrokenSpaceEntities(
              dofs_to_remove, prb_name, ROW, piolaStress,
              (*bcSpatialFreeTractionVecPtr)[d], SPACE_DIM, d, d);
      // remove piola dofs, i.e. traction free boundary
      CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, ROW,
                                                                dofs_to_remove);
      CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, COL,
                                                                dofs_to_remove);
    }
    MoFEMFunctionReturn(0);
  };
  CHKERR remove_dofs_on_broken_skin("ESHELBY_PLASTICITY");
 
  // Create elastic sub-problem
  dmElastic = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
  CHKERR DMMoFEMSetDestroyProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, piolaStress);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis);
  if (!noStretch) {
    CHKERR DMMoFEMAddSubFieldRow(dmElastic, stretchTensor);
  }
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, hybridSpatialDisp);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, contactDisp);
  CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName);
  CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement);
  CHKERR DMMoFEMAddElement(dmElastic, skeletonElement);
  CHKERR DMMoFEMAddElement(dmElastic, contactElement);
  CHKERR DMMoFEMAddElement(dmElastic, skinElement);
  CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
  CHKERR DMSetUp(dmElastic);
 
  // dmMaterial = createDM(mField.get_comm(), "DMMOFEM");
  // CHKERR DMMoFEMCreateSubDM(dmMaterial, dM, "MATERIAL_PROBLEM");
  // CHKERR DMMoFEMSetDestroyProblem(dmMaterial, PETSC_TRUE);
  // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, eshelbyStress);
  // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, materialL2Disp);
  // CHKERR DMMoFEMAddElement(dmMaterh elementVolumeName);
  // CHKERR DMMoFEMAddElement(dmMaterial, naturalBcElement);
  // CHKERR DMMoFEMAddElement(dmMaterial, skinElement);
  // CHKERR DMMoFEMSetSquareProblem(dmMaterial, PETSC_TRUE);
  // CHKERR DMMoFEMSetIsPartitioned(dmMaterial, PETSC_TRUE);
  // CHKERR DMSetUp(dmMaterial);
 
  auto set_zero_block = [&]() {
    MoFEMFunctionBegin;
    if (!noStretch) {
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", spatialL2Disp, stretchTensor);
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", stretchTensor, spatialL2Disp);
    }
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", spatialL2Disp, rotAxis);
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", rotAxis, spatialL2Disp);
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", spatialL2Disp, bubbleField);
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", bubbleField, spatialL2Disp);
    if (!noStretch) {
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", bubbleField, bubbleField);
      CHKERR
      mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", piolaStress, piolaStress);
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", bubbleField, piolaStress);
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", piolaStress, bubbleField);
    }
 
    solTSStep = createDMVector(dmElastic);
    MoFEMFunctionReturn(0);
  };
 
  auto set_section = [&]() {
    MoFEMFunctionBegin;
    PetscSection section;
    CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
                                                           &section);
    CHKERR DMSetSection(dmElastic, section);
    CHKERR DMSetGlobalSection(dmElastic, section);
    CHKERR PetscSectionDestroy(&section);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_zero_block();
  CHKERR set_section();
 
  dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
  CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, spatialH1Disp);
  CHKERR DMMoFEMAddElement(dmPrjSpatial, elementVolumeName);
  CHKERR DMMoFEMSetSquareProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMSetUp(dmPrjSpatial);
 
  // CHKERR mField.getInterface<BcManager>()
  //     ->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
  //         "PROJECT_SPATIAL", spatialH1Disp, true, false);
 
  MoFEMFunctionReturn(0);
}
 
BcDisp::BcDisp(std::string name, std::vector<double> attr, Range faces)
    : blockName(name), faces(faces) {
  vals.resize(3, false);
  flags.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    vals[ii] = attr[ii];
    flags[ii] = static_cast<int>(attr[ii + 3]);
  }
 
  MOFEM_LOG("EP", Sev::inform) << "Add BCDisp " << name;
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCDisp vals " << vals[0] << " " << vals[1] << " " << vals[2];
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCDisp flags " << flags[0] << " " << flags[1] << " " << flags[2];
  MOFEM_LOG("EP", Sev::inform) << "Add BCDisp nb. of faces " << faces.size();
}
 
BcRot::BcRot(std::string name, std::vector<double> attr, Range faces)
    : blockName(name), faces(faces) {
  vals.resize(attr.size(), false);
  for (int ii = 0; ii != attr.size(); ++ii) {
    vals[ii] = attr[ii];
  }
  theta = attr[3];
}
 
TractionBc::TractionBc(std::string name, std::vector<double> attr, Range faces)
    : blockName(name), faces(faces) {
  vals.resize(3, false);
  flags.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    vals[ii] = attr[ii];
    flags[ii] = static_cast<int>(attr[ii + 3]);
  }
 
  MOFEM_LOG("EP", Sev::inform) << "Add BCForce " << name;
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCForce vals " << vals[0] << " " << vals[1] << " " << vals[2];
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCForce flags " << flags[0] << " " << flags[1] << " " << flags[2];
  MOFEM_LOG("EP", Sev::inform) << "Add BCForce nb. of faces " << faces.size();
}
 
NormalDisplacementBc::NormalDisplacementBc(std::string name,
                                           std::vector<double> attr,
                                           Range faces)
    : blockName(name), faces(faces) {
 
  blockName = name;
  if (attr.size() < 1) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong size of normal displacement BC");
  }
 
  val = attr[0];
 
  MOFEM_LOG("EP", Sev::inform) << "Add NormalDisplacementBc " << name;
  MOFEM_LOG("EP", Sev::inform) << "Add NormalDisplacementBc val " << val;
  MOFEM_LOG("EP", Sev::inform)
      << "Add NormalDisplacementBc nb. of faces " << faces.size();
}
 
PressureBc::PressureBc(std::string name, std::vector<double> attr, Range faces)
    : blockName(name), faces(faces) {
 
  blockName = name;
  if (attr.size() < 1) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong size of normal displacement BC");
  }
 
  val = attr[0];
 
  MOFEM_LOG("EP", Sev::inform) << "Add PressureBc " << name;
  MOFEM_LOG("EP", Sev::inform) << "Add PressureBc val " << val;
  MOFEM_LOG("EP", Sev::inform)
      << "Add PressureBc nb. of faces " << faces.size();
}
 
 
ExternalStrain::ExternalStrain(std::string name, std::vector<double> attr,
                           Range ents)
    : blockName(name), ents(ents) {
 
  blockName = name;
  if (attr.size() < 2) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong size of external strain attribute");
  }
 
  val = attr[0];
  bulkModulusK = attr[1]; 
 
  MOFEM_LOG("EP", Sev::inform) << "Add ExternalStrain " << name;
  MOFEM_LOG("EP", Sev::inform) << "Add ExternalStrain val " << val;
  MOFEM_LOG("EP", Sev::inform) << "Add ExternalStrain bulk modulus K "
                               << bulkModulusK;
  MOFEM_LOG("EP", Sev::inform)
      << "Add ExternalStrain nb. of tets " << ents.size();
 
 
}
 
AnalyticalDisplacementBc::AnalyticalDisplacementBc(std::string name,
                                                   std::vector<double> attr,
                                                   Range faces)
    : blockName(name), faces(faces) {
 
  blockName = name;
  if (attr.size() < 3) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong size of analytical displacement BC");
  }
 
  flags.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    flags[ii] = attr[ii];
  }
 
  MOFEM_LOG("EP", Sev::inform) << "Add AnalyticalDisplacementBc " << name;
  MOFEM_LOG("EP", Sev::inform)
      << "Add AnalyticalDisplacementBc flags " << flags[0] << " " << flags[1]
      << " " << flags[2];
  MOFEM_LOG("EP", Sev::inform)
      << "Add AnalyticalDisplacementBc nb. of faces " << faces.size();
}
 
AnalyticalTractionBc::AnalyticalTractionBc(std::string name,
                                                   std::vector<double> attr,
                                                   Range faces)
    : blockName(name), faces(faces) {
 
  blockName = name;
  if (attr.size() < 3) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong size of analytical traction BC");
  }
 
  flags.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    flags[ii] = attr[ii];
  }
 
  MOFEM_LOG("EP", Sev::inform) << "Add AnalyticalTractionBc " << name;
  MOFEM_LOG("EP", Sev::inform)
      << "Add AnalyticalTractionBc flags " << flags[0] << " " << flags[1]
      << " " << flags[2];
  MOFEM_LOG("EP", Sev::inform)
      << "Add AnalyticalTractionBc nb. of faces " << faces.size();
}
 
 
MoFEMErrorCode
EshelbianCore::getTractionFreeBc(const EntityHandle meshset,
                                 boost::shared_ptr<TractionFreeBc> &bc_ptr,
                                 const std::string contact_set_name) {
  MoFEMFunctionBegin;
 
  // get skin from all tets
  Range tets;
  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
  Range tets_skin_part;
  Skinner skin(&mField.get_moab());
  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
  Range tets_skin;
  CHKERR pcomm->filter_pstatus(tets_skin_part,
                               PSTATUS_SHARED | PSTATUS_MULTISHARED,
                               PSTATUS_NOT, -1, &tets_skin);
 
  bc_ptr->resize(3);
  for (int dd = 0; dd != 3; ++dd)
    (*bc_ptr)[dd] = tets_skin;
 
  // Do not remove dofs on which traction is applied
  if (bcSpatialDispVecPtr)
    for (auto &v : *bcSpatialDispVecPtr) {
      if (v.flags[0])
        (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      if (v.flags[1])
        (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      if (v.flags[2])
        (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  // Do not remove dofs on which rotation is applied
  if (bcSpatialRotationVecPtr)
    for (auto &v : *bcSpatialRotationVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialNormalDisplacementVecPtr)
    for (auto &v : *bcSpatialNormalDisplacementVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialAnalyticalDisplacementVecPtr)
    for (auto &v : *bcSpatialAnalyticalDisplacementVecPtr) {
      if (v.flags[0])
        (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      if (v.flags[1])
        (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      if (v.flags[2])
        (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialTractionVecPtr)
    for (auto &v : *bcSpatialTractionVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialAnalyticalTractionVecPtr)
    for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
        
  if (bcSpatialPressureVecPtr)
    for (auto &v : *bcSpatialPressureVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  // remove contact
  for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
           std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
    Range faces;
    CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                              true);
    (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
    (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
    (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
  }
 
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Set integration rule on element
 * \param order on row
 * \param order on column
 * \param order on data
 *
 * Use maximal oder on data in order to determine integration rule
 *
 */
struct VolRule {
  int operator()(int p_row, int p_col, int p_data) const {
    return 2 * p_data + 1;
  }
};
 
struct FaceRule {
  int operator()(int p_row, int p_col, int p_data) const {
    return 2 * (p_data + 1);
  }
};
 
 
 
MoFEMErrorCode EshelbianCore::setBaseVolumeElementOps(
    const int tag, const bool do_rhs, const bool do_lhs, const bool calc_rates,
    SmartPetscObj<Vec> var_vec,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe) {
  MoFEMFunctionBegin;
 
  auto bubble_cache =
      boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0, MatrixDouble());
  fe->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      fe->getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions, frontAdjEdges);
 
  // set integration rule
  fe->getRuleHook = [](int, int, int) { return -1; };
  fe->setRuleHook = SetIntegrationAtFrontVolume(frontVertices, frontAdjEdges);
  // fe->getRuleHook = VolRule();
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
    dataAtPts->physicsPtr = physicalEquations;
  }
 
  // calculate fields values
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
      piolaStress, dataAtPts->getApproxPAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
      bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
      piolaStress, dataAtPts->getDivPAtPts()));
 
  if (noStretch) {
    fe->getOpPtrVector().push_back(
        physicalEquations->returnOpCalculateStretchFromStress(
            dataAtPts, physicalEquations));
  } else {
    fe->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValues<3>(
            stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
  }
 
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
 
  // H1 displacements
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
      spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
 
  // velocities
  if (calc_rates) {
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
        spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
    if (noStretch) {
    } else {
      fe->getOpPtrVector().push_back(
          new OpCalculateTensor2SymmetricFieldValuesDot<3>(
              stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
      fe->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradientDot<6, 3>(
              stretchTensor, dataAtPts->getGradLogStretchDotTensorAtPts(),
              MBTET));
    }
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
        rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradientDot<3, 3>(
        rotAxis, dataAtPts->getRotAxisGradDotAtPts(), MBTET));
 
    // acceleration
    if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
      fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
          spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
    }
  }
 
  // variations
  if (var_vec.use_count()) {
    fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
        piolaStress, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
        var_vec));
    fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
        bubbleField, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
        var_vec, MBMAXTYPE));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        rotAxis, dataAtPts->getVarRotAxisPts(), var_vec, MBTET));
    fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
        piolaStress, dataAtPts->getDivVarPiolaPts(), var_vec));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        spatialL2Disp, dataAtPts->getVarWL2Pts(), var_vec, MBTET));
 
    if (noStretch) {
      fe->getOpPtrVector().push_back(
          physicalEquations->returnOpCalculateVarStretchFromStress(
              dataAtPts, physicalEquations));
    } else {
      fe->getOpPtrVector().push_back(
          new OpCalculateTensor2SymmetricFieldValues<3>(
              stretchTensor, dataAtPts->getVarLogStreachPts(), var_vec, MBTET));
    }
  }
 
  // calculate other derived quantities
  fe->getOpPtrVector().push_back(new OpCalculateRotationAndSpatialGradient(
      dataAtPts, ((do_rhs || do_lhs) && calc_rates) ? alphaOmega : 0.));
 
  // evaluate integration points
  if (noStretch) {
  } else {
    fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
        tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setVolumeElementOps(
    const int tag, const bool add_elastic, const bool add_material,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
  MoFEMFunctionBegin;
 
  /** Contact requires that body is marked */
  auto get_body_range = [this](auto name, int dim) {
    std::map<int, Range> map;
 
    for (auto m_ptr :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % name).str()
 
                 ))
 
    ) {
      Range ents;
      CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
                                                            dim, ents, true),
                     "by dim");
      map[m_ptr->getMeshsetId()] = ents;
    }
 
    return map;
  };
 
  auto rule_contact = [](int, int, int o) { return -1; };
  auto refine = Tools::refineTriangle(contactRefinementLevels);
 
  auto set_rule_contact = [refine](
 
                              ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                              int order_col, int order_data
 
                          ) {
    MoFEMFunctionBegin;
    auto rule = 2 * order_data;
    fe_raw_ptr->gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
    MoFEMFunctionReturn(0);
  };
 
  // Right hand side
  fe_rhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
  CHKERR setBaseVolumeElementOps(tag, true, false, true, SmartPetscObj<Vec>(),
                                 fe_rhs);
 
  // elastic
  if (add_elastic) {
 
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialEquilibrium(spatialL2Disp, dataAtPts, alphaW, alphaRho));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialRotation(rotAxis, dataAtPts, alphaOmega));
    if (noStretch) {
      // do nothing - no stretch approximation
    } else {
      if (!internalStressTagName.empty()) {
        switch (internalStressInterpOrder) {
        case 0:
          fe_rhs->getOpPtrVector().push_back(
              new OpGetInternalStress<0>(dataAtPts, internalStressTagName));
          break;
        case 1:
          fe_rhs->getOpPtrVector().push_back(
              new OpGetInternalStress<1>(dataAtPts, internalStressTagName));
          break;
        default:
          SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
                   "Unsupported internal stress interpolation order %d",
                   internalStressInterpOrder);
        }
        if (internalStressVoigt) {
          fe_rhs->getOpPtrVector().push_back(
              new OpSpatialPhysicalInternalStress<true>(stretchTensor,
                                                        dataAtPts));
        } else {
          fe_rhs->getOpPtrVector().push_back(
              new OpSpatialPhysicalInternalStress<false>(stretchTensor,
                                                         dataAtPts));
        }
      }
      if (auto op = physicalEquations->returnOpSpatialPhysicalExternalStrain(
              stretchTensor, dataAtPts, externalStrainVecPtr, timeScaleMap)) {
        fe_rhs->getOpPtrVector().push_back(op);
      } else if (externalStrainVecPtr && !externalStrainVecPtr->empty()) {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
                "OpSpatialPhysicalExternalStrain not implemented for this "
                "material");
      }
 
      fe_rhs->getOpPtrVector().push_back(
          physicalEquations->returnOpSpatialPhysical(stretchTensor, dataAtPts,
                                                     alphaU));
    }
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyP(piolaStress, dataAtPts));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyBubble(bubbleField, dataAtPts));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyDivTerm(piolaStress, dataAtPts));
 
    auto set_hybridisation = [&](auto &pip) {
      MoFEMFunctionBegin;
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using EleOnSide =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
      using SideEleOp = EleOnSide::UserDataOperator;
      using BdyEleOp = BoundaryEle::UserDataOperator;
 
      // First: Iterate over skeleton FEs adjacent to Domain FEs
      // Note:  BoundaryEle, i.e. uses skeleton interation rule
      auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
          mField, skeletonElement, SPACE_DIM - 1, Sev::noisy);
      // op_loop_skeleton_side->getSideFEPtr()->getRuleHook = FaceRule();
      op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
        return -1;
      };
      op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
          SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
 
      CHKERR EshelbianPlasticity::
          AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
              op_loop_skeleton_side->getOpPtrVector(), {L2},
              materialH1Positions, frontAdjEdges);
 
      // Second: Iterate over domain FEs adjacent to skelton, particularly one
      // domain element.
      auto broken_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
      // Note: EleOnSide, i.e. uses on domain projected skeleton rule
      auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
          mField, elementVolumeName, SPACE_DIM, Sev::noisy);
      op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
          boost::make_shared<CGGUserPolynomialBase>();
      CHKERR
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
          materialH1Positions, frontAdjEdges);
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
      auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(piolaStress,
                                                               flux_mat_ptr));
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpSetFlux<SideEleOp>(broken_data_ptr, flux_mat_ptr));
 
      // Assemble on skeleton
      op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
      using OpC_dHybrid = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
          GAUSS>::OpBrokenSpaceConstrainDHybrid<SPACE_DIM>;
      using OpC_dBroken = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
          GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
      op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dHybrid(
          hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0)));
      auto hybrid_ptr = boost::make_shared<MatrixDouble>();
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(hybridSpatialDisp,
                                                      hybrid_ptr));
      op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dBroken(
          broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0)));
 
      // Add skeleton to domain pipeline
      pip.push_back(op_loop_skeleton_side);
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_contact = [&](auto &pip) {
      MoFEMFunctionBegin;
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using EleOnSide =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
      using SideEleOp = EleOnSide::UserDataOperator;
      using BdyEleOp = BoundaryEle::UserDataOperator;
 
      // First: Iterate over skeleton FEs adjacent to Domain FEs
      // Note:  BoundaryEle, i.e. uses skeleton interation rule
      auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
          mField, contactElement, SPACE_DIM - 1, Sev::noisy);
 
      op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
      op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
      CHKERR EshelbianPlasticity::
          AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
              op_loop_skeleton_side->getOpPtrVector(), {L2},
              materialH1Positions, frontAdjEdges);
 
      // Second: Iterate over domain FEs adjacent to skelton, particularly
      // one domain element.
      auto broken_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
 
      // Data storing contact fields
      auto contact_common_data_ptr =
          boost::make_shared<ContactOps::CommonData>();
 
      auto add_ops_domain_side = [&](auto &pip) {
        MoFEMFunctionBegin;
        // Note: EleOnSide, i.e. uses on domain projected skeleton rule
        auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
            mField, elementVolumeName, SPACE_DIM, Sev::noisy);
        op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
            boost::make_shared<CGGUserPolynomialBase>();
        CHKERR
        EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
            op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
            materialH1Positions, frontAdjEdges);
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpGetBrokenBaseSideData<SideEleOp>(piolaStress,
                                                   broken_data_ptr));
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
                piolaStress, contact_common_data_ptr->contactTractionPtr()));
        pip.push_back(op_loop_domain_side);
        MoFEMFunctionReturn(0);
      };
 
      auto add_ops_contact_rhs = [&](auto &pip) {
        MoFEMFunctionBegin;
        // get body id and SDF range
        auto contact_sfd_map_range_ptr =
            boost::make_shared<std::map<int, Range>>(
                get_body_range("CONTACT_SDF", SPACE_DIM - 1));
 
        pip.push_back(new OpCalculateVectorFieldValues<3>(
            contactDisp, contact_common_data_ptr->contactDispPtr()));
        auto u_h1_ptr = boost::make_shared<MatrixDouble>();
        pip.push_back(
            new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
        pip.push_back(new OpTreeSearch(
            contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
            get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
            nullptr));
        pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Rhs(
            contactDisp, contact_common_data_ptr, contactTreeRhs,
            contact_sfd_map_range_ptr));
        pip.push_back(
            new EshelbianPlasticity::OpConstrainBoundaryHDivRhs<A, GAUSS>(
                broken_data_ptr, contact_common_data_ptr, contactTreeRhs));
 
        MoFEMFunctionReturn(0);
      };
 
      // push ops to face/side pipeline
      CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
      CHKERR add_ops_contact_rhs(op_loop_skeleton_side->getOpPtrVector());
 
      // Add skeleton to domain pipeline
      pip.push_back(op_loop_skeleton_side);
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_hybridisation(fe_rhs->getOpPtrVector());
    CHKERR set_contact(fe_rhs->getOpPtrVector());
 
    // Body forces
    using BodyNaturalBC =
        NaturalBC<VolumeElementForcesAndSourcesCore::UserDataOperator>::
            Assembly<PETSC>::LinearForm<GAUSS>;
    using OpBodyForce =
        BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
 
    auto body_time_scale =
        boost::make_shared<DynamicRelaxationTimeScale>("body_force.txt");
    CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
        fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {body_time_scale},
        "BODY_FORCE", Sev::inform);
  }
 
  // Left hand side
  fe_lhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
  CHKERR setBaseVolumeElementOps(tag, true, true, true, SmartPetscObj<Vec>(),
                                 fe_lhs);
 
  // elastic
  if (add_elastic) {
 
    if (noStretch) {
      fe_lhs->getOpPtrVector().push_back(
          new OpSpatialConsistency_dP_dP(piolaStress, piolaStress, dataAtPts));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dP(
          bubbleField, piolaStress, dataAtPts));
      fe_lhs->getOpPtrVector().push_back(
          new OpSpatialConsistency_dBubble_dBubble(bubbleField, bubbleField,
                                                   dataAtPts));
    } else {
      fe_lhs->getOpPtrVector().push_back(
          physicalEquations->returnOpSpatialPhysical_du_du(
              stretchTensor, stretchTensor, dataAtPts, alphaU));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
          stretchTensor, piolaStress, dataAtPts, true));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
          stretchTensor, bubbleField, dataAtPts, true));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
          stretchTensor, rotAxis, dataAtPts,
          symmetrySelector == SYMMETRIC ? true : false));
    }
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
        spatialL2Disp, piolaStress, dataAtPts, true));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
        spatialL2Disp, spatialL2Disp, dataAtPts, alphaW, alphaRho));
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_domega(
        piolaStress, rotAxis, dataAtPts,
        symmetrySelector == SYMMETRIC ? true : false));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
        bubbleField, rotAxis, dataAtPts,
        symmetrySelector == SYMMETRIC ? true : false));
 
    if (symmetrySelector > SYMMETRIC) {
      if (!noStretch) {
        fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_du(
            rotAxis, stretchTensor, dataAtPts, false));
      }
      fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
          rotAxis, piolaStress, dataAtPts, false));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
          rotAxis, bubbleField, dataAtPts, false));
    }
    fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_domega(
        rotAxis, rotAxis, dataAtPts, alphaOmega));
    
    auto set_hybridisation = [&](auto &pip) {
      MoFEMFunctionBegin;
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using EleOnSide =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
      using SideEleOp = EleOnSide::UserDataOperator;
      using BdyEleOp = BoundaryEle::UserDataOperator;
 
      // First: Iterate over skeleton FEs adjacent to Domain FEs
      // Note:  BoundaryEle, i.e. uses skeleton interation rule
      auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
          mField, skeletonElement, SPACE_DIM - 1, Sev::noisy);
      // op_loop_skeleton_side->getSideFEPtr()->getRuleHook = FaceRule();
      op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
        return -1;
      };
      op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
          SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
      CHKERR EshelbianPlasticity::
          AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
              op_loop_skeleton_side->getOpPtrVector(), {L2},
              materialH1Positions, frontAdjEdges);
 
      // Second: Iterate over domain FEs adjacent to skelton, particularly one
      // domain element.
      auto broken_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
      // Note: EleOnSide, i.e. uses on domain projected skeleton rule
      auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
          mField, elementVolumeName, SPACE_DIM, Sev::noisy);
      op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
          boost::make_shared<CGGUserPolynomialBase>();
      CHKERR
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
          materialH1Positions, frontAdjEdges);
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
 
      op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
      using OpC = FormsIntegrators<BdyEleOp>::Assembly<A>::BiLinearForm<
          GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpC(hybridSpatialDisp, broken_data_ptr,
                  boost::make_shared<double>(1.0), true, false));
 
      pip.push_back(op_loop_skeleton_side);
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_contact = [&](auto &pip) {
      MoFEMFunctionBegin;
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using EleOnSide =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
      using SideEleOp = EleOnSide::UserDataOperator;
      using BdyEleOp = BoundaryEle::UserDataOperator;
 
      // First: Iterate over skeleton FEs adjacent to Domain FEs
      // Note:  BoundaryEle, i.e. uses skeleton interation rule
      auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
          mField, contactElement, SPACE_DIM - 1, Sev::noisy);
 
      op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
      op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
      CHKERR EshelbianPlasticity::
          AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
              op_loop_skeleton_side->getOpPtrVector(), {L2},
              materialH1Positions, frontAdjEdges);
 
      // Second: Iterate over domain FEs adjacent to skelton, particularly
      // one domain element.
      auto broken_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
 
      // Data storing contact fields
      auto contact_common_data_ptr =
          boost::make_shared<ContactOps::CommonData>();
 
      auto add_ops_domain_side = [&](auto &pip) {
        MoFEMFunctionBegin;
        // Note: EleOnSide, i.e. uses on domain projected skeleton rule
        auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
            mField, elementVolumeName, SPACE_DIM, Sev::noisy);
        op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
            boost::make_shared<CGGUserPolynomialBase>();
        CHKERR
        EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
            op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
            materialH1Positions, frontAdjEdges);
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpGetBrokenBaseSideData<SideEleOp>(piolaStress,
                                                   broken_data_ptr));
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
                piolaStress, contact_common_data_ptr->contactTractionPtr()));
        pip.push_back(op_loop_domain_side);
        MoFEMFunctionReturn(0);
      };
 
      auto add_ops_contact_lhs = [&](auto &pip) {
        MoFEMFunctionBegin;
        pip.push_back(new OpCalculateVectorFieldValues<3>(
            contactDisp, contact_common_data_ptr->contactDispPtr()));
        auto u_h1_ptr = boost::make_shared<MatrixDouble>();
        pip.push_back(
            new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
        pip.push_back(new OpTreeSearch(
            contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
            get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
            nullptr));
 
        // get body id and SDF range
        auto contact_sfd_map_range_ptr =
            boost::make_shared<std::map<int, Range>>(
                get_body_range("CONTACT_SDF", SPACE_DIM - 1));
 
        pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU(
            contactDisp, contactDisp, contact_common_data_ptr, contactTreeRhs,
            contact_sfd_map_range_ptr));
        pip.push_back(
            new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP<A, GAUSS>(
                contactDisp, broken_data_ptr, contact_common_data_ptr,
                contactTreeRhs, contact_sfd_map_range_ptr));
        pip.push_back(
            new EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU<A, GAUSS>(
                broken_data_ptr, contactDisp, contact_common_data_ptr,
                contactTreeRhs));
 
        MoFEMFunctionReturn(0);
      };
 
      // push ops to face/side pipeline
      CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
      CHKERR add_ops_contact_lhs(op_loop_skeleton_side->getOpPtrVector());
 
      // Add skeleton to domain pipeline
      pip.push_back(op_loop_skeleton_side);
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_hybridisation(fe_lhs->getOpPtrVector());
    CHKERR set_contact(fe_lhs->getOpPtrVector());
  }
 
  if (add_material) {
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setFaceElementOps(
    const bool add_elastic, const bool add_material,
    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
  MoFEMFunctionBegin;
 
  fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
  fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
 
  // set integration rule
  // fe_rhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
  // fe_lhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
  fe_rhs->getRuleHook = [](int, int, int) { return -1; };
  fe_lhs->getRuleHook = [](int, int, int) { return -1; };
  fe_rhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
  fe_lhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
 
  CHKERR
  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      fe_rhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
  CHKERR
  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      fe_lhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
 
  if (add_elastic) {
 
    auto get_broken_op_side = [this](auto &pip) {
      using EleOnSide =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
      using SideEleOp = EleOnSide::UserDataOperator;
      // Iterate over domain FEs adjacent to boundary.
      auto broken_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
      // Note: EleOnSide, i.e. uses on domain projected skeleton rule
      auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
          mField, elementVolumeName, SPACE_DIM, Sev::noisy);
      op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
          boost::make_shared<CGGUserPolynomialBase>();
      CHKERR
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
          materialH1Positions, frontAdjEdges);
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
      boost::shared_ptr<double> piola_scale_ptr(new double);
      op_loop_domain_side->getOpPtrVector().push_back(
          physicalEquations->returnOpSetScale(piola_scale_ptr,
                                              physicalEquations));
 
      auto piola_stress_mat_ptr = boost::make_shared<MatrixDouble>();
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateHVecTensorField<3, 3>(piolaStress,
                                               piola_stress_mat_ptr));
      pip.push_back(op_loop_domain_side);
      return std::make_tuple(broken_data_ptr, piola_scale_ptr,
                             piola_stress_mat_ptr);
    };
 
    auto set_rhs = [&]() {
      MoFEMFunctionBegin;
 
      auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
          get_broken_op_side(fe_rhs->getOpPtrVector());
 
      fe_rhs->getOpPtrVector().push_back(
          new OpDispBc(broken_data_ptr, bcSpatialDispVecPtr, timeScaleMap));
      fe_rhs->getOpPtrVector().push_back(new OpAnalyticalDispBc(
          broken_data_ptr, bcSpatialAnalyticalDisplacementVecPtr,
          timeScaleMap));
      fe_rhs->getOpPtrVector().push_back(new OpRotationBc(
          broken_data_ptr, bcSpatialRotationVecPtr, timeScaleMap));
 
      fe_rhs->getOpPtrVector().push_back(
          new OpBrokenTractionBc(hybridSpatialDisp, bcSpatialTractionVecPtr,
                                 piola_scale_ptr, timeScaleMap));
      auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
      // if you push gradient of L2 base to physical element, it will not work.
      fe_rhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM - 1>(
              hybridSpatialDisp, hybrid_grad_ptr));
      fe_rhs->getOpPtrVector().push_back(new OpBrokenPressureBc(
          hybridSpatialDisp, bcSpatialPressureVecPtr, piola_scale_ptr,
          hybrid_grad_ptr, timeScaleMap));
      fe_rhs->getOpPtrVector().push_back(new OpBrokenAnalyticalTractionBc(
          hybridSpatialDisp, bcSpatialAnalyticalTractionVecPtr, piola_scale_ptr,
          timeScaleMap));
 
      auto hybrid_ptr = boost::make_shared<MatrixDouble>();
      fe_rhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(hybridSpatialDisp,
                                                      hybrid_ptr));
      fe_rhs->getOpPtrVector().push_back(new OpNormalDispRhsBc(
          hybridSpatialDisp, hybrid_ptr, piola_stress_mat_ptr,
          bcSpatialNormalDisplacementVecPtr, timeScaleMap));
 
      auto get_normal_disp_bc_faces = [&]() {
        auto faces =
            get_range_from_block(mField, "NORMAL_DISPLACEMENT", SPACE_DIM - 1);
        return boost::make_shared<Range>(faces);
      };
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using BdyEleOp = BoundaryEle::UserDataOperator;
      using OpC_dBroken = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
          GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
      fe_rhs->getOpPtrVector().push_back(new OpC_dBroken(
          broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0),
          get_normal_disp_bc_faces()));
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_lhs = [&]() {
      MoFEMFunctionBegin;
 
      auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
          get_broken_op_side(fe_lhs->getOpPtrVector());
 
      fe_lhs->getOpPtrVector().push_back(new OpNormalDispLhsBc_dU(
          hybridSpatialDisp, bcSpatialNormalDisplacementVecPtr, timeScaleMap));
      fe_lhs->getOpPtrVector().push_back(new OpNormalDispLhsBc_dP(
          hybridSpatialDisp, broken_data_ptr, bcSpatialNormalDisplacementVecPtr,
          timeScaleMap));
 
      auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
      // if you push gradient of L2 base to physical element, it will not work.
      fe_rhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM - 1>(
              hybridSpatialDisp, hybrid_grad_ptr));
      fe_lhs->getOpPtrVector().push_back(new OpBrokenPressureBcLhs_dU(
          hybridSpatialDisp, bcSpatialPressureVecPtr, hybrid_grad_ptr,
          timeScaleMap));
 
      auto get_normal_disp_bc_faces = [&]() {
        auto faces =
            get_range_from_block(mField, "NORMAL_DISPLACEMENT", SPACE_DIM - 1);
        return boost::make_shared<Range>(faces);
      };
      
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using BdyEleOp = BoundaryEle::UserDataOperator;
      using OpC = FormsIntegrators<BdyEleOp>::Assembly<A>::BiLinearForm<
          GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
      fe_lhs->getOpPtrVector().push_back(new OpC(
          hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0),
          true, true, get_normal_disp_bc_faces()));
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_rhs();
    CHKERR set_lhs();
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setContactElementRhsOps(
 
    boost::shared_ptr<ContactTree> &fe_contact_tree
 
) {
  MoFEMFunctionBegin;
 
  /** Contact requires that body is marked */
  auto get_body_range = [this](auto name, int dim, auto sev) {
    std::map<int, Range> map;
 
    for (auto m_ptr :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % name).str()
 
                 ))
 
    ) {
      Range ents;
      CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
                                                            dim, ents, true),
                     "by dim");
      map[m_ptr->getMeshsetId()] = ents;
      MOFEM_LOG("EPSYNC", sev) << "Meshset: " << m_ptr->getMeshsetId() << " "
                               << ents.size() << " entities";
    }
 
    MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), sev);
    return map;
  };
 
  auto get_map_skin = [this](auto &&map) {
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
 
    Skinner skin(&mField.get_moab());
    for (auto &m : map) {
      Range skin_faces;
      CHKERR skin.find_skin(0, m.second, false, skin_faces);
      CHK_MOAB_THROW(pcomm->filter_pstatus(skin_faces,
                                           PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                           PSTATUS_NOT, -1, nullptr),
                     "filter");
      m.second.swap(skin_faces);
    }
    return map;
  };
 
  /* The above code is written in C++ and it appears to be defining and using
  various operations on boundary elements and side elements. */
  using BoundaryEle = FaceElementForcesAndSourcesCore;
  using BoundaryEleOp = BoundaryEle::UserDataOperator;
 
  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
 
  auto calcs_side_traction = [&](auto &pip) {
    MoFEMFunctionBegin;
    using EleOnSide =
        PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
    using SideEleOp = EleOnSide::UserDataOperator;
    auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
        mField, elementVolumeName, SPACE_DIM, Sev::noisy);
    op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
        boost::make_shared<CGGUserPolynomialBase>();
    EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
        materialH1Positions, frontAdjEdges);
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
            piolaStress, contact_common_data_ptr->contactTractionPtr(),
            boost::make_shared<double>(1.0)));
    pip.push_back(op_loop_domain_side);
    MoFEMFunctionReturn(0);
  };
 
  auto add_contact_three = [&]() {
    MoFEMFunctionBegin;
    auto tree_moab_ptr = boost::make_shared<moab::Core>();
    fe_contact_tree = boost::make_shared<ContactTree>(
        mField, tree_moab_ptr, spaceOrder,
        get_body_range("CONTACT", SPACE_DIM - 1, Sev::inform));
    fe_contact_tree->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(
            contactDisp, contact_common_data_ptr->contactDispPtr()));
    CHKERR calcs_side_traction(fe_contact_tree->getOpPtrVector());
    auto u_h1_ptr = boost::make_shared<MatrixDouble>();
    fe_contact_tree->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
    fe_contact_tree->getOpPtrVector().push_back(
        new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_contact_three();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setElasticElementOps(const int tag) {
  MoFEMFunctionBegin;
 
  // Add contact operators. Note that only for rhs. THe lhs is assembled with
  // volume element, to enable schur complement evaluation.
  CHKERR setContactElementRhsOps(contactTreeRhs);
 
  CHKERR setVolumeElementOps(tag, true, false, elasticFeRhs, elasticFeLhs);
  CHKERR setFaceElementOps(true, false, elasticBcRhs, elasticBcLhs);
 
  auto adj_cache =
      boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
 
  auto get_op_contact_bc = [&]() {
    using SideEle = FaceElementForcesAndSourcesCore;
    auto op_loop_side = new OpLoopSide<SideEle>(
        mField, contactElement, SPACE_DIM - 1, Sev::noisy, adj_cache);
    return op_loop_side;
  };
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setElasticElementToTs(DM dm) {
  MoFEMFunctionBegin;
  boost::shared_ptr<FEMethod> null;
 
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
 
    CHKERR DMMoFEMTSSetI2Function(dm, elementVolumeName, elasticFeRhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Function(dm, naturalBcElement, elasticBcRhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
                                  null);
 
  } else {
    CHKERR DMMoFEMTSSetIFunction(dm, elementVolumeName, elasticFeRhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIFunction(dm, naturalBcElement, elasticBcRhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
                                 null);
  }
 
  MoFEMFunctionReturn(0);
}
 
#include "impl/EshelbianMonitor.cpp"
#include "impl/SetUpSchurImpl.cpp"
#include "impl/TSElasticPostStep.cpp"
#include "impl/CGGUserPolynomialBase.cpp"
 
struct solve_elastic_setup {
 
  inline static auto setup(EshelbianCore *ep_ptr, TS ts, Vec x,
                           bool set_ts_monitor) {
 
#ifdef ENABLE_PYTHON_BINDING
    auto setup_sdf = [&]() {
      boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
 
      auto file_exists = [](std::string myfile) {
        std::ifstream file(myfile.c_str());
        if (file) {
          return true;
        }
        return false;
      };
 
      char sdf_file_name[255] = "sdf.py";
      CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-sdf_file",
                                   sdf_file_name, 255, PETSC_NULLPTR);
 
      if (file_exists(sdf_file_name)) {
        MOFEM_LOG("EP", Sev::inform) << sdf_file_name << " file found";
        sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
        CHKERR sdf_python_ptr->sdfInit(sdf_file_name);
        ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
        MOFEM_LOG("EP", Sev::inform) << "SdfPython initialized";
      } else {
        MOFEM_LOG("EP", Sev::warning) << sdf_file_name << " file NOT found";
      }
      return sdf_python_ptr;
    };
    auto sdf_python_ptr = setup_sdf();
#endif
 
    auto setup_ts_monitor = [&]() {
      boost::shared_ptr<TsCtx> ts_ctx;
      DMMoFEMGetTsCtx(ep_ptr->dmElastic, ts_ctx);
      if (set_ts_monitor) {
        CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULLPTR);
        auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*ep_ptr);
        ts_ctx->getLoopsMonitor().push_back(
            TsCtx::PairNameFEMethodPtr(ep_ptr->elementVolumeName, monitor_ptr));
      }
      MOFEM_LOG("EP", Sev::inform) << "TS monitor setup";
      return std::make_tuple(ts_ctx);
    };
 
    auto setup_snes_monitor = [&]() {
      MoFEMFunctionBeginHot;
      SNES snes;
      CHKERR TSGetSNES(ts, &snes);
      auto snes_ctx = getDMSnesCtx(ep_ptr->dmElastic);
      CHKERR SNESMonitorSet(snes,
                            (MoFEMErrorCode (*)(SNES, PetscInt, PetscReal,
                                                void *))MoFEMSNESMonitorEnergy,
                            (void *)(snes_ctx.get()), PETSC_NULLPTR);
      MOFEM_LOG("EP", Sev::inform) << "SNES monitor setup";
      MoFEMFunctionReturnHot(0);
    };
 
    auto setup_snes_conergence_test = [&]() {
      MoFEMFunctionBegin;
 
      auto snes_convergence_test = [](SNES snes, PetscInt it, PetscReal xnorm,
                                      PetscReal snorm, PetscReal fnorm,
                                      SNESConvergedReason *reason, void *cctx) {
        MoFEMFunctionBegin;
        // EshelbianCore *ep_ptr = (EshelbianCore *)cctx;
        CHKERR SNESConvergedDefault(snes, it, xnorm, snorm, fnorm, reason,
                                    PETSC_NULLPTR);
 
        Vec x_update, r;
        CHKERR SNESGetSolutionUpdate(snes, &x_update);
        CHKERR SNESGetFunction(snes, &r, PETSC_NULLPTR, PETSC_NULLPTR);
 
        // *reason = SNES_CONVERGED_ITERATING;
        // if (!it) {
        //   /* set parameter for default relative tolerance convergence test */
        //   snes->ttol = fnorm * snes->rtol;
        //   snes->rnorm0 = fnorm;
        // }
        // if (PetscIsInfOrNanReal(fnorm)) {
        //   MOFEM_LOG("EP", Sev::error)
        //       << "SNES convergence test: function norm is NaN";
        //   *reason = SNES_DIVERGED_FNORM_NAN;
        // } else if (fnorm < snes->abstol && (it || !snes->forceiteration)) {
        //   MOFEM_LOG("EP", Sev::inform)
        //       << "SNES convergence test: Converged due to function norm "
        //       << std::setw(14) << std::setprecision(12) << std::scientific
        //       << fnorm << " < " << std::setw(14) << std::setprecision(12)
        //       << std::scientific << snes->abstol;
 
        //   PetscCall(PetscInfo(
        //       snes, "Converged due to function norm %14.12e < %14.12e\n",
        //       (double)fnorm, (double)snes->abstol));
        //   *reason = SNES_CONVERGED_FNORM_ABS;
        // } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
        //   PetscCall(PetscInfo(
        //       snes,
        //       "Exceeded maximum number of function evaluations: %" PetscInt_FMT
        //       " > %" PetscInt_FMT "\n",
        //       snes->nfuncs, snes->max_funcs));
        //   *reason = SNES_DIVERGED_FUNCTION_COUNT;
        // }
 
        // if (it && !*reason) {
        //   if (fnorm <= snes->ttol) {
        //     PetscCall(PetscInfo(snes,
        //                         "Converged due to function norm %14.12e < "
        //                         "%14.12e (relative tolerance)\n",
        //                         (double)fnorm, (double)snes->ttol));
        //     *reason = SNES_CONVERGED_FNORM_RELATIVE;
        //   } else if (snorm < snes->stol * xnorm) {
        //     PetscCall(PetscInfo(snes,
        //                         "Converged due to small update length: %14.12e "
        //                         "< %14.12e * %14.12e\n",
        //                         (double)snorm, (double)snes->stol,
        //                         (double)xnorm));
        //     *reason = SNES_CONVERGED_SNORM_RELATIVE;
        //   } else if (snes->divtol != PETSC_UNLIMITED &&
        //              (fnorm > snes->divtol * snes->rnorm0)) {
        //     PetscCall(PetscInfo(snes,
        //                         "Diverged due to increase in function norm: "
        //                         "%14.12e > %14.12e * %14.12e\n",
        //                         (double)fnorm, (double)snes->divtol,
        //                         (double)snes->rnorm0));
        //     *reason = SNES_DIVERGED_DTOL;
        //   }
        // }
 
        MoFEMFunctionReturn(0);
      };
 
      // SNES snes;
      // CHKERR TSGetSNES(ts, &snes);
      // CHKERR SNESSetConvergenceTest(snes, snes_convergence_test, ep_ptr,
      //                               PETSC_NULLPTR);
      // MOFEM_LOG("EP", Sev::inform) << "SNES convergence test setup";
      MoFEMFunctionReturn(0);
    };
 
    auto setup_section = [&]() {
      PetscSection section_raw;
      CHKERR DMGetSection(ep_ptr->dmElastic, &section_raw);
      int num_fields;
      CHKERR PetscSectionGetNumFields(section_raw, &num_fields);
      for (int ff = 0; ff != num_fields; ff++) {
        const char *field_name;
        CHKERR PetscSectionGetFieldName(section_raw, ff, &field_name);
        MOFEM_LOG_C("EP", Sev::inform, "Field %d name %s", ff, field_name);
      }
      return section_raw;
    };
 
    auto set_vector_on_mesh = [&]() {
      MoFEMFunctionBeginHot;
      CHKERR DMoFEMMeshToLocalVector(ep_ptr->dmElastic, x, INSERT_VALUES,
                                     SCATTER_FORWARD);
      CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
      MOFEM_LOG("EP", Sev::inform) << "Vector set on mesh";
      MoFEMFunctionReturnHot(0);
    };
 
    auto setup_schur_block_solver = [&]() {
      MOFEM_LOG("EP", Sev::inform) << "Setting up Schur block solver";
      CHKERR TSAppendOptionsPrefix(ts, "elastic_");
      CHKERR TSSetFromOptions(ts);
      CHKERR TSSetDM(ts, ep_ptr->dmElastic);
      // Adding field split solver
      boost::shared_ptr<EshelbianCore::SetUpSchur> schur_ptr;
      if constexpr (A == AssemblyType::BLOCK_MAT) {
        schur_ptr =
            EshelbianCore::SetUpSchur::createSetUpSchur(ep_ptr->mField, ep_ptr);
        CHK_THROW_MESSAGE(schur_ptr->setUp(ts), "setup schur");
      }
      MOFEM_LOG("EP", Sev::inform) << "Setting up Schur block solver done";
      return schur_ptr;
    };
 
    // Warning: sequence of construction is not guaranteed for tuple. You have
    // to enforce order by proper packaging.
 
#ifdef ENABLE_PYTHON_BINDING
    return std::make_tuple(setup_sdf(), setup_ts_monitor(),
                           setup_snes_monitor(), setup_snes_conergence_test(),
                           setup_section(), set_vector_on_mesh(),
                           setup_schur_block_solver());
#else
    return std::make_tuple(setup_ts_monitor(), setup_snes_monitor(),
                           setup_snes_conergence_test(), setup_section(),
                           set_vector_on_mesh(), setup_schur_block_solver());
#endif
  }
};
 
MoFEMErrorCode EshelbianCore::solveElastic(TS ts, Vec x) {
  MoFEMFunctionBegin;
 
  PetscBool debug_model = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-debug_model", &debug_model,
                             PETSC_NULLPTR);
 
  if (debug_model == PETSC_TRUE) {
    auto ts_ctx_ptr = getDMTsCtx(dmElastic);
    auto post_proc = [&](TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, Vec F,
                         void *ctx) {
      MoFEMFunctionBegin;
 
      SNES snes;
      CHKERR TSGetSNES(ts, &snes);
      int it;
      CHKERR SNESGetIterationNumber(snes, &it);
      std::string file_name = "snes_iteration_" + std::to_string(it) + ".h5m";
      CHKERR postProcessResults(1, file_name, F, u_t);
      std::string file_skel_name =
          "snes_iteration_skel_" + std::to_string(it) + ".h5m";
 
      auto get_material_force_tag = [&]() {
        auto &moab = mField.get_moab();
        Tag tag;
        CHK_MOAB_THROW(moab.tag_get_handle("MaterialForce", tag),
                       "can't get tag");
        return tag;
      };
 
      CHKERR calculateFaceMaterialForce(1, ts);
      CHKERR postProcessSkeletonResults(1, file_skel_name, F,
                                        {get_material_force_tag()});
 
      MoFEMFunctionReturn(0);
    };
    ts_ctx_ptr->tsDebugHook = post_proc;
  }
 
  auto storage = solve_elastic_setup::setup(this, ts, x, true);
 
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
    Vec xx;
    CHKERR VecDuplicate(x, &xx);
    CHKERR VecZeroEntries(xx);
    CHKERR TS2SetSolution(ts, x, xx);
    CHKERR VecDestroy(&xx);
  } else {
    CHKERR TSSetSolution(ts, x);
  }
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
  CHKERR TSSetUp(ts);
  CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
  CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  CHKERR TSElasticPostStep::postStepInitialise(this);
  CHKERR TSSolve(ts, PETSC_NULLPTR);
  CHKERR TSElasticPostStep::postStepDestroy();
  TetPolynomialBase::switchCacheBaseOff<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
 
  SNES snes;
  CHKERR TSGetSNES(ts, &snes);
  int lin_solver_iterations;
  CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
  MOFEM_LOG("EP", Sev::inform)
      << "Number of linear solver iterations " << lin_solver_iterations;
 
  PetscBool test_cook_flg = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_cook", &test_cook_flg,
                             PETSC_NULLPTR);
  if (test_cook_flg) {
    constexpr int expected_lin_solver_iterations = 11;
    if (lin_solver_iterations > expected_lin_solver_iterations)
      SETERRQ(
          PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
          "Expected number of iterations is different than expected %d > %d",
          lin_solver_iterations, expected_lin_solver_iterations);
  }
 
  PetscBool test_sslv116_flag = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_sslv116",
                             &test_sslv116_flag, PETSC_NULLPTR);
 
  if (test_sslv116_flag) {
    double max_val = 0.0;
    double min_val = 0.0;
    auto field_min_max = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
      MoFEMFunctionBeginHot;
      auto ent_type = ent_ptr->getEntType();
      if (ent_type == MBVERTEX) {
        max_val = std::max(ent_ptr->getEntFieldData()[SPACE_DIM - 1], max_val);
        min_val = std::min(ent_ptr->getEntFieldData()[SPACE_DIM - 1], min_val);
      }
      MoFEMFunctionReturnHot(0);
    };
    CHKERR mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
        field_min_max, spatialH1Disp);
 
    double global_max_val = 0.0;
    double global_min_val = 0.0;
    MPI_Allreduce(&max_val, &global_max_val, 1, MPI_DOUBLE, MPI_MAX,
                  mField.get_comm());
    MPI_Allreduce(&min_val, &global_min_val, 1, MPI_DOUBLE, MPI_MIN,
                  mField.get_comm());
    MOFEM_LOG("EP", Sev::inform)
        << "Max " << spatialH1Disp << " value: " << global_max_val;
    MOFEM_LOG("EP", Sev::inform)
        << "Min " << spatialH1Disp << " value: " << global_min_val;
 
    double ref_max_val = 0.00767;
    double ref_min_val = -0.00329;
    if (std::abs(global_max_val - ref_max_val) > 1e-5) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
               "Incorrect max value of the displacement field: %f != %f",
               global_max_val, ref_max_val);
    }
    if (std::abs(global_min_val - ref_min_val) > 4e-5) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
               "Incorrect min value of the displacement field: %f != %f",
               global_min_val, ref_min_val);
    }
  }
 
  CHKERR gettingNorms();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::solveDynamicRelaxation(TS ts, Vec x,
                                                     int start_step,
                                                     double start_time) {
  MoFEMFunctionBegin;
 
  auto storage = solve_elastic_setup::setup(this, ts, x, false);
 
  double final_time = 1;
  double delta_time = 0.1;
  int max_it = 10;
  PetscBool ts_h1_update = PETSC_FALSE;
 
  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options",
                           "none");
 
  CHKERR PetscOptionsScalar("-dynamic_final_time",
                            "dynamic relaxation final time", "", final_time,
                            &final_time, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-dynamic_delta_time",
                            "dynamic relaxation final time", "", delta_time,
                            &delta_time, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-dynamic_max_it", "dynamic relaxation iterations", "",
                         max_it, &max_it, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-dynamic_h1_update", "update each ts step", "",
                          ts_h1_update, &ts_h1_update, PETSC_NULLPTR);
 
  PetscOptionsEnd();
 
  auto setup_ts_monitor = [&]() {
    auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
    return monitor_ptr;
  };
  auto monitor_ptr = setup_ts_monitor();
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
  CHKERR TSSetUp(ts);
  CHKERR TSElasticPostStep::postStepInitialise(this);
 
  double ts_delta_time;
  CHKERR TSGetTimeStep(ts, &ts_delta_time);
 
  if (ts_h1_update) {
    CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
    CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  }
 
  CHKERR TSElasticPostStep::preStepFun(ts);
  CHKERR TSElasticPostStep::postStepFun(ts);
 
  dynamicTime = start_time;
  dynamicStep = start_step;
  monitor_ptr->ts_u = PETSC_NULLPTR;
  monitor_ptr->ts_t = dynamicTime;
  monitor_ptr->ts_step = dynamicStep;
  CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
  MOFEM_LOG("EP", Sev::inform) << "IT " << dynamicStep << " Time "
                               << dynamicTime << " delta time " << delta_time;
  dynamicTime += delta_time;
  ++dynamicStep;
 
  for (; dynamicTime < final_time; dynamicTime += delta_time) {
    MOFEM_LOG("EP", Sev::inform) << "IT " << dynamicStep << " Time "
                                 << dynamicTime << " delta time " << delta_time;
 
    CHKERR TSSetStepNumber(ts, 0);
    CHKERR TSSetTime(ts, 0);
    CHKERR TSSetTimeStep(ts, ts_delta_time);
    if (!ts_h1_update) {
      CHKERR TSElasticPostStep::preStepFun(ts);
    }
    CHKERR TSSolve(ts, PETSC_NULLPTR);
    if (!ts_h1_update) {
      CHKERR TSElasticPostStep::postStepFun(ts);
    }
 
    monitor_ptr->ts_u = PETSC_NULLPTR;
    monitor_ptr->ts_t = dynamicTime;
    monitor_ptr->ts_step = dynamicStep;
    CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
    ++dynamicStep;
    if (dynamicStep > max_it)
      break;
  }
 
  CHKERR TSElasticPostStep::postStepDestroy();
  TetPolynomialBase::switchCacheBaseOff<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setBlockTagsOnSkin() {
  MoFEMFunctionBegin;
 
  auto set_block = [&](auto name, int dim) {
    std::map<int, Range> map;
    auto set_tag_impl = [&](auto name) {
      MoFEMFunctionBegin;
      auto mesh_mng = mField.getInterface<MeshsetsManager>();
      auto bcs = mesh_mng->getCubitMeshsetPtr(
 
          std::regex((boost::format("%s(.*)") % name).str())
 
      );
      for (auto bc : bcs) {
        Range r;
        CHKERR bc->getMeshsetIdEntitiesByDimension(mField.get_moab(), dim, r,
                                                   true);
        map[bc->getMeshsetId()] = r;
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_tag_impl(name);
 
    return std::make_pair(name, map);
  };
 
  auto set_skin = [&](auto &&map) {
    for (auto &m : map.second) {
      auto s = filter_true_skin(mField, get_skin(mField, m.second));
      m.second.swap(s);
    }
    return map;
  };
 
  auto set_tag = [&](auto &&map) {
    Tag th;
    auto name = map.first;
    int def_val[] = {-1};
    CHK_MOAB_THROW(
        mField.get_moab().tag_get_handle(name, 1, MB_TYPE_INTEGER, th,
                                         MB_TAG_SPARSE | MB_TAG_CREAT, def_val),
        "create tag");
    for (auto &m : map.second) {
      int id = m.first;
      CHK_MOAB_THROW(mField.get_moab().tag_clear_data(th, m.second, &id),
                     "clear tag");
    }
    return th;
  };
 
  listTagsToTransfer.push_back(set_tag(set_skin(set_block("BODY", 3))));
  listTagsToTransfer.push_back(set_tag(set_skin(set_block("MAT_ELASTIC", 3))));
  listTagsToTransfer.push_back(
      set_tag(set_skin(set_block("MAT_NEOHOOKEAN", 3))));
  listTagsToTransfer.push_back(set_tag(set_block("CONTACT", 2)));
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EshelbianCore::postProcessResults(const int tag, const std::string file,
                                  Vec f_residual, Vec var_vector,
                                  std::vector<Tag> tags_to_transfer) {
  MoFEMFunctionBegin;
 
  // mark crack surface
  if (crackingOn) {
    Tag th;
    rval = mField.get_moab().tag_get_handle("CRACK", th);
    if (rval == MB_SUCCESS) {
      CHKERR mField.get_moab().tag_delete(th);
    }
    int def_val[] = {0};
    CHKERR mField.get_moab().tag_get_handle(
        "CRACK", 1, MB_TYPE_INTEGER, th, MB_TAG_SPARSE | MB_TAG_CREAT, def_val);
    int mark[] = {1};
    CHKERR mField.get_moab().tag_clear_data(th, *crackFaces, mark);
    tags_to_transfer.push_back(th);
 
    auto get_tag = [&](auto name, auto dim) {
      auto &mob = mField.get_moab();
      Tag tag;
      double def_val[] = {0., 0., 0.};
      CHK_MOAB_THROW(mob.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
                                        MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
                     "create tag");
      return tag;
    };
 
    tags_to_transfer.push_back(get_tag("MaterialForce", 3));
    // tags_to_transfer.push_back(get_tag("GriffithForce", 1));
  }
 
  // add tags to transfer
  for (auto t : listTagsToTransfer) {
    tags_to_transfer.push_back(t);
  }
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
  }
 
  struct exclude_sdf {
    exclude_sdf(Range &&r) : map(r) {}
    bool operator()(FEMethod *fe_method_ptr) {
      auto ent = fe_method_ptr->getFEEntityHandle();
      if (map.find(ent) != map.end()) {
        return false;
      }
      return true;
    }
 
  private:
    Range map;
  };
 
  contactTreeRhs->exeTestHook =
      exclude_sdf(get_range_from_block(mField, "CONTACT_SDF", SPACE_DIM - 1));
 
  CHKERR DMoFEMLoopFiniteElements(dM, contactElement, contactTreeRhs);
 
  auto get_post_proc = [&](auto &post_proc_mesh, auto sense) {
    using FaceEle = FaceElementForcesAndSourcesCore;
    auto post_proc_ptr =
        boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
            mField, post_proc_mesh);
    EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
        frontAdjEdges);
 
    auto domain_ops = [&](auto &fe, int sense) {
      MoFEMFunctionBegin;
      fe.getUserPolynomialBase() =
          boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          fe.getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions,
          frontAdjEdges);
      auto piola_scale_ptr = boost::make_shared<double>(1.0);
      fe.getOpPtrVector().push_back(physicalEquations->returnOpSetScale(
          piola_scale_ptr, physicalEquations));
      fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
          piolaStress, dataAtPts->getApproxPAtPts(), piola_scale_ptr));
      fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
          bubbleField, dataAtPts->getApproxPAtPts(), piola_scale_ptr,
          SmartPetscObj<Vec>(), MBMAXTYPE));
      if (noStretch) {
        fe.getOpPtrVector().push_back(
            physicalEquations->returnOpCalculateStretchFromStress(
                dataAtPts, physicalEquations));
      } else {
        fe.getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<3>(
                stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
      }
      if (var_vector) {
        auto vec = SmartPetscObj<Vec>(var_vector, true);
        fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
            piolaStress, dataAtPts->getVarPiolaPts(),
            boost::make_shared<double>(1), vec));
        fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
            bubbleField, dataAtPts->getVarPiolaPts(),
            boost::make_shared<double>(1), vec, MBMAXTYPE));
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            rotAxis, dataAtPts->getVarRotAxisPts(), vec, MBTET));
        if (noStretch) {
          fe.getOpPtrVector().push_back(
              physicalEquations->returnOpCalculateVarStretchFromStress(
                  dataAtPts, physicalEquations));
        } else {
          fe.getOpPtrVector().push_back(
              new OpCalculateTensor2SymmetricFieldValues<3>(
                  stretchTensor, dataAtPts->getVarLogStreachPts(), vec, MBTET));
        }
      }
 
      fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
      fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
 
      fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
      fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
      fe.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
          spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
      // evaluate derived quantities
      fe.getOpPtrVector().push_back(
          new OpCalculateRotationAndSpatialGradient(dataAtPts));
 
      // evaluate integration points
      fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
          tag, true, false, dataAtPts, physicalEquations));
      if (auto op =
              physicalEquations->returnOpCalculateEnergy(dataAtPts, nullptr)) {
        fe.getOpPtrVector().push_back(op);
        fe.getOpPtrVector().push_back(new OpCalculateEshelbyStress(dataAtPts));
      }
 
      // // post-proc
      using OpPPMap = OpPostProcMapInMoab<
          SPACE_DIM, SPACE_DIM,
          VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator>;
 
      struct OpSidePPMap : public OpPPMap {
        OpSidePPMap(moab::Interface &post_proc_mesh,
                    std::vector<EntityHandle> &map_gauss_pts,
                    DataMapVec data_map_scalar, DataMapMat data_map_vec,
                    DataMapMat data_map_mat, DataMapMat data_symm_map_mat,
                    int sense)
            : OpPPMap(post_proc_mesh, map_gauss_pts, data_map_scalar,
                      data_map_vec, data_map_mat, data_symm_map_mat),
              tagSense(sense) {}
 
        MoFEMErrorCode doWork(int side, EntityType type,
                              EntitiesFieldData::EntData &data) {
          MoFEMFunctionBegin;
 
          if (tagSense != OpPPMap::getSkeletonSense())
            MoFEMFunctionReturnHot(0);
 
          CHKERR OpPPMap::doWork(side, type, data);
          MoFEMFunctionReturn(0);
        }
 
      private:
        int tagSense;
      };
 
      OpPPMap::DataMapMat vec_fields;
      vec_fields["SpatialDisplacementL2"] = dataAtPts->getSmallWL2AtPts();
      vec_fields["SpatialDisplacementH1"] = dataAtPts->getSmallWH1AtPts();
      vec_fields["Omega"] = dataAtPts->getRotAxisAtPts();
      vec_fields["ContactDisplacement"] = dataAtPts->getContactL2AtPts();
      vec_fields["AngularMomentum"] = dataAtPts->getLeviKirchhoffAtPts();
      vec_fields["X"] = dataAtPts->getLargeXH1AtPts();
      if (!noStretch) {
        vec_fields["EiegnLogStreach"] = dataAtPts->getEigenValsAtPts();
      }
      if (var_vector) {
        auto vec = SmartPetscObj<Vec>(var_vector, true);
        vec_fields["VarOmega"] = dataAtPts->getVarRotAxisPts();
        vec_fields["VarSpatialDisplacementL2"] =
            boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            spatialL2Disp, vec_fields["VarSpatialDisplacementL2"], vec, MBTET));
      }
      if (f_residual) {
        auto vec = SmartPetscObj<Vec>(f_residual, true);
        vec_fields["ResSpatialDisplacementL2"] =
            boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            spatialL2Disp, vec_fields["ResSpatialDisplacementL2"], vec, MBTET));
        vec_fields["ResOmega"] = boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            rotAxis, vec_fields["ResOmega"], vec, MBTET));
      }
 
      OpPPMap::DataMapMat mat_fields;
      mat_fields["PiolaStress"] = dataAtPts->getApproxPAtPts();
      if (var_vector) {
        mat_fields["VarPiolaStress"] = dataAtPts->getVarPiolaPts();
      }
      if (f_residual) {
        auto vec = SmartPetscObj<Vec>(f_residual, true);
        mat_fields["ResPiolaStress"] = boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
            piolaStress, mat_fields["ResPiolaStress"],
            boost::make_shared<double>(1), vec));
        fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
            bubbleField, mat_fields["ResPiolaStress"],
            boost::make_shared<double>(1), vec, MBMAXTYPE));
      }
      if (!internalStressTagName.empty()) {
        mat_fields[internalStressTagName] = dataAtPts->getInternalStressAtPts();
        switch (internalStressInterpOrder) {
        case 0:
          fe.getOpPtrVector().push_back(
              new OpGetInternalStress<0>(dataAtPts, internalStressTagName));
          break;
        case 1:
          fe.getOpPtrVector().push_back(
              new OpGetInternalStress<1>(dataAtPts, internalStressTagName));
          break;
        default:
          SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
                   "Unsupported internal stress interpolation order %d",
                   internalStressInterpOrder);
        }
      }
 
      OpPPMap::DataMapMat mat_fields_symm;
      mat_fields_symm["LogSpatialStretch"] =
          dataAtPts->getLogStretchTensorAtPts();
      mat_fields_symm["SpatialStretch"] = dataAtPts->getStretchTensorAtPts();
      if (var_vector) {
        mat_fields_symm["VarLogSpatialStretch"] =
            dataAtPts->getVarLogStreachPts();
      }
      if (f_residual) {
        auto vec = SmartPetscObj<Vec>(f_residual, true);
        if (!noStretch) {
          mat_fields_symm["ResLogSpatialStretch"] =
              boost::make_shared<MatrixDouble>();
          fe.getOpPtrVector().push_back(
              new OpCalculateTensor2SymmetricFieldValues<3>(
                  stretchTensor, mat_fields_symm["ResLogSpatialStretch"], vec,
                  MBTET));
        }
      }
 
      fe.getOpPtrVector().push_back(
 
          new OpSidePPMap(
 
              post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
              {},
 
              vec_fields,
 
              mat_fields,
 
              mat_fields_symm,
 
              sense
 
              )
 
      );
 
      fe.getOpPtrVector().push_back(new OpPostProcDataStructure(
          post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
          dataAtPts, sense));
 
      MoFEMFunctionReturn(0);
    };
 
    post_proc_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(contactDisp,
                                            dataAtPts->getContactL2AtPts()));
 
    auto X_h1_ptr = boost::make_shared<MatrixDouble>();
    // H1 material positions
    post_proc_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(materialH1Positions,
                                            dataAtPts->getLargeXH1AtPts()));
 
    // domain
    auto op_loop_side = new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
        mField, elementVolumeName, SPACE_DIM);
    domain_ops(*(op_loop_side->getSideFEPtr()), sense);
    post_proc_ptr->getOpPtrVector().push_back(op_loop_side);
 
    return post_proc_ptr;
  };
 
  // contact
  auto calcs_side_traction_and_displacements = [&](auto &post_proc_ptr,
                                                   auto &pip) {
    MoFEMFunctionBegin;
    auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
    // evaluate traction
    using EleOnSide =
        PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
    using SideEleOp = EleOnSide::UserDataOperator;
    auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
        mField, elementVolumeName, SPACE_DIM, Sev::noisy);
    op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
        boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
    EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
        materialH1Positions, frontAdjEdges);
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
            piolaStress, contact_common_data_ptr->contactTractionPtr(),
            boost::make_shared<double>(1.0)));
    pip.push_back(op_loop_domain_side);
    // evaluate contact displacement and contact conditions
    auto u_h1_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
        contactDisp, contact_common_data_ptr->contactDispPtr()));
    pip.push_back(new OpTreeSearch(
        contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
        get_range_from_block(mField, "CONTACT", SPACE_DIM - 1),
        &post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
 
    if (f_residual) {
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      auto op_this = new OpLoopThis<BoundaryEle>(mField, contactElement);
      pip.push_back(op_this);
      auto contact_residual = boost::make_shared<MatrixDouble>();
      op_this->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(
              contactDisp, contact_residual,
              SmartPetscObj<Vec>(f_residual, true)));
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
      op_this->getOpPtrVector().push_back(
 
          new OpPPMap(
 
              post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
              {},
 
              {{"res_contact", contact_residual}},
 
              {},
 
              {}
 
              )
 
      );
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto post_proc_mesh = boost::make_shared<moab::Core>();
  auto post_proc_ptr = get_post_proc(post_proc_mesh, 1);
  auto post_proc_negative_sense_ptr = get_post_proc(post_proc_mesh, -1);
  CHKERR calcs_side_traction_and_displacements(post_proc_ptr,
                                               post_proc_ptr->getOpPtrVector());
 
  auto own_tets =
      CommInterface::getPartEntities(mField.get_moab(), mField.get_comm_rank())
          .subset_by_dimension(SPACE_DIM);
  Range own_faces;
  CHKERR mField.get_moab().get_adjacencies(own_tets, SPACE_DIM - 1, true,
                                           own_faces, moab::Interface::UNION);
 
  auto get_post_negative = [&](auto &&ents) {
    auto crack_faces_pos = ents;
    auto crack_faces_neg = crack_faces_pos;
    auto skin = get_skin(mField, own_tets);
    auto crack_on_proc_skin = intersect(crack_faces_pos, skin);
    for (auto f : crack_on_proc_skin) {
      Range tet;
      CHKERR mField.get_moab().get_adjacencies(&f, 1, SPACE_DIM, false, tet);
      tet = intersect(tet, own_tets);
      int side_number, sense, offset;
      CHKERR mField.get_moab().side_number(tet[0], f, side_number, sense,
                                           offset);
      if (sense == 1) {
        crack_faces_neg.erase(f);
      } else {
        crack_faces_pos.erase(f);
      }
    }
    return std::make_pair(crack_faces_pos, crack_faces_neg);
  };
 
  auto get_crack_faces = [&](auto crack_faces) {
    auto get_adj = [&](auto e, auto dim) {
      Range adj;
      CHKERR mField.get_moab().get_adjacencies(e, dim, true, adj,
                                               moab::Interface::UNION);
      return adj;
    };
    auto tets = get_adj(crack_faces, 3);
    auto faces = subtract(get_adj(tets, 2), crack_faces);
    tets = subtract(tets, get_adj(faces, 3));
    return subtract(crack_faces, get_adj(tets, 2));
  };
 
  auto [crack_faces_pos, crack_faces_neg] =
      get_post_negative(intersect(own_faces, get_crack_faces(*crackFaces)));
 
  auto only_crack_faces = [&](FEMethod *fe_method_ptr) {
    auto ent = fe_method_ptr->getFEEntityHandle();
    if (crack_faces_pos.find(ent) == crack_faces_pos.end()) {
      return false;
    }
    return true;
  };
 
  auto only_negative_crack_faces = [&](FEMethod *fe_method_ptr) {
    auto ent = fe_method_ptr->getFEEntityHandle();
    if (crack_faces_neg.find(ent) == crack_faces_neg.end()) {
      return false;
    }
    return true;
  };
 
  auto get_adj_front = [&]() {
    auto skeleton_faces = *skeletonFaces;
    Range adj_front;
    CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2, true, adj_front,
                                             moab::Interface::UNION);
 
    adj_front = intersect(adj_front, skeleton_faces);
    adj_front = subtract(adj_front, *crackFaces);
    adj_front = intersect(own_faces, adj_front);
    return adj_front;
  };
 
  post_proc_ptr->setTagsToTransfer(tags_to_transfer);
  post_proc_negative_sense_ptr->setTagsToTransfer(tags_to_transfer);
 
  auto post_proc_begin =
      PostProcBrokenMeshInMoabBaseBegin(mField, post_proc_mesh);
  CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin.getFEMethod());
  CHKERR DMoFEMLoopFiniteElements(dM, skinElement, post_proc_ptr);
  post_proc_ptr->exeTestHook = only_crack_faces;
  post_proc_negative_sense_ptr->exeTestHook = only_negative_crack_faces;
  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
      dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, skeletonElement,
                                              post_proc_negative_sense_ptr, 0,
                                              mField.get_comm_size());
 
  constexpr bool debug = false;
  if (debug) {
 
    auto [adj_front_pos, adj_front_neg] =
        get_post_negative(filter_owners(mField, get_adj_front()));
 
    auto only_front_faces_pos = [adj_front_pos](FEMethod *fe_method_ptr) {
      auto ent = fe_method_ptr->getFEEntityHandle();
      if (adj_front_pos.find(ent) == adj_front_pos.end()) {
        return false;
      }
      return true;
    };
 
    auto only_front_faces_neg = [adj_front_neg](FEMethod *fe_method_ptr) {
      auto ent = fe_method_ptr->getFEEntityHandle();
      if (adj_front_neg.find(ent) == adj_front_neg.end()) {
        return false;
      }
      return true;
    };
 
    post_proc_ptr->exeTestHook = only_front_faces_pos;
    CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
        dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
    post_proc_negative_sense_ptr->exeTestHook = only_front_faces_neg;
    CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, skeletonElement,
                                                post_proc_negative_sense_ptr, 0,
                                                mField.get_comm_size());
  }
  auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
  CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
 
  CHKERR post_proc_end.writeFile(file.c_str());
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EshelbianCore::postProcessSkeletonResults(const int tag, const std::string file,
                                          Vec f_residual,
                                          std::vector<Tag> tags_to_transfer) {
  MoFEMFunctionBegin;
 
  using FaceEle = FaceElementForcesAndSourcesCore;
 
  auto post_proc_mesh = boost::make_shared<moab::Core>();
  auto post_proc_ptr =
      boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
          mField, post_proc_mesh);
  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM - 1, SPACE_DIM>::add(
      post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
      frontAdjEdges);
 
  auto hybrid_disp = boost::make_shared<MatrixDouble>();
  post_proc_ptr->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<3>(hybridSpatialDisp, hybrid_disp));
  post_proc_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);
  post_proc_ptr->getOpPtrVector().push_back(op_loop_domain_side);
 
  // evaluated in side domain, that is op_loop_domain_side
  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>();
  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
      materialH1Positions, frontAdjEdges);
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
          piolaStress, dataAtPts->getApproxPAtPts()));
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpCalculateHTensorTensorField<SPACE_DIM, SPACE_DIM>(
          bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(
          rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  if (noStretch) {
    op_loop_domain_side->getOpPtrVector().push_back(
        physicalEquations->returnOpCalculateStretchFromStress(
            dataAtPts, physicalEquations));
  } else {
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
            stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
  }
 
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
  OpPPMap::DataMapMat vec_fields;
  vec_fields["HybridDisplacement"] = hybrid_disp;
  vec_fields["Omega"] = dataAtPts->getRotAxisAtPts();
  OpPPMap::DataMapMat mat_fields;
  mat_fields["PiolaStress"] = dataAtPts->getApproxPAtPts();
  mat_fields["HybridDisplacementGradient"] =
      dataAtPts->getGradHybridDispAtPts();
  OpPPMap::DataMapMat mat_fields_symm;
  mat_fields_symm["LogSpatialStretch"] = dataAtPts->getLogStretchTensorAtPts();
 
  post_proc_ptr->getOpPtrVector().push_back(
 
      new OpPPMap(
 
          post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
          {},
 
          vec_fields,
 
          mat_fields,
 
          mat_fields_symm
 
          )
 
  );
 
  if (f_residual) {
    auto hybrid_res = boost::make_shared<MatrixDouble>();
    post_proc_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(
            hybridSpatialDisp, hybrid_res,
            SmartPetscObj<Vec>(f_residual, true)));
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
    post_proc_ptr->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
            {},
 
            {{"res_hybrid", hybrid_res}},
 
            {},
 
            {}
 
            )
 
    );
  }
 
  post_proc_ptr->setTagsToTransfer(tags_to_transfer);
 
  auto post_proc_begin =
      PostProcBrokenMeshInMoabBaseBegin(mField, post_proc_mesh);
  CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin.getFEMethod());
  CHKERR DMoFEMLoopFiniteElements(dM, skeletonElement, post_proc_ptr);
  auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
  CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
 
  CHKERR post_proc_end.writeFile(file.c_str());
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::calculateFaceMaterialForce(const int tag, TS ts) {
  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, AREAGROWTH, GRIFFITHFORCE, FACEPRESSURE };
 
  auto tags = get_tags_vec({{"MaterialForce", 3},
                            {"AreaGrowth", 3},
                            {"GriffithForce", 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);
 
      // hybrid disp, evalated 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>();
 
      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 OpCalculateHTensorTensorField<SPACE_DIM, SPACE_DIM>(
              bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
      // op_loop_domain_side->getOpPtrVector().push_back(
      //     new OpCalculateHVecTensorGradient<3, 9, 3>(
      //         piolaStress, dataAtPts->getGradPAtPts()));
      // op_loop_domain_side->getOpPtrVector().push_back(
      //     new OpCalculateHVecTensorGradient<9, 9, 3>(
      //         bubbleField, dataAtPts->getGradPAtPts(), MBMAXTYPE));
 
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(
              rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
      if (noStretch) {
        op_loop_domain_side->getOpPtrVector().push_back(
            physicalEquations->returnOpCalculateStretchFromStress(
                dataAtPts, physicalEquations));
      } else {
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
                stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
      }
 
      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 calculate_front_material_force = [&](auto nb_J_integral_contours) {
    MoFEMFunctionBegin;
    FTENSOR_INDEX(SPACE_DIM, I);
 
    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_material_force = [&](auto r, auto th) {
      MatrixDouble material_forces(r.size(), 3, false);
      CHK_MOAB_THROW(
          mField.get_moab().tag_get_data(th, r, material_forces.data().data()),
          "get data");
      return material_forces;
    };
 
    if (mField.get_comm_rank() == 0) {
 
      auto crack_edges = get_adj_range(*crackFaces, 1);
      auto front_nodes = get_conn_range(*frontEdges);
      auto body_edges = get_range_from_block(mField, "EDGES", 1);
      // auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
      // front_block_edges = subtract(front_block_edges, crack_edges);
      // auto front_block_edges_conn = get_conn_range(front_block_edges);
 
      // #ifndef NDEBUG
      Range all_skin_faces;
      Range all_front_faces;
      // #endif
 
      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;
      };
 
      // take bubble tets at node, and then avarage over the edges
      auto calculate_force_through_node = [&]() {
        MoFEMFunctionBegin;
 
        FTENSOR_INDEX(SPACE_DIM, I);
        FTENSOR_INDEX(SPACE_DIM, J);
        FTENSOR_INDEX(SPACE_DIM, K);
        FTENSOR_INDEX(SPACE_DIM, L);
 
        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);
 
        // calculate nodal material force
        for (auto n : front_nodes) {
          auto adj_tets = get_adj(n, 3);
          for (int ll = 0; ll < nb_J_integral_contours; ++ll) {
            auto conn = get_conn_range(adj_tets);
            adj_tets = get_adj_range(conn, 3);
          }
          auto skin_faces = get_skin(mField, adj_tets);
          auto material_forces = get_material_force(skin_faces, tags[0]);
 
          #ifndef NDEBUG
          if (debug) {
            all_skin_faces.merge(skin_faces);
          }
          #endif
 
          auto calculate_node_material_force = [&]() {
            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, 3), 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);
            }
 
            return t_node_force;
          };
 
          auto calculate_crack_area_growth_direction =
              [&](auto n, auto &t_node_force) {
                // 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, 2), 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);
                }
 
                auto adj_faces = intersect(get_adj(n, 2), *crackFaces);
                if (adj_faces.empty()) {
                  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;
                }
 
                // calculate direction
                auto front_edges = get_adj(n, 1);
                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 t_node_force = calculate_node_material_force();
          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");
 
          auto get_area_dir = [&]() {
            FTensor::Tensor1<double, SPACE_DIM> t_area_dir{0., 0., 0.};
            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);
            t_area_dir(I) = 0;
            for (auto sn : seed_n) {
              auto t_area_dir_sn =
                  calculate_crack_area_growth_direction(sn, t_node_force);
              t_area_dir(I) += t_area_dir_sn(I);
            }
            for (auto sn : skin_adj_edges) {
              auto t_area_dir_sn =
                  calculate_crack_area_growth_direction(sn, t_node_force);
              t_area_dir(I) += t_area_dir_sn(I) / 2;
            }
            return t_area_dir;
          };
 
          auto t_area_dir = get_area_dir();
 
          CHK_MOAB_THROW(
              mField.get_moab().tag_set_data(tags[ExhangeTags::AREAGROWTH], &n,
                                             1, &t_area_dir(0)),
              "set data");
          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[ExhangeTags::GRIFFITHFORCE],
                                             &n, 1, &griffith),
              "set data");
        }
 
        // iterate over edges, and calculate average edge material force
        auto ave_node_force = [&](auto th) {
          MoFEMFunctionBegin;
 
          for (auto e : *frontEdges) {
 
            auto conn = get_conn(e);
            auto data = get_material_force(conn, th);
            auto t_node = getFTensor1FromPtr<SPACE_DIM>(data.data().data());
            FTensor::Tensor1<double, SPACE_DIM> t_edge{0., 0., 0.};
            for (auto n : conn) {
              NOT_USED(n);
              t_edge(I) += t_node(I);
              ++t_node;
            }
            t_edge(I) /= conn.size();
 
            FTensor::Tensor1<double, SPACE_DIM> t_edge_direction =
                calculate_edge_direction(e);
            t_edge_direction.normalize();
 
            FTENSOR_INDEX(SPACE_DIM, I);
            FTENSOR_INDEX(SPACE_DIM, J);
            FTENSOR_INDEX(SPACE_DIM, K);
            FTensor::Tensor1<double, SPACE_DIM> t_cross;
            t_cross(K) =
                FTensor::levi_civita(I, J, K) * t_edge_direction(I) * t_edge(J);
            t_edge(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(0));
          }
          MoFEMFunctionReturn(0);
        };
 
        // iterate over edges, and calculate average edge griffith energy
        auto ave_node_griffith_energy = [&](auto th) {
          MoFEMFunctionBegin;
          for (auto e : *frontEdges) {
            FTensor::Tensor1<double, SPACE_DIM> t_edge_force;
            CHKERR mField.get_moab().tag_get_data(
                tags[ExhangeTags::MATERIALFORCE], &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(th, &e, 1, &griffith_energy);
          }
          MoFEMFunctionReturn(0);
        };
 
        CHKERR ave_node_force(tags[ExhangeTags::MATERIALFORCE]);
        CHKERR ave_node_force(tags[ExhangeTags::AREAGROWTH]);
        CHKERR ave_node_griffith_energy(tags[ExhangeTags::GRIFFITHFORCE]);
 
        MoFEMFunctionReturn(0);
      };
 
      CHKERR calculate_force_through_node();
 
      // calculate face cross
      for (auto e : *frontEdges) {
        auto adj_faces = get_adj(e, 2);
        auto crack_face = intersect(get_adj(e, 2), *crackFaces);
 
        // #ifndef NDEBUG
        if (debug) {
          all_front_faces.merge(adj_faces);
        }
        // #endif
 
        FTensor::Tensor1<double, SPACE_DIM> t_edge_force;
        CHKERR mField.get_moab().tag_get_data(tags[ExhangeTags::MATERIALFORCE],
                                              &e, 1, &t_edge_force(0));
        FTensor::Tensor1<double, SPACE_DIM> t_edge_direction =
            calculate_edge_direction(e);
        t_edge_direction.normalize();
        FTENSOR_INDEX(SPACE_DIM, I);
        FTENSOR_INDEX(SPACE_DIM, J);
        FTENSOR_INDEX(SPACE_DIM, K);
        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[ExhangeTags::GRIFFITHFORCE], &f, 1, &dot),
                         "set data");
        }
      }
 
      #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);
        CHKERR save_range(mField.get_moab(),
                          "front_skin_faces_material_force_" +
                              std::to_string(ts_step) + ".vtk",
                          all_skin_faces);
        CHKERR save_range(mField.get_moab(),
                          "front_faces_material_force_" +
                              std::to_string(ts_step) + ".vtk",
                          all_front_faces);
      }
      #endif
    }
 
    auto edge_exchange = CommInterface::createEntitiesPetscVector(
        mField.get_comm(), mField.get_moab(), 1, 3, Sev::inform);
    CHKERR CommInterface::updateEntitiesPetscVector(
        mField.get_moab(), edge_exchange, tags[ExhangeTags::MATERIALFORCE]);
    CHKERR CommInterface::updateEntitiesPetscVector(
        mField.get_moab(), edge_exchange, tags[ExhangeTags::AREAGROWTH]);
    CHKERR CommInterface::updateEntitiesPetscVector(
        mField.get_moab(), edgeExchange, tags[ExhangeTags::GRIFFITHFORCE]);
 
    MoFEMFunctionReturn(0);
  };
 
  auto print_results = [&](auto nb_J_integral_conturs) {
    MoFEMFunctionBegin;
 
    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);
        auto material_force =
            get_tag_data(conn, tags[ExhangeTags::MATERIALFORCE], 3);
        auto area_growth = get_tag_data(conn, tags[ExhangeTags::AREAGROWTH], 3);
        auto griffith_force =
            get_tag_data(conn, tags[ExhangeTags::GRIFFITHFORCE], 1);
        std::vector<double> coords(conn.size() * 3);
        CHK_MOAB_THROW(mField.get_moab().get_coords(conn, coords.data()),
                       "get coords");
        if (conn.size())
          MOFEM_LOG("EPSELF", Sev::inform) << "Material force at nodes";
        for (size_t i = 0; i < conn.size(); ++i) {
          MOFEM_LOG("EPSELF", Sev::inform)
              << "Node " << conn[i] << " coords " << coords[i * 3 + 0] << " "
              << coords[i * 3 + 1] << " " << coords[i * 3 + 2]
              << " material force " << material_force[i * 3 + 0] << " "
              << material_force[i * 3 + 1] << " " << material_force[i * 3 + 2]
              << " area growth " << area_growth[i * 3 + 0] << " "
              << area_growth[i * 3 + 1] << " " << area_growth[i * 3 + 2]
              << " griffith force " << std::setprecision(12)
              << griffith_force[i] << " contour " << nb_J_integral_conturs;
        }
        MoFEMFunctionReturn(0);
      };
 
      at_nodes();
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR calculate_material_forces();
 
  PetscBool all_conturs = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "",
                             "-calculate_J_integral_all_levels", &all_conturs,
                             PETSC_NULLPTR); // for backward compatibility
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "",
                             "-calculate_J_integral_all_contours", &all_conturs,
                             PETSC_NULLPTR); // new name
 
 
  if (all_conturs == PETSC_TRUE) {
    for (int l = 0; l < nbJIntegralContours; ++l) {
      CHKERR calculate_front_material_force(l);
      CHKERR print_results(l);
    }
  }
 
  CHKERR calculate_front_material_force(nbJIntegralContours);
  CHKERR print_results(nbJIntegralContours);
 
 
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::calculateOrientation(const int tag,
                                                   bool set_orientation) {
  MoFEMFunctionBegin;
 
  // constexpr bool debug = false;
  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_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;
    };
 
    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) {
      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);
        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;
      };
 
      // calulate 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::inform)
                            << "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;
            sort_by_energy[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}) {
 
          // 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 < 0) {
                  MOFEM_LOG("EPSELF", Sev::inform)
                      << "Skip negative face " << max_face << " with energy "
                      << energy << " and pressure " << face_pressure;
                  continue;
                }
              }
            }
 
            MOFEM_LOG("EPSELF", Sev::inform)
                << "Check face " << max_face << " edge " << max_edge
                << " energy " << energy << " optimal angle " << opt_angle
                << " face pressure " << face_pressure;
 
            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 (!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);
    };
 
    MOFEM_LOG("EP", sev) << "Front edges " << frontEdges->size();
    CHKERR evaluate_face_energy_and_set_orientation(
        *frontEdges, get_adj_front(false), sides_pair, th_front_position);
  }
 
  // 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()) {
      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::saveOrgCoords() {
  MoFEMFunctionBegin;
  auto crack_faces =
      get_range_from_block(mField, "CRACK_COMPUTED", SPACE_DIM - 1);
  Range conn;
  CHKERR mField.get_moab().get_connectivity(crack_faces, conn, true);
  Range verts;
  CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, verts);
  verts = subtract(verts, conn);
  std::vector<double> coords(3 * verts.size());
  CHKERR mField.get_moab().get_coords(verts, coords.data());
  double def_coords[] = {0., 0., 0.};
  Tag th_org_coords;
  CHKERR mField.get_moab().tag_get_handle(
      "ORG_COORDS", 3, MB_TYPE_DOUBLE, th_org_coords,
      MB_TAG_CREAT | MB_TAG_DENSE, def_coords);
  CHKERR mField.get_moab().tag_set_data(th_org_coords, verts, coords.data());
  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);
};
 
//! [Getting norms]
MoFEMErrorCode EshelbianCore::gettingNorms() {
  MoFEMFunctionBegin;
 
  auto post_proc_norm_fe =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
 
  auto post_proc_norm_rule_hook = [](int, int, int p) -> int {
    return 2 * (p);
  };
  post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
 
  post_proc_norm_fe->getUserPolynomialBase() =
      boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
 
  CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV}, materialH1Positions,
      frontAdjEdges);
 
  enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
  auto norms_vec =
      createVectorMPI(mField.get_comm(), LAST_NORM, PETSC_DETERMINE);
  CHKERR VecZeroEntries(norms_vec);
 
  auto u_l2_ptr = boost::make_shared<MatrixDouble>();
  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp, u_l2_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(spatialH1Disp, u_h1_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_NORM_L2));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_h1_ptr, norms_vec, U_NORM_H1));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_ERROR_L2,
                                         u_h1_ptr));
 
  auto piola_ptr = boost::make_shared<MatrixDouble>();
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateHVecTensorField<3, 3>(piolaStress, piola_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateHTensorTensorField<3, 3>(bubbleField, piola_ptr,
                                              MBMAXTYPE));
 
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor2<3, 3>(piola_ptr, norms_vec, PIOLA_NORM));
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>({post_proc_norm_fe.get()});
  CHKERR mField.loop_finite_elements("ELASTIC_PROBLEM", elementVolumeName,
                                     *post_proc_norm_fe);
  TetPolynomialBase::switchCacheBaseOff<HDIV>({post_proc_norm_fe.get()});
 
  CHKERR VecAssemblyBegin(norms_vec);
  CHKERR VecAssemblyEnd(norms_vec);
  const double *norms;
  CHKERR VecGetArrayRead(norms_vec, &norms);
  MOFEM_LOG("EP", Sev::inform) << "norm_u: " << std::sqrt(norms[U_NORM_L2]);
  MOFEM_LOG("EP", Sev::inform) << "norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
  MOFEM_LOG("EP", Sev::inform)
      << "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
  MOFEM_LOG("EP", Sev::inform)
      << "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
  CHKERR VecRestoreArrayRead(norms_vec, &norms);
 
  MoFEMFunctionReturn(0);
}
//! [Getting norms]
 
MoFEMErrorCode EshelbianCore::getSpatialDispBc() {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      "", piolaStress, false, false);
 
  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    if (auto disp_bc = bc.second->dispBcPtr) {
 
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      MOFEM_LOG("EP", Sev::noisy) << "Displacement BC: " << *disp_bc;
 
      std::vector<double> block_attributes(6, 0.);
      if (disp_bc->data.flag1 == 1) {
        block_attributes[0] = disp_bc->data.value1;
        block_attributes[3] = 1;
      }
      if (disp_bc->data.flag2 == 1) {
        block_attributes[1] = disp_bc->data.value2;
        block_attributes[4] = 1;
      }
      if (disp_bc->data.flag3 == 1) {
        block_attributes[2] = disp_bc->data.value3;
        block_attributes[5] = 1;
      }
      auto faces = bc.second->bcEnts.subset_by_dimension(2);
      bcSpatialDispVecPtr->emplace_back(block_name, block_attributes, faces);
    }
  }
  // old way of naming blocksets for displacement BCs
  CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
 
  bcSpatialNormalDisplacementVecPtr =
      boost::make_shared<NormalDisplacementBcVec>();
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)NORMAL_DISPLACEMENT(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialNormalDisplacementVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
  bcSpatialAnalyticalDisplacementVecPtr =
      boost::make_shared<AnalyticalDisplacementBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)ANALYTICAL_DISPLACEMENT(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialAnalyticalDisplacementVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
  auto ts_displacement =
      boost::make_shared<DynamicRelaxationTimeScale>("disp_history.txt");
  for (auto &bc : *bcSpatialDispVecPtr) {
    MOFEM_LOG("EP", Sev::noisy)
        << "Add time scaling displacement BC: " << bc.blockName;
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_displacement, "disp_history", ".txt", bc.blockName);
  }
 
  auto ts_normal_displacement =
      boost::make_shared<DynamicRelaxationTimeScale>("normal_disp_history.txt");
  for (auto &bc : *bcSpatialNormalDisplacementVecPtr) {
    MOFEM_LOG("EP", Sev::noisy)
        << "Add time scaling normal displacement BC: " << bc.blockName;
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_normal_displacement, "normal_disp_history", ".txt",
            bc.blockName);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::getSpatialTractionBc() {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<ForceCubitBcData>("", piolaStress,
                                                          false, false);
 
  bcSpatialTractionVecPtr = boost::make_shared<TractionBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    if (auto force_bc = bc.second->forceBcPtr) {
 
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      MOFEM_LOG("EP", Sev::noisy) << "Force BC: " << *force_bc;
 
      std::vector<double> block_attributes(6, 0.);
      block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
      block_attributes[3] = 1;
      block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
      block_attributes[4] = 1;
      block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
      block_attributes[5] = 1;
      auto faces = bc.second->bcEnts.subset_by_dimension(2);
      bcSpatialTractionVecPtr->emplace_back(block_name, block_attributes,
                                            faces);
    }
  }
  CHKERR getBc(bcSpatialTractionVecPtr, "SPATIAL_TRACTION_BC", 6);
 
  bcSpatialPressureVecPtr = boost::make_shared<PressureBcVec>();
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)PRESSURE(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
 
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialPressureVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
bcSpatialAnalyticalTractionVecPtr =
      boost::make_shared<AnalyticalTractionBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)ANALYTICAL_TRACTION(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialAnalyticalTractionVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
  auto ts_traction =
      boost::make_shared<DynamicRelaxationTimeScale>("traction_history.txt");
  for (auto &bc : *bcSpatialTractionVecPtr) {
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_traction, "traction_history", ".txt", bc.blockName);
  }
 
  auto ts_pressure =
      boost::make_shared<DynamicRelaxationTimeScale>("pressure_history.txt");
  for (auto &bc : *bcSpatialPressureVecPtr) {
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_pressure, "pressure_history", ".txt", bc.blockName);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::getExternalStrain() {
  MoFEMFunctionBegin;
 
  auto getExternalStrain = [&](boost::shared_ptr<ExternalStrainVec> &ext_strain_vec_ptr,
                             const std::string block_name,
                             const int nb_attributes) {
    MoFEMFunctionBegin;
    for (auto it : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
             std::regex((boost::format("%s(.*)") % block_name).str()))
    ) {
      std::vector<double> block_attributes;
      CHKERR it->getAttributes(block_attributes);
      if (block_attributes.size() < nb_attributes) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "In block %s expected %d attributes, but given %ld",
                 it->getName().c_str(), nb_attributes, block_attributes.size());
      }
 
      auto get_block_ents = [&]() {
        Range ents;
        CHKERR mField.get_moab().get_entities_by_handle(it->meshset, ents,
                                                        true);
        return ents;
      };
      auto Ents = get_block_ents();
      ext_strain_vec_ptr->emplace_back(it->getName(), block_attributes,
                               get_block_ents());
    }
    MoFEMFunctionReturn(0);
  };
 
  externalStrainVecPtr = boost::make_shared<ExternalStrainVec>();
  CHKERR getExternalStrain(externalStrainVecPtr, "EXTERNALSTRAIN", 2);
 
  auto ts_pre_stretch =
      boost::make_shared<DynamicRelaxationTimeScale>("externalstrain_history.txt");
  for (auto &ext_strain_block: *externalStrainVecPtr) {
    MOFEM_LOG("EP", Sev::noisy)
        << "Add time scaling external strain: " << ext_strain_block.blockName;
    timeScaleMap[ext_strain_block.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_pre_stretch, "externalstrain_history", ".txt", ext_strain_block.blockName);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::createExchangeVectors(Sev sev) {
  MoFEMFunctionBegin;
 
  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);
  };
 
  volumeExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 3, 1, sev);
  CHKERR print_loc_size(volumeExchange, "volumeExchange", sev);
  faceExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 2, 1, Sev::inform);
  CHKERR print_loc_size(faceExchange, "faceExchange", sev);
  edgeExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 1, 1, Sev::inform);
  CHKERR print_loc_size(edgeExchange, "edgeExchange", sev);
  vertexExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 0, 3, Sev::inform);
  CHKERR print_loc_size(vertexExchange, "vertexExchange", sev);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EshelbianCore::calculateCrackArea(boost::shared_ptr<double> area_ptr) {
  MoFEMFunctionBegin;
 
  if (!area_ptr) {
    CHK_THROW_MESSAGE(MOFEM_INVALID_DATA, "area_ptr is null");
  }
 
  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
 
#include <impl/EshelbianContact.cpp>