EshelbianPlasticity.cpp

Eshelbian plasticity implementation.

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_SFD
#include <boost/python.hpp>
#include <boost/python/def.hpp>
namespace bp = boost::python;
 
#pragma message "With PYTHON_SFD"
 
#else
 
#pragma message "Without PYTHON_SFD"
 
#endif
 
#include <EshelbianContact.hpp>
 
namespace EshelbianPlasticity {
 
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"), materialL2Disp("W"),
      streachTensor("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"};
  PetscInt choice_rot = rotSelector;
 
  const char *list_stretches[] = {"linear", "log"};
  PetscInt choice_stretch = StretchSelector::LOG;
 
  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
                           "none");
 
  spaceOrder = 1;
  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
                         spaceOrder, &spaceOrder, PETSC_NULL);
  materialOrder = 1;
  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
                         "", materialOrder, &materialOrder, PETSC_NULL);
 
  alphaU = 0;
  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
                            &alphaU, PETSC_NULL);
 
  alphaW = 0;
  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
                            &alphaW, PETSC_NULL);
 
  alphaRho = 0;
  CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
                            &alphaRho, PETSC_NULL);
 
  precEps = 0;
  CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
                            precEps, &precEps, PETSC_NULL);
 
  precEpsOmega = 0;
  CHKERR PetscOptionsScalar("-preconditioner_eps_omega", "preconditioner_eps",
                            "", precEpsOmega, &precEpsOmega, PETSC_NULL);
  precEpsW = 0;
  CHKERR PetscOptionsScalar("-preconditioner_eps_w", "preconditioner_eps", "",
                            precEpsW, &precEpsW, PETSC_NULL);
 
  CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
                           LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
                           PETSC_NULL);
  CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
                            &exponentBase, PETSC_NULL);
  CHKERR PetscOptionsEList(
      "-stretches", "stretches", "", list_stretches, StretchSelector::LOG + 1,
      list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
 
  ierr = PetscOptionsEnd();
  CHKERRG(ierr);
 
  rotSelector = static_cast<RotSelector>(choice_rot);
  switch (choice_stretch) {
  case StretchSelector::LINEAR:
    f = f_linear;
    d_f = d_f_linear;
    dd_f = dd_f_linear;
    break;
  case StretchSelector::LOG:
    f = f_log;
    d_f = d_f_log;
    dd_f = 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) << "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[rotSelector];
  if (exponentBase != exp(1))
    MOFEM_LOG("EP", Sev::inform) << "Base exponent " << exponentBase;
  else
    MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
  MOFEM_LOG("EP", Sev::inform) << "Stretch " << list_stretches[choice_stretch];
 
  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 msh_mng = mField.getInterface<MeshsetsManager>();
        auto reg_exp = std::regex((boost::format("%s(.*)") % block_name).str());
        for (auto m : msh_mng->getCubitMeshsetPtr(reg_exp)) {
          Range faces;
          CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(),
                                                    SPACE_DIM - 1, faces, true);
          tets_skin = subtract(tets_skin, faces);
          MOFEM_LOG("EP", Sev::inform)
              << "Subtracting " << m->getName() << " " << faces.size();
        }
        MoFEMFunctionReturn(0);
      };
  CHKERR subtract_faces_where_displacements_are_applied("CONTACT");
 
  Range faces;
  CHKERR mField.get_moab().get_adjacencies(tets, 2, 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);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    CHKERR mField.set_field_order(faces_not_on_the_skin, field_name, order);
    CHKERR mField.set_field_order(tets_skin, field_name, 0);
    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_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(streachTensor, spaceOrder, 6);
 
  // 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);
 
  // Add history filedes
  // CHKERR add_l2_field(materialGradient + "0", materialOrder - 1, 9);
 
  // spatial displacement
  CHKERR add_h1_field(spatialH1Disp, spaceOrder, 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, streachTensor);
    CHKERR add_field_to_fe(elementVolumeName, rotAxis);
    CHKERR add_field_to_fe(elementVolumeName, spatialL2Disp);
    CHKERR add_field_to_fe(elementVolumeName, streachTensor);
    CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
    // 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;
 
  auto bc_elements_add_to_range = [&](const std::string block_name, Range &r) {
    MoFEMFunctionBegin;
 
    auto mesh_mng = mField.getInterface<MeshsetsManager>();
    auto bcs = mesh_mng->getCubitMeshsetPtr(
 
        std::regex((boost::format("%s(.*)") % block_name).str())
 
    );
 
    for (auto bc : bcs) {
      Range faces;
      CHKERR bc->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                                 true);
      r.merge(faces);
    }
 
    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);
    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, eshelbyStress);
  CHKERR mField.build_finite_elements(skinElement);
 
  Range contact_range;
  CHKERR bc_elements_add_to_range("CONTACT", contact_range);
  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 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);
        MoFEMFunctionReturn(0);
      };
  CHKERR remove_dofs_on_essential_spatial_stress_boundary("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.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, streachTensor.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp.c_str());
  CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName);
  CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement);
  CHKERR DMMoFEMAddElement(dmElastic, essentialBcElement);
  CHKERR DMMoFEMAddElement(dmElastic, contactElement);
  CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
  CHKERR DMSetUp(dmElastic);
 
  {
    PetscSection section;
    CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
                                                           &section);
    CHKERR DMSetSection(dmElastic, section);
    CHKERR DMSetGlobalSection(dmElastic, section);
    CHKERR PetscSectionDestroy(&section);
  }
 
  dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
  CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, spatialH1Disp.c_str());
  CHKERR DMMoFEMAddElement(dmPrjSpatial, elementVolumeName);
  CHKERR DMMoFEMSetSquareProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMSetUp(dmPrjSpatial);
 
  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 {
    return p_data + p_data + (p_data - 1);
  }
};
 
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;
 
  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();
    // number of base functions
 
    // 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);
 
    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()));
  // fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
  //     eshelbyStress, dataAtPts->getApproxSigmaAtPts()));
  // fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
  //     eshelbyStress, dataAtPts->getDivSigmaAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
      streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  // fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
  //     materialGradient, dataAtPts->getBigGAtPts(), MBTET));
  // fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
  //     materialGradient + "0", dataAtPts->getBigG0AtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialL2Disp, dataAtPts->getSmallWAtPts(), MBTET));
 
  // velocities
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
      spatialL2Disp, dataAtPts->getSmallWDotAtPts(), MBTET));
  fe->getOpPtrVector().push_back(
      new OpCalculateTensor2SymmetricFieldValuesDot<3>(
          streachTensor, dataAtPts->getLogStreachDotTensorAtPts(), 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->getSmallWDotDotAtPts(), MBTET));
  }
 
  // calculate other derived quantities
  fe->getOpPtrVector().push_back(
      new OpCalculateRotationAndSpatialGradient(dataAtPts));
 
  // evaluate integration points
  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));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialPhysical(streachTensor, 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, "w", {time_scale}, "BODY_FORCE",
        Sev::inform);
  }
 
  // Left hand side
  CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
 
  // elastic
  if (add_elastic) {
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
        streachTensor, streachTensor, dataAtPts, alphaU));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
        streachTensor, bubbleField, dataAtPts, true));
 
    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, false));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
        bubbleField, rotAxis, dataAtPts, false));
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
        streachTensor, piolaStress, dataAtPts, true));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
        streachTensor, rotAxis, 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));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_du(
        rotAxis, streachTensor, dataAtPts, false));
 
    // Stabilisation
    if constexpr (A == AssemblyType::SCHUR) {
      // Note that we assemble to AMat, however Assembly<SCHUR> assemble by
      // default to PMat
      using OpMassStab = FormsIntegrators<VolUserDataOperator>::Assembly<
          SCHUR>::BiLinearForm<GAUSS>::OpMass<1, SPACE_DIM>;
      fe_lhs->getOpPtrVector().push_back(
          new OpMassStab(rotAxis, rotAxis, [this](double, double, double) {
            return precEpsOmega;
          }));
      if (std::abs(alphaRho + alphaW) <
          std::numeric_limits<double>::epsilon()) {
 
        fe_lhs->getOpPtrVector().push_back(new OpMassStab(
            spatialL2Disp, spatialL2Disp,
            [this](double, double, double) { return 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));
    fe_rhs->getOpPtrVector().push_back(
        new OpRotationBc(piolaStress, dataAtPts, bcSpatialRotationVecPtr));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setContactElementOps(
    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
  MoFEMFunctionBegin;
 
  /* 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_rhs = boost::make_shared<BoundaryEle>(mField);
  fe_lhs = boost::make_shared<BoundaryEle>(mField);
  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});
 
  auto adj_cache =
      boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
 
  auto get_op_side = [&](auto disp_ptr, auto col_ent_data_ptr) {
    using SideEle = VolumeElementForcesAndSourcesCoreOnSide;
    auto op_loop_side = new OpLoopSide<SideEle>(
        mField, elementVolumeName, SPACE_DIM, Sev::noisy, adj_cache);
    op_loop_side->getSideFEPtr()->getUserPolynomialBase() =
        boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        op_loop_side->getOpPtrVector(), {HDIV, H1, L2});
    op_loop_side->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp, disp_ptr));
    if (col_ent_data_ptr != nullptr) {
      op_loop_side->getOpPtrVector().push_back(
          new EshelbianPlasticity::OpLoopSideGetDataForSideEle(
              spatialL2Disp, col_ent_data_ptr));
    }
    return op_loop_side;
  };
 
  auto add_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    auto col_ent_data_ptr = boost::make_shared<EntitiesFieldData::EntData>();
    auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
    pip.push_back(
        get_op_side(common_data_ptr->contactDispPtr(), col_ent_data_ptr));
    pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
        piolaStress, common_data_ptr->contactTractionPtr()));
    pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryLhs_dU(
        piolaStress, common_data_ptr, col_ent_data_ptr));
    pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryLhs_dP(
        piolaStress, common_data_ptr));
    MoFEMFunctionReturn(0);
  };
 
  auto add_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
    auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
    pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
        piolaStress, common_data_ptr->contactTractionPtr()));
    pip.push_back(get_op_side(common_data_ptr->contactDispPtr(), nullptr));
    pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryRhs(
        piolaStress, common_data_ptr));
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_ops_lhs(fe_lhs->getOpPtrVector());
  CHKERR add_ops_rhs(fe_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(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);
  }
  MoFEMFunctionReturn(0);
}
 
SmartPetscObj<KSP>
    prjKsp; ///< KSP for projection spatial displacement on H1 space
SmartPetscObj<Vec> prjD; ///< Unknown vector for projection spatial displacement
SmartPetscObj<Vec> prjF; ///< RHS for projection spatial displacement on H1
                         ///< space
SmartPetscObj<DM> prjDM;
 
MoFEMErrorCode EshelbianCore::solveElastic(TS ts, Mat m, Vec f, Vec x) {
  MoFEMFunctionBegin;
 
#ifdef PYTHON_SFD
 
  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);
 
  struct Monitor : public FEMethod {
 
    using Ele = ForcesAndSourcesCore;
    using VolEle = VolumeElementForcesAndSourcesCore;
    using VolOp = VolumeElementForcesAndSourcesCore::UserDataOperator;
    using SetPtsData = FieldEvaluatorInterface::SetPtsData;
 
    EshelbianCore &eP;
    boost::shared_ptr<SetPtsData> dataFieldEval;
    boost::shared_ptr<VolEle> volPostProcEnergy;
    boost::shared_ptr<double> gEnergy;
 
    Monitor(EshelbianCore &ep)
        : eP(ep),
          dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
                            ->getData<VolEle>()),
          volPostProcEnergy(new VolEle(ep.mField)), gEnergy(new double) {
      ierr = ep.mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
          dataFieldEval, "EP");
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
 
      auto no_rule = [](int, int, int) { return -1; };
 
      auto set_element_for_field_eval = [&]() {
        boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr.lock());
        vol_ele->getRuleHook = no_rule;
        vol_ele->getUserPolynomialBase() =
            boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
        AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
            vol_ele->getOpPtrVector(), {HDIV, H1, L2});
 
        vol_ele->getOpPtrVector().push_back(
            new OpCalculateHVecTensorField<3, 3>(
                eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
        vol_ele->getOpPtrVector().push_back(
            new OpCalculateHTensorTensorField<3, 3>(
                eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
        vol_ele->getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<3>(
                eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
                MBTET));
        vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
        // vol_ele->getOpPtrVector().push_back(
        //     new OpCalculateTensor2FieldValues<3, 3>(
        //         eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
        vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            eP.spatialL2Disp, eP.dataAtPts->getSmallWAtPts(), MBTET));
        vol_ele->getOpPtrVector().push_back(
            new OpCalculateRotationAndSpatialGradient(eP.dataAtPts));
      };
 
      auto set_element_for_post_process = [&]() {
        volPostProcEnergy->getRuleHook = VolRule();
        // Right hand side
        // CHKERR eP.setBaseVolumeElementOps(tag, true, false,
        // volPostProcEnergy);
 
        volPostProcEnergy->getUserPolynomialBase() =
            boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
        AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
            volPostProcEnergy->getOpPtrVector(), {HDIV, H1, L2});
 
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateHVecTensorField<3, 3>(
                eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateHTensorTensorField<3, 3>(
                eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<3>(
                eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
                MBTET));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValues<3>(
                eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
        // volPostProcEnergy->getOpPtrVector().push_back(
        //   new OpCalculateTensor2FieldValues<3, 3>(
        //         eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValues<3>(
                eP.spatialL2Disp, eP.dataAtPts->getSmallWAtPts(), MBTET));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateRotationAndSpatialGradient(eP.dataAtPts));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateStrainEnergy(eP.spatialL2Disp, eP.dataAtPts,
                                        gEnergy));
      };
 
      set_element_for_field_eval();
      set_element_for_post_process();
    }
 
    MoFEMErrorCode preProcess() { return 0; }
 
    MoFEMErrorCode operator()() { return 0; }
 
    MoFEMErrorCode postProcess() {
      MoFEMFunctionBegin;
 
      // auto get_str_time = [](auto ts_t) {
      //   std::ostringstream ss;
      //   ss << boost::str(boost::format("%d") %
      //                    static_cast<int>(std::ceil(ts_t * 1e6)));
      //   std::string s = ss.str();
      //   return s;
      // };
 
      auto get_step = [](auto ts_step) {
        std::ostringstream ss;
        ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
        std::string s = ss.str();
        return s;
      };
 
      PetscViewer viewer;
      CHKERR PetscViewerBinaryOpen(
          PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
          FILE_MODE_WRITE, &viewer);
      CHKERR VecView(ts_u, viewer);
      CHKERR PetscViewerDestroy(&viewer);
 
      CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_step(ts_step) +
                                          ".h5m");
 
      // Loop boundary elements with traction boundary conditions
      *gEnergy = 0;
      CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
                                            *volPostProcEnergy);
 
      double body_energy = 0;
      MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
                    eP.mField.get_comm());
      MOFEM_LOG_C("EP", Sev::inform, "Step %d time %3.4g strain energy %3.6e",
                  ts_step, ts_t, body_energy);
 
      auto post_proc_at_points = [&](std::array<double, 3> point,
                                     std::string str) {
        MoFEMFunctionBegin;
 
        dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
 
        struct OpPrint : public VolOp {
 
          EshelbianCore &eP;
          std::array<double, 3> point;
          std::string str;
 
          OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
                  std::string &str)
              : VolOp(ep.spatialL2Disp, VolOp::OPROW), eP(ep), point(point),
                str(str) {}
 
          MoFEMErrorCode doWork(int side, EntityType type,
                                EntitiesFieldData::EntData &data) {
            MoFEMFunctionBegin;
            if (type == MBTET) {
              if (getGaussPts().size2()) {
 
                auto t_h = getFTensor2FromMat<3, 3>(eP.dataAtPts->hAtPts);
                auto t_approx_P =
                    getFTensor2FromMat<3, 3>(eP.dataAtPts->approxPAtPts);
 
                FTensor::Index<'i', 3> i;
                FTensor::Index<'j', 3> j;
                FTensor::Index<'k', 3> k;
                const double jac = determinantTensor3by3(t_h);
                FTensor::Tensor2<double, 3, 3> t_cauchy;
                t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
 
                auto add = [&]() {
                  std::ostringstream s;
                  s << str << " elem " << getFEEntityHandle() << " ";
                  return s.str();
                };
 
                auto print_tensor = [](auto &t) {
                  std::ostringstream s;
                  s << t;
                  return s.str();
                };
 
                std::ostringstream print;
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "comm rank " << eP.mField.get_comm_rank();
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "point " << getVectorAdaptor(point.data(), 3);
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "coords at gauss pts " << getCoordsAtGaussPts();
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "w " << eP.dataAtPts->wAtPts;
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "Piola " << eP.dataAtPts->approxPAtPts;
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "Cauchy " << print_tensor(t_cauchy);
              }
            }
            MoFEMFunctionReturn(0);
          }
        };
 
        if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
 
          fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
          CHKERR eP.mField.getInterface<FieldEvaluatorInterface>()
              ->evalFEAtThePoint3D(
                  point.data(), 1e-12, problemPtr->getName(), "EP",
                  dataFieldEval, eP.mField.get_comm_rank(),
                  eP.mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
          fe_ptr->getOpPtrVector().pop_back();
        }
 
        MoFEMFunctionReturn(0);
      };
 
      // Points for Cook beam
      std::array<double, 3> pointA = {48., 60., 0.};
      CHKERR post_proc_at_points(pointA, "Point A");
      MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
 
      std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
      CHKERR post_proc_at_points(pointB, "Point B");
      MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
 
      std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
      CHKERR post_proc_at_points(pointC, "Point C");
      MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
 
      MoFEMFunctionReturn(0);
    }
  };
 
  boost::shared_ptr<FEMethod> monitor_ptr(new Monitor(*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::SCHUR) {
 
    auto p = matDuplicate(m, MAT_DO_NOT_COPY_VALUES);
 
    // If density is larger than zero, use dynamic time solver
    if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
      CHKERR TSSetI2Function(ts, f, PETSC_NULL, PETSC_NULL);
      CHKERR TSSetI2Jacobian(ts, m, p, PETSC_NULL, PETSC_NULL);
    } else {
      CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
      CHKERR TSSetIJacobian(ts, m, p, PETSC_NULL, PETSC_NULL);
    }
 
    KSP ksp;
    CHKERR SNESGetKSP(snes, &ksp);
    PC pc;
    CHKERR KSPGetPC(ksp, &pc);
    schur_ptr = SetUpSchur::createSetUpSchur(
        mField, SmartPetscObj<Mat>(m, true), p, &*this);
    CHKERR schur_ptr->setUp(ksp);
 
    elasticFeLhs->preProcessHook = [&]() {
      MoFEMFunctionBegin;
      *(elasticFeLhs->matAssembleSwitch) = false;
      CHKERR schur_ptr->preProc();
      MoFEMFunctionReturn(0);
    };
    elasticFeLhs->postProcessHook = [&]() {
      MoFEMFunctionBegin;
      MoFEMFunctionReturn(0);
    };
    elasticBcLhs->preProcessHook = [&]() {
      MoFEMFunctionBegin;
      MoFEMFunctionReturn(0);
    };
    elasticBcLhs->preProcessHook = [&]() {
      MoFEMFunctionBegin;
      *(elasticBcLhs->matAssembleSwitch) = false;
      CHKERR schur_ptr->postProc();
      MoFEMFunctionReturn(0);
    };
  }
 
  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);
  }
 
  auto create_post_step_ksp = [this]() {
    auto ksp = createKSP(mField.get_comm());
 
    auto set_up = [&]() {
      MoFEMFunctionBegin;
      using DomainEle = VolumeElementForcesAndSourcesCore;
      using DomainEleOp = DomainEle::UserDataOperator;
      using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
          PETSC>::BiLinearForm<GAUSS>::OpMass<1, 3>;
      using OpRhs = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
          GAUSS>::OpBaseTimesVector<1, 3, 1>;
      auto fe_lhs = boost::make_shared<DomainEle>(mField);
      auto fe_rhs = boost::make_shared<DomainEle>(mField);
 
      fe_lhs->getUserPolynomialBase() =
          boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
      fe_rhs->getUserPolynomialBase() =
          boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
 
      fe_lhs->getOpPtrVector().push_back(
          new OpDomainMass(spatialH1Disp, spatialH1Disp));
      auto w_ptr = boost::make_shared<MatrixDouble>();
      fe_rhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>(spatialL2Disp, w_ptr));
      fe_rhs->getOpPtrVector().push_back(new OpRhs(spatialH1Disp, w_ptr));
 
      CHKERR DMMoFEMKSPSetComputeOperators(dmPrjSpatial, elementVolumeName,
                                           fe_lhs, nullptr, nullptr);
      CHKERR DMMoFEMKSPSetComputeRHS(dmPrjSpatial, elementVolumeName, fe_rhs,
                                     nullptr, nullptr);
 
      CHKERR KSPAppendOptionsPrefix(ksp, "prjspatial_");
      CHKERR KSPSetFromOptions(ksp);
      CHKERR KSPSetDM(ksp, dmPrjSpatial);
      CHKERR KSPSetUp(ksp);
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(set_up(), "set up");
 
    return ksp;
  };
 
  prjKsp = create_post_step_ksp();
  prjD = createDMVector(dmPrjSpatial);
  prjF = vectorDuplicate(prjD);
  prjDM = dmPrjSpatial;
 
  auto post_step_fun = [](TS ts) {
    MoFEMFunctionBegin;
    MOFEM_LOG("EP", Sev::inform) << "Post step";
    CHKERR VecZeroEntries(prjF);
    CHKERR KSPSolve(prjKsp, prjF, prjD);
    CHKERR VecGhostUpdateBegin(prjD, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(prjD, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(prjDM, prjD, INSERT_VALUES, SCATTER_REVERSE);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR TSSetUp(ts);
  CHKERR TSSetPostStep(ts, post_step_fun);
  CHKERR TSSolve(ts, PETSC_NULL);
 
  prjKsp.reset();
  prjD.reset();
  prjF.reset();
  prjDM.reset();
 
  // 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);
 
  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 OpCalculateHTensorTensorField<3, 3>(
        bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
    fe.getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
        streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
    // fe.getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
    //     materialGradient, dataAtPts->getBigGAtPts(), MBTET));
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        spatialL2Disp, dataAtPts->getSmallWAtPts(), MBTET));
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        spatialH1Disp, dataAtPts->getUAtPtr()));
 
    // evaluate derived quantities
    fe.getOpPtrVector().push_back(
        new OpCalculateRotationAndSpatialGradient(dataAtPts));
 
    // evaluate integration points
    fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
        tag, true, false, dataAtPts, physicalEquations));
 
    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 OpPostProcDataStructure(
      post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), spatialL2Disp,
      dataAtPts));
 
  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));
 
  CHKERR mField.loop_finite_elements("ELASTIC_PROBLEM", elementVolumeName,
                                     *post_proc_norm_fe);
 
  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]
 
struct SetUpSchurImpl : public EshelbianCore::SetUpSchur {
 
  SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj<Mat> m,
                 SmartPetscObj<Mat> p, EshelbianCore *ep_core_ptr)
      : SetUpSchur(), mField(m_field), M(m), P(p), epCorePtr(ep_core_ptr) {
    if (S) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Is expected that schur matrix is not allocated. This is "
          "possible only is PC is set up twice");
    }
  }
  virtual ~SetUpSchurImpl() { S.reset(); }
 
  MoFEMErrorCode setUp(KSP solver);
  MoFEMErrorCode preProc();
  MoFEMErrorCode postProc();
 
private:
  SmartPetscObj<Mat> M;
  SmartPetscObj<Mat> P;
  SmartPetscObj<Mat> S;
  SmartPetscObj<AO> aoUp;
  SmartPetscObj<IS> isA00;
  SmartPetscObj<IS> isShift;
 
  MoFEM::Interface &mField;
  EshelbianCore *epCorePtr;
};
 
MoFEMErrorCode SetUpSchurImpl::setUp(KSP solver) {
  MoFEMFunctionBegin;
 
  auto get_ents_by_dim = [&](int dim) {
    Range ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents);
    return ents;
  };
 
  std::vector<std::string> schur_field_list{
      epCorePtr->streachTensor, epCorePtr->bubbleField, epCorePtr->rotAxis,
      epCorePtr->spatialL2Disp, epCorePtr->piolaStress};
  std::vector<std::string> field_list{epCorePtr->piolaStress};
 
  auto create_schur_dm = [&](SmartPetscObj<DM> &dm_sub) {
    MoFEMFunctionBegin;
 
    dm_sub = createDM(mField.get_comm(), "DMMOFEM");
    CHKERR DMMoFEMCreateSubDM(dm_sub, epCorePtr->dmElastic, "SUB_SCHUR");
    CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
    CHKERR DMMoFEMSetIsPartitioned(dm_sub, PETSC_TRUE);
    CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->elementVolumeName);
    CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->naturalBcElement);
 
    auto faces_ptr = boost::make_shared<Range>(get_ents_by_dim(SPACE_DIM - 1));
    std::vector<boost::shared_ptr<Range>> dm_range_list{faces_ptr};
 
    if (field_list.size() != dm_range_list.size()) {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
              "Both ranges should have the same size");
    }
 
    int r_idx = 0;
    for (auto f : field_list) {
      MOFEM_LOG("EP", Sev::inform) << "Add schur field: " << f;
      CHKERR DMMoFEMAddSubFieldRow(dm_sub, f.c_str(), dm_range_list[r_idx]);
      CHKERR DMMoFEMAddSubFieldCol(dm_sub, f.c_str(), dm_range_list[r_idx]);
      ++r_idx;
    }
    CHKERR DMSetUp(dm_sub);
    MoFEMFunctionReturn(0);
  };
 
  auto create_a00_is = [&](SmartPetscObj<IS> &is_a00) {
    MoFEMFunctionBegin;
    auto vols = get_ents_by_dim(SPACE_DIM);
    std::vector<SmartPetscObj<IS>> is_vec;
    std::vector<Range *> range_list_ptr(schur_field_list.size(), nullptr);
    range_list_ptr.back() = &vols;
    int r_idx = 0;
    for (auto f : schur_field_list) {
      SmartPetscObj<IS> is;
      CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
          "ELASTIC_PROBLEM", ROW, f, 0, MAX_DOFS_ON_ENTITY, is,
          range_list_ptr[r_idx]);
      is_vec.push_back(is);
      ++r_idx;
    }
    if (is_vec.size()) {
      is_a00 = is_vec[0];
      for (auto i = 1; i < is_vec.size(); ++i) {
        IS is_union_raw;
        CHKERR ISExpand(is_a00, is_vec[i], &is_union_raw);
        is_a00 = SmartPetscObj<IS>(is_union_raw);
      }
      CHKERR ISSort(is_a00);
    }
    MoFEMFunctionReturn(0);
  };
 
  PC pc;
  CHKERR KSPSetFromOptions(solver);
  CHKERR KSPGetPC(solver, &pc);
  PetscBool is_pcfs = PETSC_FALSE;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
  if (is_pcfs) {
 
    SmartPetscObj<DM> schur_dm, a00_dm;
    CHKERR create_schur_dm(schur_dm);
    CHKERR create_a00_is(isA00);
 
    if (S) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Is expected that schur matrix is not allocated. This is "
          "possible only is PC is set up twice");
    }
 
    S = createDMMatrix(schur_dm);
 
    auto set_ops = [&]() {
      MoFEMFunctionBegin;
      std::vector<boost::shared_ptr<Range>> ranges_list(schur_field_list.size(),
                                                        nullptr);
      ranges_list.back() =
          boost::make_shared<Range>(get_ents_by_dim(SPACE_DIM));
      std::vector<SmartPetscObj<AO>> ao_list(schur_field_list.size(),
                                             SmartPetscObj<AO>());
      std::vector<SmartPetscObj<Mat>> mat_list(schur_field_list.size(),
                                               SmartPetscObj<Mat>());
      std::vector<bool> symm_list(schur_field_list.size(), false);
 
      auto dm_is = getDMSubData(schur_dm)->getSmartRowIs();
      aoUp = createAOMappingIS(dm_is, PETSC_NULL);
      ao_list.back() = aoUp;
      mat_list.back() = S;
 
      epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
          new OpSchurAssembleBegin());
      epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
          new OpSchurAssembleEnd<SCHUR_DGESV>(schur_field_list, ranges_list,
                                              ao_list, mat_list, symm_list,
                                              false));
      epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
          new OpSchurAssembleBegin());
      epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
          new OpSchurAssembleEnd<SCHUR_DGESV>(schur_field_list, ranges_list,
                                              ao_list, mat_list, symm_list,
                                              false));
      MoFEMFunctionReturn(0);
    };
 
    auto set_pc = [&]() {
      MoFEMFunctionBegin;
      CHKERR PCFieldSplitSetIS(pc, NULL, isA00);
      CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_ops();
    CHKERR set_pc();
 
  } else {
    epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
        new OpSchurAssembleBegin());
    epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
        new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
    epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
        new OpSchurAssembleBegin());
    epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
        new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode SetUpSchurImpl::preProc() {
  MoFEMFunctionBegin;
  MOFEM_LOG("EP", Sev::verbose) << "Zero Schur";
  if (SetUpSchurImpl::S) {
    CHKERR MatZeroEntries(S);
  }
  CHKERR MatZeroEntries(M);
  CHKERR MatZeroEntries(P);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode SetUpSchurImpl::postProc() {
  MoFEMFunctionBegin;
  if (S) {
    MOFEM_LOG("EP", Sev::verbose) << "Assemble Schur";
    CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
  }
  CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
  CHKERR MatAXPY(P, 1, M, SAME_NONZERO_PATTERN);
 
  MoFEMFunctionReturn(0);
}
 
boost::shared_ptr<EshelbianCore::SetUpSchur>
EshelbianCore::SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
                                            SmartPetscObj<Mat> m,
                                            SmartPetscObj<Mat> p,
                                            EshelbianCore *ep_core_ptr) {
  return boost::shared_ptr<SetUpSchur>(
      new SetUpSchurImpl(m_field, m, p, ep_core_ptr));
}
 
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>