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>
#include <IntegrationRules.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;
#endif
 
#include <EshelbianAux.hpp>
#include <EshelbianCohesive.hpp>
#include <EshelbianContact.hpp>
#include <EshelbianTopologicalDerivative.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_entities_by_handle(MoFEM::Interface &m_field,
                                   const std::string block_name) {
  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) {
    auto meshset = bc->getMeshset();
    CHK_MOAB_THROW(m_field.get_moab().get_entities_by_handle(meshset, r, true),
                   "get meshset ents");
  }
 
  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 {
 
auto vol_rule(int o) { return 2 * (o + 1); };
auto face_rule(int o) { return 2 * (o + 1); };
 
struct SetIntegrationAtFrontVolume {
 
  using FunRule = boost::function<int(int)>;
  FunRule funRule;
 
  SetIntegrationAtFrontVolume(
      boost::shared_ptr<Range> front_nodes,
      boost::shared_ptr<Range> front_edges,
      boost::shared_ptr<CGGUserPolynomialBase::CachePhi> cache_phi = nullptr)
      : funRule(vol_rule), 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)
      : funRule(fun_rule), frontNodes(front_nodes), frontEdges(front_edges),
        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;
      if (!funRule) {
        funRule = vol_rule;
      }
      const int rule = funRule(order_data);
      const auto xiao_rule =
          IntRules::XiaoGimbutas::getTetrahedronRule(rule);
      if (!xiao_rule) {
        SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "Xiao--Gimbutas tetrahedron rule is available for polynomial "
                "orders 0 to %d; requested %d",
                IntRules::XiaoGimbutas::tetrahedronRuleCount, rule);
      }
      if (xiao_rule->numBarycentricCoordinates != 4) {
        SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "wrong number of tetrahedron barycentric coordinates");
      }
 
      const size_t nb_gauss_pts = xiao_rule->numPoints;
      auto &gauss_pts = fe_ptr->gaussPts;
      gauss_pts.resize(4, nb_gauss_pts, false);
      cblas_dcopy(nb_gauss_pts, &xiao_rule->points[1], 4, &gauss_pts(0, 0),
                  1);
      cblas_dcopy(nb_gauss_pts, &xiao_rule->points[2], 4, &gauss_pts(1, 0),
                  1);
      cblas_dcopy(nb_gauss_pts, &xiao_rule->points[3], 4, &gauss_pts(2, 0),
                  1);
      cblas_dcopy(nb_gauss_pts, xiao_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, xiao_rule->points, 1, shape_ptr, 1);
      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 {
 
  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,
                            int (*)(int))
      : frontNodes(front_nodes), frontEdges(front_edges) {};
 
  MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                            int order_col, int order_data) {
    MoFEMFunctionBegin;
 
    constexpr bool debug = false;
 
    constexpr int numNodes = 3;
    constexpr int numEdges = 3;
    constexpr int refinementLevels = 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;
      const int rule = face_rule(order_data);
      const auto xiao_rule = IntRules::XiaoGimbutas::getTriangleRule(rule);
      if (!xiao_rule) {
        SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "Xiao--Gimbutas triangle rule is available for polynomial "
                "orders 0 to %d; requested %d",
                IntRules::XiaoGimbutas::triangleRuleCount, rule);
      }
      if (xiao_rule->numBarycentricCoordinates != 3) {
        SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "wrong number of triangle barycentric coordinates");
      }
 
      const size_t nb_gauss_pts = xiao_rule->numPoints;
      auto &gauss_pts = fe_ptr->gaussPts;
      gauss_pts.resize(3, nb_gauss_pts, false);
      cblas_dcopy(nb_gauss_pts, &xiao_rule->points[1], 3, &gauss_pts(0, 0),
                  1);
      cblas_dcopy(nb_gauss_pts, &xiao_rule->points[2], 3, &gauss_pts(1, 0),
                  1);
      cblas_dcopy(nb_gauss_pts, xiao_rule->weights, 1, &gauss_pts(2, 0), 1);
      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);
        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));
            }
 
            const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
            MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
            MatrixDouble shape_n(nb_gauss_pts, 3, false);
            CHKERR ShapeMBTRI(&shape_n(0, 0), &fe_ptr->gaussPts(0, 0),
                              &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
            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;
};
 
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_release[] = {"griffith_force", "griffith_skeleton"};
  const char *list_stretches[] = {"linear", "log", "log_quadratic"};
  const char *list_broken_hdiv_bases[] = {"demkowicz", "ainsworth"};
  PetscInt choice_rot = EshelbianCore::rotSelector;
  PetscInt choice_grad = EshelbianCore::gradApproximator;
  PetscInt choice_release = EshelbianCore::energyReleaseSelector;
  PetscInt choice_stretch = StretchSelector::LOG;
  PetscInt choice_solver = SolverType::TimeSolver;
  PetscInt choice_broken_hdiv_base = 0;
  PetscBool l2_user_base_scale_set = PETSC_FALSE;
  switch (EshelbianCore::brokenHdivBase) {
  case DEMKOWICZ_JACOBI_BASE:
    choice_broken_hdiv_base = 0;
    break;
  case AINSWORTH_LEGENDRE_BASE:
    choice_broken_hdiv_base = 1;
    break;
  default:
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
            "Unsupported broken HDIV base %s",
            ApproximationBaseNames[EshelbianCore::brokenHdivBase]);
  }
  char analytical_expr_file_name[255] = "analytical_expr.py";
  PetscBool no_stretch = isNoStretch() ? PETSC_TRUE : PETSC_FALSE;
 
  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",
                         "", materialH1Order, &materialH1Order, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
                            &alphaU, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
                            &alphaW, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-alpha_omega", "rot H1 penalty", "", alphaOmega,
                            &alphaOmega, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-alpha_omega0", "rot H1 penalty scaled by |M0|",
                            "", alphaOmega0, &alphaOmega0, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-alpha_r", "rot L2 penalty", "", alphaR, &alphaR,
                            PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-alpha_r0", "rot L2 penalty scaled by |M0|", "",
                            alphaR0, &alphaR0, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_omega", "rot viscosity", "",
                            alphaViscousOmega, &alphaViscousOmega,
                            PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_omega0",
                            "rot viscosity scaled by |M0|", "",
                            alphaViscousOmega0, &alphaViscousOmega0,
                            PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_r", "rot L2 viscosity", "",
                            alphaViscousR, &alphaViscousR, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_r0",
                            "rot L2 viscosity scaled by |M0|", "",
                            alphaViscousR0, &alphaViscousR0, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
                            &alphaRho, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-alpha_tau", "tau", "", alphaTau, &alphaTau,
                            PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-alpha_tau0", "tau0", "", alphaTau0, &alphaTau0,
                            PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-alpha_tau_bc_disp", "tau for displacement BC", "",
                            alphaTauBcDisp, &alphaTauBcDisp, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-alpha_tau0_bc_disp",
                            "tau0 for displacement BC", "", alphaTauBcDisp0,
                            &alphaTauBcDisp0, 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("-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", "",
                          no_stretch, &no_stretch, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-set_singularity", "set singularity", "",
                          setSingularity, &setSingularity, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-l2_user_base_scale", "streach scale", "",
                          l2UserBaseScale, &l2UserBaseScale,
                          &l2_user_base_scale_set);
  CHKERR PetscOptionsEList(
      "-broken_hdiv_base", "broken HDIV stress approximation base", "",
      list_broken_hdiv_bases, 2,
      list_broken_hdiv_bases[choice_broken_hdiv_base],
      &choice_broken_hdiv_base, PETSC_NULLPTR);
 
  // dynamic relaxation
 
  // @deprecate this option
  CHKERR PetscOptionsBool("-dynamic_relaxation", "dynamic time relaxation", "",
                          physicalTimeFlg, &physicalTimeFlg, PETSC_NULLPTR);
  CHKERR PetscOptionsEList(
      "-solver_type", "solver type", "", EshelbianCore::listSolvers,
      EshelbianCore::SolverType::LastSolver,
      EshelbianCore::listSolvers[choice_solver], &choice_solver, PETSC_NULLPTR);
 
  if (choice_solver != SolverType::TimeSolver) {
    CHKERR PetscOptionsScalar("-physical_final_time", "physical final time", "",
                              finalPhysicalTime,
                              &EshelbianCore::finalPhysicalTime, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-physical_delta_time", "physical delta time", "",
                              physicalDt, &EshelbianCore::physicalDt,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsInt("-physical_max_steps", "physical max iterations", "",
                           physicalMaxSteps, &EshelbianCore::physicalMaxSteps,
                           PETSC_NULLPTR);
    CHKERR PetscOptionsBool(
        "-physical_h1_update", "update each physicalsolver step", "",
        physicalH1Update, &EshelbianCore::physicalH1Update, PETSC_NULLPTR);
  }
 
  // contact parameters
  CHKERR PetscOptionsInt("-contact_max_post_proc_ref_level", "refinement level",
                         "", contactRefinementLevels, &contactRefinementLevels,
                         PETSC_NULLPTR);
  // cohesive interface
  CHKERR PetscOptionsBool("-cohesive_interface_on", "cohesive interface ON", "",
                          interfaceCrack, &interfaceCrack, PETSC_NULLPTR);
  CHKERR PetscOptionsInt(
      "-cohesive_interface_remove_level", "cohesive interface remove level", "",
      interfaceRemoveLevel, &interfaceRemoveLevel, PETSC_NULLPTR);
 
  // cracking parameters
  CHKERR PetscOptionsBool("-cracking_on", "cracking ON", "", crackingOn,
                          &crackingOn, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-cracking_add_time", "cracking add time", "",
                            crackingAddTime, &crackingAddTime, 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 PetscOptionsScalar("-cracking_rtol", "Cracking relative tolerance", "",
                            crackingRtol, &crackingRtol, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-cracking_atol", "Cracking absolute tolerance", "",
                            crackingAtol, &crackingAtol, 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 PetscOptionsBool("-internal_stress_voigt", "Voigt index notation", "",
                          internalStressVoigt, &internalStressVoigt,
                          PETSC_NULLPTR);
 
  // Heterogenous Young's modulus
  char tag_heterogeneous_youngs_modulus_name[255] = "";
  CHKERR PetscOptionsString(
      "-heterogeneous_youngs_modulus", "heterogeneous Young's modulus tag name",
      "", "", tag_heterogeneous_youngs_modulus_name, 255, PETSC_NULLPTR);
  heterogeneousYoungModTagName = string(tag_heterogeneous_youngs_modulus_name);
 
  PetscBool has_analytical_expr_file_option = PETSC_FALSE;
  CHKERR PetscOptionsGetString(
      PETSC_NULLPTR, PETSC_NULLPTR, "-analytical_expr_file",
      analytical_expr_file_name, 255, &has_analytical_expr_file_option);
  if (!has_analytical_expr_file_option) {
    const auto analytical_expr_script =
        mField.getInterface<JsonConfigManager>()->getPythonScriptByKey(
            "analytical_expr");
    if (!analytical_expr_script.empty()) {
      CHKERR PetscStrncpy(analytical_expr_file_name,
                          analytical_expr_script.c_str(),
                          sizeof(analytical_expr_file_name));
      MOFEM_LOG("EP", Sev::inform)
          << "Using Python script 'analytical_expr' from JSON config: "
          << analytical_expr_file_name;
    }
  }
 
  PetscOptionsEnd();
 
  alphaTauBcDisp = std::max(alphaTauBcDisp, alphaTau);
  alphaTauBcDisp0 = std::max(alphaTauBcDisp0, alphaTau0);
 
  PetscOptionsBegin(PETSC_COMM_WORLD, "mesh_transfer_", "mesh data transfer",
                    "none");
  char tag_mesh_transfer_source_file_name[255] = "";
  CHKERR PetscOptionsString("-source_file", "source mesh file name", "",
                            "source.h5m", tag_mesh_transfer_source_file_name,
                            255, &meshTransferSourceMeshFileSpecified);
  meshTransferSourceMeshFileName = string(tag_mesh_transfer_source_file_name);
  CHKERR PetscOptionsInt("-interp_order", "interpolation order", "", 0,
                         &meshTransferInterpOrder, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-hybrid_interp", "use hybrid interpolation", "",
                          meshTransferHybridInterp, &meshTransferHybridInterp,
                          PETSC_NULLPTR);
  PetscOptionsEnd();
 
  if (meshTransferInterpOrder != 0 && meshTransferInterpOrder != 1) {
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
            "Unsupported mesh transfer interpolation order %d",
            meshTransferInterpOrder);
  }
  if (!internalStressTagName.empty())
    listTagsToProject.push_back(internalStressTagName);
  if (!heterogeneousYoungModTagName.empty())
    listTagsToProject.push_back(heterogeneousYoungModTagName);
 
  const PetscBool l2_user_base_scale_option = l2UserBaseScale;
  if (setSingularity && !l2_user_base_scale_set) {
    l2UserBaseScale = PETSC_TRUE;
  }
 
  EshelbianCore::stretchHandling =
      no_stretch ? NO_STREACH_LINEAR : FULL_ORDER_STREACH;
  EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
  EshelbianCore::gradApproximator = static_cast<RotSelector>(choice_grad);
  EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
  EshelbianCore::energyReleaseSelector =
      static_cast<EnergyReleaseSelector>(choice_release);
  switch (choice_broken_hdiv_base) {
  case 0:
    EshelbianCore::brokenHdivBase = DEMKOWICZ_JACOBI_BASE;
    break;
  case 1:
    EshelbianCore::brokenHdivBase = AINSWORTH_LEGENDRE_BASE;
    break;
  default:
    SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
            "Unknown broken HDIV base option");
  }
 
  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:
    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_e;
    EshelbianCore::inv_d_f = inv_d_f_log_e;
    EshelbianCore::inv_dd_f = inv_dd_f_log_e;
    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 = inv_f_log_e_quadratic;
    EshelbianCore::inv_d_f = inv_d_f_log_e_quadratic;
    EshelbianCore::inv_dd_f = inv_dd_f_log_e_quadratic;
    break;
  default:
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
    break;
  };
 
  const PetscBool dynamic_relaxation_option = physicalTimeFlg;
  if (physicalTimeFlg) {
    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);
    physicalTimeFlg = PETSC_TRUE;
    break;
  case SolverType::Cohesive:
    EshelbianCore::solverType = static_cast<SolverType>(SolverType::Cohesive);
    break;
  case SolverType::LoadFactor:
    EshelbianCore::solverType = static_cast<SolverType>(SolverType::LoadFactor);
    physicalTimeFlg = PETSC_TRUE;
    break;
  case SolverType::ShapeOptimization:
    EshelbianCore::solverType =
        static_cast<SolverType>(SolverType::ShapeOptimization);
    break;
  case SolverType::TestTopologicalDerivative:
    EshelbianCore::solverType =
        static_cast<SolverType>(SolverType::TestTopologicalDerivative);
    physicalTimeFlg = PETSC_TRUE;
    break;
  default:
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown solver");
    break;
  };
 
  // start cracking after adding crack elements
  crackingStartTime = std::max(crackingStartTime, crackingAddTime);
 
  const auto yes_no = [](auto flag) { return flag ? "yes" : "no"; };
 
  MOFEM_LOG("EP", Sev::inform) << "spaceOrder: -space_order " << spaceOrder;
  MOFEM_LOG("EP", Sev::inform)
      << "spaceH1Order: -space_h1_order " << spaceH1Order;
  MOFEM_LOG("EP", Sev::inform)
      << "materialH1Order: -material_order " << materialH1Order;
  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: -alpha_omega " << alphaOmega;
  MOFEM_LOG("EP", Sev::inform)
      << "alphaOmega0: -alpha_omega0 " << alphaOmega0;
  MOFEM_LOG("EP", Sev::inform) << "alphaR: -alpha_r " << alphaR;
  MOFEM_LOG("EP", Sev::inform) << "alphaR0: -alpha_r0 " << alphaR0;
  MOFEM_LOG("EP", Sev::inform)
      << "alphaViscousOmega: -viscosity_alpha_omega "
      << alphaViscousOmega;
  MOFEM_LOG("EP", Sev::inform)
      << "alphaViscousOmega0: -viscosity_alpha_omega0 "
      << alphaViscousOmega0;
  MOFEM_LOG("EP", Sev::inform)
      << "alphaViscousR: -viscosity_alpha_r " << alphaViscousR;
  MOFEM_LOG("EP", Sev::inform)
      << "alphaViscousR0: -viscosity_alpha_r0 " << alphaViscousR0;
  MOFEM_LOG("EP", Sev::inform) << "alphaRho: -density_alpha_rho " << alphaRho;
  MOFEM_LOG("EP", Sev::inform) << "alphaTau: -alpha_tau " << alphaTau;
  MOFEM_LOG("EP", Sev::inform) << "alphaTau0: -alpha_tau0 " << alphaTau0;
  MOFEM_LOG("EP", Sev::inform)
      << "alphaTauBcDisp: -alpha_tau_bc_disp " << alphaTauBcDisp;
  MOFEM_LOG("EP", Sev::inform)
      << "alphaTauBcDisp0: -alpha_tau0_bc_disp " << alphaTauBcDisp0;
  MOFEM_LOG("EP", Sev::inform)
      << "Rotations: -rotations " << list_rots[EshelbianCore::rotSelector];
  MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation: -grad "
                               << list_rots[EshelbianCore::gradApproximator];
  MOFEM_LOG("EP", Sev::inform)
      << "Stretch: -stretches " << list_stretches[choice_stretch];
  MOFEM_LOG("EP", Sev::inform)
      << "No stretch: -no_stretch "
      << (isNoStretch() ? "yes" : "no");
 
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation: -dynamic_relaxation "
      << yes_no(dynamic_relaxation_option);
  MOFEM_LOG("EP", Sev::inform) << "Solver type: -solver_type "
                               << EshelbianCore::listSolvers[choice_solver];
  if (choice_solver != SolverType::TimeSolver) {
    MOFEM_LOG("EP", Sev::inform)
        << "Physical final time: -physical_final_time " << finalPhysicalTime;
    MOFEM_LOG("EP", Sev::inform)
        << "Physical delta time: -physical_delta_time " << physicalDt;
    MOFEM_LOG("EP", Sev::inform)
        << "Physical max steps: -physical_max_steps " << physicalMaxSteps;
    MOFEM_LOG("EP", Sev::inform)
        << "Physical H1 update: -physical_h1_update "
        << yes_no(physicalH1Update);
  }
  MOFEM_LOG("EP", Sev::inform)
      << "Singularity: -set_singularity " << yes_no(setSingularity);
  MOFEM_LOG("EP", Sev::inform)
      << "L2 user base scale: -l2_user_base_scale "
      << yes_no(l2_user_base_scale_option);
  if (l2UserBaseScale != l2_user_base_scale_option) {
    MOFEM_LOG("EP", Sev::inform)
        << "Effective L2 user base scale after option processing "
        << yes_no(l2UserBaseScale) << " (auto-enabled by -set_singularity)";
  }
  MOFEM_LOG("EP", Sev::inform)
      << "Broken HDIV base: -broken_hdiv_base "
      << list_broken_hdiv_bases[choice_broken_hdiv_base];
  MOFEM_LOG("EP", Sev::inform)
      << "Contact max post-proc ref level: -contact_max_post_proc_ref_level "
      << contactRefinementLevels;
 
  MOFEM_LOG("EP", Sev::inform)
      << "Cracking on: -cracking_on " << yes_no(crackingOn);
  MOFEM_LOG("EP", Sev::inform)
      << "Cracking add time: -cracking_add_time " << crackingAddTime;
  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)
      << "Cracking relative tolerance: -cracking_rtol " << crackingRtol;
  MOFEM_LOG("EP", Sev::inform)
      << "Cracking absolute tolerance: -cracking_atol " << crackingAtol;
  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 / "
         "-nb_J_integral_levels "
      << nbJIntegralContours;
  MOFEM_LOG("EP", Sev::inform)
      << "Cohesive interface on: -cohesive_interface_on "
      << ((interfaceCrack == PETSC_TRUE) ? "yes" : "no");
  MOFEM_LOG("EP", Sev::inform)
      << "Cohesive interface remove level: -cohesive_interface_remove_level "
      << interfaceRemoveLevel;
  MOFEM_LOG("EP", Sev::inform)
      << "Internal stress tag name: -internal_stress_tag_name "
      << internalStressTagName;
  MOFEM_LOG("EP", Sev::inform)
      << "Internal stress Voigt notation: -internal_stress_voigt "
      << yes_no(internalStressVoigt);
  MOFEM_LOG("EP", Sev::inform)
      << "Heterogeneous Young's modulus: -heterogeneous_youngs_modulus "
      << heterogeneousYoungModTagName;
  MOFEM_LOG("EP", Sev::inform)
      << "Analytical expression file: -analytical_expr_file "
      << analytical_expr_file_name;
  if (meshTransferSourceMeshFileSpecified) {
    MOFEM_LOG("EP", Sev::inform)
        << "Mesh transfer source file: -mesh_transfer_source_file "
        << meshTransferSourceMeshFileName;
  } else {
    MOFEM_LOG("EP", Sev::inform)
        << "Mesh transfer source file: -mesh_transfer_source_file <not set>";
  }
  MOFEM_LOG("EP", Sev::inform)
      << "Mesh transfer interpolation order: -mesh_transfer_interp_order "
      << meshTransferInterpOrder;
  MOFEM_LOG("EP", Sev::inform)
      << "Mesh transfer hybrid interpolation: -mesh_transfer_hybrid_interp "
      << yes_no(meshTransferHybridInterp);
 
#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,
                                        const bool add_bubble) {
  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 (bcSpatialSpringVecPtr)
      for (auto &v : *bcSpatialSpringVecPtr) {
        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
 
  const FieldApproximationBase broken_hdiv_base =
      EshelbianCore::brokenHdivBase;
 
  auto add_broken_hdiv_field = [this, meshset, broken_hdiv_base](
                                   const std::string field_name,
                                   const int order) {
    MoFEMFunctionBegin;
 
    const FieldApproximationBase base = broken_hdiv_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);
  };
 
  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));
    }
    if (bcSpatialSpringVecPtr)
      for (auto &bc : *bcSpatialSpringVecPtr) {
        faces.merge(intersect(bc.faces, skin));
      }
    return faces;
  };
 
  auto add_spatial_fields = [&]<FieldApproximationBase Base>() {
    MoFEMFunctionBegin;
    using Orders = EshelbianCore::FieldOrders<Base>;
    CHKERR add_broken_hdiv_field(piolaStress, Orders::stress(spaceOrder));
    if (add_bubble) {
      CHKERR add_bubble_field(bubbleField, Orders::bubble(spaceOrder), 1);
    }
    CHKERR add_l2_field(spatialL2Disp, Orders::disp(spaceOrder), 3);
    CHKERR add_user_l2_field(rotAxis, Orders::rot(spaceOrder), 3);
    CHKERR add_user_l2_field(stretchTensor,
                             !isNoStretch()
                                 ? Orders::stretch(spaceOrder)
                                 : -1,
                             6);
    CHKERR add_l2_field_by_range(hybridSpatialDisp,
                                 Orders::hybrid(spaceOrder), 2, 3,
                                 get_hybridised_disp());
    CHKERR add_l2_field_by_range(contactDisp, Orders::hybrid(spaceOrder), 2, 3,
                                 Range(*contactFaces));
    MoFEMFunctionReturn(0);
  };
 
  CHKERR withFieldOrders(add_spatial_fields);
 
  // spatial displacement
  CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
  // material positions
  CHKERR add_h1_field(materialH1Positions, materialH1Order, 3);
 
  CHKERR mField.build_fields();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::projectGeometry(const EntityHandle meshset,
                                              double time) {
  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));
  };
 
  if ((time - crackingAddTime) > std::numeric_limits<double>::epsilon()) {
    crackFaces = boost::make_shared<Range>(
        get_range_from_block(mField, "CRACK", SPACE_DIM - 1));
  } else {
    crackFaces = boost::make_shared<Range>();
  }
  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();
 
  auto get_interface_from_block = [&](auto block_name) {
    auto vol_eles = get_range_from_block(mField, block_name, SPACE_DIM);
    auto skin = filter_true_skin(mField, get_skin(mField, vol_eles));
    Range faces;
    CHKERR mField.get_moab().get_adjacencies(vol_eles, SPACE_DIM - 1, true,
                                             faces, moab::Interface::UNION);
    faces = subtract(faces, skin);
    MOFEM_LOG("EP", Sev::inform)
        << "Number of vol interface elements: " << vol_eles.size()
        << " and faces: " << faces.size();
    return faces;
  };
 
  interfaceFaces->merge(get_interface_from_block("VOLUME_INTERFACE"));
 
  auto remove_interface_from_block = [&](auto block_name, auto level) {
    MoFEMFunctionBegin;
    Range intreface_faces;
    if (mField.get_comm_rank() == 0) {
      auto ents = get_entities_by_handle(mField, block_name);
      for (auto l = 0; l < level; ++l) {
        Range adj_tets;
        CHKERR mField.get_moab().get_adjacencies(
            ents, SPACE_DIM, true, adj_tets, moab::Interface::UNION);
        Range adj_tets_faces;
        CHKERR mField.get_moab().get_adjacencies(adj_tets, SPACE_DIM - 1, true,
                                                 adj_tets_faces,
                                                 moab::Interface::UNION);
        ents.merge(adj_tets_faces);
      }
      auto faces = ents.subset_by_dimension(SPACE_DIM - 1);
      if (faces.size()) {
        MOFEM_LOG("EP", Sev::inform)
            << "Removed ents " << faces.size()
            << " interface faces: " << interfaceFaces->size();
      }
      intreface_faces = subtract(*interfaceFaces, faces);
      MOFEM_LOG("EP", Sev::noisy)
          << "Interface faces after remove " << intreface_faces;
    }
    auto intreface_faces_global = send_type(mField, intreface_faces, MBTRI);
    interfaceFaces->swap(intreface_faces_global);
    MoFEMFunctionReturn(0);
  };
  CHKERR remove_interface_from_block("REMOVE_INTERFACE", interfaceRemoveLevel);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::projectMaterialTags(const EntityHandle meshset) {
  MoFEMFunctionBegin;
#ifdef INCLUDE_MBCOUPLER
 
  double toler = 5.e-10;
  MOFEM_LOG_CHANNEL("WORLD");
  MOFEM_LOG_TAG("WORLD", "mesh_data_transfer");
  if (!meshTransferSourceMeshFileSpecified) {
    MOFEM_LOG("WORLD", Sev::verbose)
        << "No source mesh specified. Skipping projection";
    MoFEMFunctionReturnHot(0);
  }
  MOFEM_LOG("WORLD", Sev::inform)
      << "Projecting from source mesh: " << meshTransferSourceMeshFileName;
  MOFEM_LOG("WORLD", Sev::verbose)
      << "Interpolation Stress tag name: " << internalStressTagName;
  MOFEM_LOG("WORLD", Sev::verbose) << "Interpolation Young's modulus tag name: "
                                   << heterogeneousYoungModTagName;
  MOFEM_LOG("WORLD", Sev::verbose)
      << "Interpolation order: " << meshTransferInterpOrder;
  MOFEM_LOG("WORLD", Sev::verbose) << "Using hybrid interpolation: "
                                   << (meshTransferHybridInterp ? "yes" : "no");
 
  auto &moab = mField.get_moab();
 
  // check if tag exists
  for (const auto &tag_name : listTagsToProject) {
    Tag old_interp_tag;
    auto rval_check_tag = moab.tag_get_handle(tag_name.c_str(), old_interp_tag);
    if (rval_check_tag == MB_SUCCESS) {
      MOFEM_LOG("WORLD", Sev::inform)
          << "Deleting existing tag on target mesh: " << 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 (meshTransferInterpOrder) {
  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
  std::vector<Tag> interp_tags;
  std::vector<int> tag_length;
  std::vector<DataType> dtype;
  std::vector<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(
        meshTransferSourceMeshFileName.c_str(), &source_root, "");
    if (rval_source_mesh != MB_SUCCESS) {
      MOFEM_LOG("WORLD", Sev::warning) << "Error loading source mesh file: "
                                       << meshTransferSourceMeshFileName;
    }
    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);
 
    for (const auto &iterp_tag_name : listTagsToProject) {
      std::string tag_to_use = iterp_tag_name;
 
      Tag interp_tag;
      CHKERR moab.tag_get_handle(tag_to_use.c_str(), 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.push_back(interp_tag_len);
      dtype.push_back(DataType());
      storage.push_back(TagType());
      interp_tags.push_back(interp_tag);
      CHKERR moab.tag_get_data_type(interp_tag, dtype.back());
      CHKERR moab.tag_get_type(interp_tag, storage.back());
 
      // 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 (meshTransferInterpOrder == 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(tag_to_use) + "_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 (meshTransferInterpOrder == 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);
      };
 
      MOFEM_LOG("WORLD", Sev::inform)
          << "Performing interpolation for tag: " << tag_to_use;
      MOFEM_LOG("WORLD", Sev::inform)
          << "Number of target points to interpolate: " << num_pts;
      MOFEM_LOG("WORLD", Sev::inform)
          << "Interpolation method: "
          << (method == Coupler::CONSTANT ? "constant" : "linear FE");
      MOFEM_LOG("WORLD", Sev::inform)
          << "Number of components in tag: " << interp_tag_len;
 
      MOFEM_LOG("WORLD", Sev::inform)
          << "Source tag data range: ["
          << *std::min_element(source_data.begin(), source_data.end()) << ", "
          << *std::max_element(source_data.begin(), source_data.end()) << "]";
 
      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() && (meshTransferInterpOrder == 1) &&
          meshTransferHybridInterp) {
        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 after projection but keep the tags info for broadcast
    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
  int tag_size = tag_length.size();
  MPI_Bcast(&tag_size, 1, MPI_INT, 0, PETSC_COMM_WORLD);
  if (rank != 0) {
    interp_tags.resize(tag_size);
    tag_length.resize(tag_size);
    dtype.resize(tag_size);
    storage.resize(tag_size);
  }
  MPI_Bcast(interp_tags.data(), tag_size, MPI_INT, 0, PETSC_COMM_WORLD);
  MPI_Bcast(tag_length.data(), tag_size, MPI_INT, 0, PETSC_COMM_WORLD);
  MPI_Bcast(dtype.data(), tag_size, MPI_INT, 0, PETSC_COMM_WORLD);
  MPI_Bcast(storage.data(), tag_size, MPI_INT, 0, PETSC_COMM_WORLD);
 
  // create new tag on other processors
  // loop over tag index to support multiple tags projection in one run
 
  for (size_t index = 0; index < interp_tags.size(); index++) {
    // check if tag exists first
    if (world_rank) {
      Tag old_interp_tag;
      auto rval_check_tag =
          moab.tag_get_handle(listTagsToProject[index].c_str(), old_interp_tag);
      if (rval_check_tag == MB_SUCCESS) {
        MOFEM_LOG("WORLD", Sev::verbose)
            << "Deleting existing tag on target mesh (post-projection): "
            << listTagsToProject[index];
        CHKERR moab.tag_delete(old_interp_tag);
      }
    }
    Tag interp_tag_all;
    unsigned flags =
        MB_TAG_CREAT | storage[index]; // e.g., MB_TAG_DENSE or MB_TAG_SPARSE
    std::vector<double> def_val(tag_length[index], 0.);
    auto rval = moab.tag_get_handle(listTagsToProject[index].c_str(),
                                    tag_length[index], dtype[index],
                                    interp_tag_all, flags, def_val.data());
    if (rval != MB_SUCCESS && world_rank) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unable to create projection tag %s",
              listTagsToProject[index].c_str());
    }
 
    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[index],
        Sev::inform);
    auto volume_exchange = CommInterface::createEntitiesPetscVector(
        mField.get_comm(), mField.get_moab(), 3, tag_length[index],
        Sev::inform);
 
    CHKERR CommInterface::updateEntitiesPetscVector(
        mField.get_moab(), vertex_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,
                                      const bool add_bubble) {
  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);
    if (add_bubble) {
      CHKERR add_field_to_fe(elementVolumeName, bubbleField);
    }
    if (!isNoStretch())
      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 add_field_to_fe(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_row(fe, materialH1Positions);
    CHKERR mField.modify_finite_element_add_field_col(fe, materialH1Positions);
    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 (bcSpatialSpringVecPtr) {
      for (auto &v : *bcSpatialSpringVecPtr) {
        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 (!isNoStretch()) {
    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 DMMoFEMSetSquareProblem(dmMaterial, PETSC_FALSE);
  CHKERR DMMoFEMSetIsPartitioned(dmMaterial, PETSC_TRUE);
  CHKERR DMMoFEMSetDestroyProblem(dmMaterial, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmMaterial, materialH1Positions);
  CHKERR DMMoFEMAddSubFieldCol(dmMaterial, piolaStress);
  CHKERR DMMoFEMAddSubFieldCol(dmMaterial, bubbleField);
  CHKERR DMMoFEMAddSubFieldCol(dmMaterial, spatialL2Disp);
  CHKERR DMMoFEMAddSubFieldCol(dmMaterial, rotAxis);
  if (!isNoStretch()) {
    CHKERR DMMoFEMAddSubFieldCol(dmMaterial, stretchTensor);
  }
  CHKERR DMMoFEMAddSubFieldCol(dmMaterial, hybridSpatialDisp);
  CHKERR DMMoFEMAddSubFieldCol(dmMaterial, contactDisp);
  CHKERR DMMoFEMAddElement(dmMaterial, elementVolumeName);
  CHKERR DMMoFEMAddElement(dmMaterial, naturalBcElement);
  CHKERR DMMoFEMAddElement(dmMaterial, skeletonElement);
  CHKERR DMSetUp(dmMaterial);
 
  auto set_zero_block = [&]() {
    MoFEMFunctionBegin;
    if (!isNoStretch()) {
      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 (!isNoStretch()) {
      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);
    }
 
    auto zero_kinetic_constraints_block = [&]() {
      MoFEMFunctionBegin;
      // we shoudl have sparet bloc names for this. TOPO_FIX_X, TOPO_FIX_Y,
      // TOPO_FIX_Z, TOPO_FIX_ALL
      auto bc_mng = mField.getInterface<BcManager>();
      CHKERR bc_mng->removeBlockDOFsOnEntities("MATERIAL_PROBLEM", "REMOVE_X",
                                               materialH1Positions, 0, 0);
      CHKERR bc_mng->removeBlockDOFsOnEntities("MATERIAL_PROBLEM", "REMOVE_Y",
                                               materialH1Positions, 1, 1);
      CHKERR bc_mng->removeBlockDOFsOnEntities("MATERIAL_PROBLEM", "REMOVE_Z",
                                               materialH1Positions, 2, 2);
      CHKERR bc_mng->removeBlockDOFsOnEntities("MATERIAL_PROBLEM", "REMOVE_ALL",
                                               materialH1Positions, 0, 3);
      CHKERR bc_mng->removeBlockDOFsOnEntities("MATERIAL_PROBLEM", "FIX_X",
                                               materialH1Positions, 0, 0);
      CHKERR bc_mng->removeBlockDOFsOnEntities("MATERIAL_PROBLEM", "FIX_Y",
                                               materialH1Positions, 1, 1);
      CHKERR bc_mng->removeBlockDOFsOnEntities("MATERIAL_PROBLEM", "FIX_Z",
                                               materialH1Positions, 2, 2);
      CHKERR bc_mng->removeBlockDOFsOnEntities("MATERIAL_PROBLEM", "FIX_ALL",
                                               materialH1Positions, 0, 3);
      MoFEMFunctionReturn(0);
    };
 
    // CHKERR zero_kinetic_constraints_block();
 
    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,
               std::string load_history_file)
    : blockName(name), loadHistoryFile(load_history_file), 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,
             std::string load_history_file)
    : blockName(name), loadHistoryFile(load_history_file), 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,
                       std::string load_history_file)
    : blockName(name), loadHistoryFile(load_history_file), 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,
                                           std::string load_history_file)
    : blockName(name), loadHistoryFile(load_history_file), 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();
}
 
SpringBc::SpringBc(std::string name, std::vector<double> attr, Range faces)
    : blockName(name), faces(faces) {
 
  blockName = name;
  if (attr.size() < 2) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Wrong size of spring BC attributes");
  }
 
  normalStiffness = attr[0];
  tangentialStiffness = attr[1];
 
  MOFEM_LOG("EP", Sev::inform) << "Add SpringBc " << name;
  MOFEM_LOG("EP", Sev::inform) << "Add SpringBc kn " << normalStiffness;
  MOFEM_LOG("EP", Sev::inform) << "Add SpringBc kt " << tangentialStiffness;
  MOFEM_LOG("EP", Sev::inform) << "Add SpringBc nb. of faces " << faces.size();
}
 
PressureBc::PressureBc(std::string name, std::vector<double> attr, Range faces,
                       std::string load_history_file)
    : blockName(name), loadHistoryFile(load_history_file), 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, std::string load_history_file)
    : blockName(name), loadHistoryFile(load_history_file), 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 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,
    std::string load_history_file)
    : blockName(name), faces(faces) {
  (void)load_history_file;
  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,
                                           std::string load_history_file)
    : blockName(name), faces(faces) {
  (void)load_history_file;
  flags.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    flags[ii] = attr.size() < 3 ? 1 : 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 (bcSpatialSpringVecPtr)
    for (auto &v : *bcSpatialSpringVecPtr) {
      (*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,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe,
    const bool add_bubble) {
  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 = (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()));
  if (add_bubble) {
    fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
        bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  }
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
      piolaStress, dataAtPts->getDivPAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
 
  if (isNoStretch()) {
    pushOpCalculateStretchFromStress(
        fe->getOpPtrVector(), physicalEquations, dataAtPts,
        externalStrainVecPtr, timeScaleMap);
  } else {
    fe->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValues<3>(
            stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
  }
 
  CHKERR VecSetDM(solTSStep, PETSC_NULLPTR);
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
      piolaStress, dataAtPts->getApproxP0AtPts(), nullptr, solTSStep));
  if (!isNoStretch()) {
    fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
        stretchTensor, dataAtPts->getLogStretchTensor0AtPts(), solTSStep,
        MBTET));
  }
 
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
      rotAxis, dataAtPts->getRotAxisGradAtPts(), 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 (isNoStretch()) {
    } 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));
    }
  }
 
  // calculate other derived quantities
  fe->getOpPtrVector().push_back(
      new OpCalculateRotationAndSpatialGradient(dataAtPts));
 
  // evaluate integration points
  if (isNoStretch()) {
  } else {
    fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
        do_rhs, do_lhs, dataAtPts, physicalEquations));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::pushNoStretchVolumeA00Ops(
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs) {
  MoFEMFunctionBegin;
 
  bool has_nonhomogeneous_mat_block =
      EshelbianCore::hasNonHomogeneousMaterialBlock;
  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_dP(
      piolaStress, piolaStress, dataAtPts, has_nonhomogeneous_mat_block));
  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dP(
      bubbleField, piolaStress, dataAtPts, has_nonhomogeneous_mat_block));
  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dBubble(
      bubbleField, bubbleField, dataAtPts, has_nonhomogeneous_mat_block));
 
  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) {
    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, alphaR, alphaR0, alphaOmega, alphaOmega0,
      alphaViscousR, alphaViscousR0, alphaViscousOmega,
      alphaViscousOmega0));
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::pushStretchVolumeA00Ops(
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs) {
  MoFEMFunctionBegin;
 
  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) {
    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, alphaR, alphaR0, alphaOmega, alphaOmega0,
      alphaViscousR, alphaViscousR0, alphaViscousOmega,
      alphaViscousOmega0));
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::pushStressGramOps(
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs) {
  MoFEMFunctionBegin;
  CHKERR pushPiolaStressGramOps(fe_lhs);
  fe_lhs->getOpPtrVector().push_back(
      new OpStressGram_dBubble_dP(bubbleField, piolaStress, dataAtPts));
  fe_lhs->getOpPtrVector().push_back(new OpStressGram_dBubble_dBubble(
      bubbleField, bubbleField));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::pushPiolaStressGramOps(
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs) {
  MoFEMFunctionBegin;
  fe_lhs->getOpPtrVector().push_back(
      new OpStressGram_dP_dP(piolaStress, piolaStress));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setVolumeElementOps(
    const int tag, const bool add_elastic, const bool add_material,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
  MoFEMFunctionBegin;
 
  /** Contact requires that body is marked */
  auto get_body_range = [this](auto name, int dim) {
    std::map<int, Range> map;
 
    for (auto m_ptr :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % name).str()
 
                 ))
 
    ) {
      Range ents;
      CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
                                                            dim, ents, true),
                     "by dim");
      map[m_ptr->getMeshsetId()] = ents;
    }
 
    return map;
  };
 
  auto local_tau_sacale = boost::make_shared<double>(1.0);
  using BoundaryEle =
      PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
  using BdyEleOp = BoundaryEle::UserDataOperator;
  struct OpSetTauScale : public BdyEleOp {
    OpSetTauScale(boost::shared_ptr<double> local_tau_sacale, double alphaTau,
                  double alphaTau0,
                  boost::shared_ptr<MatrixDouble> flux_mat_ptr)
        : BdyEleOp(NOSPACE, BdyEleOp::OPSPACE),
          localTauSacale(local_tau_sacale), alphaTau(alphaTau),
          alphaTau0(alphaTau0), fluxMatPtr(flux_mat_ptr) {}
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) override {
      MoFEMFunctionBegin;
      auto &coords = BdyEleOp::getCoords();
      auto [centre, barycenter, h] =
          Tools::getTricircumcenter3d(coords.data().data());
 
      FTENSOR_INDEXES(3, i, J);
      auto t_P = getFTensor2FromMat<3, 3>(fluxMatPtr);
      auto t_normal = getFTensor1NormalsAtGaussPts();
      auto t_w = getFTensor0IntegrationWeight();
      auto nb_gauss_pts = getGaussPts().size2(); 
      double norm = 0;
      for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
        FTensor::Tensor1<double, 3> t_t;
        t_t(i) = t_P(i, J) * t_normal(J);
        norm += t_w * sqrt(t_t(i) * t_t(i));
        ++t_w;
        ++t_normal;
        ++t_P;
      }
 
      *localTauSacale = (alphaTau / h) + alphaTau0 * norm;
 
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<double> localTauSacale;
    boost::shared_ptr<MatrixDouble> fluxMatPtr;
    double alphaTau;
    double alphaTau0;
  };
 
  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, 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, alphaR, alphaR0, alphaOmega,
                              alphaOmega0, alphaViscousR, alphaViscousR0,
                              alphaViscousOmega, alphaViscousOmega0));
    if (isNoStretch()) {
      // do nothing - no stretch approximation
    } else {
      if (!internalStressTagName.empty()) {
        switch (meshTransferInterpOrder) {
        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 mesh transfer interpolation order %d, for "
                  "internal stress",
                  meshTransferInterpOrder);
        }
        // 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>(nullptr, true);
      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);
      op_loop_skeleton_side->getSideFEPtr()->exeTestHook = not_interface_face;
      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>(nullptr, true);
      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>();
      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));
      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, boost::make_shared<double>(1.0),
              solTSStep));
      op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpSetTauScale(local_tau_sacale, alphaTau, alphaTau0,
                            flux_mat_ptr));
 
      // Add stabilization Ugamma Ugamma skeleton
      auto hybrid_ptr = boost::make_shared<MatrixDouble>();
      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_tau_stabilsation_disp_bc_rhs = [&](auto &pip, auto side_fe_name) {
      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);
      op_loop_skeleton_side->getSideFEPtr()->exeTestHook = not_interface_face;
      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>(nullptr, true);
      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>();
      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);
      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, boost::make_shared<double>(1.0),
              solTSStep));
      op_loop_skeleton_side->getOpPtrVector().push_back(new OpSetTauScale(
          local_tau_sacale, alphaTauBcDisp, alphaTauBcDisp0, flux_mat_ptr));
 
      // Diag L2 - L2 volumes
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpTauStabilizationDispRhsBc(
              broken_disp_data_ptr, bcSpatialDispVecPtr, timeScaleMap,
              [local_tau_sacale](double, double, double) {
                return (*local_tau_sacale);
              }));
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpTauStabilizationOpAnalyticalDispBc(
              broken_disp_data_ptr, bcSpatialAnalyticalDisplacementVecPtr,
              timeScaleMap, [local_tau_sacale](double, double, double) {
                return (*local_tau_sacale);
              }));
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpTauStabilizationRotationRhsBc(
              broken_disp_data_ptr, bcSpatialRotationVecPtr, timeScaleMap,
              [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),
          interfaceFaces, pip);
    };
 
    CHKERR set_hybridisation_rhs(fe_rhs->getOpPtrVector());
    CHKERR set_contact_rhs(fe_rhs->getOpPtrVector());
    if (alphaTau > 0.0 || alphaTau0 > 0.0) {
      CHKERR set_tau_stabilsation_rhs(fe_rhs->getOpPtrVector(), skeletonElement,
                                      hybridSpatialDisp);
    }
    if (alphaTauBcDisp > 0.0 || alphaTauBcDisp0 > 0.0) {
      CHKERR set_tau_stabilsation_disp_bc_rhs(fe_rhs->getOpPtrVector(),
                                              naturalBcElement);
    }
    if (interfaceCrack == 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>;
 
    std::string body_force_history;
    CHKERR getStringArgumentFromJsonBlocksets("BODY_FORCE", "load_history",
                                              body_force_history);
    if (body_force_history.empty()) {
      body_force_history = "body_force.txt";
    } else {
      MOFEM_LOG("EP", Sev::inform)
          << "Body force load history from JSON: " << body_force_history;
    }
    auto body_time_scale =
        boost::make_shared<DynamicRelaxationTimeScale>(body_force_history);
    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, fe_lhs);
 
  // elastic
  if (add_elastic) {
 
    if (isNoStretch()) {
      CHKERR pushNoStretchVolumeA00Ops(fe_lhs);
    } else {
      CHKERR pushStretchVolumeA00Ops(fe_lhs);
    }
 
    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>(nullptr, true);
      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);
      op_loop_skeleton_side->getSideFEPtr()->exeTestHook = not_interface_face;
      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>(nullptr, true);
      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));
      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, boost::make_shared<double>(1.0),
              solTSStep));
      op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpSetTauScale(local_tau_sacale, alphaTau, alphaTau0,
                            flux_mat_ptr));
 
      // 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);
          }));
      // 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);
              }));
      // 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));
      // 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));
 
      pip.push_back(op_loop_skeleton_side);
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_tau_stabilsation_disp_bc_lhs = [&](auto &pip, auto side_fe_name) {
      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);
      op_loop_skeleton_side->getSideFEPtr()->exeTestHook = not_interface_face;
      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>(nullptr, true);
      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));
      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, boost::make_shared<double>(1.0),
              solTSStep));
      op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
      op_loop_skeleton_side->getOpPtrVector().push_back(new OpSetTauScale(
          local_tau_sacale, alphaTauBcDisp, alphaTauBcDisp0, flux_mat_ptr));
 
      // Diag L2-L2 volumes
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpTauStabilizationDispLhsBc(
              broken_disp_data_ptr, bcSpatialDispVecPtr,
              [local_tau_sacale](double, double, double) {
                return (*local_tau_sacale);
              }));
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpTauStabilizationOpAnalyticalDispLhsBc(
              broken_disp_data_ptr, bcSpatialAnalyticalDisplacementVecPtr,
              [local_tau_sacale](double, double, double) {
                return (*local_tau_sacale);
              }));
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpTauStabilizationRotationLhsBc(
              broken_disp_data_ptr, bcSpatialRotationVecPtr,
              [local_tau_sacale](double, double, double) {
                return (*local_tau_sacale);
              }));
 
      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),
          interfaceFaces, pip);
    };
 
    CHKERR set_hybridisation_lhs(fe_lhs->getOpPtrVector());
    CHKERR set_contact_lhs(fe_lhs->getOpPtrVector());
    if (alphaTau > 0.0 || alphaTau0 > 0.0) {
      CHKERR set_tau_stabilsation_lhs(fe_lhs->getOpPtrVector(), skeletonElement,
                                      hybridSpatialDisp);
    }
    if (alphaTauBcDisp > 0.0 || alphaTauBcDisp0 > 0.0) {
      CHKERR set_tau_stabilsation_disp_bc_lhs(fe_lhs->getOpPtrVector(),
                                              naturalBcElement);
    }
    if (interfaceCrack == 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>(nullptr, true);
      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));
      pip.push_back(op_loop_domain_side);
      return broken_data_ptr;
    };
 
    auto set_rhs = [&]() {
      MoFEMFunctionBegin;
 
      auto broken_data_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));
 
      auto piola_scale_ptr = boost::make_shared<double>(1.0);
      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));
      fe_rhs->getOpPtrVector().push_back(
          new OpSpringRhsBc(hybridSpatialDisp, hybrid_ptr, broken_data_ptr,
                            bcSpatialSpringVecPtr));
 
      auto get_normal_disp_bc_faces = [&]() {
        auto faces =
            get_range_from_block(mField, "NORMAL_DISPLACEMENT", SPACE_DIM - 1);
        return boost::make_shared<Range>(faces);
      };
 
      auto get_spring_bc_faces = [&]() {
        auto faces = get_range_from_block(mField, "SPRING_BC", 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()));
      fe_rhs->getOpPtrVector().push_back(new OpC_dBroken(
          broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0),
          get_spring_bc_faces()));
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_lhs = [&]() {
      MoFEMFunctionBegin;
 
      auto broken_data_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));
      fe_lhs->getOpPtrVector().push_back(
          new OpSpringLhsBc_dU(hybridSpatialDisp, bcSpatialSpringVecPtr));
      fe_lhs->getOpPtrVector().push_back(new OpSpringLhsBc_dP(
          hybridSpatialDisp, broken_data_ptr, bcSpatialSpringVecPtr));
 
      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);
      };
 
      auto get_spring_bc_faces = [&]() {
        auto faces = get_range_from_block(mField, "SPRING_BC", 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()));
      fe_lhs->getOpPtrVector().push_back(new OpC(
          hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0),
          true, true, get_spring_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/EshelbianTestingMonitor.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(ep_ptr->mField); };
#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);
        auto testing_monitor_ptr =
            boost::make_shared<EshelbianTestingMonitor>(*ep_ptr, monitor_ptr);
        ts_ctx->getLoopsMonitor().push_back(
            TsCtx::PairNameFEMethodPtr(ep_ptr->elementVolumeName, monitor_ptr));
 
        PetscBool test_cook_flg = PETSC_FALSE;
        PetscBool test_cook_pts_flg = PETSC_FALSE;
        PetscInt atom_test = 0;
        CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_cook",
                                   &test_cook_flg, PETSC_NULLPTR);
        CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_cook_pts",
                                   &test_cook_pts_flg, PETSC_NULLPTR);
        CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
                                  PETSC_NULLPTR);
        if (atom_test || test_cook_flg || test_cook_pts_flg) {
          ts_ctx->getLoopsMonitor().push_back(TsCtx::PairNameFEMethodPtr(
              ep_ptr->elementVolumeName, testing_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, PETSC_NULLPTR, {}, ts);
      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()}, ts);
 
      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) {
    PetscInt expected_lin_solver_iterations = 11;
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "",
                              "-test_cook_max_linear_iterations",
                              &expected_lin_solver_iterations, PETSC_NULLPTR);
    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);
 
  // Deprecated options
  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options", "none");
 
  CHKERR PetscOptionsScalar(
      "-dynamic_final_time", "dynamic relaxation final time", "",
      finalPhysicalTime, &finalPhysicalTime, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-dynamic_delta_time",
                            "dynamic relaxation final time", "", physicalDt,
                            &physicalDt, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-dynamic_max_it", "dynamic relaxation iterations", "",
                         physicalMaxSteps, &physicalMaxSteps, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-dynamic_h1_update", "update each ts step", "",
                          physicalH1Update, &physicalH1Update, PETSC_NULLPTR);
 
  PetscOptionsEnd();
 
  MOFEM_LOG("EP", Sev::warning)
      << "Following options are deprecated, use -physical prefix options "
         "instead";
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation final time -dynamic_final_time = "
      << finalPhysicalTime;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation delta time -dynamic_delta_time = " << physicalDt;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation max iterations -dynamic_max_it = "
      << physicalMaxSteps;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation H1 update each step -dynamic_h1_update = "
      << (physicalH1Update ? "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);
  CHKERR TSSetSolution(ts, x);
 
  if (physicalH1Update) {
    CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
    CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  } else {
    CHKERR TSSetPreStep(ts, PETSC_NULLPTR);
    CHKERR TSSetPostStep(ts, PETSC_NULLPTR);
  }
  
  CHKERR TSElasticPostStep::preStepFun(ts);
  CHKERR TSElasticPostStep::postStepFun(ts);
 
  currentPhysicalTime = start_time;
  physicalStepNumber = start_step;
  monitor_ptr->ts = PETSC_NULLPTR;
  monitor_ptr->ts_u = PETSC_NULLPTR;
  monitor_ptr->ts_t = currentPhysicalTime;
  monitor_ptr->ts_step = physicalStepNumber;
  CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
  if (physicalDt <= 0.) {
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "physicalDt must be positive, got %g", physicalDt);
  }
  for (; currentPhysicalTime < finalPhysicalTime;) {
    MOFEM_LOG("EP", Sev::inform)
        << "Load step " << physicalStepNumber << " Time " << currentPhysicalTime
        << " delta time " << physicalDt;
 
    CHKERR TSSetStepNumber(ts, 0);
    CHKERR TSSetTime(ts, 0);
    CHKERR TSSetTimeStep(ts, ts_delta_time);
    CHKERR TSSetSolution(ts, x);
    if (!physicalH1Update) {
      CHKERR TSElasticPostStep::preStepFun(ts);
    }
    CHKERR TSSolve(ts, PETSC_NULLPTR);
    if (!physicalH1Update) {
      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 = PETSC_NULLPTR;
    monitor_ptr->ts_u = x;
    monitor_ptr->ts_t = currentPhysicalTime;
    monitor_ptr->ts_step = physicalStepNumber;
    CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
    ++physicalStepNumber;
    if (physicalStepNumber > physicalMaxSteps)
      break;
 
    const double remainingPhysicalTime =
        finalPhysicalTime - currentPhysicalTime;
    if (physicalDt >= remainingPhysicalTime) {
      currentPhysicalTime = finalPhysicalTime;
    } else {
      currentPhysicalTime += physicalDt;
    }
  }
 
  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::postProcessRestartMesh(const int tag, const std::string file,
                                      std::vector<Tag> tags_to_transfer) {
  MoFEMFunctionBegin;
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
  // write file with only crack surfaces and full mesh
  if (crackingOn && !pcomm->rank()) {
    auto meshsets_mng = mField.getInterface<MeshsetsManager>();
 
    std::vector<boost::shared_ptr<TempMeshset>> meshsets_tmp_list;
    auto &list = meshsets_mng->getMeshsetsMultindex();
    std::vector<Tag> tags_list;
 
    auto meshset_ptr = get_temp_meshset_ptr(mField.get_moab());
 
    for (auto &m : list) {
      meshsets_tmp_list.push_back(get_temp_meshset_ptr(mField.get_moab()));
      EntityHandle new_meshset = *meshsets_tmp_list.back();
      auto meshset = m.getMeshset();
      std::vector<Tag> tmp_tags_list;
      CHKERR mField.get_moab().tag_get_tags_on_entity(meshset, tmp_tags_list);
      Range ents;
      CHKERR mField.get_moab().get_entities_by_handle(meshset, ents, true);
      CHKERR mField.get_moab().add_entities(new_meshset, ents);
      for (auto t : tmp_tags_list) {
        void *tag_vals[1];
        int tag_size[1];
        CHKERR mField.get_moab().tag_get_by_ptr(
            t, &meshset, 1, (const void **)tag_vals, tag_size);
        CHKERR mField.get_moab().tag_set_by_ptr(t, &new_meshset, 1, tag_vals,
                                                tag_size);
      }
      std::vector<std::string> remove_tags;
      remove_tags.push_back("AKDTree_coord_norm");
      remove_tags.push_back("__PARALLEL_");
      remove_tags.push_back("_RefBitLevel");
 
      for (auto t : tmp_tags_list) {
        std::string tag_name;
        CHKERR mField.get_moab().tag_get_name(t, tag_name);
        bool add = true;
 
        for (auto &p : remove_tags) {
          if (tag_name.compare(0, p.size(), p) == 0) {
            add = false;
            break;
          }
        }
 
        if (add)
          tags_list.push_back(t);
      }
    }
 
    for (auto &m_ptr : meshsets_tmp_list) {
      EntityHandle m = *m_ptr;
      CHKERR mField.get_moab().add_entities(*meshset_ptr, &m, 1);
    }
 
    // meshsets_tmp_list has all meshsets to write
    std::sort(tags_list.begin(), tags_list.end());
    auto new_end = std::unique(tags_list.begin(), tags_list.end());
    tags_list.resize(std::distance(tags_list.begin(), new_end));
 
    EntityHandle save_meshset = *meshset_ptr;
    CHKERR mField.get_moab().write_file(file.c_str(), "MOAB", "", &save_meshset,
                                        1, &tags_list[0], tags_list.size());
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EshelbianCore::postProcessResults(const int tag, const std::string file,
                                  Vec f_residual, Vec var_vector, Vec gradient,
                                  std::vector<Tag> tags_to_transfer, TS ts) {
  MoFEMFunctionBegin;
 
  SmartPetscObj<Vec> f_r, v_v;
  if (f_residual != PETSC_NULLPTR || var_vector != PETSC_NULLPTR) {
    SmartPetscObj<VecScatter> scatter;
    SmartPetscObj<Vec> xout;
    xout = createDMVector(dM);
    auto xin = f_residual != PETSC_NULLPTR ? f_residual : var_vector;
    CHKERR mField.getInterface<VecManager>()->vecScatterCreate(
        xin, "ELASTIC_PROBLEM", RowColData::ROW, xout, "ESHELBY_PLASTICITY",
        RowColData::ROW, scatter);
    if (f_residual) {
      f_r = vectorDuplicate(xout);
      CHKERR VecScatterBegin(scatter, f_residual, f_r, INSERT_VALUES,
                             SCATTER_FORWARD);
      CHKERR VecScatterEnd(scatter, f_residual, f_r, INSERT_VALUES,
                           SCATTER_FORWARD);
      CHKERR VecGhostUpdateBegin(f_r, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(f_r, INSERT_VALUES, SCATTER_FORWARD);
    }
    if (var_vector) {
      v_v = createDMVector(dM);
      CHKERR VecScatterBegin(scatter, var_vector, v_v, INSERT_VALUES,
                             SCATTER_FORWARD);
      CHKERR VecScatterEnd(scatter, var_vector, v_v, INSERT_VALUES,
                           SCATTER_FORWARD);
      CHKERR VecGhostUpdateBegin(v_v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(v_v, INSERT_VALUES, SCATTER_FORWARD);
    }
  }
 
  SmartPetscObj<Vec> g;
  if (gradient) {
    SmartPetscObj<VecScatter> scatter;
    g = createDMVector(dM);
    CHKERR mField.getInterface<VecManager>()->vecScatterCreate(
        gradient, "MATERIAL_PROBLEM", RowColData::ROW, g, "ESHELBY_PLASTICITY",
        RowColData::ROW, scatter);
    CHKERR VecScatterBegin(scatter, gradient, g, INSERT_VALUES,
                           SCATTER_FORWARD);
    CHKERR VecScatterEnd(scatter, gradient, g, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateBegin(g, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(g, INSERT_VALUES, SCATTER_FORWARD);
  }
 
  // 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) {
    std::string name;
    CHKERR mField.get_moab().tag_get_name(t, name);
    MOFEM_LOG("EP", Sev::verbose)
        << "Adding tag " << name << " to transfer list for post-processing";
    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);
 
    if (ts != PETSC_NULLPTR) {
      post_proc_ptr->data_ctx |= PetscData::CTX_SET_TIME;
      CHKERR TSGetTime(ts, &(post_proc_ptr->ts_t));
      CHKERR TSGetTimeStep(ts, &(post_proc_ptr->ts_dt));
    }
 
    auto domain_ops = [&](auto &fe, int sense) {
      MoFEMFunctionBegin;
 
      auto bubble_cache = boost::make_shared<CGGUserPolynomialBase::CachePhi>(
          0, 0, MatrixDouble());
      fe.getUserPolynomialBase() = boost::shared_ptr<BaseFunction>(
          new CGGUserPolynomialBase(bubble_cache));
      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(new OpCalculateHVecTensorField<3, 3>(
          piolaStress, dataAtPts->getApproxPAtPts(), piola_scale_ptr));
      constexpr bool add_bubble = true;
      if (add_bubble) {
        fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
            bubbleField, dataAtPts->getApproxPAtPts(), piola_scale_ptr,
            SmartPetscObj<Vec>(), MBMAXTYPE));
      }
      fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
      if (isNoStretch()) {
        pushOpCalculateStretchFromStress(
            fe.getOpPtrVector(), physicalEquations, dataAtPts,
            externalStrainVecPtr, timeScaleMap);
      } else {
        fe.getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<3>(
                stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
      }
      CHKERR VecSetDM(solTSStep, PETSC_NULLPTR);
      fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
          piolaStress, dataAtPts->getApproxP0AtPts(), nullptr, solTSStep));
      if (add_bubble) {
        fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
            bubbleField, dataAtPts->getApproxP0AtPts(), nullptr, solTSStep,
            MBMAXTYPE));
      }
      if (!isNoStretch()) {
        fe.getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<3>(
                stretchTensor, dataAtPts->getLogStretchTensor0AtPts(),
                solTSStep, MBTET));
      }
      if (var_vector) {
        fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
            piolaStress, dataAtPts->getVarPiolaPts(),
            boost::make_shared<double>(1), v_v));
        fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
            bubbleField, dataAtPts->getVarPiolaPts(),
            boost::make_shared<double>(1), v_v, MBMAXTYPE));
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            rotAxis, dataAtPts->getVarRotAxisPts(), v_v, MBTET));
        if (isNoStretch()) {
          fe.getOpPtrVector().push_back(
              physicalEquations->returnOpCalculateVarStretchFromStress(
                  dataAtPts, physicalEquations));
        } else {
          fe.getOpPtrVector().push_back(
              new OpCalculateTensor2SymmetricFieldValues<3>(
                  stretchTensor, dataAtPts->getVarLogStreachPts(), v_v, MBTET));
        }
      }
      if (gradient) {
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            materialH1Positions, dataAtPts->getGradientAtPts(), g));
      }
 
      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(
          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 (!isNoStretch()) {
        vec_fields["EiegnLogStreach"] = dataAtPts->getEigenVals();
      }
      if (var_vector) {
        vec_fields["VarOmega"] = dataAtPts->getVarRotAxisPts();
        vec_fields["VarSpatialDisplacementL2"] =
            boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            spatialL2Disp, vec_fields["VarSpatialDisplacementL2"], v_v, MBTET));
      }
      if (f_residual) {
        vec_fields["ResSpatialDisplacementL2"] =
            boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            spatialL2Disp, vec_fields["ResSpatialDisplacementL2"], f_r, MBTET));
        vec_fields["ResOmega"] = boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            rotAxis, vec_fields["ResOmega"], f_r, MBTET));
      }
      if (gradient) {
        vec_fields["Gradient"] = dataAtPts->getGradientAtPts();
      }
 
      OpPPMap::DataMapMat mat_fields;
      mat_fields["PiolaStress"] = dataAtPts->getApproxPAtPts();
      if (var_vector) {
        mat_fields["VarPiolaStress"] = dataAtPts->getVarPiolaPts();
      }
      if (f_residual) {
        mat_fields["ResPiolaStress"] = boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
            piolaStress, mat_fields["ResPiolaStress"],
            boost::make_shared<double>(1), f_r));
        fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
            bubbleField, mat_fields["ResPiolaStress"],
            boost::make_shared<double>(1), f_r, MBMAXTYPE));
      }
      if (!internalStressTagName.empty()) {
        mat_fields[internalStressTagName] = dataAtPts->getInternalStressAtPts();
        switch (meshTransferInterpOrder) {
        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 mesh transfer interpolation order %d, for "
                  "internal stress",
                  meshTransferInterpOrder);
        }
      }
 
      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) {
        if (!isNoStretch()) {
          mat_fields_symm["ResLogSpatialStretch"] =
              boost::make_shared<MatrixDouble>();
          fe.getOpPtrVector().push_back(
              new OpCalculateTensor2SymmetricFieldValues<3>(
                  stretchTensor, mat_fields_symm["ResLogSpatialStretch"], f_r,
                  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(nullptr, true));
    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);
    pip.push_back(op_this);
 
    op_this->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
            {},
 
            {{"ContactDisplacement", dataAtPts->getContactL2AtPts()}},
 
            {},
 
            {}
 
            )
 
    );
 
    if (f_residual) {
 
      auto contact_residual = boost::make_shared<MatrixDouble>();
      op_this->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(
              contactDisp, contact_residual, f_r, MBTET));
      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 get_interface_from_block = [&](auto block_name) {
    auto vol_eles = get_range_from_block(mField, block_name, SPACE_DIM);
    auto skin = filter_true_skin(mField, get_skin(mField, vol_eles));
    Range faces;
    CHKERR mField.get_moab().get_adjacencies(vol_eles, SPACE_DIM - 1, true,
                                             faces, moab::Interface::UNION);
    faces = subtract(faces, skin);
    return faces;
  };
 
  auto crack_faces = unite(get_crack_faces(*crackFaces), *interfaceFaces);
  // VOLUME_INTERFACE faces were already merged into interfaceFaces in
  // projectGeometry(), after applying REMOVE_INTERFACE exclusions.
  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,
                                          TS ts) {
  MoFEMFunctionBegin;
 
  SmartPetscObj<Vec> f_r;
  if (f_residual != PETSC_NULLPTR) {
    SmartPetscObj<VecScatter> scatter;
    f_r = createDMVector(dM);
    CHKERR mField.getInterface<VecManager>()->vecScatterCreate(
        f_residual, "ELASTIC_PROBLEM", RowColData::ROW, f_r,
        "ESHELBY_PLASTICITY", RowColData::ROW, scatter);
    CHKERR VecScatterBegin(scatter, f_residual, f_r, INSERT_VALUES,
                           SCATTER_FORWARD);
    CHKERR VecScatterEnd(scatter, f_residual, f_r, INSERT_VALUES,
                         SCATTER_FORWARD);
  }
 
  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);
  if (ts != PETSC_NULLPTR) {
    post_proc_ptr->data_ctx |= PetscData::CtxSetTime;
    CHKERR TSGetTime(ts, &post_proc_ptr->ts_t);
    CHKERR TSGetTimeStep(ts, &post_proc_ptr->ts_dt);
  }
  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>(nullptr, true);
  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 OpCalculateVectorFieldValues<SPACE_DIM>(
          rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<3>(
          spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
 
  if (isNoStretch()) {
    pushOpCalculateStretchFromStress(
        op_loop_domain_side->getOpPtrVector(), physicalEquations, dataAtPts,
        externalStrainVecPtr, timeScaleMap);
  } 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,
                                            f_r));
    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 bubble_cache =
      boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0, MatrixDouble());
  post_proc_norm_fe->getUserPolynomialBase() =
      boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase(bubble_cache));
  post_proc_norm_fe->getRuleHook = [](int, int, int) { return -1; };
  post_proc_norm_fe->setRuleHook = SetIntegrationAtFrontVolume(
      frontVertices, frontAdjEdges, vol_rule, bubble_cache);
  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<BcMeshsetType<DISPLACEMENTSET>>(
      "", piolaStress, false, false);
  CHKERR bc_mng->pushMarkDOFsOnEntities<BcDisplacementMeshsetType<BLOCKSET>>(
      "", piolaStress, false, false);
 
  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
  auto get_fix_load_history = [&](const std::string &block_name) {
    for (const auto type_name : {"FIX_X", "FIX_Y", "FIX_Z", "FIX_ALL"}) {
      for (auto it : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
               std::regex(
 
                   (boost::format("%s(.*)") % type_name).str()
 
                       ))
 
      ) {
        if (it->getName() == block_name) {
          return getStringArgumentFromJsonBlockset(
              type_name, it->getMeshsetId(), "load_history");
        }
      }
    }
    return std::string();
  };
 
  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,
                                        get_fix_load_history(block_name));
    }
  }
  // old way of naming blocksets for displacement BCs
  CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
 
  bcSpatialNormalDisplacementVecPtr =
      boost::make_shared<NormalDisplacementBcVec>();
  CHKERR getBc(bcSpatialNormalDisplacementVecPtr, "NORMAL_DISPLACEMENT", 1);
 
  bcSpatialSpringVecPtr = boost::make_shared<SpringBcVec>();
  auto mesh_mng = mField.getInterface<MeshsetsManager>();
  for (auto it : mesh_mng->getCubitMeshsetPtr(
           std::regex((boost::format("(.*)%s(.*)") % "SPRING_BC").str()))) {
    std::vector<double> block_attributes;
    CHKERR it->getAttributes(block_attributes);
    if (block_attributes.size() < 2) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "In block %s expected 2 attributes, but given %ld",
              it->getName().c_str(), block_attributes.size());
    }
    Range faces;
    CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                               true);
    MOFEM_LOG("EP", Sev::inform)
        << "Found spring BC on block " << it->getName();
    MOFEM_LOG("EP", Sev::inform)
        << "  kn = " << block_attributes[0] << ", kt = " << block_attributes[1];
    MOFEM_LOG("EP", Sev::inform) << "  nb. of faces " << faces.size();
    bcSpatialSpringVecPtr->emplace_back(it->getName(), block_attributes, faces);
  }
 
  bcSpatialAnalyticalDisplacementVecPtr =
      boost::make_shared<AnalyticalDisplacementBcVec>();
  CHKERR getBc(bcSpatialAnalyticalDisplacementVecPtr, "ANALYTICAL_DISPLACEMENT",
               3);
 
  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;
    if (!bc.loadHistoryFile.empty()) {
      MOFEM_LOG("EP", Sev::inform)
          << "Displacement load history from JSON for " << bc.blockName << ": "
          << bc.loadHistoryFile;
      timeScaleMap[bc.blockName] =
          boost::make_shared<DynamicRelaxationTimeScale>(bc.loadHistoryFile);
    } else {
      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;
    if (!bc.loadHistoryFile.empty()) {
      MOFEM_LOG("EP", Sev::inform)
          << "Normal displacement load history from JSON for " << bc.blockName
          << ": " << bc.loadHistoryFile;
      timeScaleMap[bc.blockName] =
          boost::make_shared<DynamicRelaxationTimeScale>(bc.loadHistoryFile);
    } else {
      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>();
  CHKERR getBc(bcSpatialPressureVecPtr, "PRESSURE", 1);
 
  bcSpatialAnalyticalTractionVecPtr =
      boost::make_shared<AnalyticalTractionBcVec>();
  CHKERR getBc(bcSpatialAnalyticalTractionVecPtr, "ANALYTICAL_TRACTION", 3);
 
  auto ts_traction =
      boost::make_shared<DynamicRelaxationTimeScale>("traction_history.txt");
  for (auto &bc : *bcSpatialTractionVecPtr) {
    if (!bc.loadHistoryFile.empty()) {
      MOFEM_LOG("EP", Sev::inform)
          << "Traction load history from JSON for " << bc.blockName << ": "
          << bc.loadHistoryFile;
      timeScaleMap[bc.blockName] =
          boost::make_shared<DynamicRelaxationTimeScale>(bc.loadHistoryFile);
    } else {
      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) {
    if (!bc.loadHistoryFile.empty()) {
      MOFEM_LOG("EP", Sev::inform)
          << "Pressure load history from JSON for " << bc.blockName << ": "
          << bc.loadHistoryFile;
      timeScaleMap[bc.blockName] =
          boost::make_shared<DynamicRelaxationTimeScale>(bc.loadHistoryFile);
    } else {
      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;
      const bool analytical_external_strain = std::regex_match(
          it->getName(), std::regex("(.*)ANALYTICAL_EXTERNALSTRAIN(.*)"));
      const std::string json_block_name =
          analytical_external_strain ? "ANALYTICAL_EXTERNALSTRAIN" : block_name;
 
      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;
      };
 
      std::string load_history;
      if (!analytical_external_strain) {
        load_history = getStringArgumentFromJsonBlockset(
            json_block_name, it->getMeshsetId(), "load_history");
      }
      ext_strain_vec_ptr->emplace_back(it->getName(), block_attributes,
                                       get_block_ents(), load_history);
    }
    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;
    if (!ext_strain_block.loadHistoryFile.empty()) {
      MOFEM_LOG("EP", Sev::inform)
          << "External strain load history from JSON for "
          << ext_strain_block.blockName << ": "
          << ext_strain_block.loadHistoryFile;
      timeScaleMap[ext_strain_block.blockName] =
          boost::make_shared<DynamicRelaxationTimeScale>(
              ext_strain_block.loadHistoryFile);
    } else {
      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),
      SmartPetscObj<TS>(ts, true));
 
  // Deprecated options
  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options", "none");
 
  CHKERR PetscOptionsScalar(
      "-dynamic_final_time", "dynamic relaxation final time", "",
      finalPhysicalTime, &finalPhysicalTime, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-dynamic_delta_time",
                            "dynamic relaxation final time", "", physicalDt,
                            &physicalDt, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-dynamic_max_it", "dynamic relaxation iterations", "",
                         physicalMaxSteps, &physicalMaxSteps, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-dynamic_h1_update", "update each ts step", "",
                          physicalH1Update, &physicalH1Update, PETSC_NULLPTR);
 
  PetscOptionsEnd();
 
  EshelbianCore::physicalTimeFlg = PETSC_TRUE;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation final time -dynamic_final_time = "
      << finalPhysicalTime;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation delta time -dynamic_delta_time = " << physicalDt;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation max iterations -dynamic_max_it = "
      << physicalMaxSteps;
  MOFEM_LOG("EP", Sev::inform)
      << "Dynamic relaxation H1 update each step -dynamic_h1_update = "
      << (physicalH1Update ? "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 (physicalH1Update) {
    CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
    CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  }
 
  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());
 
  currentPhysicalTime = start_time;
  physicalStepNumber = start_step;
  monitor_ptr->ts = PETSC_NULLPTR;
  monitor_ptr->ts_u = PETSC_NULLPTR;
  monitor_ptr->ts_t = currentPhysicalTime;
  monitor_ptr->ts_step = physicalStepNumber;
  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);
 
  if (physicalDt <= 0.) {
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "physicalDt must be positive, got %g", physicalDt);
  }
  for (; currentPhysicalTime < finalPhysicalTime;) {
    MOFEM_LOG("EP", Sev::inform)
        << "Load step " << physicalStepNumber << " Time " << currentPhysicalTime
        << " delta time " << physicalDt;
 
    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);
 
    if (!physicalH1Update && physicalStepNumber > start_step) {
      CHKERR TSElasticPostStep::preStepFun(ts);
    }
    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 = PETSC_NULLPTR;
    monitor_ptr->ts_u = x;
    monitor_ptr->ts_t = currentPhysicalTime;
    monitor_ptr->ts_step = physicalStepNumber;
    CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
    if (!physicalH1Update) {
      CHKERR TSElasticPostStep::postStepFun(ts);
    }
 
    ++physicalStepNumber;
    if (physicalStepNumber > physicalMaxSteps)
      break;
 
    const double remainingPhysicalTime =
        finalPhysicalTime - currentPhysicalTime;
    if (physicalDt >= remainingPhysicalTime) {
      currentPhysicalTime = finalPhysicalTime;
    } else {
      currentPhysicalTime += physicalDt;
    }
  }
 
  CHKERR TSElasticPostStep::postStepDestroy();
  TetPolynomialBase::switchCacheBaseOff<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::solveLoadFactor(TS ts, Vec x, int start_step,
                                              double start_time) {
  MoFEMFunctionBegin;
 
  loadFactorTSSolveExecuted = PETSC_FALSE;
 
  auto storage = solve_elastic_setup::setup(this, ts, x, 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();
 
  auto test_monitor_ptr =
      boost::make_shared<EshelbianTestingMonitor>(*this, monitor_ptr);
 
  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 (physicalH1Update) {
    CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
    CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  }
 
  CHKERR TSElasticPostStep::preStepFun(ts);
  CHKERR TSElasticPostStep::postStepFun(ts);
 
  double load_factor_change_clip = 0.1;
 
  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Load Factor Options", "none");
 
  CHKERR PetscOptionsScalar("-initial_load_factor", "Initial load factor", "",
                            loadFactor, &loadFactor, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar(
      "-max_crack_ext_area", "Maximum crack extension area", "",
      maxCrackExtension, &maxCrackExtension, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar(
      "-clip_load_factor_percent", "Upper bound for load factor change", "",
      load_factor_change_clip, &load_factor_change_clip, PETSC_NULLPTR);
  PetscOptionsEnd();
 
  currentPhysicalTime = start_time;
  physicalStepNumber = start_step;
  monitor_ptr->ts = ts;
  monitor_ptr->ts_u = PETSC_NULLPTR;
  monitor_ptr->ts_t = currentPhysicalTime;
  monitor_ptr->ts_step = physicalStepNumber;
  CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
  PetscBool test_cook_flg = PETSC_FALSE;
  PetscInt atom_test = 0;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_cook", &test_cook_flg,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
                            PETSC_NULLPTR);
  if (atom_test || test_cook_flg) {
    test_monitor_ptr->ts = ts;
    test_monitor_ptr->ts_u = PETSC_NULLPTR;
    test_monitor_ptr->ts_t = currentPhysicalTime;
    test_monitor_ptr->ts_step = physicalStepNumber;
 
    CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName,
                                    test_monitor_ptr);
  }
 
  MOFEM_LOG("EP", Sev::inform)
      << "Initial crack area: " << *currentCrackAreaPtr;
  MOFEM_LOG("EP", Sev::inform) << "Initial load factor: " << loadFactor;
  MOFEM_LOG("EP", Sev::inform)
      << "Initial crack front energy: " << avgGriffithsEnergy;
 
  for (; currentPhysicalTime < finalPhysicalTime;) {
    MOFEM_LOG("EP", Sev::inform)
        << "Load step " << physicalStepNumber << " Load Factor "
        << currentPhysicalTime << " delta load factor " << physicalDt;
 
    oldCrackArea = *currentCrackAreaPtr;
    oldStrainEnergy = strainEnergy;
    oldLoadFactor = loadFactor;
 
    CHKERR TSSetStepNumber(ts, 0);
    CHKERR TSSetTime(ts, 0);
    CHKERR TSSetTimeStep(ts, ts_delta_time);
    if (!physicalH1Update) {
      CHKERR TSElasticPostStep::preStepFun(ts);
    }
    CHKERR TSSetSolution(ts, x);
    CHKERR TSSolve(ts, PETSC_NULLPTR);
    loadFactorTSSolveExecuted = PETSC_TRUE;
    if (!physicalH1Update) {
      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 = ts;
    monitor_ptr->ts_u = x;
    monitor_ptr->ts_t = currentPhysicalTime;
    monitor_ptr->ts_step = physicalStepNumber;
    CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
    if (atom_test || test_cook_flg) {
      test_monitor_ptr->ts = ts;
      test_monitor_ptr->ts_u = x;
      test_monitor_ptr->ts_t = currentPhysicalTime;
      test_monitor_ptr->ts_step = physicalStepNumber;
      CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName,
                                      test_monitor_ptr);
    }
 
    if (mField.get_comm_rank() == 0) {
      const double delta_area = *currentCrackAreaPtr - oldCrackArea;
      const bool has_crack_extension = delta_area > 0.0;
 
      if (has_crack_extension) {
        const double denom = 0.5 * std::abs(avgGriffithsEnergy);
        if (denom > 0.0) {
          const double updated_load_factor =
              oldLoadFactor * std::sqrt(griffithEnergy / denom);
          loadFactor = std::max(updated_load_factor, 1.0e-6);
        }
      }
 
      const bool is_first_step = physicalStepNumber == start_step;
      const double initial_step_range = 5;
      const double min_load_factor = 1.0e-6;
      const double max_load_factor =
          oldLoadFactor * (1.0 + load_factor_change_clip);
 
      if (physicalStepNumber >= start_step + initial_step_range) {
        loadFactor = std::clamp(loadFactor, min_load_factor, max_load_factor);
        MOFEM_LOG("EP", Sev::inform)
            << "Allowable range for load factor [" << min_load_factor << ", "
            << max_load_factor << "]";
      }
 
      const double previous_load_factor = is_first_step ? 0. : oldLoadFactor;
      physicalDt = loadFactor - previous_load_factor;
 
      MOFEM_LOG("EP", Sev::inform)
          << "Setting new load factor to: " << loadFactor;
    }
    double load_control_data[] = {physicalDt, loadFactor};
    CHKERR MPI_Bcast(load_control_data, 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    physicalDt = load_control_data[0];
    loadFactor = load_control_data[1];
 
    ++physicalStepNumber;
    if (physicalStepNumber > physicalMaxSteps)
      break;
 
    const double remainingPhysicalTime =
        finalPhysicalTime - currentPhysicalTime;
    if (physicalDt >= remainingPhysicalTime) {
      currentPhysicalTime = finalPhysicalTime;
    } else {
      currentPhysicalTime += physicalDt;
    }
  }
 
  CHKERR TSElasticPostStep::postStepDestroy();
  TetPolynomialBase::switchCacheBaseOff<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
  MOFEM_LOG("EP", Sev::inform) << "Final load factor: " << loadFactor;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::solveSchapeOptimisation(TS ts, Vec x,
                                                      int start_step,
                                                      double start_time) {
  MoFEMFunctionBegin;
 
  auto storage = solve_elastic_setup::setup(this, ts, x, false);
 
  auto topological_tao_ctx = createTopologicalTAOCtx(
      this, SetIntegrationAtFrontVolume(frontVertices, frontAdjEdges),
      SetIntegrationAtFrontFace(frontVertices, frontAdjEdges),
      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();
 
  EshelbianCore::physicalTimeFlg = PETSC_TRUE;
  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);
 
  auto tao = createTao(mField.get_comm());
  CHKERR TaoSetType(tao, TAOLMVM);
  auto g = createDMVector(dmMaterial, RowColData::ROW);
  CHKERR TaoSetObjectiveAndGradient(tao, g,
                                    topologicalEvaluateObjectiveAndGradient,
                                    (void *)topological_tao_ctx.get());
 
  currentPhysicalTime = start_time;
  physicalStepNumber = start_step;
  monitor_ptr->ts = PETSC_NULLPTR;
  monitor_ptr->ts_u = PETSC_NULLPTR;
  monitor_ptr->ts_t = currentPhysicalTime;
  monitor_ptr->ts_step = physicalStepNumber;
  CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
  auto tao_sol0 = createDMVector(dmMaterial, RowColData::ROW);
  CHKERR DMoFEMMeshToLocalVector(dmMaterial, tao_sol0, INSERT_VALUES,
                                 SCATTER_FORWARD, RowColData::ROW);
  CHKERR VecGhostUpdateBegin(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
 
  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)
      << "Toplogical data vector size " << tao_sol_size << " local size "
      << tao_sol_loc_size << " number of interface faces "
      << interfaceFaces->size();
 
  CHKERR TaoSetFromOptions(tao);
 
  if (delta_time <= 0.) {
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "delta_time must be positive, got %g", delta_time);
  }
  for (; currentPhysicalTime < final_time;) {
    MOFEM_LOG("EP", Sev::inform)
        << "Load step " << physicalStepNumber << " Time " << currentPhysicalTime
        << " 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 = topological_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 = PETSC_NULLPTR;
    monitor_ptr->ts_u = x;
    monitor_ptr->ts_t = currentPhysicalTime;
    monitor_ptr->ts_step = physicalStepNumber;
    CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
    ++physicalStepNumber;
    if (physicalStepNumber > max_it)
      break;
 
    const double remainingPhysicalTime = final_time - currentPhysicalTime;
    if (delta_time >= remainingPhysicalTime) {
      currentPhysicalTime = final_time;
    } else {
      currentPhysicalTime += delta_time;
    }
  }
 
  CHKERR TSElasticPostStep::postStepDestroy();
  TetPolynomialBase::switchCacheBaseOff<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::solveTestTopologicalDerivative(TS ts, Vec x,
                                                             int start_step,
                                                             double start_time) {
  MoFEMFunctionBegin;
 
  auto storage = solve_elastic_setup::setup(this, ts, x, false);
 
  auto topological_tao_ctx = createTopologicalTAOCtx(
      this, SetIntegrationAtFrontVolume(frontVertices, frontAdjEdges),
      SetIntegrationAtFrontFace(frontVertices, frontAdjEdges),
      SmartPetscObj<TS>(ts, true));
 
  EshelbianCore::physicalTimeFlg = PETSC_TRUE;
  CHKERR addDebugModel(ts);
 
  auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
 
  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 (physicalH1Update) {
    CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
    CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  }
 
  CHKERR TSElasticPostStep::preStepFun(ts);
  CHKERR TSElasticPostStep::postStepFun(ts);
 
  const bool restart_run =
      start_step != 0 ||
      std::abs(start_time) > std::numeric_limits<double>::epsilon();
  const double test_time = restart_run ? start_time : finalPhysicalTime;
  if (!restart_run &&
      std::abs(test_time) < std::numeric_limits<double>::epsilon()) {
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Set non-zero -physical_final_time for test_topological_derivative");
  }
 
  currentPhysicalTime = test_time;
  physicalStepNumber = start_step;
  monitor_ptr->ts = PETSC_NULLPTR;
  monitor_ptr->ts_u = PETSC_NULLPTR;
  monitor_ptr->ts_t = currentPhysicalTime;
  monitor_ptr->ts_step = physicalStepNumber;
  CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
  MOFEM_LOG("EP", Sev::inform)
      << "Solving load step before topological derivative test: "
      << physicalStepNumber << " Time " << currentPhysicalTime
      << " TS delta time " << ts_delta_time;
 
  CHKERR TSSetStepNumber(ts, 0);
  CHKERR TSSetTime(ts, 0);
  CHKERR TSSetTimeStep(ts, ts_delta_time);
  if (!physicalH1Update) {
    CHKERR TSElasticPostStep::preStepFun(ts);
  }
  CHKERR TSSetSolution(ts, x);
  CHKERR TSSolve(ts, PETSC_NULLPTR);
  if (!physicalH1Update) {
    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 = PETSC_NULLPTR;
  monitor_ptr->ts_u = x;
  monitor_ptr->ts_t = currentPhysicalTime;
  monitor_ptr->ts_step = physicalStepNumber;
  CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
  auto tao_sol0 = createDMVector(dmMaterial, RowColData::ROW);
  CHKERR DMoFEMMeshToLocalVector(dmMaterial, tao_sol0, INSERT_VALUES,
                                 SCATTER_FORWARD, RowColData::ROW);
  CHKERR VecGhostUpdateBegin(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(tao_sol0, INSERT_VALUES, SCATTER_FORWARD);
 
  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)
      << "Topological data vector size " << tao_sol_size << " local size "
      << tao_sol_loc_size << " number of interface faces "
      << interfaceFaces->size();
 
  const char *list_objective_models[ObjectiveModelType::LAST_MODEL] = {
      "python_model", "hencky_model"};
#ifdef ENABLE_PYTHON_BINDING
  PetscInt choice_objective_model = ObjectiveModelType::PYTHON_MODEL;
#else
  PetscInt choice_objective_model = ObjectiveModelType::HENCKY_MODEL;
#endif
  CHKERR PetscOptionsGetEList(
      PETSC_NULLPTR, PETSC_NULLPTR, "-objective_model_type",
      list_objective_models, ObjectiveModelType::LAST_MODEL,
      &choice_objective_model, PETSC_NULLPTR);
  const auto objective_model_type =
      static_cast<ObjectiveModelType>(choice_objective_model);
  MOFEM_LOG("EP", Sev::inform)
      << "Objective model type: -objective_model_type "
      << list_objective_models[objective_model_type];
 
  auto g = createDMVector(dmMaterial, RowColData::ROW);
  PetscReal obj_value;
  CHKERR testTopologicalDerivative(topological_tao_ctx.get(), tao_sol0,
                                   &obj_value, g, objective_model_type);
 
  CHKERR TSElasticPostStep::postStepDestroy();
  TetPolynomialBase::switchCacheBaseOff<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
 
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity
 
#include <impl/EshelbianFracture.cpp>