thermoplastic.cpp

Go to the documentation of this file.

/**
 * \file thermoplastic.cpp
 * \example thermoplastic.cpp
 *
 * Thermoplasticity in 2d and 3d
 *
 */
 
/* The above code is a preprocessor directive in C++ that checks if the macro
"EXECUTABLE_DIMENSION" has been defined. If it has not been defined, it replaces
the " */
#ifndef EXECUTABLE_DIMENSION
  #define EXECUTABLE_DIMENSION 3
#endif
 
#undef ADD_CONTACT
 
#include <MoFEM.hpp>
#include <MatrixFunction.hpp>
#include <IntegrationRules.hpp>
#include <cstdlib>
extern "C" {
#include <phg-quadrule/quad.h>
}
 
using namespace MoFEM;
 
template <int DIM> struct ElementsAndOps;
 
template <> struct ElementsAndOps<2> {
  using DomainEle = PipelineManager::FaceEle;
  using BoundaryEle = PipelineManager::EdgeEle;
  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
  static constexpr FieldSpace CONTACT_SPACE = HCURL;
};
 
template <> struct ElementsAndOps<3> {
  using DomainEle = PipelineManager::VolEle;
  using BoundaryEle = PipelineManager::FaceEle;
  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
  static constexpr FieldSpace CONTACT_SPACE = HDIV;
};
 
constexpr int SPACE_DIM =
    EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
 
constexpr bool IS_LARGE_STRAINS =
    FINITE_DEFORMATION_FLAG; //< Flag to turn off/on geometric nonlinearities
 
constexpr AssemblyType AT =
    (SCHUR_ASSEMBLE) ? AssemblyType::BLOCK_SCHUR
                     : AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType IT =
    IntegrationType::GAUSS; //< selected integration type
 
constexpr FieldSpace ElementsAndOps<2>::CONTACT_SPACE;
constexpr FieldSpace ElementsAndOps<3>::CONTACT_SPACE;
constexpr FieldSpace CONTACT_SPACE = ElementsAndOps<SPACE_DIM>::CONTACT_SPACE;
 
using EntData = EntitiesFieldData::EntData;
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEle = ElementsAndOps<SPACE_DIM>::BoundaryEle;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
using SkinPostProcEle = PostProcBrokenMeshInMoab<BoundaryEle>;
using SideEle = ElementsAndOps<SPACE_DIM>::SideEle;
using SetPtsData = FieldEvaluatorInterface::SetPtsData;
 
using AssemblyDomainEleOp = FormsIntegrators<DomainEleOp>::Assembly<AT>::OpBase;
 
using ScalerFunTwoArgs = boost::function<double(const double, const double)>;
using ScalerFunThreeArgs =
    boost::function<double(const double, const double, const double)>;
 
namespace {
 
inline double triangleAreaLengthQuality(const double *coords) {
  std::array<std::array<int, 2>, 3> edge_nodes = {
      std::array<int, 2>{0, 1}, std::array<int, 2>{1, 2},
      std::array<int, 2>{2, 0}};
 
  double normal[3];
  CHK_THROW_MESSAGE(MoFEM::Tools::getTriNormal(coords, normal),
                    "calculate triangle normal");
 
  const double area =
      0.5 *
      std::sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
                normal[2] * normal[2]);
 
  double sum_edge_length_sq = 0;
  for (const auto &edge : edge_nodes) {
    for (int dd = 0; dd != 3; ++dd) {
      const double diff = coords[3 * edge[0] + dd] - coords[3 * edge[1] + dd];
      sum_edge_length_sq += diff * diff;
    }
  }
 
  if (sum_edge_length_sq <= std::numeric_limits<double>::epsilon()) {
    return 0;
  }
 
  return 4. * std::sqrt(3.) * area / sum_edge_length_sq;
}
 
inline MoFEMErrorCode getTriangleAreaAndNormal(const double *coords,
                                               double &area,
                                               std::array<double, 3> &normal) {
  MoFEMFunctionBeginHot;
  CHKERR MoFEM::Tools::getTriNormal(coords, normal.data());
  area = 0.5 * std::sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
                         normal[2] * normal[2]);
  MoFEMFunctionReturnHot(0);
}
 
inline MoFEMErrorCode minTriangleQuality(
    MoFEM::Interface &m_field, const Range &tris, double &min_quality, Tag th,
    boost::function<double(double, double)> f =
        [](double a, double b) -> double { return std::min(a, b); }) {
  MoFEMFunctionBegin;
  auto &moab = m_field.get_moab();
  const EntityHandle *conn = nullptr;
  int num_nodes = 0;
  double coords[9];
  for (auto tri : tris) {
    CHKERR moab.get_connectivity(tri, conn, num_nodes, true);
    if (th) {
      CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
    } else {
      CHKERR moab.get_coords(conn, num_nodes, coords);
    }
    double q = triangleAreaLengthQuality(coords);
    if (!std::isnormal(q))
      q = -2;
    min_quality = f(q, min_quality);
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
inline double H_thermal(double H, double omega_h, double temp_0, double temp) {
  return H * (1. - omega_h * (temp - temp_0));
}
 
inline double dH_thermal_dtemp(double H, double omega_h) {
  return -H * omega_h;
}
 
inline double y_0_thermal(double sigmaY, double omega_0, double temp_0,
                          double temp) {
  return sigmaY * (1. - omega_0 * (temp - temp_0));
}
 
inline double dy_0_thermal_dtemp(double sigmaY, double omega_0) {
  return -sigmaY * omega_0;
}
 
inline double Qinf_thermal(double Qinf, double omega_h, double temp_0,
                           double temp) {
  return Qinf * (1. - omega_h * (temp - temp_0));
}
 
inline double dQinf_thermal_dtemp(double Qinf, double omega_h) {
  return -Qinf * omega_h;
}
 
inline double iso_hardening_exp(double tau, double b_iso) {
  return std::exp(-b_iso * tau);
}
 
double exp_C = -6.284;
double exp_D = -1.024;
 
inline double temp_hat(double temp) {
  double JC_ref_temp = 0.;
  double JC_melt_temp = 1650.;
 
  return (temp - JC_ref_temp) / (JC_melt_temp - JC_ref_temp);
}
 
inline double exp_softening(double temp_hat) {
  return std::exp(exp_C * pow(temp_hat, 2) + exp_D * temp_hat);
}
 
/**
 * Isotropic hardening
 */
inline double iso_hardening(double tau, double H, double omega_0, double Qinf,
                            double omega_h, double b_iso, double sigmaY,
                            double temp_0, double temp) {
  return (sigmaY + H * tau + Qinf * (1. - iso_hardening_exp(tau, b_iso))) *
         exp_softening(temp_hat(temp));
}
 
inline double iso_hardening_dtau(double tau, double H, double omega_0,
                                 double Qinf, double omega_h, double b_iso,
                                 double sigmaY, double temp_0, double temp) {
  return (H + Qinf * b_iso * iso_hardening_exp(tau, b_iso)) *
         exp_softening(temp_hat(temp));
}
 
inline double iso_hardening_dtemp(double tau, double H, double omega_0,
                                  double Qinf, double omega_h, double b_iso,
                                  double sigmaY, double temp_0, double temp) {
  double JC_ref_temp = 0.;
  double JC_melt_temp = 1650.;
  double T_hat = temp_hat(temp);
 
  return (sigmaY + H * tau + Qinf * (1. - iso_hardening_exp(tau, b_iso))) *
         (exp_D + 2 * T_hat * exp_C) *
         std::exp(exp_D * T_hat + pow(T_hat, 2) * exp_C) /
         (JC_melt_temp - JC_ref_temp);
}
 
/**
 * Kinematic hardening
 */
template <typename T, int DIM>
inline auto
kinematic_hardening(FTensor::Tensor2_symmetric<T, DIM> &t_plastic_strain,
                    double C1_k) {
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  FTensor::Tensor2_symmetric<double, DIM> t_alpha;
  if (C1_k < std::numeric_limits<double>::epsilon()) {
    t_alpha(i, j) = 0;
    return t_alpha;
  }
  t_alpha(i, j) = C1_k * t_plastic_strain(i, j);
  return t_alpha;
}
 
template <int DIM>
inline auto kinematic_hardening_dplastic_strain(double C1_k) {
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  FTensor::Index<'k', DIM> k;
  FTensor::Index<'l', DIM> l;
  FTensor::Ddg<double, DIM, DIM> t_diff;
  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
  t_diff(i, j, k, l) = C1_k * (t_kd(i, k) ^ t_kd(j, l)) / 4.;
  return t_diff;
}
 
const bool is_large_strains =
    IS_LARGE_STRAINS;                  // PETSC_TRUE; ///< Large strains
PetscBool set_timer = PETSC_FALSE;     ///< Set timer
PetscBool do_eval_field = PETSC_FALSE; ///< Evaluate field
char restart_file[255];
PetscBool restart_flg = PETSC_TRUE;
double restart_time = 0;
int restart_step = 0;
PetscBool is_distributed_mesh = PETSC_TRUE;
 
double scale = 1.;
double default_thermal_conductivity_scale = 1.;
double default_heat_capacity_scale = 1.;
ScalerFunTwoArgs thermal_conductivity_scaling;
ScalerFunTwoArgs heat_capacity_scaling;
ScalerFunThreeArgs inelastic_heat_fraction_scaling;
 
double young_modulus = 115000; ///< Young modulus
double poisson_ratio = 0.31;   ///< Poisson ratio
double sigmaY = 936.2;         ///< Yield stress
double H = 584.3;              ///< Hardening
double visH = 0;               ///< Viscous hardening
double zeta = 5e-2;            ///< regularisation parameter
double Qinf = 174.2;           ///< Saturation yield stress
double b_iso = 63.69;          ///< Saturation exponent
double C1_k = 0;               ///< Kinematic hardening
 
double cn0 = 1;
double cn1 = 1;
 
double default_ref_temp = 20.0;
enum ICType { IC_UNIFORM, IC_GAUSSIAN, IC_LINEAR };
const char *const ICTypes[] = {"uniform", "gaussian", "linear",
                               "ICType",  "IC_",      nullptr};
// CHANGE this from PetscBool to ICType
ICType ic_type = IC_UNIFORM;
double init_temp = 20.0;
double peak_temp = 1000.0;
double width = 10.0; ///< Width of Gaussian distribution
double default_coeff_expansion = 1e-5;
double default_heat_conductivity =
    11.4; // Force / (time temperature )  or Power /
          // (length temperature) // Time unit is hour. force unit kN
double default_heat_capacity = 10*2.332; // length^2/(time^2 temperature) //
                                      // length is millimeter time is hour
 
double omega_0 = 2e-3; ///< flow stress softening
double omega_h = 2e-3; ///< hardening softening
double inelastic_heat_fraction =
    0.9;            ///< fraction of plastic dissipation converted to heat
double temp_0 = 20; ///< reference temperature for thermal softening
 
int order = 2;             ///< Order displacement
int tau_order = order - 2; ///< Order of tau field
int ep_order = order - 1;  ///< Order of ep field
int flux_order = order;    ///< Order of tau field
int T_order = order - 1;   ///< Order of ep field
int geom_order = 2;        ///< Order if fixed.
 
int atom_test = 0;
 
PetscBool order_edge = PETSC_FALSE;
PetscBool order_face = PETSC_FALSE;
PetscBool order_volume = PETSC_FALSE;
 
PetscBool is_quasi_static = PETSC_TRUE;
double rho = 0.0;
double alpha_damping = 0;
double init_dt = 0.05;
double min_dt = 1e-12;
double qual_tol = 0;
double qual_thresh = 0.1;
double edge_growth_thresh = 0.2;
int num_refinement_levels = 0;
 
enum ThermoPlasticityBits {
  COMPUTATION_BIT = 0,
  FLIPPED_BIT = 1,
  FIRST_REF_BIT = 2,
  REFINED_EDGES_BIT = BITREFLEVEL_SIZE - 4, 
  PREVIOUS_BIT = BITREFLEVEL_SIZE - 3,
  VIRGIN_BIT = BITREFLEVEL_SIZE - 2,
  STORAGE_BIT = BITREFLEVEL_SIZE - 1
};
 
BitRefLevel comp_mesh_bit = BitRefLevel().set(COMPUTATION_BIT);
BitRefLevel flipped_bit = BitRefLevel().set(FLIPPED_BIT);
BitRefLevel refined_bit = BitRefLevel().set(FIRST_REF_BIT);
BitRefLevel prev_mesh_bit = BitRefLevel().set(PREVIOUS_BIT);
BitRefLevel virgin_mesh_bit = BitRefLevel().set(VIRGIN_BIT);
BitRefLevel storage_bit = BitRefLevel().set(STORAGE_BIT);
 
auto Gaussian_distribution = [](double init_temp, double peak_temp, double x,
                                double y, double z) {
  double s =
      width /
      std::sqrt(-2. * std::log(0.01 * peak_temp / (peak_temp - init_temp)));
  return (peak_temp - init_temp) * exp(-pow(y, 2) / (2. * pow(s, 2))) +
         init_temp;
};
 
auto linear_distribution = [](double init_temp, double peak_temp, double x,
                              double y, double z) {
  return ((peak_temp - init_temp) / width) * y + (peak_temp + init_temp) / 2.0;
};
 
auto uniform_distribution = [](double init_temp, double peak_temp, double x,
                               double y, double z) { return init_temp; };
 
/**
 * @brief Initialisation function for temperature field
 */
auto init_T = [](double init_temp, double peak_temp, double x, double y,
                 double z) {
  switch (ic_type) {
  case IC_GAUSSIAN:
    return -Gaussian_distribution(init_temp, peak_temp, x, y, z);
  case IC_LINEAR:
    return -linear_distribution(init_temp, peak_temp, x, y, z);
  case IC_UNIFORM:
  default:
    return -uniform_distribution(init_temp, peak_temp, x, y, z);
  }
};
 
#include <HenckyOps.hpp>
#include <PlasticOps.hpp>
#include <PlasticNaturalBCs.hpp>
#include <ThermoElasticOps.hpp>
#include <FiniteThermalOps.hpp>
#include <ThermoPlasticOps.hpp>
#include <ThermalOps.hpp>
 
#include <EdgeFlippingOps.hpp>
#include <SolutionMapping.hpp>
 
#include <ThermalConvection.hpp>
#include <ThermalRadiation.hpp>
 
#ifdef ADD_CONTACT
  #ifdef ENABLE_PYTHON_BINDING
    #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
 
using namespace PlasticOps;
using namespace HenckyOps;
using namespace ThermoElasticOps;
using namespace FiniteThermalOps;
using namespace ThermalOps;
 
int post_processing_counter = 0;
 
auto postProcessHere(MoFEM::Interface &m_field, SmartPetscObj<DM> &dm,
                     std::string output_name, int &counter = *(new int(0))) {
  MoFEMFunctionBegin;
 
  auto create_post_process_elements = [&]() {
    MoFEMFunctionBegin;
    auto pp_fe = boost::make_shared<PostProcEle>(m_field);
 
    auto push_vol_post_proc_ops = [&](auto &pp_fe) {
      MoFEMFunctionBegin;
 
      auto &pip = pp_fe->getOpPtrVector();
 
      auto TAU_ptr = boost::make_shared<VectorDouble>();
      pip.push_back(new OpCalculateScalarFieldValues("TAU", TAU_ptr));
      auto T_ptr = boost::make_shared<VectorDouble>();
      pip.push_back(new OpCalculateScalarFieldValues("T", T_ptr));
 
      auto T_grad_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(
          new OpCalculateScalarFieldGradient<SPACE_DIM>("T", T_grad_ptr));
      auto U_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", U_ptr));
      auto FLUX_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(
          new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", FLUX_ptr));
 
      auto EP_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(
          new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>("EP", EP_ptr));
 
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      pip.push_back(
 
          new OpPPMap(
 
              pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
 
              {{"TAU", TAU_ptr}, {"T", T_ptr}},
 
              {{"U", U_ptr}, {"FLUX", FLUX_ptr}, {"T_GRAD", T_grad_ptr}},
 
              {},
 
              {{"EP", EP_ptr}}
 
              )
 
      );
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR push_vol_post_proc_ops(pp_fe);
 
    if (m_field.get_comm_size() == 1) {
 
      CHKERR DMoFEMLoopFiniteElements(dm, "dFE", pp_fe);
      CHKERR pp_fe->writeFile("out_" + std::to_string(counter) + "_" +
                              output_name + ".h5m");
    } else {
 
      CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dm, "dFE", pp_fe, 0,
                                                  m_field.get_comm_size());
      auto &post_proc_moab = pp_fe->getPostProcMesh();
      auto file_name = "out_" + std::to_string(counter) + "_" + output_name +
                       "_" + std::to_string(m_field.get_comm_rank()) + ".vtk";
      MOFEM_LOG("WORLD", Sev::inform)
          << "Writing file " << file_name << std::endl;
      CHKERR post_proc_moab.write_file(file_name.c_str(), "VTK");
    }
    counter++;
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR create_post_process_elements();
 
  MoFEMFunctionReturn(0);
}
 
auto postProcessPETScHere(MoFEM::Interface &m_field, SmartPetscObj<DM> &dm,
                          Vec sol, std::string output_name) {
  MoFEMFunctionBegin;
 
  auto smart_sol = SmartPetscObj<Vec>(sol, true);
 
  auto create_post_process_elements = [&]() {
    MoFEMFunctionBegin;
    auto pp_fe = boost::make_shared<PostProcEle>(m_field);
 
    auto push_vol_post_proc_ops = [&](auto &pp_fe) {
      MoFEMFunctionBegin;
 
      auto &pip = pp_fe->getOpPtrVector();
 
      auto TAU_ptr = boost::make_shared<VectorDouble>();
      pip.push_back(
          new OpCalculateScalarFieldValues("TAU", TAU_ptr, smart_sol));
      auto T_ptr = boost::make_shared<VectorDouble>();
      pip.push_back(new OpCalculateScalarFieldValues("T", T_ptr, smart_sol));
 
      auto U_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", U_ptr, smart_sol));
      auto FLUX_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateHVecVectorField<3, SPACE_DIM>(
          "FLUX", FLUX_ptr, smart_sol));
 
      auto EP_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
          "EP", EP_ptr, smart_sol));
 
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      pip.push_back(
 
          new OpPPMap(
 
              pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
 
              {{"TAU", TAU_ptr}, {"T", T_ptr}},
 
              {{"U", U_ptr}, {"FLUX", FLUX_ptr}},
 
              {},
 
              {{"EP", EP_ptr}}
 
              )
 
      );
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR push_vol_post_proc_ops(pp_fe);
 
    CHKERR DMoFEMLoopFiniteElements(dm, "dFE", pp_fe);
    CHKERR pp_fe->writeFile("out_" + output_name + ".h5m");
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR create_post_process_elements();
 
  MoFEMFunctionReturn(0);
}
 
auto printOnAllCores(MoFEM::Interface &m_field,
                     const std::string &out_put_string,
                     const auto &out_put_quantity) {
  MoFEMFunctionBegin;
  auto rank = m_field.get_comm_rank();
  auto size = m_field.get_comm_size();
 
  // Loop through all processes and print one by one
  for (int i = 0; i < size; ++i) {
    // If it is this process's turn, print the data
    if (i == rank) {
      std::cout << "Proc " << rank << " " + out_put_string + " "
                << out_put_quantity << std::endl;
    }
    // All processes wait here until the current one has finished
    // printing
    CHKERR PetscBarrier(NULL);
  }
  MoFEMFunctionReturn(0);
};
 
using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;
using OpDomainRhsBCs =
    DomainRhsBCs::OpFlux<PlasticOps::DomainBCs, 1, SPACE_DIM>;
using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;
using OpBoundaryRhsBCs =
    BoundaryRhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::BiLinearForm<IT>;
using OpBoundaryLhsBCs =
    BoundaryLhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
 
//! [Body and heat source]
using DomainNaturalBCRhs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;
using OpBodyForce =
    DomainNaturalBCRhs::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
using OpHeatSource =
    DomainNaturalBCRhs::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 1>;
using DomainNaturalBCLhs =
    NaturalBC<DomainEleOp>::Assembly<AT>::BiLinearForm<IT>;
//! [Body and heat source]
 
//! [Natural boundary conditions]
using BoundaryNaturalBC =
    NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;
using OpForce = BoundaryNaturalBC::OpFlux<NaturalForceMeshsets, 1, SPACE_DIM>;
using OpTemperatureBC =
    BoundaryNaturalBC::OpFlux<NaturalTemperatureMeshsets, 3, SPACE_DIM>;
//! [Natural boundary conditions]
 
//! [Essential boundary conditions (Least square approach)]
using OpEssentialFluxRhs = EssentialBC<BoundaryEleOp>::Assembly<AT>::LinearForm<
    IT>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
using OpEssentialFluxLhs = EssentialBC<BoundaryEleOp>::Assembly<
    AT>::BiLinearForm<IT>::OpEssentialLhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
//! [Essential boundary conditions (Least square approach)]
 
using OpSetTemperatureRhs =
    DomainNaturalBCRhs::OpFlux<SetTargetTemperature, 1, 1>;
using OpSetTemperatureLhs =
    DomainNaturalBCLhs::OpFlux<SetTargetTemperature, 1, 1>;
 
auto save_range = [](moab::Interface &moab, const std::string name,
                     const Range r) {
  MoFEMFunctionBegin;
  auto out_meshset = get_temp_meshset_ptr(moab);
  CHKERR moab.add_entities(*out_meshset, r);
  CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
  MoFEMFunctionReturn(0);
};
 
auto get_string_from_vector = [](auto vec) {
  std::string output_string = "";
  for (auto el : vec) {
    output_string += std::to_string(el) + " ";
  }
  return output_string;
};
 
struct Example;
 
/**
 * @brief Set of functions called by PETSc solver used to refine and update
 * mesh.
 *
 * @note Currently theta method is only handled by this code.
 *
 */
struct TSPrePostProc {
 
  TSPrePostProc() = default;
  virtual ~TSPrePostProc() = default;
 
  /**
   * @brief Used to setup TS solver
   *
   * @param ts
   * @return MoFEMErrorCode
   */
  MoFEMErrorCode tsSetUp(TS ts);
 
  Example *ExRawPtr;
 
private:
  /**
   * @brief Pre process time step
   *
   * Refine mesh and update fields
   *
   * @param ts
   * @return MoFEMErrorCode
   */
  static MoFEMErrorCode tsPreProc(TS ts);
 
  /**
   * @brief Post process time step.
   *
   * Currently that function do not make anything major
   *
   * @param ts
   * @return MoFEMErrorCode
   */
  static MoFEMErrorCode tsPostProc(TS ts);
 
  static MoFEMErrorCode tsPreStage(TS ts);
};
 
static boost::weak_ptr<TSPrePostProc> tsPrePostProc;
 
struct ResizeCtx {
  MoFEM::Interface *m_field;
  // any flags / lists indicating which elements flipped or were added
  PetscBool mesh_changed; // set by your code that flips elements
  std::multimap<double, EntityHandle> el_q_map;
  Range flipped_els;
  std::vector<EntityHandle> new_connectivity;
  Tag th_spatial_coords;
  boost::weak_ptr<TSPrePostProc> ptr;
  // store whatever mapping info you need
};
 
static std::unordered_map<TS, ResizeCtx *> ts_to_rctx;
 
SetEnrichedIntegration::SetEnrichedIntegration() {}
 
MoFEMErrorCode
SetEnrichedIntegration::operator()(ForcesAndSourcesCore *fe_raw_ptr,
                                   int order_row, int order_col,
                                   int order_data) {
  MoFEMFunctionBegin;
 
  constexpr int numNodes = 3;
  constexpr int numEdges = 3;
 
  auto &m_field = fe_raw_ptr->mField;
  auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
  auto fe_handle = fe_ptr->getFEEntityHandle();
 
  auto set_base_quadrature = [&]() {
    MoFEMFunctionBegin;
    int rule = 2 * (order_data + 1);
    if (rule < QUAD_2D_TABLE_SIZE) {
      if (QUAD_2D_TABLE[rule]->dim != 2) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
      }
      if (QUAD_2D_TABLE[rule]->order < rule) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
      }
      const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
      fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
      cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
                  &fe_ptr->gaussPts(0, 0), 1);
      cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
                  &fe_ptr->gaussPts(1, 0), 1);
      cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
                  &fe_ptr->gaussPts(2, 0), 1);
      auto &data = fe_ptr->dataOnElement[H1];
      data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
                                                            false);
      double *shape_ptr =
          &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
      cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
                  1);
      data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
      std::copy(
          Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
          data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
 
    } else {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_base_quadrature();
 
  auto set_gauss_pts = [&](auto &ref_gauss_pts) {
    MoFEMFunctionBegin;
    fe_ptr->gaussPts.swap(ref_gauss_pts);
    const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
    auto &data = fe_ptr->dataOnElement[H1];
    data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3);
    double *shape_ptr =
        &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
    CHKERR ShapeMBTRI(shape_ptr, &fe_ptr->gaussPts(0, 0),
                      &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
    MoFEMFunctionReturn(0);
  };
 
  auto refine_quadrature = [&]() {
    MoFEMFunctionBegin;
 
    moab::Core moab_ref;
    double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
    EntityHandle nodes[numNodes];
    for (int nn = 0; nn != numNodes; nn++)
      CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
    EntityHandle tri;
    CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
    MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
    MoFEM::Interface &m_field_ref = mofem_ref_core;
 
    Range ref_tri(tri, tri);
    Range edges;
    CHKERR m_field_ref.get_moab().get_adjacencies(ref_tri, 1, true, edges,
                                                  moab::Interface::UNION);
    CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
        ref_tri, BitRefLevel().set(0), false, VERBOSE);
 
    Range nodes_at_front;
    for (int nn = 0; nn != numNodes; nn++) {
      EntityHandle ent;
      CHKERR moab_ref.side_element(tri, 0, nn, ent);
      nodes_at_front.insert(ent);
    }
 
    EntityHandle meshset;
    CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
    for (int ee = 0; ee != numEdges; ee++) {
      EntityHandle ent;
      CHKERR moab_ref.side_element(tri, 1, ee, ent);
      CHKERR moab_ref.add_entities(meshset, &ent, 1);
    }
 
    // refine mesh
    auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
    Range ref_edges;
    CHKERR m_field_ref.getInterface<BitRefManager>()
        ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(0),
                                       BitRefLevel().set(), MBEDGE, ref_edges);
 
    // Range tris_to_refine;
    // CHKERR m_field_ref.getInterface<BitRefManager>()
    //     ->getEntitiesByTypeAndRefLevel(
    //         BitRefLevel().set(0), BitRefLevel().set(), MBTRI,
    //         tris_to_refine);
    CHKERR m_ref->addVerticesInTheMiddleOfEdges(ref_edges,
                                                BitRefLevel().set(1));
    CHKERR m_ref->refineTris(ref_tri, BitRefLevel().set(1));
    CHKERR m_field_ref.getInterface<BitRefManager>()
        ->updateMeshsetByEntitiesChildren(meshset, BitRefLevel().set(1),
                                          meshset, MBEDGE, true);
    CHKERR m_field_ref.getInterface<BitRefManager>()
        ->updateMeshsetByEntitiesChildren(meshset, BitRefLevel().set(1),
                                          meshset, MBTRI, true);
 
    // get ref coords
    Range output_ref_tris;
    CHKERR m_field_ref.getInterface<BitRefManager>()
        ->getEntitiesByTypeAndRefLevel(
            BitRefLevel().set(1), BitRefLevel().set(), MBTRI, output_ref_tris);
 
#ifndef NDEBUG
    CHKERR save_range(moab_ref, "ref_tris.vtk", output_ref_tris);
#endif
 
    MatrixDouble ref_coords(output_ref_tris.size(), 9, false);
    int tt = 0;
    for (Range::iterator tit = output_ref_tris.begin();
         tit != output_ref_tris.end(); tit++, tt++) {
      int num_nodes;
      const EntityHandle *conn;
      CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
      CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
    }
 
    auto &data = fe_ptr->dataOnElement[H1];
    const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
    MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
    MatrixDouble &shape_n = data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
    int gg = 0;
    for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
      double *tri_coords = &ref_coords(tt, 0);
      FTensor::Tensor1<double, 3> t_normal;
      CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
      auto det = t_normal.l2();
      for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
        for (int dd = 0; dd != 2; dd++) {
          ref_gauss_pts(dd, gg) = shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
                                  shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
                                  shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
        }
        MOFEM_LOG("REMESHING", Sev::noisy)
            << ref_gauss_pts(0, gg) << ", " << ref_gauss_pts(1, gg);
        ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
      }
    }
 
    int num_nodes;
    const EntityHandle *conn;
    CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
                                               true);
    std::bitset<numNodes> all_nodes;
    for (auto nn = 0; nn != numNodes; ++nn) {
      all_nodes.set(nn);
    }
 
    mapRefCoords[all_nodes.to_ulong()].swap(ref_gauss_pts);
    CHKERR set_gauss_pts(mapRefCoords[all_nodes.to_ulong()]);
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR refine_quadrature();
 
  MoFEMFunctionReturn(0);
}
 
extern PetscErrorCode MyTSResizeTransfer(TS, PetscInt, Vec[], Vec[], void *);
 
struct Example {
 
  Example(MoFEM::Interface &m_field) : mField(m_field) {}
 
  MoFEMErrorCode runProblem();
 
  MoFEM::Interface &mField;
 
private:
  friend PetscErrorCode ::MyTSResizeTransfer(TS, PetscInt, Vec[], Vec[],
                                             void *);
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode createCommonData();
  MoFEMErrorCode initialConditions();
  MoFEMErrorCode mechanicalBC(BitRefLevel bit, BitRefLevel mask);
  MoFEMErrorCode thermalBC(BitRefLevel bit, BitRefLevel mask);
  MoFEMErrorCode
  getElementQuality(std::multimap<double, EntityHandle> &el_q_map,
                    Range &flipped_els,
                    std::vector<EntityHandle> &new_connectivity,
                    bool &do_refine, Tag &th_spatial_coords);
  MoFEMErrorCode edgeFlips(BitRefLevel parent_bit, BitRefLevel child_bit);
  MoFEMErrorCode doEdgeFlips(std::multimap<double, EntityHandle> &el_q_map,
                             Range &flipped_els, Tag &th_spatial_coords,
                             std::vector<EntityHandle> &new_connectivity);
  MoFEMErrorCode doEdgeSplits(bool &refine, bool add_ents);
  MoFEMErrorCode OPs();
  MoFEMErrorCode tsSolve();
 
  boost::shared_ptr<DomainEle> reactionFe;
 
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
 
  boost::shared_ptr<VectorDouble> tempFieldPtr =
      boost::make_shared<VectorDouble>();
  boost::shared_ptr<MatrixDouble> fluxFieldPtr =
      boost::make_shared<MatrixDouble>();
  boost::shared_ptr<MatrixDouble> dispFieldPtr =
      boost::make_shared<MatrixDouble>();
  boost::shared_ptr<MatrixDouble> dispGradPtr =
      boost::make_shared<MatrixDouble>();
  boost::shared_ptr<MatrixDouble> strainFieldPtr =
      boost::make_shared<MatrixDouble>();
  boost::shared_ptr<MatrixDouble> stressFieldPtr =
      boost::make_shared<MatrixDouble>();
  boost::shared_ptr<VectorDouble> plasticMultiplierFieldPtr =
      boost::make_shared<VectorDouble>();
  boost::shared_ptr<MatrixDouble> plasticStrainFieldPtr =
      boost::make_shared<MatrixDouble>();
 
  struct ScaledTimeScale : public MoFEM::TimeScale {
    using MoFEM::TimeScale::TimeScale;
    double getScale(const double time) {
      return scale * MoFEM::TimeScale::getScale(time);
    };
  };
 
#ifdef ADD_CONTACT
  #ifdef ENABLE_PYTHON_BINDING
  boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
  #endif
#endif // ADD_CONTACT
 
  template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
  MoFEMErrorCode opThermoPlasticFactoryDomainRhs(
      MoFEM::Interface &m_field, std::string block_name,
      std::string thermal_block_name, std::string thermoelastic_block_name,
      std::string thermoplastic_block_name, Pip &pip, std::string u,
      std::string ep, std::string tau, std::string temperature) {
    MoFEMFunctionBegin;
 
    using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
        A>::template LinearForm<I>;
    using OpInternalForceCauchy =
        typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
    using OpInternalForcePiola =
        typename B::template OpGradTimesTensor<1, DIM, DIM>;
 
    using P = PlasticityIntegrators<DomainEleOp>;
 
    using H = HenckyOps::HenckyIntegrators<DomainEleOp>;
 
    auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
          common_thermoelastic_ptr, common_thermoplastic_ptr] =
        createCommonThermoPlasticOps<DIM, I, DomainEleOp>(
            m_field, block_name, thermal_block_name, thermoelastic_block_name,
            thermoplastic_block_name, pip, u, ep, tau, temperature, scale,
            thermal_conductivity_scaling, heat_capacity_scaling,
            inelastic_heat_fraction_scaling, Sev::inform);
 
    auto m_D_ptr = common_hencky_ptr->matDPtr;
 
    pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
        ep, common_plastic_ptr->getPlasticStrainDotPtr()));
    pip.push_back(new OpCalculateScalarFieldValuesDot(
        tau, common_plastic_ptr->getPlasticTauDotPtr()));
    pip.push_back(new OpCalculateScalarFieldValues(
        temperature, common_thermoplastic_ptr->getTempPtr()));
    pip.push_back(new OpCalculateHVecVectorField<3, SPACE_DIM>(
        "FLUX", common_thermoplastic_ptr->getHeatFluxPtr()));
    pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
        u, common_plastic_ptr, m_D_ptr, common_thermoplastic_ptr));
 
    pip.push_back(
        new
        typename H::template OpCalculateHenckyThermoPlasticStress<DIM, I, 0>(
            u, common_thermoplastic_ptr->getTempPtr(), common_hencky_ptr,
            common_thermoelastic_ptr->getCoeffExpansionPtr(),
            common_thermoelastic_ptr->getRefTempPtr()));
 
    // Calculate internal forces
    if (common_hencky_ptr) {
      pip.push_back(new OpInternalForcePiola(
          u, common_hencky_ptr->getMatFirstPiolaStress()));
    } else {
      pip.push_back(
          new OpInternalForceCauchy(u, common_plastic_ptr->mStressPtr));
    }
 
    auto inelastic_heat_frac_ptr =
        common_thermoplastic_ptr->getInelasticHeatFractionPtr();
    auto inelastic_heat_fraction =
        [inelastic_heat_frac_ptr](const double, const double, const double) {
          return (*inelastic_heat_frac_ptr);
        };
 
    // TODO: add scenario for when not using Hencky material
    pip.push_back(new ThermoPlasticOps::OpCalculateAdiabaticHeatingRhs<
                  SPACE_DIM, IT, AssemblyDomainEleOp, IS_LARGE_STRAINS>(
        temperature, common_hencky_ptr->getMatHenckyStress(),
        common_plastic_ptr->getPlasticStrainDotPtr(), inelastic_heat_fraction));
 
    pip.push_back(
        new
        typename P::template Assembly<A>::template OpCalculateConstraintsRhs<I>(
            tau, common_plastic_ptr, m_D_ptr));
    pip.push_back(
        new
        typename P::template Assembly<A>::template OpCalculatePlasticFlowRhs<
            DIM, I>(ep, common_plastic_ptr, m_D_ptr));
 
    MoFEMFunctionReturn(0);
  }
 
  template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
  MoFEMErrorCode opThermoPlasticFactoryDomainLhs(
      MoFEM::Interface &m_field, std::string block_name,
      std::string thermal_block_name, std::string thermoelastic_block_name,
      std::string thermoplastic_block_name, Pip &pip, std::string u,
      std::string ep, std::string tau, std::string temperature) {
    MoFEMFunctionBegin;
 
    using namespace HenckyOps;
    using H = HenckyOps::HenckyIntegrators<DomainEleOp>;
 
    using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
        A>::template BiLinearForm<I>;
    using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
    using OpKCauchy = typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
 
    using P = PlasticityIntegrators<DomainEleOp>;
 
    auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
          common_thermoelastic_ptr, common_thermoplastic_ptr] =
        createCommonThermoPlasticOps<DIM, I, DomainEleOp>(
            m_field, block_name, thermal_block_name, thermoelastic_block_name,
            thermoplastic_block_name, pip, u, ep, tau, temperature, scale,
            thermal_conductivity_scaling, heat_capacity_scaling,
            inelastic_heat_fraction_scaling, Sev::inform);
 
    auto m_D_ptr = common_hencky_ptr->matDPtr;
 
    pip.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
        ep, common_plastic_ptr->getPlasticStrainDotPtr()));
    pip.push_back(new OpCalculateScalarFieldValuesDot(
        tau, common_plastic_ptr->getPlasticTauDotPtr()));
    pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
        u, common_plastic_ptr, m_D_ptr, common_thermoplastic_ptr));
 
    if (common_hencky_ptr) {
      pip.push_back(new typename H::template OpHenckyTangent<DIM, I, 0>(
          u, common_hencky_ptr, m_D_ptr));
      pip.push_back(new OpKPiola(u, u, common_hencky_ptr->getMatTangent()));
      pip.push_back(
          new typename P::template Assembly<A>::
              template OpCalculatePlasticInternalForceLhs_LogStrain_dEP<DIM, I>(
                  u, ep, common_plastic_ptr, common_hencky_ptr, m_D_ptr));
    } else {
      pip.push_back(new OpKCauchy(u, u, m_D_ptr));
      pip.push_back(new typename P::template Assembly<A>::
                        template OpCalculatePlasticInternalForceLhs_dEP<DIM, I>(
                            u, ep, common_plastic_ptr, m_D_ptr));
    }
 
    if (common_hencky_ptr) {
      pip.push_back(
          new typename P::template Assembly<A>::
              template OpCalculateConstraintsLhs_LogStrain_dU<DIM, I>(
                  tau, u, common_plastic_ptr, common_hencky_ptr, m_D_ptr));
      pip.push_back(
          new typename P::template Assembly<A>::
              template OpCalculatePlasticFlowLhs_LogStrain_dU<DIM, I>(
                  ep, u, common_plastic_ptr, common_hencky_ptr, m_D_ptr));
    } else {
      pip.push_back(new typename P::template Assembly<A>::
                        template OpCalculateConstraintsLhs_dU<DIM, I>(
                            tau, u, common_plastic_ptr, m_D_ptr));
      pip.push_back(new typename P::template Assembly<A>::
                        template OpCalculatePlasticFlowLhs_dU<DIM, I>(
                            ep, u, common_plastic_ptr, m_D_ptr));
    }
 
    pip.push_back(new typename P::template Assembly<A>::
                      template OpCalculatePlasticFlowLhs_dEP<DIM, I>(
                          ep, ep, common_plastic_ptr, m_D_ptr));
    pip.push_back(new typename P::template Assembly<A>::
                      template OpCalculatePlasticFlowLhs_dTAU<DIM, I>(
                          ep, tau, common_plastic_ptr, m_D_ptr));
    pip.push_back(
        new typename ThermoPlasticOps::OpCalculatePlasticFlowLhs_dT<DIM, I>(
            ep, temperature, common_thermoplastic_ptr));
    pip.push_back(new typename P::template Assembly<A>::
                      template OpCalculateConstraintsLhs_dEP<DIM, I>(
                          tau, ep, common_plastic_ptr, m_D_ptr));
    pip.push_back(
        new typename P::template Assembly<
            A>::template OpCalculateConstraintsLhs_dTAU<I>(tau, tau,
                                                           common_plastic_ptr));
 
    // TODO: add scenario for when not using Hencky material
    pip.push_back(new typename H::template OpCalculateHenckyThermalStressdT<
                  DIM, I, AssemblyDomainEleOp, 0>(
        u, temperature, common_hencky_ptr,
        common_thermoelastic_ptr->getCoeffExpansionPtr()));
 
    auto inelastic_heat_frac_ptr =
        common_thermoplastic_ptr->getInelasticHeatFractionPtr();
    auto inelastic_heat_fraction =
        [inelastic_heat_frac_ptr](const double, const double, const double) {
          return (*inelastic_heat_frac_ptr);
        };
 
    // TODO: add scenario for when not using Hencky material
    pip.push_back(new ThermoPlasticOps::OpCalculateAdiabaticHeatingLhsdT<
                  SPACE_DIM, IT, AssemblyDomainEleOp, IS_LARGE_STRAINS>(
        temperature, temperature, common_hencky_ptr, common_plastic_ptr,
        inelastic_heat_fraction,
        common_thermoelastic_ptr->getCoeffExpansionPtr()));
 
    // TODO: add scenario for when not using Hencky material
    pip.push_back(new ThermoPlasticOps::OpCalculateAdiabaticHeatingLhsdEP<
                  SPACE_DIM, IT, AssemblyDomainEleOp, IS_LARGE_STRAINS>(
        temperature, ep, common_hencky_ptr, common_plastic_ptr,
        inelastic_heat_fraction));
 
    // TODO: add scenario for when not using Hencky material
    pip.push_back(new ThermoPlasticOps::OpCalculateAdiabaticHeatingLhsdU<
                  SPACE_DIM, IT, AssemblyDomainEleOp, IS_LARGE_STRAINS>(
        temperature, u, common_hencky_ptr, common_plastic_ptr,
        inelastic_heat_fraction));
 
    pip.push_back(new ThermoPlasticOps::OpCalculateConstraintsLhs_dT<I>(
        tau, temperature, common_thermoplastic_ptr));
 
    MoFEMFunctionReturn(0);
  }
 
  template <int DIM, IntegrationType I, typename DomainEleOp>
  auto createCommonThermoPlasticOps(
      MoFEM::Interface &m_field, std::string plastic_block_name,
      std::string thermal_block_name, std::string thermoelastic_block_name,
      std::string thermoplastic_block_name,
      boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
      std::string u, std::string ep, std::string tau, std::string temperature,
      double scale, ScalerFunTwoArgs thermal_conductivity_scaling,
      ScalerFunTwoArgs heat_capacity_scaling,
      ScalerFunThreeArgs inelastic_heat_fraction_scaling, Sev sev,
      bool with_rates = true) {
 
    using P = PlasticityIntegrators<DomainEleOp>;
 
    auto common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
    auto common_thermal_ptr =
        boost::make_shared<ThermoElasticOps::BlockedThermalParameters>();
    auto common_thermoelastic_ptr =
        boost::make_shared<ThermoElasticOps::BlockedThermoElasticParameters>();
    auto common_thermoplastic_ptr =
        boost::make_shared<ThermoPlasticOps::ThermoPlasticBlockedParameters>();
 
    constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
    auto make_d_mat = []() {
      return boost::make_shared<MatrixDouble>(size_symm * size_symm, 1);
    };
 
    common_plastic_ptr->mDPtr = boost::make_shared<MatrixDouble>();
    common_plastic_ptr->mGradPtr = boost::make_shared<MatrixDouble>();
    common_plastic_ptr->mStrainPtr = boost::make_shared<MatrixDouble>();
    common_plastic_ptr->mStressPtr = boost::make_shared<MatrixDouble>();
 
    auto m_D_ptr = common_plastic_ptr->mDPtr;
 
    CHK_THROW_MESSAGE(PlasticOps::addMatBlockOps<DIM>(
                          m_field, plastic_block_name, pip, m_D_ptr,
                          common_plastic_ptr->getParamsPtr(), scale, sev),
                      "add mat block plastic ops");
    CHK_THROW_MESSAGE(ThermoElasticOps::addMatThermalBlockOps(
                          m_field, pip, thermal_block_name, common_thermal_ptr,
                          default_heat_conductivity, default_heat_capacity,
                          default_thermal_conductivity_scale,
                          default_heat_capacity_scale, sev,
                          thermal_conductivity_scaling, heat_capacity_scaling),
                      "add mat block thermal ops");
    CHK_THROW_MESSAGE(ThermoElasticOps::addMatThermoElasticBlockOps(
                          m_field, pip, thermoelastic_block_name,
                          common_thermoelastic_ptr, default_coeff_expansion,
                          default_ref_temp, sev),
                      "add mat block thermal ops");
    CHK_THROW_MESSAGE(ThermoPlasticOps::addMatThermoPlasticBlockOps(
                          m_field, pip, thermoplastic_block_name,
                          common_thermoplastic_ptr, common_thermal_ptr, sev,
                          inelastic_heat_fraction_scaling),
                      "add mat block thermoplastic ops");
    auto common_hencky_ptr =
        HenckyOps::commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
            mField, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
 
    common_plastic_ptr->mDPtr = common_hencky_ptr->matDPtr;
 
    pip.push_back(new OpCalculateScalarFieldValues(
        tau, common_plastic_ptr->getPlasticTauPtr()));
    pip.push_back(new OpCalculateTensor2SymmetricFieldValues<DIM>(
        ep, common_plastic_ptr->getPlasticStrainPtr()));
    pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(
        u, common_plastic_ptr->mGradPtr));
 
    pip.push_back(new OpCalculateScalarFieldValues(
        temperature, common_thermoplastic_ptr->getTempPtr()));
    pip.push_back(new OpCalculateHVecVectorField<3, SPACE_DIM>(
        "FLUX", common_thermoplastic_ptr->getHeatFluxPtr()));
 
    common_plastic_ptr->mGradPtr = common_hencky_ptr->matGradPtr;
    common_plastic_ptr->mDPtr = common_hencky_ptr->matDPtr;
    common_hencky_ptr->matLogCPlastic =
        common_plastic_ptr->getPlasticStrainPtr();
    common_plastic_ptr->mStrainPtr = common_hencky_ptr->getMatLogC();
    common_plastic_ptr->mStressPtr = common_hencky_ptr->getMatHenckyStress();
 
    using H = HenckyOps::HenckyIntegrators<DomainEleOp>;
 
    pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
        u, common_hencky_ptr));
    pip.push_back(
        new typename H::template OpCalculateLogC<DIM, I>(u, common_hencky_ptr));
    pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
        u, common_hencky_ptr));
 
    pip.push_back(
        new
        typename H::template OpCalculateHenckyThermoPlasticStress<DIM, I, 0>(
            u, common_thermoplastic_ptr->getTempPtr(), common_hencky_ptr,
            common_thermoelastic_ptr->getCoeffExpansionPtr(),
            common_thermoelastic_ptr->getRefTempPtr()));
    pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I, 0>(
        u, common_hencky_ptr));
 
    pip.push_back(new typename P::template OpCalculatePlasticSurface<DIM, I>(
        u, common_plastic_ptr));
 
    pip.push_back(new typename P::template OpCalculatePlasticHardening<DIM, I>(
        u, common_plastic_ptr, common_thermoplastic_ptr));
 
    return std::make_tuple(common_plastic_ptr, common_hencky_ptr,
                           common_thermal_ptr, common_thermoelastic_ptr,
                           common_thermoplastic_ptr);
  }
 
  friend struct TSPrePostProc;
};
 
MoFEMErrorCode TSPrePostProc::tsPreProc(TS ts) {
  MoFEMFunctionBegin;
 
  if (auto ptr = tsPrePostProc.lock()) {
    MOFEM_LOG("THERMOPLASTICITY", Sev::inform) << "Run step pre proc";
 
    auto &m_field = ptr->ExRawPtr->mField;
    auto simple = m_field.getInterface<Simple>();
 
    ResizeCtx *ctx = nullptr;
    CHKERR TSGetApplicationContext(ts, (void **)&ctx);
    if (!ctx)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "ResizeCtx missing");
 
    if (ctx->mesh_changed == PETSC_TRUE) { // need to be set for new meshsets
      ptr->ExRawPtr->mechanicalBC(comp_mesh_bit, BitRefLevel().set());       // if remeshed
      ptr->ExRawPtr->thermalBC(comp_mesh_bit, BitRefLevel().set());
      // #ifdef ADD_CONTACT
      //       auto pip_mng = m_field.getInterface<PipelineManager>();
      //       pip_mng->getOpDomainLhsPipeline().clear();
      //       pip_mng->getOpDomainRhsPipeline().clear();
      //       pip_mng->getOpBoundaryLhsPipeline().clear();
      //       pip_mng->getOpBoundaryRhsPipeline().clear();
      //       CHKERR ptr->ExRawPtr->OPs();
      // #endif // ADD_CONTACT
    }
 
    // get vector norm
    auto get_norm = [&](auto x) {
      double nrm;
      CHKERR VecNorm(x, NORM_2, &nrm);
      return nrm;
    };
 
    auto save_range = [](moab::Interface &moab, const std::string name,
                         const Range r) {
      MoFEMFunctionBegin;
      auto out_meshset = get_temp_meshset_ptr(moab);
      CHKERR moab.add_entities(*out_meshset, r);
      CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(),
                             1);
      MoFEMFunctionReturn(0);
    };
 
    auto dm = simple->getDM();
    auto d_vector = createDMVector(dm);
    CHKERR DMoFEMMeshToLocalVector(dm, d_vector, INSERT_VALUES,
                                   SCATTER_FORWARD);
    CHKERR VecGhostUpdateBegin(d_vector, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(d_vector, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToGlobalVector(dm, d_vector, INSERT_VALUES,
                                    SCATTER_REVERSE);
 
    Tag th_spatial_coords;
    std::multimap<double, EntityHandle> el_q_map;
    std::vector<EntityHandle> new_connectivity;
    Range flipped_els;
    bool do_refine = false;
 
    auto check_element_quality = [&]() {
      MoFEMFunctionBegin;
 
      // Get Range of all elements not to be flipped in the mesh
      Range virgin_entities;
      CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
          virgin_mesh_bit, BitRefLevel().set(), virgin_entities);
 
      // Calculates the element quality
      CHKERR ptr->ExRawPtr->getElementQuality(el_q_map, flipped_els,
                                              new_connectivity, do_refine,
                                              th_spatial_coords);
 
      MoFEMFunctionReturn(0);
    };
 
    auto solve_equil_sub_problem = [&]() {
      MoFEMFunctionBegin;
      if (el_q_map.size() == 0) {
        MoFEMFunctionReturnHot(0);
      }
 
      auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order; };
 
      auto simple = m_field.getInterface<Simple>();
 
      using UDO = ForcesAndSourcesCore::UserDataOperator;
 
      auto U_ptr = boost::make_shared<MatrixDouble>();
      auto EP_ptr = boost::make_shared<MatrixDouble>();
      auto TAU_ptr = boost::make_shared<VectorDouble>();
      auto T_ptr = boost::make_shared<VectorDouble>();
      // auto FLUX_ptr = boost::make_shared<MatrixDouble>();
      // auto T_grad_ptr = boost::make_shared<MatrixDouble>();
 
      auto post_proc = [&](auto dm) {
        MoFEMFunctionBegin;
 
        auto post_proc_fe = boost::make_shared<PostProcEle>(m_field);
 
        CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
            post_proc_fe->getOpPtrVector(), {H1, L2, HDIV});
 
        using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
        // post_proc_fe->getOpPtrVector().push_back(
        //     new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX",
        //     FLUX_ptr));
        post_proc_fe->getOpPtrVector().push_back(
            new OpCalculateScalarFieldValues("T", T_ptr));
        // post_proc_fe->getOpPtrVector().push_back(
        //     new OpCalculateScalarFieldGradient<SPACE_DIM>("T",
        //     T_grad_ptr));
        post_proc_fe->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValues<SPACE_DIM>("U", U_ptr));
        post_proc_fe->getOpPtrVector().push_back(
            new OpCalculateScalarFieldValues("TAU", TAU_ptr));
        post_proc_fe->getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>("EP",
                                                                  EP_ptr));
        post_proc_fe->getOpPtrVector().push_back(
 
            new OpPPMap(post_proc_fe->getPostProcMesh(),
                        post_proc_fe->getMapGaussPts(),
 
                        {
                            {"TAU", TAU_ptr},
                            {"T", T_ptr},
                        },
 
                        {{"U", U_ptr}},
 
                        {},
 
                        {{"EP", EP_ptr}}
 
                        )
 
        );
 
        CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
        CHKERR post_proc_fe->writeFile("out_equilibrate.h5m");
 
        MoFEMFunctionReturn(0);
      };
 
      auto solve_equilibrate_solution = [&]() {
        MoFEMFunctionBegin;
 
        auto set_domain_rhs = [&](auto &pip) {
          MoFEMFunctionBegin;
 
          CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip,
                                                                {H1, HDIV});
 
          using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
              AT>::template LinearForm<IT>;
          using OpInternalForceCauchy =
              typename B::template OpGradTimesSymTensor<1, SPACE_DIM,
                                                        SPACE_DIM>;
          using OpInternalForcePiola =
              typename B::template OpGradTimesTensor<1, SPACE_DIM, SPACE_DIM>;
 
          using P = PlasticityIntegrators<DomainEleOp>;
 
          auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
                common_thermoelastic_ptr, common_thermoplastic_ptr] =
              ptr->ExRawPtr
                  ->createCommonThermoPlasticOps<SPACE_DIM, IT, DomainEleOp>(
                      m_field, "MAT_PLASTIC", "MAT_THERMAL",
                      "MAT_THERMOELASTIC", "MAT_THERMOPLASTIC", pip, "U", "EP",
                      "TAU", "T", scale, thermal_conductivity_scaling,
                      heat_capacity_scaling, inelastic_heat_fraction_scaling,
                      Sev::inform);
 
          auto m_D_ptr = common_plastic_ptr->mDPtr;
 
          pip.push_back(new OpCalculateScalarFieldValues(
              "T", common_thermoplastic_ptr->getTempPtr()));
          pip.push_back(
              new
              typename P::template OpCalculatePlasticityWithoutRates<SPACE_DIM,
                                                                     IT>(
                  "U", common_plastic_ptr, m_D_ptr, common_thermoplastic_ptr));
 
          // Calculate internal forces
          if (common_hencky_ptr) {
            pip.push_back(new OpInternalForcePiola(
                "U", common_hencky_ptr->getMatFirstPiolaStress()));
          } else {
            pip.push_back(
                new OpInternalForceCauchy("U", common_plastic_ptr->mStressPtr));
          }
 
          MoFEMFunctionReturn(0);
        };
 
        auto set_domain_lhs = [&](auto &pip) {
          MoFEMFunctionBegin;
 
          CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip,
                                                                {H1, L2, HDIV});
 
          using namespace HenckyOps;
 
          using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
              AT>::template BiLinearForm<IT>;
          using OpKPiola =
              typename B::template OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
          using OpKCauchy =
              typename B::template OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM,
                                                       0>;
 
          using P = PlasticityIntegrators<DomainEleOp>;
 
          auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
                common_thermoelastic_ptr, common_thermoplastic_ptr] =
              ptr->ExRawPtr
                  ->createCommonThermoPlasticOps<SPACE_DIM, IT, DomainEleOp>(
                      m_field, "MAT_PLASTIC", "MAT_THERMAL",
                      "MAT_THERMOELASTIC", "MAT_THERMOPLASTIC", pip, "U", "EP",
                      "TAU", "T", scale, thermal_conductivity_scaling,
                      heat_capacity_scaling, inelastic_heat_fraction_scaling,
                      Sev::inform);
 
          auto m_D_ptr = common_plastic_ptr->mDPtr;
 
          pip.push_back(
              new
              typename P::template OpCalculatePlasticityWithoutRates<SPACE_DIM,
                                                                     IT>(
                  "U", common_plastic_ptr, m_D_ptr, common_thermoplastic_ptr));
 
          if (common_hencky_ptr) {
            using H = HenckyOps::HenckyIntegrators<DomainEleOp>;
            pip.push_back(
                new typename H::template OpHenckyTangent<SPACE_DIM, IT, 0>(
                    "U", common_hencky_ptr, m_D_ptr));
            pip.push_back(
                new OpKPiola("U", "U", common_hencky_ptr->getMatTangent()));
            // pip.push_back(
            //     new typename P::template Assembly<AT>::
            //         template
            //         OpCalculatePlasticInternalForceLhs_LogStrain_dEP<
            //             SPACE_DIM, IT>("U", "EP", common_plastic_ptr,
            //                     common_hencky_ptr, m_D_ptr));
          } else {
            pip.push_back(new OpKCauchy("U", "U", m_D_ptr));
            // pip.push_back(
            //     new typename P::template Assembly<AT>::
            //         template
            //         OpCalculatePlasticInternalForceLhs_dEP<SPACE_DIM,
            //         IT>(
            //             "U", "EP", common_plastic_ptr, m_D_ptr));
          }
 
          // TODO: add scenario for when not using Hencky material
          // pip.push_back(new
          // HenckyOps::OpCalculateHenckyThermalStressdT<
          //               SPACE_DIM, IT, AssemblyDomainEleOp>(
          //     "U", "T", common_hencky_ptr,
          //     common_thermal_ptr->getCoeffExpansionPtr()));
 
          MoFEMFunctionReturn(0);
        };
 
        auto create_sub_dm = [&](SmartPetscObj<DM> base_dm) {
          auto dm_sub = createDM(m_field.get_comm(), "DMMOFEM");
          CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "EQUIL_DM");
          CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
          CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
          for (auto f : {"U", "EP", "TAU", "T"}) {
            CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
            CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
          }
          CHKERR DMSetUp(dm_sub);
          return dm_sub;
        };
 
        auto sub_dm = create_sub_dm(simple->getDM());
 
        auto fe_rhs = boost::make_shared<DomainEle>(m_field);
        auto fe_lhs = boost::make_shared<DomainEle>(m_field);
 
        fe_rhs->getRuleHook = vol_rule;
        fe_lhs->getRuleHook = vol_rule;
        CHKERR set_domain_rhs(fe_rhs->getOpPtrVector());
        CHKERR set_domain_lhs(fe_lhs->getOpPtrVector());
 
        auto bc_mng = m_field.getInterface<BcManager>();
 
        CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
            "EQUIL_DM", "U");
 
        auto &bc_map = bc_mng->getBcMapByBlockName();
        for (auto bc : bc_map)
          MOFEM_LOG("PLASTICITY", Sev::verbose) << "Marker " << bc.first;
 
        auto time_scale = boost::make_shared<TimeScale>(
            "", false, [](const double) { return 1; });
        auto def_time_scale = [time_scale](const double time) {
          return (!time_scale->argFileScale) ? time : 1;
        };
        auto def_file_name = [time_scale](const std::string &&name) {
          return (!time_scale->argFileScale) ? name : "";
        };
        auto time_displacement_scaling = boost::make_shared<TimeScale>(
            def_file_name("displacement_bc_scale.txt"), false, def_time_scale);
 
        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>();
 
        PetscReal current_time;
        TSGetTime(ts, &current_time);
        pre_proc_ptr->ts_t = current_time;
 
        // Set BCs by eliminating degrees of freedom
        auto get_bc_hook_rhs = [&]() {
          EssentialPreProc<DisplacementCubitBcData> hook(
              ptr->ExRawPtr->mField, pre_proc_ptr,
              {time_scale, time_displacement_scaling});
          return hook;
        };
        auto get_bc_hook_lhs = [&]() {
          EssentialPreProc<DisplacementCubitBcData> hook(
              ptr->ExRawPtr->mField, pre_proc_ptr,
              {time_scale, time_displacement_scaling});
          return hook;
        };
 
        pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
        pre_proc_ptr->preProcessHook = get_bc_hook_lhs();
 
        auto get_post_proc_hook_rhs = [&]() {
          MoFEMFunctionBegin;
          CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
              ptr->ExRawPtr->mField, post_proc_rhs_ptr, nullptr,
              Sev::verbose)();
          CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
              ptr->ExRawPtr->mField, post_proc_rhs_ptr, 1.)();
          MoFEMFunctionReturn(0);
        };
        auto get_post_proc_hook_lhs = [&]() {
          return EssentialPostProcLhs<DisplacementCubitBcData>(
              ptr->ExRawPtr->mField, post_proc_lhs_ptr, 1.);
        };
        post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
        post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
 
        auto snes_ctx_ptr = getDMSnesCtx(sub_dm);
        snes_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_ptr);
        snes_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_ptr);
        snes_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs_ptr);
        snes_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs_ptr);
 
        auto null_fe = boost::shared_ptr<FEMethod>();
        // CHKERR DMMoFEMKSPSetComputeOperators(
        //     sub_dm, simple->getDomainFEName(), fe_lhs, null_fe,
        //     null_fe);
        // CHKERR DMMoFEMKSPSetComputeRHS(sub_dm,
        // simple->getDomainFEName(),
        //                                fe_rhs, null_fe, null_fe);
        CHKERR DMMoFEMSNESSetFunction(sub_dm, simple->getDomainFEName(), fe_rhs,
                                      null_fe, null_fe);
        CHKERR DMMoFEMSNESSetJacobian(sub_dm, simple->getDomainFEName(), fe_lhs,
                                      null_fe, null_fe);
 
        // auto ksp = MoFEM::createKSP(m_field.get_comm());
        // CHKERR KSPSetDM(ksp, sub_dm);
        // CHKERR KSPSetFromOptions(ksp);
 
        auto D = createDMVector(sub_dm);
 
        auto snes = MoFEM::createSNES(m_field.get_comm());
        CHKERR SNESSetDM(snes, sub_dm);
        CHKERR SNESSetFromOptions(snes);
        CHKERR SNESSetUp(snes);
 
        CHKERR DMoFEMMeshToLocalVector(sub_dm, D, INSERT_VALUES,
                                       SCATTER_FORWARD);
        CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
        // SmartPetscObj<Vec> global_rhs;
        // CHKERR DMCreateGlobalVector_MoFEM(
        //     sub_dm, global_rhs);
        // CHKERR KSPSolve(ksp, PETSC_NULLPTR, D);
        CHKERR SNESSolve(snes, PETSC_NULLPTR, D);
 
        CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR DMoFEMMeshToLocalVector(sub_dm, D, INSERT_VALUES,
                                       SCATTER_REVERSE);
 
        MoFEMFunctionReturn(0);
      };
 
      CHKERR solve_equilibrate_solution();
      CHKERR post_proc(simple->getDM());
 
      MOFEM_LOG("REMESHING", Sev::inform)
          << "Equilibrated solution after mapping";
 
      MoFEMFunctionReturn(0);
    };
 
    if (ctx->mesh_changed == PETSC_FALSE) // protects for rollback
      CHKERR check_element_quality();     // improve element quality
    ctx->mesh_changed = PETSC_FALSE;      // needs reset if true
    if (el_q_map.size() > 0 || do_refine) {
      // Now you can safely update the context fields
      ctx->ptr = ptr;
      ctx->el_q_map = el_q_map;
      ctx->flipped_els = flipped_els;
      ctx->new_connectivity = new_connectivity;
      ctx->th_spatial_coords = th_spatial_coords;
      ctx->mesh_changed = PETSC_TRUE;
 
      // Then trigger the resize process
      // CHKERR TSResize(ts);
    }
 
    // Need barriers, some functions in TS solver need are called
    // collectively and requite the same state of variables
    CHKERR PetscBarrier((PetscObject)ts);
 
    MOFEM_LOG_CHANNEL("SYNC");
    MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "REMESHING") << "PreProc done";
    MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode TSPrePostProc::tsPostProc(TS ts) {
  if (auto ptr = tsPrePostProc.lock()) {
    auto &m_field = ptr->ExRawPtr->mField;
    MOFEM_LOG_CHANNEL("SYNC");
    MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "THERMOPLASTICITY")
        << "PostProc done";
    MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
  }
  return 0;
}
 
MoFEMErrorCode TSPrePostProc::tsSetUp(TS ts) {
 
  MoFEMFunctionBegin;
 
  CHKERR TSSetPreStep(ts, tsPreProc);
  CHKERR TSSetPostStep(ts, tsPostProc);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ThermoPlasticOps::addMatThermoPlasticBlockOps(
    MoFEM::Interface &m_field,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
    std::string thermoplastic_block_name,
    boost::shared_ptr<ThermoPlasticOps::ThermoPlasticBlockedParameters>
        blockedParamsPtr,
    boost::shared_ptr<ThermoElasticOps::BlockedThermalParameters>
        &blockedThermalParamsPtr,
    Sev sev, ScalerFunThreeArgs inelastic_heat_fraction_scaling_func) {
  MoFEMFunctionBegin;
 
  struct OpMatThermoPlasticBlocks : public DomainEleOp {
    OpMatThermoPlasticBlocks(
        boost::shared_ptr<double> omega_0_ptr,
        boost::shared_ptr<double> omega_h_ptr,
        boost::shared_ptr<double> inelastic_heat_fraction_ptr,
        boost::shared_ptr<double> temp_0_ptr,
        boost::shared_ptr<double> heat_conductivity,
        boost::shared_ptr<double> heat_capacity,
        boost::shared_ptr<double> heat_conductivity_scaling,
        boost::shared_ptr<double> heat_capacity_scaling,
        ScalerFunThreeArgs inelastic_heat_fraction_scaling,
        MoFEM::Interface &m_field, Sev sev,
        std::vector<const CubitMeshSets *> meshset_vec_ptr)
        : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), omega0Ptr(omega_0_ptr),
          omegaHPtr(omega_h_ptr),
          inelasticHeatFractionPtr(inelastic_heat_fraction_ptr),
          temp0Ptr(temp_0_ptr), heatConductivityPtr(heat_conductivity),
          heatCapacityPtr(heat_capacity),
          heatConductivityScalingPtr(heat_conductivity_scaling),
          heatCapacityScalingPtr(heat_capacity_scaling),
          inelasticHeatFractionScaling(inelastic_heat_fraction_scaling) {
      CHK_THROW_MESSAGE(
          extractThermoPlasticBlockData(m_field, meshset_vec_ptr, sev),
          "Can not get data from block");
    }
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
 
      auto scale_param = [this](auto parameter, double scaling) {
        return parameter * scaling;
      };
 
      for (auto &b : blockData) {
 
        if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
          *omega0Ptr = b.omega_0;
          *omegaHPtr = b.omega_h;
          *inelasticHeatFractionPtr = scale_param(
              b.inelastic_heat_fraction,
              inelasticHeatFractionScaling(
                  young_modulus,
                  *heatConductivityPtr / (*heatConductivityScalingPtr),
                  *heatCapacityPtr / (*heatCapacityScalingPtr)));
          *temp0Ptr = b.temp_0;
          MoFEMFunctionReturnHot(0);
        }
      }
 
      *omega0Ptr = omega_0;
      *omegaHPtr = omega_h;
      *inelasticHeatFractionPtr =
          scale_param(inelastic_heat_fraction,
                      inelasticHeatFractionScaling(
                          young_modulus,
                          *heatConductivityPtr / (*heatConductivityScalingPtr),
                          *heatCapacityPtr / (*heatCapacityScalingPtr)));
      *temp0Ptr = temp_0;
 
      MoFEMFunctionReturn(0);
    }
 
  private:
    ScalerFunThreeArgs inelasticHeatFractionScaling;
    boost::shared_ptr<double> heatConductivityPtr;
    boost::shared_ptr<double> heatCapacityPtr;
    boost::shared_ptr<double> heatConductivityScalingPtr;
    boost::shared_ptr<double> heatCapacityScalingPtr;
    struct BlockData {
      double omega_0;
      double omega_h;
      double inelastic_heat_fraction;
      double temp_0;
      Range blockEnts;
    };
 
    std::vector<BlockData> blockData;
 
    MoFEMErrorCode extractThermoPlasticBlockData(
        MoFEM::Interface &m_field,
        std::vector<const CubitMeshSets *> meshset_vec_ptr, Sev sev) {
      MoFEMFunctionBegin;
 
      for (auto m : meshset_vec_ptr) {
        MOFEM_TAG_AND_LOG("WORLD", sev, "Mat ThermoPlastic Block") << *m;
        std::vector<double> block_data;
        CHKERR m->getAttributes(block_data);
        if (block_data.size() < 3) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Expected that block has four attributes");
        }
        auto get_block_ents = [&]() {
          Range ents;
          CHKERR
          m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
          return ents;
        };
 
        MOFEM_TAG_AND_LOG("WORLD", sev, "Mat ThermoPlastic Block")
            << m->getName() << ": omega_0 = " << block_data[0]
            << " omega_h = " << block_data[1]
            << " inelastic_heat_fraction = " << block_data[2] << " temp_0 "
            << block_data[3];
 
        blockData.push_back({block_data[0], block_data[1], block_data[2],
                             block_data[3], get_block_ents()});
      }
      MOFEM_LOG_CHANNEL("WORLD");
      MoFEMFunctionReturn(0);
    }
 
    boost::shared_ptr<double> omega0Ptr;
    boost::shared_ptr<double> omegaHPtr;
    boost::shared_ptr<double> inelasticHeatFractionPtr;
    boost::shared_ptr<double> temp0Ptr;
  };
 
  pipeline.push_back(new OpMatThermoPlasticBlocks(
      blockedParamsPtr->getOmega0Ptr(), blockedParamsPtr->getOmegaHPtr(),
      blockedParamsPtr->getInelasticHeatFractionPtr(),
      blockedParamsPtr->getTemp0Ptr(),
      blockedThermalParamsPtr->getHeatConductivityPtr(),
      blockedThermalParamsPtr->getHeatCapacityPtr(),
      blockedThermalParamsPtr->getHeatConductivityScalingPtr(),
      blockedThermalParamsPtr->getHeatCapacityScalingPtr(),
      inelastic_heat_fraction_scaling_func, m_field, sev,
 
      // Get blockset using regular expression
      m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
          (boost::format("%s(.*)") % thermoplastic_block_name).str()
 
              ))
 
          ));
 
  MoFEMFunctionReturn(0);
}
 
//! [Get element quality]
MoFEMErrorCode
Example::getElementQuality(std::multimap<double, EntityHandle> &el_q_map,
                           Range &flipped_els,
                           std::vector<EntityHandle> &new_connectivity,
                           bool &do_refine, Tag &th_spatial_coords) {
  MoFEMFunctionBegin;
 
  Range verts;
  CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, verts);
  std::vector<double> coords(verts.size() * 3);
  CHKERR mField.get_moab().get_coords(verts, coords.data());
  auto t_x = getFTensor1FromPtr<3>(coords.data());
 
  auto save_tag = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
    MoFEMFunctionBeginHot;
    FTENSOR_INDEX(3, i)
    auto field_data = ent_ptr->getEntFieldData();
    FTensor::Tensor1<double, 3> t_u;
    if (SPACE_DIM == 2) {
      t_u = {field_data[0], field_data[1], 0.0};
    } else {
      t_u = {field_data[0], field_data[1], field_data[2]};
    }
 
    t_x(i) += t_u(i);
    ++t_x;
    MoFEMFunctionReturnHot(0);
  };
 
  mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(save_tag, "U",
                                                          &verts);
  double def_coord[3] = {0.0, 0.0, 0.0};
  CHKERR mField.get_moab().tag_get_handle(
      "SpatialCoord", 3, MB_TYPE_DOUBLE, th_spatial_coords,
      MB_TAG_DENSE | MB_TAG_CREAT, def_coord);
  CHKERR mField.get_moab().tag_set_data(th_spatial_coords, verts,
                                        coords.data());
 
  // get a map which has the element quality for each tag
 
  Range all_els;
  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, all_els,
                                                     true);
  CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
      virgin_mesh_bit, BitRefLevel().set(), all_els);
 
  std::multimap<double, EntityHandle> candidate_el_q_map;
  double spatial_coords[9], material_coords[9];
  const EntityHandle *conn;
  int num_nodes;
  Range::iterator nit = all_els.begin();
  for (int gg = 0; nit != all_els.end(); nit++, gg++) {
    CHKERR mField.get_moab().get_connectivity(*nit, conn, num_nodes, true);
 
    CHKERR mField.get_moab().get_coords(conn, num_nodes, material_coords);
    CHKERR mField.get_moab().tag_get_data(th_spatial_coords, conn, num_nodes,
                                          spatial_coords);
 
    double q = triangleAreaLengthQuality(spatial_coords);
    if (!std::isnormal(q))
      SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
              "Calculated quality of element is not "
              "as expected: %f",
              q);
 
    if (q < qual_tol && q > 0.)
      candidate_el_q_map.insert(std::pair<double, EntityHandle>(q, *nit));
  }
 
  double min_q = 1;
  CHKERR minTriangleQuality(mField, all_els, min_q, th_spatial_coords);
  MOFEM_LOG("REMESHING", Sev::inform)
      << "Old minimum element quality: " << min_q;
  // Get the first element using begin()
  auto pair = candidate_el_q_map.begin();
  MOFEM_LOG("REMESHING", Sev::inform)
      << "New minimum element quality: " << pair->first;
 
  for (auto pair = candidate_el_q_map.begin(); pair != candidate_el_q_map.end();
       ++pair) {
    double quality = pair->first;
    Range element;
    element.insert(pair->second);
    if (!flipped_els.contains(element)) {
      // Get edges of element
      Range edges;
      CHKERR mField.get_moab().get_adjacencies(element, 1, false, edges);
 
      // Get the longest edge
      EntityHandle longest_edge;
      double longest_edge_length = 0;
      std::vector<std::pair<double, EntityHandle>> edge_lengths;
      edge_lengths.reserve(edges.size());
      for (auto edge : edges) {
        edge_lengths.emplace_back(
            mField.getInterface<Tools>()->getEdgeLength(edge), edge);
      }
      if (!edge_lengths.empty()) {
        const auto it = std::max_element(
            edge_lengths.begin(), edge_lengths.end(),
            [](const auto &a, const auto &b) { return a.first < b.first; });
        longest_edge_length = it->first;
        longest_edge = it->second;
      } else {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "Unable to calculate edge lengths to find longest edge.");
      }
 
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "Edge flip longest edge length: " << longest_edge_length
          << " for edge: " << longest_edge;
 
      auto get_skin = [&]() {
        Range body_ents;
        CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                           body_ents);
        Skinner skin(&mField.get_moab());
        Range skin_ents;
        CHKERR skin.find_skin(0, body_ents, false, skin_ents);
        return skin_ents;
      };
 
      auto filter_true_skin = [&](auto skin) {
        Range boundary_ents;
        ParallelComm *pcomm =
            ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
        CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                     PSTATUS_NOT, -1, &boundary_ents);
        return boundary_ents;
      };
 
      Range boundary_ents = get_skin();
 
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "Boundary entities: " << boundary_ents;
 
      // Get neighbouring element with the longest edge
      Range flip_candidate_els;
      CHKERR mField.get_moab().get_adjacencies(&longest_edge, 1, SPACE_DIM,
                                               false, flip_candidate_els);
      CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
          virgin_mesh_bit, BitRefLevel().set(), flip_candidate_els);
 
      Range neighbouring_el = subtract(flip_candidate_els, element);
      CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
          virgin_mesh_bit, BitRefLevel().set(), neighbouring_el);
 
      if (boundary_ents.contains(Range(longest_edge, longest_edge))) {
        continue;
      }
 
      if (flipped_els.contains(neighbouring_el))
        continue; // Already flipped
 
#ifndef NDEBUG
      if (neighbouring_el.size() != 1) {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "Should be 1 neighbouring element to bad element for edge "
                "flip. Instead, there are %zu",
                neighbouring_el.size());
      }
#endif
 
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "flip_candidate_els: " << flip_candidate_els;
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "Neighbouring element: " << neighbouring_el;
 
      // Check the quality of the neighbouring element
      std::vector<EntityHandle> neighbouring_nodes;
      CHKERR mField.get_moab().get_connectivity(&neighbouring_el.front(), 1,
                                                neighbouring_nodes, true);
 
      std::vector<EntityHandle> element_nodes;
      CHKERR mField.get_moab().get_connectivity(&element.front(), 1,
                                                element_nodes, true);
 
      CHKERR mField.get_moab().tag_get_data(
          th_spatial_coords, &neighbouring_nodes.front(), 3, spatial_coords);
 
      double neighbour_qual = triangleAreaLengthQuality(spatial_coords);
      if (!std::isnormal(neighbour_qual))
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "Calculated quality of neighbouring element is not as "
                "expected: %f",
                neighbour_qual);
 
      // Get the nodes of flip_candidate_els as well as the nodes of
      // longest_edge
 
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "Element nodes before swap: "
          << get_string_from_vector(element_nodes);
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "Neighbouring nodes before swap: "
          << get_string_from_vector(neighbouring_nodes);
 
      CHKERR save_range(mField.get_moab(), "bad_element.vtk", element);
      CHKERR save_range(mField.get_moab(), "bad_element_neighbour.vtk",
                        neighbouring_el);
 
      // Get new canonical ordering of new nodes and new neighbouring
      // nodes
      std::vector<EntityHandle> reversed_neighbouring_nodes =
          neighbouring_nodes;
      std::reverse(reversed_neighbouring_nodes.begin(),
                   reversed_neighbouring_nodes.end());
 
      int num_matches = 0;
      std::vector<bool> mismatch_mask(element_nodes.size());
      int loop_counter = 0; // To prevent infinite loop
      while (num_matches != 2) {
        // Permute first element to the end
        std::rotate(reversed_neighbouring_nodes.begin(),
                    reversed_neighbouring_nodes.begin() + 1,
                    reversed_neighbouring_nodes.end());
        // Create a boolean mask
        std::transform(element_nodes.begin(), element_nodes.end(),
                       reversed_neighbouring_nodes.begin(),
                       mismatch_mask.begin(), std::equal_to<EntityHandle>());
        // Check if common edge is found
        num_matches =
            std::count(mismatch_mask.begin(), mismatch_mask.end(), true);
 
        ++loop_counter;
        if (loop_counter > 3) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW,
                  "Not found matching nodes for edge flipping");
        }
      }
 
      // Get matching nodes
      std::vector<EntityHandle> matched_elements(element_nodes.size());
      std::transform(element_nodes.begin(), element_nodes.end(),
                     mismatch_mask.begin(), matched_elements.begin(),
                     [](EntityHandle el, bool match) {
                       return match ? el : -1; // Or some other "null" value
                     });
 
      // Remove zero (or "null") elements
      matched_elements.erase(
          std::remove(matched_elements.begin(), matched_elements.end(), -1),
          matched_elements.end());
 
      // Get mismatching nodes
      std::vector<EntityHandle> mismatched_elements(element_nodes.size()),
          neighbouring_mismatched_elements(neighbouring_nodes.size());
      std::transform(element_nodes.begin(), element_nodes.end(),
                     mismatch_mask.begin(), mismatched_elements.begin(),
                     [](EntityHandle el, bool match) {
                       return match ? -1 : el; // Or some other "null" value
                     });
      std::transform(reversed_neighbouring_nodes.begin(),
                     reversed_neighbouring_nodes.end(), mismatch_mask.begin(),
                     neighbouring_mismatched_elements.begin(),
                     [](EntityHandle el, bool match) {
                       return match ? -1 : el; // Or some other "null" value
                     });
 
      // Remove zero (or "null") elements
      mismatched_elements.erase(std::remove(mismatched_elements.begin(),
                                            mismatched_elements.end(), -1),
                                mismatched_elements.end());
      neighbouring_mismatched_elements.erase(
          std::remove(neighbouring_mismatched_elements.begin(),
                      neighbouring_mismatched_elements.end(), -1),
          neighbouring_mismatched_elements.end());
 
      mismatched_elements.insert(mismatched_elements.end(),
                                 neighbouring_mismatched_elements.begin(),
                                 neighbouring_mismatched_elements.end());
 
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "Reversed neighbouring nodes: "
          << get_string_from_vector(reversed_neighbouring_nodes);
 
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "mismatch mask: " << get_string_from_vector(mismatch_mask);
 
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "Old nodes are: " << get_string_from_vector(matched_elements);
 
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "New nodes are: " << get_string_from_vector(mismatched_elements);
 
      auto replace_correct_nodes = [](std::vector<EntityHandle> &ABC,
                                      std::vector<EntityHandle> &DBA,
                                      const std::vector<EntityHandle> &AB,
                                      const std::vector<EntityHandle> &CD) {
        MoFEMFunctionBegin;
        std::vector<std::vector<EntityHandle>> results;
        // Assume AB.size() == 2, CD.size() == 2 for tris
        for (int i = 0; i < 2; ++i) {   // i: 0 for A, 1 for B
          for (int j = 0; j < 2; ++j) { // j: 0 for C, 1 for D
            // Only try to replace if CD[j] is not already in ABC
            if (std::find(ABC.begin(), ABC.end(), CD[j]) == ABC.end()) {
              std::vector<EntityHandle> tmp = ABC;
              // Find AB[i] in ABC and replace with CD[j]
              auto it = std::find(tmp.begin(), tmp.end(), AB[i]);
              if (it != tmp.end()) {
                *it = CD[j];
                results.push_back(tmp);
              }
            }
          }
        }
 
        if (results.size() != 2) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW,
                  "Failed to find two valid vertex replacements for edge "
                  "flip");
        }
 
        ABC = results[0];
        DBA = results[1];
 
        MoFEMFunctionReturn(0);
      };
 
      CHKERR replace_correct_nodes(element_nodes, neighbouring_nodes,
                                   matched_elements, mismatched_elements);
 
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "Element nodes after swap: "
          << get_string_from_vector(element_nodes);
      MOFEM_LOG("REMESHING", Sev::verbose)
          << "Neighbouring nodes after swap: "
          << get_string_from_vector(neighbouring_nodes);
 
      // Calculate the quality of the new elements
      CHKERR mField.get_moab().tag_get_data(
          th_spatial_coords, &element_nodes.front(), 3, spatial_coords);
 
      double new_qual = triangleAreaLengthQuality(spatial_coords);
      if (new_qual < 0.0 || !std::isfinite(new_qual))
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "Calculated quality of new element is not as expected: %f",
                new_qual);
 
#ifndef NDEBUG
      auto check_normal_direction = [&](double qual_val) {
        MoFEMFunctionBegin;
        FTensor::Index<'i', 3> i;
        auto t_correct_normal = FTensor::Tensor1<double, 3>(0.0, 0.0, 1.0);
        double new_area = 0;
        std::array<double, 3> new_normal;
        CHKERR getTriangleAreaAndNormal(spatial_coords, new_area, new_normal);
        auto t_new_normal = FTensor::Tensor1<double, 3>(
            new_normal[0], new_normal[1], new_normal[2]);
        auto t_diff = FTensor::Tensor1<double, 3>();
        t_diff(i) = t_new_normal(i) - t_correct_normal(i);
        if (qual_val > 1e-6 && t_diff(i) * t_diff(i) > 1e-6) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW,
                  "Direction of element to be created is wrong orientation");
        }
        MoFEMFunctionReturn(0);
      };
 
      CHKERR check_normal_direction(new_qual);
#endif
 
      CHKERR mField.get_moab().tag_get_data(
          th_spatial_coords, &neighbouring_nodes.front(), 3, spatial_coords);
 
      double new_neighbour_qual = triangleAreaLengthQuality(spatial_coords);
      if (new_neighbour_qual < 0.0 || !std::isfinite(new_neighbour_qual))
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "Calculated quality of new neighbouring element is not "
                "as expected: %f",
                new_neighbour_qual);
 
#ifndef NDEBUG
      CHKERR check_normal_direction(new_neighbour_qual);
#endif
 
      // If the minimum element quality has improved, do the flip
      if (std::min(new_qual, new_neighbour_qual) >
          (1. + qual_thresh) * std::min(quality, neighbour_qual)) {
        MOFEM_LOG("REMESHING", Sev::inform)
            << "Element quality improved from " << quality << " and "
            << neighbour_qual << " to " << new_qual << " and "
            << new_neighbour_qual << " for elements" << element << " and "
            << neighbouring_el;
 
        // then push to creation "pipeline".
        flipped_els.merge(flip_candidate_els);
        el_q_map.insert(
            std::pair<double, EntityHandle>(pair->first, pair->second));
        new_connectivity.insert(new_connectivity.end(), element_nodes.begin(),
                                element_nodes.end());
        new_connectivity.insert(new_connectivity.end(),
                                neighbouring_nodes.begin(),
                                neighbouring_nodes.end());
      }
    }
  }
 
  if (el_q_map.size() > 0) {
    MOFEM_LOG("REMESHING", Sev::verbose) << "Flipped elements: " << flipped_els;
    MOFEM_LOG("REMESHING", Sev::verbose)
        << "New connectivity: " << get_string_from_vector(new_connectivity);
  }
 
  Range edges;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
      virgin_mesh_bit, BitRefLevel().set(), SPACE_DIM - 1, edges);
  CHKERR mField.get_moab().tag_get_handle(
      "SpatialCoord", 3, MB_TYPE_DOUBLE, th_spatial_coords,
      MB_TAG_DENSE | MB_TAG_CREAT, def_coord);
 
  Range prev_ref_ents;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      BitRefLevel().set(REFINED_EDGES_BIT), BitRefLevel().set(), MBEDGE,
      prev_ref_ents);
 
  for (auto e : edges) {
    const EntityHandle *conn;
    int num_nodes;
    CHKERR mField.get_moab().get_connectivity(e, conn, num_nodes, true);
    std::array<double, 6> ref_coords, spatial_coords;
    CHKERR mField.get_moab().get_coords(conn, num_nodes, ref_coords.data());
    CHKERR mField.get_moab().tag_get_data(th_spatial_coords, conn, num_nodes,
                                          spatial_coords.data());
    auto get_length = [](auto &a) {
      FTENSOR_INDEX(SPACE_DIM, i);
      FTensor::Tensor1<double, 3> p0{a[0], a[1], a[2]};
      FTensor::Tensor1<double, 3> p1{a[3], a[4], a[5]};
      p1(i) = p1(i) - p0(i);
      return p1.l2();
    };
    auto ref_edge_length = get_length(ref_coords);
    auto spatial_edge_length = get_length(spatial_coords);
    auto change = spatial_edge_length / ref_edge_length;
    if ((change >= 1. + edge_growth_thresh) && (num_refinement_levels > 0) &&
        !prev_ref_ents.contains(Range(e, e))) {
      do_refine = true;
    }
  }
 
  MoFEMFunctionReturn(0);
};
//! [Get element quality]
 
//! [Edge flips]
MoFEMErrorCode Example::edgeFlips(BitRefLevel parent_bit,
                                  BitRefLevel child_bit) {
  MoFEMFunctionBegin;
 
  moab::Interface &moab = mField.get_moab();
 
  auto make_edge_flip = [&](auto edge, auto adj_faces, Range &new_tris) {
    MoFEMFunctionBegin;
 
    auto get_conn = [&](EntityHandle e, EntityHandle *conn_cpy) {
      MoFEMFunctionBegin;
      const EntityHandle *conn;
      int num_nodes;
      CHKERR moab.get_connectivity(e, conn, num_nodes, true);
      std::copy(conn, conn + num_nodes, conn_cpy);
      MoFEMFunctionReturn(0);
    };
 
    auto get_tri_normals = [&](auto &conn) {
      std::array<double, 18> coords;
      CHKERR moab.get_coords(conn.data(), 6, coords.data());
      std::array<FTensor::Tensor1<double, 3>, 2> tri_normals;
      for (int t = 0; t != 2; ++t) {
        CHKERR Tools::getTriNormal(&coords[9 * t], &tri_normals[t](0));
      }
      return tri_normals;
    };
 
    auto test_flip = [&](auto &&t_normals) {
      FTENSOR_INDEX(3, i);
      if (t_normals[0](i) * t_normals[1](i) <
          std::numeric_limits<float>::epsilon())
        return false;
      return true;
    };
 
    std::array<EntityHandle, 6> adj_conn;
    CHKERR get_conn(adj_faces[0], &adj_conn[0]);
    CHKERR get_conn(adj_faces[1], &adj_conn[3]);
    std::array<EntityHandle, 2> edge_conn;
    CHKERR get_conn(edge, edge_conn.data());
    std::array<EntityHandle, 2> new_edge_conn;
 
    int j = 1;
    for (int i = 0; i != 6; ++i) {
      if (adj_conn[i] != edge_conn[0] && adj_conn[i] != edge_conn[1]) {
        new_edge_conn[j] = adj_conn[i];
        --j;
      }
    }
 
    auto &new_conn = adj_conn; //< just alias this
    for (int t = 0; t != 2; ++t) {
      for (int i = 0; i != 3; ++i) {
        if (
 
            (adj_conn[3 * t + i % 3] == edge_conn[0] &&
             adj_conn[3 * t + (i + 1) % 3] == edge_conn[1])
 
            ||
 
            (adj_conn[3 * t + i % 3] == edge_conn[1] &&
             adj_conn[3 * t + (i + 1) % 3] == edge_conn[0])
 
        ) {
          new_conn[3 * t + (i + 1) % 3] = new_edge_conn[t];
          break;
        }
      }
    }
 
    if (test_flip(get_tri_normals(new_conn))) {
      for (int t = 0; t != 2; ++t) {
        Range rtri;
        CHKERR moab.get_adjacencies(&new_conn[3 * t], SPACE_DIM + 1, SPACE_DIM,
                                    false, rtri);
        if (!rtri.size()) {
          EntityHandle tri;
          CHKERR moab.create_element(MBTRI, &new_conn[3 * t], SPACE_DIM + 1,
                                     tri);
          new_tris.insert(tri);
        } else {
#ifndef NDEBUG
          if (rtri.size() != 1) {
            MOFEM_LOG("SELF", Sev::error)
                << "Multiple tries created during edge flip for edge " << edge
                << " adjacent faces " << std::endl
                << rtri;
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "Multiple tries created during edge flip");
          }
#endif // NDEBUG
          new_tris.merge(rtri);
        }
      }
 
      Range new_edges;
      CHKERR moab.get_adjacencies(new_tris, SPACE_DIM - 1, true, new_edges,
                                  moab::Interface::UNION);
    } else {
 
      MOFEM_LOG_CHANNEL("SELF");
      MOFEM_LOG("SELF", Sev::warning) << "Edge flip rejected for edge " << edge
                                      << " adjacent faces " << adj_faces;
    }
 
    MoFEMFunctionReturn(0);
  };
 
  Range tris;
  CHKERR moab.get_entities_by_dimension(0, SPACE_DIM, tris);
  CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
      parent_bit, BitRefLevel().set(), tris);
  Skinner skin(&moab);
  Range skin_edges;
  CHKERR skin.find_skin(0, tris, false, skin_edges);
 
  Range edges;
  CHKERR moab.get_entities_by_dimension(0, SPACE_DIM - 1, edges);
  edges = subtract(edges, skin_edges);
  CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
      parent_bit, BitRefLevel().set(), edges);
 
  Range new_tris, flipped_tris, forbidden_tris;
  int flip_count = 0;
  for (auto edge : edges) {
    Range adjacent_tris;
    CHKERR moab.get_adjacencies(&edge, 1, SPACE_DIM, true, adjacent_tris);
 
    adjacent_tris = intersect(adjacent_tris, tris);
    adjacent_tris = subtract(adjacent_tris, forbidden_tris);
    if (adjacent_tris.size() == 2) {
 
#ifndef NDEBUG
      int side_number0, sense0, offset0;
      CHKERR mField.get_moab().side_number(adjacent_tris[0], edge, side_number0,
                                           sense0, offset0);
      int side_number1, sense1, offset1;
      CHKERR mField.get_moab().side_number(adjacent_tris[1], edge, side_number1,
                                           sense1, offset1);
      if (sense0 * sense1 > 0)
        SETERRQ(
            PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Cannot flip edge with same orientation in both adjacent faces");
#endif // NDEBUG
 
      Range new_flipped_tris;
      CHKERR make_edge_flip(edge, adjacent_tris, new_flipped_tris);
      if (new_flipped_tris.size()) {
        flipped_tris.merge(adjacent_tris);
        forbidden_tris.merge(adjacent_tris);
        new_tris.merge(new_flipped_tris);
 
#ifndef NDEBUG
        CHKERR save_range(moab,
                          "flipped_tris_" + std::to_string(flip_count) + ".vtk",
                          adjacent_tris);
        CHKERR save_range(
            moab, "new_flipped_tris_" + std::to_string(flip_count) + ".vtk",
            new_flipped_tris);
 
#endif // NDEBUG
 
        ++flip_count;
      }
    }
  }
 
  Range all_tris;
  CHKERR moab.get_entities_by_dimension(0, SPACE_DIM, all_tris);
  Range not_flipped_tris = subtract(all_tris, flipped_tris);
 
  MOFEM_LOG("SELF", Sev::noisy)
      << "Flipped " << flip_count << " edges with two adjacent faces.";
  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevel(not_flipped_tris,
                                                              child_bit);
  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevel(new_tris,
                                                              child_bit);
  CHKERR mField.getInterface<BitRefManager>()->writeBitLevel(
      child_bit, BitRefLevel().set(), "edge_flips_before_refinement.vtk", "VTK",
      "");
 
  MoFEMFunctionReturn(0);
}
//! [Edge flips]
 
//! [Do Edge Flips]
MoFEMErrorCode
Example::doEdgeFlips(std::multimap<double, EntityHandle> &el_q_map,
                     Range &flipped_els, Tag &th_spatial_coords,
                     std::vector<EntityHandle> &new_connectivity) {
  MoFEMFunctionBegin;
 
  ReadUtilIface *read_util;
  CHKERR mField.get_moab().query_interface(read_util);
 
  const int num_ele = el_q_map.size() * 2; // each flip creates 2 elements
  int num_nod_per_ele;
  EntityType ent_type;
 
  if (SPACE_DIM == 2) {
    num_nod_per_ele = 3;
    ent_type = MBTRI;
  } else {
    num_nod_per_ele = 4;
    ent_type = MBTET;
  }
 
  Range new_elements;
  auto new_conn = new_connectivity.begin();
  for (auto e = 0; e != num_ele; ++e) {
    EntityHandle conn[num_nod_per_ele];
    for (int n = 0; n < num_nod_per_ele; ++n) {
      conn[n] = *new_conn;
      ++new_conn;
    }
    EntityHandle new_ele;
    Range adj_ele;
    CHKERR mField.get_moab().get_adjacencies(conn, num_nod_per_ele, SPACE_DIM,
                                             false, adj_ele);
    if (adj_ele.size()) {
      if (adj_ele.size() != 1) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW,
                "Element duplication");
      } else {
        new_ele = adj_ele.front();
      }
    } else {
      CHKERR mField.get_moab().create_element(ent_type, conn, num_nod_per_ele,
                                              new_ele);
    }
    new_elements.insert(new_ele);
  }
 
  MOFEM_LOG("REMESHING", Sev::verbose)
      << "New elements from edge flipping: " << new_elements;
 
  auto reset_flip_bit = [](EntityHandle ent, BitRefLevel &bit) {
    bit.set(STORAGE_BIT, true);
    bit.set(FLIPPED_BIT, false);
  };
  CHKERR mField.getInterface<BitRefManager>()->lambdaBitRefLevel(
      reset_flip_bit);
 
  auto get_adj = [&](auto ents) {
    Range adj;
    for (auto d = 0; d != SPACE_DIM; ++d) {
      CHK_THROW_MESSAGE(
          mField.get_moab().get_adjacencies(ents, d, true, adj,
                                            moab::Interface::UNION),
          "Getting adjacencies of dimension " + std::to_string(d) + " failed");
    }
    return adj;
  };
 
  Range non_flipped_range;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
      virgin_mesh_bit, BitRefLevel().set(), SPACE_DIM, non_flipped_range);
 
  non_flipped_range = subtract(non_flipped_range, flipped_els);
  auto flip_bit_ents = new_elements;
  flip_bit_ents.merge(non_flipped_range);
  auto adj = get_adj(new_elements);
  // flip_bit_ents.merge(get_adj(new_elements));
 
  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevel(
      flip_bit_ents, flipped_bit, false);
 
  CHKERR mField.getInterface<BitRefManager>()->writeBitLevelByDim(
      BitRefLevel().set(FLIPPED_BIT), BitRefLevel().set(), 2,
      "new_elements_after_edge_flips.vtk", "VTK", "");
 
  MoFEMFunctionReturn(0);
};
//! [Do Edge Flips]
 
//! [Refine edges]
MoFEMErrorCode Example::doEdgeSplits(bool &refined, bool add_ents) {
  MoFEMFunctionBegin;
 
  Tag th_spatial_coords;
  double def_coord[3] = {0.0, 0.0, 0.0};
 
  auto refined_mesh = [&](auto level) {
    MoFEMFunctionBeginHot;
    auto bit = BitRefLevel().set(FLIPPED_BIT + level);
    Range edges;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
        bit, BitRefLevel().set(), SPACE_DIM - 1, edges);
    CHKERR mField.get_moab().tag_get_handle(
        "SpatialCoord", 3, MB_TYPE_DOUBLE, th_spatial_coords,
        MB_TAG_DENSE | MB_TAG_CREAT, def_coord);
 
    Range refine_edges;
    if (add_ents)
      for (auto e : edges) {
        const EntityHandle *conn;
        int num_nodes;
        CHKERR mField.get_moab().get_connectivity(e, conn, num_nodes, true);
        std::array<double, 6> ref_coords, spatial_coords;
        CHKERR mField.get_moab().get_coords(conn, num_nodes, ref_coords.data());
        CHKERR mField.get_moab().tag_get_data(th_spatial_coords, conn,
                                              num_nodes, spatial_coords.data());
        auto get_length = [](auto &a) {
          FTENSOR_INDEX(SPACE_DIM, i);
          FTensor::Tensor1<double, 3> p0{a[0], a[1], a[2]};
          FTensor::Tensor1<double, 3> p1{a[3], a[4], a[5]};
          p1(i) = p1(i) - p0(i);
          return p1.l2();
        };
        auto ref_edge_length = get_length(ref_coords);
        auto spatial_edge_length = get_length(spatial_coords);
        auto change = spatial_edge_length / ref_edge_length;
        if (change >= 1. + edge_growth_thresh) {
          refine_edges.insert(e);
        }
      }
 
    if (refine_edges.size())
      refined = true;
 
    Range prev_level_ents;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        bit, BitRefLevel().set(), MBEDGE, prev_level_ents);
 
    Range prev_ref_ents;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        BitRefLevel().set(REFINED_EDGES_BIT), BitRefLevel().set(), MBEDGE,
        prev_ref_ents);
 
    refine_edges.merge(intersect(prev_ref_ents, prev_level_ents));
    CHKERR mField.getInterface<BitRefManager>()->setNthBitRefLevel(
        refine_edges, REFINED_EDGES_BIT, true);
 
    auto refine = mField.getInterface<MeshRefinement>();
    auto refined_bit = BitRefLevel().set(FIRST_REF_BIT + level);
    CHKERR refine->addVerticesInTheMiddleOfEdges(refine_edges, refined_bit);
    Range tris_to_refine;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
        bit, BitRefLevel().set(), SPACE_DIM, tris_to_refine);
    CHKERR refine->refineTris(tris_to_refine, refined_bit, QUIET, false);
 
    auto set_spatial_coords = [&]() {
      MoFEMFunctionBegin;
      Range new_vertices;
      CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
          refined_bit, refined_bit, MBVERTEX, new_vertices);
      auto th_parent_handle =
          mField.getInterface<BitRefManager>()->get_th_RefParentHandle();
      for (auto v : new_vertices) {
        EntityHandle parent;
        CHKERR mField.get_moab().tag_get_data(th_parent_handle, &v, 1, &parent);
        const EntityHandle *conn;
        int num_nodes;
        CHKERR mField.get_moab().get_connectivity(parent, conn, num_nodes,
                                                  true);
        std::array<double, 6> spatial_coords;
        CHKERR mField.get_moab().tag_get_data(th_spatial_coords, conn,
                                              num_nodes, spatial_coords.data());
        for (auto d = 0; d < 3; ++d) {
          spatial_coords[d] += spatial_coords[3 + d];
        }
        for (auto d = 0; d < 3; ++d) {
          spatial_coords[d] /= 2.0;
        }
        CHKERR mField.get_moab().tag_set_data(th_spatial_coords, &v, 1,
                                              spatial_coords.data());
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_spatial_coords();
 
    CHKERR mField.getInterface<BitRefManager>()->writeBitLevelByDim(
        BitRefLevel().set(FIRST_REF_BIT + level), BitRefLevel().set(), 2,
        ("A_refined_" + boost::lexical_cast<std::string>(level) + ".vtk")
            .c_str(),
        "VTK", "");
 
    MoFEMFunctionReturnHot(0);
  };
 
  auto reset_ref_bits = [](EntityHandle ent, BitRefLevel &bit) {
    bit.set(STORAGE_BIT, true);
    for (int l = 0; l < num_refinement_levels; ++l) {
      bit.set(FIRST_REF_BIT + l, false);
    }
  };
  CHKERR mField.getInterface<BitRefManager>()->lambdaBitRefLevel(
      reset_ref_bits);
 
  for (auto l = 0; l < num_refinement_levels; ++l) {
    CHKERR refined_mesh(l);
  };
 
  CHKERR mField.getInterface<BitRefManager>()->writeBitLevelByDim(
      BitRefLevel().set(FLIPPED_BIT + num_refinement_levels),
      BitRefLevel().set(), 2, "new_elements_after_edge_splits.vtk", "VTK", "");
 
  MoFEMFunctionReturn(0);
};
//! [Refine edges]
 
//! [Run problem]
MoFEMErrorCode Example::runProblem() {
  MoFEMFunctionBegin;
  CHKERR createCommonData();
  CHKERR setupProblem();
  if (!restart_flg)
    CHKERR initialConditions();
  CHKERR mechanicalBC(comp_mesh_bit, BitRefLevel().set());
  CHKERR thermalBC(comp_mesh_bit, BitRefLevel().set());
  CHKERR OPs();
  CHKERR tsSolve();
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Set up problem]
MoFEMErrorCode Example::setupProblem() {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  Range domain_ents;
  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
                                                     true);
  auto get_ents_by_dim = [&](const auto dim) {
    if (dim == SPACE_DIM) {
      return domain_ents;
    } else {
      Range ents;
      if (dim == 0)
        CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
      else
        CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
      return ents;
    }
  };
 
  auto get_base = [&]() {
    auto domain_ents = get_ents_by_dim(SPACE_DIM);
    if (domain_ents.empty())
      CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
    const auto type = type_from_handle(domain_ents[0]);
    switch (type) {
    case MBQUAD:
      return DEMKOWICZ_JACOBI_BASE;
    case MBHEX:
      return DEMKOWICZ_JACOBI_BASE;
    case MBTRI:
      return AINSWORTH_LEGENDRE_BASE;
    case MBTET:
      return AINSWORTH_LEGENDRE_BASE;
    default:
      CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
    }
    return NOBASE;
  };
 
  const auto base = get_base();
  MOFEM_LOG("PLASTICITY", Sev::inform)
      << "Base " << ApproximationBaseNames[base];
 
  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
  CHKERR simple->addDomainField("EP", L2, base, size_symm);
  CHKERR simple->addDomainField("TAU", L2, base, 1);
  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
 
  constexpr auto flux_space = (SPACE_DIM == 2) ? HCURL : HDIV;
  CHKERR simple->addDomainField("T", L2, base, 1);
  CHKERR simple->addDomainField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
  CHKERR simple->addBoundaryField("FLUX", flux_space, DEMKOWICZ_JACOBI_BASE, 1);
  CHKERR simple->addBoundaryField("TBC", L2, base, 1);
 
  // CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
 
  CHKERR simple->setFieldOrder("U", order);
  CHKERR simple->setFieldOrder("EP", ep_order);
  CHKERR simple->setFieldOrder("TAU", tau_order);
  CHKERR simple->setFieldOrder("FLUX", flux_order);
  CHKERR simple->setFieldOrder("T", T_order);
  CHKERR simple->setFieldOrder("TBC", T_order);
 
  // CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
 
#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; ///< blocs interation is collective, so that is set irrespective
                ///< if there are entities 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()));
  CHKERR simple->setFieldOrder("SIGMA", 0);
  CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
#endif
 
  CHKERR simple->setUp(is_distributed_mesh);
  // Cannot resetup after this block has been defined
  // as empty when remeshing
  if (qual_tol < std::numeric_limits<double>::epsilon() &&
      num_refinement_levels == 0)
    CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
 
  // auto project_ho_geometry = [&]() {
  //   Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
  //   return mField.loop_dofs("GEOMETRY", ent_method);
  // };
  // PetscBool project_geometry = PETSC_TRUE;
  // CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-project_geometry",
  //                            &project_geometry, PETSC_NULLPTR);
  // if (project_geometry) {
  //   CHKERR project_ho_geometry();
  // }
 
  MoFEMFunctionReturn(0);
}
//! [Set up problem]
 
//! [Create common data]
MoFEMErrorCode Example::createCommonData() {
  MoFEMFunctionBegin;
 
  auto get_command_line_parameters = [&]() {
    MoFEMFunctionBegin;
 
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-quasi_static",
                               &is_quasi_static, PETSC_NULLPTR);
    MOFEM_LOG("PLASTICITY", Sev::inform)
        << "Is quasi static: " << (is_quasi_static ? "true" : "false");
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-ts_dt", &init_dt,
                                 PETSC_NULLPTR);
    MOFEM_LOG("THERMOPLASTICITY", Sev::inform) << "Initial dt: " << init_dt;
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-ts_adapt_dt_min", &min_dt,
                                 PETSC_NULLPTR);
    MOFEM_LOG("THERMOPLASTICITY", Sev::inform) << "Minimum dt: " << min_dt;
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-min_quality", &qual_tol,
                                 PETSC_NULLPTR);
    MOFEM_LOG("REMESHING", Sev::inform)
        << "Minumum quality for edge flipping: " << qual_tol;
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "",
                                 "-quality_improvement_threshold", &qual_thresh,
                                 PETSC_NULLPTR);
    MOFEM_LOG("REMESHING", Sev::inform)
        << "Quality improvement threshold for edge flipping: " << qual_thresh;
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-edge_growth_threshold",
                                 &edge_growth_thresh, PETSC_NULLPTR);
    MOFEM_LOG("REMESHING", Sev::inform)
        << "Edge growth threshold for edge splitting: " << edge_growth_thresh;
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-num_refinement_levels",
                              &num_refinement_levels, PETSC_NULLPTR);
    MOFEM_LOG("REMESHING", Sev::inform)
        << "Number of refinement levels for edge splitting: "
        << num_refinement_levels;
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-scale", &scale,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus",
                                 &young_modulus, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio",
                                 &poisson_ratio, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-hardening", &H,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-hardening_viscous", &visH,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-yield_stress", &sigmaY,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn0", &cn0,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn1", &cn1,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-zeta", &zeta,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-Qinf", &Qinf,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-b_iso", &b_iso,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-C1_k", &C1_k,
                                 PETSC_NULLPTR);
    // CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-large_strains",
    //                            &is_large_strains, PETSC_NULLPTR);
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-set_timer", &set_timer,
                               PETSC_NULLPTR);
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-coeff_expansion",
                                 &default_coeff_expansion, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-ref_temp",
                                 &default_ref_temp, PETSC_NULLPTR);
    CHKERR PetscOptionsGetEList("", "-IC_type", ICTypes, 3,
                                (PetscInt *)&ic_type, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-init_temp", &init_temp,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-peak_temp", &peak_temp,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-width", &width,
                                 PETSC_NULLPTR);
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-capacity",
                                 &default_heat_capacity, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-conductivity",
                                 &default_heat_conductivity, PETSC_NULLPTR);
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-omega_0", &omega_0,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-omega_h", &omega_h,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-inelastic_heat_fraction",
                                 &inelastic_heat_fraction, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-temp_0", &temp_0,
                                 PETSC_NULLPTR);
 
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order,
                              PETSC_NULLPTR);
    PetscBool tau_order_is_set; ///< true if tau order is set
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-tau_order", &tau_order,
                              &tau_order_is_set);
    PetscBool ep_order_is_set; ///< true if tau order is set
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-ep_order", &ep_order,
                              &ep_order_is_set);
    PetscBool flux_order_is_set; ///< true if flux order is set
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-flux_order", &flux_order,
                              &flux_order_is_set);
    PetscBool T_order_is_set; ///< true if T order is set
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-T_order", &T_order,
                              &T_order_is_set);
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-geom_order", &geom_order,
                              PETSC_NULLPTR);
 
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
                              PETSC_NULLPTR);
 
    CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-restart", restart_file,
                                 255, &restart_flg);
    sscanf(restart_file, "restart_%d.dat", &restart_step);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-restart_time",
                                 &restart_time, PETSC_NULLPTR);
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-rho", &rho,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-alpha_damping",
                                 &alpha_damping, PETSC_NULLPTR);
 
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Young modulus " << young_modulus;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Poisson ratio " << poisson_ratio;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Yield stress " << sigmaY;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Hardening " << H;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Viscous hardening " << visH;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation yield stress " << Qinf;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation exponent " << b_iso;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Kinematic hardening " << C1_k;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "cn0 " << cn0;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "cn1 " << cn1;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "zeta " << zeta;
 
    MOFEM_LOG("THERMAL", Sev::inform)
        << "default_ref_temp " << default_ref_temp;
    MOFEM_LOG("THERMAL", Sev::inform) << "init_temp " << init_temp;
    MOFEM_LOG("THERMAL", Sev::inform) << "peak_temp " << peak_temp;
    MOFEM_LOG("THERMAL", Sev::inform) << "width " << width;
    MOFEM_LOG("THERMAL", Sev::inform)
        << "default_coeff_expansion " << default_coeff_expansion;
    MOFEM_LOG("THERMAL", Sev::inform)
        << "heat_conductivity " << default_heat_conductivity;
    MOFEM_LOG("THERMAL", Sev::inform)
        << "heat_capacity " << default_heat_capacity;
 
    MOFEM_LOG("THERMOPLASTICITY", Sev::inform)
        << "inelastic_heat_fraction " << inelastic_heat_fraction;
    MOFEM_LOG("THERMOPLASTICITY", Sev::inform) << "omega_0 " << omega_0;
    MOFEM_LOG("THERMOPLASTICITY", Sev::inform) << "omega_h " << omega_h;
 
    if (tau_order_is_set == PETSC_FALSE)
      tau_order = order - 2;
    if (ep_order_is_set == PETSC_FALSE)
      ep_order = order - 1;
 
    if (flux_order_is_set == PETSC_FALSE)
      flux_order = order;
    if (T_order_is_set == PETSC_FALSE)
      T_order = order - 1;
 
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Approximation order " << order;
    MOFEM_LOG("PLASTICITY", Sev::inform)
        << "Ep approximation order " << ep_order;
    MOFEM_LOG("PLASTICITY", Sev::inform)
        << "Tau approximation order " << tau_order;
    MOFEM_LOG("THERMAL", Sev::inform)
        << "FLUX approximation order " << flux_order;
    MOFEM_LOG("THERMAL", Sev::inform) << "T approximation order " << T_order;
    // MOFEM_LOG("PLASTICITY", Sev::inform)
    //     << "Geometry approximation order " << geom_order;
 
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Density " << rho;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "alpha_damping " << alpha_damping;
 
    PetscBool is_scale = PETSC_TRUE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_scale", &is_scale,
                               PETSC_NULLPTR);
    MOFEM_LOG("THERMOPLASTICITY", Sev::inform) << "Scaling used? " << is_scale;
 
    switch (is_scale) {
    case PETSC_TRUE:
      MOFEM_LOG("THERMOPLASTICITY", Sev::warning)
          << "Scaling does not yet work with radiation and convection BCs";
      // FIXME: Need to properly sort scaling as well as when using convection
      // and radiation BCs
      scale /= young_modulus;
      thermal_conductivity_scaling = [&](double thermal_cond, double heat_cap) {
        if (heat_cap != 0.0 && thermal_cond != 0.0)
          return 1. / (thermal_cond * heat_cap);
        else
          return 1. / thermal_cond;
      };
      heat_capacity_scaling = [&](double thermal_cond, double heat_cap) {
        if (heat_cap != 0.0)
          return 1. / (thermal_cond * heat_cap);
        else
          return 1.0;
      };
      inelastic_heat_fraction_scaling =
          [&](double young_modulus, double thermal_cond, double heat_cap) {
            if (heat_cap != 0.0)
              return young_modulus / (thermal_cond * heat_cap);
            else
              return young_modulus / thermal_cond;
          };
      break;
    default:
      thermal_conductivity_scaling = [&](double thermal_cond, double heat_cap) {
        return 1.0;
      };
      heat_capacity_scaling = [&](double thermal_cond, double heat_cap) {
        return 1.0;
      };
      inelastic_heat_fraction_scaling = [&](double young_modulus,
                                            double thermal_cond,
                                            double heat_cap) { return 1.0; };
      break;
    }
 
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Scale " << scale;
    MOFEM_LOG("THERMAL", Sev::inform)
        << "Thermal conductivity scale "
        << thermal_conductivity_scaling(default_heat_conductivity,
                                        default_heat_capacity);
    MOFEM_LOG("THERMAL", Sev::inform)
        << "Thermal capacity scale "
        << heat_capacity_scaling(default_heat_conductivity,
                                 default_heat_capacity);
    MOFEM_LOG("THERMOPLASTICITY", Sev::inform)
        << "Inelastic heat fraction scale "
        << inelastic_heat_fraction_scaling(
               young_modulus, default_heat_conductivity, default_heat_capacity);
 
#ifdef ADD_CONTACT
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn_contact",
                                 &ContactOps::cn_contact, PETSC_NULLPTR);
    MOFEM_LOG("CONTACT", Sev::inform)
        << "cn_contact " << ContactOps::cn_contact;
#endif // ADD_CONTACT
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR get_command_line_parameters();
 
#ifdef ADD_CONTACT
  #ifdef ENABLE_PYTHON_BINDING
  sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
  CHKERR sdfPythonPtr->sdfInit("sdf.py");
  ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
  #endif
#endif // ADD_CONTACT
 
  MoFEMFunctionReturn(0);
}
//! [Create common data]
 
//! [Initial conditions]
MoFEMErrorCode Example::initialConditions() {
  MoFEMFunctionBegin;
 
  auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order; };
 
  // #ifdef PYTHON_INIT_SURFACE
  //   auto get_py_surface_init = []() {
  //     auto py_surf_init = boost::make_shared<SurfacePython>();
  //     CHKERR py_surf_init->surfaceInit("surface.py");
  //     surfacePythonWeakPtr = py_surf_init;
  //     return py_surf_init;
  //   };
  //   auto py_surf_init = get_py_surface_init();
  // #endif
 
  auto simple = mField.getInterface<Simple>();
 
  //   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);
  //   CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
  //                                            "REMOVE_ALL", "U", 0, 3);
 
  // #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");
 
  using UDO = ForcesAndSourcesCore::UserDataOperator;
 
  auto T_ptr = boost::make_shared<VectorDouble>();
 
  auto post_proc = [&](auto dm) {
    MoFEMFunctionBegin;
 
    auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("T", T_ptr));
    // post_proc_fe->getOpPtrVector().push_back(
    //     new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
    if (atom_test == 8) {
 
      auto TAU_ptr = boost::make_shared<VectorDouble>();
      auto EP_ptr = boost::make_shared<MatrixDouble>();
 
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateScalarFieldValues("TAU", TAU_ptr));
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>("EP", EP_ptr));
 
      post_proc_fe->getOpPtrVector().push_back(
 
          new OpPPMap(post_proc_fe->getPostProcMesh(),
                      post_proc_fe->getMapGaussPts(),
 
                      {{"T", T_ptr}, {"TAU", TAU_ptr}},
 
                      {},
 
                      {},
 
                      {{"EP", EP_ptr}}
 
                      )
 
      );
    } else {
      post_proc_fe->getOpPtrVector().push_back(
 
          new OpPPMap(post_proc_fe->getPostProcMesh(),
                      post_proc_fe->getMapGaussPts(),
 
                      {{"T", T_ptr}},
 
                      {},
 
                      {},
 
                      {}
 
                      )
 
      );
    }
 
    CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
    CHKERR post_proc_fe->writeFile("out_init.h5m");
 
    MoFEMFunctionReturn(0);
  };
 
  auto solve_init = [&]() {
    MoFEMFunctionBegin;
 
    auto set_domain_rhs = [&](auto &pip) {
      MoFEMFunctionBegin;
 
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV});
 
      pip.push_back(new OpCalculateScalarFieldValues("T", T_ptr));
      pip.push_back(new OpRhsSetInitT<AT, IT>(
          "T", nullptr, T_ptr, nullptr, boost::make_shared<double>(init_temp),
          boost::make_shared<double>(peak_temp)));
 
      if (atom_test == 8) {
        auto TAU_ptr = boost::make_shared<VectorDouble>();
        auto EP_ptr = boost::make_shared<MatrixDouble>();
 
        pip.push_back(new OpCalculateScalarFieldValues("TAU", TAU_ptr));
        auto min_tau = boost::make_shared<double>(1.0);
        auto max_tau = boost::make_shared<double>(2.0);
        pip.push_back(new OpRhsSetInitT<AT, IT>("TAU", nullptr, TAU_ptr,
                                                nullptr, min_tau, max_tau));
 
        pip.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
            "EP", EP_ptr));
        auto min_EP = boost::make_shared<double>(0.0);
        auto max_EP = boost::make_shared<double>(0.01);
        pip.push_back(new EdgeFlippingOps::OpRhsSetInitEP<AT, IT, SPACE_DIM>(
            "EP", nullptr, EP_ptr, nullptr, min_EP, max_EP));
      }
 
      // using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
      //     AT>::template LinearForm<IT>;
      // using OpInternalForceCauchy =
      //     typename B::template OpGradTimesSymTensor<1, SPACE_DIM, SPACE_DIM>;
      // using OpInternalForcePiola =
      //     typename B::template OpGradTimesTensor<1, SPACE_DIM, SPACE_DIM>;
 
      // using P = PlasticityIntegrators<DomainEleOp>;
 
      // auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
      // common_thermoelastic_ptr, common_thermoplastic_ptr] =
      //     createCommonThermoPlasticOps<SPACE_DIM, IT, DomainEleOp>(
      //         mField, "MAT_ELASTIC", "MAT_THERMAL", "MAT_THERMOPLASTIC", pip,
      //         "U", "EP", "TAU", "T", scale, thermal_conductivity_scale,
      //         heat_capacity_scale, inelastic_heat_fraction_scale,
      //         Sev::inform);
      // auto m_D_ptr = common_hencky_ptr->matDPtr;
 
      // using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
      //     AT>::template LinearForm<IT>;
      // using H = HenckyOps::HenckyIntegrators<DomainEleOp>;
      // auto coeff_expansion_ptr = common_thermal_ptr->getCoeffExpansionPtr();
      // auto ref_temp_ptr = common_thermal_ptr->getRefTempPtr();
      // pip.push_back(new
      //               typename H::template
      //               OpCalculateHenckyThermalStress<SPACE_DIM, IT>(
      //                   "U", T_ptr, common_hencky_ptr,
      //                   coeff_expansion_ptr, ref_temp_ptr));
      // pip.push_back(new typename H::template
      // OpCalculatePiolaStress<SPACE_DIM, IT>(
      //     "U", common_hencky_ptr));
      // using OpInternalForcePiola =
      //     typename B::template OpGradTimesTensor<1, SPACE_DIM, SPACE_DIM>;
      // pip.push_back(new OpInternalForcePiola(
      //     "U", common_hencky_ptr->getMatFirstPiolaStress()));
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_domain_lhs = [&](auto &pip) {
      MoFEMFunctionBegin;
 
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV});
 
      using OpLhsScalarLeastSquaresProj = FormsIntegrators<
          DomainEleOp>::Assembly<AT>::BiLinearForm<IT>::OpMass<1, 1>;
      pip.push_back(new OpLhsScalarLeastSquaresProj("T", "T"));
      if (atom_test == 8) {
        pip.push_back(new OpLhsScalarLeastSquaresProj("TAU", "TAU"));
        pip.push_back(
            new EdgeFlippingOps::OpLhsSymTensorLeastSquaresProj<SPACE_DIM, IT>(
                "EP", "EP"));
      }
 
      // auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
      // common_thermoelastic_ptr, common_thermoplastic_ptr] =
      //     createCommonThermoPlasticOps<SPACE_DIM, IT, DomainEleOp>(
      //         mField, "MAT_ELASTIC", "MAT_THERMAL", "MAT_THERMOPLASTIC", pip,
      //         "U", "EP", "TAU", "T", scale, thermal_conductivity_scale,
      //         heat_capacity_scale, inelastic_heat_fraction_scale,
      //         Sev::inform);
 
      // auto m_D_ptr = common_hencky_ptr->matDPtr;
 
      // using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
      //     AT>::template BiLinearForm<IT>;
      // using OpKPiola = typename B::template OpGradTensorGrad<1, SPACE_DIM,
      // SPACE_DIM, 1>;
 
      // using H = HenckyIntegrators<DomainEleOp>;
      // pip.push_back(new OpCalculateScalarFieldValues("T", T_ptr));
      // auto coeff_expansion_ptr = common_thermal_ptr->getCoeffExpansionPtr();
      // auto ref_temp_ptr = common_thermal_ptr->getRefTempPtr();
      // pip.push_back(new
      //               typename H::template
      //               OpCalculateHenckyThermalStress<SPACE_DIM, IT>(
      //                   "U", T_ptr, common_hencky_ptr,
      //                   coeff_expansion_ptr, ref_temp_ptr));
      // pip.push_back(new typename H::template
      // OpCalculatePiolaStress<SPACE_DIM, IT>(
      //     "U", common_hencky_ptr));
      // pip.push_back(new typename H::template OpHenckyTangent<SPACE_DIM, IT>(
      //     "U", common_hencky_ptr));
      // pip.push_back(new OpKPiola("U", "U",
      //                            common_hencky_ptr->getMatTangent()));
      // pip.push_back(new typename HenckyOps::OpCalculateHenckyThermalStressdT<
      //               SPACE_DIM, IT, AssemblyDomainEleOp>("U", "T",
      //                                            common_hencky_ptr,
      //                                            coeff_expansion_ptr));
 
      MoFEMFunctionReturn(0);
    };
 
    auto create_sub_dm = [&](SmartPetscObj<DM> base_dm) {
      auto dm_sub = createDM(mField.get_comm(), "DMMOFEM");
      CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "INIT_DM");
      CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
      CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
      for (auto f : {"T"}) {
        CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
        CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
      }
      if (atom_test == 8) {
        for (auto f : {"TAU", "EP"}) {
          CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
          CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
        }
      }
      CHKERR DMSetUp(dm_sub);
      return dm_sub;
    };
 
    auto fe_rhs = boost::make_shared<DomainEle>(mField);
    auto fe_lhs = boost::make_shared<DomainEle>(mField);
    fe_rhs->getRuleHook = vol_rule;
    fe_lhs->getRuleHook = vol_rule;
    CHKERR set_domain_rhs(fe_rhs->getOpPtrVector());
    CHKERR set_domain_lhs(fe_lhs->getOpPtrVector());
 
    auto sub_dm = create_sub_dm(simple->getDM());
 
    auto null_fe = boost::shared_ptr<FEMethod>();
    CHKERR DMMoFEMKSPSetComputeOperators(sub_dm, simple->getDomainFEName(),
                                         fe_lhs, null_fe, null_fe);
    CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(), fe_rhs,
                                   null_fe, null_fe);
 
    auto ksp = MoFEM::createKSP(mField.get_comm());
    CHKERR KSPSetDM(ksp, sub_dm);
    CHKERR KSPSetFromOptions(ksp);
 
    auto D = createDMVector(sub_dm);
 
    CHKERR KSPSolve(ksp, PETSC_NULLPTR, D);
 
    CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToGlobalVector(sub_dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR solve_init();
  CHKERR post_proc(simple->getDM());
 
  MOFEM_LOG("THERMAL", Sev::inform) << "Set thermoelastic initial conditions";
 
  MoFEMFunctionReturn(0);
}
//! [Initial conditions]
 
//! [Mechanical boundary conditions]
MoFEMErrorCode Example::mechanicalBC(BitRefLevel bit, BitRefLevel mask) {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
 
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
                                           "U", 0, 0, true,
                                           is_distributed_mesh);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
                                           "U", 1, 1, true,
                                           is_distributed_mesh);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
                                           "U", 2, 2, true,
                                           is_distributed_mesh);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                           "REMOVE_ALL", "U", 0, 3, true,
                                           is_distributed_mesh);
 
#ifdef ADD_CONTACT
  for (auto b : {"FIX_X", "REMOVE_X"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(
        simple->getProblemName(), b, "SIGMA", 0, 0, false, is_distributed_mesh);
  for (auto b : {"FIX_Y", "REMOVE_Y"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(
        simple->getProblemName(), b, "SIGMA", 1, 1, false, is_distributed_mesh);
  for (auto b : {"FIX_Z", "REMOVE_Z"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(
        simple->getProblemName(), b, "SIGMA", 2, 2, false, is_distributed_mesh);
  for (auto b : {"FIX_ALL", "REMOVE_ALL"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(
        simple->getProblemName(), b, "SIGMA", 0, 3, false, is_distributed_mesh);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                           "NO_CONTACT", "SIGMA", 0, 3, false,
                                           is_distributed_mesh);
#endif
 
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "U");
 
  auto &bc_map = bc_mng->getBcMapByBlockName();
  for (auto bc : bc_map)
    MOFEM_LOG("PLASTICITY", Sev::verbose) << "Marker " << bc.first;
 
  MoFEMFunctionReturn(0);
}
//! [Mechanical boundary conditions]
 
//! [Thermal boundary conditions]
MoFEMErrorCode Example::thermalBC(BitRefLevel bit, BitRefLevel mask) {
  MoFEMFunctionBegin;
 
  MOFEM_LOG("SYNC", Sev::noisy) << "bC";
  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
 
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
 
  auto get_skin = [&]() {
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
    CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
        bit, BitRefLevel().set(), 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_flux_blocks = [&](auto skin, bool temp_bc = false) {
    auto remove_cubit_blocks = [&](auto c) {
      MoFEMFunctionBegin;
      for (auto m :
 
           mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(c)
 
      ) {
        Range ents;
        CHKERR mField.get_moab().get_entities_by_dimension(
            m->getMeshset(), SPACE_DIM - 1, ents, true);
        skin = subtract(skin, ents);
      }
      MoFEMFunctionReturn(0);
    };
 
    auto remove_named_blocks = [&](auto n) {
      MoFEMFunctionBegin;
      for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
               std::regex(
 
                   (boost::format("%s(.*)") % n).str()
 
                       ))
 
      ) {
        Range ents;
        CHKERR mField.get_moab().get_entities_by_dimension(
            m->getMeshset(), SPACE_DIM - 1, ents, true);
        skin = subtract(skin, ents);
      }
      MoFEMFunctionReturn(0);
    };
    if (!temp_bc) {
      CHK_THROW_MESSAGE(remove_cubit_blocks(NODESET | TEMPERATURESET),
                        "remove_cubit_blocks");
      CHK_THROW_MESSAGE(remove_named_blocks("TEMPERATURE"),
                        "remove_named_blocks");
    }
    CHK_THROW_MESSAGE(remove_cubit_blocks(SIDESET | HEATFLUXSET),
                      "remove_cubit_blocks");
    CHK_THROW_MESSAGE(remove_named_blocks("HEATFLUX"), "remove_named_blocks");
    CHK_THROW_MESSAGE(remove_named_blocks("CONVECTION"), "remove_named_blocks");
    CHK_THROW_MESSAGE(remove_named_blocks("RADIATION"), "remove_named_blocks");
    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 remove_flux_ents = filter_true_skin(filter_flux_blocks(get_skin()));
  auto remove_temp_bc_ents =
      filter_true_skin(filter_flux_blocks(get_skin(), true));
 
  // CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
  //     remove_flux_ents);
  // CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
  //     remove_temp_bc_ents);
 
  MOFEM_LOG("SYNC", Sev::noisy) << remove_flux_ents << endl;
  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
 
  MOFEM_LOG("SYNC", Sev::noisy) << remove_temp_bc_ents << endl;
  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
 
#ifndef NDEBUG
 
  CHKERR save_range(
      mField.get_moab(),
      (boost::format("flux_remove_%d.vtk") % mField.get_comm_rank()).str(),
      remove_flux_ents);
 
#endif
 
  if (is_distributed_mesh == PETSC_TRUE) {
    CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
        simple->getProblemName(), "FLUX", remove_flux_ents);
    CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
        simple->getProblemName(), "TBC", remove_temp_bc_ents);
  } else {
    CHKERR mField.getInterface<ProblemsManager>()
        ->removeDofsOnEntitiesNotDistributed(simple->getProblemName(), "FLUX",
                                             remove_flux_ents);
    CHKERR mField.getInterface<ProblemsManager>()
        ->removeDofsOnEntitiesNotDistributed(simple->getProblemName(), "TBC",
                                             remove_temp_bc_ents);
  }
 
  // auto set_init_temp = [](boost::shared_ptr<FieldEntity> field_entity_ptr) {
  //   field_entity_ptr->getEntFieldData()[0] = init_temp;
  //   return 0;
  // };
  // CHKERR
  // mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(set_init_temp,
  //                                                                "T");
 
  CHKERR bc_mng->pushMarkDOFsOnEntities<HeatFluxCubitBcData>(
      simple->getProblemName(), "FLUX", false);
 
  MoFEMFunctionReturn(0);
}
//! [Thermal Boundary conditions]
 
//! [Push operators to pipeline]
MoFEMErrorCode Example::OPs() {
  MoFEMFunctionBegin;
  auto pip_mng = mField.getInterface<PipelineManager>();
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
 
  auto boundary_marker =
      bc_mng->getMergedBlocksMarker(vector<string>{"HEATFLUX"});
 
  auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
 
  auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order; };
 
  CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
  CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
 
  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
 
  auto thermal_block_params =
      boost::make_shared<ThermoElasticOps::BlockedThermalParameters>();
  auto heat_conductivity_ptr = thermal_block_params->getHeatConductivityPtr();
  auto heat_capacity_ptr = thermal_block_params->getHeatCapacityPtr();
 
  auto thermoelastic_block_params =
      boost::make_shared<ThermoElasticOps::BlockedThermoElasticParameters>();
  auto coeff_expansion_ptr = thermoelastic_block_params->getCoeffExpansionPtr();
  auto ref_temp_ptr = thermoelastic_block_params->getRefTempPtr();
 
  // Default time scaling of BCs and sources, from command line arguments
  auto time_scale =
      boost::make_shared<TimeScale>("", false, [](const double) { return 1; });
  auto def_time_scale = [time_scale](const double time) {
    return (!time_scale->argFileScale) ? time : 1;
  };
  auto def_file_name = [time_scale](const std::string &&name) {
    return (!time_scale->argFileScale) ? name : "";
  };
 
  // Files which define scaling for separate variables, if provided
  auto time_bodyforce_scaling = boost::make_shared<TimeScale>(
      def_file_name("bodyforce_scale.txt"), false, def_time_scale);
  auto time_heatsource_scaling = boost::make_shared<TimeScale>(
      def_file_name("heatsource_scale.txt"), false, def_time_scale);
  auto time_temperature_scaling = boost::make_shared<TimeScale>(
      def_file_name("temperature_bc_scale.txt"), false, def_time_scale);
  auto time_displacement_scaling = boost::make_shared<TimeScale>(
      def_file_name("displacement_bc_scale.txt"), false, def_time_scale);
  auto time_flux_scaling = boost::make_shared<TimeScale>(
      def_file_name("flux_bc_scale.txt"), false, def_time_scale);
  auto time_force_scaling = boost::make_shared<TimeScale>(
      def_file_name("force_bc_scale.txt"), false, def_time_scale);
 
  auto add_boundary_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    pip.push_back(new OpSetBc("FLUX", true, boundary_marker));
 
    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV});
    pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
 
    // Add Natural BCs to LHS
    CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
        pip, mField, "U", Sev::inform);
 
#ifdef ADD_CONTACT
    CHKERR
    ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
        pip, "SIGMA", "U");
    CHKERR
    ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
        mField, pip, simple->getDomainFEName(), "SIGMA", "U", "", vol_rule);
#endif // ADD_CONTACT
 
    using T =
        NaturalBC<BoundaryEleOp>::template Assembly<PETSC>::BiLinearForm<GAUSS>;
    using OpConvectionLhsBC =
        T::OpFlux<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1>;
    using OpRadiationLhsBC =
        T::OpFlux<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1>;
    auto temp_bc_ptr = boost::make_shared<VectorDouble>();
    pip.push_back(new OpCalculateScalarFieldValues("TBC", temp_bc_ptr));
    T::AddFluxToPipeline<OpConvectionLhsBC>::add(pip, mField, "TBC", "TBC",
                                                 "CONVECTION", Sev::verbose);
    T::AddFluxToPipeline<OpRadiationLhsBC>::add(
        pip, mField, "TBC", "TBC", temp_bc_ptr, "RADIATION", Sev::verbose);
 
    // This is temporary implementation. It be moved to LinearFormsIntegrators,
    // however for hybridised case, what is goal of this changes, such function
    // is implemented for fluxes in broken space. Thus ultimately this operator
    // would be not needed.
 
    struct OpTBCTimesNormalFlux
        : public FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase {
 
      using OP = FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase;
 
      OpTBCTimesNormalFlux(const std::string row_field_name,
                           const std::string col_field_name,
                           boost::shared_ptr<Range> ents_ptr = nullptr)
          : OP(row_field_name, col_field_name, OP::OPROWCOL, ents_ptr) {
        this->sYmm = false;
        this->assembleTranspose = true;
        this->onlyTranspose = false;
      }
 
      MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                               EntitiesFieldData::EntData &col_data) {
        MoFEMFunctionBegin;
 
        FTENSOR_INDEX(SPACE_DIM, i);
 
        // get integration weights
        auto t_w = OP::getFTensor0IntegrationWeight();
        // get base function values on rows
        auto t_row_base = row_data.getFTensor0N();
        // get normal vector
        auto t_normal = OP::getFTensor1NormalsAtGaussPts();
        // loop over integration points
        for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
          // loop over rows base functions
          auto a_mat_ptr = &*OP::locMat.data().begin();
          int rr = 0;
          for (; rr != OP::nbRows; rr++) {
            // take into account Jacobian
            const double alpha = t_w * t_row_base;
            // get column base functions values at gauss point gg
            auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
            // loop over columns
            for (int cc = 0; cc != OP::nbCols; cc++) {
              // calculate element of local matrix
              // using L2 norm of normal vector for convenient area calculation
              // for quads, tris etc.
              *a_mat_ptr += alpha * (t_col_base(i) * t_normal(i));
              ++t_col_base;
              ++a_mat_ptr;
            }
            ++t_row_base;
          }
          for (; rr < OP::nbRowBaseFunctions; ++rr)
            ++t_row_base;
          ++t_normal;
          ++t_w; // move to another integration weight
        }
        EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
        if (fe_type == MBTRI) {
          OP::locMat /= 2;
        }
        MoFEMFunctionReturn(0);
      }
    };
 
    pip.push_back(new OpTBCTimesNormalFlux("TBC", "FLUX"));
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_boundary_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV});
    pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
 
    CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
        pip, mField, "U", {time_scale, time_force_scaling}, "FORCE",
        "PRESSURE", Sev::inform);
 
    // Temperature BC
 
    using OpTemperatureBC =
        BoundaryNaturalBC::OpFlux<NaturalTemperatureMeshsets, 3, SPACE_DIM>;
    CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpTemperatureBC>::add(
        pip, mField, "FLUX", {time_scale, time_temperature_scaling},
        "TEMPERATURE", Sev::inform);
 
    // Note: fluxes for temperature should be aggregated. Here separate is
    // NaturalMeshsetType<HEATFLUXSET>, we should add version with BLOCKSET,
    // convection, see example tutorials/vec-0/src/NaturalBoundaryBC.hpp.
    // and vec-0/elastic.cpp
 
    using OpFluxBC =
        BoundaryNaturalBC::OpFlux<NaturalMeshsetType<HEATFLUXSET>, 1, 1>;
    CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpFluxBC>::add(
        pip, mField, "TBC", {time_scale, time_flux_scaling}, "FLUX",
        Sev::inform);
 
    using T = NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
    using OpConvectionRhsBC =
        T::OpFlux<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1>;
    using OpRadiationRhsBC =
        T::OpFlux<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1>;
    auto temp_bc_ptr = boost::make_shared<VectorDouble>();
    pip.push_back(new OpCalculateScalarFieldValues("TBC", temp_bc_ptr));
    T::AddFluxToPipeline<OpConvectionRhsBC>::add(
        pip, mField, "TBC", temp_bc_ptr, "CONVECTION", Sev::inform);
    T::AddFluxToPipeline<OpRadiationRhsBC>::add(pip, mField, "TBC", temp_bc_ptr,
                                                "RADIATION", Sev::inform);
 
    auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(
        new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
 
    // This is temporary implementation. It be moved to LinearFormsIntegrators,
    // however for hybridised case, what is goal of this changes, such function
    // is implemented for fluxes in broken space. Thus ultimately this operator
    // would be not needed.
 
    struct OpTBCTimesNormalFlux
        : public FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase {
 
      using OP = FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase;
 
      OpTBCTimesNormalFlux(const std::string field_name,
                           boost::shared_ptr<MatrixDouble> vec,
                           boost::shared_ptr<Range> ents_ptr = nullptr)
          : OP(field_name, field_name, OP::OPROW, ents_ptr), sourceVec(vec) {}
 
      MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
        MoFEMFunctionBegin;
        FTENSOR_INDEX(SPACE_DIM, i);
        // get integration weights
        auto t_w = OP::getFTensor0IntegrationWeight();
        // get base function values on rows
        auto t_row_base = row_data.getFTensor0N();
        // get normal vector
        auto t_normal = OP::getFTensor1NormalsAtGaussPts();
        // get vector values
        auto t_vec = getFTensor1FromMat<
            SPACE_DIM, 1, DataLayoutTraits<DataLayout::CoeffsByGauss>>(
            *sourceVec);
        // loop over integration points
        for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
          // take into account Jacobian
          const double alpha = t_w * (t_vec(i) * t_normal(i));
          // loop over rows base functions
          int rr = 0;
          for (; rr != OP::nbRows; ++rr) {
            OP::locF[rr] += alpha * t_row_base;
            ++t_row_base;
          }
          for (; rr < OP::nbRowBaseFunctions; ++rr)
            ++t_row_base;
          ++t_w; // move to another integration weight
          ++t_vec;
          ++t_normal;
        }
        EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
        if (fe_type == MBTRI) {
          OP::locF /= 2;
        }
        MoFEMFunctionReturn(0);
      }
 
    protected:
      boost::shared_ptr<MatrixDouble> sourceVec;
    };
    pip.push_back(new OpTBCTimesNormalFlux("TBC", mat_flux_ptr));
 
    struct OpBaseNormalTimesTBC
        : public FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase {
 
      using OP = FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase;
 
      OpBaseNormalTimesTBC(const std::string field_name,
                           boost::shared_ptr<VectorDouble> vec,
                           boost::shared_ptr<Range> ents_ptr = nullptr)
          : OP(field_name, field_name, OP::OPROW, ents_ptr), sourceVec(vec) {}
 
      MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
        MoFEMFunctionBegin;
        FTENSOR_INDEX(SPACE_DIM, i);
        // get integration weights
        auto t_w = OP::getFTensor0IntegrationWeight();
        // get base function values on rows
        auto t_row_base = row_data.getFTensor1N<3>();
        // get normal vector
        auto t_normal = OP::getFTensor1NormalsAtGaussPts();
        // get vector values
        auto t_vec = getFTensor0FromVec(*sourceVec);
        // loop over integration points
        for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
          // take into account Jacobian
          const double alpha = t_w * t_vec;
          // loop over rows base functions
          int rr = 0;
          for (; rr != OP::nbRows; ++rr) {
            OP::locF[rr] += alpha * (t_row_base(i) * t_normal(i));
            ++t_row_base;
          }
          for (; rr < OP::nbRowBaseFunctions; ++rr)
            ++t_row_base;
          ++t_w; // move to another integration weight
          ++t_vec;
          ++t_normal;
        }
        EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
        if (fe_type == MBTRI) {
          OP::locF /= 2;
        }
        MoFEMFunctionReturn(0);
      }
 
    protected:
      boost::shared_ptr<VectorDouble> sourceVec;
    };
 
    pip.push_back(new OpBaseNormalTimesTBC("FLUX", temp_bc_ptr));
 
#ifdef ADD_CONTACT
    CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
        pip, "SIGMA", "U");
#endif // ADD_CONTACT
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    pip.push_back(new OpSetBc("FLUX", true, boundary_marker));
 
    CHKERR ThermoElasticOps::addMatThermalBlockOps(
        mField, pip, "MAT_THERMAL", thermal_block_params,
        default_heat_conductivity, default_heat_capacity,
        default_thermal_conductivity_scale, default_heat_capacity_scale,
        Sev::inform, thermal_conductivity_scaling, heat_capacity_scaling);
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV});
 
    if (is_quasi_static == PETSC_FALSE) {
 
      //! [Only used for dynamics]
      using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
          IT>::OpMass<1, SPACE_DIM>;
      //! [Only used for dynamics]
 
      auto get_inertia_and_mass_damping = [this](const double, const double,
                                                 const double) {
        auto *pip = mField.getInterface<PipelineManager>();
        auto &fe_domain_lhs = pip->getDomainLhsFE();
        return (rho / scale) * fe_domain_lhs->ts_aa +
               (alpha_damping / scale) * fe_domain_lhs->ts_a;
      };
      pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
    }
 
    CHKERR opThermoPlasticFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
        mField, "MAT_PLASTIC", "MAT_THERMAL", "MAT_THERMOELASTIC",
        "MAT_THERMOPLASTIC", pip, "U", "EP", "TAU", "T");
 
    auto hencky_common_data_ptr =
        HenckyOps::commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
            mField, pip, "U", "MAT_PLASTIC", Sev::inform, scale);
    auto mat_D_ptr = hencky_common_data_ptr->matDPtr;
    auto mat_grad_ptr = hencky_common_data_ptr->matGradPtr;
 
    pip.push_back(new OpSetBc("FLUX", true, boundary_marker));
 
    auto resistance = [heat_conductivity_ptr](const double, const double,
                                              const double) {
      return (1. / (*heat_conductivity_ptr));
    };
    auto capacity = [heat_capacity_ptr](const double, const double,
                                        const double) {
      return -(*heat_capacity_ptr);
    };
 
    pip.push_back(new OpHdivHdiv("FLUX", "FLUX", resistance, mat_grad_ptr));
    pip.push_back(
        new OpHdivT("FLUX", "T", []() constexpr { return -1; }, true));
 
    auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(
        new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
 
    pip.push_back(
        new OpHdivU("FLUX", "U", mat_flux_ptr, resistance, mat_grad_ptr));
 
    auto op_capacity = new OpCapacity("T", "T", capacity);
    op_capacity->feScalingFun = [](const FEMethod *fe_ptr) {
      return fe_ptr->ts_a;
    };
    pip.push_back(op_capacity);
 
    pip.push_back(new OpUnSetBc("FLUX"));
 
    auto vec_temp_ptr = boost::make_shared<VectorDouble>();
    pip.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
    CHKERR DomainNaturalBCLhs::AddFluxToPipeline<OpSetTemperatureLhs>::add(
        pip, mField, "T", vec_temp_ptr, "SETTEMP", Sev::verbose);
 
    pip.push_back(new OpUnSetBc("FLUX"));
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    pip.push_back(new OpSetBc("FLUX", true, boundary_marker));
 
    CHKERR ThermoElasticOps::addMatThermalBlockOps(
        mField, pip, "MAT_THERMAL", thermal_block_params,
        default_heat_conductivity, default_heat_capacity,
        default_thermal_conductivity_scale, default_heat_capacity_scale,
        Sev::inform, thermal_conductivity_scaling, heat_capacity_scaling);
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV});
 
    auto hencky_common_data_ptr =
        HenckyOps::commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
            mField, pip, "U", "MAT_PLASTIC", Sev::inform, scale);
    auto mat_D_ptr = hencky_common_data_ptr->matDPtr;
    auto mat_grad_ptr = hencky_common_data_ptr->matGradPtr;
    auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
    auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
 
    auto vec_temp_ptr = boost::make_shared<VectorDouble>();
    auto vec_temp_dot_ptr = boost::make_shared<VectorDouble>();
    auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
    auto vec_temp_div_ptr = boost::make_shared<VectorDouble>();
 
    pip.push_back(new OpCalculateScalarFieldValues("T", vec_temp_ptr));
    pip.push_back(new OpCalculateScalarFieldValuesDot("T", vec_temp_dot_ptr));
    pip.push_back(new OpCalculateHdivVectorDivergence<3, SPACE_DIM>(
        "FLUX", vec_temp_div_ptr));
    pip.push_back(
        new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", mat_flux_ptr));
 
    auto resistance = [heat_conductivity_ptr](const double, const double,
                                              const double) {
      return (1. / (*heat_conductivity_ptr));
    };
    // negative value is consequence of symmetric system
    auto capacity = [heat_capacity_ptr](const double, const double,
                                        const double) {
      return -(*heat_capacity_ptr);
    };
    auto unity = [](const double, const double, const double) constexpr {
      return -1.;
    };
    pip.push_back(
        new OpHdivFlux("FLUX", mat_flux_ptr, resistance, mat_grad_ptr));
    pip.push_back(new OpHDivTemp("FLUX", vec_temp_ptr, unity));
    pip.push_back(new OpBaseDivFlux("T", vec_temp_div_ptr, unity));
    pip.push_back(new OpBaseDotT("T", vec_temp_dot_ptr, capacity));
 
    CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
        pip, mField, "U",
        {boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
        Sev::inform);
 
    // only in case of dynamics
    if (is_quasi_static == PETSC_FALSE) {
 
      //! [Only used for dynamics]
      using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
          AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
      //! [Only used for dynamics]
 
      auto mat_acceleration = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
          "U", mat_acceleration));
      pip.push_back(
          new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
            return rho / scale;
          }));
      if (alpha_damping > 0) {
        auto mat_velocity = boost::make_shared<MatrixDouble>();
        pip.push_back(
            new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
        pip.push_back(
            new OpInertiaForce("U", mat_velocity, [](double, double, double) {
              return alpha_damping / scale;
            }));
      }
    }
 
    CHKERR opThermoPlasticFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
        mField, "MAT_PLASTIC", "MAT_THERMAL", "MAT_THERMOELASTIC",
        "MAT_THERMOPLASTIC", pip, "U", "EP", "TAU", "T");
 
#ifdef ADD_CONTACT
    CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
        pip, "SIGMA", "U");
#endif // ADD_CONTACT
 
    pip.push_back(new OpUnSetBc("FLUX"));
 
    MoFEMFunctionReturn(0);
  };
 
  auto create_reaction_pipeline = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
    CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
        mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
    MoFEMFunctionReturn(0);
  };
 
  reactionFe = boost::make_shared<DomainEle>(mField);
  reactionFe->getRuleHook = vol_rule;
  CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
  reactionFe->postProcessHook =
      EssentialPreProcReaction<DisplacementCubitBcData>(mField, reactionFe);
 
  // Set BCs by eliminating degrees of freedom
  auto get_bc_hook_rhs = [&]() {
    EssentialPreProc<DisplacementCubitBcData> hook(
        mField, pip_mng->getDomainRhsFE(),
        {time_scale, time_displacement_scaling});
    return hook;
  };
  auto get_bc_hook_lhs = [&]() {
    EssentialPreProc<DisplacementCubitBcData> hook(
        mField, pip_mng->getDomainLhsFE(),
        {time_scale, time_displacement_scaling});
    return hook;
  };
 
  pip_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
  pip_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
 
  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
 
  // Boundary
  CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
  CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
 
  MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
 
//! [Solve]
struct SetUpSchur {
 
  /**
   * @brief Create data structure for handling Schur complement
   *
   * @param m_field
   * @param sub_dm  Schur complement sub dm
   * @param field_split_it IS of Schur block
   * @param ao_map AO map from sub dm to main problem
   * @return boost::shared_ptr<SetUpSchur>
   */
  static boost::shared_ptr<SetUpSchur> createSetUpSchur(
 
      MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
      SmartPetscObj<IS> field_split_it, SmartPetscObj<AO> ao_map
 
  );
  virtual MoFEMErrorCode setUp(TS solver) = 0;
 
protected:
  SetUpSchur() = default;
};
 
struct MyTsCtx {
  boost::shared_ptr<MoFEM::TsCtx> dmts_ctx;
  PetscInt snes_iter_counter = 0;
};
 
static std::unordered_map<TS, MyTsCtx *> ts_to_ctx;
 
MoFEMErrorCode TSIterationPreStage(TS ts, PetscReal stagetime) {
  MoFEMFunctionBegin;
  MyTsCtx *my_ctx = ts_to_ctx[ts];
  my_ctx->snes_iter_counter = 0;
  MoFEMFunctionReturn(0);
}
 
// Monitor function
MoFEMErrorCode SNESIterationMonitor(SNES snes, PetscInt its, PetscReal norm,
                                    void *ctx) {
  MoFEMFunctionBegin;
  MyTsCtx *my_ctx = static_cast<MyTsCtx *>(ctx);
  my_ctx->snes_iter_counter++;
  MoFEMFunctionReturn(0);
}
 
PetscErrorCode MyTSResizeSetup(TS ts, PetscInt nstep, PetscReal time, Vec sol,
                               PetscBool *resize, void *ctx) {
  PetscFunctionBegin;
  ResizeCtx *rctx = static_cast<ResizeCtx *>(ctx);
  *resize = rctx->mesh_changed;
  PetscFunctionReturn(0);
}
 
PetscErrorCode MyTSResizeTransfer(TS ts, PetscInt nv, Vec ts_vecsin[],
                                  Vec ts_vecsout[], void *ctx) {
  PetscFunctionBegin;
  ResizeCtx *rctx = static_cast<ResizeCtx *>(ctx);
 
  if (auto ptr = rctx->ptr.lock()) {
 
    MoFEM::Interface &m_field = *(rctx->m_field);
 
    auto bit_mng = m_field.getInterface<BitRefManager>();
 
    auto set_prev_mesh_bit = [](EntityHandle ent, BitRefLevel &bit) {
      if (bit[COMPUTATION_BIT]) {
        bit.set(COMPUTATION_BIT, false);
        bit.set(PREVIOUS_BIT, true);
      }
    };
 
    CHKERR bit_mng->lambdaBitRefLevel(set_prev_mesh_bit);
 
    std::string improved_reasons = "";
 
    Range edges_to_not_refine;
 
    CHKERR ptr->ExRawPtr->doEdgeFlips(rctx->el_q_map, rctx->flipped_els,
                                      rctx->th_spatial_coords,
                                      rctx->new_connectivity);
    if (rctx->el_q_map.size() > 0) {
      improved_reasons += ", Edge flipping";
    }
 
    bool refined = false;
    int step;
    CHKERR TSGetStepNumber(ts, &step);
    CHKERR ptr->ExRawPtr->doEdgeSplits(refined, (step) ? true : false);
    if (refined) {
      improved_reasons += ", Edge splitting";
    }
 
    // Erase leading ", " from improved_reasons
    improved_reasons.erase(0, 2);
 
    if (improved_reasons != "")
      MOFEM_LOG("REMESHING", Sev::inform)
          << "Improved mesh quality by: " << improved_reasons;
 
#ifndef NDEBUG
 
    // comp_mesh_bit = BitRefLevel().set(COMPUTATION_BIT);
    // flipped_bit = BitRefLevel().set(FLIPPED_BIT);
    // refined_bit = BitRefLevel().set(FIRST_REF_BIT);
    // prev_mesh_bit = BitRefLevel().set(PREVIOUS_BIT);
    // virgin_mesh_bit = BitRefLevel().set(VIRGIN_BIT);
    // storage_bit = BitRefLevel().set(STORAGE_BIT);
 
    if (m_field.get_comm_rank() == 0) {
      CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByDim(
          comp_mesh_bit, BitRefLevel().set(), 2,
          "before_mapping_comp_mesh_bit.vtk", "VTK", "");
      CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByDim(
          flipped_bit, BitRefLevel().set(), 2, "before_mapping_flipped_bit.vtk",
          "VTK", "");
      CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByDim(
          refined_bit, BitRefLevel().set(), 2, "before_mapping_refined_bit.vtk",
          "VTK", "");
      CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByDim(
          prev_mesh_bit, BitRefLevel().set(), 2,
          "before_mapping_prev_mesh_bit.vtk", "VTK", "");
      CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByDim(
          virgin_mesh_bit, BitRefLevel().set(), 2,
          "before_mapping_virgin_mesh_bit.vtk", "VTK", "");
    }
#endif // NDEBUG
 
    MOFEM_LOG("REMESHING", Sev::verbose) << "number of vectors to map = " << nv;
 
    for (PetscInt i = 0; i < nv; ++i) {
      double nrm;
      CHKERR VecNorm(ts_vecsin[i], NORM_2, &nrm);
      MOFEM_LOG("REMESHING", Sev::warning)
          << "Before remeshing: ts_vecsin[" << i << "] norm = " << nrm;
    }
 
    auto scatter_mng = m_field.getInterface<VecManager>();
    auto simple = m_field.getInterface<Simple>();
 
    // CHKERR m_field.getInterface<BitRefManager>()->writeBitLevel(
    //     BitRefLevel().set(1), BitRefLevel().set(), "before_resettingup.vtk",
    //     "VTK", "");
 
    // Rebuilding DM with old and new mesh entities
    auto bit = BitRefLevel().set();
    CHKERR SolutionMapping(m_field).reSetUp(
        {"TAU", "EP", "T"}, {tau_order, ep_order, T_order}, {"U"},
        {order}, {"FLUX"}, {flux_order}, bit);
 
    // CHKERR m_field.getInterface<BitRefManager>()->writeBitLevel(
    //     BitRefLevel().set(1), BitRefLevel().set(), "after_resettingup.vtk",
    //     "VTK", "");
 
    auto intermediate_dm = createDM(m_field.get_comm(), "DMMOFEM");
    CHKERR DMMoFEMCreateMoFEM(
        intermediate_dm, &m_field, "INTERMEDIATE_DM",
        prev_mesh_bit | BitRefLevel().set(FLIPPED_BIT + num_refinement_levels),
        BitRefLevel().set());
    CHKERR DMMoFEMSetDestroyProblem(intermediate_dm, PETSC_TRUE);
    CHKERR DMSetFromOptions(intermediate_dm);
    CHKERR DMMoFEMAddElement(intermediate_dm, simple->getDomainFEName());
    CHKERR DMMoFEMSetSquareProblem(intermediate_dm, PETSC_TRUE);
    CHKERR DMMoFEMSetIsPartitioned(intermediate_dm, is_distributed_mesh);
    CHKERR DMSetUp(intermediate_dm);
 
    Vec vec_in[nv], vec_out[nv];
    for (PetscInt i = 0; i < nv; ++i) {
      CHKERR DMCreateGlobalVector(intermediate_dm, &vec_in[i]);
      CHKERR VecDuplicate(vec_in[i], &vec_out[i]);
    }
 
    VecScatter scatter_to_intermediate;
 
    for (PetscInt i = 0; i < nv; ++i) {
      CHKERR scatter_mng->vecScatterCreate(
          ts_vecsin[i], simple->getProblemName(), RowColData::ROW, vec_in[i],
          "INTERMEDIATE_DM", RowColData::ROW, &scatter_to_intermediate);
      CHKERR VecScatterBegin(scatter_to_intermediate, ts_vecsin[i], vec_in[i],
                             INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecScatterEnd(scatter_to_intermediate, ts_vecsin[i], vec_in[i],
                           INSERT_VALUES, SCATTER_FORWARD);
    }
 
    CHKERR VecScatterDestroy(&scatter_to_intermediate);
 
    CHKERR SolutionMapping(m_field).mapFields(
      intermediate_dm, prev_mesh_bit,
      BitRefLevel().set(FLIPPED_BIT + num_refinement_levels), nv, vec_in,
      vec_out);
 
    // Build problem on new bit ref level
    auto ts_problem_name = simple->getProblemName();
    CHKERR m_field.modify_problem_ref_level_set_bit(
        ts_problem_name,
        BitRefLevel().set(FLIPPED_BIT + num_refinement_levels));
    CHKERR m_field.modify_problem_mask_ref_level_set_bit(ts_problem_name,
                                                         BitRefLevel().set());
    CHKERR DMMoFEMSetIsPartitioned(simple->getDM(), is_distributed_mesh);
    PetscBarrier(nullptr);
    CHKERR DMSetUp_MoFEM(simple->getDM());
 
    // Update meshsets for correct BCs
    MeshsetsManager *meshsets_manager_ptr;
    CHKERR m_field.getInterface(meshsets_manager_ptr);
    CHKERR meshsets_manager_ptr->updateAllMeshsetsByEntitiesChildren(
        BitRefLevel().set());
 
    CHKERR ptr->ExRawPtr->mechanicalBC(
        BitRefLevel().set(FLIPPED_BIT + num_refinement_levels),
        BitRefLevel().set());
    CHKERR ptr->ExRawPtr->thermalBC(
        BitRefLevel().set(FLIPPED_BIT + num_refinement_levels),
        BitRefLevel().set());
 
    VecScatter scatter_to_final;
 
    for (PetscInt i = 0; i < nv; ++i) {
      CHKERR DMCreateGlobalVector(simple->getDM(), &ts_vecsout[i]);
      CHKERR scatter_mng->vecScatterCreate(
          ts_vecsout[i], simple->getProblemName(), RowColData::ROW, vec_out[i],
          "INTERMEDIATE_DM", RowColData::ROW, &scatter_to_final);
      CHKERR VecScatterBegin(scatter_to_final, vec_out[i], ts_vecsout[i],
                             INSERT_VALUES, SCATTER_REVERSE);
      CHKERR VecScatterEnd(scatter_to_final, vec_out[i], ts_vecsout[i],
                           INSERT_VALUES, SCATTER_REVERSE);
      CHKERR VecScatterDestroy(&scatter_to_final);
    }
 
    for (PetscInt i = 0; i < nv; ++i) {
      CHKERR VecDestroy(&vec_in[i]);
      CHKERR VecDestroy(&vec_out[i]);
    }
 
    auto set_all_bit = [&](EntityHandle ent, BitRefLevel &bit) {
      for (auto l = 0; l <= FLIPPED_BIT + num_refinement_levels; ++l) {
        bit.set(STORAGE_BIT, true);
      }
      if (bit[VIRGIN_BIT]) {
        bit.set(VIRGIN_BIT, false);
      }
      if (bit[FLIPPED_BIT]) {
        bit.set(VIRGIN_BIT, true);
      }
      if (bit[PREVIOUS_BIT]) {
        bit.set(PREVIOUS_BIT, false);
      }
      if (bit[FLIPPED_BIT + num_refinement_levels]) {
        bit.set(COMPUTATION_BIT, true);
      }
    };
    CHKERR bit_mng->lambdaBitRefLevel(set_all_bit);
 
#ifndef NDEBUG
 
    // comp_mesh_bit = BitRefLevel().set(COMPUTATION_BIT);
    // flipped_bit = BitRefLevel().set(FLIPPED_BIT);
    // refined_bit = BitRefLevel().set(FIRST_REF_BIT);
    // prev_mesh_bit = BitRefLevel().set(PREVIOUS_BIT);
    // virgin_mesh_bit = BitRefLevel().set(VIRGIN_BIT);
    // storage_bit = BitRefLevel().set(STORAGE_BIT);
 
    if (m_field.get_comm_rank() == 0) {
      CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByDim(
          comp_mesh_bit, BitRefLevel().set(), 2,
          "after_mapping_comp_mesh_bit.vtk", "VTK", "");
      CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByDim(
          flipped_bit, BitRefLevel().set(), 2, "after_mapping_flipped_bit.vtk",
          "VTK", "");
      CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByDim(
          refined_bit, BitRefLevel().set(), 2, "after_mapping_refined_bit.vtk",
          "VTK", "");
      CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByDim(
          prev_mesh_bit, BitRefLevel().set(), 2,
          "after_mapping_prev_mesh_bit.vtk", "VTK", "");
      CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByDim(
          virgin_mesh_bit, BitRefLevel().set(), 2,
          "after_mapping_virgin_mesh_bit.vtk", "VTK", "");
    }
#endif // NDEBUG
 
    CHKERR m_field.modify_problem_ref_level_set_bit(ts_problem_name,
                                                    comp_mesh_bit);
    CHKERR m_field.modify_problem_mask_ref_level_set_bit(ts_problem_name,
                                                         BitRefLevel().set());
 
    CHKERR TSSetDM(ts, simple->getDM());
 
    auto B = createDMMatrix(simple->getDM());
    CHKERR TSSetIJacobian(ts, B, B, nullptr, nullptr);
  }
 
  for (PetscInt i = 0; i < nv; ++i) {
    double nrm;
    CHKERR VecNorm(ts_vecsout[i], NORM_2, &nrm);
    MOFEM_LOG("REMESHING", Sev::warning)
        << "After remeshing: ts_vecsout[" << i << "] norm = " << nrm;
  }
 
  PetscFunctionReturn(0);
}
 
Example *thermoplastic_raw_ptr = nullptr;
 
MoFEMErrorCode Example::tsSolve() {
  MoFEMFunctionBegin;
 
  Simple *simple = mField.getInterface<Simple>();
  PipelineManager *pip_mng = mField.getInterface<PipelineManager>();
  ISManager *is_manager = mField.getInterface<ISManager>();
 
  auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
 
  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 create_post_process_elements = [&]() {
    auto pp_fe = boost::make_shared<PostProcEle>(mField);
 
    auto push_vol_ops = [this](auto &pip) {
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV});
 
      ScalerFunTwoArgs unity_2_args = [&](const double, const double) {
        return 1.;
      };
      ScalerFunThreeArgs unity_3_args = [&](const double, const double,
                                            const double) { return 1.; };
 
      auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
            common_thermoelastic_ptr, common_thermoplastic_ptr] =
          createCommonThermoPlasticOps<SPACE_DIM, IT, DomainEleOp>(
              mField, "MAT_PLASTIC", "MAT_THERMAL", "MAT_THERMOELASTIC",
              "MAT_THERMOPLASTIC", pip, "U", "EP", "TAU", "T", 1., unity_2_args,
              unity_2_args, unity_3_args, Sev::inform);
 
      if (common_hencky_ptr) {
        if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
          CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
      }
 
      return std::make_tuple(common_plastic_ptr, common_hencky_ptr,
                             common_thermoplastic_ptr);
    };
 
    auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
      MoFEMFunctionBegin;
 
      auto &pip = pp_fe->getOpPtrVector();
 
      auto [common_plastic_ptr, common_hencky_ptr, common_thermoplastic_ptr] =
          p;
 
      auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
      auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
 
      auto vec_temp_ptr = boost::make_shared<VectorDouble>();
      auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
 
      auto u_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
 
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      if (is_large_strains) {
 
        pip.push_back(
 
            new OpPPMap(
 
                pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
 
                {{"ISOTROPIC_HARDENING",
                  common_plastic_ptr->getPlasticIsoHardeningPtr()},
                 {"PLASTIC_SURFACE",
                  common_plastic_ptr->getPlasticSurfacePtr()},
                 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()},
                 {"T", common_thermoplastic_ptr->getTempPtr()}},
 
                {{"U", u_ptr},
                 {"FLUX", common_thermoplastic_ptr->getHeatFluxPtr()}},
 
                {{"GRAD", common_hencky_ptr->matGradPtr},
                 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
 
                {{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
                 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()},
                 {"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
                 {"HENCKY_STRESS", common_hencky_ptr->getMatHenckyStress()}}
 
                )
 
        );
 
      } else {
 
        pip.push_back(
 
            new OpPPMap(
 
                pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
 
                {{"PLASTIC_SURFACE",
                  common_plastic_ptr->getPlasticSurfacePtr()},
                 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()},
                 {"T", common_thermoplastic_ptr->getTempPtr()}},
 
                {{"U", u_ptr},
                 {"FLUX", common_thermoplastic_ptr->getHeatFluxPtr()}},
 
                {},
 
                {{"STRAIN", common_plastic_ptr->mStrainPtr},
                 {"STRESS", common_plastic_ptr->mStressPtr},
                 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
                 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
 
                )
 
        );
      }
 
      MoFEMFunctionReturn(0);
    };
 
    // Defaults to outputting volume for 2D and skin for 3D
    PetscBool post_proc_vol = PETSC_FALSE;
    PetscBool post_proc_skin = PETSC_TRUE;
    switch (SPACE_DIM) {
    case 2:
      post_proc_vol = PETSC_TRUE;
      post_proc_skin = PETSC_FALSE;
      break;
    case 3:
      post_proc_vol = PETSC_FALSE;
      post_proc_skin = PETSC_TRUE;
      break;
    default:
      CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "SPACE_DIM neither 2 nor 3");
      break;
    }
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_vol",
                               &post_proc_vol, PETSC_NULLPTR);
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin",
                               &post_proc_skin, PETSC_NULLPTR);
 
    auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
                          &post_proc_vol]() {
      if (post_proc_vol == PETSC_FALSE)
        return boost::shared_ptr<PostProcEle>();
      auto pp_fe = boost::make_shared<PostProcEle>(mField);
      CHK_MOAB_THROW(
          push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
          "push_vol_post_proc_ops");
      return pp_fe;
    };
 
    auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
                           &post_proc_skin]() {
      if (post_proc_skin == PETSC_FALSE)
        return boost::shared_ptr<SkinPostProcEle>();
 
      auto simple = mField.getInterface<Simple>();
      auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
      auto op_side = new OpLoopSide<SideEle>(mField, simple->getDomainFEName(),
                                             SPACE_DIM, Sev::verbose);
      pp_fe->getOpPtrVector().push_back(op_side);
      CHK_MOAB_THROW(push_vol_post_proc_ops(
                         pp_fe, push_vol_ops(op_side->getOpPtrVector())),
                     "push_vol_post_proc_ops");
      return pp_fe;
    };
 
    return std::make_pair(vol_post_proc(), skin_post_proc());
  };
 
  auto scatter_create = [&](auto D, auto coeff) {
    SmartPetscObj<IS> is;
    CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
                                                   ROW, "U", coeff, coeff, is);
    int loc_size;
    CHKERR ISGetLocalSize(is, &loc_size);
    Vec v;
    CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
    VecScatter scatter;
    CHKERR VecScatterCreate(D, is, v, PETSC_NULLPTR, &scatter);
    return std::make_tuple(SmartPetscObj<Vec>(v),
                           SmartPetscObj<VecScatter>(scatter));
  };
 
  boost::shared_ptr<SetPtsData> field_eval_data;
  boost::shared_ptr<MatrixDouble> u_field_ptr;
 
  std::array<double, 3> field_eval_coords = {0., 0., 0.};
  int coords_dim = 3;
  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
                                  field_eval_coords.data(), &coords_dim,
                                  &do_eval_field);
 
  boost::shared_ptr<std::map<std::string, boost::shared_ptr<VectorDouble>>>
      scalar_field_ptrs = boost::make_shared<
          std::map<std::string, boost::shared_ptr<VectorDouble>>>();
  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
      vector_field_ptrs = boost::make_shared<
          std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
      sym_tensor_field_ptrs = boost::make_shared<
          std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
      tensor_field_ptrs = boost::make_shared<
          std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
 
  auto test_monitor_ptr = boost::make_shared<FEMethod>();
 
  if (do_eval_field && atom_test == 0) {
    field_eval_data =
        mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
    CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
        field_eval_data, simple->getDomainFEName());
 
    field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
    auto no_rule = [](int, int, int) { return -1; };
    auto field_eval_fe_ptr = field_eval_data->feMethodPtr;
    field_eval_fe_ptr->getRuleHook = no_rule;
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV});
 
    ScalerFunTwoArgs unity_2_args = [&](const double, const double) {
      return 1.;
    };
    ScalerFunThreeArgs unity_3_args = [&](const double, const double,
                                          const double) { return 1.; };
 
    auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
          common_thermoelastic_ptr, common_thermoplastic_ptr] =
        createCommonThermoPlasticOps<SPACE_DIM, IT, DomainEleOp>(
            mField, "MAT_PLASTIC", "MAT_THERMAL", "MAT_THERMOELASTIC",
            "MAT_THERMOPLASTIC", field_eval_fe_ptr->getOpPtrVector(), "U", "EP",
            "TAU", "T", 1., unity_2_args, unity_2_args, unity_3_args,
            Sev::inform);
 
    auto u_field_ptr = boost::make_shared<MatrixDouble>();
    field_eval_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_field_ptr));
    auto T_field_ptr = boost::make_shared<VectorDouble>();
    field_eval_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("T", T_field_ptr));
    auto FLUX_field_ptr = boost::make_shared<MatrixDouble>();
    field_eval_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", FLUX_field_ptr));
 
    if ((common_plastic_ptr) && (common_hencky_ptr) && (scalar_field_ptrs)) {
      if (is_large_strains) {
        scalar_field_ptrs->insert(std::make_pair("TEMPERATURE", T_field_ptr));
        scalar_field_ptrs->insert(std::make_pair(
            "PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()));
        scalar_field_ptrs->insert(std::make_pair(
            "PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()));
        vector_field_ptrs->insert(std::make_pair("U", u_field_ptr));
        vector_field_ptrs->insert(std::make_pair("FLUX", FLUX_field_ptr));
        sym_tensor_field_ptrs->insert(std::make_pair(
            "PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()));
        sym_tensor_field_ptrs->insert(std::make_pair(
            "PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()));
        sym_tensor_field_ptrs->insert(std::make_pair(
            "HENCKY_STRESS", common_hencky_ptr->getMatHenckyStress()));
        sym_tensor_field_ptrs->insert(
            std::make_pair("HENCKY_STRAIN", common_hencky_ptr->getMatLogC()));
        tensor_field_ptrs->insert(
            std::make_pair("GRAD", common_hencky_ptr->matGradPtr));
        tensor_field_ptrs->insert(std::make_pair(
            "FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()));
      } else {
        scalar_field_ptrs->insert(std::make_pair("TEMPERATURE", T_field_ptr));
        scalar_field_ptrs->insert(std::make_pair(
            "PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()));
        scalar_field_ptrs->insert(std::make_pair(
            "PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()));
        vector_field_ptrs->insert(std::make_pair("U", u_field_ptr));
        vector_field_ptrs->insert(std::make_pair("FLUX", FLUX_field_ptr));
        sym_tensor_field_ptrs->insert(
            std::make_pair("STRAIN", common_plastic_ptr->mStrainPtr));
        sym_tensor_field_ptrs->insert(
            std::make_pair("STRESS", common_plastic_ptr->mStressPtr));
        sym_tensor_field_ptrs->insert(std::make_pair(
            "PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()));
        sym_tensor_field_ptrs->insert(std::make_pair(
            "PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()));
      }
    }
  } else if (atom_test > 0) {
    field_eval_data =
        mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
 
    CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
        field_eval_data, simple->getDomainFEName());
 
    field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
    auto no_rule = [](int, int, int) { return -1; };
 
    auto field_eval_fe_ptr = field_eval_data->feMethodPtr;
    field_eval_fe_ptr->getRuleHook = no_rule;
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV});
 
    ScalerFunTwoArgs unity_2_args = [&](const double, const double) {
      return 1.;
    };
    ScalerFunThreeArgs unity_3_args = [&](const double, const double,
                                          const double) { return 1.; };
 
    auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
          common_thermoelastic_ptr, common_thermoplastic_ptr] =
        createCommonThermoPlasticOps<SPACE_DIM, IT, DomainEleOp>(
            mField, "MAT_PLASTIC", "MAT_THERMAL", "MAT_THERMOELASTIC",
            "MAT_THERMOPLASTIC", field_eval_fe_ptr->getOpPtrVector(), "U", "EP",
            "TAU", "T", 1, unity_2_args, unity_2_args, unity_3_args,
            Sev::inform);
 
    auto dispGradPtr = common_hencky_ptr->matGradPtr;
    auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
 
    auto coeff_expansion_ptr = common_thermoelastic_ptr->getCoeffExpansionPtr();
    auto ref_temp_ptr = common_thermoelastic_ptr->getRefTempPtr();
 
    field_eval_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("T", tempFieldPtr));
    field_eval_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", fluxFieldPtr));
    field_eval_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("U", dispFieldPtr));
    field_eval_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
                                                                 dispGradPtr));
    plasticMultiplierFieldPtr = common_plastic_ptr->getPlasticTauPtr();
    plasticStrainFieldPtr = common_plastic_ptr->getPlasticStrainPtr();
 
    using H = HenckyOps::HenckyIntegrators<DomainEleOp>;
 
    field_eval_fe_ptr->getOpPtrVector().push_back(
        new typename H::template OpCalculateHenckyThermoPlasticStress<SPACE_DIM,
                                                                      IT, 0>(
            "U", tempFieldPtr, common_hencky_ptr, coeff_expansion_ptr,
            ref_temp_ptr));
    if (!IS_LARGE_STRAINS) {
      field_eval_fe_ptr->getOpPtrVector().push_back(
          new OpSymmetrizeTensor<SPACE_DIM>(
              common_hencky_ptr->getMatFirstPiolaStress(), stressFieldPtr));
      field_eval_fe_ptr->getOpPtrVector().push_back(
          new OpSymmetrizeTensor<SPACE_DIM>(dispGradPtr, strainFieldPtr));
    } else {
      field_eval_fe_ptr->getOpPtrVector().push_back(
          new typename H::template OpCalculateLogC<SPACE_DIM, IT>(
              "U", common_hencky_ptr));
      stressFieldPtr = common_hencky_ptr->getMatFirstPiolaStress();
      strainFieldPtr = common_hencky_ptr->getMatLogC();
    };
  }
 
  auto set_time_monitor = [&](auto dm, auto solver) {
    MoFEMFunctionBegin;
    if (atom_test == 0) {
      boost::shared_ptr<Monitor<SPACE_DIM>> monitor_ptr(new Monitor<SPACE_DIM>(
          dm, create_post_process_elements(), reactionFe, uXScatter, uYScatter,
          uZScatter, field_eval_coords, field_eval_data, scalar_field_ptrs,
          vector_field_ptrs, sym_tensor_field_ptrs, tensor_field_ptrs));
      boost::shared_ptr<ForcesAndSourcesCore> null;
      CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
                                 monitor_ptr, null, null);
    } else {
 
      test_monitor_ptr->preProcessHook = [&]() {
        MoFEMFunctionBegin;
 
        CHKERR mField.getInterface<FieldEvaluatorInterface>()
            ->evalFEAtThePoint<SPACE_DIM>(
                field_eval_coords.data(), 1e-12, simple->getProblemName(),
                simple->getDomainFEName(), field_eval_data,
                mField.get_comm_rank(), mField.get_comm_rank(), nullptr,
                MF_EXIST, QUIET);
 
        auto eval_num_vec = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
        CHKERR VecZeroEntries(eval_num_vec);
        if (tempFieldPtr->size()) {
          CHKERR VecSetValue(eval_num_vec, 0, 1, ADD_VALUES);
        }
        CHKERR VecAssemblyBegin(eval_num_vec);
        CHKERR VecAssemblyEnd(eval_num_vec);
 
        double eval_num;
        CHKERR VecSum(eval_num_vec, &eval_num);
        if (!(int)eval_num) {
          SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                  "atom test %d failed: did not find elements to evaluate "
                  "the field, check the coordinates",
                  atom_test);
        }
 
        if (tempFieldPtr->size()) {
          auto t_temp = getFTensor0FromVec(*tempFieldPtr);
          MOFEM_LOG("THERMOPLASTICITY", Sev::inform)
              << "Eval point T magnitude: " << t_temp;
          if (atom_test && fabs(test_monitor_ptr->ts_t - 10) < 1e-12) {
            if (atom_test <= 3 && fabs(t_temp - 554.48) > 1e-2) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong temperature value",
                      atom_test);
            }
            if (atom_test == 4 && fabs(t_temp - 250) > 1e-2) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong temperature value",
                      atom_test);
            }
            if (atom_test == 5 && fabs(t_temp - 1) > 1e-2) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong temperature value",
                      atom_test);
            }
          }
          if (atom_test == 8 && fabs(t_temp - 0.5) > 1e-12) {
            SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                    "atom test %d failed: wrong temperature value", atom_test);
          }
        }
        if (fluxFieldPtr->size1()) {
          FTensor::Index<'i', SPACE_DIM> i;
          auto t_flux = getFTensor1FromMat<SPACE_DIM>(*fluxFieldPtr);
          auto flux_mag = sqrt(t_flux(i) * t_flux(i));
          MOFEM_LOG("THERMOPLASTICITY", Sev::inform)
              << "Eval point FLUX magnitude: " << flux_mag;
          if (atom_test && fabs(test_monitor_ptr->ts_t - 10) < 1e-12) {
            if (atom_test <= 3 && fabs(flux_mag - 27008.0) > 2e1) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong flux value", atom_test);
            }
            if (atom_test == 4 && fabs(flux_mag - 150e3) > 1e-1) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong flux value", atom_test);
            }
            if (atom_test == 5 && fabs(flux_mag) > 1e-6) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong flux value", atom_test);
            }
          }
        }
        if (dispFieldPtr->size1()) {
          FTensor::Index<'i', SPACE_DIM> i;
          auto t_disp = getFTensor1FromMat<SPACE_DIM>(*dispFieldPtr);
          auto disp_mag = sqrt(t_disp(i) * t_disp(i));
          MOFEM_LOG("THERMOPLASTICITY", Sev::inform)
              << "Eval point U magnitude: " << disp_mag;
          if (atom_test && fabs(test_monitor_ptr->ts_t - 10) < 1e-12) {
            if (atom_test == 1 && fabs(disp_mag - 0.00345) > 1e-5) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong displacement value",
                      atom_test);
            }
            if ((atom_test == 2 || atom_test == 3) &&
                fabs(disp_mag - 0.00265) > 1e-5) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong displacement value",
                      atom_test);
            }
            if ((atom_test == 5) &&
                fabs(t_disp(0) - (std::sqrt(std::exp(0.2)) - 1)) > 1e-5 &&
                fabs(t_disp(1) - (std::sqrt(std::exp(0.2)) - 1)) > 1e-5) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong displacement value",
                      atom_test);
            }
          }
        }
        if (strainFieldPtr->size1()) {
          FTensor::Index<'i', SPACE_DIM> i;
          auto t_strain =
              getFTensor2SymmetricFromMat<SPACE_DIM>(*strainFieldPtr);
          auto t_strain_trace = t_strain(i, i);
          if (atom_test && fabs(test_monitor_ptr->ts_t - 10) < 1e-12) {
            if (atom_test == 1 && fabs(t_strain_trace - 0.00679) > 1e-5) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong strain value", atom_test);
            }
            if ((atom_test == 2 || atom_test == 3) &&
                fabs(t_strain_trace - 0.00522) > 1e-5) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong strain value", atom_test);
            }
            if ((atom_test == 5) && fabs(t_strain_trace - 0.2) > 1e-5) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong strain value", atom_test);
            }
          }
        }
        if (stressFieldPtr->size1()) {
          double von_mises_stress;
          auto getVonMisesStress = [&](auto t_stress) {
            MoFEMFunctionBegin;
            von_mises_stress = std::sqrt(
                0.5 *
                ((t_stress(0, 0) - t_stress(1, 1)) *
                     (t_stress(0, 0) - t_stress(1, 1)) +
                 (SPACE_DIM == 3 ? (t_stress(1, 1) - t_stress(2, 2)) *
                                       (t_stress(1, 1) - t_stress(2, 2))
                                 : 0) +
                 (SPACE_DIM == 3 ? (t_stress(2, 2) - t_stress(0, 0)) *
                                       (t_stress(2, 2) - t_stress(0, 0))
                                 : 0) +
                 6 * (t_stress(0, 1) * t_stress(0, 1) +
                      (SPACE_DIM == 3 ? t_stress(1, 2) * t_stress(1, 2) : 0) +
                      (SPACE_DIM == 3 ? t_stress(2, 0) * t_stress(2, 0) : 0))));
            MoFEMFunctionReturn(0);
          };
          if (!IS_LARGE_STRAINS) {
            auto t_stress =
                getFTensor2SymmetricFromMat<SPACE_DIM>(*stressFieldPtr);
            CHKERR getVonMisesStress(t_stress);
          } else {
            auto t_stress =
                getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*stressFieldPtr);
            CHKERR getVonMisesStress(t_stress);
          }
          MOFEM_LOG("THERMOPLASTICITY", Sev::inform)
              << "Eval point von Mises Stress: " << von_mises_stress;
          if (atom_test && fabs(test_monitor_ptr->ts_t - 10) < 1e-12) {
            if (atom_test == 1 && fabs(von_mises_stress - 523.0) > 5e-1) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong von Mises stress value",
                      atom_test);
            }
            if (atom_test == 2 && fabs(von_mises_stress - 16.3) > 5e-2) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong von Mises stress value",
                      atom_test);
            }
            if (atom_test == 3 && fabs(von_mises_stress - 14.9) > 5e-2) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong von Mises stress value",
                      atom_test);
            }
            if (atom_test == 5 && fabs(von_mises_stress) > 5e-2) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong von Mises stress value",
                      atom_test);
            }
          }
        }
        if (plasticMultiplierFieldPtr->size()) {
          auto t_plastic_multiplier =
              getFTensor0FromVec(*plasticMultiplierFieldPtr);
          MOFEM_LOG("THERMOPLASTICITY", Sev::inform)
              << "Eval point TAU magnitude: " << t_plastic_multiplier;
          if (atom_test == 8 && fabs(t_plastic_multiplier - 1.5) > 1e-12) {
            SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                    "atom test %d failed: wrong plastic multiplier value",
                    atom_test);
          }
        }
        if (plasticStrainFieldPtr->size1()) {
          FTENSOR_INDEX(SPACE_DIM, i);
          FTENSOR_INDEX(SPACE_DIM, j);
          auto t_plastic_strain =
              getFTensor2SymmetricFromMat<SPACE_DIM>(*plasticStrainFieldPtr);
          MOFEM_LOG("THERMOPLASTICITY", Sev::inform)
              << "Eval point EP(0,0) = " << t_plastic_strain(0, 0)
              << ", EP(0,1) = " << t_plastic_strain(0, 1)
              << ", EP(1,1) = " << t_plastic_strain(1, 1);
          if (atom_test == 8) {
            if (fabs(t_plastic_strain(0, 0) - 0.005) > 1e-12) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong EP(0,0) value", atom_test);
            }
            if (fabs(t_plastic_strain(0, 1) - 0.025) > 1e-12) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong EP(0,1) value", atom_test);
            }
            if (fabs(t_plastic_strain(1, 1) - 0.015) > 1e-12) {
              SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                      "atom test %d failed: wrong EP(1,1) value", atom_test);
            }
          }
        }
 
        MOFEM_LOG_SYNCHRONISE(mField.get_comm());
        MoFEMFunctionReturn(0);
      };
      auto null = boost::shared_ptr<FEMethod>();
      CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(), null,
                                 test_monitor_ptr, null);
    };
 
    MoFEMFunctionReturn(0);
  };
 
  auto set_schur_pc = [&](auto solver,
                          boost::shared_ptr<SetUpSchur> &schur_ptr) {
    MoFEMFunctionBeginHot;
 
    auto name_prb = simple->getProblemName();
 
    // create sub dm for Schur complement
    auto create_schur_dm = [&](SmartPetscObj<DM> base_dm,
                               SmartPetscObj<DM> &dm_sub) {
      MoFEMFunctionBegin;
      dm_sub = createDM(mField.get_comm(), "DMMOFEM");
      CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SCHUR");
      CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
      CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
      CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
      for (auto f : {"U", "FLUX"}) {
        CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
        CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
      }
      CHKERR DMSetUp(dm_sub);
 
      MoFEMFunctionReturn(0);
    };
 
    auto create_block_dm = [&](SmartPetscObj<DM> base_dm,
                               SmartPetscObj<DM> &dm_sub) {
      MoFEMFunctionBegin;
      dm_sub = createDM(mField.get_comm(), "DMMOFEM");
      CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "BLOCK");
      CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
      CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
      CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
#ifdef ADD_CONTACT
      for (auto f : {"SIGMA", "EP", "TAU", "T"}) {
        CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
        CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
      }
#else
      for (auto f : {"EP", "TAU", "T"}) {
        CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
        CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
      }
#endif
      CHKERR DMSetUp(dm_sub);
      MoFEMFunctionReturn(0);
    };
 
    // Create nested (sub BC) Schur DM
    if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
 
      SmartPetscObj<DM> dm_schur;
      CHKERR create_schur_dm(simple->getDM(), dm_schur);
      SmartPetscObj<DM> dm_block;
      CHKERR create_block_dm(simple->getDM(), dm_block);
 
#ifdef ADD_CONTACT
 
      auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
        auto block_mat_data = createBlockMatStructure(
            simple->getDM(),
 
            {
 
                {simple->getDomainFEName(),
 
                 {{"U", "U"},
                  {"SIGMA", "SIGMA"},
                  {"U", "SIGMA"},
                  {"SIGMA", "U"},
                  {"EP", "EP"},
                  {"TAU", "TAU"},
                  {"U", "EP"},
                  {"EP", "U"},
                  {"EP", "TAU"},
                  {"TAU", "EP"},
                  {"TAU", "U"},
                  {"T", "T"},
                  {"U", "T"},
                  {"T", "U"},
                  {"EP", "T"},
                  {"T", "EP"},
                  {"TAU", "T"},
                  {"T", "TAU"}
 
                 }},
 
                {simple->getBoundaryFEName(),
 
                 {{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
 
                 }}
 
            }
 
        );
 
        return createSchurNestedMatrixStruture(
 
            {dm_schur, dm_block}, block_mat_data,
 
            {"SIGMA", "EP", "TAU", "T"}, {nullptr, nullptr, nullptr, nullptr},
            true
 
        );
      };
 
#else
 
      auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
        auto block_mat_data =
            createBlockMatStructure(simple->getDM(),
 
                                    {{simple->getDomainFEName(),
 
                                      {{"U", "U"},
                                       {"EP", "EP"},
                                       {"TAU", "TAU"},
                                       {"U", "EP"},
                                       {"EP", "U"},
                                       {"EP", "TAU"},
                                       {"TAU", "U"},
                                       {"TAU", "EP"},
                                       {"T", "T"},
                                       {"T", "U"},
                                       {"T", "FLUX"},
                                       {"T", "EP"},
                                       {"FLUX", "T"},
                                       {"FLUX", "U"},
                                       {"U", "T"},
                                       {"EP", "T"},
                                       {"TAU", "T"}
 
                                      }}}
 
            );
 
        return createSchurNestedMatrixStruture(
 
            {dm_schur, dm_block}, block_mat_data,
 
            {"EP", "TAU", "T"}, {nullptr, nullptr, nullptr}, false
 
        );
      };
 
#endif
 
      auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
      CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
 
      auto block_is = getDMSubData(dm_block)->getSmartRowIs();
      auto ao_schur = getDMSubData(dm_schur)->getSmartRowMap();
 
      // Indices has to be map fro very to level, while assembling Schur
      // complement.
      schur_ptr =
          SetUpSchur::createSetUpSchur(mField, dm_schur, block_is, ao_schur);
      CHKERR schur_ptr->setUp(solver);
    }
 
    MoFEMFunctionReturnHot(0);
  };
 
  auto dm = simple->getDM();
  auto D = createDMVector(dm);
  auto DD = vectorDuplicate(D);
  uXScatter = scatter_create(D, 0);
  uYScatter = scatter_create(D, 1);
  if constexpr (SPACE_DIM == 3)
    uZScatter = scatter_create(D, 2);
 
  auto create_solver = [pip_mng]() {
    if (is_quasi_static == PETSC_TRUE)
      return pip_mng->createTSIM();
    else
      return pip_mng->createTSIM2();
  };
 
  SmartPetscObj<Vec> solution_vector;
  CHKERR DMCreateGlobalVector_MoFEM(dm, solution_vector);
 
  if (restart_flg) {
    PetscViewer viewer;
    CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, restart_file, FILE_MODE_READ,
                                 &viewer);
    CHKERR VecLoad(solution_vector, viewer);
    CHKERR PetscViewerDestroy(&viewer);
    CHKERR DMoFEMMeshToLocalVector(dm, solution_vector, INSERT_VALUES,
                                   SCATTER_REVERSE);
  }
 
  auto solver = create_solver();
 
  auto active_pre_lhs = []() {
    MoFEMFunctionBegin;
    std::fill(PlasticOps::CommonData::activityData.begin(),
              PlasticOps::CommonData::activityData.end(), 0);
    MoFEMFunctionReturn(0);
  };
 
  auto active_post_lhs = [&]() {
    MoFEMFunctionBegin;
    auto get_iter = [&]() {
      SNES snes;
      CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
      int iter;
      CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
                        "Can not get iter");
      return iter;
    };
 
    auto iter = get_iter();
    if (iter >= 0) {
 
      std::array<int, 5> activity_data;
      std::fill(activity_data.begin(), activity_data.end(), 0);
      MPI_Allreduce(PlasticOps::CommonData::activityData.data(),
                    activity_data.data(), activity_data.size(), MPI_INT,
                    MPI_SUM, mField.get_comm());
 
      int &active_points = activity_data[0];
      int &avtive_full_elems = activity_data[1];
      int &avtive_elems = activity_data[2];
      int &nb_points = activity_data[3];
      int &nb_elements = activity_data[4];
 
      if (nb_points) {
 
        double proc_nb_points =
            100 * static_cast<double>(active_points) / nb_points;
        double proc_nb_active =
            100 * static_cast<double>(avtive_elems) / nb_elements;
        double proc_nb_full_active = 100;
        if (avtive_elems)
          proc_nb_full_active =
              100 * static_cast<double>(avtive_full_elems) / avtive_elems;
 
        MOFEM_LOG_C("PLASTICITY", Sev::inform,
                    "Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active "
                    "elements %d "
                    "(%3.3f\%) nb full active elems %d (%3.3f\%)",
                    iter, nb_points, active_points, proc_nb_points,
                    avtive_elems, proc_nb_active, avtive_full_elems,
                    proc_nb_full_active, iter);
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_active_dofs_elem = [&](auto dm) {
    MoFEMFunctionBegin;
    auto fe_pre_proc = boost::make_shared<FEMethod>();
    fe_pre_proc->preProcessHook = active_pre_lhs;
    auto fe_post_proc = boost::make_shared<FEMethod>();
    fe_post_proc->postProcessHook = active_post_lhs;
    auto ts_ctx_ptr = getDMTsCtx(dm);
    ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
    ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
    MoFEMFunctionReturn(0);
  };
 
  auto set_essential_bc = [&](auto dm, auto solver) {
    MoFEMFunctionBegin;
    // This is low level pushing finite elements (pipelines) to solver
 
    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>();
 
    // Add boundary condition scaling
    auto disp_time_scale = boost::make_shared<TimeScale>();
 
    auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
      EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
                                                     {disp_time_scale}, false);
      return hook;
    };
    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]() {
      return EssentialPostProcLhs<DisplacementCubitBcData>(
          mField, post_proc_lhs_ptr, 1.);
    };
    post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
 
    auto ts_ctx_ptr = getDMTsCtx(dm);
    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);
    post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
    ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
 
  if (is_quasi_static == PETSC_TRUE) {
    CHKERR TSSetSolution(solver, D);
  } else {
    CHKERR TS2SetSolution(solver, D, DD);
  }
 
  auto set_up_adapt = [&](auto solver) {
    MoFEMFunctionBegin;
    CHKERR TSAdaptRegister(TSADAPTMOFEM, TSAdaptCreateMoFEM);
    TSAdapt adapt;
    CHKERR TSGetAdapt(solver, &adapt);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_section_monitor(solver);
  CHKERR set_time_monitor(dm, solver);
  CHKERR set_up_adapt(solver);
  CHKERR TSSetFromOptions(solver);
 
  thermoplastic_raw_ptr = this;
 
  if (restart_flg) {
    CHKERR TSSetTime(solver, restart_time);
    CHKERR TSSetStepNumber(solver, restart_step);
  }
 
  CHKERR add_active_dofs_elem(dm);
  boost::shared_ptr<SetUpSchur> schur_ptr;
  CHKERR set_schur_pc(solver, schur_ptr);
  CHKERR set_essential_bc(dm, solver);
 
  auto my_ctx = boost::make_shared<MyTsCtx>();
 
  ts_to_ctx[solver] = my_ctx.get();
 
  SNES snes;
  CHKERR TSGetSNES(solver, &snes);
 
  CHKERR SNESMonitorSet(snes, SNESIterationMonitor, my_ctx.get(), NULL);
 
  CHKERR TSSetPreStage(solver, TSIterationPreStage);
 
  auto my_rhs = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx) {
    MyTsCtx *my_ts_ctx = static_cast<MyTsCtx *>(ctx);
    auto ts_ctx = my_ts_ctx->dmts_ctx;
    auto &m_field = ts_ctx->mField;
    auto simple = m_field.getInterface<Simple>();
 
    // Get the iteration number of the snes solver
    SNES snes;
    CHK_THROW_MESSAGE(TSGetSNES(ts, &snes), "Cannot get SNES");
 
    double time_increment;
    CHK_THROW_MESSAGE(TSGetTimeStep(ts, &time_increment), "Cannot get iter");
    if (time_increment < min_dt) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW,
              "Minimum time increment exceeded");
    }
 
    if (my_ts_ctx->snes_iter_counter == 0) {
      int error = system("rm ./out_snes_residual_iter_*.h5m > /dev/null 2>&1");
    }
 
    // Get the time step
    double ts_dt;
    CHK_THROW_MESSAGE(TSGetTimeStep(ts, &ts_dt),
                      "Cannot get timestep from TS object");
 
    if (ts_dt < min_dt) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW,
              "Minimum timestep exceeded");
    }
 
    // Post process the residuals
    auto post_proc = [&](auto dm, auto f_res, auto sol, auto sol_dot,
                         auto out_name) {
      MoFEMFunctionBegin;
 
      auto smart_f_res = SmartPetscObj<Vec>(f_res, true);
      auto smart_sol = SmartPetscObj<Vec>(sol, true);
      auto smart_sol_dot = SmartPetscObj<Vec>(sol_dot, true);
 
      auto post_proc_fe =
          boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
              ts_ctx->mField);
      auto &pip = post_proc_fe->getOpPtrVector();
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV});
 
      // pointers for residuals
      auto disp_res_mat = boost::make_shared<MatrixDouble>();
      auto tau_res_vec = boost::make_shared<VectorDouble>();
      auto plastic_strain_res_mat = boost::make_shared<MatrixDouble>();
      auto flux_res_mat = boost::make_shared<MatrixDouble>();
      auto temp_res_vec = boost::make_shared<VectorDouble>();
 
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
          "U", disp_res_mat, smart_f_res));
      pip.push_back(
          new OpCalculateScalarFieldValues("TAU", tau_res_vec, smart_f_res));
      pip.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
          "EP", plastic_strain_res_mat, smart_f_res));
      // pip.push_back(
      //     new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", flux_res_mat,
      //                                                 smart_f_res));
      pip.push_back(
          new OpCalculateScalarFieldValues("T", temp_res_vec, smart_f_res));
 
      // pointers for fields
      auto disp_mat = boost::make_shared<MatrixDouble>();
      auto tau_vec = boost::make_shared<VectorDouble>();
      auto plastic_strain_mat = boost::make_shared<MatrixDouble>();
      auto flux_mat = boost::make_shared<MatrixDouble>();
      auto temp_vec = boost::make_shared<VectorDouble>();
 
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", disp_mat));
      pip.push_back(new OpCalculateScalarFieldValues("TAU", tau_vec));
      pip.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
          "EP", plastic_strain_mat));
      // pip.push_back(
      //     new OpCalculateVectorFieldValues<SPACE_DIM>("FLUX", flux_mat,
      //                                                 smart_sol));
      pip.push_back(new OpCalculateScalarFieldValues("T", temp_vec));
 
      // pointers for rates of fields
      auto disp_dot_mat = boost::make_shared<MatrixDouble>();
      auto tau_dot_vec = boost::make_shared<VectorDouble>();
      auto plastic_strain_dot_mat = boost::make_shared<MatrixDouble>();
      auto flux_dot_mat = boost::make_shared<MatrixDouble>();
      auto temp_dot_vec = boost::make_shared<VectorDouble>();
 
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
          "U", disp_dot_mat, smart_sol_dot));
      pip.push_back(
          new OpCalculateScalarFieldValues("TAU", tau_dot_vec, smart_sol_dot));
      pip.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
          "EP", plastic_strain_dot_mat, smart_sol_dot));
      // pip.push_back(
      //     new OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", flux_dot_mat,
      //                                                 smart_sol_dot));
      pip.push_back(
          new OpCalculateScalarFieldValues("T", temp_dot_vec, smart_sol_dot));
 
      constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
      auto make_d_mat = [&]() {
        return boost::make_shared<MatrixDouble>(size_symm * size_symm, 1);
      };
 
      auto push_vol_ops = [&](auto &pip) {
        ScalerFunTwoArgs unity_2_args = [&](const double, const double) {
          return 1.;
        };
        ScalerFunThreeArgs unity_3_args = [&](const double, const double,
                                              const double) { return 1.; };
 
        auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
              common_thermoelastic_ptr, common_thermoplastic_ptr] =
            thermoplastic_raw_ptr
                ->createCommonThermoPlasticOps<SPACE_DIM, IT, DomainEleOp>(
                    m_field, "MAT_PLASTIC", "MAT_THERMAL", "MAT_THERMOELASTIC",
                    "MAT_THERMOPLASTIC", pip, "U", "EP", "TAU", "T", 1.,
                    unity_2_args, unity_2_args, unity_3_args, Sev::inform);
 
        if (common_hencky_ptr) {
          if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
            CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                              "Wrong pointer for grad");
        }
 
        return std::make_tuple(common_plastic_ptr, common_hencky_ptr,
                               common_thermoplastic_ptr);
      };
 
      auto push_vol_post_proc_ops = [&](auto &pp_fe, auto &&p) {
        MoFEMFunctionBegin;
 
        auto &pip = pp_fe->getOpPtrVector();
 
        auto [common_plastic_ptr, common_hencky_ptr, common_thermoplastic_ptr] =
            p;
 
        auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
        auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
 
        auto vec_temp_ptr = boost::make_shared<VectorDouble>();
        auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
 
        auto u_ptr = boost::make_shared<MatrixDouble>();
        pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
 
        using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
        if (is_large_strains) {
 
          pip.push_back(
 
              new OpPPMap(
 
                  pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
 
                  {{"RESIDUAL_TAU", tau_res_vec},
                   {"RESIDUAL_T", temp_res_vec},
                   {"TAU", tau_vec},
                   {"T", temp_vec},
                   {"RATE_TAU", tau_dot_vec},
                   {"RATE_T", temp_dot_vec},
                   {"ISOTROPIC_HARDENING",
                    common_plastic_ptr->getPlasticIsoHardeningPtr()},
                   {"PLASTIC_SURFACE",
                    common_plastic_ptr->getPlasticSurfacePtr()},
                   {"PLASTIC_MULTIPLIER",
                    common_plastic_ptr->getPlasticTauPtr()},
                   {"YIELD_FUNCTION",
                    common_plastic_ptr->getPlasticSurfacePtr()},
                   {"T", common_thermoplastic_ptr->getTempPtr()}},
 
                  {{"RESIDUAL_U", disp_res_mat},
                   {"RATE_U", disp_dot_mat},
                   {"U", u_ptr},
                   {"FLUX", common_thermoplastic_ptr->getHeatFluxPtr()}},
 
                  {{"GRAD", common_hencky_ptr->matGradPtr},
                   {"FIRST_PIOLA",
                    common_hencky_ptr->getMatFirstPiolaStress()}},
 
                  {{"RESIDUAL_EP", plastic_strain_res_mat},
                   {"RATE_EP", plastic_strain_dot_mat},
                   {"PLASTIC_STRAIN",
                    common_plastic_ptr->getPlasticStrainPtr()},
                   {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()},
                   {"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
                   {"HENCKY_STRESS", common_hencky_ptr->getMatHenckyStress()}}
 
                  )
 
          );
 
        } else {
 
          pip.push_back(
 
              new OpPPMap(
 
                  pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
 
                  {{"RESIDUAL_TAU", tau_res_vec},
                   {"RESIDUAL_T", temp_res_vec},
                   {"TAU", tau_vec},
                   {"T", temp_vec},
                   {"RATE_TAU", tau_dot_vec},
                   {"RATE_T", temp_dot_vec},
                   {"ISOTROPIC_HARDENING",
                    common_plastic_ptr->getPlasticIsoHardeningPtr()},
                   {"PLASTIC_SURFACE",
                    common_plastic_ptr->getPlasticSurfacePtr()},
                   {"PLASTIC_MULTIPLIER",
                    common_plastic_ptr->getPlasticTauPtr()},
                   {"YIELD_FUNCTION",
                    common_plastic_ptr->getPlasticSurfacePtr()},
                   {"T", common_thermoplastic_ptr->getTempPtr()}},
 
                  {{"RESIDUAL_U", disp_res_mat},
                   //  {"RESIDUAL_FLUX", flux_res_mat},
                   {"U", disp_mat},
                   //  {"FLUX", flux_mat},
                   {"RATE_U", disp_dot_mat},
                   {"U", u_ptr},
                   {"FLUX", common_thermoplastic_ptr->getHeatFluxPtr()}},
 
                  {},
 
                  {{"RESIDUAL_PLASTIC_STRAIN", plastic_strain_res_mat},
                   {"PLASTIC_STRAIN", plastic_strain_mat},
                   {"RATE_PLASTIC_STRAIN", plastic_strain_dot_mat},
                   {"STRAIN", common_plastic_ptr->mStrainPtr},
                   {"STRESS", common_plastic_ptr->mStressPtr},
                   {"PLASTIC_STRAIN",
                    common_plastic_ptr->getPlasticStrainPtr()},
                   {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
 
                  )
 
          );
        }
 
        MoFEMFunctionReturn(0);
      };
 
      CHK_MOAB_THROW(push_vol_post_proc_ops(post_proc_fe, push_vol_ops(pip)),
                     "push_vol_post_proc_ops");
 
      CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                      post_proc_fe);
 
      post_proc_fe->writeFile(out_name);
      MoFEMFunctionReturn(0);
    };
 
    // Output the RHS
 
    auto set_RHS = TsSetIFunction(ts, t, u, u_t, F, ts_ctx.get());
    CHKERR post_proc(simple->getDM(),
                     //
                     F, u, u_t,
                     //
                     "out_snes_residual_iter_" +
                         std::to_string(my_ts_ctx->snes_iter_counter) + ".h5m");
    return set_RHS;
  };
 
  PetscBool is_output_residual_fields = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-output_residual_fields",
                             &is_output_residual_fields, PETSC_NULLPTR);
 
  if (is_output_residual_fields == PETSC_TRUE) {
    my_ctx->dmts_ctx = getDMTsCtx(simple->getDM());
    // auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
    CHKERR TSSetIFunction(solver, PETSC_NULLPTR, my_rhs, my_ctx.get());
  }
 
  CHKERR TSSetDM(solver, dm);
 
  auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
  tsPrePostProc = ts_pre_post_proc;
 
  if (auto ptr = tsPrePostProc.lock()) {
    ptr->ExRawPtr = this;
    std::multimap<double, EntityHandle> no_el_q_map;
    Range no_flipped_els;
    std::vector<EntityHandle> no_new_connectivity;
    Tag no_th_spatial_coords;
 
    ResizeCtx *ctx = new ResizeCtx{
        &mField,        PETSC_FALSE,         no_el_q_map,
        no_flipped_els, no_new_connectivity, no_th_spatial_coords};
    PetscBool rollback =
        PETSC_TRUE; // choose true if you want to restart current step
    CHKERR TSSetResize(solver, rollback, MyTSResizeSetup, MyTSResizeTransfer,
                       ctx);
    CHKERR TSSetApplicationContext(solver, (void *)ctx);
    // ptr->globSol = createDMVector(dm);
    // CHKERR DMoFEMMeshToLocalVector(dm, ptr->globSol, INSERT_VALUES,
    //                                SCATTER_FORWARD);
    // CHKERR VecAssemblyBegin(ptr->globSol);
    // CHKERR VecAssemblyEnd(ptr->globSol);
    MOFEM_LOG_CHANNEL("TIMER");
    MOFEM_LOG_TAG("TIMER", "timer");
    if (set_timer)
      BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
    MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";
    CHKERR ptr->tsSetUp(solver);
    CHKERR TSSetUp(solver);
    MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";
    MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";
    CHKERR TSSolve(solver, PETSC_NULLPTR);
    MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";
  }
 
  MoFEMFunctionReturn(0);
}
//! [Solve]
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
#ifdef ADD_CONTACT
  #ifdef ENABLE_PYTHON_BINDING
  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(), "PLASTICITY"));
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "THERMAL"));
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "THERMOPLASTICITY"));
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "REMESHING"));
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
  LogManager::setLog("PLASTICITY");
  LogManager::setLog("THERMAL");
  LogManager::setLog("THERMOPLASTICITY");
  LogManager::setLog("REMESHING");
  MOFEM_LOG_TAG("PLASTICITY", "Plasticity");
  MOFEM_LOG_TAG("THERMAL", "Thermal");
  MOFEM_LOG_TAG("THERMOPLASTICITY", "Thermoplasticity");
  MOFEM_LOG_TAG("REMESHING", "Remeshing");
 
#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]
 
    //! [Load mesh]
    Simple *simple = m_field.getInterface<Simple>();
    CHKERR simple->getOptions();
 
    simple->getBitRefLevel() = comp_mesh_bit;
    simple->getBitRefLevelMask() = BitRefLevel().set();
 
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_distributed_mesh",
                               &is_distributed_mesh, PETSC_NULLPTR);
    if (is_distributed_mesh == PETSC_TRUE) {
      CHKERR simple->loadFile();
    } else {
      char meshFileName[255];
      CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-file_name",
                                   meshFileName, 255, PETSC_NULLPTR);
      CHKERR simple->loadFile("", meshFileName);
    }
 
    Range ents;
    CHKERR m_field.get_moab().get_entities_by_dimension(0, SPACE_DIM, ents,
                                                        true);
    CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
        ents, comp_mesh_bit | prev_mesh_bit | virgin_mesh_bit | storage_bit,
        false);
    //! [Load mesh]
 
    //! [Example]
    Example ex(m_field);
    CHKERR ex.runProblem();
    //! [Example]
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
 
#ifdef ADD_CONTACT
  #ifdef ENABLE_PYTHON_BINDING
  if (Py_FinalizeEx() < 0) {
    exit(120);
  }
  #endif
#endif // ADD_CONTACT
 
  return 0;
}
 
struct SetUpSchurImpl : public SetUpSchur {
 
  SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj<DM> sub_dm,
                 SmartPetscObj<IS> field_split_is, SmartPetscObj<AO> ao_up)
      : SetUpSchur(), mField(m_field), subDM(sub_dm),
        fieldSplitIS(field_split_is), aoSchur(ao_up) {
    if (S) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Is expected that schur matrix is not "
                        "allocated. This is "
                        "possible only is PC is set up twice");
    }
  }
  virtual ~SetUpSchurImpl() { S.reset(); }
 
  MoFEMErrorCode setUp(TS solver);
  MoFEMErrorCode preProc();
  MoFEMErrorCode postProc();
 
private:
  SmartPetscObj<Mat> S;
 
  MoFEM::Interface &mField;
  SmartPetscObj<DM> subDM;        ///< field split sub dm
  SmartPetscObj<IS> fieldSplitIS; ///< IS for split Schur block
  SmartPetscObj<AO> aoSchur;      ///> main DM to subDM
};
 
MoFEMErrorCode SetUpSchurImpl::setUp(TS solver) {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pip_mng = mField.getInterface<PipelineManager>();
 
  SNES snes;
  CHKERR TSGetSNES(solver, &snes);
  KSP ksp;
  CHKERR SNESGetKSP(snes, &ksp);
  CHKERR KSPSetFromOptions(ksp);
 
  PC pc;
  CHKERR KSPGetPC(ksp, &pc);
  PetscBool is_pcfs = PETSC_FALSE;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
  if (is_pcfs) {
    if (S) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Is expected that schur matrix is not "
                        "allocated. This is "
                        "possible only is PC is set up twice");
    }
 
    S = createDMMatrix(subDM);
    CHKERR MatSetBlockSize(S, SPACE_DIM);
    
 
    // Set DM to use shell block matrix
    DM solver_dm;
    CHKERR TSGetDM(solver, &solver_dm);
    CHKERR DMSetMatType(solver_dm, MATSHELL);
 
    auto ts_ctx_ptr = getDMTsCtx(solver_dm);
    auto A = createDMBlockMat(simple->getDM());
    auto P = createDMNestSchurMat(simple->getDM());
 
    if (is_quasi_static == PETSC_TRUE) {
      auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
                              Mat A, Mat B, void *ctx) {
        return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
      };
      CHKERR TSSetIJacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
    } else {
      auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
                              PetscReal a, PetscReal aa, Mat A, Mat B,
                              void *ctx) {
        return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
      };
      CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
    }
    CHKERR KSPSetOperators(ksp, A, P);
 
    auto set_ops = [&]() {
      MoFEMFunctionBegin;
      auto pip_mng = mField.getInterface<PipelineManager>();
 
#ifndef ADD_CONTACT
      // Boundary
      pip_mng->getOpBoundaryLhsPipeline().push_front(
          createOpSchurAssembleBegin());
      pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
 
          {"EP", "TAU", "T"}, {nullptr, nullptr, nullptr}, aoSchur, S, false,
          false
 
          ));
      // Domain
      pip_mng->getOpDomainLhsPipeline().push_front(
          createOpSchurAssembleBegin());
      pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
 
          {"EP", "TAU", "T"}, {nullptr, nullptr, nullptr}, aoSchur, S, false,
          false
 
          ));
#else
 
      double eps_stab = 1e-4;
      CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-eps_stab", &eps_stab,
                                   PETSC_NULLPTR);
 
      using B = FormsIntegrators<BoundaryEleOp>::Assembly<
          BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<IT>;
      using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
 
      // Boundary
      pip_mng->getOpBoundaryLhsPipeline().push_front(
          createOpSchurAssembleBegin());
      pip_mng->getOpBoundaryLhsPipeline().push_back(
          new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
            return eps_stab;
          }));
      pip_mng->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd(
 
          {"SIGMA", "EP", "TAU", "T"}, {nullptr, nullptr, nullptr, nullptr},
          aoSchur, S, false, false
 
          ));
      // Domain
      pip_mng->getOpDomainLhsPipeline().push_front(
          createOpSchurAssembleBegin());
      pip_mng->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd(
 
          {"SIGMA", "EP", "TAU", "T"}, {nullptr, nullptr, nullptr, nullptr},
          aoSchur, S, false, false
 
          ));
#endif // ADD_CONTACT
      MoFEMFunctionReturn(0);
    };
 
    auto set_assemble_elems = [&]() {
      MoFEMFunctionBegin;
      auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
      schur_asmb_pre_proc->preProcessHook = [this]() {
        MoFEMFunctionBegin;
        CHKERR MatZeroEntries(S);
        MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
        MoFEMFunctionReturn(0);
      };
      auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
 
      schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
        MoFEMFunctionBegin;
        MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
 
        // Apply essential constrains to Schur complement
        CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
        CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
        CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
            mField, schur_asmb_post_proc, 1, S, aoSchur)();
 
        MoFEMFunctionReturn(0);
      };
      auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
      ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
      ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
      MoFEMFunctionReturn(0);
    };
 
    auto set_pc = [&]() {
      MoFEMFunctionBegin;
      CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
      CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
      MoFEMFunctionReturn(0);
    };
 
    auto set_diagonal_pc = [&]() {
      MoFEMFunctionBegin;
      KSP *subksp;
      CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
      auto get_pc = [](auto ksp) {
        PC pc_raw;
        CHKERR KSPGetPC(ksp, &pc_raw);
        return SmartPetscObj<PC>(pc_raw,
                                 true); // bump reference
      };
      CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
      CHKERR PetscFree(subksp);
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_ops();
    CHKERR set_pc();
    CHKERR set_assemble_elems();
 
    CHKERR TSSetUp(solver);
    CHKERR KSPSetUp(ksp);
    CHKERR set_diagonal_pc();
 
  } else {
    pip_mng->getOpBoundaryLhsPipeline().push_front(
        createOpSchurAssembleBegin());
    pip_mng->getOpBoundaryLhsPipeline().push_back(
        createOpSchurAssembleEnd({}, {}));
    pip_mng->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
    pip_mng->getOpDomainLhsPipeline().push_back(
        createOpSchurAssembleEnd({}, {}));
  }
 
  // fieldSplitIS.reset();
  // aoSchur.reset();
  MoFEMFunctionReturn(0);
}
 
boost::shared_ptr<SetUpSchur>
SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
                             SmartPetscObj<DM> sub_dm, SmartPetscObj<IS> is_sub,
                             SmartPetscObj<AO> ao_up) {
  return boost::shared_ptr<SetUpSchur>(
      new SetUpSchurImpl(m_field, sub_dm, is_sub, ao_up));
}