adolc_plasticity.cpp

\brife Implementation of plasticity problem using ADOLC-C

/**
 * \file adolc_plasticity.cpp
 * \ingroup adoc_plasticity
 * \example adolc_plasticity.cpp
 *
 * \brife Implementation of plasticity problem using ADOLC-C
 *
 */
#include <MoFEM.hpp>
#include <MatrixFunction.hpp>
 
using namespace MoFEM;
 
constexpr int SPACE_DIM =
    EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
 
using EntData = EntitiesFieldData::EntData;
 
/** Select finite element type for integration on domain based on problem
 * dimension
 */
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
 
/** Select finite element type for integrate on boundary based on problem
 * dimension
 */
using BoundaryEle =
    PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
/** Operators used to assemble domain integrals
 */
using DomainEleOp = DomainEle::UserDataOperator;
 
/** Operators used to assemble boundary integrals
 */
using BoundaryEleOp = BoundaryEle::UserDataOperator;
 
/** Linear forms used to integrate body forces
 */
using DomainNaturalBC =
    NaturalBC<DomainEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
 
/** Select linear froms reading data from blockest (e.g. "BODY_FORCE") and
 * applying body force.
 */
using OpBodyForce =
    DomainNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
 
// Select natural BC boundary condition
 
// Use natural boundary conditions implemented with adv-1 example which includes
// spring BC
#include <ContactNaturalBC.hpp>
 
/** Use Gauss quadrature rule and PETSc assembly to integrate neural boundary
 * conditions. Select linear forms operators.
 */
using BoundaryRhsBCs =
    NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
 
/* Use specialization from adv-1 integrating boundary conditions on forces and
 * with springs. OP is used to integrate the right hand side
 */
using OpBoundaryRhsBCs =
    BoundaryRhsBCs::OpFlux<ContactOps::BoundaryBCs, 1, SPACE_DIM>;
 
/** Use Gauss quadrature rule and PETSc assembly to integrate neural boundary
 * conditions. Select bi-linear forms operators.
 */
using BoundaryLhsBCs =
    NaturalBC<BoundaryEleOp>::Assembly<PETSC>::BiLinearForm<GAUSS>;
 
/** Use specialization from adv-1 integrating boundary conditions on forces and
 * with springs
 */
using OpBoundaryLhsBCs =
    BoundaryLhsBCs::OpFlux<ContactOps::BoundaryBCs, 1, SPACE_DIM>;
 
// Note: We create only skin post-processing mesh.
 
template <int DIM>
struct ElementsAndOps; //< Select finite elements and operators used for
                       // postprocessing and contact space
 
template <> struct ElementsAndOps<2> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBase<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBase<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
  static constexpr FieldSpace CONTACT_SPACE = HCURL;
};
 
template <> struct ElementsAndOps<3> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBase<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBase<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
  static constexpr FieldSpace CONTACT_SPACE = HDIV;
};
// Define contact space by dimension
constexpr FieldSpace ElementsAndOps<2>::CONTACT_SPACE;
constexpr FieldSpace ElementsAndOps<3>::CONTACT_SPACE;
constexpr FieldSpace CONTACT_SPACE = ElementsAndOps<SPACE_DIM>::CONTACT_SPACE;
 
using PostProcEleDomain = ElementsAndOps<SPACE_DIM>::PostProcEleDomain;
using SideEle = ElementsAndOps<SPACE_DIM>::SideEle;
using PostProcEleBdy = ElementsAndOps<SPACE_DIM>::PostProcEleBdy;
 
double scale = 1.;
#include <HenckyOps.hpp>
#include <ADOLCPlasticity.hpp>
#include <ADOLCPlasticityLargeStrain.hpp>
#include <ADOLCPlasticityMaterialModels.hpp>
using namespace HenckyOps;
using namespace ADOLCPlasticity;
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
namespace ContactOps {
 
double cn_contact = 0.1;
 
}; // namespace ContactOps
 
#include <ContactOps.hpp>
#endif // ADD_CONTACT
 
struct PlasticProblem {
 
  PlasticProblem(MoFEM::Interface &m_field) : mField(m_field) {}
 
  MoFEMErrorCode runProblem();
 
private:
  MoFEM::Interface &mField;
 
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode boundaryCondition();
  MoFEMErrorCode assembleSystem();
  MoFEMErrorCode solveSystem();
  MoFEMErrorCode gettingNorms();
  MoFEMErrorCode outputResults();
  MoFEMErrorCode checkResults();
 
  int approxOrder = 2;
  int geomOrder = 2;
  PetscBool largeStrain = PETSC_FALSE;  ///< Large strain flag
  FieldApproximationBase approxBase = AINSWORTH_LEGENDRE_BASE;
  boost::shared_ptr<ClosestPointProjection> cpPtr; ///< Closest point
                                                   ///< projection
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
  boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
#endif
#endif // ADD_CONTACT
};
 
//! [Run problem]
MoFEMErrorCode PlasticProblem::runProblem() {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR gettingNorms();
  CHKERR outputResults();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Read mesh]
MoFEMErrorCode PlasticProblem::readMesh() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
 
  if (simple->getDim() != SPACE_DIM)
    SETERRQ2(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
             "Wrong dimension of mesh %d != %d", simple->getDim(), SPACE_DIM);
 
  // Check if mesh has boundary conditions tag and set yield strength as elastic
  Tag th_boundary_conditions;
  ErrorCode rval_bc = mField.get_moab().tag_get_handle(
      "BOUNDARY_CONDITIONS", 1, MB_TYPE_INTEGER, th_boundary_conditions,
      MB_TAG_SPARSE);
 
  // Set yield strength as elastic if boundary conditions tag is set
  if (rval_bc == MB_SUCCESS) {
    
    auto yield_strength_compressive = std::numeric_limits<double>::max();
    auto yield_strength_tension = std::numeric_limits<double>::max();
 
    Tag th_compressive_yield_stress;
    Tag th_tension_yield_stress;
 
    ErrorCode rval_compressive = mField.get_moab().tag_get_handle(
        "Yield_Strength_C", 1, MB_TYPE_DOUBLE, th_compressive_yield_stress,
        MB_TAG_SPARSE);
 
    ErrorCode rval_tension = mField.get_moab().tag_get_handle(
        "Yield_Strength_T", 1, MB_TYPE_DOUBLE, th_tension_yield_stress,
        MB_TAG_SPARSE);
 
    Range fe_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, fe_ents);
 
    for (auto fe_ent : fe_ents) {
      Range nodes;
      // get nodes of the element
      CHKERR mField.get_moab().get_adjacencies(&fe_ent, 1, 0, false, nodes);
 
      std::vector<double> v_bc_val;
      for (auto node : nodes) {
        int bc_val;
        CHKERR mField.get_moab().tag_get_data(th_boundary_conditions, &node, 1,
                                              static_cast<void *>(&bc_val));
        v_bc_val.push_back(bc_val);
      }
 
      // Set yield strength as elastic if boundary conditions is FIX_X or
      // CONTACT
      if (std::find_if(v_bc_val.begin(), v_bc_val.end(), [](int val) {
            return val == 5 || val == 20;
          }) != v_bc_val.end()) {
 
        if (rval_compressive == MB_SUCCESS) {
 
          CHKERR mField.get_moab().tag_set_data(th_compressive_yield_stress,
                                                &fe_ent, 1,
                                                &yield_strength_compressive);
        } else {
          MOFEM_LOG("ADOLCPlasticityWold", Sev::warning)
              << "Yield_Strength_C tag does not exist using default value";
        }
        if (rval_tension == MB_SUCCESS) {
 
          CHKERR mField.get_moab().tag_set_data(
              th_tension_yield_stress, &fe_ent, 1, &yield_strength_tension);
        } else {
          MOFEM_LOG("ADOLCPlasticityWold", Sev::warning)
              << "Yield_Strength_T tag does not exist using default value";
        }
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
//! [Read mesh]
 
//! [Set up problem]
MoFEMErrorCode PlasticProblem::setupProblem() {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
  PetscInt choice_base_value = AINSWORTH;
  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
                              LASBASETOPT, &choice_base_value, PETSC_NULL);
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &approxOrder, PETSC_NULL);
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geomOrder,
                            PETSC_NULL);
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strain", &largeStrain,
                            PETSC_NULL);
 
  switch (choice_base_value) {
  case AINSWORTH:
    approxBase = AINSWORTH_LEGENDRE_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set AINSWORTH_LEGENDRE_BASE for displacements";
    break;
  case DEMKOWICZ:
    approxBase = DEMKOWICZ_JACOBI_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set DEMKOWICZ_JACOBI_BASE for displacements";
    break;
  default:
    approxBase = LASTBASE;
    break;
  }
 
  // Add field
  CHKERR simple->addDomainField("U", H1, approxBase, SPACE_DIM);
  CHKERR simple->addBoundaryField("U", H1, approxBase, SPACE_DIM);
  CHKERR simple->addDataField("GEOMETRY", H1, approxBase, SPACE_DIM);
 
#ifdef ADD_CONTACT
  CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
                                SPACE_DIM);
  CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
                                  SPACE_DIM);
  auto get_skin = [&]() {
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
    Skinner skin(&mField.get_moab());
    Range skin_ents;
    CHKERR skin.find_skin(0, body_ents, false, skin_ents);
    return skin_ents;
  };
 
  auto filter_blocks = [&](auto skin) {
    bool is_contact_block = false;
    Range contact_range;
    for (auto m :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % "CONTACT").str()
 
                 ))
 
    ) {
      is_contact_block =
          true; ///< bloks interation is collectibe, so that is set irrespective
                ///< if there are enerities in given rank or not in the block
      MOFEM_LOG("CONTACT", Sev::inform)
          << "Find contact block set:  " << m->getName();
      auto meshset = m->getMeshset();
      Range contact_meshset_range;
      CHKERR mField.get_moab().get_entities_by_dimension(
          meshset, SPACE_DIM - 1, contact_meshset_range, true);
 
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          contact_meshset_range);
      contact_range.merge(contact_meshset_range);
    }
    if (is_contact_block) {
      MOFEM_LOG("SYNC", Sev::inform)
          << "Nb entities in contact surface: " << contact_range.size();
      MOFEM_LOG_SYNCHRONISE(mField.get_comm());
      skin = intersect(skin, contact_range);
    }
    return skin;
  };
 
  auto filter_true_skin = [&](auto skin) {
    Range boundary_ents;
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &boundary_ents);
    return boundary_ents;
  };
 
  auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
#ifdef ADD_CONTACT
  CHKERR simple->setFieldOrder("SIGMA", 0);
  CHKERR simple->setFieldOrder("SIGMA", approxOrder - 1, &boundary_ents);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn_contact",
                               &ContactOps::cn_contact, PETSC_NULL);
  MOFEM_LOG("CONTACT", Sev::inform) << "cn_contact " << ContactOps::cn_contact;
 
#endif
 
#ifdef PYTHON_SDF
  sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
  CHKERR sdfPythonPtr->sdfInit("sdf.py");
  ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
#endif
#endif
 
  CHKERR simple->setFieldOrder("U", approxOrder);
  CHKERR simple->setFieldOrder("GEOMETRY", geomOrder);
  CHKERR simple->setUp();
 
  auto project_ho_geometry = [&]() {
    Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
    return mField.loop_dofs("GEOMETRY", ent_method);
  };
  PetscBool project_geometry = PETSC_TRUE;
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
                             &project_geometry, PETSC_NULL);
  if (project_geometry) {
    CHKERR project_ho_geometry();
  }
 
  //! [Material models selection]
 
  /**
   * FIXME: Here only array of material models is initalized. Each model has
   * unique set of the ADOL-C tags. Pointer is attached based on block name to
   * which entity belongs. That will enable heterogeneity of the model, in
   * addition of heterogeneity of the material properties.
   */
 
  enum MaterialModel {
    VonMisses,
    VonMissesPlaneStress,
    Paraboloidal,
    HeterogeneousParaboloidal,
    LastModel
  };
  const char *list_materials[LastModel] = {"VonMisses", "VonMissesPlaneStress",
                                           "Paraboloidal","HeterogeneousParaboloidal"};
  PetscInt choice_material = VonMisses;
  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-material", list_materials,
                              LastModel, &choice_material, PETSC_NULL);
 
  switch (choice_material) {
  case VonMisses:
    cpPtr = createMaterial<J2Plasticity<3>>();
    break;
  case VonMissesPlaneStress:
    if (SPACE_DIM != 2)
      SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
              "VonMissesPlaneStrain is only for 2D case");
    cpPtr = createMaterial<J2Plasticity<2>>();
    break;
  case Paraboloidal:
    cpPtr = createMaterial<ParaboloidalPlasticity>();
    break;
  case HeterogeneousParaboloidal:
    cpPtr = createMaterial<HeterogeneousParaboloidalPlasticity>();
    break;
  default:
    SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Wrong material model");
  }
 
  if (approxBase == DEMKOWICZ_JACOBI_BASE && SPACE_DIM == 2)
    cpPtr->integrationRule = [](int, int, int p) { return 2 * (p - 2); };
 
  //! [Material models selection]
 
  MoFEMFunctionReturn(0);
}
//! [Set up problem]
 
//! [Boundary condition]
MoFEMErrorCode PlasticProblem::boundaryCondition() {
  MoFEMFunctionBegin;
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
  auto time_scale = boost::make_shared<TimeScale>();
 
  auto rule = [](int, int, int p) { return 2 * p; };
  CHKERR pipeline_mng->setDomainRhsIntegrationRule(cpPtr->integrationRule);
  CHKERR pipeline_mng->setDomainLhsIntegrationRule(cpPtr->integrationRule);
  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule);
  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(rule);
 
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pipeline_mng->getOpBoundaryLhsPipeline(), {HDIV}, "GEOMETRY");
  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
      new OpSetHOWeightsOnSubDim<SPACE_DIM>());
 
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pipeline_mng->getOpBoundaryRhsPipeline(), {HDIV}, "GEOMETRY");
  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
      new OpSetHOWeightsOnSubDim<SPACE_DIM>());
 
  // Add Natural BCs to RHS
  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
      pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
      Sev::inform);
  // Add Natural BCs to LHS
  CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
      pipeline_mng->getOpBoundaryLhsPipeline(), mField, "U", Sev::inform);
 
#ifdef ADD_CONTACT
  CHKERR
  ContactOps::opFactoryBoundaryLhs<SPACE_DIM, PETSC, GAUSS, BoundaryEleOp>(
      pipeline_mng->getOpBoundaryLhsPipeline(), "SIGMA", "U");
  CHKERR
  ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEle>(
      mField, pipeline_mng->getOpBoundaryLhsPipeline(),
      simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
      cpPtr->integrationRule);
 
  CHKERR
  ContactOps::opFactoryBoundaryRhs<SPACE_DIM, PETSC, GAUSS, BoundaryEleOp>(
      pipeline_mng->getOpBoundaryRhsPipeline(), "SIGMA", "U");
#endif // ADD_CONTACT
 
  //! Add body froces
  CHKERR DomainNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
      pipeline_mng->getOpDomainRhsPipeline(), mField, "U", {time_scale},
      "BODY_FORCE", Sev::inform);
 
  // Essential BC
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
                                           "U", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
                                           "U", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
                                           "U", 2, 2);
#ifdef ADD_CONTACT
  for (auto b : {"FIX_X", "REMOVE_X"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                             "SIGMA", 0, 0, false, true);
  for (auto b : {"FIX_Y", "REMOVE_Y"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                             "SIGMA", 1, 1, false, true);
  for (auto b : {"FIX_Z", "REMOVE_Z"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                             "SIGMA", 2, 2, false, true);
  for (auto b : {"FIX_ALL", "REMOVE_ALL"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                             "SIGMA", 0, 3, false, true);
  CHKERR bc_mng->removeBlockDOFsOnEntities(
      simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
#endif
 
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "U");
 
  MoFEMFunctionReturn(0);
}
//! [Boundary condition]
 
//! [Push operators to pipeline]
MoFEMErrorCode PlasticProblem::assembleSystem() {
  MoFEMFunctionBegin;
  auto *simple = mField.getInterface<Simple>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  // assemble operator to the right hand side
  auto add_domain_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
    // push forward finite element bases from reference to physical element
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
                                                          "GEOMETRY");
    // create local common data
    auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
 
    if (largeStrain) {
      CHKERR
      opFactoryDomainHenckyStrainLhs<SPACE_DIM, GAUSS, DomainEleOp>(
          mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr);
      using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
          PETSC>::template BiLinearForm<GAUSS>;
 
      using OpKPiola =
          typename B::template OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
 
      pip.push_back(
          new OpKPiola("U", "U", common_data_ptr->getMatTangentPtr()));
    } else {
      pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
          "U", common_data_ptr->getGradAtGaussPtsPtr()));
      CHKERR opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
          mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
    // push forward finite element bases from reference to physical element
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
                                                          "GEOMETRY");
    // create local common data
    auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
 
    if (largeStrain) {
      CHKERR
      opFactoryDomainHenckyStrainRhs<SPACE_DIM, GAUSS, DomainEleOp>(
          mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr);
      using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
          PETSC>::template LinearForm<GAUSS>;
      using OpInternalForcePiola =
          typename B::template OpGradTimesTensor<1, SPACE_DIM, SPACE_DIM>;
      pip.push_back(
          new OpInternalForcePiola("U", common_data_ptr->getStressMatrixPtr()));
    } else {
      pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
          "U", common_data_ptr->getGradAtGaussPtsPtr()));
      CHKERR opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
          mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr);
    }
 
#ifdef ADD_CONTACT
    CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
        pip, "SIGMA", "U");
#endif // ADD_CONTACT
    MoFEMFunctionReturn(0);
  };
 
  // Push operators to the left hand side pipeline. Indicate that domain (i.e.
  // volume/interior) element is used.
  CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
  // Push operators to the right hand side pipeline.
  CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
 
  MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
 
/**
 * @brief Monitor solution
 *
 * This functions is called by TS solver at the end of each step. It is used
 * to output results to the hard drive.
 */
struct Monitor : public FEMethod {
  Monitor(SmartPetscObj<DM> dm,
          std::pair<boost::shared_ptr<PostProcEleDomain>,
                    boost::shared_ptr<PostProcEleBdy>>
              pair_post_proc_fe,
          boost::shared_ptr<DomainEle> reaction_fe,
          std::vector<boost::shared_ptr<ScalingMethod>> smv)
      : dM(dm), reactionFE(reaction_fe), vecOfTimeScalingMethods(smv) {
    fRes = createDMVector(dM);
    domainPostProcFe = pair_post_proc_fe.first;
    skinPostProcFe = pair_post_proc_fe.second;
 
    MoFEM::Interface *m_field_ptr;
    CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
#ifdef ADD_CONTACT
    auto get_integrate_traction = [&]() {
      auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
      auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
      CHK_THROW_MESSAGE(
          (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
              integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
          "Apply transform");
      // We have to integrate on curved face geometry, thus integration weight
      // have to adjusted.
      integrate_traction->getOpPtrVector().push_back(
          new OpSetHOWeightsOnSubDim<SPACE_DIM>());
      integrate_traction->getRuleHook = [](int, int, int approx_order) {
        return 2 * approx_order + 2 - 1;
      };
 
      CHK_THROW_MESSAGE(
          (ContactOps::opFactoryCalculateTraction<SPACE_DIM, GAUSS,
                                                  BoundaryEleOp>(
              integrate_traction->getOpPtrVector(), "SIGMA", 0)),
          "push operators to calculate traction");
 
      return integrate_traction;
    };
 
    integrateTraction = get_integrate_traction();
#endif
  }
  MoFEMErrorCode postProcess() {
    MoFEMFunctionBegin;
    int save_every_nth_step = 1;
    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every_nth_step",
                              &save_every_nth_step, PETSC_NULL);
#ifdef ADD_CONTACT
    auto print_traction = [&](const std::string msg) {
      MoFEMFunctionBegin;
      MoFEM::Interface *m_field_ptr;
      CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
      if (!m_field_ptr->get_comm_rank()) {
        const double *t_ptr;
        CHKERR VecGetArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
        MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e %3.4e %3.4e %3.4e",
                    msg.c_str(), ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
        CHKERR VecRestoreArrayRead(ContactOps::CommonData::totalTraction,
                                   &t_ptr);
      }
      MoFEMFunctionReturn(0);
    };
#endif
 
    auto make_vtk = [&]() {
      MoFEMFunctionBegin;
      if (domainPostProcFe) {
        CHKERR DMoFEMLoopFiniteElements(dM, "dFE", domainPostProcFe,
                                        getCacheWeakPtr());
        CHKERR domainPostProcFe->writeFile(
            "out_plastic_" + boost::lexical_cast<std::string>(ts_step) +
            ".h5m");
      }
      if (skinPostProcFe) {
        CHKERR DMoFEMLoopFiniteElements(dM, "bFE", skinPostProcFe,
                                        getCacheWeakPtr());
        CHKERR skinPostProcFe->writeFile(
            "out_skin_plastic_" + boost::lexical_cast<std::string>(ts_step) +
            ".h5m");
      }
      MoFEMFunctionReturn(0);
    };
 
    if (!(ts_step % save_every_nth_step)) {
      CHKERR make_vtk();
    }
    if (reactionFE) {
      CHKERR VecZeroEntries(fRes);
      reactionFE->f = fRes;
      CHKERR DMoFEMLoopFiniteElements(dM, "dFE", reactionFE, getCacheWeakPtr());
      CHKERR VecAssemblyBegin(fRes);
      CHKERR VecAssemblyEnd(fRes);
      CHKERR VecGhostUpdateBegin(fRes, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(fRes, ADD_VALUES, SCATTER_REVERSE);
 
      MoFEM::Interface *m_field_ptr;
      CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
      CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
          *m_field_ptr, reactionFE, fRes)();
 
      double nrm;
      CHKERR VecNorm(fRes, NORM_2, &nrm);
      MOFEM_LOG("PlasticPrb", Sev::verbose)
          << "Residual norm " << nrm << " at time step " << ts_step;
    }
 
#ifdef ADD_CONTACT
    auto calculate_traction = [&] {
      MoFEMFunctionBegin;
      CHKERR VecZeroEntries(ContactOps::CommonData::totalTraction);
      CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateTraction);
      CHKERR VecAssemblyBegin(ContactOps::CommonData::totalTraction);
      CHKERR VecAssemblyEnd(ContactOps::CommonData::totalTraction);
      MoFEMFunctionReturn(0);
    };
#endif
 
    auto get_min_max_displacement = [&]() {
      MoFEMFunctionBegin;
      MoFEM::Interface *m_field_ptr;
      CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
 
      std::array<double, 4> a_min = {DBL_MAX, DBL_MAX, DBL_MAX, 0};
      std::array<double, 4> a_max = {-DBL_MAX, -DBL_MAX, -DBL_MAX, 0};
 
      auto get_min_max = [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
        MoFEMFunctionBegin;
        int d = 0;
        for (auto v : field_entity_ptr->getEntFieldData()) {
          a_min[d] = std::min(a_min[d], v);
          a_max[d] = std::max(a_max[d], v);
          ++d;
        }
        MoFEMFunctionReturn(0);
      };
 
      a_min[SPACE_DIM] = 0;
      a_max[SPACE_DIM] = 0;
 
      Range verts;
      CHKERR m_field_ptr->get_field_entities_by_type("U", MBVERTEX, verts);
      CHKERR m_field_ptr->getInterface<FieldBlas>()->fieldLambdaOnEntities(
          get_min_max, "U", &verts);
 
      auto mpi_reduce = [&](auto &a, auto op) {
        std::array<double, 3> a_mpi = {0, 0, 0};
        MPI_Allreduce(a.data(), a_mpi.data(), 3, MPI_DOUBLE, op,
                      m_field_ptr->get_comm());
        return a_mpi;
      };
 
      auto a_min_mpi = mpi_reduce(a_min, MPI_MIN);
      auto a_max_mpi = mpi_reduce(a_max, MPI_MAX);
 
      MOFEM_LOG("PlasticPrb", Sev::inform)
          << "Min displacement " << a_min_mpi[0] << " " << a_min_mpi[1] << " "
          << a_min_mpi[2];
      MOFEM_LOG("PlasticPrb", Sev::inform)
          << "Max displacement " << a_max_mpi[0] << " " << a_max_mpi[1] << " "
          << a_max_mpi[2];
      MoFEMFunctionReturn(0);
    };
    CHKERR get_min_max_displacement();
#ifdef ADD_CONTACT
    CHKERR calculate_traction();
    CHKERR print_traction("Contact force");
#endif
    double scale = 1;
    for (auto s : vecOfTimeScalingMethods) {
      scale *= s->getScale(this->ts_t);
    }
 
    MOFEM_LOG("PlasticPrb", Sev::inform)
        << "Time: " << this->ts_t << " scale: " << scale;
 
    MoFEMFunctionReturn(0);
  }
 
private:
  SmartPetscObj<DM> dM;
  boost::shared_ptr<PostProcEleDomain> domainPostProcFe;
  boost::shared_ptr<PostProcEleBdy> skinPostProcFe;
  boost::shared_ptr<FEMethod> reactionFE;
 
  boost::shared_ptr<BoundaryEle> integrateTraction;
 
  SmartPetscObj<Vec> fRes;
  VecOfTimeScalingMethods vecOfTimeScalingMethods;
};
 
static boost::shared_ptr<TSUpdate> ts_update_ptr = nullptr;
 
//! [Solve]
MoFEMErrorCode PlasticProblem::solveSystem() {
  MoFEMFunctionBegin;
  auto *simple = mField.getInterface<Simple>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto dm = simple->getDM();
  auto time_scale = boost::make_shared<TimeScale>();
  auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
 
  // Setup postprocessing
  auto create_post_proc_fe = [dm, this, simple]() {
    //! [DG projection]
    // Post-process domain, i.e. volume elements. If 3d case creates stresses
    // are evaluated at points which are trace of the volume element on boundary
    auto post_proc_ele_domain = [this](auto &pip_domain, auto &fe_name) {
      // Push bases on reference element to the phyiscal element
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          pip_domain, {H1, HDIV}, "GEOMETRY");
 
      // calculate displacement gradients at nodes of post processing mesh. For
      // gradient DG projection is obsolete, since displacements can be
      // evaluated at arbitrary points.
      auto grad_ptr = boost::make_shared<MatrixDouble>();
      pip_domain.push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
                                                                   grad_ptr));
 
      // Create operator to run (neted) pipeline in operator. InteriorElem
      // element will have iteration rule and bases as domain element used to
      // solve problem, and remember history variables.
      using InteriorElem = DomainEle;
      auto op_this = new OpLoopThis<InteriorElem>(mField, fe_name, Sev::noisy);
      // add interior element to post-processing pipeline
      pip_domain.push_back(op_this);
 
      // get pointer to interior element, which is in essence pipeline which
      // pipeline.
      auto fe_physics = op_this->getThisFEPtr();
 
      auto evaluate_stress_on_physical_element = [&]() {
        // evaluate stress and plastic strain at Gauss points
        fe_physics->getRuleHook = cpPtr->integrationRule;
        CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
            fe_physics->getOpPtrVector(), {H1});
        auto common_data_ptr =
            boost::make_shared<ADOLCPlasticity::CommonData>();
        // note that gradient of displacements is evaluate again, at
        // physical nodes
        if (largeStrain) {
          CHKERR
          opFactoryDomainHenckyStrainRhs<SPACE_DIM, GAUSS, DomainEleOp>(
              mField, "U", fe_physics->getOpPtrVector(), "ADOLCMAT", common_data_ptr, cpPtr);
              
        } else {
          fe_physics->getOpPtrVector().push_back(
              new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
                  "U", common_data_ptr->getGradAtGaussPtsPtr()));
          CHKERR cpPtr->addMatBlockOps(mField, fe_physics->getOpPtrVector(), "ADOLCMAT", Sev::inform);
          fe_physics->getOpPtrVector().push_back(getRawPtrOpCalculateStress<SPACE_DIM, SMALL_STRAIN>(
              mField, common_data_ptr, cpPtr, false));
        }
        return common_data_ptr;
      };
 
      auto dg_projection_froward_and_back = [&](auto &common_data_ptr) {
        // here we do actual projection of stress and plastic strain to DG space
        int dg_order = approxOrder - 1;
        auto entity_data_l2 =
            boost::make_shared<EntitiesFieldData>(MBENTITYSET);
        auto mass_ptr =
            boost::make_shared<MatrixDouble>(); //< projection operator (mass
                                                // matrix)
        fe_physics->getOpPtrVector().push_back(new OpDGProjectionMassMatrix(
            dg_order, mass_ptr, entity_data_l2, approxBase, L2));
 
        auto coeffs_ptr_stress = boost::make_shared<MatrixDouble>();
        // project stress on DG space on physical element
        fe_physics->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
            common_data_ptr->getStressMatrixPtr(), coeffs_ptr_stress, mass_ptr,
            entity_data_l2, approxBase, L2));
        // project strains plastic strains DG space on physical element
        auto coeffs_ptr_plastic_strain = boost::make_shared<MatrixDouble>();
        fe_physics->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
            common_data_ptr->getPlasticStrainMatrixPtr(),
            coeffs_ptr_plastic_strain, mass_ptr, entity_data_l2, approxBase,
            L2));
 
        // project stress and plastic strain for DG space on post-process
        // element
        pip_domain.push_back(new OpDGProjectionEvaluation(
            common_data_ptr->getStressMatrixPtr(), coeffs_ptr_stress,
            entity_data_l2, approxBase, L2));
        pip_domain.push_back(new OpDGProjectionEvaluation(
            common_data_ptr->getPlasticStrainMatrixPtr(),
            coeffs_ptr_plastic_strain, entity_data_l2, approxBase, L2));
      };
 
      auto common_data_ptr = evaluate_stress_on_physical_element();
      dg_projection_froward_and_back(common_data_ptr);
 
      return boost::make_tuple(grad_ptr, common_data_ptr->getStressMatrixPtr(),
                               common_data_ptr->getPlasticStrainMatrixPtr());
    };
    //! [DG projection]
 
    // Create tags on post-processing mesh, i.e. those tags are visible in
    // Preview
    auto add_post_proc_map = [&](auto post_proc_fe, auto u_ptr, auto grad_ptr,
                                 auto stress_ptr, auto plastic_strain_ptr,
                                 auto contact_stress_ptr, auto X_ptr) {
      using OpPPMapSPACE_DIM = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
      post_proc_fe->getOpPtrVector().push_back(
 
          new OpPPMapSPACE_DIM(
 
              post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
              {},
 
              {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
 
              {{"GRAD", grad_ptr}, {"SIGMA", contact_stress_ptr}},
 
              {}
 
              )
 
      );
 
      using OpPPMap3D = OpPostProcMapInMoab<3, 3>;
      if (largeStrain) {
        post_proc_fe->getOpPtrVector().push_back(
 
            new OpPPMap3D(
 
                post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
                {},
 
                {},
 
                {{"PK1", stress_ptr}},
 
                {{"PLASTIC_STRAIN", plastic_strain_ptr}}
 
                )
 
        );
      } else {
        post_proc_fe->getOpPtrVector().push_back(
 
            new OpPPMap3D(
 
                post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
                {},
 
                {},
 
                {},
 
                {{"STRESS", stress_ptr}, {"PLASTIC_STRAIN", plastic_strain_ptr}}
 
                )
 
        );
      }
 
      return post_proc_fe;
    };
 
    auto vol_post_proc = [this, simple, post_proc_ele_domain,
                          add_post_proc_map]() {
      PetscBool post_proc_vol = PETSC_FALSE;
      if (SPACE_DIM == 2)
        post_proc_vol = PETSC_TRUE;
 
      CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",
                                 &post_proc_vol, PETSC_NULL);
      if (post_proc_vol == PETSC_FALSE)
        return boost::shared_ptr<PostProcEleDomain>();
      auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
      auto u_ptr = boost::make_shared<MatrixDouble>();
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
      auto X_ptr = boost::make_shared<MatrixDouble>();
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
      auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
#ifdef ADD_CONTACT
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
              "SIGMA", contact_stress_ptr));
#else
      contact_stress_ptr = nullptr;
#endif
      auto [grad_ptr, stress_ptr, plastic_strain_ptr] = post_proc_ele_domain(
          post_proc_fe->getOpPtrVector(), simple->getDomainFEName());
 
      return add_post_proc_map(post_proc_fe, u_ptr, grad_ptr, stress_ptr,
                               plastic_strain_ptr, contact_stress_ptr, X_ptr);
    };
 
    auto skin_post_proc = [this, simple, post_proc_ele_domain,
                           add_post_proc_map]() {
      // create skin of the volume mesh for post-processing,
      // i.e. boundary of the volume mesh
      PetscBool post_proc_skin = PETSC_TRUE;
      if (SPACE_DIM == 2)
        post_proc_skin = PETSC_FALSE;
 
      CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
                                 &post_proc_skin, PETSC_NULL);
 
      if (post_proc_skin == PETSC_FALSE)
        return boost::shared_ptr<PostProcEleBdy>();
 
      auto post_proc_fe = boost::make_shared<PostProcEleBdy>(mField);
      auto u_ptr = boost::make_shared<MatrixDouble>();
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
      auto X_ptr = boost::make_shared<MatrixDouble>();
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
      auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
      auto op_loop_side =
          new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
      auto [grad_ptr, stress_ptr, plastic_strain_ptr] = post_proc_ele_domain(
          op_loop_side->getOpPtrVector(), simple->getDomainFEName());
#ifdef ADD_CONTACT
      op_loop_side->getOpPtrVector().push_back(
          new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
              "SIGMA", contact_stress_ptr));
#else
      contact_stress_ptr = nullptr;
#endif
      post_proc_fe->getOpPtrVector().push_back(op_loop_side);
 
      return add_post_proc_map(post_proc_fe, u_ptr, grad_ptr, stress_ptr,
                               plastic_strain_ptr, contact_stress_ptr, X_ptr);
    };
 
    return std::make_pair(vol_post_proc(), skin_post_proc());
  };
 
  auto create_reaction_fe = [&]() {
    auto fe_ptr = boost::make_shared<DomainEle>(mField);
    fe_ptr->getRuleHook = cpPtr->integrationRule;
 
    auto &pip = fe_ptr->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
    auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
 
    if (largeStrain) {
      CHKERR
      opFactoryDomainHenckyStrainRhs<SPACE_DIM, GAUSS, DomainEleOp>(
          mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr);
      using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
          PETSC>::template LinearForm<GAUSS>;
      using OpInternalForcePiola =
          typename B::template OpGradTimesTensor<1, SPACE_DIM, SPACE_DIM>;
      pip.push_back(
          new OpInternalForcePiola("U", common_data_ptr->getStressMatrixPtr()));
    } else {
      pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
          "U", common_data_ptr->getGradAtGaussPtsPtr()));
      CHKERR opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
          mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr);
    }
 
    return fe_ptr;
  };
 
  auto add_extra_finite_elements_to_ksp_solver_pipelines = [&]() {
    MoFEMFunctionBegin;
 
    auto pre_proc_ptr = boost::make_shared<FEMethod>();
    auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
    auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
 
    auto get_bc_hook_rhs = [this, pre_proc_ptr, time_scale]() {
      MoFEMFunctionBeginHot;
      CHKERR EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_ptr,
                                                       {time_scale}, false)();
      MoFEMFunctionReturnHot(0);
    };
 
    pre_proc_ptr->preProcessHook = get_bc_hook_rhs;
 
    auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
      MoFEMFunctionBegin;
      CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
          mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
          mField, post_proc_rhs_ptr, 1.)();
      MoFEMFunctionReturn(0);
    };
    auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
      MoFEMFunctionBegin;
      CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
          mField, post_proc_lhs_ptr, 1.)();
      MoFEMFunctionReturn(0);
    };
    post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
    post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs;
 
    // This is low level pushing finite elements (pipelines) to solver
    auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
    ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
    ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
    ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
    ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
    MoFEMFunctionReturn(0);
  };
 
  // Add extra finite elements to SNES solver pipelines to resolve essential
  // boundary conditions
  CHKERR add_extra_finite_elements_to_ksp_solver_pipelines();
 
  auto create_monitor_fe = [dm, time_scale](auto &&post_proc_fe,
                                            auto &&reaction_fe) {
    return boost::make_shared<Monitor>(
        dm, post_proc_fe, reaction_fe,
        std::vector<boost::shared_ptr<ScalingMethod>>{time_scale});
  };
 
  //! [Set up post-step]
  auto set_up_post_step = [&](auto ts) {
    MoFEMFunctionBegin;
 
    // create finite element (pipeline) and populate it with operators to
    // update history variables
    auto create_update_ptr = [&]() {
      auto fe_ptr = boost::make_shared<DomainEle>(mField);
      fe_ptr->getRuleHook = cpPtr->integrationRule;
      AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(fe_ptr->getOpPtrVector(),
                                                     {H1});
      auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
 
      if (largeStrain) {
        CHKERR opFactoryDomainHenckyStrainRhs<SPACE_DIM, GAUSS,
                                              DomainEleOp>(
            mField, "U", fe_ptr->getOpPtrVector(), "ADOLCMAT", common_data_ptr,
            cpPtr);
      } else {
        fe_ptr->getOpPtrVector().push_back(
            new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
                "U", common_data_ptr->getGradAtGaussPtsPtr()));
        CHKERR cpPtr->addMatBlockOps(mField, fe_ptr->getOpPtrVector(),
                                     "ADOLCMAT", Sev::noisy);
        fe_ptr->getOpPtrVector().push_back(
            getRawPtrOpCalculateStress<SPACE_DIM, SMALL_STRAIN>(mField, common_data_ptr,
                                                  cpPtr, false));
      }
 
      CHK_THROW_MESSAGE(
          opFactoryDomainUpdate<SPACE_DIM>(mField, fe_ptr->getOpPtrVector(),
                                           "ADOLCMAT", common_data_ptr, cpPtr),
          "push update ops");
      return fe_ptr;
    };
 
    // ts_update_ptr is global static variable
    ts_update_ptr =
        createTSUpdate(simple->getDomainFEName(), create_update_ptr());
 
    //! [TS post-step function]
    // This is pure "C" function which we can to the TS solver
    auto ts_step_post_proc = [](TS ts) {
      MoFEMFunctionBegin;
      if (ts_update_ptr)
        CHKERR ts_update_ptr->postProcess(ts);
      MoFEMFunctionReturn(0);
    };
    //! [TS post-step function]
 
    // finally set up post-step
    CHKERR TSSetPostStep(ts, ts_step_post_proc);
    MoFEMFunctionReturn(0);
  };
  //! [Set up post-step]
 
  // Set monitor which postprocessing results and saves them to the hard drive
  auto set_up_monitor = [&](auto ts) {
    MoFEMFunctionBegin;
    boost::shared_ptr<FEMethod> null_fe;
    auto monitor_ptr =
        create_monitor_fe(create_post_proc_fe(), create_reaction_fe());
    CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
                               null_fe, monitor_ptr);
    MoFEMFunctionReturn(0);
  };
 
  auto set_section_monitor = [&](auto solver) {
    MoFEMFunctionBegin;
    SNES snes;
    CHKERR TSGetSNES(solver, &snes);
    CHKERR SNESMonitorSet(snes,
                          (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
                                             void *))MoFEMSNESMonitorFields,
                          (void *)(snes_ctx_ptr.get()), nullptr);
    MoFEMFunctionReturn(0);
  };
 
  auto set_up_adapt = [&](auto ts) {
    MoFEMFunctionBegin;
    CHKERR TSAdaptRegister(TSADAPTMOFEM, TSAdaptCreateMoFEM);
    TSAdapt adapt;
    CHKERR TSGetAdapt(ts, &adapt);
    MoFEMFunctionReturn(0);
  };
 
  //! [Create TS]
  auto ts = pipeline_mng->createTSIM();
 
  // Set time solver
  double ftime = 1;
  CHKERR TSSetMaxTime(ts, ftime);
  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
  auto D = createDMVector(simple->getDM());
#ifdef ADD_CONTACT
  ContactOps::CommonData::createTotalTraction(mField);
#endif
  CHKERR TSSetSolution(ts, D);
  CHKERR set_section_monitor(ts);
 
  // Set monitor, step adaptivity, and post-step to update history variables
  CHKERR set_up_monitor(ts);
  CHKERR set_up_post_step(ts);
  CHKERR set_up_adapt(ts);
  CHKERR TSSetFromOptions(ts);
 
  CHKERR TSSolve(ts, NULL);
  //! [Create TS]
 
  CHKERR TSGetTime(ts, &ftime);
 
  PetscInt steps, snesfails, rejects, nonlinits, linits;
  CHKERR TSGetStepNumber(ts, &steps);
  CHKERR TSGetSNESFailures(ts, &snesfails);
  CHKERR TSGetStepRejections(ts, &rejects);
  CHKERR TSGetSNESIterations(ts, &nonlinits);
  CHKERR TSGetKSPIterations(ts, &linits);
  MOFEM_LOG_C("PlasticPrb", Sev::inform,
              "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
              "%d, linits %d",
              steps, rejects, snesfails, ftime, nonlinits, linits);
 
  MoFEMFunctionReturn(0);
}
//! [Solve]
 
//! [Getting norms]
MoFEMErrorCode PlasticProblem::gettingNorms() {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  auto dm = simple->getDM();
 
  auto T = createDMVector(simple->getDM());
  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
                                 SCATTER_FORWARD);
  double nrm2;
  CHKERR VecNorm(T, NORM_2, &nrm2);
  MOFEM_LOG("PlasticPrb", Sev::inform) << "Solution norm " << nrm2;
 
  auto post_proc_norm_fe = boost::make_shared<DomainEle>(mField);
 
  post_proc_norm_fe->getRuleHook = cpPtr->integrationRule;
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      post_proc_norm_fe->getOpPtrVector(), {H1});
 
  enum NORMS { U_NORM_L2 = 0, PIOLA_NORM, LAST_NORM };
  auto norms_vec =
      createVectorMPI(mField.get_comm(),
                      (mField.get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
  CHKERR VecZeroEntries(norms_vec);
 
  auto u_ptr = boost::make_shared<MatrixDouble>();
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
 
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_ptr, norms_vec, U_NORM_L2));
 
  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                  post_proc_norm_fe);
 
  CHKERR VecAssemblyBegin(norms_vec);
  CHKERR VecAssemblyEnd(norms_vec);
 
  MOFEM_LOG_CHANNEL("SELF"); // Clear channel from old tags
  if (mField.get_comm_rank() == 0) {
    const double *norms;
    CHKERR VecGetArrayRead(norms_vec, &norms);
    MOFEM_TAG_AND_LOG("SELF", Sev::inform, "example")
        << "norm_u: " << std::scientific << std::sqrt(norms[U_NORM_L2]);
    CHKERR VecRestoreArrayRead(norms_vec, &norms);
  }
 
  MoFEMFunctionReturn(0);
}
//! [Getting norms]
 
//! [Postprocessing results]
MoFEMErrorCode PlasticProblem::outputResults() {
  MoFEMFunctionBegin;
  PetscInt test_nb = 0;
  PetscBool test_flg = PETSC_FALSE;
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-test", &test_nb, &test_flg);
 
  if (test_flg) {
    auto simple = mField.getInterface<Simple>();
    auto T = createDMVector(simple->getDM());
    CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
                                   SCATTER_FORWARD);
    double nrm2;
    CHKERR VecNorm(T, NORM_2, &nrm2);
    MOFEM_LOG("PlasticPrb", Sev::verbose) << "Regression norm " << nrm2;
    double regression_value = 0;
    switch (test_nb) {
    default:
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong test number.");
      break;
    }
    if (fabs(nrm2 - regression_value) > 1e-2)
      SETERRQ2(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
               "Regression test field; wrong norm value. %6.4e != %6.4e", nrm2,
               regression_value);
  }
  MoFEMFunctionReturn(0);
}
//! [Postprocessing results]
 
//! [Check]
MoFEMErrorCode PlasticProblem::checkResults() {
  MoFEMFunctionBegin;
#ifdef ADD_CONTACT
  ContactOps::CommonData::totalTraction.reset();
#endif
  MoFEMFunctionReturn(0);
}
//! [Check]
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
  Py_Initialize();
  np::initialize();
#endif
#endif // ADD_CONTACT
 
  // Initialisation of MoFEM/PETSc and MOAB data structures
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  // Add logging channel for example
  auto core_log = logging::core::get();
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "PlasticPrb"));
  LogManager::setLog("PlasticPrb");
  MOFEM_LOG_TAG("PlasticPrb", "PlasticPrb");
  MOFEM_LOG("PlasticPrb", Sev::inform) << "SPACE_DIM " << SPACE_DIM;
#ifdef ADD_CONTACT
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
  LogManager::setLog("CONTACT");
  MOFEM_LOG_TAG("CONTACT", "Contact");
#endif // ADD_CONTACT
 
  try {
 
    //! [Register MoFEM discrete manager in PETSc]
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    //! [Register MoFEM discrete manager in PETSc
 
    //! [Create MoAB]
    moab::Core mb_instance;              ///< mesh database
    moab::Interface &moab = mb_instance; ///< mesh database interface
    //! [Create MoAB]
 
    //! [Create MoFEM]
    MoFEM::Core core(moab);           ///< finite element database
    MoFEM::Interface &m_field = core; ///< finite element database interface
    //! [Create MoFEM]
 
    //! [PlasticProblem]
    PlasticProblem ex(m_field);
    CHKERR ex.runProblem();
    //! [PlasticProblem]
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
  if (Py_FinalizeEx() < 0) {
    exit(120);
  }
#endif
#endif // ADD_CONTACT
  return 0;
}