EshelbianPlasticity.cpp

Go to the documentation of this file.

/**
 * \file EshelbianPlasticity.cpp
 * \example
 * mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp
 *
 * \brief Eshelbian plasticity implementation
 *
 * \copyright 2024. Various authors, some of them anonymous contributors under
 * MiT core contributors license agreement.
 */
 
#define SINGULARITY
#include <MoFEM.hpp>
 
#ifdef INCLUDE_MBCOUPLER
  #include <mbcoupler/Coupler.hpp>
#endif
using namespace MoFEM;
 
#include <CGGTonsorialBubbleBase.hpp>
 
#include <EshelbianPlasticity.hpp>
#include <boost/math/constants/constants.hpp>
 
#include <cholesky.hpp>
#ifdef ENABLE_PYTHON_BINDING
  #include <boost/python.hpp>
  #include <boost/python/def.hpp>
  #include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
 
  #pragma message "With ENABLE_PYTHON_BINDING"
 
#else
 
  #pragma message "Without ENABLE_PYTHON_BINDING"
 
#endif
 
#include <EshelbianAux.hpp>
#include <EshelbianContact.hpp>
#include <EshlabianCohesive.hpp>
#include <TSElasticPostStep.hpp>
 
extern "C" {
#include <phg-quadrule/quad.h>
}
 
#include <queue>
 
namespace EshelbianPlasticity {
struct VolUserDataOperatorStabAssembly
    : public VolumeElementForcesAndSourcesCore::UserDataOperator {
  using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
struct FaceUserDataOperatorStabAssembly
    : public FaceElementForcesAndSourcesCore::UserDataOperator {
  using FaceElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
 
} // namespace EshelbianPlasticity
 
static auto send_type(MoFEM::Interface &m_field, Range r,
                      const EntityType type) {
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
 
  auto dim = CN::Dimension(type);
 
  std::vector<int> sendcounts(pcomm->size());
  std::vector<int> displs(pcomm->size());
  std::vector<int> sendbuf(r.size());
  if (pcomm->rank() == 0) {
    for (auto p = 1; p != pcomm->size(); p++) {
      auto part_ents = m_field.getInterface<CommInterface>()
                           ->getPartEntities(m_field.get_moab(), p)
                           .subset_by_dimension(SPACE_DIM);
      Range faces;
      CHKERR m_field.get_moab().get_adjacencies(part_ents, dim, true, faces,
                                                moab::Interface::UNION);
      faces = intersect(faces, r);
      sendcounts[p] = faces.size();
      displs[p] = sendbuf.size();
      for (auto f : faces) {
        auto id = id_from_handle(f);
        sendbuf.push_back(id);
      }
    }
  }
 
  int recv_data;
  MPI_Scatter(sendcounts.data(), 1, MPI_INT, &recv_data, 1, MPI_INT, 0,
              pcomm->comm());
  std::vector<int> recvbuf(recv_data);
  MPI_Scatterv(sendbuf.data(), sendcounts.data(), displs.data(), MPI_INT,
               recvbuf.data(), recv_data, MPI_INT, 0, pcomm->comm());
 
  if (pcomm->rank() > 0) {
    Range r;
    for (auto &f : recvbuf) {
      r.insert(ent_form_type_and_id(type, f));
    }
    return r;
  }
 
  return r;
}
 
static auto get_range_from_block(MoFEM::Interface &m_field,
                                 const std::string block_name, int dim) {
  Range r;
 
  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
  auto bcs = mesh_mng->getCubitMeshsetPtr(
 
      std::regex((boost::format("%s(.*)") % block_name).str())
 
  );
 
  for (auto bc : bcs) {
    Range faces;
    CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
                                                       faces, true),
                   "get meshset ents");
    r.merge(faces);
  }
 
  return r;
};
 
static auto get_range_from_block_map(MoFEM::Interface &m_field,
                                     const std::string block_name, int dim) {
  std::map<std::string, Range> r;
 
  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
  auto bcs = mesh_mng->getCubitMeshsetPtr(
 
      std::regex((boost::format("%s(.*)") % block_name).str())
 
  );
 
  for (auto bc : bcs) {
    Range faces;
    CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
                                                       faces, true),
                   "get meshset ents");
    r[bc->getName()] = faces;
  }
 
  return r;
}
 
static auto get_block_meshset(MoFEM::Interface &m_field, const int ms_id,
                              const unsigned int cubit_bc_type) {
  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
  EntityHandle meshset;
  CHKERR mesh_mng->getMeshset(ms_id, cubit_bc_type, meshset);
  return meshset;
};
 
static auto save_range(moab::Interface &moab, const std::string name,
                       const Range r, std::vector<Tag> tags = {}) {
  MoFEMFunctionBegin;
  auto out_meshset = get_temp_meshset_ptr(moab);
  CHKERR moab.add_entities(*out_meshset, r);
  if (r.size()) {
    CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1,
                           tags.data(), tags.size());
  } else {
    MOFEM_LOG("SELF", Sev::warning) << "Empty range for " << name;
  }
  MoFEMFunctionReturn(0);
};
 
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin) {
  Range boundary_ents;
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
  CHK_MOAB_THROW(pcomm->filter_pstatus(skin,
                                       PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                       PSTATUS_NOT, -1, &boundary_ents),
                 "filter_pstatus");
  return boundary_ents;
};
 
static auto filter_owners(MoFEM::Interface &m_field, Range skin) {
  Range owner_ents;
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
  CHK_MOAB_THROW(pcomm->filter_pstatus(skin, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1,
                                       &owner_ents),
                 "filter_pstatus");
  return owner_ents;
};
 
static auto get_skin(MoFEM::Interface &m_field, Range body_ents) {
  Skinner skin(&m_field.get_moab());
  Range skin_ents;
  CHK_MOAB_THROW(skin.find_skin(0, body_ents, false, skin_ents), "find_skin");
  return skin_ents;
};
 
static auto get_crack_front_edges(MoFEM::Interface &m_field,
                                  Range crack_faces) {
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
  auto &moab = m_field.get_moab();
  Range crack_skin_without_bdy;
  if (pcomm->rank() == 0) {
    Range crack_edges;
    CHKERR moab.get_adjacencies(crack_faces, 1, true, crack_edges,
                                moab::Interface::UNION);
    auto crack_skin = get_skin(m_field, crack_faces);
    Range body_ents;
    CHK_MOAB_THROW(
        m_field.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents),
        "get_entities_by_dimension");
    auto body_skin = get_skin(m_field, body_ents);
    Range body_skin_edges;
    CHK_MOAB_THROW(moab.get_adjacencies(body_skin, 1, true, body_skin_edges,
                                        moab::Interface::UNION),
                   "get_adjacencies");
    crack_skin_without_bdy = subtract(crack_skin, body_skin_edges);
    auto front_edges_map = get_range_from_block_map(m_field, "FRONT", 1);
    for (auto &m : front_edges_map) {
      auto add_front = subtract(m.second, crack_edges);
      auto i = intersect(m.second, crack_edges);
      if (i.empty()) {
        crack_skin_without_bdy.merge(add_front);
      } else {
        auto i_skin = get_skin(m_field, i);
        Range adj_i_skin;
        CHKERR moab.get_adjacencies(i_skin, 1, true, adj_i_skin,
                                    moab::Interface::UNION);
        adj_i_skin = subtract(intersect(adj_i_skin, m.second), crack_edges);
        crack_skin_without_bdy.merge(adj_i_skin);
      }
    }
  }
  return send_type(m_field, crack_skin_without_bdy, MBEDGE);
}
 
static auto get_two_sides_of_crack_surface(MoFEM::Interface &m_field,
                                           Range crack_faces) {
 
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
 
  MOFEM_LOG("EP", Sev::noisy) << "get_two_sides_of_crack_surface";
 
  if (!pcomm->rank()) {
 
    auto impl = [&](auto &saids) {
      MoFEMFunctionBegin;
 
      auto &moab = m_field.get_moab();
 
      auto get_adj = [&](auto e, auto dim) {
        Range adj;
        CHK_MOAB_THROW(m_field.get_moab().get_adjacencies(
                           e, dim, true, adj, moab::Interface::UNION),
                       "get adj");
        return adj;
      };
 
      auto get_conn = [&](auto e) {
        Range conn;
        CHK_MOAB_THROW(m_field.get_moab().get_connectivity(e, conn, true),
                       "get connectivity");
        return conn;
      };
 
      constexpr bool debug = false;
      Range body_ents;
      CHKERR m_field.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                          body_ents);
      auto body_skin = get_skin(m_field, body_ents);
      auto body_skin_edges = get_adj(body_skin, 1);
 
      auto crack_skin =
          subtract(get_skin(m_field, crack_faces), body_skin_edges);
      auto crack_skin_conn = get_conn(crack_skin);
      auto crack_skin_conn_edges = get_adj(crack_skin_conn, 1);
      auto crack_edges = get_adj(crack_faces, 1);
      crack_edges = subtract(crack_edges, crack_skin);
      auto all_tets = get_adj(crack_edges, 3);
      crack_edges = subtract(crack_edges, crack_skin_conn_edges);
      auto crack_conn = get_conn(crack_edges);
      all_tets.merge(get_adj(crack_conn, 3));
 
      if (debug) {
        CHKERR save_range(m_field.get_moab(), "crack_faces.vtk", crack_faces);
        CHKERR save_range(m_field.get_moab(), "all_crack_tets.vtk", all_tets);
        CHKERR save_range(m_field.get_moab(), "crack_edges_all.vtk",
                          crack_edges);
      }
 
      if (crack_faces.size()) {
        auto grow = [&](auto r) {
          auto crack_faces_conn = get_conn(crack_faces);
          Range v;
          auto size_r = 0;
          while (size_r != r.size() && r.size() > 0) {
            size_r = r.size();
            CHKERR moab.get_connectivity(r, v, true);
            v = subtract(v, crack_faces_conn);
            if (v.size()) {
              CHKERR moab.get_adjacencies(v, SPACE_DIM, true, r,
                                          moab::Interface::UNION);
              r = intersect(r, all_tets);
            }
            if (r.empty()) {
              break;
            }
          }
          return r;
        };
 
        Range all_tets_ord = all_tets;
        while (all_tets.size()) {
          Range faces = get_adj(unite(saids.first, saids.second), 2);
          faces = subtract(crack_faces, faces);
          if (faces.size()) {
            Range tets;
            auto fit = faces.begin();
            for (; fit != faces.end(); ++fit) {
              tets = intersect(get_adj(Range(*fit, *fit), 3), all_tets);
              if (tets.size() == 2) {
                break;
              }
            }
            if (tets.empty()) {
              break;
            } else {
              saids.first.insert(tets[0]);
              saids.first = grow(saids.first);
              all_tets = subtract(all_tets, saids.first);
              if (tets.size() == 2) {
                saids.second.insert(tets[1]);
                saids.second = grow(saids.second);
                all_tets = subtract(all_tets, saids.second);
              }
            }
          } else {
            break;
          }
        }
 
        saids.first = subtract(all_tets_ord, saids.second);
        saids.second = subtract(all_tets_ord, saids.first);
      }
 
      MoFEMFunctionReturn(0);
    };
 
    std::pair<Range, Range> saids;
    if (crack_faces.size())
      CHK_THROW_MESSAGE(impl(saids), "get crack both sides");
    return saids;
  }
 
  MOFEM_LOG("EP", Sev::noisy) << "get_two_sides_of_crack_surface <- done";
 
  return std::pair<Range, Range>();
}
 
namespace EshelbianPlasticity {
 
struct SetIntegrationAtFrontVolume {
 
  using FunRule = boost::function<int(int)>;
  FunRule funRule = [](int p) { return 2 * p + 1; };
 
  SetIntegrationAtFrontVolume(
      boost::shared_ptr<Range> front_nodes,
      boost::shared_ptr<Range> front_edges,
      boost::shared_ptr<CGGUserPolynomialBase::CachePhi> cache_phi = nullptr)
      : frontNodes(front_nodes), frontEdges(front_edges), cachePhi(cache_phi){};
 
  SetIntegrationAtFrontVolume(
      boost::shared_ptr<Range> front_nodes,
      boost::shared_ptr<Range> front_edges, FunRule fun_rule,
      boost::shared_ptr<CGGUserPolynomialBase::CachePhi> cache_phi = nullptr)
      : frontNodes(front_nodes), frontEdges(front_edges), funRule(fun_rule),
        cachePhi(cache_phi){};
 
  MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                            int order_col, int order_data) {
    MoFEMFunctionBegin;
 
    constexpr bool debug = false;
 
    constexpr int numNodes = 4;
    constexpr int numEdges = 6;
    constexpr int refinementLevels = 6;
 
    auto &m_field = fe_raw_ptr->mField;
    auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
    auto fe_handle = fe_ptr->getFEEntityHandle();
 
    auto set_base_quadrature = [&]() {
      MoFEMFunctionBegin;
      int rule = funRule(order_data);
      if (rule < QUAD_3D_TABLE_SIZE) {
        if (QUAD_3D_TABLE[rule]->dim != 3) {
          SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
                  "wrong dimension");
        }
        if (QUAD_3D_TABLE[rule]->order < rule) {
          SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
                  "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
        }
        const size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
        auto &gauss_pts = fe_ptr->gaussPts;
        gauss_pts.resize(4, nb_gauss_pts, false);
        cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
                    &gauss_pts(0, 0), 1);
        cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
                    &gauss_pts(1, 0), 1);
        cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
                    &gauss_pts(2, 0), 1);
        cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
                    &gauss_pts(3, 0), 1);
        auto &data = fe_ptr->dataOnElement[H1];
        data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
                                                              false);
        double *shape_ptr =
            &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
        cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
                    1);
      } else {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_base_quadrature();
 
    if (frontNodes->size() && EshelbianCore::setSingularity) {
 
      auto get_singular_nodes = [&]() {
        int num_nodes;
        const EntityHandle *conn;
        CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
                                                   true);
        std::bitset<numNodes> singular_nodes;
        for (auto nn = 0; nn != numNodes; ++nn) {
          if (frontNodes->find(conn[nn]) != frontNodes->end()) {
            singular_nodes.set(nn);
          } else {
            singular_nodes.reset(nn);
          }
        }
        return singular_nodes;
      };
 
      auto get_singular_edges = [&]() {
        std::bitset<numEdges> singular_edges;
        for (int ee = 0; ee != numEdges; ee++) {
          EntityHandle edge;
          CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
          if (frontEdges->find(edge) != frontEdges->end()) {
            singular_edges.set(ee);
          } else {
            singular_edges.reset(ee);
          }
        }
        return singular_edges;
      };
 
      auto set_gauss_pts = [&](auto &ref_gauss_pts) {
        MoFEMFunctionBegin;
        fe_ptr->gaussPts.swap(ref_gauss_pts);
        const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
        auto &data = fe_ptr->dataOnElement[H1];
        data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
        double *shape_ptr =
            &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
        CHKERR ShapeMBTET(shape_ptr, &fe_ptr->gaussPts(0, 0),
                          &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
                          nb_gauss_pts);
        MoFEMFunctionReturn(0);
      };
 
      auto singular_nodes = get_singular_nodes();
      if (singular_nodes.count()) {
        auto it_map_ref_coords = mapRefCoords.find(singular_nodes.to_ulong());
        if (it_map_ref_coords != mapRefCoords.end()) {
          CHKERR set_gauss_pts(it_map_ref_coords->second);
          MoFEMFunctionReturnHot(0);
        } else {
 
          auto refine_quadrature = [&]() {
            MoFEMFunctionBegin;
 
            const int max_level = refinementLevels;
            EntityHandle tet;
 
            moab::Core moab_ref;
            double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
            EntityHandle nodes[4];
            for (int nn = 0; nn != 4; nn++)
              CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
            CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
            MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
            MoFEM::Interface &m_field_ref = mofem_ref_core;
            {
              Range tets(tet, tet);
              Range edges;
              CHKERR m_field_ref.get_moab().get_adjacencies(
                  tets, 1, true, edges, moab::Interface::UNION);
              CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
                  tets, BitRefLevel().set(0), false, VERBOSE);
            }
 
            Range nodes_at_front;
            for (int nn = 0; nn != numNodes; nn++) {
              if (singular_nodes[nn]) {
                EntityHandle ent;
                CHKERR moab_ref.side_element(tet, 0, nn, ent);
                nodes_at_front.insert(ent);
              }
            }
 
            auto singular_edges = get_singular_edges();
 
            EntityHandle meshset;
            CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
            for (int ee = 0; ee != numEdges; ee++) {
              if (singular_edges[ee]) {
                EntityHandle ent;
                CHKERR moab_ref.side_element(tet, 1, ee, ent);
                CHKERR moab_ref.add_entities(meshset, &ent, 1);
              }
            }
 
            // refine mesh
            auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
            for (int ll = 0; ll != max_level; ll++) {
              Range edges;
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                                 BitRefLevel().set(), MBEDGE,
                                                 edges);
              Range ref_edges;
              CHKERR moab_ref.get_adjacencies(
                  nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
              ref_edges = intersect(ref_edges, edges);
              Range ents;
              CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
              ref_edges = intersect(ref_edges, ents);
              Range tets;
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->getEntitiesByTypeAndRefLevel(
                      BitRefLevel().set(ll), BitRefLevel().set(), MBTET, tets);
              CHKERR m_ref->addVerticesInTheMiddleOfEdges(
                  ref_edges, BitRefLevel().set(ll + 1));
              CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->updateMeshsetByEntitiesChildren(meshset,
                                                    BitRefLevel().set(ll + 1),
                                                    meshset, MBEDGE, true);
            }
 
            // get ref coords
            Range tets;
            CHKERR m_field_ref.getInterface<BitRefManager>()
                ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
                                               BitRefLevel().set(), MBTET,
                                               tets);
 
            if (debug) {
              CHKERR save_range(moab_ref, "ref_tets.vtk", tets);
            }
 
            MatrixDouble ref_coords(tets.size(), 12, false);
            int tt = 0;
            for (Range::iterator tit = tets.begin(); tit != tets.end();
                 tit++, tt++) {
              int num_nodes;
              const EntityHandle *conn;
              CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
              CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
            }
 
            auto &data = fe_ptr->dataOnElement[H1];
            const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
            MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
            MatrixDouble &shape_n =
                data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
            int gg = 0;
            for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
              double *tet_coords = &ref_coords(tt, 0);
              double det = Tools::tetVolume(tet_coords);
              det *= 6;
              for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
                for (int dd = 0; dd != 3; dd++) {
                  ref_gauss_pts(dd, gg) =
                      shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
                      shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
                      shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
                      shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
                }
                ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
              }
            }
 
            mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
            CHKERR set_gauss_pts(mapRefCoords[singular_nodes.to_ulong()]);
 
            // clear cache bubble
            cachePhi->get<0>() = 0;
            cachePhi->get<1>() = 0;
            // tet base cache
            TetPolynomialBase::switchCacheBaseOff<HDIV>({fe_raw_ptr});
            TetPolynomialBase::switchCacheBaseOn<HDIV>({fe_raw_ptr});
 
            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;
 
  boost::shared_ptr<CGGUserPolynomialBase::CachePhi> cachePhi;
 
  static inline std::map<long int, MatrixDouble> mapRefCoords;
};
 
struct SetIntegrationAtFrontFace {
 
  using FunRule = boost::function<int(int)>;
  FunRule funRule = [](int p) { return 2 * (p + 1) + 1; };
 
  SetIntegrationAtFrontFace(boost::shared_ptr<Range> front_nodes,
                            boost::shared_ptr<Range> front_edges)
      : frontNodes(front_nodes), frontEdges(front_edges){};
 
  SetIntegrationAtFrontFace(boost::shared_ptr<Range> front_nodes,
                            boost::shared_ptr<Range> front_edges,
                            FunRule fun_rule)
      : frontNodes(front_nodes), frontEdges(front_edges), funRule(fun_rule){};
 
  MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                            int order_col, int order_data) {
    MoFEMFunctionBegin;
 
    constexpr bool debug = false;
 
    constexpr int numNodes = 3;
    constexpr int numEdges = 3;
    constexpr int refinementLevels = 6;
 
    auto &m_field = fe_raw_ptr->mField;
    auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
    auto fe_handle = fe_ptr->getFEEntityHandle();
 
    auto set_base_quadrature = [&]() {
      MoFEMFunctionBegin;
      int rule = funRule(order_data);
      if (rule < QUAD_2D_TABLE_SIZE) {
        if (QUAD_2D_TABLE[rule]->dim != 2) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
        }
        if (QUAD_2D_TABLE[rule]->order < rule) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
        }
        const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
        fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
        cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
                    &fe_ptr->gaussPts(0, 0), 1);
        cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
                    &fe_ptr->gaussPts(1, 0), 1);
        cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
                    &fe_ptr->gaussPts(2, 0), 1);
        auto &data = fe_ptr->dataOnElement[H1];
        data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
                                                              false);
        double *shape_ptr =
            &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
        cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
                    1);
        data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
        std::copy(
            Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
            data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
 
      } else {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_base_quadrature();
 
    if (frontNodes->size() && EshelbianCore::setSingularity) {
 
      auto get_singular_nodes = [&]() {
        int num_nodes;
        const EntityHandle *conn;
        CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
                                                   true);
        std::bitset<numNodes> singular_nodes;
        for (auto nn = 0; nn != numNodes; ++nn) {
          if (frontNodes->find(conn[nn]) != frontNodes->end()) {
            singular_nodes.set(nn);
          } else {
            singular_nodes.reset(nn);
          }
        }
        return singular_nodes;
      };
 
      auto get_singular_edges = [&]() {
        std::bitset<numEdges> singular_edges;
        for (int ee = 0; ee != numEdges; ee++) {
          EntityHandle edge;
          CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
          if (frontEdges->find(edge) != frontEdges->end()) {
            singular_edges.set(ee);
          } else {
            singular_edges.reset(ee);
          }
        }
        return singular_edges;
      };
 
      auto set_gauss_pts = [&](auto &ref_gauss_pts) {
        MoFEMFunctionBegin;
        fe_ptr->gaussPts.swap(ref_gauss_pts);
        const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
        auto &data = fe_ptr->dataOnElement[H1];
        data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
        double *shape_ptr =
            &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
        CHKERR ShapeMBTRI(shape_ptr, &fe_ptr->gaussPts(0, 0),
                          &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
        MoFEMFunctionReturn(0);
      };
 
      auto singular_nodes = get_singular_nodes();
      if (singular_nodes.count()) {
        auto it_map_ref_coords = mapRefCoords.find(singular_nodes.to_ulong());
        if (it_map_ref_coords != mapRefCoords.end()) {
          CHKERR set_gauss_pts(it_map_ref_coords->second);
          MoFEMFunctionReturnHot(0);
        } else {
 
          auto refine_quadrature = [&]() {
            MoFEMFunctionBegin;
 
            const int max_level = refinementLevels;
 
            moab::Core moab_ref;
            double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
            EntityHandle nodes[numNodes];
            for (int nn = 0; nn != numNodes; nn++)
              CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
            EntityHandle tri;
            CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
            MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
            MoFEM::Interface &m_field_ref = mofem_ref_core;
            {
              Range tris(tri, tri);
              Range edges;
              CHKERR m_field_ref.get_moab().get_adjacencies(
                  tris, 1, true, edges, moab::Interface::UNION);
              CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
                  tris, BitRefLevel().set(0), false, VERBOSE);
            }
 
            Range nodes_at_front;
            for (int nn = 0; nn != numNodes; nn++) {
              if (singular_nodes[nn]) {
                EntityHandle ent;
                CHKERR moab_ref.side_element(tri, 0, nn, ent);
                nodes_at_front.insert(ent);
              }
            }
 
            auto singular_edges = get_singular_edges();
 
            EntityHandle meshset;
            CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
            for (int ee = 0; ee != numEdges; ee++) {
              if (singular_edges[ee]) {
                EntityHandle ent;
                CHKERR moab_ref.side_element(tri, 1, ee, ent);
                CHKERR moab_ref.add_entities(meshset, &ent, 1);
              }
            }
 
            // refine mesh
            auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
            for (int ll = 0; ll != max_level; ll++) {
              Range edges;
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                                 BitRefLevel().set(), MBEDGE,
                                                 edges);
              Range ref_edges;
              CHKERR moab_ref.get_adjacencies(
                  nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
              ref_edges = intersect(ref_edges, edges);
              Range ents;
              CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
              ref_edges = intersect(ref_edges, ents);
              Range tris;
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->getEntitiesByTypeAndRefLevel(
                      BitRefLevel().set(ll), BitRefLevel().set(), MBTRI, tris);
              CHKERR m_ref->addVerticesInTheMiddleOfEdges(
                  ref_edges, BitRefLevel().set(ll + 1));
              CHKERR m_ref->refineTris(tris, BitRefLevel().set(ll + 1));
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->updateMeshsetByEntitiesChildren(meshset,
                                                    BitRefLevel().set(ll + 1),
                                                    meshset, MBEDGE, true);
            }
 
            // get ref coords
            Range tris;
            CHKERR m_field_ref.getInterface<BitRefManager>()
                ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
                                               BitRefLevel().set(), MBTRI,
                                               tris);
 
            if (debug) {
              CHKERR save_range(moab_ref, "ref_tris.vtk", tris);
              SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "debug");
            }
 
            MatrixDouble ref_coords(tris.size(), 9, false);
            int tt = 0;
            for (Range::iterator tit = tris.begin(); tit != tris.end();
                 tit++, tt++) {
              int num_nodes;
              const EntityHandle *conn;
              CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
              CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
            }
 
            auto &data = fe_ptr->dataOnElement[H1];
            const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
            MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
            MatrixDouble &shape_n =
                data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
            int gg = 0;
            for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
              double *tri_coords = &ref_coords(tt, 0);
              FTensor::Tensor1<double, 3> t_normal;
              CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
              auto det = t_normal.l2();
              for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
                for (int dd = 0; dd != 2; dd++) {
                  ref_gauss_pts(dd, gg) =
                      shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
                      shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
                      shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
                }
                ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
              }
            }
 
            mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
            CHKERR set_gauss_pts(mapRefCoords[singular_nodes.to_ulong()]);
 
            MoFEMFunctionReturn(0);
          };
 
          CHKERR refine_quadrature();
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  struct Fe : public ForcesAndSourcesCore {
    using ForcesAndSourcesCore::dataOnElement;
 
  private:
    using ForcesAndSourcesCore::ForcesAndSourcesCore;
  };
 
  boost::shared_ptr<Range> frontNodes;
  boost::shared_ptr<Range> frontEdges;
 
  static inline std::map<long int, MatrixDouble> mapRefCoords;
};
 
double EshelbianCore::exponentBase = exp(1);
boost::function<double(const double)> EshelbianCore::f = EshelbianCore::f_log_e;
boost::function<double(const double)> EshelbianCore::d_f =
    EshelbianCore::d_f_log_e;
boost::function<double(const double)> EshelbianCore::dd_f =
    EshelbianCore::dd_f_log_e;
boost::function<double(const double)> EshelbianCore::inv_f =
    EshelbianCore::inv_f_log_e;
boost::function<double(const double)> EshelbianCore::inv_d_f =
    EshelbianCore::inv_d_f_log_e;
boost::function<double(const double)> EshelbianCore::inv_dd_f =
    EshelbianCore::inv_dd_f_log_e;
 
MoFEMErrorCode
EshelbianCore::query_interface(boost::typeindex::type_index type_index,
                               UnknownInterface **iface) const {
  *iface = const_cast<EshelbianCore *>(this);
  return 0;
}
 
MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
 
  if (evalRhs)
    CHKERR evaluateRhs(data);
 
  if (evalLhs)
    CHKERR evaluateLhs(data);
 
  MoFEMFunctionReturn(0);
}
 
EshelbianCore::EshelbianCore(MoFEM::Interface &m_field) : mField(m_field) {
  CHK_THROW_MESSAGE(getOptions(), "getOptions failed");
}
 
EshelbianCore::~EshelbianCore() = default;
 
MoFEMErrorCode EshelbianCore::getOptions() {
  MoFEMFunctionBegin;
  const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
  const char *list_symm[] = {"symm", "not_symm"};
  const char *list_release[] = {"griffith_force", "griffith_skeleton"};
  const char *list_stretches[] = {"linear", "log", "log_quadratic"};
  const char *list_solvers[] = {"time_solver", "dynamic_relaxation",
                                "cohesive"};
  PetscInt choice_rot = EshelbianCore::rotSelector;
  PetscInt choice_grad = EshelbianCore::gradApproximator;
  PetscInt choice_symm = EshelbianCore::symmetrySelector;
  PetscInt choice_release = EshelbianCore::energyReleaseSelector;
  PetscInt choice_stretch = StretchSelector::LOG;
  PetscInt choice_solver = SolverType::TimeSolver;
  char analytical_expr_file_name[255] = "analytical_expr.py";
 
  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity", "none");
  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
                         spaceOrder, &spaceOrder, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
                         spaceH1Order, &spaceH1Order, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
                         "", materialOrder, &materialOrder, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
                            &alphaU, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
                            &alphaW, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_omega", "rot viscosity", "",
                            alphaOmega, &alphaOmega, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
                            &alphaRho, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-alpha_tau", "tau", "", alphaTau, &alphaTau,
                            PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
                           LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
                           PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-grad", "gradient of defamation approximate", "",
                           list_rots, NO_H1_CONFIGURATION + 1,
                           list_rots[choice_grad], &choice_grad, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-symm", "symmetric variant", "", list_symm, 2,
                           list_symm[choice_symm], &choice_symm, PETSC_NULLPTR);
 
  CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
                            &EshelbianCore::exponentBase, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-stretches", "stretches", "", list_stretches,
                           StretchSelector::STRETCH_SELECTOR_LAST,
                           list_stretches[choice_stretch], &choice_stretch,
                           PETSC_NULLPTR);
 
  CHKERR PetscOptionsBool("-no_stretch", "do not solve for stretch", "",
                          noStretch, &noStretch, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-set_singularity", "set singularity", "",
                          setSingularity, &setSingularity, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-l2_user_base_scale", "streach scale", "",
                          l2UserBaseScale, &l2UserBaseScale, PETSC_NULLPTR);
 
  // dynamic relaxation
 
  // @deprecate this option
  CHKERR PetscOptionsBool("-dynamic_relaxation", "dynamic time relaxation", "",
                          dynamicRelaxation, &dynamicRelaxation, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-solver_type", "solver type", "", list_solvers,
                           SolverType::LastSolver, list_solvers[choice_solver],
                           &choice_solver, PETSC_NULLPTR);
 
  // contact parameters
  CHKERR PetscOptionsInt("-contact_max_post_proc_ref_level", "refinement level",
                         "", contactRefinementLevels, &contactRefinementLevels,
                         PETSC_NULLPTR);
  // cohesive interfae
  CHKERR PetscOptionsBool("-cohesive_interface_on", "cohesive interface ON", "",
                          intefaceCrack, &intefaceCrack, PETSC_NULLPTR);
 
  // cracking parameters
  CHKERR PetscOptionsBool("-cracking_on", "cracking ON", "", crackingOn,
                          &crackingOn, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-cracking_start_time", "cracking start time", "",
                            crackingStartTime, &crackingStartTime,
                            PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-griffith_energy", "Griffith energy", "",
                            griffithEnergy, &griffithEnergy, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-energy_release_variant", "energy release variant",
                           "", list_release, 2, list_release[choice_release],
                           &choice_release, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-nb_J_integral_levels", "Number of J integarl levels",
                         "", nbJIntegralContours, &nbJIntegralContours,
                         PETSC_NULLPTR); // backward compatibility
  CHKERR PetscOptionsInt(
      "-nb_J_integral_contours", "Number of J integral contours", "",
      nbJIntegralContours, &nbJIntegralContours, PETSC_NULLPTR);
 
  // internal stress
  char tag_name[255] = "";
  CHKERR PetscOptionsString("-internal_stress_tag_name",
                            "internal stress tag name", "", "", tag_name, 255,
                            PETSC_NULLPTR);
  internalStressTagName = string(tag_name);
  CHKERR PetscOptionsInt(
      "-internal_stress_interp_order", "internal stress interpolation order",
      "", internalStressInterpOrder, &internalStressInterpOrder, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-internal_stress_voigt", "Voigt index notation", "",
                          internalStressVoigt, &internalStressVoigt,
                          PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR,
                               "-analytical_expr_file",
                               analytical_expr_file_name, 255, PETSC_NULLPTR);
 
  PetscOptionsEnd();
 
  if (internalStressInterpOrder != 0 && internalStressInterpOrder != 1) {
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
            "Unsupported internal stress interpolation order %d",
            internalStressInterpOrder);
  }
 
  if (setSingularity) {
    l2UserBaseScale = PETSC_TRUE;
  }
 
  EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
  EshelbianCore::gradApproximator = static_cast<RotSelector>(choice_grad);
  EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
  EshelbianCore::symmetrySelector = static_cast<SymmetrySelector>(choice_symm);
  EshelbianCore::energyReleaseSelector =
      static_cast<EnergyReleaseSelector>(choice_release);
 
  switch (EshelbianCore::stretchSelector) {
  case StretchSelector::LINEAR:
    EshelbianCore::f = f_linear;
    EshelbianCore::d_f = d_f_linear;
    EshelbianCore::dd_f = dd_f_linear;
    EshelbianCore::inv_f = inv_f_linear;
    EshelbianCore::inv_d_f = inv_d_f_linear;
    EshelbianCore::inv_dd_f = inv_dd_f_linear;
    break;
  case StretchSelector::LOG:
    if (std::fabs(EshelbianCore::exponentBase - exp(1)) >
        std::numeric_limits<float>::epsilon()) {
      EshelbianCore::f = f_log;
      EshelbianCore::d_f = d_f_log;
      EshelbianCore::dd_f = dd_f_log;
      EshelbianCore::inv_f = inv_f_log;
      EshelbianCore::inv_d_f = inv_d_f_log;
      EshelbianCore::inv_dd_f = inv_dd_f_log;
    } else {
      EshelbianCore::f = f_log_e;
      EshelbianCore::d_f = d_f_log_e;
      EshelbianCore::dd_f = dd_f_log_e;
      EshelbianCore::inv_f = inv_f_log;
      EshelbianCore::inv_d_f = inv_d_f_log;
      EshelbianCore::inv_dd_f = inv_dd_f_log;
    }
    break;
  case StretchSelector::LOG_QUADRATIC:
    EshelbianCore::f = f_log_e_quadratic;
    EshelbianCore::d_f = d_f_log_e_quadratic;
    EshelbianCore::dd_f = dd_f_log_e_quadratic;
    EshelbianCore::inv_f = [](const double x) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "No logarithmic quadratic stretch for this case");
      return 0;
    };
    EshelbianCore::inv_d_f = EshelbianCore::inv_f;
    EshelbianCore::inv_dd_f = EshelbianCore::inv_f;
    break; // no stretch, do not use stretch functions
  default:
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
    break;
  };
 
  if (dynamicRelaxation) {
    MOFEM_LOG("EP", Sev::warning)
        << "-dynamic_relaxation option is deprecated, use -solver_type "
           "dynamic_relaxation instead.";
    choice_solver = SolverType::DynamicRelaxation;
  }
 
  switch (choice_solver) {
  case SolverType::TimeSolver:
    EshelbianCore::solverType = static_cast<SolverType>(SolverType::TimeSolver);
    break;
  case SolverType::DynamicRelaxation:
    EshelbianCore::solverType =
        static_cast<SolverType>(SolverType::DynamicRelaxation);
    dynamicRelaxation = PETSC_TRUE;
    break;
  case SolverType::Cohesive:
    EshelbianCore::solverType = static_cast<SolverType>(SolverType::Cohesive);
    break;
  default:
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown solver");
    break;
  };
 
  MOFEM_LOG("EP", Sev::inform) << "spaceOrder: -space_order " << spaceOrder;
  MOFEM_LOG("EP", Sev::inform)
      << "spaceH1Order: -space_h1_order " << spaceH1Order;
  MOFEM_LOG("EP", Sev::inform)
      << "materialOrder: -material_order " << materialOrder;
  MOFEM_LOG("EP", Sev::inform) << "alphaU: -viscosity_alpha_u " << alphaU;
  MOFEM_LOG("EP", Sev::inform) << "alphaW: -viscosity_alpha_w " << alphaW;
  MOFEM_LOG("EP", Sev::inform)
      << "alphaOmega: -viscosity_alpha_omega " << alphaOmega;
  MOFEM_LOG("EP", Sev::inform) << "alphaRho: -density_alpha_rho " << alphaRho;
  MOFEM_LOG("EP", Sev::inform) << "alphaTau: -alpha_tau " << alphaTau;
  MOFEM_LOG("EP", Sev::inform)
      << "Rotations: -rotations " << list_rots[EshelbianCore::rotSelector];
  MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
                               << list_rots[EshelbianCore::gradApproximator];
  MOFEM_LOG("EP", Sev::inform)
      << "Symmetry: -symm " << list_symm[EshelbianCore::symmetrySelector];
  if (exponentBase != exp(1))
    MOFEM_LOG("EP", Sev::inform)
        << "Base exponent: -exponent_base " << EshelbianCore::exponentBase;
  else
    MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
  MOFEM_LOG("EP", Sev::inform)
      << "Stretch: -stretches " << list_stretches[choice_stretch];
  MOFEM_LOG("EP", Sev::inform) << "No stretch: -no_stretch " << (noStretch)
      ? "yes"
      : "no";
 
  MOFEM_LOG("EP", Sev::inform)
          << "Dynamic relaxation: -dynamic_relaxation " << (dynamicRelaxation)
      ? "yes"
      : "no";
  MOFEM_LOG("EP", Sev::inform)
      << "Solver type: -solver_type " << list_solvers[choice_solver];
  MOFEM_LOG("EP", Sev::inform)
          << "Singularity: -set_singularity " << (setSingularity)
      ? "yes"
      : "no";
  MOFEM_LOG("EP", Sev::inform)
          << "L2 user base scale: -l2_user_base_scale " << (l2UserBaseScale)
      ? "yes"
      : "no";
 
  MOFEM_LOG("EP", Sev::inform) << "Cracking on: -cracking_on " << (crackingOn)
      ? "yes"
      : "no";
  MOFEM_LOG("EP", Sev::inform)
      << "Cracking start time: -cracking_start_time " << crackingStartTime;
  MOFEM_LOG("EP", Sev::inform)
      << "Griffith energy: -griffith_energy " << griffithEnergy;
  MOFEM_LOG("EP", Sev::inform)
      << "Energy release variant: -energy_release_variant "
      << list_release[EshelbianCore::energyReleaseSelector];
  MOFEM_LOG("EP", Sev::inform)
      << "Number of J integral contours: -nb_J_integral_contours "
      << nbJIntegralContours;
  MOFEM_LOG("EP", Sev::inform)
          << "Cohesive interface on: -cohesive_interface_on "
          << (intefaceCrack == PETSC_TRUE)
      ? "yes"
      : "no";
 
#ifdef ENABLE_PYTHON_BINDING
  auto file_exists = [](std::string myfile) {
    std::ifstream file(myfile.c_str());
    if (file) {
      return true;
    }
    return false;
  };
 
  if (file_exists(analytical_expr_file_name)) {
    MOFEM_LOG("EP", Sev::inform) << analytical_expr_file_name << " file found";
 
    AnalyticalExprPythonPtr = boost::make_shared<AnalyticalExprPython>();
    CHKERR AnalyticalExprPythonPtr->analyticalExprInit(
        analytical_expr_file_name);
    AnalyticalExprPythonWeakPtr = AnalyticalExprPythonPtr;
  } else {
    MOFEM_LOG("EP", Sev::warning)
        << analytical_expr_file_name << " file NOT found";
  }
#endif
 
  if (spaceH1Order == -1)
    spaceH1Order = spaceOrder;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::addFields(const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  auto get_tets = [&]() {
    Range tets;
    CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
    return tets;
  };
 
  auto get_tets_skin = [&]() {
    Range tets_skin_part;
    Skinner skin(&mField.get_moab());
    CHKERR skin.find_skin(0, get_tets(), false, tets_skin_part);
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    Range tets_skin;
    CHKERR pcomm->filter_pstatus(tets_skin_part,
                                 PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &tets_skin);
    return tets_skin;
  };
 
  auto subtract_boundary_conditions = [&](auto &&tets_skin) {
    // That mean, that hybrid field on all faces on which traction is applied,
    // on other faces, or enforcing displacements as
    // natural boundary condition.
    if (bcSpatialTractionVecPtr)
      for (auto &v : *bcSpatialTractionVecPtr) {
        tets_skin = subtract(tets_skin, v.faces);
      }
    if (bcSpatialAnalyticalTractionVecPtr)
      for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
        tets_skin = subtract(tets_skin, v.faces);
      }
 
    if (bcSpatialPressureVecPtr)
      for (auto &v : *bcSpatialPressureVecPtr) {
        tets_skin = subtract(tets_skin, v.faces);
      }
 
    return tets_skin;
  };
 
  auto add_blockset = [&](auto block_name, auto &&tets_skin) {
    auto crack_faces =
        get_range_from_block(mField, "block_name", SPACE_DIM - 1);
    tets_skin.merge(crack_faces);
    return tets_skin;
  };
 
  auto subtract_blockset = [&](auto block_name, auto &&tets_skin) {
    auto contact_range =
        get_range_from_block(mField, block_name, SPACE_DIM - 1);
    tets_skin = subtract(tets_skin, contact_range);
    return tets_skin;
  };
 
  auto get_stress_trace_faces = [&](auto &&tets_skin) {
    Range faces;
    CHKERR mField.get_moab().get_adjacencies(get_tets(), SPACE_DIM - 1, true,
                                             faces, moab::Interface::UNION);
    Range trace_faces = subtract(faces, tets_skin);
    return trace_faces;
  };
 
  auto tets = get_tets();
 
  // remove also contact faces, i.e. that is also kind of hybrid field but
  // named but used to enforce contact conditions
  auto trace_faces = get_stress_trace_faces(
 
      subtract_blockset("CONTACT",
                        subtract_boundary_conditions(get_tets_skin()))
 
  );
 
  contactFaces = boost::make_shared<Range>(intersect(
      trace_faces, get_range_from_block(mField, "CONTACT", SPACE_DIM - 1)));
  skeletonFaces =
      boost::make_shared<Range>(subtract(trace_faces, *contactFaces));
 
#ifndef NDEBUG
  if (contactFaces->size())
    CHKERR save_range(mField.get_moab(),
                      "contact_faces_" +
                          std::to_string(mField.get_comm_rank()) + ".vtk",
                      *contactFaces);
  if (skeletonFaces->size())
    CHKERR save_range(mField.get_moab(),
                      "skeleton_faces_" +
                          std::to_string(mField.get_comm_rank()) + ".vtk",
                      *skeletonFaces);
#endif
 
  auto add_broken_hdiv_field = [this, meshset](const std::string field_name,
                                               const int order) {
    MoFEMFunctionBegin;
 
    constexpr FieldApproximationBase base = DEMKOWICZ_JACOBI_BASE;
 
    auto get_side_map_hdiv = [&]() {
      return std::vector<
 
          std::pair<EntityType,
                    std::function<MoFEMErrorCode(BaseFunction::DofsSideMap &)>
 
                    >>{
 
          {MBTET,
           [&](BaseFunction::DofsSideMap &dofs_side_map) -> MoFEMErrorCode {
             return TetPolynomialBase::setDofsSideMap(HDIV, DISCONTINUOUS, base,
                                                      dofs_side_map);
           }}
 
      };
    };
 
    CHKERR mField.add_broken_field(field_name, HDIV, base, SPACE_DIM,
                                   get_side_map_hdiv(), MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_l2_field = [this, meshset](const std::string field_name,
                                      const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_h1_field = [this, meshset](const std::string field_name,
                                      const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
    CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_l2_field_by_range = [this](const std::string field_name,
                                      const int order, const int dim,
                                      const int field_dim, Range &&r) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, field_dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(r);
    CHKERR mField.add_ents_to_field_by_dim(r, dim, field_name);
    CHKERR mField.set_field_order(r, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_bubble_field = [this, meshset](const std::string field_name,
                                          const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
                            MF_ZERO);
    // Modify field
    auto field_ptr = mField.get_field_structure(field_name);
    auto field_order_table =
        const_cast<Field *>(field_ptr)->getFieldOrderTable();
    auto get_cgg_bubble_order_zero = [](int p) { return 0; };
    auto get_cgg_bubble_order_tet = [](int p) {
      return NBVOLUMETET_CCG_BUBBLE(p);
    };
    field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
    field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
    field_order_table[MBTRI] = get_cgg_bubble_order_zero;
    field_order_table[MBTET] = get_cgg_bubble_order_tet;
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_user_l2_field = [this, meshset](const std::string field_name,
                                           const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, L2, USER_BASE, dim, MB_TAG_DENSE,
                            MF_ZERO);
    // Modify field
    auto field_ptr = mField.get_field_structure(field_name);
    auto field_order_table =
        const_cast<Field *>(field_ptr)->getFieldOrderTable();
    auto zero_dofs = [](int p) { return 0; };
    auto dof_l2_tet = [](int p) { return NBVOLUMETET_L2(p); };
    field_order_table[MBVERTEX] = zero_dofs;
    field_order_table[MBEDGE] = zero_dofs;
    field_order_table[MBTRI] = zero_dofs;
    field_order_table[MBTET] = dof_l2_tet;
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  // spatial fields
  CHKERR add_broken_hdiv_field(piolaStress, spaceOrder);
  CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
  CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
  CHKERR add_user_l2_field(rotAxis, spaceOrder - 1, 3);
  CHKERR add_user_l2_field(stretchTensor, noStretch ? -1 : spaceOrder, 6);
 
  if (!skeletonFaces)
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
  if (!contactFaces)
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No contact faces");
 
  auto get_hybridised_disp = [&]() {
    auto faces = *skeletonFaces;
    auto skin = subtract_boundary_conditions(get_tets_skin());
    for (auto &bc : *bcSpatialNormalDisplacementVecPtr) {
      faces.merge(intersect(bc.faces, skin));
    }
    return faces;
  };
 
  CHKERR add_l2_field_by_range(hybridSpatialDisp, spaceOrder - 1, 2, 3,
                               get_hybridised_disp());
  CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
                               Range(*contactFaces));
 
  // spatial displacement
  CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
  // material positions
  CHKERR add_h1_field(materialH1Positions, 2, 3);
 
  // Eshelby stress
  // CHKERR add_broken_hdiv_field(eshelbyStress, spaceOrder);
  // CHKERR add_l2_field(materialL2Disp, spaceOrder - 1, 3);
  // CHKERR add_l2_field_by_range(hybridMaterialDisp, spaceOrder - 1, 2, 3,
  //                              Range(*materialSkeletonFaces));
 
  CHKERR mField.build_fields();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::projectGeometry(const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  Range meshset_ents;
  CHKERR mField.get_moab().get_entities_by_handle(meshset, meshset_ents);
 
  auto project_ho_geometry = [&](auto field) {
    Projection10NodeCoordsOnField ent_method(mField, materialH1Positions);
    return mField.loop_dofs(field, ent_method);
  };
  CHKERR project_ho_geometry(materialH1Positions);
 
  auto get_adj_front_edges = [&](auto &front_edges) {
    Range front_crack_nodes;
    Range crack_front_edges_with_both_nodes_not_at_front;
 
    if (mField.get_comm_rank() == 0) {
      auto &moab = mField.get_moab();
      CHK_MOAB_THROW(
          moab.get_connectivity(front_edges, front_crack_nodes, true),
          "get_connectivity failed");
      Range crack_front_edges;
      CHK_MOAB_THROW(moab.get_adjacencies(front_crack_nodes, SPACE_DIM - 2,
                                          false, crack_front_edges,
                                          moab::Interface::UNION),
                     "get_adjacencies failed");
      Range crack_front_edges_nodes;
      CHK_MOAB_THROW(moab.get_connectivity(crack_front_edges,
                                           crack_front_edges_nodes, true),
                     "get_connectivity failed");
      // those nodes are hannging nodes
      crack_front_edges_nodes =
          subtract(crack_front_edges_nodes, front_crack_nodes);
      Range crack_front_edges_with_both_nodes_not_at_front;
      CHK_MOAB_THROW(
          moab.get_adjacencies(crack_front_edges_nodes, 1, false,
                               crack_front_edges_with_both_nodes_not_at_front,
                               moab::Interface::UNION),
          "get_adjacencies failed");
      // those edges are have one node not at the crack front
      crack_front_edges_with_both_nodes_not_at_front = intersect(
          crack_front_edges, crack_front_edges_with_both_nodes_not_at_front);
    }
 
    front_crack_nodes = send_type(mField, front_crack_nodes, MBVERTEX);
    crack_front_edges_with_both_nodes_not_at_front = send_type(
        mField, crack_front_edges_with_both_nodes_not_at_front, MBEDGE);
 
    return std::make_pair(boost::make_shared<Range>(front_crack_nodes),
                          boost::make_shared<Range>(
                              crack_front_edges_with_both_nodes_not_at_front));
  };
 
  crackFaces = boost::make_shared<Range>(
      get_range_from_block(mField, "CRACK", SPACE_DIM - 1));
  frontEdges =
      boost::make_shared<Range>(get_crack_front_edges(mField, *crackFaces));
  auto [front_vertices, front_adj_edges] = get_adj_front_edges(*frontEdges);
  frontVertices = front_vertices;
  frontAdjEdges = front_adj_edges;
 
  MOFEM_LOG("EP", Sev::inform)
      << "Number of crack faces: " << crackFaces->size();
  MOFEM_LOG("EP", Sev::inform)
      << "Number of front edges: " << frontEdges->size();
  MOFEM_LOG("EP", Sev::inform)
      << "Number of front vertices: " << frontVertices->size();
  MOFEM_LOG("EP", Sev::inform)
      << "Number of front adjacent edges: " << frontAdjEdges->size();
 
#ifndef NDEBUG
  if (crackingOn) {
    auto rank = mField.get_comm_rank();
    // CHKERR save_range(mField.get_moab(),
    //                   (boost::format("meshset_ents_%d.vtk") % rank).str(),
    //                   meshset_ents);
    CHKERR save_range(mField.get_moab(),
                      (boost::format("crack_faces_%d.vtk") % rank).str(),
                      *crackFaces);
    CHKERR save_range(mField.get_moab(),
                      (boost::format("front_edges_%d.vtk") % rank).str(),
                      *frontEdges);
    // CHKERR save_range(mField.get_moab(),
    //                   (boost::format("front_vertices_%d.vtk") % rank).str(),
    //                   *frontVertices);
    // CHKERR save_range(mField.get_moab(),
    //                   (boost::format("front_adj_edges_%d.vtk") % rank).str(),
    //                   *frontAdjEdges);
  }
#endif // NDEBUG
 
  auto set_singular_dofs = [&](auto &front_adj_edges, auto &front_vertices) {
    MoFEMFunctionBegin;
    auto &moab = mField.get_moab();
 
    double eps = 1;
    double beta = 0;
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-singularity_eps", &beta,
                                 PETSC_NULLPTR);
    MOFEM_LOG("EP", Sev::inform) << "Singularity eps " << beta;
    eps -= beta;
 
    auto field_blas = mField.getInterface<FieldBlas>();
    auto lambda =
        [&](boost::shared_ptr<FieldEntity> field_entity_ptr) -> MoFEMErrorCode {
      MoFEMFunctionBegin;
      FTENSOR_INDEX(3, i);
      FTENSOR_INDEX(3, j);
 
      auto nb_dofs = field_entity_ptr->getEntFieldData().size();
      if (nb_dofs == 0) {
        MoFEMFunctionReturnHot(0);
      }
 
#ifndef NDEBUG
      if (field_entity_ptr->getNbOfCoeffs() != 3)
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "Expected 3 coefficients per edge");
      if (nb_dofs % 3 != 0)
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "Expected multiple of 3 coefficients per edge");
#endif // NDEBUG
 
      auto get_conn = [&]() {
        int num_nodes;
        const EntityHandle *conn;
        CHKERR moab.get_connectivity(field_entity_ptr->getEnt(), conn,
                                     num_nodes, false);
        return std::make_pair(conn, num_nodes);
      };
 
      auto get_dir = [&](auto &&conn_p) {
        auto [conn, num_nodes] = conn_p;
        double coords[6];
        CHKERR moab.get_coords(conn, num_nodes, coords);
        FTensor::Tensor1<double, 3> t_edge_dir{coords[3] - coords[0],
                                               coords[4] - coords[1],
                                               coords[5] - coords[2]};
        return t_edge_dir;
      };
 
      auto get_singularity_dof = [&](auto &&conn_p, auto &&t_edge_dir) {
        auto [conn, num_nodes] = conn_p;
        FTensor::Tensor1<double, 3> t_singularity_dof{0., 0., 0.};
        if (front_vertices.find(conn[0]) != front_vertices.end()) {
          t_singularity_dof(i) = t_edge_dir(i) * (-eps);
        } else if (front_vertices.find(conn[1]) != front_vertices.end()) {
          t_singularity_dof(i) = t_edge_dir(i) * eps;
        }
        return t_singularity_dof;
      };
 
      auto t_singularity_dof =
          get_singularity_dof(get_conn(), get_dir(get_conn()));
 
      auto field_data = field_entity_ptr->getEntFieldData();
      FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_dof{
          &field_data[0], &field_data[1], &field_data[2]};
 
      t_dof(i) = t_singularity_dof(i);
      ++t_dof;
      for (auto n = 1; n < field_data.size() / 3; ++n) {
        t_dof(i) = 0;
        ++t_dof;
      }
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR field_blas->fieldLambdaOnEntities(lambda, materialH1Positions,
                                             &front_adj_edges);
 
    MoFEMFunctionReturn(0);
  };
 
  if (setSingularity)
    CHKERR set_singular_dofs(*frontAdjEdges, *frontVertices);
 
  interfaceFaces = boost::make_shared<Range>(
      get_range_from_block(mField, "INTERFACE", SPACE_DIM - 1));
  MOFEM_LOG("EP", Sev::inform)
      << "Number of interface elements: " << interfaceFaces->size();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EshelbianCore::projectInternalStress(const EntityHandle meshset) {
  MoFEMFunctionBegin;
#ifdef INCLUDE_MBCOUPLER
 
  char mesh_source_file[255] = "source.h5m";
  char iterp_tag_name[255] = "INTERNAL_STRESS";
 
  int interp_order = 1;
  PetscBool hybrid_interp = PETSC_TRUE;
  PetscBool project_internal_stress = PETSC_FALSE;
 
  double toler = 5.e-10;
  PetscOptionsBegin(PETSC_COMM_WORLD, "mesh_transfer_", "mesh data transfer",
                    "none");
  CHKERR PetscOptionsString("-source_file", "source mesh file name", "",
                            "source.h5m", mesh_source_file, 255,
                            &project_internal_stress);
  CHKERR PetscOptionsString("-interp_tag", "interpolation tag name", "",
                            "INTERNAL_STRESS", iterp_tag_name, 255,
                            PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-interp_order", "interpolation order", "", 0,
                         &interp_order, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-hybrid_interp", "use hybrid interpolation", "",
                          hybrid_interp, &hybrid_interp, PETSC_NULLPTR);
  PetscOptionsEnd();
 
  MOFEM_LOG_CHANNEL("WORLD");
  MOFEM_LOG_TAG("WORLD", "mesh_data_transfer");
  if (!project_internal_stress) {
    MOFEM_LOG("WORLD", Sev::inform)
        << "No internal stress source mesh specified. Skipping projection of "
           "internal stress";
    MoFEMFunctionReturnHot(0);
  }
  MOFEM_LOG("WORLD", Sev::inform)
      << "Projecting internal stress from source mesh: " << mesh_source_file;
 
  auto &moab = mField.get_moab();
 
  // check if tag exists
  Tag old_interp_tag;
  auto rval_check_tag = moab.tag_get_handle(iterp_tag_name, old_interp_tag);
  if (rval_check_tag == MB_SUCCESS) {
    MOFEM_LOG("WORLD", Sev::inform)
        << "Deleting existing tag on target mesh: " << iterp_tag_name;
    CHKERR moab.tag_delete(old_interp_tag);
  }
 
  // make a size-1 communicator for the coupler (rank 0 only)
  int world_rank = -1, world_size = -1;
  MPI_Comm_rank(PETSC_COMM_WORLD, &world_rank);
  MPI_Comm_size(PETSC_COMM_WORLD, &world_size);
 
  Range original_meshset_ents;
  CHKERR moab.get_entities_by_handle(0, original_meshset_ents);
 
  MPI_Comm comm_coupler;
  if (world_rank == 0) {
    MPI_Comm_split(PETSC_COMM_WORLD, 0, 0, &comm_coupler);
  } else {
    MPI_Comm_split(PETSC_COMM_WORLD, MPI_UNDEFINED, world_rank, &comm_coupler);
  }
 
  // build a separate ParallelComm for the coupler (rank 0 only)
  ParallelComm *pcomm0 = nullptr;
  int pcomm0_id = -1;
  if (world_rank == 0) {
    pcomm0 = new ParallelComm(&moab, comm_coupler, &pcomm0_id);
  }
 
  Coupler::Method method;
  switch (interp_order) {
  case 0:
    method = Coupler::CONSTANT;
    break;
  case 1:
    method = Coupler::LINEAR_FE;
    break;
  default:
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
            "Unsupported interpolation order");
  }
 
  int nprocs, rank;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD, &nprocs);
  CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  CHKERRQ(ierr);
 
  // std::string read_opts, write_opts;
  // read_opts = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_"
  //             "DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS";
  // if (world_size > 1)
  //   read_opts += ";PARALLEL_GHOSTS=3.0.1";
  // write_opts = (world_size > 1) ? "PARALLEL=WRITE_PART" : "";
 
  // create target mesh from existing meshset
  EntityHandle target_root;
  CHKERR moab.create_meshset(MESHSET_SET, target_root);
  MOFEM_LOG("WORLD", Sev::inform)
      << "Creating target mesh from existing meshset";
  Range target_meshset_ents;
  CHKERR moab.get_entities_by_handle(0, target_meshset_ents);
  CHKERR moab.add_entities(target_root, target_meshset_ents);
 
  // variables for tags to be broadcast later
  int tag_length;
  DataType dtype;
  TagType storage;
 
  // load source mesh
  Range targ_verts, targ_elems;
  if (world_rank == 0) {
 
    EntityHandle source_root;
    CHKERR moab.create_meshset(MESHSET_SET, source_root);
 
    MOFEM_LOG("WORLD", Sev::inform) << "Loading source mesh on rank 0";
    auto rval_source_mesh = moab.load_file(mesh_source_file, &source_root, "");
    if (rval_source_mesh != MB_SUCCESS) {
      MOFEM_LOG("WORLD", Sev::warning)
          << "Error loading source mesh file: " << mesh_source_file;
    }
    MOFEM_LOG("WORLD", Sev::inform) << "Source mesh loaded.";
 
    Range src_elems;
    CHKERR moab.get_entities_by_dimension(source_root, 3, src_elems);
 
    EntityHandle part_set;
    CHKERR pcomm0->create_part(part_set);
    CHKERR moab.add_entities(part_set, src_elems);
 
    Range src_elems_part;
    CHKERR pcomm0->get_part_entities(src_elems_part, 3);
 
    Tag interp_tag;
    CHKERR moab.tag_get_handle(iterp_tag_name, interp_tag);
 
    int interp_tag_len;
    CHKERR moab.tag_get_length(interp_tag, interp_tag_len);
 
    if (interp_tag_len != 1 && interp_tag_len != 3 && interp_tag_len != 9) {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
              "Unsupported interpolation tag length: %d", interp_tag_len);
    }
 
    // store tag info for later broadcast
    tag_length = interp_tag_len;
    CHKERR moab.tag_get_data_type(interp_tag, dtype);
    CHKERR moab.tag_get_type(interp_tag, storage);
 
    // coupler is collective
    Coupler mbc(&moab, pcomm0, src_elems_part, 0, true);
 
    std::vector<double> vpos; // the positions we are interested in
    int num_pts = 0;
 
    Range tmp_verts;
 
    // First get all vertices adj to partition entities in target mesh
    CHKERR moab.get_entities_by_dimension(target_root, 3, targ_elems);
 
    if (interp_order == 0) {
      targ_verts = targ_elems;
    } else {
      CHKERR moab.get_adjacencies(targ_elems, 0, false, targ_verts,
                                  moab::Interface::UNION);
    }
 
    // Then get non-owned verts and subtract
    CHKERR pcomm0->get_pstatus_entities(0, PSTATUS_NOT_OWNED, tmp_verts);
    targ_verts = subtract(targ_verts, tmp_verts);
 
    // get position of these entities; these are the target points
    num_pts = (int)targ_verts.size();
    vpos.resize(3 * targ_verts.size());
    CHKERR moab.get_coords(targ_verts, &vpos[0]);
 
    // Locate those points in the source mesh
    boost::shared_ptr<TupleList> tl_ptr;
    tl_ptr = boost::make_shared<TupleList>();
    CHKERR mbc.locate_points(&vpos[0], num_pts, 0, toler, tl_ptr.get(), false);
 
    // If some points were not located, we need to process them
    auto find_missing_points = [&](Range &targ_verts, int &num_pts,
                                   std::vector<double> &vpos,
                                   Range &missing_verts) {
      MoFEMFunctionBegin;
      int missing_pts_num = 0;
      int i = 0;
      auto vit = targ_verts.begin();
      for (; vit != targ_verts.end(); i++) {
        if (tl_ptr->vi_rd[3 * i + 1] == -1) {
          missing_verts.insert(*vit);
          vit = targ_verts.erase(vit);
          missing_pts_num++;
        } else {
          vit++;
        }
      }
 
      int missing_pts_num_global = 0;
      // MPI_Allreduce(&missing_pts_num, &missing_pts_num_global, 1, MPI_INT,
      //               MPI_SUM, pcomm0);
      if (missing_pts_num_global) {
        MOFEM_LOG("WORLD", Sev::warning)
            << missing_pts_num_global
            << " points in target mesh were not located in source mesh. ";
      }
 
      if (missing_pts_num) {
        num_pts = (int)targ_verts.size();
        vpos.resize(3 * targ_verts.size());
        CHKERR moab.get_coords(targ_verts, &vpos[0]);
        tl_ptr->reset();
        CHKERR mbc.locate_points(&vpos[0], num_pts, 0, toler, tl_ptr.get(),
                                 false);
      }
      MoFEMFunctionReturn(0);
    };
 
    Range missing_verts;
    CHKERR find_missing_points(targ_verts, num_pts, vpos, missing_verts);
 
    std::vector<double> source_data(interp_tag_len * src_elems.size(), 0.0);
    std::vector<double> target_data(interp_tag_len * num_pts, 0.0);
 
    CHKERR moab.tag_get_data(interp_tag, src_elems, &source_data[0]);
 
    Tag scalar_tag, adj_count_tag;
    double def_scl = 0;
    string scalar_tag_name = string(iterp_tag_name) + "_COMP";
    CHKERR moab.tag_get_handle(scalar_tag_name.c_str(), 1, MB_TYPE_DOUBLE,
                               scalar_tag, MB_TAG_CREAT | MB_TAG_DENSE,
                               &def_scl);
 
    string adj_count_tag_name = "ADJ_COUNT";
    double def_adj = 0;
    CHKERR moab.tag_get_handle(adj_count_tag_name.c_str(), 1, MB_TYPE_DOUBLE,
                               adj_count_tag, MB_TAG_CREAT | MB_TAG_DENSE,
                               &def_adj);
 
    // MBCoupler functionality supports only scalar tags. For the case of
    // vector or tensor tags we need to save each component as a scalar tag
    auto create_scalar_tags = [&](const Range &src_elems,
                                  const std::vector<double> &source_data,
                                  int itag) {
      MoFEMFunctionBegin;
 
      std::vector<double> source_data_scalar(src_elems.size());
      // Populate source_data_scalar
      for (int ielem = 0; ielem < src_elems.size(); ielem++) {
        source_data_scalar[ielem] = source_data[itag + ielem * interp_tag_len];
      }
 
      // Set data on the scalar tag
      CHKERR moab.tag_set_data(scalar_tag, src_elems, &source_data_scalar[0]);
 
      if (interp_order == 1) {
        // Linear interpolation: compute average value of data on vertices
        Range src_verts;
        CHKERR moab.get_connectivity(src_elems, src_verts, true);
 
        CHKERR moab.tag_clear_data(scalar_tag, src_verts, &def_scl);
        CHKERR moab.tag_clear_data(adj_count_tag, src_verts, &def_adj);
 
        for (auto &tet : src_elems) {
          double tet_data = 0;
          CHKERR moab.tag_get_data(scalar_tag, &tet, 1, &tet_data);
 
          Range adj_verts;
          CHKERR moab.get_connectivity(&tet, 1, adj_verts, true);
 
          std::vector<double> adj_vert_data(adj_verts.size(), 0.0);
          std::vector<double> adj_vert_count(adj_verts.size(), 0.0);
 
          CHKERR moab.tag_get_data(scalar_tag, adj_verts, &adj_vert_data[0]);
          CHKERR moab.tag_get_data(adj_count_tag, adj_verts,
                                   &adj_vert_count[0]);
 
          for (int ivert = 0; ivert < adj_verts.size(); ivert++) {
            adj_vert_data[ivert] += tet_data;
            adj_vert_count[ivert] += 1;
          }
 
          CHKERR moab.tag_set_data(scalar_tag, adj_verts, &adj_vert_data[0]);
          CHKERR moab.tag_set_data(adj_count_tag, adj_verts,
                                   &adj_vert_count[0]);
        }
 
        // Reduce tags for the parallel case
        std::vector<Tag> tags = {scalar_tag, adj_count_tag};
        pcomm0->reduce_tags(tags, tags, MPI_SUM, src_verts);
 
        std::vector<double> src_vert_data(src_verts.size(), 0.0);
        std::vector<double> src_vert_adj_count(src_verts.size(), 0.0);
 
        CHKERR moab.tag_get_data(scalar_tag, src_verts, &src_vert_data[0]);
        CHKERR moab.tag_get_data(adj_count_tag, src_verts,
                                 &src_vert_adj_count[0]);
 
        for (int ivert = 0; ivert < src_verts.size(); ivert++) {
          src_vert_data[ivert] /= src_vert_adj_count[ivert];
        }
        CHKERR moab.tag_set_data(scalar_tag, src_verts, &src_vert_data[0]);
      }
      MoFEMFunctionReturn(0);
    };
 
    for (int itag = 0; itag < interp_tag_len; itag++) {
 
      CHKERR create_scalar_tags(src_elems, source_data, itag);
 
      std::vector<double> target_data_scalar(num_pts, 0.0);
      CHKERR mbc.interpolate(method, scalar_tag_name, &target_data_scalar[0],
                             tl_ptr.get());
 
      for (int ielem = 0; ielem < num_pts; ielem++) {
        target_data[itag + ielem * interp_tag_len] = target_data_scalar[ielem];
      }
    }
 
    // Use original tag
    CHKERR moab.tag_set_data(interp_tag, targ_verts, &target_data[0]);
 
    if (missing_verts.size() && (interp_order == 1) && hybrid_interp) {
      MOFEM_LOG("WORLD", Sev::warning) << "Using hybrid interpolation for "
                                          "missing points in the target mesh.";
      Range missing_adj_elems;
      CHKERR moab.get_adjacencies(missing_verts, 3, false, missing_adj_elems,
                                  moab::Interface::UNION);
 
      int num_adj_elems = (int)missing_adj_elems.size();
      std::vector<double> vpos_adj_elems;
 
      vpos_adj_elems.resize(3 * missing_adj_elems.size());
      CHKERR moab.get_coords(missing_adj_elems, &vpos_adj_elems[0]);
 
      // Locate those points in the source mesh
      tl_ptr->reset();
      CHKERR mbc.locate_points(&vpos_adj_elems[0], num_adj_elems, 0, toler,
                               tl_ptr.get(), false);
 
      Range missing_tets;
      CHKERR find_missing_points(missing_adj_elems, num_adj_elems,
                                 vpos_adj_elems, missing_tets);
      if (missing_tets.size()) {
        MOFEM_LOG("WORLD", Sev::warning)
            << missing_tets.size()
            << "points in target mesh were not located in source mesh. ";
      }
 
      std::vector<double> target_data_adj_elems(interp_tag_len * num_adj_elems,
                                                0.0);
 
      for (int itag = 0; itag < interp_tag_len; itag++) {
        CHKERR create_scalar_tags(src_elems, source_data, itag);
 
        std::vector<double> target_data_adj_elems_scalar(num_adj_elems, 0.0);
        CHKERR mbc.interpolate(method, scalar_tag_name,
                               &target_data_adj_elems_scalar[0], tl_ptr.get());
 
        for (int ielem = 0; ielem < num_adj_elems; ielem++) {
          target_data_adj_elems[itag + ielem * interp_tag_len] =
              target_data_adj_elems_scalar[ielem];
        }
      }
 
      CHKERR moab.tag_set_data(interp_tag, missing_adj_elems,
                               &target_data_adj_elems[0]);
 
      // FIXME: add implementation for parallel case
      for (auto &vert : missing_verts) {
        Range adj_elems;
        CHKERR moab.get_adjacencies(&vert, 1, 3, false, adj_elems,
                                    moab::Interface::UNION);
 
        std::vector<double> adj_elems_data(adj_elems.size() * interp_tag_len,
                                           0.0);
        CHKERR moab.tag_get_data(interp_tag, adj_elems, &adj_elems_data[0]);
 
        std::vector<double> vert_data(interp_tag_len, 0.0);
        for (int itag = 0; itag < interp_tag_len; itag++) {
          for (int i = 0; i < adj_elems.size(); i++) {
            vert_data[itag] += adj_elems_data[i * interp_tag_len + itag];
          }
          vert_data[itag] /= adj_elems.size();
        }
        CHKERR moab.tag_set_data(interp_tag, &vert, 1, &vert_data[0]);
      }
    }
 
    CHKERR moab.tag_delete(scalar_tag);
    CHKERR moab.tag_delete(adj_count_tag);
    // delete source mesh
    Range src_mesh_ents;
    CHKERR moab.get_entities_by_handle(source_root, src_mesh_ents);
    CHKERR moab.delete_entities(&source_root, 1);
    CHKERR moab.delete_entities(src_mesh_ents);
    CHKERR moab.delete_entities(&part_set, 1);
  }
 
  // broadcast tag info to other processors
  MPI_Bcast(&tag_length, 1, MPI_INT, 0, PETSC_COMM_WORLD);
  MPI_Bcast(&dtype, 1, MPI_INT, 0, PETSC_COMM_WORLD);
  MPI_Bcast(&storage, 1, MPI_INT, 0, PETSC_COMM_WORLD);
 
  // create tag on other processors
  Tag interp_tag_all;
  unsigned flags =
      MB_TAG_CREAT | storage; // e.g., MB_TAG_DENSE or MB_TAG_SPARSE
  std::vector<double> def_val(tag_length, 0.);
  CHKERR moab.tag_get_handle(iterp_tag_name, tag_length, dtype, interp_tag_all,
                             flags, def_val.data());
  MPI_Barrier(PETSC_COMM_WORLD);
 
  // exchange data for all entity types across all processors
  auto vertex_exchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 0, tag_length, Sev::inform);
  auto edge_exchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 1, tag_length, Sev::inform);
  auto face_exchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 2, tag_length, Sev::inform);
  auto volume_exchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 3, tag_length, Sev::inform);
 
  CHKERR CommInterface::updateEntitiesPetscVector(
      mField.get_moab(), vertex_exchange, interp_tag_all);
  CHKERR CommInterface::updateEntitiesPetscVector(
      mField.get_moab(), edge_exchange, interp_tag_all);
  CHKERR CommInterface::updateEntitiesPetscVector(
      mField.get_moab(), face_exchange, interp_tag_all);
  CHKERR CommInterface::updateEntitiesPetscVector(
      mField.get_moab(), volume_exchange, interp_tag_all);
 
  // delete target meshset but not the entities
  CHKERR moab.delete_entities(&target_root, 1);
 
#endif // INCLUDE_MBCOUPLER
  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);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EshelbianCore::addBoundaryFiniteElement(const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  Range meshset_ents;
  CHKERR mField.get_moab().get_entities_by_handle(meshset, meshset_ents);
 
  auto set_fe_adjacency = [&](auto fe_name) {
    MoFEMFunctionBegin;
    parentAdjSkeletonFunctionDim2 =
        boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
            bitAdjParent, bitAdjParentMask, bitAdjEnt, bitAdjEntMask);
    CHKERR mField.modify_finite_element_adjacency_table(
        fe_name, MBTRI, *parentAdjSkeletonFunctionDim2);
    MoFEMFunctionReturn(0);
  };
 
  // set finite element fields
  auto add_field_to_fe = [this](const std::string fe,
                                const std::string field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, materialH1Positions);
    MoFEMFunctionReturn(0);
  };
 
  if (!mField.check_finite_element(naturalBcElement)) {
 
    Range natural_bc_elements;
    if (bcSpatialDispVecPtr) {
      for (auto &v : *bcSpatialDispVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialRotationVecPtr) {
      for (auto &v : *bcSpatialRotationVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialNormalDisplacementVecPtr) {
      for (auto &v : *bcSpatialNormalDisplacementVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialAnalyticalDisplacementVecPtr) {
      for (auto &v : *bcSpatialAnalyticalDisplacementVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialTractionVecPtr) {
      for (auto &v : *bcSpatialTractionVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialAnalyticalTractionVecPtr) {
      for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialPressureVecPtr) {
      for (auto &v : *bcSpatialPressureVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    natural_bc_elements = intersect(natural_bc_elements, meshset_ents);
 
    CHKERR mField.add_finite_element(naturalBcElement, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
                                                     naturalBcElement);
    CHKERR add_field_to_fe(naturalBcElement, piolaStress);
    CHKERR add_field_to_fe(naturalBcElement, hybridSpatialDisp);
    CHKERR set_fe_adjacency(naturalBcElement);
    CHKERR mField.build_finite_elements(naturalBcElement);
  }
 
  auto get_skin = [&](auto &body_ents) {
    Skinner skin(&mField.get_moab());
    Range skin_ents;
    CHKERR skin.find_skin(0, body_ents, false, skin_ents);
    return skin_ents;
  };
 
  auto filter_true_skin = [&](auto &&skin) {
    Range boundary_ents;
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &boundary_ents);
    return boundary_ents;
  };
 
  if (!mField.check_finite_element(skinElement)) {
 
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM,
                                                       body_ents);
    auto skin = filter_true_skin(get_skin(body_ents));
 
    CHKERR mField.add_finite_element(skinElement, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(skin, MBTRI, skinElement);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       spatialH1Disp);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       materialH1Positions);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       contactDisp);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       hybridSpatialDisp);
 
    CHKERR 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, contactDisp);
      CHKERR add_field_to_fe(contactElement, spatialL2Disp);
      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, spatialL2Disp);
    CHKERR add_field_to_fe(skeletonElement, spatialH1Disp);
    CHKERR set_fe_adjacency(skeletonElement);
    CHKERR mField.build_finite_elements(skeletonElement);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::addDMs(const BitRefLevel bit,
                                     const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  // find adjacencies between finite elements and dofs
  CHKERR mField.build_adjacencies(bit, QUIET);
 
  // Create coupled problem
  dM = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
                            BitRefLevel().set());
  CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
  CHKERR DMMoFEMAddElement(dM, elementVolumeName);
  CHKERR DMMoFEMAddElement(dM, naturalBcElement);
  CHKERR DMMoFEMAddElement(dM, skeletonElement);
  CHKERR DMMoFEMAddElement(dM, contactElement);
  CHKERR DMMoFEMAddElement(dM, skinElement);
 
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
  CHKERR DMSetUp(dM);
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
 
  auto remove_dofs_on_broken_skin = [&](const std::string prb_name) {
    MoFEMFunctionBegin;
    for (int d : {0, 1, 2}) {
      std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
      CHKERR mField.getInterface<ProblemsManager>()
          ->getSideDofsOnBrokenSpaceEntities(
              dofs_to_remove, prb_name, ROW, piolaStress,
              (*bcSpatialFreeTractionVecPtr)[d], SPACE_DIM, d, d);
      // remove piola dofs, i.e. traction free boundary
      CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, ROW,
                                                                dofs_to_remove);
      CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, COL,
                                                                dofs_to_remove);
    }
    MoFEMFunctionReturn(0);
  };
  CHKERR remove_dofs_on_broken_skin("ESHELBY_PLASTICITY");
 
  // Create elastic sub-problem
  dmElastic = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
  CHKERR DMMoFEMSetDestroyProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, piolaStress);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis);
  if (!noStretch) {
    CHKERR DMMoFEMAddSubFieldRow(dmElastic, stretchTensor);
  }
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, hybridSpatialDisp);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, contactDisp);
  CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName);
  CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement);
  CHKERR DMMoFEMAddElement(dmElastic, skeletonElement);
  CHKERR DMMoFEMAddElement(dmElastic, contactElement);
  CHKERR DMMoFEMAddElement(dmElastic, skinElement);
  CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
  CHKERR DMSetUp(dmElastic);
 
  // dmMaterial = createDM(mField.get_comm(), "DMMOFEM");
  // CHKERR DMMoFEMCreateSubDM(dmMaterial, dM, "MATERIAL_PROBLEM");
  // CHKERR DMMoFEMSetDestroyProblem(dmMaterial, PETSC_TRUE);
  // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, eshelbyStress);
  // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, materialL2Disp);
  // CHKERR DMMoFEMAddElement(dmMaterh elementVolumeName);
  // CHKERR DMMoFEMAddElement(dmMaterial, naturalBcElement);
  // CHKERR DMMoFEMAddElement(dmMaterial, skinElement);
  // CHKERR DMMoFEMSetSquareProblem(dmMaterial, PETSC_TRUE);
  // CHKERR DMMoFEMSetIsPartitioned(dmMaterial, PETSC_TRUE);
  // CHKERR DMSetUp(dmMaterial);
 
  auto set_zero_block = [&]() {
    MoFEMFunctionBegin;
    if (!noStretch) {
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", spatialL2Disp, stretchTensor);
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", stretchTensor, spatialL2Disp);
    }
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", spatialL2Disp, rotAxis);
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", rotAxis, spatialL2Disp);
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", spatialL2Disp, bubbleField);
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", bubbleField, spatialL2Disp);
    if (!noStretch) {
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", bubbleField, bubbleField);
      CHKERR
      mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", piolaStress, piolaStress);
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", bubbleField, piolaStress);
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", piolaStress, bubbleField);
    }
 
    solTSStep = createDMVector(dmElastic);
    MoFEMFunctionReturn(0);
  };
 
  auto set_section = [&]() {
    MoFEMFunctionBegin;
    PetscSection section;
    CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
                                                           &section);
    CHKERR DMSetSection(dmElastic, section);
    CHKERR DMSetGlobalSection(dmElastic, section);
    CHKERR PetscSectionDestroy(&section);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_zero_block();
  CHKERR set_section();
 
  dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
  CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, spatialH1Disp);
  CHKERR DMMoFEMAddElement(dmPrjSpatial, elementVolumeName);
  CHKERR DMMoFEMSetSquareProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMSetUp(dmPrjSpatial);
 
  // CHKERR mField.getInterface<BcManager>()
  //     ->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
  //         "PROJECT_SPATIAL", spatialH1Disp, true, false);
 
  MoFEMFunctionReturn(0);
}
 
BcDisp::BcDisp(std::string name, std::vector<double> attr, Range faces)
    : blockName(name), faces(faces) {
  vals.resize(3, false);
  flags.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    vals[ii] = attr[ii];
    flags[ii] = static_cast<int>(attr[ii + 3]);
  }
 
  MOFEM_LOG("EP", Sev::inform) << "Add BCDisp " << name;
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCDisp vals " << vals[0] << " " << vals[1] << " " << vals[2];
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCDisp flags " << flags[0] << " " << flags[1] << " " << flags[2];
  MOFEM_LOG("EP", Sev::inform) << "Add BCDisp nb. of faces " << faces.size();
}
 
BcRot::BcRot(std::string name, std::vector<double> attr, Range faces)
    : blockName(name), faces(faces) {
  vals.resize(attr.size(), false);
  for (int ii = 0; ii != attr.size(); ++ii) {
    vals[ii] = attr[ii];
  }
  theta = attr[3];
}
 
TractionBc::TractionBc(std::string name, std::vector<double> attr, Range faces)
    : blockName(name), faces(faces) {
  vals.resize(3, false);
  flags.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    vals[ii] = attr[ii];
    flags[ii] = static_cast<int>(attr[ii + 3]);
  }
 
  MOFEM_LOG("EP", Sev::inform) << "Add BCForce " << name;
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCForce vals " << vals[0] << " " << vals[1] << " " << vals[2];
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCForce flags " << flags[0] << " " << flags[1] << " " << flags[2];
  MOFEM_LOG("EP", Sev::inform) << "Add BCForce nb. of faces " << faces.size();
}
 
NormalDisplacementBc::NormalDisplacementBc(std::string name,
                                           std::vector<double> attr,
                                           Range faces)
    : blockName(name), faces(faces) {
 
  blockName = name;
  if (attr.size() < 1) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong size of normal displacement BC");
  }
 
  val = attr[0];
 
  MOFEM_LOG("EP", Sev::inform) << "Add NormalDisplacementBc " << name;
  MOFEM_LOG("EP", Sev::inform) << "Add NormalDisplacementBc val " << val;
  MOFEM_LOG("EP", Sev::inform)
      << "Add NormalDisplacementBc nb. of faces " << faces.size();
}
 
PressureBc::PressureBc(std::string name, std::vector<double> attr, Range faces)
    : blockName(name), faces(faces) {
 
  blockName = name;
  if (attr.size() < 1) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong size of normal displacement BC");
  }
 
  val = attr[0];
 
  MOFEM_LOG("EP", Sev::inform) << "Add PressureBc " << name;
  MOFEM_LOG("EP", Sev::inform) << "Add PressureBc val " << val;
  MOFEM_LOG("EP", Sev::inform)
      << "Add PressureBc nb. of faces " << faces.size();
}
 
ExternalStrain::ExternalStrain(std::string name, std::vector<double> attr,
                               Range ents)
    : blockName(name), ents(ents) {
 
  blockName = name;
  if (attr.size() < 2) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong size of external strain attribute");
  }
 
  val = attr[0];
  bulkModulusK = attr[1];
 
  MOFEM_LOG("EP", Sev::inform) << "Add ExternalStrain " << name;
  MOFEM_LOG("EP", Sev::inform) << "Add ExternalStrain val " << val;
  MOFEM_LOG("EP", Sev::inform)
      << "Add ExternalStrain bulk modulus K " << bulkModulusK;
  MOFEM_LOG("EP", Sev::inform)
      << "Add ExternalStrain nb. of tets " << ents.size();
}
 
AnalyticalDisplacementBc::AnalyticalDisplacementBc(std::string name,
                                                   std::vector<double> attr,
                                                   Range faces)
    : blockName(name), faces(faces) {
 
  blockName = name;
  if (attr.size() < 3) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong size of analytical displacement BC");
  }
 
  flags.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    flags[ii] = attr[ii];
  }
 
  MOFEM_LOG("EP", Sev::inform) << "Add AnalyticalDisplacementBc " << name;
  MOFEM_LOG("EP", Sev::inform)
      << "Add AnalyticalDisplacementBc flags " << flags[0] << " " << flags[1]
      << " " << flags[2];
  MOFEM_LOG("EP", Sev::inform)
      << "Add AnalyticalDisplacementBc nb. of faces " << faces.size();
}
 
AnalyticalTractionBc::AnalyticalTractionBc(std::string name,
                                           std::vector<double> attr,
                                           Range faces)
    : blockName(name), faces(faces) {
 
  blockName = name;
  if (attr.size() < 3) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong size of analytical traction BC");
  }
 
  flags.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    flags[ii] = attr[ii];
  }
 
  MOFEM_LOG("EP", Sev::inform) << "Add AnalyticalTractionBc " << name;
  MOFEM_LOG("EP", Sev::inform) << "Add AnalyticalTractionBc flags " << flags[0]
                               << " " << flags[1] << " " << flags[2];
  MOFEM_LOG("EP", Sev::inform)
      << "Add AnalyticalTractionBc nb. of faces " << faces.size();
}
 
MoFEMErrorCode
EshelbianCore::getTractionFreeBc(const EntityHandle meshset,
                                 boost::shared_ptr<TractionFreeBc> &bc_ptr,
                                 const std::string contact_set_name) {
  MoFEMFunctionBegin;
 
  // get skin from all tets
  Range tets;
  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
  Range tets_skin_part;
  Skinner skin(&mField.get_moab());
  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
  Range tets_skin;
  CHKERR pcomm->filter_pstatus(tets_skin_part,
                               PSTATUS_SHARED | PSTATUS_MULTISHARED,
                               PSTATUS_NOT, -1, &tets_skin);
 
  bc_ptr->resize(3);
  for (int dd = 0; dd != 3; ++dd)
    (*bc_ptr)[dd] = tets_skin;
 
  // Do not remove dofs on which traction is applied
  if (bcSpatialDispVecPtr)
    for (auto &v : *bcSpatialDispVecPtr) {
      if (v.flags[0])
        (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      if (v.flags[1])
        (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      if (v.flags[2])
        (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  // Do not remove dofs on which rotation is applied
  if (bcSpatialRotationVecPtr)
    for (auto &v : *bcSpatialRotationVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialNormalDisplacementVecPtr)
    for (auto &v : *bcSpatialNormalDisplacementVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialAnalyticalDisplacementVecPtr)
    for (auto &v : *bcSpatialAnalyticalDisplacementVecPtr) {
      if (v.flags[0])
        (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      if (v.flags[1])
        (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      if (v.flags[2])
        (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialTractionVecPtr)
    for (auto &v : *bcSpatialTractionVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialAnalyticalTractionVecPtr)
    for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialPressureVecPtr)
    for (auto &v : *bcSpatialPressureVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  // remove contact
  for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
           std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
    Range faces;
    CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                              true);
    (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
    (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
    (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
  }
 
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Set integration rule on element
 * \param order on row
 * \param order on column
 * \param order on data
 *
 * Use maximal oder on data in order to determine integration rule
 *
 */
struct VolRule {
  int operator()(int p_row, int p_col, int p_data) const {
    return 2 * p_data + 1;
  }
};
 
struct FaceRule {
  int operator()(int p_row, int p_col, int p_data) const {
    return 2 * (p_data + 1);
  }
};
 
MoFEMErrorCode EshelbianCore::setBaseVolumeElementOps(
    const int tag, const bool do_rhs, const bool do_lhs, const bool calc_rates,
    SmartPetscObj<Vec> var_vec,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe) {
  MoFEMFunctionBegin;
 
  auto bubble_cache =
      boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0, MatrixDouble());
  fe->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      fe->getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions, frontAdjEdges);
 
  // set integration rule
  fe->getRuleHook = [](int, int, int) { return -1; };
  auto vol_rule_lin = [](int o) { return 2 * o + 1; };
  auto vol_rule_no_lin = [](int o) { return 2 * o + (o - 1) + 1; };
  auto vol_rule = (SMALL_ROT > 0) ? vol_rule_lin : vol_rule_no_lin;
  fe->setRuleHook = SetIntegrationAtFrontVolume(frontVertices, frontAdjEdges,
                                                vol_rule, bubble_cache);
  // fe->getRuleHook = VolRule();
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
    dataAtPts->physicsPtr = physicalEquations;
  }
 
  // calculate fields values
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
      piolaStress, dataAtPts->getApproxPAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
      bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
      piolaStress, dataAtPts->getDivPAtPts()));
 
  if (noStretch) {
    fe->getOpPtrVector().push_back(
        physicalEquations->returnOpCalculateStretchFromStress(
            dataAtPts, physicalEquations));
  } else {
    fe->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValues<3>(
            stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
  }
 
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
 
  // H1 displacements
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
      spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
 
  // velocities
  if (calc_rates) {
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
        spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
    if (noStretch) {
    } else {
      fe->getOpPtrVector().push_back(
          new OpCalculateTensor2SymmetricFieldValuesDot<3>(
              stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
      fe->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradientDot<6, 3>(
              stretchTensor, dataAtPts->getGradLogStretchDotTensorAtPts(),
              MBTET));
    }
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
        rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradientDot<3, 3>(
        rotAxis, dataAtPts->getRotAxisGradDotAtPts(), MBTET));
 
    // acceleration
    if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
      fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
          spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
    }
  }
 
  // variations
  if (var_vec.use_count()) {
    fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
        piolaStress, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
        var_vec));
    fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
        bubbleField, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
        var_vec, MBMAXTYPE));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        rotAxis, dataAtPts->getVarRotAxisPts(), var_vec, MBTET));
    fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
        piolaStress, dataAtPts->getDivVarPiolaPts(), var_vec));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        spatialL2Disp, dataAtPts->getVarWL2Pts(), var_vec, MBTET));
 
    if (noStretch) {
      fe->getOpPtrVector().push_back(
          physicalEquations->returnOpCalculateVarStretchFromStress(
              dataAtPts, physicalEquations));
    } else {
      fe->getOpPtrVector().push_back(
          new OpCalculateTensor2SymmetricFieldValues<3>(
              stretchTensor, dataAtPts->getVarLogStreachPts(), var_vec, MBTET));
    }
  }
 
  // calculate other derived quantities
  fe->getOpPtrVector().push_back(new OpCalculateRotationAndSpatialGradient(
      dataAtPts, ((do_rhs || do_lhs) && calc_rates) ? alphaOmega : 0.));
 
  // evaluate integration points
  if (noStretch) {
  } else {
    fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
        tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setVolumeElementOps(
    const int tag, const bool add_elastic, const bool add_material,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
  MoFEMFunctionBegin;
 
  /** Contact requires that body is marked */
  auto get_body_range = [this](auto name, int dim) {
    std::map<int, Range> map;
 
    for (auto m_ptr :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % name).str()
 
                 ))
 
    ) {
      Range ents;
      CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
                                                            dim, ents, true),
                     "by dim");
      map[m_ptr->getMeshsetId()] = ents;
    }
 
    return map;
  };
 
  constexpr bool stab_tau_dot_variant = false;
 
  auto local_tau_sacale = boost::make_shared<double>(1.0);
  using BoundaryEle =
      PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
  using BdyEleOp = BoundaryEle::UserDataOperator;
  auto set_scale_op = new BdyEleOp(NOSPACE, BdyEleOp::OPSPACE);
  set_scale_op->doWorkRhsHook =
      [this,
       local_tau_sacale](DataOperator *raw_op_ptr, int side, EntityType type,
                         EntitiesFieldData::EntData &data) -> MoFEMErrorCode {
    MoFEMFunctionBegin;
    auto op_ptr = static_cast<BdyEleOp *>(raw_op_ptr);
    auto h = std::sqrt(op_ptr->getFTensor1Normal().l2());
    *local_tau_sacale = (alphaTau / h);
    MoFEMFunctionReturn(0);
  };
 
  auto not_interface_face = [this](FEMethod *fe_method_ptr) {
    auto ent = fe_method_ptr->getFEEntityHandle();
    if (
 
        (interfaceFaces->find(ent) != interfaceFaces->end())
 
        || (crackFaces->find(ent) != crackFaces->end())
 
    ) {
      return false;
    };
    return true;
  };
 
  // Right hand side
  fe_rhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
  CHKERR setBaseVolumeElementOps(tag, true, false, true, SmartPetscObj<Vec>(),
                                 fe_rhs);
 
  // elastic
  if (add_elastic) {
 
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialEquilibrium(spatialL2Disp, dataAtPts, alphaW, alphaRho));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialRotation(rotAxis, dataAtPts, alphaOmega));
    if (noStretch) {
      // do nothing - no stretch approximation
    } else {
      if (!internalStressTagName.empty()) {
        switch (internalStressInterpOrder) {
        case 0:
          fe_rhs->getOpPtrVector().push_back(
              new OpGetInternalStress<0>(dataAtPts, internalStressTagName));
          break;
        case 1:
          fe_rhs->getOpPtrVector().push_back(
              new OpGetInternalStress<1>(dataAtPts, internalStressTagName));
          break;
        default:
          SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
                  "Unsupported internal stress interpolation order %d",
                  internalStressInterpOrder);
        }
        // set default time scaling for interal stresses to constant
        TimeScale::ScalingFun def_scaling_fun = [](double time) { return 1; };
        auto ts_internal_stress =
            boost::make_shared<DynamicRelaxationTimeScale>(
                "internal_stress_history.txt", false, def_scaling_fun);
        if (internalStressVoigt) {
          fe_rhs->getOpPtrVector().push_back(
              new OpSpatialPhysicalInternalStress<true>(
                  stretchTensor, dataAtPts, ts_internal_stress));
        } else {
          fe_rhs->getOpPtrVector().push_back(
              new OpSpatialPhysicalInternalStress<false>(
                  stretchTensor, dataAtPts, ts_internal_stress));
        }
      }
      if (auto op = physicalEquations->returnOpSpatialPhysicalExternalStrain(
              stretchTensor, dataAtPts, externalStrainVecPtr, timeScaleMap)) {
        fe_rhs->getOpPtrVector().push_back(op);
      } else if (externalStrainVecPtr && !externalStrainVecPtr->empty()) {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
                "OpSpatialPhysicalExternalStrain not implemented for this "
                "material");
      }
 
      fe_rhs->getOpPtrVector().push_back(
          physicalEquations->returnOpSpatialPhysical(stretchTensor, dataAtPts,
                                                     alphaU));
    }
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyP(piolaStress, dataAtPts));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyBubble(bubbleField, dataAtPts));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyDivTerm(piolaStress, dataAtPts));
 
    auto set_hybridisation_rhs = [&](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_tau_stabilsation_rhs = [&](auto &pip, auto side_fe_name,
                                        auto hybrid_field) {
      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, side_fe_name, 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);
 
      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);
 
      // Add stabilization operator
      auto broken_disp_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpGetBrokenBaseSideData<SideEleOp>(spatialL2Disp,
                                                 broken_disp_data_ptr));
      auto disp_mat_ptr = boost::make_shared<MatrixDouble>();
      if (stab_tau_dot_variant) {
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValuesDot<SPACE_DIM>(spatialL2Disp,
                                                           disp_mat_ptr));
      } else {
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp,
                                                        disp_mat_ptr));
      }
      // Set diag fluxes on skeleton side
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpSetFlux<SideEleOp>(broken_disp_data_ptr, disp_mat_ptr));
 
      op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
      op_loop_skeleton_side->getOpPtrVector().push_back(set_scale_op);
 
      // Add stabilization Ugamma Ugamma skeleton
      auto hybrid_ptr = boost::make_shared<MatrixDouble>();
      if (stab_tau_dot_variant) {
        op_loop_skeleton_side->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValuesDot<SPACE_DIM>(hybrid_field,
                                                           hybrid_ptr));
      } else {
        op_loop_skeleton_side->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValues<SPACE_DIM>(hybrid_field,
                                                        hybrid_ptr));
      }
 
      // Diag u_gamma - u_gamma faces
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpBaseTimesVectorFace(
              hybrid_field, hybrid_ptr,
              [local_tau_sacale, broken_disp_data_ptr](double, double, double) {
                return broken_disp_data_ptr->size() * (*local_tau_sacale);
              }));
      // Diag L2 - L2 volumes
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpBrokenBaseTimesBrokenDisp(
              broken_disp_data_ptr, [local_tau_sacale](double, double, double) {
                return (*local_tau_sacale);
              }));
      // Off-diag Ugamma - L2
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpHybridBaseTimesBrokenDisp(
              hybrid_field, broken_disp_data_ptr,
              [local_tau_sacale](double, double, double) {
                return -(*local_tau_sacale);
              }));
      // Off-diag L2 - Ugamma
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpBrokenBaseTimesHybridDisp(
              broken_disp_data_ptr, hybrid_ptr,
              [local_tau_sacale](double, double, double) {
                return -(*local_tau_sacale);
              }));
 
      // Add skeleton to domain pipeline
      pip.push_back(op_loop_skeleton_side);
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_contact_rhs = [&](auto &pip) {
      return pushContactOpsRhs(*this, contactTreeRhs, pip);
    };
 
    auto set_cohesive_rhs = [&](auto &pip) {
      return pushCohesiveOpsRhs(
          *this,
          SetIntegrationAtFrontFace(frontVertices, frontAdjEdges,
                                    [](int p) { return 2 * (p + 1) + 1; }),
          interfaceFaces, pip);
    };
 
    CHKERR set_hybridisation_rhs(fe_rhs->getOpPtrVector());
    CHKERR set_contact_rhs(fe_rhs->getOpPtrVector());
    if (alphaTau > 0.0) {
      CHKERR set_tau_stabilsation_rhs(fe_rhs->getOpPtrVector(), skeletonElement,
                                      hybridSpatialDisp);
      CHKERR set_tau_stabilsation_rhs(fe_rhs->getOpPtrVector(), contactElement,
                                      contactDisp);
    }
    if (intefaceCrack == PETSC_TRUE) {
      CHKERR set_cohesive_rhs(fe_rhs->getOpPtrVector());
    }
 
    // Body forces
    using BodyNaturalBC =
        NaturalBC<VolumeElementForcesAndSourcesCore::UserDataOperator>::
            Assembly<PETSC>::LinearForm<GAUSS>;
    using OpBodyForce =
        BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
 
    auto body_time_scale =
        boost::make_shared<DynamicRelaxationTimeScale>("body_force.txt");
    CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
        fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {body_time_scale},
        "BODY_FORCE", Sev::inform);
  }
 
  // Left hand side
  fe_lhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
  CHKERR setBaseVolumeElementOps(tag, true, true, true, SmartPetscObj<Vec>(),
                                 fe_lhs);
 
  // elastic
  if (add_elastic) {
 
    if (noStretch) {
      fe_lhs->getOpPtrVector().push_back(
          new OpSpatialConsistency_dP_dP(piolaStress, piolaStress, dataAtPts));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dP(
          bubbleField, piolaStress, dataAtPts));
      fe_lhs->getOpPtrVector().push_back(
          new OpSpatialConsistency_dBubble_dBubble(bubbleField, bubbleField,
                                                   dataAtPts));
    } else {
      fe_lhs->getOpPtrVector().push_back(
          physicalEquations->returnOpSpatialPhysical_du_du(
              stretchTensor, stretchTensor, dataAtPts, alphaU));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
          stretchTensor, piolaStress, dataAtPts, true));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
          stretchTensor, bubbleField, dataAtPts, true));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
          stretchTensor, rotAxis, dataAtPts,
          symmetrySelector == SYMMETRIC ? true : false));
    }
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
        spatialL2Disp, piolaStress, dataAtPts, true));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
        spatialL2Disp, spatialL2Disp, dataAtPts, alphaW, alphaRho));
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_domega(
        piolaStress, rotAxis, dataAtPts,
        symmetrySelector == SYMMETRIC ? true : false));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
        bubbleField, rotAxis, dataAtPts,
        symmetrySelector == SYMMETRIC ? true : false));
 
    if (symmetrySelector > SYMMETRIC) {
      if (!noStretch) {
        fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_du(
            rotAxis, stretchTensor, dataAtPts, false));
      }
      fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
          rotAxis, piolaStress, dataAtPts, false));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
          rotAxis, bubbleField, dataAtPts, false));
    }
    fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_domega(
        rotAxis, rotAxis, dataAtPts, alphaOmega));
 
    auto set_hybridisation_lhs = [&](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 = [](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_tau_stabilsation_lhs = [&](auto &pip, auto side_fe_name,
                                        auto hybrid_field) {
      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, side_fe_name, SPACE_DIM - 1, Sev::noisy);
      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);
 
      // 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);
 
      auto broken_disp_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpGetBrokenBaseSideData<SideEleOp>(spatialL2Disp,
                                                 broken_disp_data_ptr));
      op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
      op_loop_skeleton_side->getOpPtrVector().push_back(set_scale_op);
 
      auto time_shift = [](const FEMethod *fe_ptr) { return fe_ptr->ts_a; };
 
      // Diag Ugamma-Ugamma skeleton
      op_loop_skeleton_side->getOpPtrVector().push_back(new OpMassVectorFace(
          hybrid_field, hybrid_field,
          [local_tau_sacale, broken_disp_data_ptr](double, double, double) {
            return broken_disp_data_ptr->size() * (*local_tau_sacale);
          }));
      if (stab_tau_dot_variant) {
        static_cast<OpMassVectorFace &>(
            op_loop_skeleton_side->getOpPtrVector().back())
            .feScalingFun = time_shift;
      }
      // Diag L2-L2 volumes
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpBrokenBaseBrokenBase(
              broken_disp_data_ptr, [local_tau_sacale](double, double, double) {
                return (*local_tau_sacale);
              }));
      if (stab_tau_dot_variant) {
        static_cast<OpBrokenBaseBrokenBase &>(
            op_loop_skeleton_side->getOpPtrVector().back())
            .feScalingFun = time_shift;
      }
      // Off-diag Ugamma - L2
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpHyrbridBaseBrokenBase(
              hybrid_field, broken_disp_data_ptr,
              [local_tau_sacale](double, double, double) {
                return -(*local_tau_sacale);
              },
              false, false));
      if (stab_tau_dot_variant) {
        static_cast<OpHyrbridBaseBrokenBase &>(
            op_loop_skeleton_side->getOpPtrVector().back())
            .feScalingFun = time_shift;
      }
 
      // Off-diag L2 - Ugamma
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpHyrbridBaseBrokenBase(
              hybrid_field, broken_disp_data_ptr,
              [local_tau_sacale](double, double, double) {
                return -(*local_tau_sacale);
              },
              true, true));
      if (stab_tau_dot_variant) {
        static_cast<OpHyrbridBaseBrokenBase &>(
            op_loop_skeleton_side->getOpPtrVector().back())
            .feScalingFun = time_shift;
      }
 
      pip.push_back(op_loop_skeleton_side);
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_contact_lhs = [&](auto &pip) {
      return pushContactOpsLhs(*this, contactTreeRhs, pip);
    };
 
    auto set_cohesive_lhs = [&](auto &pip) {
      return pushCohesiveOpsLhs(
          *this,
          SetIntegrationAtFrontFace(frontVertices, frontAdjEdges,
                                    [](int p) { return 2 * (p + 1) + 1; }),
          interfaceFaces, pip);
    };
 
    CHKERR set_hybridisation_lhs(fe_lhs->getOpPtrVector());
    CHKERR set_contact_lhs(fe_lhs->getOpPtrVector());
    if (alphaTau > 0.0) {
      CHKERR set_tau_stabilsation_lhs(fe_lhs->getOpPtrVector(), skeletonElement,
                                      hybridSpatialDisp);
      CHKERR set_tau_stabilsation_lhs(fe_lhs->getOpPtrVector(), contactElement,
                                      contactDisp);
    }
    if (intefaceCrack == PETSC_TRUE) {
      CHKERR set_cohesive_lhs(fe_lhs->getOpPtrVector());
    }
  }
 
  if (add_material) {
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setFaceElementOps(
    const bool add_elastic, const bool add_material,
    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
  MoFEMFunctionBegin;
 
  fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
  fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
 
  // set integration rule
  // fe_rhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
  // fe_lhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
  fe_rhs->getRuleHook = [](int, int, int) { return -1; };
  fe_lhs->getRuleHook = [](int, int, int) { return -1; };
  fe_rhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
  fe_lhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
 
  CHKERR
  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      fe_rhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
  CHKERR
  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      fe_lhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
 
  if (add_elastic) {
 
    auto get_broken_op_side = [this](auto &pip) {
      using EleOnSide =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
      using SideEleOp = EleOnSide::UserDataOperator;
      // Iterate over domain FEs adjacent to boundary.
      auto broken_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
      // Note: EleOnSide, i.e. uses on domain projected skeleton rule
      auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
          mField, elementVolumeName, SPACE_DIM, Sev::noisy);
      op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
          boost::make_shared<CGGUserPolynomialBase>();
      CHKERR
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
          materialH1Positions, frontAdjEdges);
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
      boost::shared_ptr<double> piola_scale_ptr(new double);
      op_loop_domain_side->getOpPtrVector().push_back(
          physicalEquations->returnOpSetScale(piola_scale_ptr,
                                              physicalEquations));
      auto 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));
      pip.push_back(op_loop_domain_side);
      return std::make_tuple(broken_data_ptr, piola_scale_ptr);
    };
 
    auto set_rhs = [&]() {
      MoFEMFunctionBegin;
 
      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, timeScaleMap));
      fe_rhs->getOpPtrVector().push_back(new OpAnalyticalDispBc(
          broken_data_ptr, bcSpatialAnalyticalDisplacementVecPtr,
          timeScaleMap));
      fe_rhs->getOpPtrVector().push_back(new OpRotationBc(
          broken_data_ptr, bcSpatialRotationVecPtr, timeScaleMap));
 
      fe_rhs->getOpPtrVector().push_back(
          new OpBrokenTractionBc(hybridSpatialDisp, bcSpatialTractionVecPtr,
                                 piola_scale_ptr, timeScaleMap));
      auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
      // if you push gradient of L2 base to physical element, it will not work.
      fe_rhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM - 1>(
              hybridSpatialDisp, hybrid_grad_ptr));
      fe_rhs->getOpPtrVector().push_back(new OpBrokenPressureBc(
          hybridSpatialDisp, bcSpatialPressureVecPtr, piola_scale_ptr,
          hybrid_grad_ptr, timeScaleMap));
      fe_rhs->getOpPtrVector().push_back(new OpBrokenAnalyticalTractionBc(
          hybridSpatialDisp, bcSpatialAnalyticalTractionVecPtr, piola_scale_ptr,
          timeScaleMap));
 
      auto hybrid_ptr = boost::make_shared<MatrixDouble>();
      fe_rhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(hybridSpatialDisp,
                                                      hybrid_ptr));
      fe_rhs->getOpPtrVector().push_back(new OpNormalDispRhsBc(
          hybridSpatialDisp, hybrid_ptr, broken_data_ptr,
          bcSpatialNormalDisplacementVecPtr, timeScaleMap));
 
      auto get_normal_disp_bc_faces = [&]() {
        auto faces =
            get_range_from_block(mField, "NORMAL_DISPLACEMENT", SPACE_DIM - 1);
        return boost::make_shared<Range>(faces);
      };
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using BdyEleOp = BoundaryEle::UserDataOperator;
      using OpC_dBroken = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
          GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
      fe_rhs->getOpPtrVector().push_back(new OpC_dBroken(
          broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0),
          get_normal_disp_bc_faces()));
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_lhs = [&]() {
      MoFEMFunctionBegin;
 
      auto [broken_data_ptr, piola_scale_ptr] =
          get_broken_op_side(fe_lhs->getOpPtrVector());
 
      fe_lhs->getOpPtrVector().push_back(new OpNormalDispLhsBc_dU(
          hybridSpatialDisp, bcSpatialNormalDisplacementVecPtr, timeScaleMap));
      fe_lhs->getOpPtrVector().push_back(new OpNormalDispLhsBc_dP(
          hybridSpatialDisp, broken_data_ptr, bcSpatialNormalDisplacementVecPtr,
          timeScaleMap));
 
      auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
      // if you push gradient of L2 base to physical element, it will not work.
      fe_lhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM - 1>(
              hybridSpatialDisp, hybrid_grad_ptr));
      fe_lhs->getOpPtrVector().push_back(new OpBrokenPressureBcLhs_dU(
          hybridSpatialDisp, bcSpatialPressureVecPtr, hybrid_grad_ptr,
          timeScaleMap));
 
      auto get_normal_disp_bc_faces = [&]() {
        auto faces =
            get_range_from_block(mField, "NORMAL_DISPLACEMENT", SPACE_DIM - 1);
        return boost::make_shared<Range>(faces);
      };
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using BdyEleOp = BoundaryEle::UserDataOperator;
      using OpC = FormsIntegrators<BdyEleOp>::Assembly<A>::BiLinearForm<
          GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
      fe_lhs->getOpPtrVector().push_back(new OpC(
          hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0),
          true, true, get_normal_disp_bc_faces()));
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_rhs();
    CHKERR set_lhs();
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setFaceInterfaceOps(
    const bool add_elastic, const bool add_material,
    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setContactElementRhsOps(
 
    boost::shared_ptr<ForcesAndSourcesCore> &fe_contact_tree
 
) {
  MoFEMFunctionBegin;
  fe_contact_tree = createContactDetectionFiniteElement(*this);
  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/CGGUserPolynomialBase.cpp"
#include "impl/EshelbianMonitor.cpp"
#include "impl/SetUpSchurImpl.cpp"
 
struct solve_elastic_setup {
 
  inline static auto setup(EshelbianCore *ep_ptr, TS ts, Vec x,
                           bool set_ts_monitor) {
 
#ifdef ENABLE_PYTHON_BINDING
    auto setup_sdf = [&]() { return setupContactSdf(); };
#endif
 
    auto setup_ts_monitor = [&]() {
      boost::shared_ptr<TsCtx> ts_ctx;
      CHK_THROW_MESSAGE(DMMoFEMGetTsCtx(ep_ptr->dmElastic, ts_ctx),
                        "get TS ctx");
      if (set_ts_monitor) {
        CHK_THROW_MESSAGE(
            TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULLPTR),
            "TS monitor set");
        auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*ep_ptr);
        ts_ctx->getLoopsMonitor().push_back(
            TsCtx::PairNameFEMethodPtr(ep_ptr->elementVolumeName, monitor_ptr));
      }
      MOFEM_LOG("EP", Sev::inform) << "TS monitor setup";
      return std::make_tuple(ts_ctx);
    };
 
    auto setup_snes_monitor = [&]() {
      MoFEMFunctionBeginHot;
      SNES snes;
      CHKERR TSGetSNES(ts, &snes);
      auto snes_ctx = getDMSnesCtx(ep_ptr->dmElastic);
      CHKERR SNESMonitorSet(snes,
                            (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
                                               void *))MoFEMSNESMonitorEnergy,
                            (void *)(snes_ctx.get()), PETSC_NULLPTR);
      MOFEM_LOG("EP", Sev::inform) << "SNES monitor setup";
      MoFEMFunctionReturnHot(0);
    };
 
    auto setup_snes_conergence_test = [&]() {
      MoFEMFunctionBegin;
 
      auto snes_convergence_test = [](SNES snes, PetscInt it, PetscReal xnorm,
                                      PetscReal snorm, PetscReal fnorm,
                                      SNESConvergedReason *reason, void *cctx) {
        MoFEMFunctionBegin;
        // EshelbianCore *ep_ptr = (EshelbianCore *)cctx;
        CHKERR SNESConvergedDefault(snes, it, xnorm, snorm, fnorm, reason,
                                    PETSC_NULLPTR);
 
        Vec x_update, r;
        CHKERR SNESGetSolutionUpdate(snes, &x_update);
        CHKERR SNESGetFunction(snes, &r, PETSC_NULLPTR, PETSC_NULLPTR);
 
        MoFEMFunctionReturn(0);
      };
 
      // SNES snes;
      // CHKERR TSGetSNES(ts, &snes);
      // CHKERR SNESSetConvergenceTest(snes, snes_convergence_test, ep_ptr,
      //                               PETSC_NULLPTR);
      // MOFEM_LOG("EP", Sev::inform) << "SNES convergence test setup";
      MoFEMFunctionReturn(0);
    };
 
    auto setup_section = [&]() {
      PetscSection section_raw;
      CHK_THROW_MESSAGE(DMGetSection(ep_ptr->dmElastic, &section_raw),
                        "get DM section");
      int num_fields;
      CHK_THROW_MESSAGE(PetscSectionGetNumFields(section_raw, &num_fields),
                        "get num fields");
      for (int ff = 0; ff != num_fields; ff++) {
        const char *field_name;
        CHK_THROW_MESSAGE(
            PetscSectionGetFieldName(section_raw, ff, &field_name),
            "get field name");
        MOFEM_LOG_C("EP", Sev::inform, "Field %d name %s", ff, field_name);
      }
      return SmartPetscObj<PetscSection>(section_raw, true);
    };
 
    auto set_vector_on_mesh = [&]() {
      MoFEMFunctionBeginHot;
      CHKERR DMoFEMMeshToLocalVector(ep_ptr->dmElastic, x, INSERT_VALUES,
                                     SCATTER_FORWARD);
      CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
      MOFEM_LOG("EP", Sev::inform) << "Vector set on mesh";
      MoFEMFunctionReturnHot(0);
    };
 
    auto setup_schur_block_solver = [&]() {
      MOFEM_LOG("EP", Sev::inform) << "Setting up Schur block solver";
      CHK_THROW_MESSAGE(TSAppendOptionsPrefix(ts, "elastic_"),
                        "append options prefix");
      CHK_THROW_MESSAGE(TSSetFromOptions(ts), "set from options");
      CHK_THROW_MESSAGE(TSSetDM(ts, ep_ptr->dmElastic), "set DM");
      // Adding field split solver
      boost::shared_ptr<EshelbianCore::SetUpSchur> schur_ptr;
      if constexpr (A == AssemblyType::BLOCK_MAT) {
        schur_ptr =
            EshelbianCore::SetUpSchur::createSetUpSchur(ep_ptr->mField, ep_ptr);
        CHK_THROW_MESSAGE(schur_ptr->setUp(ts), "setup schur");
      }
      MOFEM_LOG("EP", Sev::inform) << "Setting up Schur block solver done";
      return schur_ptr;
    };
 
    // Warning: sequence of construction is not guaranteed for tuple. You have
    // to enforce order by proper packaging.
 
#ifdef ENABLE_PYTHON_BINDING
    return std::make_tuple(setup_sdf(), setup_ts_monitor(),
                           setup_snes_monitor(), setup_snes_conergence_test(),
                           setup_section(), set_vector_on_mesh(),
                           setup_schur_block_solver());
#else
    return std::make_tuple(setup_ts_monitor(), setup_snes_monitor(),
                           setup_snes_conergence_test(), setup_section(),
                           set_vector_on_mesh(), setup_schur_block_solver());
#endif
  }
};
 
MoFEMErrorCode EshelbianCore::addDebugModel(TS ts) {
  MoFEMFunctionBegin;
 
  PetscBool debug_model = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-debug_model", &debug_model,
                             PETSC_NULLPTR);
  MOFEM_LOG("EP", Sev::inform)
      << "Debug model flag is " << (debug_model ? "ON" : "OFF");
 
  if (debug_model == PETSC_TRUE) {
    auto ts_ctx_ptr = getDMTsCtx(dmElastic);
    auto post_proc = [&](TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, Vec F,
                         void *ctx) {
      MoFEMFunctionBegin;
 
      SNES snes;
      CHKERR TSGetSNES(ts, &snes);
      int it;
      CHKERR SNESGetIterationNumber(snes, &it);
      std::string file_name = "snes_iteration_" + std::to_string(it) + ".h5m";
      CHKERR postProcessResults(1, file_name, F, u_t);
      std::string file_skel_name =
          "snes_iteration_skel_" + std::to_string(it) + ".h5m";
 
      auto get_material_force_tag = [&]() {
        auto &moab = mField.get_moab();
        Tag tag;
        CHK_MOAB_THROW(moab.tag_get_handle("MaterialForce", tag),
                       "can't get tag");
        return tag;
      };
 
      CHKERR calculateFaceMaterialForce(1, ts);
      CHKERR postProcessSkeletonResults(1, file_skel_name, F,
                                        {get_material_force_tag()});
 
      MoFEMFunctionReturn(0);
    };
    ts_ctx_ptr->tsDebugHook = post_proc;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::solveElastic(TS ts, Vec x) {
  MoFEMFunctionBegin;
 
  CHKERR addDebugModel(ts);
 
  auto storage = solve_elastic_setup::setup(this, ts, x, true);
 
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
    Vec xx;
    CHKERR VecDuplicate(x, &xx);
    CHKERR VecZeroEntries(xx);
    CHKERR TS2SetSolution(ts, x, xx);
    CHKERR VecDestroy(&xx);
  } else {
    CHKERR TSSetSolution(ts, x);
  }
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
  CHKERR TSSetUp(ts);
  CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
  CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  CHKERR TSElasticPostStep::postStepInitialise(this);
  CHKERR TSSolve(ts, PETSC_NULLPTR);
  CHKERR TSElasticPostStep::postStepDestroy();
  TetPolynomialBase::switchCacheBaseOff<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
 
#ifndef NDEBUG
  // Make graph
  if (mField.get_comm_rank() == 0) {
    auto ts_ctx_ptr = getDMTsCtx(dmElastic);
    CHKERR PipelineGraph::writeTSGraphGraphviz(ts_ctx_ptr.get(),
                                               "solve_elastic_graph.dot");
  }
#endif
 
  SNES snes;
  CHKERR TSGetSNES(ts, &snes);
  int lin_solver_iterations;
  CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
  MOFEM_LOG("EP", Sev::inform)
      << "Number of linear solver iterations " << lin_solver_iterations;
 
  PetscBool test_cook_flg = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_cook", &test_cook_flg,
                             PETSC_NULLPTR);
  if (test_cook_flg) {
    constexpr int expected_lin_solver_iterations = 11;
    if (lin_solver_iterations > expected_lin_solver_iterations)
      SETERRQ(
          PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
          "Expected number of iterations is different than expected %d > %d",
          lin_solver_iterations, expected_lin_solver_iterations);
  }
 
  PetscBool test_sslv116_flag = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_sslv116",
                             &test_sslv116_flag, PETSC_NULLPTR);
 
  if (test_sslv116_flag) {
    double max_val = 0.0;
    double min_val = 0.0;
    auto field_min_max = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
      MoFEMFunctionBeginHot;
      auto ent_type = ent_ptr->getEntType();
      if (ent_type == MBVERTEX) {
        max_val = std::max(ent_ptr->getEntFieldData()[SPACE_DIM - 1], max_val);
        min_val = std::min(ent_ptr->getEntFieldData()[SPACE_DIM - 1], min_val);
      }
      MoFEMFunctionReturnHot(0);
    };
    CHKERR mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
        field_min_max, spatialH1Disp);
 
    double global_max_val = 0.0;
    double global_min_val = 0.0;
    MPI_Allreduce(&max_val, &global_max_val, 1, MPI_DOUBLE, MPI_MAX,
                  mField.get_comm());
    MPI_Allreduce(&min_val, &global_min_val, 1, MPI_DOUBLE, MPI_MIN,
                  mField.get_comm());
    MOFEM_LOG("EP", Sev::inform)
        << "Max " << spatialH1Disp << " value: " << global_max_val;
    MOFEM_LOG("EP", Sev::inform)
        << "Min " << spatialH1Disp << " value: " << global_min_val;
 
    double ref_max_val = 0.00767;
    double ref_min_val = -0.00329;
    if (std::abs(global_max_val - ref_max_val) > 1e-5) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Incorrect max value of the displacement field: %f != %f",
              global_max_val, ref_max_val);
    }
    if (std::abs(global_min_val - ref_min_val) > 4e-5) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Incorrect min value of the displacement field: %f != %f",
              global_min_val, ref_min_val);
    }
  }
 
  CHKERR gettingNorms();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::solveDynamicRelaxation(TS ts, Vec x,
                                                     int start_step,
                                                     double start_time) {
  MoFEMFunctionBegin;
 
  auto storage = solve_elastic_setup::setup(this, ts, x, false);
 
  double final_time = 1;
  double delta_time = 0.1;
  int max_it = 10;
  PetscBool ts_h1_update = PETSC_FALSE;
 
  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options", "none");
 
  CHKERR PetscOptionsScalar("-dynamic_final_time",
                            "dynamic relaxation final time", "", final_time,
                            &final_time, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-dynamic_delta_time",
                            "dynamic relaxation final time", "", delta_time,
                            &delta_time, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-dynamic_max_it", "dynamic relaxation iterations", "",
                         max_it, &max_it, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-dynamic_h1_update", "update each ts step", "",
                          ts_h1_update, &ts_h1_update, PETSC_NULLPTR);
 
  PetscOptionsEnd();
 
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation final time -dynamic_final_time = " << final_time;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation delta time -dynamic_delta_time = " << delta_time;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation max iterations -dynamic_max_it = " << max_it;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation H1 update each step -dynamic_h1_update = "
      << (ts_h1_update ? "TRUE" : "FALSE");
 
  CHKERR addDebugModel(ts);
 
  auto setup_ts_monitor = [&]() {
    auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
    return monitor_ptr;
  };
  auto monitor_ptr = setup_ts_monitor();
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
  CHKERR TSSetUp(ts);
  CHKERR TSElasticPostStep::postStepInitialise(this);
 
  double ts_delta_time;
  CHKERR TSGetTimeStep(ts, &ts_delta_time);
 
  if (ts_h1_update) {
    CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
    CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  }
 
  CHKERR TSElasticPostStep::preStepFun(ts);
  CHKERR TSElasticPostStep::postStepFun(ts);
 
  dynamicTime = start_time;
  dynamicStep = start_step;
  monitor_ptr->ts_u = PETSC_NULLPTR;
  monitor_ptr->ts_t = dynamicTime;
  monitor_ptr->ts_step = dynamicStep;
  CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
  for (; dynamicTime < final_time; dynamicTime += delta_time) {
    MOFEM_LOG("EP", Sev::inform) << "Load step " << 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 TSSetSolution(ts, x);
    CHKERR TSSolve(ts, PETSC_NULLPTR);
    if (!ts_h1_update) {
      CHKERR TSElasticPostStep::postStepFun(ts);
    }
 
    CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES,
                                   SCATTER_FORWARD);
    CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
 
    monitor_ptr->ts_u = x;
    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;
        MOFEM_LOG("EP", Sev::inform)
            << "Block " << name << " id " << bc->getMeshsetId() << " has "
            << r.size() << " entities";
      }
      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);
      MOFEM_LOG("EP", Sev::inform)
          << "Skin for block " << map.first << " id " << m.first << " has "
          << m.second.size() << " entities";
    }
    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) {
    auto get_tag = [&](auto name, auto dim) {
      auto &mob = mField.get_moab();
      Tag tag;
      double def_val[] = {0., 0., 0.};
      CHK_MOAB_THROW(mob.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
                                        MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
                     "create tag");
      return tag;
    };
    tags_to_transfer.push_back(get_tag("MaterialForce", 3));
  }
 
  {
    auto get_crack_tag = [&]() {
      Tag th;
      rval = mField.get_moab().tag_get_handle("CRACK", th);
      if (rval == MB_SUCCESS) {
        MOAB_THROW(mField.get_moab().tag_delete(th));
      }
      int def_val[] = {0};
      MOAB_THROW(mField.get_moab().tag_get_handle(
          "CRACK", 1, MB_TYPE_INTEGER, th, MB_TAG_SPARSE | MB_TAG_CREAT,
          def_val));
      return th;
    };
 
    Tag th = get_crack_tag();
    tags_to_transfer.push_back(th);
    int mark[] = {1};
    Range mark_faces;
    if (crackFaces)
      mark_faces.merge(*crackFaces);
    if (interfaceFaces)
      mark_faces.merge(*interfaceFaces);
    CHKERR mField.get_moab().tag_clear_data(th, mark_faces, mark);
  }
 
  // 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 != 0) {
            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["AngularMomentum"] = dataAtPts->getLeviKirchhoffAtPts();
      vec_fields["X"] = dataAtPts->getLargeXH1AtPts();
      if (!noStretch) {
        vec_fields["EiegnLogStreach"] = dataAtPts->getEigenValsAtPts();
      }
      if (var_vector) {
        auto vec = SmartPetscObj<Vec>(var_vector, true);
        vec_fields["VarOmega"] = dataAtPts->getVarRotAxisPts();
        vec_fields["VarSpatialDisplacementL2"] =
            boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            spatialL2Disp, vec_fields["VarSpatialDisplacementL2"], vec, MBTET));
      }
      if (f_residual) {
        auto vec = SmartPetscObj<Vec>(f_residual, true);
        vec_fields["ResSpatialDisplacementL2"] =
            boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            spatialL2Disp, vec_fields["ResSpatialDisplacementL2"], vec, MBTET));
        vec_fields["ResOmega"] = boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            rotAxis, vec_fields["ResOmega"], vec, MBTET));
      }
 
      OpPPMap::DataMapMat mat_fields;
      mat_fields["PiolaStress"] = dataAtPts->getApproxPAtPts();
      if (var_vector) {
        mat_fields["VarPiolaStress"] = dataAtPts->getVarPiolaPts();
      }
      if (f_residual) {
        auto vec = SmartPetscObj<Vec>(f_residual, true);
        mat_fields["ResPiolaStress"] = boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
            piolaStress, mat_fields["ResPiolaStress"],
            boost::make_shared<double>(1), vec));
        fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
            bubbleField, mat_fields["ResPiolaStress"],
            boost::make_shared<double>(1), vec, MBMAXTYPE));
      }
      if (!internalStressTagName.empty()) {
        mat_fields[internalStressTagName] = dataAtPts->getInternalStressAtPts();
        switch (internalStressInterpOrder) {
        case 0:
          fe.getOpPtrVector().push_back(
              new OpGetInternalStress<0>(dataAtPts, internalStressTagName));
          break;
        case 1:
          fe.getOpPtrVector().push_back(
              new OpGetInternalStress<1>(dataAtPts, internalStressTagName));
          break;
        default:
          SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
                  "Unsupported internal stress interpolation order %d",
                  internalStressInterpOrder);
        }
      }
 
      OpPPMap::DataMapMat mat_fields_symm;
      mat_fields_symm["LogSpatialStretch"] =
          dataAtPts->getLogStretchTensorAtPts();
      mat_fields_symm["SpatialStretch"] = dataAtPts->getStretchTensorAtPts();
      if (var_vector) {
        mat_fields_symm["VarLogSpatialStretch"] =
            dataAtPts->getVarLogStreachPts();
      }
      if (f_residual) {
        auto vec = SmartPetscObj<Vec>(f_residual, true);
        if (!noStretch) {
          mat_fields_symm["ResLogSpatialStretch"] =
              boost::make_shared<MatrixDouble>();
          fe.getOpPtrVector().push_back(
              new OpCalculateTensor2SymmetricFieldValues<3>(
                  stretchTensor, mat_fields_symm["ResLogSpatialStretch"], vec,
                  MBTET));
        }
      }
 
      fe.getOpPtrVector().push_back(
 
          new OpSidePPMap(
 
              post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
              {},
 
              vec_fields,
 
              mat_fields,
 
              mat_fields_symm,
 
              sense
 
              )
 
      );
 
      fe.getOpPtrVector().push_back(new OpPostProcDataStructure(
          post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
          dataAtPts, sense));
 
      MoFEMFunctionReturn(0);
    };
 
    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;
    // 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);
    auto traction_ptr = boost::make_shared<MatrixDouble>();
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
            piolaStress, traction_ptr, boost::make_shared<double>(1.0)));
 
    pip.push_back(new OpCalculateVectorFieldValues<3>(
        contactDisp, dataAtPts->getContactL2AtPts()));
    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(getOpContactDetection(
        *this, contactTreeRhs, u_h1_ptr, traction_ptr,
        get_range_from_block(mField, "CONTACT", SPACE_DIM - 1),
        &post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
    using BoundaryEle =
        PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
    auto op_this = new OpLoopThis<BoundaryEle>(mField, contactElement);
    op_this->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
            {},
 
            {{"ContactDisplacement", dataAtPts->getContactL2AtPts()}},
 
            {},
 
            {}
 
            )
 
    );
 
    if (f_residual) {
 
      pip.push_back(op_this);
      auto contact_residual = boost::make_shared<MatrixDouble>();
      op_this->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(
              contactDisp, contact_residual,
              SmartPetscObj<Vec>(f_residual, true)));
      op_this->getOpPtrVector().push_back(
 
          new OpPPMap(
 
              post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
              {},
 
              {{"res_contact", contact_residual}},
 
              {},
 
              {}
 
              )
 
      );
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto post_proc_mesh = boost::make_shared<moab::Core>();
  auto post_proc_ptr = get_post_proc(post_proc_mesh, /*positive sense*/ 1);
  auto post_proc_negative_sense_ptr =
      get_post_proc(post_proc_mesh, /*negative sense*/ -1);
  auto skin_post_proc_ptr = get_post_proc(post_proc_mesh, /*positive sense*/ 1);
  CHKERR calcs_side_traction_and_displacements(
      skin_post_proc_ptr, skin_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_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;
    };
    // this removes faces
    auto tets = get_adj(crack_faces, 3);
    // faces adjacent to tets not in crack_faces
    auto faces = subtract(get_adj(tets, 2), crack_faces);
    // what is left from below, are tets fully inside crack_faces
    tets = subtract(tets, get_adj(faces, 3));
    return subtract(crack_faces, get_adj(tets, 2));
  };
 
  auto side_one_faces = [&](auto &faces) {
    std::pair<Range, Range> sides;
    for (auto f : faces) {
      Range adj;
      MOAB_THROW(mField.get_moab().get_adjacencies(&f, 1, 3, false, adj));
      adj = intersect(own_tets, adj);
      for (auto t : adj) {
        int side, sense, offset;
        MOAB_THROW(mField.get_moab().side_number(t, f, side, sense, offset));
        if (sense == 1) {
          sides.first.insert(f);
        } else {
          sides.second.insert(f);
        }
      }
    }
    return sides;
  };
 
  auto crack_faces =
      unite(get_crack_faces(*crackFaces), get_crack_faces(*interfaceFaces));
  auto crack_side_faces = side_one_faces(crack_faces);
  auto side_one_crack_faces = [crack_side_faces](FEMethod *fe_method_ptr) {
    auto ent = fe_method_ptr->getFEEntityHandle();
    if (crack_side_faces.first.find(ent) == crack_side_faces.first.end()) {
      return false;
    }
    return true;
  };
  auto side_minus_crack_faces = [crack_side_faces](FEMethod *fe_method_ptr) {
    auto ent = fe_method_ptr->getFEEntityHandle();
    if (crack_side_faces.second.find(ent) == crack_side_faces.second.end()) {
      return false;
    }
    return true;
  };
 
  skin_post_proc_ptr->setTagsToTransfer(tags_to_transfer);
  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, skin_post_proc_ptr);
  post_proc_ptr->exeTestHook = side_one_crack_faces;
  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
      dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
  post_proc_negative_sense_ptr->exeTestHook = side_minus_crack_faces;
  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, skeletonElement,
                                              post_proc_negative_sense_ptr, 0,
                                              mField.get_comm_size());
 
  constexpr bool debug = false;
  if (debug) {
 
    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;
    };
 
    auto adj_front = filter_owners(mField, get_adj_front());
    auto only_front_faces = [adj_front](FEMethod *fe_method_ptr) {
      auto ent = fe_method_ptr->getFEEntityHandle();
      if (adj_front.find(ent) == adj_front.end()) {
        return false;
      }
      return true;
    };
 
    post_proc_ptr->exeTestHook = only_front_faces;
    CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
        dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
    post_proc_negative_sense_ptr->exeTestHook = only_front_faces;
    CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, skeletonElement,
                                                post_proc_negative_sense_ptr, 0,
                                                mField.get_comm_size());
  }
  auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
  CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
 
  CHKERR post_proc_end.writeFile(file.c_str());
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EshelbianCore::postProcessSkeletonResults(const int tag, const std::string file,
                                          Vec f_residual,
                                          std::vector<Tag> tags_to_transfer) {
  MoFEMFunctionBegin;
 
  using FaceEle = FaceElementForcesAndSourcesCore;
 
  auto post_proc_mesh = boost::make_shared<moab::Core>();
  auto post_proc_ptr =
      boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
          mField, post_proc_mesh);
  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM - 1, SPACE_DIM>::add(
      post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
      frontAdjEdges);
 
  auto hybrid_disp = boost::make_shared<MatrixDouble>();
  post_proc_ptr->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<3>(hybridSpatialDisp, hybrid_disp));
  post_proc_ptr->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
          hybridSpatialDisp, dataAtPts->getGradHybridDispAtPts()));
 
  auto op_loop_domain_side =
      new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
          mField, elementVolumeName, SPACE_DIM, Sev::noisy);
  post_proc_ptr->getOpPtrVector().push_back(op_loop_domain_side);
 
  // evaluated in side domain, that is op_loop_domain_side
  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>();
  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
      materialH1Positions, frontAdjEdges);
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
          piolaStress, dataAtPts->getApproxPAtPts()));
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpCalculateHTensorTensorField<SPACE_DIM, SPACE_DIM>(
          bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(
          rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<3>(
          spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
 
  if (noStretch) {
    op_loop_domain_side->getOpPtrVector().push_back(
        physicalEquations->returnOpCalculateStretchFromStress(
            dataAtPts, physicalEquations));
  } else {
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
            stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
  }
 
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
  OpPPMap::DataMapMat vec_fields;
  vec_fields["HybridDisplacement"] = hybrid_disp;
  // note that grad and omega have not trace, so this is only other side value
  vec_fields["spatialL2Disp"] = dataAtPts->getSmallWL2AtPts();
  vec_fields["Omega"] = dataAtPts->getRotAxisAtPts();
  OpPPMap::DataMapMat mat_fields;
  mat_fields["PiolaStress"] = dataAtPts->getApproxPAtPts();
  mat_fields["HybridDisplacementGradient"] =
      dataAtPts->getGradHybridDispAtPts();
  OpPPMap::DataMapMat mat_fields_symm;
  mat_fields_symm["LogSpatialStretch"] = dataAtPts->getLogStretchTensorAtPts();
 
  post_proc_ptr->getOpPtrVector().push_back(
 
      new OpPPMap(
 
          post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
          {},
 
          vec_fields,
 
          mat_fields,
 
          mat_fields_symm
 
          )
 
  );
 
  if (f_residual) {
    auto hybrid_res = boost::make_shared<MatrixDouble>();
    post_proc_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(
            hybridSpatialDisp, hybrid_res,
            SmartPetscObj<Vec>(f_residual, true)));
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
    post_proc_ptr->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
            {},
 
            {{"res_hybrid", hybrid_res}},
 
            {},
 
            {}
 
            )
 
    );
  }
 
  post_proc_ptr->setTagsToTransfer(tags_to_transfer);
 
  auto post_proc_begin =
      PostProcBrokenMeshInMoabBaseBegin(mField, post_proc_mesh);
  CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin.getFEMethod());
  CHKERR DMoFEMLoopFiniteElements(dM, skeletonElement, post_proc_ptr);
  auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
  CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
 
  CHKERR post_proc_end.writeFile(file.c_str());
 
  MoFEMFunctionReturn(0);
}
 
//! [Getting norms]
MoFEMErrorCode EshelbianCore::gettingNorms() {
  MoFEMFunctionBegin;
 
  auto post_proc_norm_fe =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
 
  auto post_proc_norm_rule_hook = [](int, int, int p) -> int {
    return 2 * (p);
  };
  post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
 
  post_proc_norm_fe->getUserPolynomialBase() =
      boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
 
  CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV}, materialH1Positions,
      frontAdjEdges);
 
  enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
  auto norms_vec =
      createVectorMPI(mField.get_comm(), LAST_NORM, PETSC_DETERMINE);
  CHKERR VecZeroEntries(norms_vec);
 
  auto u_l2_ptr = boost::make_shared<MatrixDouble>();
  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp, u_l2_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(spatialH1Disp, u_h1_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_NORM_L2));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_h1_ptr, norms_vec, U_NORM_H1));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_ERROR_L2,
                                         u_h1_ptr));
 
  auto piola_ptr = boost::make_shared<MatrixDouble>();
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateHVecTensorField<3, 3>(piolaStress, piola_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateHTensorTensorField<3, 3>(bubbleField, piola_ptr,
                                              MBMAXTYPE));
 
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor2<3, 3>(piola_ptr, norms_vec, PIOLA_NORM));
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>({post_proc_norm_fe.get()});
  CHKERR mField.loop_finite_elements("ELASTIC_PROBLEM", elementVolumeName,
                                     *post_proc_norm_fe);
  TetPolynomialBase::switchCacheBaseOff<HDIV>({post_proc_norm_fe.get()});
 
  CHKERR VecAssemblyBegin(norms_vec);
  CHKERR VecAssemblyEnd(norms_vec);
  const double *norms;
  CHKERR VecGetArrayRead(norms_vec, &norms);
  MOFEM_LOG("EP", Sev::inform) << "norm_u: " << std::sqrt(norms[U_NORM_L2]);
  MOFEM_LOG("EP", Sev::inform) << "norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
  MOFEM_LOG("EP", Sev::inform)
      << "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
  MOFEM_LOG("EP", Sev::inform)
      << "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
  CHKERR VecRestoreArrayRead(norms_vec, &norms);
 
  MoFEMFunctionReturn(0);
}
//! [Getting norms]
 
MoFEMErrorCode EshelbianCore::getSpatialDispBc() {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      "", piolaStress, false, false);
 
  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    if (auto disp_bc = bc.second->dispBcPtr) {
 
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      MOFEM_LOG("EP", Sev::noisy) << "Displacement BC: " << *disp_bc;
 
      std::vector<double> block_attributes(6, 0.);
      if (disp_bc->data.flag1 == 1) {
        block_attributes[0] = disp_bc->data.value1;
        block_attributes[3] = 1;
      }
      if (disp_bc->data.flag2 == 1) {
        block_attributes[1] = disp_bc->data.value2;
        block_attributes[4] = 1;
      }
      if (disp_bc->data.flag3 == 1) {
        block_attributes[2] = disp_bc->data.value3;
        block_attributes[5] = 1;
      }
      auto faces = bc.second->bcEnts.subset_by_dimension(2);
      bcSpatialDispVecPtr->emplace_back(block_name, block_attributes, faces);
    }
  }
  // old way of naming blocksets for displacement BCs
  CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
 
  bcSpatialNormalDisplacementVecPtr =
      boost::make_shared<NormalDisplacementBcVec>();
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)NORMAL_DISPLACEMENT(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialNormalDisplacementVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
  bcSpatialAnalyticalDisplacementVecPtr =
      boost::make_shared<AnalyticalDisplacementBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)ANALYTICAL_DISPLACEMENT(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialAnalyticalDisplacementVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
  auto ts_displacement =
      boost::make_shared<DynamicRelaxationTimeScale>("disp_history.txt");
  for (auto &bc : *bcSpatialDispVecPtr) {
    MOFEM_LOG("EP", Sev::noisy)
        << "Add time scaling displacement BC: " << bc.blockName;
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_displacement, "disp_history", ".txt", bc.blockName);
  }
 
  auto ts_normal_displacement =
      boost::make_shared<DynamicRelaxationTimeScale>("normal_disp_history.txt");
  for (auto &bc : *bcSpatialNormalDisplacementVecPtr) {
    MOFEM_LOG("EP", Sev::noisy)
        << "Add time scaling normal displacement BC: " << bc.blockName;
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_normal_displacement, "normal_disp_history", ".txt",
            bc.blockName);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::getSpatialTractionBc() {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<ForceCubitBcData>("", piolaStress,
                                                          false, false);
 
  bcSpatialTractionVecPtr = boost::make_shared<TractionBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    if (auto force_bc = bc.second->forceBcPtr) {
 
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      MOFEM_LOG("EP", Sev::noisy) << "Force BC: " << *force_bc;
 
      std::vector<double> block_attributes(6, 0.);
      block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
      block_attributes[3] = 1;
      block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
      block_attributes[4] = 1;
      block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
      block_attributes[5] = 1;
      auto faces = bc.second->bcEnts.subset_by_dimension(2);
      bcSpatialTractionVecPtr->emplace_back(block_name, block_attributes,
                                            faces);
    }
  }
  CHKERR getBc(bcSpatialTractionVecPtr, "SPATIAL_TRACTION_BC", 6);
 
  bcSpatialPressureVecPtr = boost::make_shared<PressureBcVec>();
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)PRESSURE(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
 
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialPressureVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
  bcSpatialAnalyticalTractionVecPtr =
      boost::make_shared<AnalyticalTractionBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)ANALYTICAL_TRACTION(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialAnalyticalTractionVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
  auto ts_traction =
      boost::make_shared<DynamicRelaxationTimeScale>("traction_history.txt");
  for (auto &bc : *bcSpatialTractionVecPtr) {
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_traction, "traction_history", ".txt", bc.blockName);
  }
 
  auto ts_pressure =
      boost::make_shared<DynamicRelaxationTimeScale>("pressure_history.txt");
  for (auto &bc : *bcSpatialPressureVecPtr) {
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_pressure, "pressure_history", ".txt", bc.blockName);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::getExternalStrain() {
  MoFEMFunctionBegin;
 
  auto getExternalStrain = [&](boost::shared_ptr<ExternalStrainVec>
                                   &ext_strain_vec_ptr,
                               const std::string block_name,
                               const int nb_attributes) {
    MoFEMFunctionBegin;
    for (auto it : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
             std::regex((boost::format("(.*)%s(.*)") % block_name).str()))) {
      std::vector<double> block_attributes;
      CHKERR it->getAttributes(block_attributes);
      if (block_attributes.size() < nb_attributes) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "In block %s expected %d attributes, but given %ld",
                it->getName().c_str(), nb_attributes, block_attributes.size());
      }
 
      auto get_block_ents = [&]() {
        Range ents;
        CHKERR mField.get_moab().get_entities_by_handle(it->meshset, ents,
                                                        true);
        return ents;
      };
      auto Ents = get_block_ents();
      ext_strain_vec_ptr->emplace_back(it->getName(), block_attributes,
                                       get_block_ents());
    }
    MoFEMFunctionReturn(0);
  };
 
  externalStrainVecPtr = boost::make_shared<ExternalStrainVec>();
 
  CHKERR getExternalStrain(externalStrainVecPtr, "EXTERNALSTRAIN", 2);
 
  auto ts_pre_stretch = boost::make_shared<DynamicRelaxationTimeScale>(
      "externalstrain_history.txt");
  for (auto &ext_strain_block : *externalStrainVecPtr) {
    MOFEM_LOG("EP", Sev::noisy)
        << "Add time scaling external strain: " << ext_strain_block.blockName;
    timeScaleMap[ext_strain_block.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_pre_stretch, "externalstrain_history", ".txt",
            ext_strain_block.blockName);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::createExchangeVectors(Sev sev) {
  MoFEMFunctionBegin;
 
  auto print_loc_size = [this](auto v, auto str, auto sev) {
    MoFEMFunctionBegin;
    int size;
    CHKERR VecGetLocalSize(v.second, &size);
    int low, high;
    CHKERR VecGetOwnershipRange(v.second, &low, &high);
    MOFEM_LOG("EPSYNC", sev) << str << " local size " << size << " ( " << low
                             << " " << high << " ) ";
    MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), sev);
    MoFEMFunctionReturn(0);
  };
 
  volumeExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 3, 1, sev);
  CHKERR print_loc_size(volumeExchange, "volumeExchange", sev);
  faceExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 2, 1, Sev::inform);
  CHKERR print_loc_size(faceExchange, "faceExchange", sev);
  edgeExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 1, 1, Sev::inform);
  CHKERR print_loc_size(edgeExchange, "edgeExchange", sev);
  vertexExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 0, 3, Sev::inform);
  CHKERR print_loc_size(vertexExchange, "vertexExchange", sev);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::solveCohesiveCrackGrowth(TS ts, Vec x,
                                                       int start_step,
                                                       double start_time) {
  MoFEMFunctionBegin;
 
  auto storage = solve_elastic_setup::setup(this, ts, x, false);
 
  auto cohesive_tao_ctx = createCohesiveTAOCtx(
      this,
      SetIntegrationAtFrontFace(frontVertices, frontAdjEdges,
                                [](int p) { return 2 * (p + 1) + 1; }),
      SmartPetscObj<TS>(ts, true));
 
  double final_time = 1;
  double delta_time = 0.1;
  int max_it = 10;
  PetscBool ts_h1_update = PETSC_FALSE;
 
  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options", "none");
 
  CHKERR PetscOptionsScalar("-dynamic_final_time",
                            "dynamic relaxation final time", "", final_time,
                            &final_time, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-dynamic_delta_time",
                            "dynamic relaxation final time", "", delta_time,
                            &delta_time, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-dynamic_max_it", "dynamic relaxation iterations", "",
                         max_it, &max_it, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-dynamic_h1_update", "update each ts step", "",
                          ts_h1_update, &ts_h1_update, PETSC_NULLPTR);
 
  PetscOptionsEnd();
 
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation final time -dynamic_final_time = " << final_time;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation delta time -dynamic_delta_time = " << delta_time;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation max iterations -dynamic_max_it = " << max_it;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation H1 update each step -dynamic_h1_update = "
      << (ts_h1_update ? "TRUE" : "FALSE");
 
  CHKERR initializeCohesiveKappaField(*this);
  CHKERR addDebugModel(ts);
 
  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);
 
  auto tao = createTao(mField.get_comm());
  CHKERR TaoSetType(tao, TAOLMVM);
  auto g = cohesive_tao_ctx->duplicateGradientVec();
  CHKERR TaoSetObjectiveAndGradient(tao, g,
                                    cohesiveEvaluateObjectiveAndGradient,
                                    (void *)cohesive_tao_ctx.get());
 
  dynamicTime = start_time;
  dynamicStep = start_step;
  monitor_ptr->ts_u = PETSC_NULLPTR;
  monitor_ptr->ts_t = dynamicTime;
  monitor_ptr->ts_step = dynamicStep;
  CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
  auto tao_sol0 = cohesive_tao_ctx->duplicateKappaVec();
  int tao_sol_size, tao_sol_loc_size;
  CHKERR VecGetSize(tao_sol0, &tao_sol_size);
  CHKERR VecGetLocalSize(tao_sol0, &tao_sol_loc_size);
  MOFEM_LOG("EP", Sev::inform)
      << "Cohesive crack growth initial kappa vector size " << tao_sol_size
      << " local size " << tao_sol_loc_size << " number of interface faces "
      << interfaceFaces->size();
 
  CHKERR TaoSetFromOptions(tao);
 
  auto xl = vectorDuplicate(tao_sol0);
  auto xu = vectorDuplicate(tao_sol0);
  CHKERR VecSet(xl, 0.0);
  CHKERR VecSet(xu, PETSC_INFINITY);
  CHKERR TaoSetVariableBounds(tao, xl, xu);
 
  for (; dynamicTime < final_time; dynamicTime += delta_time) {
    MOFEM_LOG("EP", Sev::inform) << "Load step " << dynamicStep << " Time "
                                 << dynamicTime << " delta time " << delta_time;
 
    CHKERR VecZeroEntries(tao_sol0);
    CHKERR VecGhostUpdateBegin(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR TaoSetSolution(tao, tao_sol0);
    CHKERR TaoSolve(tao);
    Vec tao_sol;
    CHKERR TaoGetSolution(tao, &tao_sol);
 
    // add solution increment to kappa vec/tags
    auto &kappa_vec = cohesive_tao_ctx->getKappaVec();
    CHKERR CommInterface::setVectorFromTag(mField.get_moab(), kappa_vec,
                                           get_kappa_tag(mField.get_moab()));
    CHKERR VecAXPY(kappa_vec.second, 1.0, tao_sol);
    CHKERR VecGhostUpdateBegin(kappa_vec.second, INSERT_VALUES,
                               SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(kappa_vec.second, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR CommInterface::setTagFromVector(mField.get_moab(), kappa_vec,
                                           get_kappa_tag(mField.get_moab()));
 
    CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES,
                                   SCATTER_FORWARD);
    CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
    monitor_ptr->ts_u = x;
    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);
}
 
} // namespace EshelbianPlasticity
 
#include <impl/EshelbianFracture.cpp>