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);

  alpha_u = 0;
  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alpha_u,
                            &alpha_u, PETSC_NULL);

  alpha_w = 0;
  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alpha_w,
                            &alpha_w, PETSC_NULL);

  preconditioner_eps = 1e-3;
  CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
                            preconditioner_eps, &preconditioner_eps,
                            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 OpCalculateVectorFieldValuesDot<3>(
      rotAxis, dataAtPts->getRotAxisDotAtPts(), 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));

  // 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, alpha_w));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialRotation(rotAxis, dataAtPts));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialPhysical(streachTensor, dataAtPts, alpha_u));
    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, alpha_u));
    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, alpha_w));

    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 (alpha_w < 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, preconditioner_eps));

    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));

  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);

  schurAssembly = boost::make_shared<EpFEMethod>();
  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);
  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, Mat S, 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, Mat S, 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, Mat S, 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, Mat S, 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, Mat S, 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, Mat S, 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, Mat S, 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, Mat S, 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) {
    Mat Suu;
    CHKERR DMCreateMatrix(dmElasticSchurStreach, &Suu);
    AO aoSuu;
    CHKERR sub_data->getRowMap(&aoSuu);

    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);

    CHKERR MatDestroy(&Suu);
    CHKERR AODestroy(&aoSuu);

    const MoFEM::Problem *schur_bubble_prb_ptr;
    CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_prb_ptr);
    if (auto bubble_data = schur_bubble_prb_ptr->subProblemData) {
      Mat SBubble;
      CHKERR DMCreateMatrix(dmElasticSchurBubble, &SBubble);
      AO aoSBubble;
      CHKERR bubble_data->getRowMap(&aoSBubble);

      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);

      CHKERR MatDestroy(&SBubble);
      CHKERR AODestroy(&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) {
        Mat SOmega;
        CHKERR DMCreateMatrix(dmElasticSchurOmega, &SOmega);
        AO aoSOmega;
        CHKERR tet_stress_data->getRowMap(&aoSOmega);

        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) {

          Mat Sw;
          CHKERR DMCreateMatrix(dmElasticSchurSpatialDisp, &Sw);
          AO aoSw;
          CHKERR spatial_disp_data->getRowMap(&aoSw);

          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);

          CHKERR MatDestroy(&Sw);
          CHKERR AODestroy(&aoSw);
        } else
          SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                  "Problem does not have sub-problem data");

        CHKERR MatDestroy(&SOmega);
        CHKERR AODestroy(&aoSOmega);

      } 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> volPostProcEle;
    boost::shared_ptr<double> gEnergy;

    Monitor(EshelbianCore &ep)
        : eP(ep),
          dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
                            ->getData<VolEle>()),
          volPostProcEle(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 = [&]() {
        volPostProcEle->getRuleHook = VolRule();
        volPostProcEle->getUserPolynomialBase() =
            boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
        volPostProcEle->getOpPtrVector().push_back(new OpL2Transform());

        volPostProcEle->getOpPtrVector().push_back(
            new OpCalculateHVecTensorField<3, 3>(
                eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
        volPostProcEle->getOpPtrVector().push_back(
            new OpCalculateHTensorTensorField<3, 3>(
                eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
        volPostProcEle->getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<3>(
                eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(), MBTET));
        volPostProcEle->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValues<3>(
                eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
        volPostProcEle->getOpPtrVector().push_back(
            new OpCalculateTensor2FieldValues<3, 3>(
                eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
        volPostProcEle->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValues<3>(
                eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
        volPostProcEle->getOpPtrVector().push_back(
            new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
                                                      eP.dataAtPts));
        volPostProcEle->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;
      };

      PetscViewer viewer;
      CHKERR PetscViewerBinaryOpen(
          PETSC_COMM_WORLD, ("restart_" + get_str_time(ts_t) + ".dat").c_str(),
          FILE_MODE_WRITE, &viewer);
      CHKERR VecView(ts_u, viewer);
      CHKERR PetscViewerDestroy(&viewer);

      CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_str_time(ts_t) +
                                          ".h5m");

      // Loop boundary elements with traction boundary conditions
      *gEnergy = 0;
      CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
                                            *volPostProcEle,nullptr);

      double body_energy;
      MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
                    eP.mField.get_comm());
      CHKERR PetscPrintf(eP.mField.get_comm(),
                         "Step %d time %3.4g strain energy %3.6e\n", 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();
                };

                std::ostringstream print;
                print << add() << "comm rank " << eP.mField.get_comm_rank()
                      << std::endl;
                print << add() << "point " << getVectorAdaptor(point.data(), 3)
                      << std::endl;
                print << add() << "coords at gauss pts " << getCoordsAtGaussPts()
                      << std::endl;
                print << add() << "w " << eP.dataAtPts->wAtPts << std::endl;
                print << add() << "Piola " << eP.dataAtPts->approxPAtPts
                      << std::endl;
                print << add() << "Cauchy " << t_cauchy << std::endl;
                print << std::endl;
                CHKERR PetscSynchronizedPrintf(eP.mField.get_comm(), "%s",
                                               print.str().c_str());
                                               
              }
            }
            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");
      PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);

      std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
      CHKERR post_proc_at_points(pointB, "Point B");
      PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);

      std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
      CHKERR post_proc_at_points(pointC, "Point C");
      PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);

      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);
    CHKERR PetscPrintf(PETSC_COMM_WORLD, "Field %d name %s\n", ff, field_name);
  }

  CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);

  // PetscRandom rctx;
  // PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
  // VecSetRandom(x, rctx);
  // VecScale(x, -1e-1);
  // PetscRandomDestroy(&rctx);
  // CHKERR VecView(x, PETSC_VIEWER_STDOUT_WORLD);

  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);

  // CHKERR VecView(x, PETSC_VIEWER_STDOUT_WORLD);

  // 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);
    }
  }

  CHKERR TSSolve(ts, x);
  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