EshelbianPlasticity.cpp

Go to the documentation of this file.

/**
 * \file EshelbianPlasticity.cpp
 * \example EshelbianPlasticity.cpp
 *
 * \brief Eshelbian plasticity implementation
 */
 
#define SINGULARITY
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <CGGTonsorialBubbleBase.hpp>
 
#include <EshelbianPlasticity.hpp>
#include <boost/math/constants/constants.hpp>
 
#include <cholesky.hpp>
#ifdef PYTHON_SDF
  #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 PYTHON_SDF"
 
#else
 
  #pragma message "Without PYTHON_SDF"
 
#endif
 
#include <EshelbianAux.hpp>
#include <EshelbianContact.hpp>
 
extern "C" {
#include <phg-quadrule/quad.h>
}
 
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_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;
  if (pcomm->rank() == 0) {
    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 = subtract(crack_skin, body_skin_edges);
  }
  return send_type(m_field, crack_skin, 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::inform) << "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::inform) << "get_two_sides_of_crack_surface <- done";
 
  return std::pair<Range, Range>();
}
 
namespace EshelbianPlasticity {
 
struct DynamicRelaxationTimeScale : public TimeScale {
 
  using TimeScale::TimeScale;
 
  double getScale(const double time) override {
    if (EshelbianCore::dynamicRelaxation) {
      return TimeScale::getScale(EshelbianCore::dynamicTime);
    } else {
      return TimeScale::getScale(time);
    }
  }
};
 
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 = 4;
 
    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) {
          SETERRQ2(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 {
        SETERRQ2(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 {
 
  SetIntegrationAtFrontFace(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 = 3;
    constexpr int numEdges = 3;
    constexpr int refinementLevels = 4;
 
    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_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) {
          SETERRQ2(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 {
        SETERRQ2(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"};
  PetscInt choice_rot = EshelbianCore::rotSelector;
  PetscInt choice_grad = EshelbianCore::gradApproximator;
  PetscInt choice_symm = EshelbianCore::symmetrySelector;
 
  const char *list_stretches[] = {"linear", "log"};
  PetscInt choice_stretch = StretchSelector::LOG;
 
  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
                           "none");
  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
                         spaceOrder, &spaceOrder, PETSC_NULL);
  CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
                         spaceH1Order, &spaceH1Order, PETSC_NULL);
  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
                         "", materialOrder, &materialOrder, PETSC_NULL);
  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
                            &alphaU, PETSC_NULL);
  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
                            &alphaW, PETSC_NULL);
  CHKERR PetscOptionsScalar("-viscosity_alpha_omega", "rot viscosity", "",
                            alphaOmega, &alphaOmega, PETSC_NULL);
  CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
                            &alphaRho, PETSC_NULL);
  CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
                           LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
                           PETSC_NULL);
  CHKERR PetscOptionsEList("-grad", "gradient of defamation approximate", "",
                           list_rots, NO_H1_CONFIGURATION + 1,
                           list_rots[choice_grad], &choice_grad, PETSC_NULL);
  CHKERR PetscOptionsEList("-symm", "symmetric variant", "", list_symm, 2,
                           list_symm[choice_symm], &choice_symm, PETSC_NULL);
 
  CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
                            &EshelbianCore::exponentBase, PETSC_NULL);
  CHKERR PetscOptionsEList(
      "-stretches", "stretches", "", list_stretches, StretchSelector::LOG + 1,
      list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
 
  CHKERR PetscOptionsBool("-no_stretch", "do not solve for stretch", "",
                          noStretch, &noStretch, PETSC_NULL);
  CHKERR PetscOptionsBool("-set_singularity", "set singularity", "",
                          setSingularity, &setSingularity, PETSC_NULL);
  CHKERR PetscOptionsBool("-l2_user_base_scale", "streach scale", "",
                          l2UserBaseScale, &l2UserBaseScale, PETSC_NULL);
 
  // dynamic relaxation
  CHKERR PetscOptionsBool("-dynamic_relaxation", "dynamic time relaxation", "",
                          dynamicRelaxation, &dynamicRelaxation, PETSC_NULL);
 
  // contact parameters
  CHKERR PetscOptionsInt("-contact_max_post_proc_ref_level", "refinement level",
                         "", contactRefinementLevels, &contactRefinementLevels,
                         PETSC_NULL);
 
  // cracking parameters
  CHKERR PetscOptionsBool("-cracking_on", "cracking ON", "", crackingOn,
                          &crackingOn, PETSC_NULL);
  CHKERR PetscOptionsScalar("-griffith_energy", "Griffith energy", "",
                            griffithEnergy, &griffithEnergy, PETSC_NULL);
 
  ierr = PetscOptionsEnd();
  CHKERRG(ierr);
 
  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);
 
  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 (EshelbianCore::exponentBase != exp(1)) {
      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;
  default:
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
    break;
  };
 
  MOFEM_LOG("EP", Sev::inform) << "spaceOrder " << spaceOrder;
  MOFEM_LOG("EP", Sev::inform) << "spaceH1Order " << spaceH1Order;
  MOFEM_LOG("EP", Sev::inform) << "materialOrder " << materialOrder;
  MOFEM_LOG("EP", Sev::inform) << "alphaU " << alphaU;
  MOFEM_LOG("EP", Sev::inform) << "alphaW " << alphaW;
  MOFEM_LOG("EP", Sev::inform) << "alphaOmega " << alphaOmega;
  MOFEM_LOG("EP", Sev::inform) << "alphaRho " << alphaRho;
  MOFEM_LOG("EP", Sev::inform)
      << "Rotations " << list_rots[EshelbianCore::rotSelector];
  MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
                               << list_rots[EshelbianCore::gradApproximator];
  MOFEM_LOG("EP", Sev::inform)
      << "Symmetry " << list_symm[EshelbianCore::symmetrySelector];
  if (exponentBase != exp(1))
    MOFEM_LOG("EP", Sev::inform)
        << "Base exponent " << EshelbianCore::exponentBase;
  else
    MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
  MOFEM_LOG("EP", Sev::inform) << "Stretch " << list_stretches[choice_stretch];
 
  MOFEM_LOG("EP", Sev::inform) << "Dynamic relaxation " << (dynamicRelaxation)
      ? "yes"
      : "no";
  MOFEM_LOG("EP", Sev::inform) << "Singularity " << (setSingularity) ? "yes"
                                                                     : "no";
  MOFEM_LOG("EP", Sev::inform) << "L2 user base scale " << (l2UserBaseScale)
      ? "yes"
      : "no";
 
  MOFEM_LOG("EP", Sev::inform) << "Griffith energy " << griffithEnergy;
 
  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 (bcSpatialTraction)
      for (auto &v : *bcSpatialTraction) {
        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.vtk", *contactFaces);
  if (skeletonFaces->size())
    CHKERR save_range(mField.get_moab(), "skeleton_faces.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_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");
 
  CHKERR add_l2_field_by_range(hybridSpatialDisp, spaceOrder - 1, 2, 3,
                               Range(*skeletonFaces));
  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) {
    auto &moab = mField.get_moab();
    Range front_crack_nodes;
    CHKERR moab.get_connectivity(front_edges, front_crack_nodes, true);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        front_crack_nodes);
    Range crack_front_edges;
    CHKERR moab.get_adjacencies(front_crack_nodes, SPACE_DIM - 2, false,
                                crack_front_edges, moab::Interface::UNION);
    crack_front_edges =
        intersect(subtract(crack_front_edges, front_edges), meshset_ents);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        crack_front_edges);
 
    Range crack_front_edges_nodes;
    CHKERR moab.get_connectivity(crack_front_edges, crack_front_edges_nodes,
                                 true);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        crack_front_edges_nodes);
    crack_front_edges_nodes =
        subtract(crack_front_edges_nodes, front_crack_nodes);
    Range crack_front_edges_with_both_nodes_not_at_front;
    CHKERR moab.get_adjacencies(crack_front_edges_nodes, SPACE_DIM - 2, false,
                                crack_front_edges_with_both_nodes_not_at_front,
                                moab::Interface::UNION);
    crack_front_edges_with_both_nodes_not_at_front =
        intersect(crack_front_edges_with_both_nodes_not_at_front, meshset_ents);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        crack_front_edges_with_both_nodes_not_at_front);
    crack_front_edges_with_both_nodes_not_at_front = intersect(
        crack_front_edges_with_both_nodes_not_at_front, crack_front_edges);
 
    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;
 
  // #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_NULL, "-singularity_eps", &beta,
                                 PETSC_NULL);
    MOFEM_LOG("EP", Sev::inform) << "Singularity eps " << beta;
    eps -= beta;
 
    for (auto edge : front_adj_edges) {
      int num_nodes;
      const EntityHandle *conn;
      CHKERR moab.get_connectivity(edge, conn, num_nodes, false);
      double coords[6];
      CHKERR moab.get_coords(conn, num_nodes, coords);
      const double dir[3] = {coords[3] - coords[0], coords[4] - coords[1],
                             coords[5] - coords[2]};
      double dof[3] = {0, 0, 0};
      if (front_vertices.find(conn[0]) != front_vertices.end()) {
        for (int dd = 0; dd != 3; dd++) {
          dof[dd] = -dir[dd] * eps;
        }
      } else if (front_vertices.find(conn[1]) != front_vertices.end()) {
        for (int dd = 0; dd != 3; dd++) {
          dof[dd] = +dir[dd] * eps;
        }
      }
      for (_IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(
               mField, materialH1Positions, edge, dit)) {
        const int idx = dit->get()->getEntDofIdx();
        if (idx > 2) {
          dit->get()->getFieldData() = 0;
        } else {
          dit->get()->getFieldData() = dof[idx];
        }
      }
    }
 
    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 (bcSpatialTraction) {
      for (auto &v : *bcSpatialTraction) {
        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, hybridSpatialDisp);
    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,
                                                       hybridSpatialDisp);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       hybridMaterialDisp);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       contactDisp);
    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,
              (*bcSpatialFreeTraction)[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(3, false);
  for (int ii = 0; ii != 3; ++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();
}
 
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 (bcSpatialTraction)
    for (auto &v : *bcSpatialTraction) {
      (*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);
  }
};
 
CGGUserPolynomialBase::CGGUserPolynomialBase(
    boost::shared_ptr<CachePhi> cache_phi_otr)
    : TetPolynomialBase(), cachePhiPtr(cache_phi_otr) {}
 
MoFEMErrorCode CGGUserPolynomialBase::query_interface(
    boost::typeindex::type_index type_index,
    BaseFunctionUnknownInterface **iface) const {
  *iface = const_cast<CGGUserPolynomialBase *>(this);
  return 0;
}
 
MoFEMErrorCode
CGGUserPolynomialBase::getValue(MatrixDouble &pts,
                                boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
  MoFEMFunctionBeginHot;
 
  this->cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
 
  int nb_gauss_pts = pts.size2();
  if (!nb_gauss_pts) {
    MoFEMFunctionReturnHot(0);
  }
 
  if (pts.size1() < 3) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Wrong dimension of pts, should be at least 3 rows with "
            "coordinates");
  }
 
  const auto base = this->cTx->bAse;
  EntitiesFieldData &data = this->cTx->dAta;
 
  switch (this->cTx->sPace) {
  case HDIV:
    CHKERR getValueHdivForCGGBubble(pts);
    break;
  case L2:
    data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4, false);
    CHKERR Tools::shapeFunMBTET(
        &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
        &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
    data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
    std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
              data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
    this->cTx->basePolynomialsType0 = Legendre_polynomials;
    CHKERR getValueL2AinsworthBase(pts);
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
  }
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode
CGGUserPolynomialBase::getValueHdivForCGGBubble(MatrixDouble &pts) {
  MoFEMFunctionBegin;
 
  // This should be used only in case USER_BASE is selected
  if (this->cTx->bAse != USER_BASE) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Wrong base, should be USER_BASE");
  }
  // get access to data structures on element
  EntitiesFieldData &data = this->cTx->dAta;
  // Get approximation order on element. Note that bubble functions are only
  // on tetrahedron.
  const int order = data.dataOnEntities[MBTET][0].getOrder();
  /// number of integration points
  const int nb_gauss_pts = pts.size2();
 
  auto check_cache = [this](int order, int nb_gauss_pts) -> bool {
    if (cachePhiPtr) {
      return cachePhiPtr->get<0>() == order &&
             cachePhiPtr->get<1>() == nb_gauss_pts;
    } else {
      return false;
    }
  };
 
  if (check_cache(order, nb_gauss_pts)) {
    auto &mat = cachePhiPtr->get<2>();
    auto &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
    phi.resize(mat.size1(), mat.size2(), false);
    noalias(phi) = mat;
 
  } else {
    // calculate shape functions, i.e. barycentric coordinates
    shapeFun.resize(nb_gauss_pts, 4, false);
    CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
                      &pts(2, 0), nb_gauss_pts);
    // derivatives of shape functions
    double diff_shape_fun[12];
    CHKERR ShapeDiffMBTET(diff_shape_fun);
 
    const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
    // get base functions and set size
    MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
    phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
    // finally calculate base functions
    FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_phi(
        &phi(0, 0), &phi(0, 1), &phi(0, 2),
 
        &phi(0, 3), &phi(0, 4), &phi(0, 5),
 
        &phi(0, 6), &phi(0, 7), &phi(0, 8));
    CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
                                nb_gauss_pts);
 
    if (cachePhiPtr) {
      cachePhiPtr->get<0>() = order;
      cachePhiPtr->get<1>() = nb_gauss_pts;
      cachePhiPtr->get<2>().resize(phi.size1(), phi.size2(), false);
      noalias(cachePhiPtr->get<2>()) = phi;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
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 OpCalculateVectorFieldValuesDot<3>(
        rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradientDot<3, 3>(
        rotAxis, dataAtPts->getRotAxisGradDotAtPts(), MBTET));
 
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradientDot<6, 3>(
        stretchTensor, dataAtPts->getGradLogStretchDotTensorAtPts(), 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);
  };
 
  auto time_scale = boost::make_shared<DynamicRelaxationTimeScale>();
 
  // 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 {
      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>;
    CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
        fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {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 ? false : true) {
      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));
      pip.push_back(op_loop_domain_side);
      return std::make_pair(broken_data_ptr, piola_scale_ptr);
    };
 
    auto [broken_data_ptr, piola_scale_ptr] =
        get_broken_op_side(fe_rhs->getOpPtrVector());
 
    fe_rhs->getOpPtrVector().push_back(new OpDispBc(
        broken_data_ptr, bcSpatialDispVecPtr,
        {
 
            boost::make_shared<DynamicRelaxationTimeScale>("disp_history.txt")
 
        }));
    fe_rhs->getOpPtrVector().push_back(
        new OpRotationBc(broken_data_ptr, bcSpatialRotationVecPtr,
 
                         {
 
                             boost::make_shared<DynamicRelaxationTimeScale>(
                                 "rotation_history.txt")
 
                         }));
 
    fe_rhs->getOpPtrVector().push_back(new OpBrokenTractionBc(
        hybridSpatialDisp, bcSpatialTraction, piola_scale_ptr,
 
        {boost::make_shared<DynamicRelaxationTimeScale>("traction_history.txt")}
 
        ));
  }
 
  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"
 
struct solve_elastic_setup {
 
  inline static auto setup(EshelbianCore *ep_ptr, TS ts, Vec x,
                           bool set_ts_monitor) {
 
#ifdef PYTHON_SDF
    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_NULL, PETSC_NULL, "-sdf_file",
                                   sdf_file_name, 255, PETSC_NULL);
 
      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;
      } 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_NULL);
        auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*ep_ptr);
        ts_ctx->getLoopsMonitor().push_back(
            TsCtx::PairNameFEMethodPtr(ep_ptr->elementVolumeName, monitor_ptr));
      }
      return std::make_tuple(ts_ctx);
    };
 
    auto setup_snes_monitor = [&]() {
      MoFEMFunctionBeginHot;
      SNES snes;
      CHKERR TSGetSNES(ts, &snes);
      PetscViewerAndFormat *vf;
      CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
                                        PETSC_VIEWER_DEFAULT, &vf);
      CHKERR SNESMonitorSet(
          snes,
          (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
                             void *))SNESMonitorFields,
          vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
      MoFEMFunctionReturnHot(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);
      MoFEMFunctionReturnHot(0);
    };
 
    auto setup_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");
      }
      return schur_ptr;
    };
 
    // Warning: sequence of construction is not guaranteed for tuple. You have
    // to enforce order by proper packaging.
 
#ifdef PYTHON_SDF
    return std::make_tuple(setup_sdf(), setup_ts_monitor(),
                           setup_snes_monitor(), setup_section(),
                           set_vector_on_mesh(), setup_schur_block_solver());
#else
    return std::make_tuple(setup_ts_monitor(), setup_snes_monitor(),
                           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_NULL, "", "-debug_model", &debug_model,
                             PETSC_NULL);
 
  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);
      std::string file_skel_name =
          "snes_iteration_skel_" + std::to_string(it) + ".h5m";
      CHKERR postProcessSkeletonResults(1, file_skel_name, F);
 
      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_NULL);
  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_NULL, "", "-test_cook", &test_cook_flg,
                             PETSC_NULL);
  if (test_cook_flg) {
    constexpr int expected_lin_solver_iterations = 11;
    if (lin_solver_iterations > expected_lin_solver_iterations)
      SETERRQ2(
          PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
          "Expected number of iterations is different than expected %d > %d",
          lin_solver_iterations, expected_lin_solver_iterations);
  }
 
  CHKERR gettingNorms();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::solveDynamicRelaxation(TS ts, Vec x) {
  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;
 
  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options",
                           "none");
 
  CHKERR PetscOptionsScalar("-dynamic_final_time",
                            "dynamic relaxation final time", "", final_time,
                            &final_time, PETSC_NULL);
  CHKERR PetscOptionsScalar("-dynamic_delta_time",
                            "dynamic relaxation final time", "", delta_time,
                            &delta_time, PETSC_NULL);
  CHKERR PetscOptionsInt("-dynamic_max_it", "dynamic relaxation iterations", "",
                         max_it, &max_it, PETSC_NULL);
  CHKERR PetscOptionsBool("-dynamic_h1_update", "update each ts step", "",
                          ts_h1_update, &ts_h1_update, PETSC_NULL);
 
  ierr = PetscOptionsEnd();
  CHKERRG(ierr);
 
  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 = 0;
  dynamicStep = 0;
  monitor_ptr->ts_u = PETSC_NULL;
  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_NULL);
    if (!ts_h1_update) {
      CHKERR TSElasticPostStep::postStepFun(ts);
    }
 
    monitor_ptr->ts_u = PETSC_NULL;
    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);
  }
 
  // add tags to transfer
  for (auto t : listTagsToTransfer) {
    tags_to_transfer.push_back(t);
  }
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
  }
 
  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 (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));
      }
 
      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()));
    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) {
  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, 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));
 
  auto hybrid_res = boost::make_shared<MatrixDouble>();
  post_proc_ptr->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<3>(hybridSpatialDisp, hybrid_disp));
 
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
  post_proc_ptr->getOpPtrVector().push_back(
 
      new OpPPMap(
 
          post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
          {},
 
          {{"hybrid_disp", hybrid_disp}},
 
          {},
 
          {}
 
          )
 
  );
 
  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}},
 
            {},
 
            {}
 
            )
 
    );
  }
 
  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::calculateEnergyRelease(const int tag, TS ts) {
  MoFEMFunctionBegin;
  MOFEM_LOG_CHANNEL("SELF");
 
  constexpr bool debug = true;
 
  SmartPetscObj<IS> streach_is;
  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
      "ELASTIC_PROBLEM", ROW, stretchTensor, 0, 6, streach_is);
  SmartPetscObj<IS> piola_is;
  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
      "ELASTIC_PROBLEM", ROW, piolaStress, 0, SPACE_DIM, piola_is);
  SmartPetscObj<IS> bubble_is;
  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
      "ELASTIC_PROBLEM", ROW, bubbleField, 0, SPACE_DIM, bubble_is);
  SmartPetscObj<IS> rot_axis_is;
  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
      "ELASTIC_PROBLEM", ROW, rotAxis, 0, SPACE_DIM, rot_axis_is);
  SmartPetscObj<IS> spatial_l2_disp_is;
  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
      "ELASTIC_PROBLEM", ROW, spatialL2Disp, 0, SPACE_DIM, spatial_l2_disp_is);
  SmartPetscObj<IS> hybrid_is;
  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
      "ELASTIC_PROBLEM", ROW, hybridSpatialDisp, 0, SPACE_DIM, hybrid_is);
 
  auto get_tags_vec = [&]() {
    std::vector<Tag> tags(1);
 
    auto create_and_clean = [&]() {
      MoFEMFunctionBegin;
      auto &moab = mField.get_moab();
      auto rval = moab.tag_get_handle("ReleaseEnergy", tags[0]);
      if (rval == MB_SUCCESS) {
        moab.tag_delete(tags[0]);
      }
      double def_val[] = {0};
      CHKERR moab.tag_get_handle("ReleaseEnergy", 1, 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 tags = get_tags_vec();
 
  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);
    return adj_front;
  };
 
  auto get_nb_front_faces_on_proc = [&](auto &&adj_front) {
    int send_size = adj_front.size();
    std::vector<int> recv_data(mField.get_comm_size());
    MPI_Allgather(&send_size, 1, MPI_INT, recv_data.data(), 1, MPI_INT,
                  mField.get_comm());
    recv_data[mField.get_comm_rank()] = send_size;
    return std::pair(recv_data, adj_front);
  };
 
  auto get_snes = [&](TS ts) {
    SNES snes;
    CHKERR TSGetSNES(ts, &snes);
    return snes;
  };
 
  auto get_ksp = [&](SNES snes) {
    KSP ksp;
    CHKERR SNESGetKSP(snes, &ksp);
    CHKERR KSPSetFromOptions(ksp);
    return ksp;
  };
 
  auto integrate_energy = [&](auto x, auto release_energy_ptr) {
    MoFEMFunctionBegin;
    auto set_volume = [&]() {
      auto fe_ptr =
          boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
      fe_ptr->ts_ctx = TSMethod::CTX_TSSETIJACOBIAN;
      CHKERR setBaseVolumeElementOps(tag, false, false, false, x, fe_ptr);
      fe_ptr->getOpPtrVector().push_back(
          new OpReleaseEnergy(dataAtPts, release_energy_ptr));
      return fe_ptr;
    };
    CHKERR DMoFEMLoopFiniteElements(dM, elementVolumeName, set_volume());
    MoFEMFunctionReturn(0);
  };
 
  auto calc_release_energy = [&](auto t, auto &&tuple_of_vectors,
                                 auto adj_front, double &energy, double &area) {
    MoFEMFunctionBegin;
    MOFEM_LOG("EP", Sev::inform) << "Calculate face release energy";
 
 
    auto copy_is = [&](auto is, auto a, auto b) {
      MoFEMFunctionBegin;
      const PetscInt *is_array;
      PetscInt is_size;
      CHKERR ISGetLocalSize(is, &is_size);
      CHKERR ISGetIndices(is, &is_array);
      double *pa;
      CHKERR VecGetArray(a, &pa);
      const double *pb;
      CHKERR VecGetArrayRead(b, &pb);
      for (int i = 0; i != is_size; ++i) {
        pa[is_array[i]] = pb[is_array[i]];
      }
      CHKERR VecRestoreArrayRead(b, &pb);
      CHKERR VecRestoreArray(a, &pa);
      CHKERR ISRestoreIndices(is, &is_array);
      MoFEMFunctionReturn(0);
    };
 
    auto axpy_is = [&](auto is, double alpha, auto a, auto b) {
      MoFEMFunctionBegin;
      const PetscInt *is_array;
      PetscInt is_size;
      CHKERR ISGetLocalSize(is, &is_size);
      CHKERR ISGetIndices(is, &is_array);
      double *pa;
      CHKERR VecGetArray(a, &pa);
      const double *pb;
      CHKERR VecGetArrayRead(b, &pb);
      for (int i = 0; i != is_size; ++i) {
        pa[is_array[i]] += alpha * pb[is_array[i]];
      }
      CHKERR VecRestoreArrayRead(b, &pb);
      CHKERR VecRestoreArray(a, &pa);
      CHKERR ISRestoreIndices(is, &is_array);
      MoFEMFunctionReturn(0);
    };
 
    auto zero_is = [&](auto is, auto a) {
      MoFEMFunctionBegin;
      const PetscInt *is_array;
      PetscInt is_size;
      CHKERR ISGetLocalSize(is, &is_size);
      CHKERR ISGetIndices(is, &is_array);
      double *pa;
      CHKERR VecGetArray(a, &pa);
      for (int i = 0; i != is_size; ++i) {
        pa[is_array[i]] = 0;
      }
      CHKERR VecRestoreArray(a, &pa);
      CHKERR ISRestoreIndices(is, &is_array);
      MoFEMFunctionReturn(0);
    };
 
    auto estimate_energy = [&](auto x) {
      auto release_energy_ptr = boost::make_shared<double>(0);
      CHK_THROW_MESSAGE(integrate_energy(x, release_energy_ptr),
                        "integrate_energy");
 
      double area = 0;
      for (auto f : filter_owners(mField, adj_front)) {
        FTensor::Tensor1<double, SPACE_DIM> t_normal;
        CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
        area += t_normal.l2() / 2;
      }
 
      std::array<double, 2> send_data{area, *release_energy_ptr};
      std::array<double, 2> recv_data{0, 0};
 
      MPI_Allreduce(send_data.data(), recv_data.data(), 2, MPI_DOUBLE, MPI_SUM,
                    mField.get_comm());
      return std::pair(recv_data[1], recv_data[0]);
    };
 
    auto set_tag = [&](auto &&release_energy, auto adj_front) {
      MoFEMFunctionBegin;
      auto own = filter_owners(mField, adj_front);
      MOFEM_LOG("EP", Sev::inform) << "Energy release " << release_energy;
      auto &moab = mField.get_moab();
      CHKERR moab.tag_clear_data(tags[0], own, &release_energy);
      MoFEMFunctionReturn(0);
    };
 
    auto solve = [&]() {
      MoFEMFunctionBegin;
 
      auto [x, f] = tuple_of_vectors;
 
      Mat mA;
      auto ksp = get_ksp(get_snes(ts));
      CHKERR KSPGetOperators(ksp, &mA, PETSC_NULL);
 
      Range faces = unite(adj_front, *crackFaces);
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(faces);
 
      std::vector<boost::weak_ptr<NumeredDofEntity>> piola_skelton_dofs_vec;
      CHKERR mField.getInterface<ProblemsManager>()
          ->getSideDofsOnBrokenSpaceEntities(
              piola_skelton_dofs_vec, "ELASTIC_PROBLEM", ROW, piolaStress,
              faces, SPACE_DIM, 0, SPACE_DIM);
      SmartPetscObj<IS> skeleton_piola_is_local;
      CHKERR mField.getInterface<ISManager>()
          ->isCreateProblemBrokenFieldAndRankLocal(piola_skelton_dofs_vec,
                                                   skeleton_piola_is_local);
      SmartPetscObj<IS> skeleton_piola_is_global;
      CHKERR mField.getInterface<ISManager>()
          ->isCreateProblemBrokenFieldAndRank(piola_skelton_dofs_vec,
                                              skeleton_piola_is_global);
 
      SmartPetscObj<IS> skeleton_hybrid_is_local;
      CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
          "ELASTIC_PROBLEM", ROW, hybridSpatialDisp, 0, SPACE_DIM,
          skeleton_hybrid_is_local, &faces);
      SmartPetscObj<IS> skeleton_hybrid_is_global;
      CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
          "ELASTIC_PROBLEM", ROW, hybridSpatialDisp, 0, SPACE_DIM,
          skeleton_hybrid_is_global, &faces);
 
      CHKERR VecZeroEntries(x);
      CHKERR copy_is(skeleton_piola_is_local, x, t);
      CHKERR MatMult(mA, x, f); // f = A*x
      CHKERR zero_is(skeleton_piola_is_local, f);
      CHKERR zero_is(skeleton_hybrid_is_local, f);
 
      auto mat_storage = getBlockMatStorageMat(mA);
      auto mat_precon_storage = getBlockMatPrconditionerStorageMat(mA);
      std::vector<double> save_mat_storage(*mat_storage);
      std::vector<double> save_mat_precon_storage(*mat_precon_storage);
 
      CHKERR MatZeroRowsColumnsIS(mA, skeleton_piola_is_global, 1, PETSC_NULL,
                                  PETSC_NULL);
      CHKERR MatZeroRowsColumnsIS(mA, skeleton_hybrid_is_global, 1, PETSC_NULL,
                                  PETSC_NULL);
 
      PetscBool factor_schur = PETSC_FALSE;
      CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-crack_factor_schur",
                                 &factor_schur, PETSC_NULL);
      if (factor_schur == PETSC_TRUE) {
        SmartPetscObj<IS> schur_skeleton_hybrid_is;
        CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
            "SUB_SCHUR", ROW, hybridSpatialDisp, 0, SPACE_DIM,
            schur_skeleton_hybrid_is, &faces);
        CHKERR MatZeroEntries(S);
        CHKERR assembleBlockMatSchur(mField, mA, S, a00FieldList, a00RangeList,
                                     SmartPetscObj<AO>(aoS, true));
        // Apply essential constrains to Schur complement
        CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
        CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
        CHKERR MatZeroRowsColumnsIS(S, schur_skeleton_hybrid_is, 1, PETSC_NULL,
                                    PETSC_NULL);
      }
 
      CHKERR copy_is(skeleton_piola_is_local, f, t);
      CHKERR VecScale(f, 1. / griffithEnergy);
      CHKERR KSPSolve(ksp, f, x);
   
      CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
 
      std::copy(save_mat_storage.begin(), save_mat_storage.end(),
                mat_storage->begin());
      std::copy(save_mat_precon_storage.begin(), save_mat_precon_storage.end(),
                mat_precon_storage->begin());
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR solve();
    auto [x, f] = tuple_of_vectors;
    std::tie(energy, area) = estimate_energy(x);
    CHKERR set_tag(energy, adj_front);
 
    MoFEMFunctionReturn(0);
  };
 
  auto run_faces = [&](auto &&p) {
    auto [recv_data, adj_front] = p;
    Vec t;
    CHKERR TSGetSolution(ts, &t);
    auto face_x = vectorDuplicate(t);
    auto face_f = vectorDuplicate(t);
 
    std::map<int, std::tuple<double, double, Range, SmartPetscObj<Vec>>>
        release_by_energy;
 
    auto solve_and_add = [&](auto &&face, auto id) {
      MoFEMFunctionBegin;
      MOFEM_LOG("SYNC", Sev::noisy) << face;
      MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(face);
      double energy, area;
      CHKERR calc_release_energy(t,
                                 std::make_tuple(face_x, face_f),
                                 face, energy, area);
 
      if (release_by_energy.find(id) == release_by_energy.end()) {
        auto vec = vectorDuplicate(face_x);
        CHKERR VecCopy(face_x, vec);
        release_by_energy[id] = std::make_tuple(energy, area, face, vec);
      } else {
        auto [e, a, faces, vec] = release_by_energy.at(id);
        if (e < energy) {
          CHKERR VecCopy(face_x, vec);
          faces.merge(face);
          release_by_energy[id] = std::make_tuple(energy, area, faces, vec);
        }
      }
 
      MoFEMFunctionReturn(0);
    };
 
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
 
    for (auto p = 0; p < mField.get_comm_rank(); ++p) {
      for (auto f = 0; f < recv_data[p]; ++f) {
        int id;
        MPI_Bcast(&id, 1, MPI_INT, p, mField.get_comm());
        MOFEM_LOG_CHANNEL("SYNC");
        MOFEM_LOG("SYNC", Sev::noisy) << "edge id " << id;
        solve_and_add(Range(), id);
      }
    }
    for (auto f = 0; f < recv_data[mField.get_comm_rank()]; ++f) {
      auto face = Range(adj_front[f], adj_front[f]);
      Range edges;
      CHK_MOAB_THROW(mField.get_moab().get_adjacencies(face, 1, true, edges,
                                                       moab::Interface::UNION),
                     "get adj");
      edges = intersect(edges, *frontEdges);
      int owner;
      EntityHandle owner_handle;
      CHK_MOAB_THROW(pcomm->get_owner_handle(edges[0], owner, owner_handle),
                     "get handle");
      int id = id_from_handle(owner_handle);
      MPI_Bcast(&id, 1, MPI_INT, mField.get_comm_rank(), mField.get_comm());
      MOFEM_LOG_CHANNEL("SYNC");
      MOFEM_LOG("SYNC", Sev::noisy) << "owner " << owner << "edge id " << id;
      solve_and_add(face, id);
    }
    for (auto p = mField.get_comm_rank() + 1; p < mField.get_comm_size(); ++p) {
      for (auto f = 0; f < recv_data[p]; ++f) {
        int id;
        MPI_Bcast(&id, 1, MPI_INT, p, mField.get_comm());
        MOFEM_LOG_CHANNEL("SYNC");
        MOFEM_LOG("SYNC", Sev::noisy) << "edge id " << id;
        solve_and_add(Range(), id);
      }
    }
    MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
 
    return release_by_energy;
  };
 
  auto release_by_energy = run_faces(
 
      get_nb_front_faces_on_proc(filter_owners(mField, get_adj_front()))
 
  );
 
  std::map<double, std::tuple<double, EntityHandle, Range, SmartPetscObj<Vec>>>
      sorted_release_by_energy;
 
  for (auto &m : release_by_energy) {
    auto [energy, area, faces, vec] = m.second;
    sorted_release_by_energy[energy] =
        std::make_tuple(area, m.first, faces, vec);
  }
 
  double total_area = 0;
  double total_energy_sum = 0;
  Vec t;
  CHKERR TSGetSolution(ts, &t);
  auto x = vectorDuplicate(t);
  CHKERR VecZeroEntries(x);
  for (auto m = sorted_release_by_energy.rbegin();
       m != sorted_release_by_energy.rend(); ++m) {
    auto [area, id, faces, vec] = m->second;
    auto edge = ent_form_type_and_id(MBEDGE, id);
    CHKERR VecAXPY(x, 1., vec);
    CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
    auto release_energy_ptr = boost::make_shared<double>(0);
    CHKERR integrate_energy(x, release_energy_ptr);
    std::array<double, 2> send_data{area, *release_energy_ptr};
    std::array<double, 2> recv_data{0, 0};
    MPI_Allreduce(send_data.data(), recv_data.data(), 2, MPI_DOUBLE, MPI_SUM,
                  mField.get_comm());
    if (recv_data[1] < total_energy_sum) {
      CHKERR VecAXPY(x, -1., vec);
      double release_energy = 0;
      CHKERR mField.get_moab().tag_clear_data(tags[0], faces, &release_energy);
    } else {
      total_energy_sum = recv_data[1];
      total_area += recv_data[0];
    }
    MOFEM_LOG("EP", Sev::inform)
        << "Aggregated energy release " << total_area << " " << recv_data[1];
  }
 
  CHKERR CommInterface::updateEntitiesPetscVector(mField.get_moab(),
                                                  volumeExchange, tags[0]);
  CHKERR CommInterface::updateEntitiesPetscVector(mField.get_moab(),
                                                  faceExchange, tags[0]);
 
  if (debug) {
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
    auto skin = filter_true_skin(mField, get_skin(mField, body_ents));
    for (auto f : skin) {
      Range tet;
      CHKERR mField.get_moab().get_adjacencies(&f, 1, SPACE_DIM, false, tet);
      if (tet.size() != 1) {
        MOFEM_LOG("EP", Sev::error) << "tet.size()!=1";
        continue;
      }
      double release_energy;
      CHKERR mField.get_moab().tag_get_data(tags[0], tet, &release_energy);
      release_energy /= total_energy_sum;
      CHKERR mField.get_moab().tag_set_data(tags[0], &f, 1, &release_energy);
    }
 
    CHKERR postProcessResults(1, "test_perturbation.h5m", PETSC_NULL, x, tags);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::calculateOrientation(const int tag,
                                                   bool set_orientation) {
  MoFEMFunctionBegin;
 
  constexpr bool debug = true;
  constexpr auto sev = Sev::inform;
 
  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);
 
  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("SELF", sev) << "Nb. layers " << layers_top.size();
      int l = 0;
      for (auto &r : layers_top) {
        MOFEM_LOG("SELF", 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("SELF", 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) {
      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);
      auto l = t_dir.l2();
      t_dir(i) /= l;
      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;
      CHKERR mField.get_moab().tag_get_handle("ReleaseEnergy", th_face_energy);
 
      /**
       * Iterate over front edges, get adjacent faces, find maximal face energy.
       * Maximal face energy is stored in the edge. The is then maximises across
       * all processors.
       *
       */
      auto find_maximal_face_energy = [&](auto front_edges, auto front_faces,
                                          auto &edge_face_max_energy_map) {
        MoFEMFunctionBegin;
        for (auto e : front_edges) {
          Range faces;
          CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
          faces = intersect(faces, front_faces);
 
          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;
          CHKERR mField.get_moab().tag_set_data(th_face_energy, &e, 1,
                                                &max_energy);
          edge_face_max_energy_map[e] =
              std::make_tuple(faces[max_energy_it - face_energy.begin()],
                              max_energy, static_cast<double>(0));
          MOFEM_LOG("SELF", sev)
              << "Edge " << e << " max energy " << max_energy;
        };
        MoFEMFunctionReturn(0);
      };
 
      /**
       * For each front edge, find maximal face energy and orientation. This is
       * done by solving quadratic equation.
       *
       */
      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);
                if (adj_tets_faces.size() != 3 && adj_tets_faces.size() != 2) {
                  MOFEM_LOG("SELF", Sev::error)
                      << "Wrong number of crack faces " << adj_tets_faces.size()
                      << " " << adj_tets.size();
                  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                                    "Wrong number of crack faces");
                }
 
                FTENSOR_INDEX(SPACE_DIM, i);
 
                // cross product of face normal and edge direction
                auto t_cross_max =
                    get_cross(calculate_edge_direction(e), 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),
                                             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 energy;
                    CHKERR mField.get_moab().tag_get_data(
                        th_face_energy, adj_tets_faces, &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(energy, angle, adj_tets_faces[0]);
 
                    } else {
                      face_angle_map_down[e] =
                          std::make_tuple(energy, -angle, adj_tets_faces[0]);
                    }
                  }
                }
              }
            }
          }
 
          MoFEMFunctionReturn(0);
        };
 
        MatrixDouble A(3, 3);
        VectorDouble b(3);
 
        auto calc_optimal_angle_impl = [&](double a0, double e0, double a1,
                                           double e1, double a2, double e2) {
          A(0, 0) = a1 * a1;
          A(1, 0) = a1;
          A(2, 0) = 1;
          A(0, 1) = a2 * a2;
          A(1, 1) = a2;
          A(2, 1) = 1;
          A(0, 2) = a0 * a0;
          A(1, 2) = a0;
          A(2, 2) = 1;
 
          b[0] = e1;
          b[1] = e2;
          b[2] = e0;
 
          CHKERR solveLinearSystem(A, b);
 
          if (b[0] > -std::numeric_limits<float>::epsilon()) {
            return std::make_pair(e0, 0.);
 
          } else {
 
            auto optimal_angle = -b[1] / (2 * b[0]);
            auto optimal_energy = b[0] * optimal_angle * optimal_angle +
                                  b[1] * optimal_angle + b[2];
 
#ifndef NDEBUG
            if (debug) {
              MOFEM_LOG("SELF", sev)
                  << "Optimal_energy " << optimal_energy << " ( " << e0 << " "
                  << (optimal_energy - e0) / e0 << "% ) " << optimal_angle
                  << " e1 : a1 " << e1 << " : " << a1 << " e2 : a2 " << e2
                  << " : " << a2;
 
              double angle = -M_PI;
              for (double a = -M_PI; a < M_PI; a += M_PI / 20.) {
                double energy = b[0] * a * a + b[1] * a + b[2];
                MOFEM_LOG("SELF", sev) << "PlotAngles " << a << " " << energy;
              }
              MOFEM_LOG("SELF", sev) << "PlotAngles " << endl;
 
              MOFEM_LOG("SELF", sev) << "AnglePts " << a1 << " " << e1;
              MOFEM_LOG("SELF", sev) << "AnglePts " << a0 << " " << e0;
              MOFEM_LOG("SELF", sev) << "AnglePts " << a2 << " " << e2;
              MOFEM_LOG("SELF", sev) << "AnglePts " << endl;
            }
#endif // NDEBUG
 
            return std::make_pair(optimal_energy, optimal_angle);
          }
        };
 
        if (debug) {
#ifndef NDEBUG
          auto [opt_e0, opt_a0] = calc_optimal_angle_impl(1, 1, 0, 0, 3, -3);
          if (std::abs(opt_e0 - 1) > std::numeric_limits<double>::epsilon()) {
            CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong optimal energy");
          }
          if (std::abs(opt_a0 - 1) > std::numeric_limits<double>::epsilon()) {
            CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong optimal angle");
          }
 
#endif // NDEBUG
        }
 
        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()) {
 
              } else {
                auto [e1, a1, face_up] = face_angle_map_up.at(e);
                auto [e2, a2, face_down] = face_angle_map_down.at(e);
 
                MOFEM_LOG("SELF", sev) << "Edge " << e;
                auto [optimal_energy, optimal_angle] =
                    calc_optimal_angle_impl(a0, e0, a1, e1, a2, e2);
 
                e0 = optimal_energy;
                a0 = optimal_angle;
              }
            }
          }
 
          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);
 
        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);
          }
        }
 
        MoFEMFunctionReturn(0);
      };
 
      auto 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_face_orientation = [&](auto &sort_by_energy, auto &layers,
                                      auto &set_vertices, double &all_quality) {
        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 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 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);
          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);
        };
 
        Range rotated_normals;
 
        // iterate over energies
        for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
             ++it) {
          auto energy = it->first;
          auto [e, max_face, opt_angle] = it->second;
 
          // calculate rotation of max energy face
          auto [t_normal, t_rotated_normal] =
              get_rotated_normal(e, max_face, opt_angle);
          rotated_normals.insert(max_face);
 
          Range front_vertex; // vertices at crack front
          CHKERR mField.get_moab().get_connectivity(&e, 1, front_vertex, true);
          Range adj_tets; // adjacent tets to max face
          CHKERR mField.get_moab().get_adjacencies(&max_face, 1, 3, false,
                                                   adj_tets);
 
#ifndef NDEBUG
          if (adj_tets.size() != 2) {
            MOFEM_LOG("SELF", sev)
                << "Wrong number of tets adjacent to max face "
                << adj_tets.size();
            CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                              "Wrong number of crack faces");
          }
#endif                           // NDEBUG
          Range adj_front_faces; // adjacent faces to front edge
          CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false,
                                                   adj_front_faces);
          adj_front_faces = subtract(adj_front_faces, *crackFaces);
          Range adj_tet_faces; // adjacent faces to adjacent tets
          CHKERR mField.get_moab().get_adjacencies(
              adj_tets, 2, false, adj_tet_faces, moab::Interface::UNION);
          adj_tet_faces =
              intersect(adj_tet_faces, adj_front_faces); // only three faces
#ifndef NDEBUG
          if (adj_tet_faces.size() != 3) {
            MOFEM_LOG("SELF", sev) << "Wrong number of faces adj. to test "
                                   << adj_tet_faces.size();
            CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                              "Wrong number of crack faces");
          }
#endif // NDEBUG
 
          Range adj_front_faces_edges; // edges adjacent to front faces
          CHKERR mField.get_moab().get_adjacencies(adj_front_faces, 1, false,
                                                   adj_front_faces_edges,
                                                   moab::Interface::UNION);
 
          // only process to faces
          std::array<EntityHandle, 3> processed_faces{0, max_face, 0};
          int set_index = 0;
          for (auto &l : layers) {
            auto adj_tets_layer = intersect(adj_tets, l);
            if (adj_tets_layer.empty()) {
              continue;
            }
            Range faces_layer; // adjacent faces to adjacent tets
            CHKERR mField.get_moab().get_adjacencies(
                adj_tets_layer, 2, false, faces_layer, moab::Interface::UNION);
            faces_layer = intersect(faces_layer, adj_front_faces);
            faces_layer = subtract(faces_layer, Range(max_face, max_face));
            if (set_index > 0) {
              faces_layer = subtract(
                  faces_layer, Range(processed_faces[0], processed_faces[0]));
            }
#ifndef NDEBUG
            if (faces_layer.size() != 1) {
              MOFEM_LOG("SELF", sev)
                  << "Wrong number of faces to move " << faces_layer.size();
              CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                                "Wrong number of crack faces");
            }
#endif // NDEBUG
            processed_faces[set_index] = faces_layer[0];
            set_index = 2;
          }
 
          // choose face, if one of the vertices is already set
          for (auto f : {0, 1, 2}) {
 
            if (processed_faces[f] == 0)
              break;
 
            Range vertex; // vertex to move
            auto rval = mField.get_moab().get_connectivity(&processed_faces[f],
                                                           1, vertex, true);
            if (rval != MB_SUCCESS) {
              MOFEM_LOG("SELF", Sev::error)
                  << "get_connectivity failed " << f << " "
                  << moab::CN::EntityTypeName(type_from_handle(f));
              SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                      "can noy get connectivity");
            }
            vertex = subtract(vertex, front_vertex);
            if (set_vertices.find(vertex[0]) != set_vertices.end()) {
              if (f == 0) {
                processed_faces[1] = 0;
                processed_faces[2] = 0;
              } else if (f == 1) {
                processed_faces[0] = 0;
                processed_faces[2] = 0;
              } else {
                processed_faces[0] = 0;
                processed_faces[1] = 0;
              }
            } 
          }
 
          // [quality] -> [vertex, face, position]
          std::map<double, std::tuple<EntityHandle, EntityHandle,
                                      FTensor::Tensor1<double, SPACE_DIM>>>
              sort_by_quality;
 
          int face_idx = 0;
          for (auto f : processed_faces) {
 
            if (f == 0)
              continue;
 
            Range vertex; // vertex to move
            auto rval = mField.get_moab().get_connectivity(&f, 1, vertex, true);
            if (rval != MB_SUCCESS) {
              MOFEM_LOG("SELF", Sev::error)
                  << "get_connectivity failed " << f << " "
                  << moab::CN::EntityTypeName(type_from_handle(f));
              SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                      "can noy get connectivity");
            }
            vertex = subtract(vertex, front_vertex);
#ifndef NDEBUG
            if (vertex.size() != 1) {
              MOFEM_LOG("SELF", sev)
                  << "Wrong number of vertex to move " << vertex.size();
            }
#endif // NDEBUG
            FTensor::Tensor1<double, 3>
                t_vertex_coords; // coords of vertex to move
            CHKERR mField.get_moab().get_coords(vertex, &t_vertex_coords(0));
 
            Range vertex_edges; // edges adjacent to vertex
            CHKERR mField.get_moab().get_adjacencies(vertex, 1, false,
                                                     vertex_edges);
            vertex_edges = subtract(vertex_edges, adj_front_faces_edges);
            if (boundary_skin_verts.size() &&
                boundary_skin_verts.find(vertex[0]) !=
                    boundary_skin_verts.end()) {
              MOFEM_LOG("SELF", 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("SELF", sev) << "Geometry edge vertex";
              vertex_edges = intersect(vertex_edges, geometry_edges);
            }
            if (crack_faces_verts.size() &&
                crack_faces_verts.find(vertex_edges[0]) !=
                    crack_faces_verts.end()) {
              MOFEM_LOG("SELF", sev) << "Crack face vertex";
              vertex_edges = intersect(vertex_edges, crack_faces_edges);
            }
 
            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);
            std::map<EntityHandle, FTensor::Tensor1<double, SPACE_DIM>>
                node_positions_map;
            for (auto e : vertex_edges) {
              Range edge_conn;
              CHKERR mField.get_moab().get_connectivity(&e, 1, edge_conn, true);
              edge_conn = subtract(edge_conn, vertex);
#ifndef NDEBUG
              if (edge_conn.size() != 1) {
                MOFEM_LOG("SELF", sev)
                    << "Wrong number of edge conn " << edge_conn.size();
              }
#endif // NDEBUG
              auto it = set_vertices.find(vertex[0]);
              if (it != set_vertices.end()) {
                node_positions_map[e](i) = std::get<2>(it->second)(i);
              } else {
                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;
                MOFEM_LOG("SELF", sev)
                    << face_idx << " edge: " << e << " gamma: " << gamma;
                if (std::isnormal(gamma) &&
                    gamma > -std::numeric_limits<double>::epsilon() &&
                    gamma < 1 - std::numeric_limits<double>::epsilon()) {
                  for (auto i : {0, 1, 2}) {
                    node_positions_map[e](i) =
                        gamma * (t_v3(i) - t_vertex_coords(i));
                  }
                }
              }
            }
            if (vertex_edges.size() == 0) {
              node_positions_map[0](i) = 0;
            }
 
            Range adj_vertex_tets;
            CHKERR mField.get_moab().get_adjacencies(vertex, 3, false,
                                                     adj_vertex_tets);
            for (auto &p : node_positions_map) {
              double min_quality = 1;
              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> coords;
                CHKERR mField.get_moab().get_coords(conn, num_nodes,
                                                    coords.data());
                int idx = 0;
                for (; idx < num_nodes; ++idx) {
                  if (conn[idx] == vertex[0]) {
                    break;
                  }
                }
 
                if (set_vertices.find(vertex[0]) == set_vertices.end()) {
                  for (auto ii : {0, 1, 2}) {
                    coords[3 * idx + ii] += p.second(ii);
                  }
                }
                for (auto vv : {0, 1, 2, 3}) {
                  auto it = set_vertices.find(conn[vv]);
                  if (it != set_vertices.end()) {
                    auto &[f, energy, t_position] = it->second;
                    for (auto ii : {0, 1, 2})
                      coords[3 * vv + ii] += t_position(ii);
                  }
                }
 
                double q = Tools::volumeLengthQuality(coords.data());
                if (!std::isnormal(q))
                  q = -2;
                min_quality = std::min(min_quality, q);
              }
              sort_by_quality[min_quality] =
                  std::make_tuple(vertex[0], f,
                                  FTensor::Tensor1<double, SPACE_DIM>{
                                      p.second(0), p.second(1), p.second(2)});
            }
 
            ++face_idx;
          }
 
          if (debug) {
            for (auto s : sort_by_quality) {
              auto &[vertex, face, t_position] = s.second;
              MOFEM_LOG("SELF", sev)
                  << "Quality: " << s.first << " vertex " << vertex << " face "
                  << face << " point: " << t_position;
            }
          }
          for (auto it = sort_by_quality.rbegin(); it != sort_by_quality.rend();
               ++it) {
            auto &[vertex, face, t_position] = it->second;
            all_quality = std::min(all_quality, it->first);
            MOFEM_LOG("SELF", sev)
                << "Set quality: " << it->first << " vertex " << vertex
                << " face " << face << " point: " << t_position;
            set_vertices[vertex] = std::make_tuple(face, energy, t_position);
            break;
          }
        }
 
        MoFEMFunctionReturn(0);
      };
 
      // map: {edge, {face, energy, optimal_angle}}
      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);
 
      auto edge_map_sort_by_energy = sort_by_energy(edge_face_max_energy_map);
 
      MOFEM_LOG("SELF", sev) << "Top";
      double top_quality = 1;
      std::map<EntityHandle, std::tuple<EntityHandle, double,
                                        FTensor::Tensor1<double, SPACE_DIM>>>
          set_vertices_top;
      CHKERR set_face_orientation(edge_map_sort_by_energy, layers_top,
                                  set_vertices_top, top_quality);
      MOFEM_LOG("SELF", sev) << "Bottom";
      double bottom_quality = 1;
      std::map<EntityHandle, std::tuple<EntityHandle, double,
                                        FTensor::Tensor1<double, SPACE_DIM>>>
          set_vertices_bottom;
      CHKERR set_face_orientation(edge_map_sort_by_energy, layers_bottom,
                                  set_vertices_bottom, bottom_quality);
 
      MOFEM_LOG("SELF", sev) << "Top quality " << top_quality;
      MOFEM_LOG("SELF", sev) << "Bottom quality " << bottom_quality;
 
      auto set_tag = [&](auto &set_vertices, auto th_position) {
        MoFEMFunctionBegin;
        for (auto m : set_vertices) {
          auto v = m.first;
          auto &[f, energy, t_p] = m.second;
          CHKERR mField.get_moab().tag_set_data(th_position[0], &v, 1, &t_p(0));
          CHKERR mField.get_moab().tag_set_data(th_max_face_energy[0], &f, 1,
                                                &energy);
        }
        MoFEMFunctionReturn(0);
      };
 
      if (top_quality > bottom_quality) {
        CHKERR set_tag(set_vertices_top, th_position);
      } else {
        CHKERR set_tag(set_vertices_bottom, th_position);
      }
 
      if (debug) {
        auto th_front_top = get_tags_vec("FrontPositionTop", 3);
        auto th_front_bottom = get_tags_vec("FrontPositionBottom", 3);
 
        Range top_moved_faces;
        for (auto m : set_vertices_top) {
          auto v = m.first;
          auto &[f, energy, t_p] = m.second;
          top_moved_faces.insert(f);
          CHKERR mField.get_moab().tag_set_data(th_front_top[0], &v, 1,
                                                &t_p(0));
        }
        Range bottom_moved_faces;
        for (auto m : set_vertices_bottom) {
          auto v = m.first;
          auto &[f, energy, t_p] = m.second;
          bottom_moved_faces.insert(f);
          CHKERR mField.get_moab().tag_set_data(th_front_bottom[0], &v, 1,
                                                &t_p(0));
        }
 
        for (auto m : set_vertices_bottom) {
          auto &[f, energy, t_p] = m.second;
          bottom_moved_faces.insert(f);
        }
 
        CHKERR save_range(mField.get_moab(), "top_moved_faces.vtk",
                          top_moved_faces);
        CHKERR save_range(mField.get_moab(), "bottom_moved_faces.vtk",
                          bottom_moved_faces);
        auto front_faces = get_adj_front(false);
        CHKERR save_range(mField.get_moab(), "front_move.vtk", front_faces);
      }
 
      MoFEMFunctionReturn(0);
    };
 
    MOFEM_LOG("SELF", sev) << "Front edges " << frontEdges->size();
    CHKERR evaluate_face_energy_and_set_orientation(
        *frontEdges, get_adj_front(false), sides_pair, th_front_position);
  }
 
  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);
  };
 
  maxMovedFaces = get_max_moved_faces();
 
  if (debug) {
    Range tets;
    CHKERR mField.get_moab().get_entities_by_dimension(0, 3, tets);
    CHKERR save_range(
        mField.get_moab(),
        "projected_tets_" +
            boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
        tets);
    CHKERR save_range(
        mField.get_moab(),
        "max_moved_faces_" +
            boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
        *maxMovedFaces);
  }
 
  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);
  }
 
  constexpr bool debug = true;
  if (debug) {
 
    CHKERR save_range(
        mField.get_moab(),
        "set_coords_faces_" +
            boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
        *maxMovedFaces);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::addCrackSurfaces() {
  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);
 
        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 =
              subtract(get_skin(mField, crack_faces_org), body_skin_edges);
          auto crack_skin_nodes = get_nodes(crack_skin);
 
          auto crack_skin_faces =
              get_adj(crack_skin, 2, moab::Interface::UNION);
          crack_skin_faces = subtract(crack_skin_faces, crack_faces_org);
 
          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);
            auto edge_conn = get_nodes(Range(edges));
            if (edges.size() == 2) {
              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);
                if (edges_conn.size() == 1) {
 
                  auto node_edges = subtract(intersect(
                      get_adj(edges_conn, 1, moab::Interface::INTERSECT),
                      crack_faces_edges), crack_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 if (edges.size() == 3) {
              crack_faces.insert(f);
            }
          }
 
          crack_faces_org = crack_faces;
 
        } while (s != crack_faces.size());
      };
 
      return send_type(mField, crack_faces, MBTRI);
    };
 
    return get_faces_of_crack_front_verts(get_crack_faces());
  };
 
#ifndef NDEBUG
  constexpr bool debug = false;
  if (debug) {
    CHKERR save_range(
        mField.get_moab(),
        "new_crack_surface_" +
            boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
        get_extended_crack_faces());
  }
#endif // NDEBUG
 
  auto new_crack_faces = subtract(get_extended_crack_faces(), *crackFaces);
  CHKERR mField.getInterface<MeshsetsManager>()->addEntitiesToMeshset(
      BLOCKSET, addCrackMeshsetId, new_crack_faces);
 
  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(BLOCKSET, addCrackMeshsetId))
    ++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 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) {
 
      MOFEM_LOG("EP", Sev::noisy) << *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(bc.first, block_attributes, faces);
    }
  }
 
  // old way of naming blocksets for displacement BCs
  CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::getSpatialTractionBc() {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<ForceCubitBcData>("", piolaStress,
                                                          false, false);
 
  bcSpatialTraction = boost::make_shared<TractionBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    if (auto force_bc = bc.second->forceBcPtr) {
      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);
      bcSpatialTraction->emplace_back(bc.first, block_attributes, faces);
    }
  }
 
  CHKERR getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::createExchangeVectors(Sev sev) {
  MoFEMFunctionBegin;
 
  volumeExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 3, 1, sev);
  faceExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 2, 1, Sev::inform);
  edgeExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 1, 1, Sev::inform);
  vertexExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 0, 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(volumeExchange, "volumeExchange", sev);
  CHKERR print_loc_size(faceExchange, "faceExchange", sev);
  CHKERR print_loc_size(edgeExchange, "edgeExchange", sev);
  CHKERR print_loc_size(vertexExchange, "vertexExchange", sev);
 
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity
 
#include <impl/EshelbianContact.cpp>