ThermoPlasticOps.hpp

Copyright

Copyright (c) 2024

/**
 * \file ThermoPlasticOps.hpp
 * \example ThermoPlasticOps.hpp
 *
 * @copyright Copyright (c) 2024
 */
 
#ifndef __THERMO_PLASTIC_OPS_HPP__
#define __THERMO_PLASTIC_OPS_HPP__
 
namespace ThermoPlasticOps {
 
struct BlockedParameters
    : public boost::enable_shared_from_this<BlockedParameters> {
  double coeffExpansion;
  double refTemp;
  double heatConductivity;
  double heatCapacity;
  double inelasticHeatFraction;
 
  inline auto getCoeffExpansionPtr();
  inline auto getRefTempPtr();
  inline auto getHeatConductivityPtr();
  inline auto getHeatCapacityPtr();
  inline auto getInelasticHeatFractionPtr();
};
 
struct ThermoPlasticBlockedParameters
    : public boost::enable_shared_from_this<ThermoPlasticBlockedParameters> {
  double omega0;
  double omegaH;
  double inelasticHeatFraction;
  double temp0;
  VectorDouble temperature;
  MatrixDouble heat_flux;
 
  VectorDouble resCdTemperature;
  MatrixDouble resFlowDtemp;
 
  inline auto getOmega0Ptr() {
    return boost::shared_ptr<double>(shared_from_this(), &omega0);
  }
 
  inline auto getOmegaHPtr() {
    return boost::shared_ptr<double>(shared_from_this(), &omegaH);
  }
 
  inline auto getInelasticHeatFractionPtr() {
    return boost::shared_ptr<double>(shared_from_this(),
                                     &inelasticHeatFraction);
  }
 
  inline auto getTemp0Ptr() {
    return boost::shared_ptr<double>(shared_from_this(), &temp0);
  }
 
  inline auto getTempPtr() {
    return boost::shared_ptr<VectorDouble>(shared_from_this(), &temperature);
  }
 
  inline auto getHeatFluxPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &heat_flux);
  }
};
 
MoFEMErrorCode addMatThermoPlasticBlockOps(
    MoFEM::Interface &m_field,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
    std::string thermoplastic_block_name,
    boost::shared_ptr<ThermoPlasticBlockedParameters> blockedParamsPtr,
    boost::shared_ptr<ThermoElasticOps::BlockedThermalParameters>
        &blockedThermalParamsPtr,
    Sev sev, ScalerFunThreeArgs inelastic_heat_fraction_scaling);
 
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);
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingRhsImpl;
 
// Calculate deltaT \zeta \bm{T} : \dot{\bm{E}}_p
// Calculate deltaT \zeta T_ij : \dot{E_ij}_p
template <int DIM, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingRhsImpl<DIM, GAUSS, AssemblyDomainEleOp>
    : public AssemblyDomainEleOp {
 
  OpCalculateAdiabaticHeatingRhsImpl(
      const std::string field_name,
      boost::shared_ptr<MatrixDouble> stress_measure,
      boost::shared_ptr<MatrixDouble> plastic_strain_rate_measure,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
        stressMeasurePtr(stress_measure),
        plasticStrainRateMeasurePtr(plastic_strain_rate_measure),
        inelasticHeatingFunction(inelastic_heating_function) {}
 
protected:
  boost::shared_ptr<MatrixDouble> stressMeasurePtr;
  boost::shared_ptr<MatrixDouble> plasticStrainRateMeasurePtr;
  ScalarFun inelasticHeatingFunction;
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
};
 
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculateAdiabaticHeatingRhsImpl<DIM, GAUSS, AssemblyDomainEleOp>::iNtegrate(
    EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
 
  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
 
  // get stress and strain values
  auto t_stress = getFTensor2SymmetricFromMat<DIM>(*stressMeasurePtr);
  auto t_plastic_strain_rate =
      getFTensor2SymmetricFromMat<DIM>(*plasticStrainRateMeasurePtr);
 
  // get element volume
  const double vol = AssemblyDomainEleOp::getMeasure();
  // get integration weights
  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
  // get base function gradient on rows
  auto t_row_base = data.getFTensor0N();
 
  // get coordinate at integration points
  auto t_coords = AssemblyDomainEleOp::getFTensor1CoordsAtGaussPts();
 
  // loop over integration points
  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
    // take into account Jacobian
    const double alpha =
        t_w * vol *
        inelasticHeatingFunction(t_coords(0), t_coords(1), t_coords(2)) *
        t_stress(i, j) * t_plastic_strain_rate(i, j);
 
    // loop over rows base functions
    size_t rr = 0;
    for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
      AssemblyDomainEleOp::locF[rr] += alpha * t_row_base;
      ++t_row_base;
    }
    for (; rr < AssemblyDomainEleOp::nbRowBaseFunctions; ++rr)
      ++t_row_base;
    ++t_w; // move to another integration weight
    ++t_coords;
    ++t_stress;
    ++t_plastic_strain_rate;
  };
  MoFEMFunctionReturn(0);
}
 
// Templated on IS_LARGE_STRAINS to allow different implementations at compile
// time
template <int SPACE_DIM, IntegrationType I, typename AssemblyDomainEleOp,
          bool IS_LARGE_STRAINS>
struct OpCalculateAdiabaticHeatingRhs;
 
template <int SPACE_DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingRhs<SPACE_DIM, I, AssemblyDomainEleOp, false>
    : public OpCalculateAdiabaticHeatingRhsImpl<SPACE_DIM, GAUSS,
                                                AssemblyDomainEleOp> {
  OpCalculateAdiabaticHeatingRhs(
      const std::string field_name,
      boost::shared_ptr<MatrixDouble> stress_measure,
      boost::shared_ptr<MatrixDouble> plastic_strain_rate_measure,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpCalculateAdiabaticHeatingRhsImpl<SPACE_DIM, GAUSS,
                                           AssemblyDomainEleOp>(
            field_name, stress_measure, plastic_strain_rate_measure,
            inelastic_heating_function, ents_ptr) {}
};
 
template <int SPACE_DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingRhs<SPACE_DIM, I, AssemblyDomainEleOp, true>
    : public OpCalculateAdiabaticHeatingRhsImpl<SPACE_DIM, GAUSS,
                                                AssemblyDomainEleOp> {
  OpCalculateAdiabaticHeatingRhs(
      const std::string field_name,
      boost::shared_ptr<MatrixDouble> stress_measure,
      boost::shared_ptr<MatrixDouble> plastic_strain_rate_measure,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpCalculateAdiabaticHeatingRhsImpl<SPACE_DIM, GAUSS,
                                           AssemblyDomainEleOp>(
            field_name, stress_measure, plastic_strain_rate_measure,
            inelastic_heating_function, ents_ptr) {}
};
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdTImpl;
 
template <int DIM, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdTImpl<DIM, GAUSS, AssemblyDomainEleOp>
    : public AssemblyDomainEleOp {
  OpCalculateAdiabaticHeatingLhsdTImpl(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<HenckyOps::CommonData> elastic_common_data_ptr,
      boost::shared_ptr<PlasticOps::CommonData> plastic_common_data_ptr,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : AssemblyDomainEleOp(row_field_name, col_field_name,
                            AssemblyDomainEleOp::OPROWCOL),
        elasticCommonDataPtr(elastic_common_data_ptr),
        plasticCommonDataPtr(plastic_common_data_ptr),
        inelasticHeatingFunction(inelastic_heating_function),
        coeffExpansionPtr(coeff_expansion_ptr), entsPtr(ents_ptr) {}
 
protected:
  boost::shared_ptr<HenckyOps::CommonData> elasticCommonDataPtr;
  boost::shared_ptr<PlasticOps::CommonData> plasticCommonDataPtr;
  ScalarFun inelasticHeatingFunction;
  boost::shared_ptr<VectorDouble> coeffExpansionPtr;
  boost::shared_ptr<Range> entsPtr;
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data);
};
 
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculateAdiabaticHeatingLhsdTImpl<DIM, GAUSS, AssemblyDomainEleOp>::
    iNtegrate(EntitiesFieldData::EntData &row_data,
              EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
 
  const auto nb_integration_pts = row_data.getN().size1();
  const auto nb_row_base_functions = row_data.getN().size2();
  auto t_w = this->getFTensor0IntegrationWeight();
 
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
  auto t_row_base = row_data.getFTensor0N();
 
  // get plastic strain rate values
  auto t_plastic_strain_rate = getFTensor2SymmetricFromMat<DIM>(
      *(plasticCommonDataPtr->getPlasticStrainDotPtr()));
 
  // get elasticity tensor
  auto t_D =
      getFTensor4DdgFromMat<DIM, DIM, 0>(*(elasticCommonDataPtr->matDPtr));
 
  // get coordinate at integration points
  auto t_coords = AssemblyDomainEleOp::getFTensor1CoordsAtGaussPts();
 
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  FTensor::Index<'k', DIM> k;
  FTensor::Index<'l', DIM> l;
 
  auto t_coeff_exp = FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM>();
  t_coeff_exp(i, j) = 0;
  for (auto d = 0; d != SPACE_DIM; ++d) {
    t_coeff_exp(d, d) = (*coeffExpansionPtr)[d];
  }
 
  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
 
    double alpha = this->getMeasure() * t_w;
    auto rr = 0;
    for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
      auto t_col_base = col_data.getFTensor0N(gg, 0);
      for (auto cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
#ifdef HENCKY_SMALL_STRAIN
        AssemblyDomainEleOp::locMat(rr, cc) -=
            alpha * t_row_base *
            inelasticHeatingFunction(t_coords(0), t_coords(1), t_coords(2)) *
            t_coeff_exp(i, j) * t_D(i, j, k, l) * t_plastic_strain_rate(k, l) *
            t_col_base;
#else
        AssemblyDomainEleOp::locMat(rr, cc) -=
            alpha * t_row_base *
            inelasticHeatingFunction(t_coords(0), t_coords(1), t_coords(2)) *
            t_coeff_exp(i, j) * t_D(i, j, k, l) * t_plastic_strain_rate(k, l) *
            t_col_base;
#endif
        ++t_col_base;
      }
 
      ++t_row_base;
    }
    for (; rr != nb_row_base_functions; ++rr)
      ++t_row_base;
 
    ++t_w;
    ++t_plastic_strain_rate;
    ++t_D;
    ++t_coords;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp,
          bool IS_LARGE_STRAINS>
struct OpCalculateAdiabaticHeatingLhsdT;
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdT<DIM, I, AssemblyDomainEleOp, false>
    : public OpCalculateAdiabaticHeatingLhsdTImpl<DIM, GAUSS,
                                                  AssemblyDomainEleOp> {
  OpCalculateAdiabaticHeatingLhsdT(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<HenckyOps::CommonData> elastic_common_data_ptr,
      boost::shared_ptr<PlasticOps::CommonData> plastic_common_data_ptr,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpCalculateAdiabaticHeatingLhsdTImpl<DIM, GAUSS, AssemblyDomainEleOp>(
            row_field_name, col_field_name, elastic_common_data_ptr,
            plastic_common_data_ptr, inelastic_heating_function,
            coeff_expansion_ptr, ents_ptr) {}
};
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdT<DIM, I, AssemblyDomainEleOp, true>
    : public OpCalculateAdiabaticHeatingLhsdTImpl<DIM, GAUSS,
                                                  AssemblyDomainEleOp> {
  OpCalculateAdiabaticHeatingLhsdT(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<HenckyOps::CommonData> elastic_common_data_ptr,
      boost::shared_ptr<PlasticOps::CommonData> plastic_common_data_ptr,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<VectorDouble> coeff_expansion_ptr,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpCalculateAdiabaticHeatingLhsdTImpl<DIM, GAUSS, AssemblyDomainEleOp>(
            row_field_name, col_field_name, elastic_common_data_ptr,
            plastic_common_data_ptr, inelastic_heating_function,
            coeff_expansion_ptr, ents_ptr) {}
};
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdEPImpl;
 
template <int DIM, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdEPImpl<DIM, GAUSS, AssemblyDomainEleOp>
    : public AssemblyDomainEleOp {
  OpCalculateAdiabaticHeatingLhsdEPImpl(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<HenckyOps::CommonData> elastic_common_data_ptr,
      boost::shared_ptr<PlasticOps::CommonData> plastic_common_data_ptr,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : AssemblyDomainEleOp(row_field_name, col_field_name,
                            AssemblyDomainEleOp::OPROWCOL),
        elasticCommonDataPtr(elastic_common_data_ptr),
        plasticCommonDataPtr(plastic_common_data_ptr),
        inelasticHeatingFunction(inelastic_heating_function),
        entsPtr(ents_ptr) {
    this->sYmm = false;
  }
 
protected:
  boost::shared_ptr<HenckyOps::CommonData> elasticCommonDataPtr;
  boost::shared_ptr<PlasticOps::CommonData> plasticCommonDataPtr;
  ScalarFun inelasticHeatingFunction;
  boost::shared_ptr<Range> entsPtr;
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data);
};
 
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculateAdiabaticHeatingLhsdEPImpl<DIM, GAUSS, AssemblyDomainEleOp>::
    iNtegrate(EntitiesFieldData::EntData &row_data,
              EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  FTensor::Index<'k', DIM> k;
  FTensor::Index<'l', DIM> l;
  FTensor::Index<'m', DIM> m;
  FTensor::Index<'n', DIM> n;
  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
  FTensor::Index<'L', size_symm> L;
 
  const auto nb_integration_pts = row_data.getN().size1();
  const auto nb_row_base_functions = row_data.getN().size2();
  auto t_w = this->getFTensor0IntegrationWeight();
 
  auto t_L = PlasticOps::symm_L_tensor(FTensor::Number<DIM>());
 
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
  auto t_row_base = row_data.getFTensor0N();
 
  auto t_stress = getFTensor2SymmetricFromMat<DIM>(
      *(elasticCommonDataPtr->getMatHenckyStress()));
 
  // get plastic strain rate values
  auto t_plastic_strain_rate = getFTensor2SymmetricFromMat<DIM>(
      *(plasticCommonDataPtr->getPlasticStrainDotPtr()));
 
  // get elasticity tensor
  auto t_D =
      getFTensor4DdgFromMat<DIM, DIM, 0>(*(elasticCommonDataPtr->matDPtr));
 
  // get coordinate at integration points
  auto t_coords = AssemblyDomainEleOp::getFTensor1CoordsAtGaussPts();
 
  auto t_diff = PlasticOps::diff_tensor(FTensor::Number<DIM>());
 
  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
    // get the local vector to store the symmetric tensor of size_symm
    auto t_vec = getFTensor1FromPtr<size_symm>(&this->locMat(0, 0));
    double alpha = this->getMeasure() * t_w;
    auto rr = 0;
    for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
      auto t_col_base = col_data.getFTensor0N(gg, 0);
      for (auto cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; cc++) {
#ifdef HENCKY_SMALL_STRAIN
        t_vec(L) -=
            alpha * t_row_base *
            inelasticHeatingFunction(t_coords(0), t_coords(1), t_coords(2)) *
            (t_plastic_strain_rate(i, j) * t_D(i, j, k, l) *
                 t_diff(k, l, m, n) -
             this->getTSa() * t_stress(i, j) * t_diff(i, j, m, n)) *
            t_col_base * t_L(m, n, L);
#else
        t_vec(L) -=
            alpha * t_row_base *
            inelasticHeatingFunction(t_coords(0), t_coords(1), t_coords(2)) *
            (t_plastic_strain_rate(i, j) * t_D(i, j, k, l) *
                 t_diff(k, l, m, n) -
             this->getTSa() * t_stress(i, j) * t_diff(i, j, m, n)) *
            t_col_base * t_L(m, n, L);
#endif
        ++t_vec;
        ++t_col_base;
      }
 
      ++t_row_base;
    }
    for (; rr != nb_row_base_functions; ++rr)
      ++t_row_base;
 
    ++t_w;
    ++t_plastic_strain_rate;
    ++t_D;
    ++t_coords;
    ++t_stress;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp,
          bool IS_LARGE_STRAINS>
struct OpCalculateAdiabaticHeatingLhsdEP;
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdEP<DIM, I, AssemblyDomainEleOp, false>
    : public OpCalculateAdiabaticHeatingLhsdEPImpl<DIM, GAUSS,
                                                   AssemblyDomainEleOp> {
  OpCalculateAdiabaticHeatingLhsdEP(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<HenckyOps::CommonData> elastic_common_data_ptr,
      boost::shared_ptr<PlasticOps::CommonData> plastic_common_data_ptr,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpCalculateAdiabaticHeatingLhsdEPImpl<DIM, GAUSS, AssemblyDomainEleOp>(
            row_field_name, col_field_name, elastic_common_data_ptr,
            plastic_common_data_ptr, inelastic_heating_function, ents_ptr) {}
};
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdEP<DIM, I, AssemblyDomainEleOp, true>
    : public OpCalculateAdiabaticHeatingLhsdEPImpl<DIM, GAUSS,
                                                   AssemblyDomainEleOp> {
  OpCalculateAdiabaticHeatingLhsdEP(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<HenckyOps::CommonData> elastic_common_data_ptr,
      boost::shared_ptr<PlasticOps::CommonData> plastic_common_data_ptr,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpCalculateAdiabaticHeatingLhsdEPImpl<DIM, GAUSS, AssemblyDomainEleOp>(
            row_field_name, col_field_name, elastic_common_data_ptr,
            plastic_common_data_ptr, inelastic_heating_function, ents_ptr) {}
};
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdUImpl;
 
template <int DIM, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdUImpl<DIM, GAUSS, AssemblyDomainEleOp>
    : public AssemblyDomainEleOp {
  OpCalculateAdiabaticHeatingLhsdUImpl(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<HenckyOps::CommonData> elastic_common_data_ptr,
      boost::shared_ptr<PlasticOps::CommonData> plastic_common_data_ptr,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : AssemblyDomainEleOp(row_field_name, col_field_name,
                            AssemblyDomainEleOp::OPROWCOL),
        elasticCommonDataPtr(elastic_common_data_ptr),
        plasticCommonDataPtr(plastic_common_data_ptr),
        inelasticHeatingFunction(inelastic_heating_function),
        entsPtr(ents_ptr) {
    this->sYmm = false;
  }
 
protected:
  boost::shared_ptr<HenckyOps::CommonData> elasticCommonDataPtr;
  boost::shared_ptr<PlasticOps::CommonData> plasticCommonDataPtr;
  ScalarFun inelasticHeatingFunction;
  boost::shared_ptr<Range> entsPtr;
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data);
};
 
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculateAdiabaticHeatingLhsdUImpl<DIM, GAUSS, AssemblyDomainEleOp>::
    iNtegrate(EntitiesFieldData::EntData &row_data,
              EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  FTensor::Index<'k', DIM> k;
  FTensor::Index<'l', DIM> l;
  FTensor::Index<'m', DIM> m;
  FTensor::Index<'n', DIM> n;
  FTensor::Index<'o', DIM> o;
  FTensor::Index<'p', DIM> p;
 
  const auto nb_integration_pts = row_data.getN().size1();
  const auto nb_row_base_functions = row_data.getN().size2();
 
  auto t_w = this->getFTensor0IntegrationWeight();
 
  auto t_row_base = row_data.getFTensor0N();
 
  auto t_logC_dC =
      getFTensor4DdgFromMat<DIM, DIM>(elasticCommonDataPtr->matLogCdC);
  auto t_C_dF =
      getFTensor4FromMat<DIM, DIM, DIM, DIM, 1>(elasticCommonDataPtr->matCdF);
 
  // get plastic strain rate values
  auto t_plastic_strain_rate = getFTensor2SymmetricFromMat<DIM>(
      *(plasticCommonDataPtr->getPlasticStrainDotPtr()));
 
  // get elasticity tensor
  auto t_D =
      getFTensor4DdgFromMat<DIM, DIM, 0>(*(elasticCommonDataPtr->matDPtr));
 
  // get coordinate at integration points
  auto t_coords = AssemblyDomainEleOp::getFTensor1CoordsAtGaussPts();
 
  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
    auto t_vec = getFTensor1FromPtr<DIM>(&this->locMat(0, 0));
    double alpha = this->getMeasure() * t_w;
 
    auto rr = 0;
    for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
      auto t_col_grad_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
      for (auto cc = 0; cc != AssemblyDomainEleOp::nbCols / DIM; cc++) {
#ifdef HENCKY_SMALL_STRAIN
        // TODO: implement correct small strain problem
        t_vec(m) +=
            0.5 * alpha * t_row_base *
            inelasticHeatingFunction(t_coords(0), t_coords(1), t_coords(2)) *
            (t_plastic_strain_rate(i, j) * t_D(i, j, k, l) *
             t_logC_dC(k, l, o, p) * t_C_dF(o, p, m, n)) *
            t_col_grad_base(n);
#else
        // Multiply by a half for dEdC from dlogCdC
        t_vec(m) +=
            0.5 * alpha * t_row_base *
            inelasticHeatingFunction(t_coords(0), t_coords(1), t_coords(2)) *
            (t_plastic_strain_rate(i, j) * t_D(i, j, k, l) *
             t_logC_dC(k, l, o, p) * t_C_dF(o, p, m, n)) *
            t_col_grad_base(n);
#endif
        ++t_col_grad_base;
        ++t_vec;
      }
 
      ++t_row_base;
    }
    for (; rr != nb_row_base_functions; ++rr)
      ++t_row_base;
 
    ++t_w;
    ++t_plastic_strain_rate;
    ++t_D;
    ++t_coords;
    ++t_logC_dC;
    ++t_C_dF;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp,
          bool IS_LARGE_STRAINS>
struct OpCalculateAdiabaticHeatingLhsdU;
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdU<DIM, I, AssemblyDomainEleOp, false>
    : public OpCalculateAdiabaticHeatingLhsdUImpl<DIM, GAUSS,
                                                  AssemblyDomainEleOp> {
  OpCalculateAdiabaticHeatingLhsdU(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<HenckyOps::CommonData> elastic_common_data_ptr,
      boost::shared_ptr<PlasticOps::CommonData> plastic_common_data_ptr,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpCalculateAdiabaticHeatingLhsdUImpl<DIM, GAUSS, AssemblyDomainEleOp>(
            row_field_name, col_field_name, elastic_common_data_ptr,
            plastic_common_data_ptr, inelastic_heating_function, ents_ptr) {}
};
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateAdiabaticHeatingLhsdU<DIM, I, AssemblyDomainEleOp, true>
    : public OpCalculateAdiabaticHeatingLhsdUImpl<DIM, GAUSS,
                                                  AssemblyDomainEleOp> {
  OpCalculateAdiabaticHeatingLhsdU(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<HenckyOps::CommonData> elastic_common_data_ptr,
      boost::shared_ptr<PlasticOps::CommonData> plastic_common_data_ptr,
      ScalarFun inelastic_heating_function,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpCalculateAdiabaticHeatingLhsdUImpl<DIM, GAUSS, AssemblyDomainEleOp>(
            row_field_name, col_field_name, elastic_common_data_ptr,
            plastic_common_data_ptr, inelastic_heating_function, ents_ptr) {}
};
 
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode opThermoPlasticFactoryDomainRhs(
    MoFEM::Interface &m_field, std::string block_name,
    std::string thermal_block_name,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string u, std::string ep, std::string tau, std::string temperature);
 
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode opThermoPlasticFactoryDomainLhs(
    MoFEM::Interface &m_field, std::string block_name,
    std::string thermal_block_name,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string u, std::string ep, std::string tau, std::string temperature);
 
template <IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_dTImpl;
 
template <IntegrationType I>
using OpCalculateConstraintsLhs_dT =
    OpCalculateConstraintsLhs_dTImpl<I, AssemblyDomainEleOp>;
 
template <typename AssemblyDomainEleOp>
struct OpCalculateConstraintsLhs_dTImpl<GAUSS, AssemblyDomainEleOp>
    : public AssemblyDomainEleOp {
  OpCalculateConstraintsLhs_dTImpl(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<ThermoPlasticOps::ThermoPlasticBlockedParameters>
          TP_common_data_ptr);
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data);
 
private:
  boost::shared_ptr<ThermoPlasticOps::ThermoPlasticBlockedParameters>
      TPcommonDataPtr;
};
 
template <typename AssemblyDomainEleOp>
OpCalculateConstraintsLhs_dTImpl<GAUSS, AssemblyDomainEleOp>::
    OpCalculateConstraintsLhs_dTImpl(
        const std::string row_field_name, const std::string col_field_name,
        boost::shared_ptr<ThermoPlasticOps::ThermoPlasticBlockedParameters>
            TP_common_data_ptr)
    : AssemblyDomainEleOp(row_field_name, col_field_name,
                          DomainEleOp::OPROWCOL),
      TPcommonDataPtr(TP_common_data_ptr) {
  AssemblyDomainEleOp::sYmm = false;
}
 
template <typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculateConstraintsLhs_dTImpl<GAUSS, AssemblyDomainEleOp>::iNtegrate(
    EntitiesFieldData::EntData &row_data,
    EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
 
  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
  const auto nb_row_base_functions = row_data.getN().size2();
 
  auto t_res_c_dtemperature =
      getFTensor0FromVec(TPcommonDataPtr->resCdTemperature);
  auto next = [&]() { ++t_res_c_dtemperature; };
 
  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
  auto t_row_base = row_data.getFTensor0N();
  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
    const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
    ++t_w;
 
    const auto res = alpha * (t_res_c_dtemperature);
    next();
 
    auto mat_ptr = AssemblyDomainEleOp::locMat.data().begin();
    size_t rr = 0;
    for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
      auto t_col_base = col_data.getFTensor0N(gg, 0);
      for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; ++cc) {
        *mat_ptr += t_row_base * t_col_base * res;
        ++t_col_base;
        ++mat_ptr;
      }
      ++t_row_base;
    }
    for (; rr < nb_row_base_functions; ++rr)
      ++t_row_base;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_dTImpl;
 
template <int DIM, IntegrationType I>
using OpCalculatePlasticFlowLhs_dT =
    OpCalculatePlasticFlowLhs_dTImpl<DIM, I, AssemblyDomainEleOp>;
 
template <int DIM, typename AssemblyDomainEleOp>
struct OpCalculatePlasticFlowLhs_dTImpl<DIM, GAUSS, AssemblyDomainEleOp>
    : public AssemblyDomainEleOp {
  OpCalculatePlasticFlowLhs_dTImpl(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<ThermoPlasticOps::ThermoPlasticBlockedParameters>
          common_TP_data_ptr);
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data);
 
private:
  boost::shared_ptr<ThermoPlasticOps::ThermoPlasticBlockedParameters>
      commonTPDataPtr;
};
 
template <int DIM, typename AssemblyDomainEleOp>
OpCalculatePlasticFlowLhs_dTImpl<DIM, GAUSS, AssemblyDomainEleOp>::
    OpCalculatePlasticFlowLhs_dTImpl(
        const std::string row_field_name, const std::string col_field_name,
        boost::shared_ptr<ThermoPlasticOps::ThermoPlasticBlockedParameters>
            common_TP_data_ptr)
    : AssemblyDomainEleOp(row_field_name, col_field_name,
                          DomainEleOp::OPROWCOL),
      commonTPDataPtr(common_TP_data_ptr) {
  AssemblyDomainEleOp::sYmm = false;
}
 
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode
OpCalculatePlasticFlowLhs_dTImpl<DIM, GAUSS, AssemblyDomainEleOp>::iNtegrate(
    EntitiesFieldData::EntData &row_data,
    EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
  FTensor::Index<'L', size_symm> L;
 
  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
  const size_t nb_row_base_functions = row_data.getN().size2();
  auto &locMat = AssemblyDomainEleOp::locMat;
 
  auto t_res_flow_dtemp =
      getFTensor2SymmetricFromMat<DIM>(commonTPDataPtr->resFlowDtemp);
 
  auto t_L = PlasticOps::symm_L_tensor(FTensor::Number<DIM>());
 
  auto next = [&]() { ++t_res_flow_dtemp; };
 
  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
  auto t_row_base = row_data.getFTensor0N();
  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
    double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
    ++t_w;
    FTensor::Tensor1<double, size_symm> t_res_vec;
    t_res_vec(L) = alpha * (t_res_flow_dtemp(i, j) * t_L(i, j, L));
    next();
 
    size_t rr = 0;
    for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
      auto t_mat = PlasticOps::get_mat_tensor_sym_dscalar(
          rr, locMat, FTensor::Number<DIM>());
      auto t_col_base = col_data.getFTensor0N(gg, 0);
      for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
        t_mat(L) += t_row_base * t_col_base * t_res_vec(L);
        ++t_mat;
        ++t_col_base;
      }
      ++t_row_base;
    }
    for (; rr != nb_row_base_functions; ++rr)
      ++t_row_base;
  }
 
  MoFEMFunctionReturn(0);
}
 
}; // namespace ThermoPlasticOps
 
#endif // __FINITE_THERMAL_OPS_HPP__