EshelbianPlasticity.cpp

Eshelbian plasticity implementation

/**
 * \file EshelbianPlasticity.cpp
 * \example EshelbianPlasticity.cpp
 *
 * \brief Eshelbian plasticity implementation
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <BasicFiniteElements.hpp>
#include <CGGTonsorialBubbleBase.hpp>
 
#include <EshelbianPlasticity.hpp>
#include <boost/math/constants/constants.hpp>
 
namespace EshelbianPlasticity {
 
MoFEMErrorCode EshelbianCore::query_interface(const MOFEMuuid &uuid,
                                              UnknownInterface **iface) const {
  MoFEMFunctionBeginHot;
  *iface = NULL;
  if (uuid == IDD_MOFEMEshelbianCrackInterface) {
    *iface = const_cast<EshelbianCore *>(this);
    MoFEMFunctionReturnHot(0);
  }
  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
  if (type != MBVERTEX)
    MoFEMFunctionReturnHot(0);
 
  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"), spatialDisp("w"),
      materialDisp("W"), streachTensor("u"), rotAxis("omega"),
      materialGradient("G"), tauField("TAU"), lambdaField("LAMBDA"),
      bubbleField("BUBBLE"), elementVolumeName("EP"),
      naturalBcElement("NATURAL_BC"), essentialBcElement("ESSENTIAL_BC") {
 
  ierr = getOptions();
  CHKERRABORT(PETSC_COMM_WORLD, ierr);
}
 
EshelbianCore::~EshelbianCore() {}
 
MoFEMErrorCode EshelbianCore::getOptions() {
  MoFEMFunctionBegin;
  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 = 1e-3;
  CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
                            precEps, &precEps, PETSC_NULL);
 
  ierr = PetscOptionsEnd();
  CHKERRG(ierr);
  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);
 
  auto subtract_faces_where_displacements_are_applied =
      [&](const std::string disp_block_set_name) {
        MoFEMFunctionBegin;
        for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
          if (it->getName().compare(0, disp_block_set_name.length(),
                                    disp_block_set_name) == 0) {
            Range faces;
            CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
                                                       faces, true);
            tets_skin = subtract(tets_skin, faces);
          }
        }
        MoFEMFunctionReturn(0);
      };
  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_DISP_BC");
  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_ROTATION_BC");
  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_TRACTION_BC");
 
  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_hdiv_rt_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_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, 0);
    CHKERR mField.set_field_order(tets_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_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(spatialDisp, 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(materialDisp, 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);
 
  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, spatialDisp);
    CHKERR add_field_to_fe(elementVolumeName, streachTensor);
    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 disp_block_set_name,
                                      Range &r) {
    MoFEMFunctionBegin;
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, disp_block_set_name.length(),
                                disp_block_set_name) == 0) {
        Range faces;
        CHKERR it->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;
  CHKERR bc_elements_add_to_range("SPATIAL_DISP_BC", natural_bc_elements);
  CHKERR bc_elements_add_to_range("SPATIAL_ROTATION_BC", natural_bc_elements);
  Range essentail_bc_elements;
  CHKERR bc_elements_add_to_range("SPATIAL_TRACTION_BC", essentail_bc_elements);
 
  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);
 
  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 = createSmartDM(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.c_str());
  CHKERR DMMoFEMAddElement(dM, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dM, essentialBcElement.c_str());
  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, NOISY,
              true);
        MoFEMFunctionReturn(0);
      };
  CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
 
  // Create elastic sub-problem
  dmElastic = createSmartDM(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, spatialDisp.c_str());
  CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName.c_str());
  CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dmElastic, essentialBcElement.c_str());
  CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
  CHKERR DMSetUp(dmElastic);
 
  // Create elastic streach-problem
  dmElasticSchurStreach = createSmartDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElasticSchurStreach, dmElastic,
                            "ELASTIC_PROBLEM_STREACH_SCHUR");
  CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurStreach, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, piolaStress.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, bubbleField.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, rotAxis.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, spatialDisp.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurStreach, elementVolumeName.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurStreach, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurStreach, essentialBcElement.c_str());
  CHKERR DMMoFEMSetSquareProblem(dmElasticSchurStreach, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurStreach, PETSC_TRUE);
  CHKERR DMSetUp(dmElasticSchurStreach);
 
  // Create elastic bubble-problem
  dmElasticSchurBubble = createSmartDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElasticSchurBubble, dmElasticSchurStreach,
                            "ELASTIC_PROBLEM_BUBBLE_SCHUR");
  CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurBubble, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, piolaStress.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, rotAxis.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, spatialDisp.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurBubble, elementVolumeName.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurBubble, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurBubble, essentialBcElement.c_str());
  CHKERR DMMoFEMSetSquareProblem(dmElasticSchurBubble, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurBubble, PETSC_TRUE);
  CHKERR DMSetUp(dmElasticSchurBubble);
 
  // Create elastic omega-problem
  dmElasticSchurOmega = createSmartDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElasticSchurOmega, dmElasticSchurBubble,
                            "ELASTIC_PROBLEM_OMEGA_SCHUR");
  CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurOmega, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurOmega, piolaStress.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurOmega, spatialDisp.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurOmega, elementVolumeName.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurOmega, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurOmega, essentialBcElement.c_str());
  CHKERR DMMoFEMSetSquareProblem(dmElasticSchurOmega, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurOmega, PETSC_TRUE);
  CHKERR DMSetUp(dmElasticSchurOmega);
 
  // Create elastic tet_stress-problem
  dmElasticSchurSpatialDisp = createSmartDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElasticSchurSpatialDisp, dmElasticSchurOmega,
                            "ELASTIC_PROBLEM_SPATIAL_DISP_SCHUR");
  CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurSpatialDisp, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurSpatialDisp, piolaStress.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp,
                           elementVolumeName.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp,
                           essentialBcElement.c_str());
  CHKERR DMMoFEMSetSquareProblem(dmElasticSchurSpatialDisp, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurSpatialDisp, PETSC_TRUE);
  CHKERR DMSetUp(dmElasticSchurSpatialDisp);
 
  {
    PetscSection section;
    CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
                                                           &section);
    CHKERR DMSetDefaultSection(dmElastic, section);
    CHKERR DMSetDefaultGlobalSection(dmElastic, section);
    CHKERR PetscSectionDestroy(&section);
  }
 
  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]);
  }
}
 
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]);
  }
}
 
MoFEMErrorCode
EshelbianCore::getTractionFreeBc(const EntityHandle meshset,
                                 boost::shared_ptr<TractionFreeBc> &bc_ptr,
                                 const std::string disp_block_set_name,
                                 const std::string rot_block_set_name) {
  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);
 
  bc_ptr->resize(3);
  for (int dd = 0; dd != 3; ++dd)
    (*bc_ptr)[dd] = tets_skin;
 
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
    if (it->getName().compare(0, disp_block_set_name.length(),
                              disp_block_set_name) == 0) {
      std::vector<double> block_attributes;
      CHKERR it->getAttributes(block_attributes);
      if (block_attributes.size() != 6) {
        SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "In block %s six attributes are required for given BC "
                 "blockset (3 values + "
                 "3 flags) != %d",
                 it->getName().c_str(), block_attributes.size());
      }
      Range faces;
      CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                                 true);
      if (block_attributes[3] != 0)
        (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
      if (block_attributes[4] != 0)
        (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
      if (block_attributes[5] != 0)
        (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
    }
    if (it->getName().compare(0, rot_block_set_name.length(),
                              rot_block_set_name) == 0) {
      Range faces;
      CHKERR it->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 2 * (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(const MOFEMuuid &uuid,
                                 BaseFunctionUnknownInterface **iface) const {
    MoFEMFunctionBegin;
    *iface = NULL;
    if (uuid == IDD_TET_BASE_FUNCTION) {
      *iface = const_cast<CGGUserPolynomialBase *>(this);
      MoFEMFunctionReturnHot(0);
    } else {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
    }
    CHKERR BaseFunction::query_interface(uuid, iface);
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode getValue(MatrixDouble &pts,
                          boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
    MoFEMFunctionBeginHot;
 
    BaseFunctionUnknownInterface *iface;
    CHKERR ctx_ptr->query_interface(IDD_TET_BASE_FUNCTION, &iface);
    cTx = reinterpret_cast<EntPolynomialBaseCtx *>(iface);
 
    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
    DataForcesAndSourcesCore &data = cTx->dAta;
    // Get approximation order on element. Note that bubble functions are only
    // on tetrahedron.
    const int order = data.dataOnEntities[MBTET][0].getDataOrder();
    /// 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<EpElement<VolumeElementForcesAndSourcesCore>> &fe) {
  MoFEMFunctionBegin;
  fe = boost::make_shared<EpElement<VolumeElementForcesAndSourcesCore>>(mField);
 
  fe->getUserPolynomialBase() =
      boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
  fe->getOpPtrVector().push_back(new OpL2Transform());
 
  // 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>(
      spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
 
  // velocities
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
      spatialDisp, 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>(
        spatialDisp, dataAtPts->getSmallWDotDotAtPts(), MBTET));
  }
 
  // calculate other derived quantities
  fe->getOpPtrVector().push_back(
      new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
 
  // evaluate integration points
  fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
      spatialDisp, tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setGenericVolumeElementOps(
    const int tag, const bool add_elastic, const bool add_material,
    boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_rhs,
    boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_lhs) {
  MoFEMFunctionBegin;
 
  // Right hand side
  CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
 
  // elastic
  if (add_elastic) {
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialEquilibrium(spatialDisp, 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));
  }
 
  // Left hand side
  CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
 
  // elastic
  if (add_elastic) {
 
    // Schur
    fe_lhs->getOpPtrVector().push_back(
        new OpSpatialSchurBegin(spatialDisp, dataAtPts));
 
    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(
        spatialDisp, piolaStress, dataAtPts, true));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
        spatialDisp, spatialDisp, dataAtPts, alphaW, alphaRho));
 
    fe_lhs->getOpPtrVector().push_back(
        new OpSpatialConsistency_dP_domega(piolaStress, rotAxis, dataAtPts));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
        bubbleField, rotAxis, dataAtPts));
 
    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));
 
    // Schur
    dataAtPts->ooMatPtr = boost::make_shared<MatrixDouble>();
    fe_lhs->getOpPtrVector().push_back(
        new OpSpatialPreconditionMass(rotAxis, dataAtPts->ooMatPtr));
    if (alphaW < std::numeric_limits<double>::epsilon() &&
        alphaRho < std::numeric_limits<double>::epsilon()) {
      dataAtPts->wwMatPtr = boost::make_shared<MatrixDouble>();
      fe_lhs->getOpPtrVector().push_back(
          new OpSpatialPreconditionMass(spatialDisp, dataAtPts->wwMatPtr));
    } else {
      dataAtPts->wwMatPtr.reset();
    }
    fe_lhs->getOpPtrVector().push_back(
        new OpSpatialSchurEnd(spatialDisp, dataAtPts, precEps));
 
    if (add_material) {
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setGenericFaceElementOps(
    const bool add_elastic, const bool add_material,
    boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_rhs,
    boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_lhs) {
  MoFEMFunctionBegin;
 
  fe_rhs =
      boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
  fe_lhs =
      boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
 
  // set integration rule
  fe_rhs->getRuleHook = FaceRule();
  fe_lhs->getRuleHook = FaceRule();
 
  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::setElasticElementOps(const int tag) {
  MoFEMFunctionBegin;
  CHKERR setGenericVolumeElementOps(tag, true, false, elasticFeRhs,
                                    elasticFeLhs);
  CHKERR setGenericFaceElementOps(true, false, elasticBcRhs, elasticBcLhs);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setElasticElementToTs(DM dm) {
  MoFEMFunctionBegin;
  boost::shared_ptr<FEMethod> null;
  boost::shared_ptr<EpElement<FeTractionBc>> spatial_traction_bc(
      new EpElement<FeTractionBc>(mField, piolaStress, bcSpatialTraction));
  schurAssembly = boost::make_shared<EpFEMethod>();
 
  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, DM_NO_ELEMENT, null, null,
                                  spatial_traction_bc);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, DM_NO_ELEMENT, null, schurAssembly, null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, DM_NO_ELEMENT, null, null, schurAssembly);
  } 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, DM_NO_ELEMENT, null, null,
                                 spatial_traction_bc);
    CHKERR DMMoFEMTSSetIJacobian(dm, DM_NO_ELEMENT, null, schurAssembly, null);
    CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIJacobian(dm, DM_NO_ELEMENT, null, null, schurAssembly);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::setUpTSElastic(TS ts, Mat m, Vec f) {
  MoFEMFunctionBegin;
  boost::shared_ptr<TsCtx> ts_ctx;
  CHKERR DMMoFEMGetTsCtx(dmElastic, ts_ctx);
 
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
    CHKERR TSSetI2Function(ts, f, PETSC_NULL, PETSC_NULL);
    CHKERR TSSetI2Jacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
  } else {
    CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
    CHKERR TSSetIJacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
  }
 
  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
 
  auto add_schur_streach_op = [this](auto &list, SmartPetscObj<Mat> S,
                                     SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
        fe_cast->addStreachSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_streach_pre = [this](auto &list, SmartPetscObj<Mat> S,
                                      SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
        fe_cast->addStreachSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_bubble_op = [this](auto &list, SmartPetscObj<Mat> S,
                                    SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
        fe_cast->addBubbleSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_bubble_pre = [this](auto &list, SmartPetscObj<Mat> S,
                                     SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
        fe_cast->addBubbleSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_spatial_disp_op = [this](auto &list, SmartPetscObj<Mat> S,
                                          SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
        fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_spatial_disp_pre = [this](auto &list, SmartPetscObj<Mat> S,
                                           SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
        fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_omega_op = [this](auto &list, SmartPetscObj<Mat> S,
                                   SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
        fe_cast->addOmegaSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_omega_pre = [this](auto &list, SmartPetscObj<Mat> S,
                                    SmartPetscObj<AO> ao) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
        fe_cast->addOmegaSchurMatrix(S, ao);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  const MoFEM::Problem *schur_streach_prb_ptr;
  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_streach_prb_ptr);
  if (auto sub_data = schur_streach_prb_ptr->subProblemData) {
    SmartPetscObj<Mat> Suu;
    CHKERR DMCreateMatrix_MoFEM(dmElasticSchurStreach, Suu);
    SmartPetscObj<AO> aoSuu = sub_data->getSmartRowMap();
 
    CHKERR add_schur_streach_op(ts_ctx->loops_to_do_IJacobian, Suu, aoSuu);
    CHKERR add_schur_streach_pre(ts_ctx->preProcess_IJacobian, Suu, aoSuu);
    CHKERR add_schur_streach_pre(ts_ctx->postProcess_IJacobian, Suu, aoSuu);
 
    const MoFEM::Problem *schur_bubble_prb_ptr;
    CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_prb_ptr);
    if (auto bubble_data = schur_bubble_prb_ptr->subProblemData) {
      SmartPetscObj<Mat> SBubble;
      CHKERR DMCreateMatrix_MoFEM(dmElasticSchurBubble, SBubble);
      SmartPetscObj<AO> aoSBubble = bubble_data->getSmartRowMap();
      CHKERR add_schur_bubble_op(ts_ctx->loops_to_do_IJacobian, SBubble,
                                 aoSBubble);
      CHKERR add_schur_bubble_pre(ts_ctx->preProcess_IJacobian, SBubble,
                                  aoSBubble);
      CHKERR add_schur_bubble_pre(ts_ctx->postProcess_IJacobian, SBubble,
                                  aoSBubble);
 
      const MoFEM::Problem *schur_omega_prb_ptr;
      CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_prb_ptr);
      if (auto tet_stress_data = schur_omega_prb_ptr->subProblemData) {
        SmartPetscObj<Mat> SOmega;
        CHKERR DMCreateMatrix_MoFEM(dmElasticSchurOmega, SOmega);
        SmartPetscObj<AO> aoSOmega = tet_stress_data->getSmartRowMap();
        CHKERR add_schur_omega_op(ts_ctx->loops_to_do_IJacobian, SOmega,
                                  aoSOmega);
        CHKERR add_schur_omega_pre(ts_ctx->preProcess_IJacobian, SOmega,
                                   aoSOmega);
        CHKERR add_schur_omega_pre(ts_ctx->postProcess_IJacobian, SOmega,
                                   aoSOmega);
 
        const MoFEM::Problem *schur_spatial_disp_prb_ptr;
        CHKERR DMMoFEMGetProblemPtr(dmElasticSchurSpatialDisp,
                                    &schur_spatial_disp_prb_ptr);
        if (auto spatial_disp_data =
                schur_spatial_disp_prb_ptr->subProblemData) {
          SmartPetscObj<Mat> Sw;
          CHKERR DMCreateMatrix_MoFEM(dmElasticSchurSpatialDisp, Sw);
          SmartPetscObj<AO> aoSw = spatial_disp_data->getSmartRowMap();
          CHKERR add_schur_spatial_disp_op(ts_ctx->loops_to_do_IJacobian, Sw,
                                           aoSw);
          CHKERR add_schur_spatial_disp_pre(ts_ctx->preProcess_IJacobian, Sw,
                                            aoSw);
          CHKERR add_schur_spatial_disp_pre(ts_ctx->postProcess_IJacobian, Sw,
                                            aoSw);
        } else
          SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                  "Problem does not have sub-problem data");
 
      } else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "Problem does not have sub-problem data");
 
    } else
      SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Problem does not have sub-problem data");
 
  } else
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Problem does not have sub-problem data");
 
  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());
        vol_ele->getOpPtrVector().push_back(new OpL2Transform());
 
        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.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
        vol_ele->getOpPtrVector().push_back(
            new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
                                                      eP.dataAtPts));
      };
 
      auto set_element_for_post_process = [&]() {
        volPostProcEnergy->getRuleHook = VolRule();
        volPostProcEnergy->getUserPolynomialBase() =
            boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
        volPostProcEnergy->getOpPtrVector().push_back(new OpL2Transform());
 
        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.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
                                                      eP.dataAtPts));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateStrainEnergy(eP.spatialDisp, 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;
      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.spatialDisp, VolOp::OPROW), eP(ep), point(point),
                str(str) {}
 
          MoFEMErrorCode doWork(int side, EntityType type,
                                DataForcesAndSourcesCore::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 = dEterminant(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(), MF_EXIST, QUIET);
          fe_ptr->getOpPtrVector().pop_back();
        }
 
        MoFEMFunctionReturn(0);
      };
 
      // Points for Cook beam
      std::array<double, 3> pointA = {48., 60., 5.};
      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->get_loops_to_do_Monitor().push_back(
      TsCtx::PairNameFEMethodPtr(elementVolumeName, monitor_ptr));
 
  CHKERR TSAppendOptionsPrefix(ts, "elastic_");
  CHKERR TSSetFromOptions(ts);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::solveElastic(TS ts, Vec x) {
  MoFEMFunctionBegin;
 
  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 DMGetDefaultSection(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
 
  KSP ksp;
  CHKERR SNESGetKSP(snes, &ksp);
  PC pc;
  CHKERR KSPGetPC(ksp, &pc);
  PetscBool is_uu_field_split;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_uu_field_split);
  if (is_uu_field_split) {
 
    const MoFEM::Problem *schur_uu_ptr;
    CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_uu_ptr);
    if (auto uu_data = schur_uu_ptr->subProblemData) {
 
      const MoFEM::Problem *prb_ptr;
      CHKERR DMMoFEMGetProblemPtr(dmElastic, &prb_ptr);
      map<std::string, IS> is_map;
      for (int ff = 0; ff != num_fields; ff++) {
        const char *field_name;
        CHKERR PetscSectionGetFieldName(section, ff, &field_name);
        CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
            prb_ptr->getName(), ROW, field_name, 0, MAX_DOFS_ON_ENTITY,
            &is_map[field_name]);
      }
      // CHKERR mField.getInterface<ISManager>()
      //     ->isCreateProblemFieldAndEntityType(
      //         prb_ptr->getName(), ROW, piolaStress, MBTET, MBTET, 0,
      //         MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TETS"]);
      // CHKERR mField.getInterface<ISManager>()
      //     ->isCreateProblemFieldAndEntityType(
      //         prb_ptr->getName(), ROW, piolaStress, MBTRI, MBTRI, 0,
      //         MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TRIS"]);
 
      CHKERR uu_data->getRowIs(&is_map["E_IS_SUU"]);
 
      CHKERR PCFieldSplitSetIS(pc, NULL, is_map[streachTensor]);
      CHKERR PCFieldSplitSetIS(pc, NULL, is_map["E_IS_SUU"]);
 
      CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
                                     schurAssembly->Suu);
 
      CHKERR PCSetUp(pc);
      PetscInt n;
      KSP *uu_ksp;
      CHKERR PCFieldSplitGetSubKSP(pc, &n, &uu_ksp);
      PC bubble_pc;
      CHKERR KSPGetPC(uu_ksp[1], &bubble_pc);
      PetscBool is_bubble_field_split;
      PetscObjectTypeCompare((PetscObject)bubble_pc, PCFIELDSPLIT,
                             &is_bubble_field_split);
      if (is_bubble_field_split) {
 
        const MoFEM::Problem *schur_bubble_ptr;
        CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_ptr);
        if (auto bubble_data = schur_bubble_ptr->subProblemData) {
 
          CHKERR bubble_data->getRowIs(&is_map["E_IS_BUBBLE"]);
 
          AO uu_ao;
          CHKERR uu_data->getRowMap(&uu_ao);
 
          CHKERR AOApplicationToPetscIS(uu_ao, is_map[bubbleField]);
          CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map[bubbleField]);
          CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map["E_IS_BUBBLE"]);
          CHKERR PCFieldSplitSetSchurPre(
              bubble_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->SBubble);
 
          CHKERR PCSetUp(bubble_pc);
          PetscInt bubble_n;
          KSP *bubble_ksp;
          CHKERR PCFieldSplitGetSubKSP(bubble_pc, &bubble_n, &bubble_ksp);
          PC omega_pc;
          CHKERR KSPGetPC(bubble_ksp[1], &omega_pc);
          PetscBool is_omega_field_split;
          PetscObjectTypeCompare((PetscObject)omega_pc, PCFIELDSPLIT,
                                 &is_omega_field_split);
 
          if (is_omega_field_split) {
 
            const MoFEM::Problem *schur_omega_ptr;
            CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_ptr);
            if (auto omega_data = schur_omega_ptr->subProblemData) {
 
              AO bubble_ao;
              CHKERR bubble_data->getRowMap(&bubble_ao);
 
              CHKERR AOApplicationToPetscIS(uu_ao, is_map[rotAxis]);
              CHKERR AOApplicationToPetscIS(bubble_ao, is_map[rotAxis]);
              CHKERR omega_data->getRowIs(&is_map["E_IS_OMEGA"]);
 
              CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map[rotAxis]);
              CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map["E_IS_OMEGA"]);
              CHKERR PCFieldSplitSetSchurPre(omega_pc,
                                             PC_FIELDSPLIT_SCHUR_PRE_USER,
                                             schurAssembly->SOmega);
 
              CHKERR PCSetUp(omega_pc);
              PetscInt omega_n;
              KSP *omega_ksp;
              CHKERR PCFieldSplitGetSubKSP(omega_pc, &omega_n, &omega_ksp);
              PC w_pc;
              CHKERR KSPGetPC(omega_ksp[1], &w_pc);
              PetscBool is_w_field_split;
              PetscObjectTypeCompare((PetscObject)w_pc, PCFIELDSPLIT,
                                     &is_w_field_split);
              if (is_w_field_split) {
 
                const MoFEM::Problem *schur_w_ptr;
                CHKERR DMMoFEMGetProblemPtr(dmElasticSchurSpatialDisp,
                                            &schur_w_ptr);
                if (auto w_data = schur_w_ptr->subProblemData) {
 
                  AO omega_ao;
                  CHKERR omega_data->getRowMap(&omega_ao);
 
                  CHKERR AOApplicationToPetscIS(uu_ao, is_map[spatialDisp]);
                  CHKERR AOApplicationToPetscIS(bubble_ao, is_map[spatialDisp]);
                  CHKERR AOApplicationToPetscIS(omega_ao, is_map[spatialDisp]);
                  CHKERR w_data->getRowIs(&is_map["E_IS_W"]);
 
                  CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map[spatialDisp]);
                  CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map["E_IS_W"]);
                  CHKERR PCFieldSplitSetSchurPre(
                      w_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->Sw);
 
                  CHKERR AODestroy(&omega_ao);
                }
              }
 
              CHKERR PetscFree(omega_ksp);
            }
          }
          CHKERR PetscFree(bubble_ksp);
          CHKERR AODestroy(&uu_ao);
        }
      }
      CHKERR PetscFree(uu_ksp);
 
      for (auto &m : is_map)
        CHKERR ISDestroy(&m.second);
    }
  }
 
  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);
  }
 
  CHKERR TSSolve(ts, PETSC_NULL);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::postProcessResults(const int tag,
                                                 const std::string file) {
  MoFEMFunctionBegin;
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
  }
 
  PostProcVolumeOnRefinedMesh post_proc(mField);
  post_proc.getUserPolynomialBase() =
      boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
  post_proc.getOpPtrVector().push_back(new OpL2Transform());
 
  CHKERR post_proc.generateReferenceElementMesh();
  post_proc.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
      piolaStress, dataAtPts->getApproxPAtPts()));
  post_proc.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
      bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  post_proc.getOpPtrVector().push_back(
      new OpCalculateTensor2SymmetricFieldValues<3>(
          streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  post_proc.getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
      materialGradient, dataAtPts->getBigGAtPts(), MBTET));
  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
 
  // evaluate derived quantities
  post_proc.getOpPtrVector().push_back(
      new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
 
  // evaluate integration points
  post_proc.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
      spatialDisp, tag, true, false, dataAtPts, physicalEquations));
  post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
      post_proc.postProcMesh, post_proc.mapGaussPts, spatialDisp, dataAtPts));
 
  CHKERR DMoFEMLoopFiniteElements(dM, elementVolumeName.c_str(), &post_proc);
  CHKERR post_proc.writeFile(file.c_str());
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity