EshelbianPlasticity.cpp

Eshelbian plasticity implementation

/**
 * \file EshelbianPlasticity.cpp
 * \example EshelbianPlasticity.cpp
 *
 * \brief Eshelbian plasticity implementation
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <CGGTonsorialBubbleBase.hpp>
 
#include <EshelbianPlasticity.hpp>
#include <boost/math/constants/constants.hpp>
 
#include <cholesky.hpp>
#ifdef PYTHON_SDF
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
 
#pragma message "With PYTHON_SDF"
 
#else
 
#pragma message "Without PYTHON_SDF"
 
#endif
 
#include <EshelbianContact.hpp>
 
namespace EshelbianPlasticity {
struct VolUserDataOperatorStabAssembly
    : public VolumeElementForcesAndSourcesCore::UserDataOperator {
  using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
struct FaceUserDataOperatorStabAssembly
    : public FaceElementForcesAndSourcesCore::UserDataOperator {
  using FaceElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
};
 
} // namespace EshelbianPlasticity
 
static auto get_range_from_block(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) {
    Range faces;
    CHK_MOAB_THROW(
        bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 2, faces, true),
        "get meshset ents");
    r.merge(faces);
  }
 
  return r;
};
 
template <>
typename MoFEM::OpBaseImpl<
    SCHUR,
    EshelbianPlasticity::VolUserDataOperatorStabAssembly>::MatSetValuesHook
    MoFEM::OpBaseImpl<SCHUR,
                      EshelbianPlasticity::VolUserDataOperatorStabAssembly>::
        matSetValuesHook = [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
                              const EntitiesFieldData::EntData &row_data,
                              const EntitiesFieldData::EntData &col_data,
                              MatrixDouble &m) {
          return MatSetValues<AssemblyTypeSelector<SCHUR>>(
              op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
        };
 
namespace EshelbianPlasticity {
 
enum RotSelector EshelbianCore::rotSelector = LARGE_ROT;
enum RotSelector EshelbianCore::gradApperoximator = LARGE_ROT;
enum StretchSelector EshelbianCore::stretchSelector = LOG;
PetscBool EshelbianCore::noStretch = PETSC_FALSE;
 
double EshelbianCore::exponentBase = exp(1);
boost::function<double(const double)> EshelbianCore::f = EshelbianCore::f_log;
boost::function<double(const double)> EshelbianCore::d_f =
    EshelbianCore::d_f_log;
boost::function<double(const double)> EshelbianCore::dd_f =
    EshelbianCore::dd_f_log;
boost::function<double(const double)> EshelbianCore::inv_f =
    EshelbianCore::inv_f_log;
boost::function<double(const double)> EshelbianCore::inv_d_f =
    EshelbianCore::inv_d_f_log;
boost::function<double(const double)> EshelbianCore::inv_dd_f =
    EshelbianCore::inv_dd_f_log;
 
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), piolaStress("P"), eshelbyStress("S"),
      spatialL2Disp("wL2"), spatialH1Disp("wH1"), contactDisp("contactDisp"),
      materialL2Disp("W"), stretchTensor("u"), rotAxis("omega"),
      materialGradient("G"), tauField("TAU"), lambdaField("LAMBDA"),
      bubbleField("BUBBLE"), elementVolumeName("EP"),
      naturalBcElement("NATURAL_BC"), essentialBcElement("ESSENTIAL_BC"),
      skinElement("SKIN_ELEMENT"), contactElement("CONTACT_ELEMENT") {
 
  ierr = getOptions();
  CHKERRABORT(PETSC_COMM_WORLD, ierr);
}
 
EshelbianCore::~EshelbianCore() = default;
 
MoFEMErrorCode EshelbianCore::getOptions() {
  MoFEMFunctionBegin;
  const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
  PetscInt choice_rot = EshelbianCore::rotSelector;
  PetscInt choice_grad = EshelbianCore::gradApperoximator;
 
  const char *list_stretches[] = {"linear", "log"};
  PetscInt choice_stretch = StretchSelector::LOG;
 
  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
                           "none");
  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
                         spaceOrder, &spaceOrder, PETSC_NULL);
  CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
                         spaceH1Order, &spaceH1Order, PETSC_NULL);
  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
                         "", materialOrder, &materialOrder, PETSC_NULL);
  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
                            &alphaU, PETSC_NULL);
  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
                            &alphaW, PETSC_NULL);
  CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
                            &alphaRho, PETSC_NULL);
  CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
                            precEps, &precEps, PETSC_NULL);
  CHKERR PetscOptionsScalar("-preconditioner_eps_omega", "preconditioner_eps",
                            "", precEpsOmega, &precEpsOmega, PETSC_NULL);
  CHKERR PetscOptionsScalar("-preconditioner_eps_w", "preconditioner_eps", "",
                            precEpsW, &precEpsW, PETSC_NULL);
  CHKERR PetscOptionsScalar("-preconditioner_eps_contact_disp",
                            "preconditioner_eps", "", precEpsContactDisp,
                            &precEpsContactDisp, PETSC_NULL);
  CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
                           LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
                           PETSC_NULL);
  CHKERR PetscOptionsEList("-grad", "gradient of defamation approximate", "",
                           list_rots, NO_H1_CONFIGURATION + 1,
                           list_rots[choice_grad], &choice_grad, PETSC_NULL);
  CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
                            &EshelbianCore::exponentBase, PETSC_NULL);
  CHKERR PetscOptionsEList(
      "-stretches", "stretches", "", list_stretches, StretchSelector::LOG + 1,
      list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
 
  CHKERR PetscOptionsBool("-no_stretch", "do not solve for stretch", "",
                          noStretch, &noStretch, PETSC_NULL);
 
  ierr = PetscOptionsEnd();
  CHKERRG(ierr);
 
  EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
  EshelbianCore::gradApperoximator = static_cast<RotSelector>(choice_grad);
  EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
 
  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;
    EshelbianCore::d_f = d_f_log;
    EshelbianCore::dd_f = dd_f_log;
    EshelbianCore::inv_f = inv_f_log;
    EshelbianCore::inv_d_f = inv_d_f_log;
    EshelbianCore::inv_dd_f = inv_dd_f_log;
    break;
  default:
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
    break;
  };
 
  precEpsOmega += precEps;
  precEpsW += precEps;
 
  MOFEM_LOG("EP", Sev::inform) << "spaceOrder " << spaceOrder;
  MOFEM_LOG("EP", Sev::inform) << "spaceH1Order " << spaceH1Order;
  MOFEM_LOG("EP", Sev::inform) << "materialOrder " << materialOrder;
  MOFEM_LOG("EP", Sev::inform) << "alphaU " << alphaU;
  MOFEM_LOG("EP", Sev::inform) << "alphaW " << alphaW;
  MOFEM_LOG("EP", Sev::inform) << "alphaRho " << alphaRho;
  MOFEM_LOG("EP", Sev::inform) << "precEps " << precEps;
  MOFEM_LOG("EP", Sev::inform) << "precEpsOmega " << precEpsOmega;
  MOFEM_LOG("EP", Sev::inform) << "precEpsW " << precEpsW;
  MOFEM_LOG("EP", Sev::inform)
      << "Rotations " << list_rots[EshelbianCore::rotSelector];
  MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
                               << list_rots[EshelbianCore::gradApperoximator];
  if (exponentBase != exp(1))
    MOFEM_LOG("EP", Sev::inform)
        << "Base exponent " << EshelbianCore::exponentBase;
  else
    MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
  MOFEM_LOG("EP", Sev::inform) << "Stretch " << list_stretches[choice_stretch];
 
  if (spaceH1Order == -1)
    spaceH1Order = spaceOrder;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::addFields(const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  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);
 
  if (bcSpatialDispVecPtr)
    for (auto &v : *bcSpatialDispVecPtr) {
      tets_skin = subtract(tets_skin, v.faces);
    }
  if (bcSpatialRotationVecPtr)
    for (auto &v : *bcSpatialRotationVecPtr) {
      tets_skin = subtract(tets_skin, v.faces);
    }
  if (bcSpatialTraction)
    for (auto &v : *bcSpatialTraction) {
      tets_skin = subtract(tets_skin, v.faces);
    }
 
  auto subtract_faces_where_displacements_are_applied =
      [&](const std::string block_name) {
        MoFEMFunctionBegin;
        auto contact_range = get_range_from_block(mField, block_name);
        tets_skin = subtract(tets_skin, contact_range);
        MoFEMFunctionReturn(0);
      };
  CHKERR subtract_faces_where_displacements_are_applied("CONTACT");
 
  Range faces;
  CHKERR mField.get_moab().get_adjacencies(tets, SPACE_DIM - 1, true, faces,
                                           moab::Interface::UNION);
  Range faces_not_on_the_skin = subtract(faces, tets_skin);
 
  auto add_hdiv_field = [&](const std::string field_name, const int order,
                            const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, HDIV, DEMKOWICZ_JACOBI_BASE, dim,
                            MB_TAG_SPARSE, MF_ZERO);
    auto field_ptr = mField.get_field_structure(field_name);
    auto meshset = field_ptr->getMeshset();
    CHKERR mField.get_moab().add_entities(meshset, tets);
    CHKERR mField.get_moab().add_entities(meshset, faces_not_on_the_skin);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    CHKERR mField.set_field_order(faces_not_on_the_skin, 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_SPARSE, 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);
  };
 
  // spatial fields
  CHKERR add_hdiv_field(piolaStress, spaceOrder, 3);
  CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
  CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
  CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
  CHKERR add_l2_field(stretchTensor, spaceOrder, 6);
  CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
                               get_range_from_block(mField, "CONTACT"));
 
  // material fields
  // CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
  // CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
  // CHKERR add_l2_field(materialL2Disp, materialOrder - 1, 3);
  // CHKERR add_l2_field(tauField, materialOrder - 1, 1);
  // CHKERR add_l2_field(lambdaField, materialOrder - 1, 1);
 
  // CHKERR add_l2_field(materialGradient + "0", materialOrder - 1, 9);
 
  // spatial displacement
  CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
 
  CHKERR mField.build_fields();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EshelbianCore::addVolumeFiniteElement(const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  // set finite element fields
  auto add_field_to_fe = [this](const std::string fe,
                                const std::string field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
    MoFEMFunctionReturn(0);
  };
 
  if (!mField.check_finite_element(elementVolumeName)) {
    CHKERR mField.add_finite_element(elementVolumeName, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(meshset, MBTET,
                                                     elementVolumeName);
 
    CHKERR add_field_to_fe(elementVolumeName, piolaStress);
    CHKERR add_field_to_fe(elementVolumeName, bubbleField);
    // CHKERR add_field_to_fe(elementVolumeName, eshelbyStress);
    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, materialGradient);
    // CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
    //                                                    materialGradient +
    //                                                    "0");
  }
 
  // build finite elements data structures
  CHKERR mField.build_finite_elements(elementVolumeName);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
EshelbianCore::addBoundaryFiniteElement(const EntityHandle meshset) {
  MoFEMFunctionBegin;
 
  // set finite element fields
  auto add_field_to_fe = [this](const std::string fe,
                                const std::string field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
    MoFEMFunctionReturn(0);
  };
 
  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);
    }
  }
  Range essentail_bc_elements;
  if (bcSpatialTraction) {
    for (auto &v : *bcSpatialTraction) {
      essentail_bc_elements.merge(v.faces);
    }
  }
 
  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, eshelbyStress);
  CHKERR mField.build_finite_elements(naturalBcElement);
 
  CHKERR mField.add_finite_element(essentialBcElement, MF_ZERO);
  CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
                                                   essentialBcElement);
  CHKERR add_field_to_fe(essentialBcElement, piolaStress);
  // CHKERR add_field_to_fe(essentialBcElement, eshelbyStress);
  CHKERR mField.build_finite_elements(essentialBcElement);
 
  auto get_skin = [&]() {
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, 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;
  };
 
  auto skin = filter_true_skin(get_skin());
 
  CHKERR mField.add_finite_element(skinElement, MF_ZERO);
  CHKERR mField.add_ents_to_finite_element_by_type(skin, MBTRI, skinElement);
  CHKERR add_field_to_fe(skinElement, piolaStress);
  CHKERR add_field_to_fe(skinElement, spatialH1Disp);
  CHKERR add_field_to_fe(skinElement, contactDisp);
  // CHKERR add_field_to_fe(skinElement, eshelbyStress);
  CHKERR mField.build_finite_elements(skinElement);
 
  auto contact_range = get_range_from_block(mField, "CONTACT");
  MOFEM_LOG("EP", Sev::inform) << "Contact elements " << contact_range.size();
 
  CHKERR mField.add_finite_element(contactElement, MF_ZERO);
  CHKERR mField.add_ents_to_finite_element_by_type(contact_range, MBTRI,
                                                   contactElement);
  CHKERR add_field_to_fe(contactElement, piolaStress);
  CHKERR add_field_to_fe(contactElement, spatialH1Disp);
  CHKERR add_field_to_fe(contactElement, contactDisp);
  CHKERR mField.build_finite_elements(contactElement);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::addDMs(const BitRefLevel bit) {
  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, essentialBcElement);
  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_essential_spatial_stress_boundary =
      [&](const std::string prb_name) {
        MoFEMFunctionBegin;
        // for (int d : {0, 1, 2})
        //   CHKERR
        //   mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
        //       prb_name, piolaStress, (*bcSpatialFreeTraction)[d], d, d, 0,
        //       MAX_DOFS_ON_ENTITY, NOISY, true);
        Range all_three;
        for (int d : {0, 1, 2})
          all_three.merge((*bcSpatialFreeTraction)[d]);
        for (int d : {0, 1, 2})
          all_three = intersect(all_three, (*bcSpatialFreeTraction)[d]);
        CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
            prb_name, piolaStress, all_three, 0, SPACE_DIM, 0,
            MAX_DOFS_ON_ENTITY, NOISY, true);
        MoFEMFunctionReturn(0);
      };
  CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
 
  auto contact_range = get_range_from_block(mField, "CONTACT");
 
  // 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);
  if (!noStretch) {
    CHKERR DMMoFEMAddSubFieldRow(dmElastic, stretchTensor);
  }
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, contactDisp);
  CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName);
  CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement);
  CHKERR DMMoFEMAddElement(dmElastic, essentialBcElement);
  CHKERR DMMoFEMAddElement(dmElastic, contactElement);
  CHKERR DMMoFEMAddElement(dmElastic, skinElement);
  CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
  CHKERR DMSetUp(dmElastic);
 
  auto set_zero_block = [&]() {
    MoFEMFunctionBegin;
    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);
    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);
    MoFEMFunctionReturn(0);
  };
 
  auto set_section = [&]() {  
    MoFEMFunctionBegin;
    PetscSection section;
    CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
                                                           &section);
    CHKERR DMSetSection(dmElastic, section);
    CHKERR DMSetGlobalSection(dmElastic, section);
    CHKERR PetscSectionDestroy(&section);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_zero_block();
  CHKERR set_section();
 
  dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
  CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, spatialH1Disp);
  CHKERR DMMoFEMAddElement(dmPrjSpatial, elementVolumeName);
  CHKERR DMMoFEMSetSquareProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMSetUp(dmPrjSpatial);
 
  CHKERR mField.getInterface<BcManager>()
      ->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
          "PROJECT_SPATIAL", spatialH1Disp, true, false);
 
  MoFEMFunctionReturn(0);
}
 
BcDisp::BcDisp(std::string name, std::vector<double> &attr, Range &faces)
    : blockName(name), faces(faces) {
  vals.resize(3, false);
  flags.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    vals[ii] = attr[ii];
    flags[ii] = static_cast<int>(attr[ii + 3]);
  }
 
  MOFEM_LOG("EP", Sev::inform) << "Add BCDisp " << name;
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCDisp vals " << vals[0] << " " << vals[1] << " " << vals[2];
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCDisp flags " << flags[0] << " " << flags[1] << " " << flags[2];
  MOFEM_LOG("EP", Sev::inform) << "Add BCDisp nb. of faces " << faces.size();
}
 
BcRot::BcRot(std::string name, std::vector<double> &attr, Range &faces)
    : blockName(name), faces(faces) {
  vals.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    vals[ii] = attr[ii];
  }
  theta = attr[3];
}
 
TractionBc::TractionBc(std::string name, std::vector<double> &attr,
                       Range &faces)
    : blockName(name), faces(faces) {
  vals.resize(3, false);
  flags.resize(3, false);
  for (int ii = 0; ii != 3; ++ii) {
    vals[ii] = attr[ii];
    flags[ii] = static_cast<int>(attr[ii + 3]);
  }
 
  MOFEM_LOG("EP", Sev::inform) << "Add BCForce " << name;
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCForce vals " << vals[0] << " " << vals[1] << " " << vals[2];
  MOFEM_LOG("EP", Sev::inform)
      << "Add BCForce flags " << flags[0] << " " << flags[1] << " " << flags[2];
  MOFEM_LOG("EP", Sev::inform) << "Add BCForce nb. of faces " << faces.size();
}
 
MoFEMErrorCode
EshelbianCore::getTractionFreeBc(const EntityHandle meshset,
                                 boost::shared_ptr<TractionFreeBc> &bc_ptr,
                                 const std::string contact_set_name) {
  MoFEMFunctionBegin;
 
  // get skin from all tets
  Range tets;
  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
  Range tets_skin_part;
  Skinner skin(&mField.get_moab());
  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
  Range tets_skin;
  CHKERR pcomm->filter_pstatus(tets_skin_part,
                               PSTATUS_SHARED | PSTATUS_MULTISHARED,
                               PSTATUS_NOT, -1, &tets_skin);
 
  bc_ptr->resize(3);
  for (int dd = 0; dd != 3; ++dd)
    (*bc_ptr)[dd] = tets_skin;
 
  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);
    }
 
  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);
    }
 
  // 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);
  }
 
  // for (int dd = 0; dd != 3; ++dd) {
  //   EntityHandle meshset;
  //   CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
  //   CHKERR mField.get_moab().add_entities(meshset, (*bc_ptr)[dd]);
  //   std::string file_name = disp_block_set_name +
  //       "_traction_free_bc_" + boost::lexical_cast<std::string>(dd) + ".vtk";
  //   CHKERR mField.get_moab().write_file(file_name.c_str(), " VTK ", "",
  //                                       &meshset, 1);
  //   CHKERR mField.get_moab().delete_entities(&meshset, 1);
  // }
 
  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 {
    if (EshelbianCore::gradApperoximator > MODERATE_ROT)
      return p_data + p_data + (p_data - 1);
    else
      return 2 * p_data;
  }
};
 
struct FaceRule {
  int operator()(int p_row, int p_col, int p_data) const { return 2 * p_data; }
};
 
struct CGGUserPolynomialBase : public BaseFunction {
 
  CGGUserPolynomialBase() {}
  ~CGGUserPolynomialBase() {}
 
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 BaseFunctionUnknownInterface **iface) const {
    *iface = const_cast<CGGUserPolynomialBase *>(this);
    return 0;
  }
 
  MoFEMErrorCode getValue(MatrixDouble &pts,
                          boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
    MoFEMFunctionBeginHot;
 
    cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
 
    int nb_gauss_pts = pts.size2();
    if (!nb_gauss_pts) {
      MoFEMFunctionReturnHot(0);
    }
 
    if (pts.size1() < 3) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Wrong dimension of pts, should be at least 3 rows with "
              "coordinates");
    }
 
    switch (cTx->sPace) {
    case HDIV:
      CHKERR getValueHdivForCGGBubble(pts);
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
    }
 
    MoFEMFunctionReturnHot(0);
  }
 
private:
  EntPolynomialBaseCtx *cTx;
 
  MatrixDouble shapeFun;
  boost::tuple<int, int, MatrixDouble> cachePhi = {0, 0, MatrixDouble()};
 
  MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts) {
    MoFEMFunctionBegin;
 
    // This should be used only in case USER_BASE is selected
    if (cTx->bAse != USER_BASE) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Wrong base, should be USER_BASE");
    }
    // get access to data structures on element
    EntitiesFieldData &data = cTx->dAta;
    // Get approximation order on element. Note that bubble functions are only
    // on tetrahedron.
    const int order = data.dataOnEntities[MBTET][0].getOrder();
    /// number of integration points
    const int nb_gauss_pts = pts.size2();
 
    if (cachePhi.get<0>() == order && cachePhi.get<1>() == nb_gauss_pts) {
      auto &mat = cachePhi.get<2>();
      auto &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
      phi.resize(mat.size1(), mat.size2(), false);
      noalias(phi) = mat;
    } else {
      // calculate shape functions, i.e. barycentric coordinates
      shapeFun.resize(nb_gauss_pts, 4, false);
      CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
                        &pts(2, 0), nb_gauss_pts);
      // direvatives of shape functions
      double diff_shape_fun[12];
      CHKERR ShapeDiffMBTET(diff_shape_fun);
 
      const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
      // get base functions and set size
      MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
      phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
      // finally calculate base functions
      FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_phi(
          &phi(0, 0), &phi(0, 1), &phi(0, 2),
 
          &phi(0, 3), &phi(0, 4), &phi(0, 5),
 
          &phi(0, 6), &phi(0, 7), &phi(0, 8));
      CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
                                  nb_gauss_pts);
 
      cachePhi.get<0>() = order;
      cachePhi.get<1>() = nb_gauss_pts;
      cachePhi.get<2>().resize(phi.size1(), phi.size2(), false);
      noalias(cachePhi.get<2>()) = phi;
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
MoFEMErrorCode EshelbianCore::setBaseVolumeElementOps(
    const int tag, const bool do_rhs, const bool do_lhs,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe) {
  MoFEMFunctionBegin;
  fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
 
  fe->getUserPolynomialBase() =
      boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(fe->getOpPtrVector(),
                                                 {HDIV, H1, L2});
 
  // set integration rule
  fe->getRuleHook = VolRule();
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
    dataAtPts->physicsPtr = physicalEquations;
  }
 
  // calculate fields values
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
      piolaStress, dataAtPts->getApproxPAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
      bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
      piolaStress, dataAtPts->getDivPAtPts()));
 
  if (noStretch) {
    fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
        tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
    fe->getOpPtrVector().push_back(new OpCalculateStretchFromStress(dataAtPts));
  } else {
    fe->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValues<3>(
            stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
  }
 
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
 
  // velocities
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
      spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
  if (!noStretch) {
    fe->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValuesDot<3>(
            stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
  }
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
      rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
 
  // acceleration
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
        spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
  }
 
  // H1 displacements
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
      spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
 
  // calculate other derived quantities
  fe->getOpPtrVector().push_back(
      new OpCalculateRotationAndSpatialGradient(dataAtPts));
 
  // evaluate integration points
  if (!noStretch) {
    fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
        tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setVolumeElementOps(
    const int tag, const bool add_elastic, const bool add_material,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
  MoFEMFunctionBegin;
 
  auto time_scale = boost::make_shared<TimeScale>();
 
  // Right hand side
  CHKERR setBaseVolumeElementOps(tag, true, false, 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));
    if (!noStretch) {
      fe_rhs->getOpPtrVector().push_back(
          new OpSpatialPhysical(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));
 
    // Body forces
    using BodyNaturalBC =
        NaturalBC<VolumeElementForcesAndSourcesCore::UserDataOperator>::
            Assembly<PETSC>::LinearForm<GAUSS>;
    using OpBodyForce =
        BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
    CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
        fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {time_scale},
        "BODY_FORCE", Sev::inform);
  }
 
  // Left hand side
  CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
 
  // elastic
  if (add_elastic) {
 
    const bool symmetric_system =
        (gradApperoximator <= MODERATE_ROT && rotSelector == SMALL_ROT);
 
    if(noStretch) {
      fe_lhs->getOpPtrVector().push_back(
          new OpSpatialConsistency_dP_dP(piolaStress, piolaStress, dataAtPts));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dP(
          bubbleField, piolaStress, dataAtPts));
      fe_lhs->getOpPtrVector().push_back(
          new OpSpatialConsistency_dBubble_dBubble(bubbleField, bubbleField,
                                                   dataAtPts));
    } else {
      fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_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));
      if (!symmetric_system) {
        fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
            stretchTensor, rotAxis, dataAtPts, 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, symmetric_system));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
        bubbleField, rotAxis, dataAtPts, symmetric_system));
 
    if (!symmetric_system) {
      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));
    }
 
    // Stabilisation
    if constexpr (A == AssemblyType::BLOCK_SCHUR) {
      using B = FormsIntegrators<
          VolumeElementForcesAndSourcesCore::UserDataOperator>::
          Assembly<BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<GAUSS>;
      using OpMassStab = B::OpMassCache<1, SPACE_DIM>;
      if (precEpsOmega > std::numeric_limits<double>::epsilon()) {
        auto cache = boost::make_shared<MatrixDouble>();
        CacheMatsTypeType cache_mats = {{{MBTET, MBTET}, cache}};
        fe_lhs->getOpPtrVector().push_back(
            new OpMassStab(cache_mats, rotAxis, rotAxis, precEpsOmega));
      }
      if (std::abs(alphaRho + alphaW) <
          std::numeric_limits<double>::epsilon()) {
        auto cache = boost::make_shared<MatrixDouble>();
        CacheMatsTypeType cache_mats = {{{MBTET, MBTET}, cache}};
        fe_lhs->getOpPtrVector().push_back(
            new OpMassStab(cache_mats, spatialL2Disp, spatialL2Disp, precEpsW));
      }
    }
 
    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 = FaceRule();
  fe_lhs->getRuleHook = FaceRule();
 
  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(fe_rhs->getOpPtrVector(),
                                                     {HDIV});
  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(fe_lhs->getOpPtrVector(),
                                                     {HDIV});
 
  if (add_elastic) {
 
    fe_rhs->getOpPtrVector().push_back(
        new OpDispBc(piolaStress, dataAtPts, bcSpatialDispVecPtr,
                     {
 
                         boost::make_shared<TimeScale>("disp_history.txt")
 
                     }));
    fe_rhs->getOpPtrVector().push_back(
        new OpRotationBc(piolaStress, dataAtPts, bcSpatialRotationVecPtr));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setContactElementOps(
 
    boost::shared_ptr<ContactTree> &fe_contact_tree,
    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_face_rhs,
    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_face_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 get_map_skin = [this](auto &&map) {
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
 
    Skinner skin(&mField.get_moab());
    for (auto &m : map) {
      Range skin_faces;
      CHKERR skin.find_skin(0, m.second, false, skin_faces);
      CHK_MOAB_THROW(pcomm->filter_pstatus(skin_faces,
                                           PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                           PSTATUS_NOT, -1, nullptr),
                     "filter");
      m.second.swap(skin_faces);
    }
    return map;
  };
 
  /* The above code is written in C++ and it appears to be defining and using
  various operations on boundary elements and side elements. */
  using BoundaryEle = FaceElementForcesAndSourcesCore;
  using BoundaryEleOp = BoundaryEle::UserDataOperator;
 
  fe_face_rhs = boost::make_shared<BoundaryEle>(mField);
  fe_face_lhs = boost::make_shared<BoundaryEle>(mField);
  auto rule = [](int, int, int o) { return -1; };
 
  int levels = 0;
  CHKERR PetscOptionsGetInt(PETSC_NULL, "-contact_max_post_proc_ref_level",
                            &levels, PETSC_NULL);
 
  auto refine = Tools::refineTriangle(levels);
 
  auto set_rule = [levels, refine](
 
                      ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                      int order_col, int order_data
 
                  ) {
    MoFEMFunctionBegin;
    auto rule = 2 * order_data;
    fe_raw_ptr->gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
    MoFEMFunctionReturn(0);
  };
 
  fe_face_rhs->getRuleHook = rule;
  fe_face_lhs->getRuleHook = rule;
  fe_face_rhs->setRuleHook = set_rule;
  fe_face_lhs->setRuleHook = set_rule;
 
  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      fe_face_rhs->getOpPtrVector(), {HDIV});
  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      fe_face_lhs->getOpPtrVector(), {HDIV});
 
  auto add_contact_three = [&]() {
    MoFEMFunctionBegin;
    auto tree_moab_ptr = boost::make_shared<moab::Core>();
    fe_contact_tree = boost::make_shared<ContactTree>(
        mField, tree_moab_ptr, spaceOrder,
        get_map_skin(get_body_range("BODY", 3)));
    auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
    AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        fe_contact_tree->getOpPtrVector(), {HDIV});
    fe_contact_tree->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(
            contactDisp, contact_common_data_ptr->contactDispPtr()));
    fe_contact_tree->getOpPtrVector().push_back(
        new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
            piolaStress, contact_common_data_ptr->contactTractionPtr()));
 
    auto u_h1_ptr = boost::make_shared<MatrixDouble>();
    fe_contact_tree->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
    fe_contact_tree->getOpPtrVector().push_back(
        new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
    MoFEMFunctionReturn(0);
  };
 
  auto add_ops_face_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
    auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
    pip.push_back(new OpCalculateVectorFieldValues<3>(
        contactDisp, contact_common_data_ptr->contactDispPtr()));
    pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
        piolaStress, contact_common_data_ptr->contactTractionPtr()));
    auto u_h1_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
    pip.push_back(new OpTreeSearch(
        fe_contact_tree, contact_common_data_ptr, u_h1_ptr,
        get_range_from_block(mField, "CONTACT"), nullptr, nullptr));
 
    // get body id and SDF range
    auto contact_sfd_map_range_ptr = boost::make_shared<std::map<int, Range>>(
        get_body_range("CONTACT_SDF", 2));
 
    pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU(
        contactDisp, contactDisp, contact_common_data_ptr, fe_contact_tree,
        contact_sfd_map_range_ptr));
    pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP(
        contactDisp, piolaStress, contact_common_data_ptr, fe_contact_tree,
        contact_sfd_map_range_ptr));
    pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU(
        piolaStress, contactDisp, contact_common_data_ptr, fe_contact_tree));
 
    // Stabilisation
    if constexpr (A == AssemblyType::BLOCK_SCHUR) {
      using B = FormsIntegrators<BoundaryEle::UserDataOperator>::Assembly<
          BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<GAUSS>;
      using OpMassStab = B::OpMass<1, SPACE_DIM>;
      if (precEpsContactDisp > std::numeric_limits<double>::epsilon()) {
        pip.push_back(new OpMassStab(
            contactDisp, contactDisp,
            [this](double, double, double) { return precEpsContactDisp; }));
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_ops_face_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
    auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
    pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
        piolaStress, contact_common_data_ptr->contactTractionPtr()));
    pip.push_back(new OpCalculateVectorFieldValues<3>(
        contactDisp, contact_common_data_ptr->contactDispPtr()));
    auto u_h1_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
    pip.push_back(new OpTreeSearch(
        fe_contact_tree, contact_common_data_ptr, u_h1_ptr,
        get_range_from_block(mField, "CONTACT"), nullptr, nullptr));
    auto contact_sfd_map_range_ptr = boost::make_shared<std::map<int, Range>>(
        get_body_range("CONTACT_SDF", 2));
    pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Rhs(
        contactDisp, contact_common_data_ptr, fe_contact_tree,
        contact_sfd_map_range_ptr));
    pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryHDivRhs(
        piolaStress, contact_common_data_ptr, fe_contact_tree));
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_contact_three();
  CHKERR add_ops_face_lhs(fe_face_lhs->getOpPtrVector());
  CHKERR add_ops_face_rhs(fe_face_rhs->getOpPtrVector());
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setElasticElementOps(const int tag) {
  MoFEMFunctionBegin;
  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;
  };
 
  auto op_contact_bc = get_op_contact_bc();
  elasticFeLhs->getOpPtrVector().push_back(op_contact_bc);
 
  CHKERR setContactElementOps(
 
      contactTreeRhs, contactRhs, op_contact_bc->getSideFEPtr()
 
  );
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setElasticElementToTs(DM dm) {
  MoFEMFunctionBegin;
  boost::shared_ptr<FEMethod> null;
  boost::shared_ptr<FeTractionBc> spatial_traction_bc(
      new FeTractionBc(mField, piolaStress, bcSpatialTraction));
 
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
 
    CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
                                  null);
    CHKERR DMMoFEMTSSetI2Function(dm, elementVolumeName, elasticFeRhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Function(dm, naturalBcElement, elasticBcRhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Function(dm, contactElement, contactRhs, null, null);
    CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, null, null,
                                  spatial_traction_bc);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
                                  null);
 
  } else {
    CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
                                 null);
    CHKERR DMMoFEMTSSetIFunction(dm, elementVolumeName, elasticFeRhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIFunction(dm, naturalBcElement, elasticBcRhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIFunction(dm, contactElement, contactRhs, null, null);
    CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, null,
                                 spatial_traction_bc);
    CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
                                 null);
  }
 
  auto set_set_essential = [&]() {
    MoFEMFunctionBegin;
    auto fe_rhs_ptr = boost::make_shared<FEMethod>();
    auto fe_lhs_ptr = boost::make_shared<FEMethod>();
 
    auto get_dofs = [&](auto fe_ptr) {
      auto prb_ptr = getProblemPtr(dmElastic);
 
      auto vec_glob_dofs =
          boost::shared_ptr<std::vector<int>>(fe_ptr, new std::vector<int>());
      const auto field_id = mField.get_field_bit_number(piolaStress);
      auto numered_dofs_ptr = prb_ptr->getNumeredRowDofsPtr();
 
      vec_glob_dofs->reserve(numered_dofs_ptr->size());
      for (auto d = 0; d != SPACE_DIM; ++d) {
        for (auto pit = (*bcSpatialFreeTraction)[d].pair_begin();
             pit != (*bcSpatialFreeTraction)[d].pair_end(); ++pit) {
          auto lo = numered_dofs_ptr->get<Unique_mi_tag>().lower_bound(
              DofEntity::getLoFieldEntityUId(field_id, pit->first));
          auto hi = numered_dofs_ptr->get<Unique_mi_tag>().upper_bound(
              DofEntity::getHiFieldEntityUId(field_id, pit->second));
          for (; lo != hi; ++lo) {
            if ((*lo)->getDofCoeffIdx() == d) {
              vec_glob_dofs->push_back((*lo)->getPetscGlobalDofIdx());
            }
          }
        }
      }
 
      return vec_glob_dofs;
    };
 
    fe_rhs_ptr->postProcessHook = [this, fe_rhs_ptr]() {
      MoFEMFunctionBegin;
 
      (*fe_rhs_ptr->vecAssembleSwitch) = false;
      CHKERR VecGhostUpdateBegin(fe_rhs_ptr->f, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(fe_rhs_ptr->f, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecAssemblyBegin(fe_rhs_ptr->f);
      CHKERR VecAssemblyEnd(fe_rhs_ptr->f);
 
      double *a;
      CHKERR VecGetArray(fe_rhs_ptr->f, &a);
 
      const auto field_id = mField.get_field_bit_number(piolaStress);
      auto numered_dofs_ptr = fe_rhs_ptr->problemPtr->getNumeredRowDofsPtr();
 
      for (auto d = 0; d != SPACE_DIM; ++d) {
        for (auto pit = (*bcSpatialFreeTraction)[d].pair_begin();
             pit != (*bcSpatialFreeTraction)[d].pair_end(); ++pit) {
          auto lo = numered_dofs_ptr->get<Unique_mi_tag>().lower_bound(
              DofEntity::getLoFieldEntityUId(field_id, pit->first));
          auto hi = numered_dofs_ptr->get<Unique_mi_tag>().upper_bound(
              DofEntity::getHiFieldEntityUId(field_id, pit->second));
          for (; lo != hi; ++lo) {
            if ((*lo)->getDofCoeffIdx() == d) {
              a[(*lo)->getPetscLocalDofIdx()] = 0.0;
            }
          }
        }
      }
 
      CHKERR VecRestoreArray(fe_rhs_ptr->f, &a);
 
      MoFEMFunctionReturn(0);
    };
 
    auto glob_dofs = get_dofs(fe_lhs_ptr);
    fe_lhs_ptr->postProcessHook = [this, fe_lhs_ptr, glob_dofs]() {
      MoFEMFunctionBegin;
 
      (*fe_lhs_ptr->matAssembleSwitch) = false;
      CHKERR MatAssemblyBegin(fe_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
      CHKERR MatAssemblyEnd(fe_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
 
      CHKERR MatZeroRowsColumns(fe_lhs_ptr->B, glob_dofs->size(),
                                &*glob_dofs->begin(), 1.0, PETSC_NULL,
                                PETSC_NULL);
 
      MoFEMFunctionReturn(0);
    };
 
    auto ts_ctx_ptr = getDMTsCtx(dmElastic);
    ts_ctx_ptr->getPostProcessIFunction().push_back(fe_rhs_ptr);
    ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_lhs_ptr);
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_set_essential();
 
  MoFEMFunctionReturn(0);
}
 
#include "impl/EshelbianMonitor.cpp"
#include "impl/SetUpSchurImpl.cpp"
#include "impl/TSElasticPostStep.cpp"
 
MoFEMErrorCode EshelbianCore::solveElastic(TS ts, Vec x) {
  MoFEMFunctionBegin;
 
#ifdef PYTHON_SDF
 
  boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
 
  auto file_exists = [](std::string myfile) {
    std::ifstream file(myfile.c_str());
    if (file) {
      return true;
    }
    return false;
  };
 
  if (file_exists("sdf.py")) {
    MOFEM_LOG("EP", Sev::inform) << "sdf.py file found";
    sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
    CHKERR sdf_python_ptr->sdfInit("sdf.py");
    ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
  } else {
    MOFEM_LOG("EP", Sev::warning) << "sdf.py file NOT found";
  }
#else
#endif
 
  boost::shared_ptr<TsCtx> ts_ctx;
  CHKERR DMMoFEMGetTsCtx(dmElastic, ts_ctx);
  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
 
  boost::shared_ptr<FEMethod> monitor_ptr(new EshelbianMonitor(*this));
  ts_ctx->getLoopsMonitor().push_back(
      TsCtx::PairNameFEMethodPtr(elementVolumeName, monitor_ptr));
 
  CHKERR TSAppendOptionsPrefix(ts, "elastic_");
  CHKERR TSSetFromOptions(ts);
 
  CHKERR TSSetDM(ts, dmElastic);
 
  SNES snes;
  CHKERR TSGetSNES(ts, &snes);
 
  PetscViewerAndFormat *vf;
  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
                                    PETSC_VIEWER_DEFAULT, &vf);
  CHKERR SNESMonitorSet(
      snes,
      (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
      vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
 
  PetscSection section;
  CHKERR DMGetSection(dmElastic, &section);
  int num_fields;
  CHKERR PetscSectionGetNumFields(section, &num_fields);
  for (int ff = 0; ff != num_fields; ff++) {
    const char *field_name;
    CHKERR PetscSectionGetFieldName(section, ff, &field_name);
    MOFEM_LOG_C("EP", Sev::inform, "Field %d name %s", ff, field_name);
  }
 
  CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
 
  // Adding field split solver
  boost::shared_ptr<SetUpSchur> schur_ptr;
  if constexpr (A == AssemblyType::BLOCK_SCHUR) {
    schur_ptr = SetUpSchur::createSetUpSchur(mField, &*this);
    CHKERR schur_ptr->setUp(ts);
  }
 
  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::swichCacheHDivBaseDemkowiczOn(
      {elasticFeLhs.get(), elasticFeRhs.get()});
  CHKERR TSSetUp(ts);
  CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
  CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  CHKERR TSElasticPostStep::postStepInitialise(this);
  CHKERR TSSolve(ts, PETSC_NULL);
  CHKERR TSElasticPostStep::postStepDestroy();
  TetPolynomialBase::swichCacheHDivBaseDemkowiczOff(
      {elasticFeLhs.get(), elasticFeRhs.get()});
 
  // CHKERR TSGetSNES(ts, &snes);
  int lin_solver_iterations;
  CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
  MOFEM_LOG("EP", Sev::inform)
      << "Number of linear solver iterations " << lin_solver_iterations;
 
  PetscBool test_cook_flg = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_cook", &test_cook_flg,
                             PETSC_NULL);
  if (test_cook_flg) {
    constexpr int expected_lin_solver_iterations = 11;
    // if (lin_solver_iterations != expected_lin_solver_iterations)
    //   SETERRQ2(
    //       PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
    //       "Expected number of iterations is different than expected %d !=
    //       %d", lin_solver_iterations, expected_lin_solver_iterations);
  }
 
  CHKERR gettingNorms();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::postProcessResults(const int tag,
                                                 const std::string file) {
  MoFEMFunctionBegin;
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
  }
 
  PostProcBrokenMeshInMoab<FaceElementForcesAndSourcesCore> post_proc(mField);
  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(post_proc.getOpPtrVector(),
                                                     {HDIV});
 
  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
 
  auto domain_ops = [&](auto &fe) {
    MoFEMFunctionBegin;
    fe.getUserPolynomialBase() =
        boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(fe.getOpPtrVector(),
                                                   {HDIV, H1, L2});
    fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
        piolaStress, dataAtPts->getApproxPAtPts()));
    fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
        bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
    fe.getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
        stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
        spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
    // evaluate derived quantities
    fe.getOpPtrVector().push_back(
        new OpCalculateRotationAndSpatialGradient(dataAtPts));
 
    // evaluate integration points
    fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
        tag, true, false, dataAtPts, physicalEquations));
 
    // contact
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
        spatialL2Disp, contact_common_data_ptr->contactDispPtr()));
 
    MoFEMFunctionReturn(0);
  };
 
  auto op_loop_side = new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
      mField, elementVolumeName, SPACE_DIM);
  CHKERR domain_ops(*(op_loop_side->getSideFEPtr()));
  post_proc.getOpPtrVector().push_back(op_loop_side);
  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      contactDisp, dataAtPts->getContactL2AtPts()));
  post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
      post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), dataAtPts));
 
  // contact
  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
  auto lambda_h1_ptr = boost::make_shared<MatrixDouble>();
  post_proc.getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
 
  using BoundaryEle = FaceElementForcesAndSourcesCore;
  using BoundaryEleOp = BoundaryEle::UserDataOperator;
 
  post_proc.getOpPtrVector().push_back(
      new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
          piolaStress, contact_common_data_ptr->contactTractionPtr()));
 
  post_proc.getOpPtrVector().push_back(new OpTreeSearch(
      contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
      get_range_from_block(mField, "CONTACT"), &post_proc.getPostProcMesh(),
      &post_proc.getMapGaussPts()));
 
  CHKERR DMoFEMLoopFiniteElements(dM, contactElement, contactTreeRhs);
  CHKERR DMoFEMLoopFiniteElements(dM, skinElement.c_str(), &post_proc);
  CHKERR post_proc.writeFile(file.c_str());
  MoFEMFunctionReturn(0);
}
 
//! [Getting norms]
MoFEMErrorCode EshelbianCore::gettingNorms() {
  MoFEMFunctionBegin;
 
  auto post_proc_norm_fe =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
 
  auto post_proc_norm_rule_hook = [](int, int, int p) -> int {
    return 2 * (p);
  };
  post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
 
  post_proc_norm_fe->getUserPolynomialBase() =
      boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV});
 
  enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
  auto norms_vec =
      createVectorMPI(mField.get_comm(), LAST_NORM, PETSC_DETERMINE);
  CHKERR VecZeroEntries(norms_vec);
 
  auto u_l2_ptr = boost::make_shared<MatrixDouble>();
  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp, u_l2_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(spatialH1Disp, u_h1_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_NORM_L2));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_h1_ptr, norms_vec, U_NORM_H1));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_ERROR_L2,
                                         u_h1_ptr));
 
  auto piola_ptr = boost::make_shared<MatrixDouble>();
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateHVecTensorField<3, 3>(piolaStress, piola_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor2<3, 3>(piola_ptr, norms_vec, PIOLA_NORM));
 
  TetPolynomialBase::swichCacheHDivBaseDemkowiczOn({post_proc_norm_fe.get()});
  CHKERR mField.loop_finite_elements("ELASTIC_PROBLEM", elementVolumeName,
                                     *post_proc_norm_fe);
  TetPolynomialBase::swichCacheHDivBaseDemkowiczOff({post_proc_norm_fe.get()});
 
  CHKERR VecAssemblyBegin(norms_vec);
  CHKERR VecAssemblyEnd(norms_vec);
  const double *norms;
  CHKERR VecGetArrayRead(norms_vec, &norms);
  MOFEM_LOG("EP", Sev::inform) << "norm_u: " << std::sqrt(norms[U_NORM_L2]);
  MOFEM_LOG("EP", Sev::inform) << "norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
  MOFEM_LOG("EP", Sev::inform)
      << "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
  MOFEM_LOG("EP", Sev::inform)
      << "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
  CHKERR VecRestoreArrayRead(norms_vec, &norms);
 
  MoFEMFunctionReturn(0);
}
//! [Getting norms]
 
MoFEMErrorCode EshelbianCore::getSpatialDispBc() {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      "", piolaStress, false, false);
 
  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    if (auto disp_bc = bc.second->dispBcPtr) {
 
      MOFEM_LOG("EP", Sev::noisy) << *disp_bc;
 
      std::vector<double> block_attributes(6, 0.);
      if (disp_bc->data.flag1 == 1) {
        block_attributes[0] = disp_bc->data.value1;
        block_attributes[3] = 1;
      }
      if (disp_bc->data.flag2 == 1) {
        block_attributes[1] = disp_bc->data.value2;
        block_attributes[4] = 1;
      }
      if (disp_bc->data.flag3 == 1) {
        block_attributes[2] = disp_bc->data.value3;
        block_attributes[5] = 1;
      }
      auto faces = bc.second->bcEnts.subset_by_dimension(2);
      bcSpatialDispVecPtr->emplace_back(bc.first, block_attributes, faces);
    }
  }
 
  // old way of naming blocksets for displacement BCs
  CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::getSpatialTractionBc() {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<ForceCubitBcData>("", piolaStress,
                                                          false, false);
 
  bcSpatialTraction = boost::make_shared<TractionBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    if (auto force_bc = bc.second->forceBcPtr) {
      std::vector<double> block_attributes(6, 0.);
      block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
      block_attributes[3] = 1;
      block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
      block_attributes[4] = 1;
      block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
      block_attributes[5] = 1;
      auto faces = bc.second->bcEnts.subset_by_dimension(2);
      bcSpatialTraction->emplace_back(bc.first, block_attributes, faces);
    }
  }
 
  CHKERR getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
 
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity
 
#include <impl/EshelbianContact.cpp>