ADOLCPlasticityMaterialModels.hpp

/** \file ADOLCPlasticityMaterialModels.hpp
 * \ingroup adoc_plasticity
 * \brief Matetial models for plasticity
 *
 * \example ADOLCPlasticityMaterialModels.hpp
 *
 */
 
#ifndef __ADOLC_MATERIAL_MODELS_HPP
#define __ADOLC_MATERIAL_MODELS_HPP
 
#ifndef WITH_ADOL_C
#error "MoFEM need to be compiled with ADOL-C"
#endif
 
namespace ADOLCPlasticity {
 
/** \brief J2 plasticity (Kinematic Isotropic (Linear) Hardening)
 * \ingroup adoc_plasticity
 */
 
/**
 * @brief J2 plasticity (Kinematic Isotropic (Linear) Hardening)
 *
 * @tparam DIM model dimension (2 - is plane stress)
 */
template <int DIM> struct J2Plasticity;
 
/**
 * @brief J2 (Von Misses) plasticity
 * @ingroup adoc_plasticity
 *
 * Free energy:
 *
 *  \f[
 * \psi(\varepsilon^e, \alpha, \beta) =
 * \frac{1}{2} \lambda \textrm{tr}[\varepsilon^e]^2
 * + \mu \varepsilon^e : \varepsilon^e
 * +
 * \sigma_y\alpha
 * +
 * \frac{1}{2} \phi H \alpha^2
 * +
 * \frac{1}{2} (1-\phi)K \beta^2
 * \f]
 *
 * Isotropic hardening variable \f$\alpha\f$ is in first index of
 * internal variables vector. Kinematic hardening variable is in
 * the remaining indices of internal variables vector.
 *
 * Yield function:
 *
 * Deviator of actual stress:
 * \f[
 * \eta=\textrm{dev}[\sigma]-\overline{\beta}
 * \f]
 * where \f$\overline{\beta}\f$ is back stress.
 *
 * Square of the deviator norm
 * \f[
 * f = \frac{3}{2} \eta:\eta
 * \f]
 *
 * Yield function:
 *  \f[
 * y = \sqrt{f} - \overline{\alpha}
 * \f]
 * where \f$f\f$ is defined in \ref eva
 *
 * This is associated potential model, so flow potential is equal to yield
 *
 */
template <> struct J2Plasticity<3> : public ClosestPointProjection {
 
  J2Plasticity() : ClosestPointProjection() {
    nbInternalVariables = 7;
    activeVariablesW.resize(12 + nbInternalVariables, false);
  }
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'Z', 6> Z;
 
  FTensor::Tensor2_symmetric<adouble, 3> tTotalStrain;
  FTensor::Tensor2_symmetric<adouble, 3> tPlasticStrain;
  FTensor::Tensor2_symmetric<adouble, 3> tElasticStrain;
  FTensor::Tensor2_symmetric<adouble, 3> tInternalTensor;
  FTensor::Tensor2_symmetric<adouble, 3> tStress;
  FTensor::Tensor2_symmetric<adouble, 3> tBackStress;
  FTensor::Tensor2_symmetric<adouble, 3> tDeviator;
 
  adouble tR;
 
  double mu;     //< Lamé parameter
  double lambda; //< Lamé parameters
  double H;      //< Isotropic hardening
  double K;      //< Kinematic hardening
  double phi;    //< Combined isotropic/kinematic hardening
  double sigmaY; //< Yield stress
 
  enum ModelParams { LAMBDA, MU, ISOH, KINK, PHI, SIGMA_Y, LAST_PARAM };
 
  using ParamsArray = std::array<double, LAST_PARAM>;
  boost::shared_ptr<ParamsArray> defaultParamsArrayPtr = nullptr;
  boost::shared_ptr<ParamsArray> paramsArrayPtr = nullptr;
  boost::shared_ptr<ParamsArray> oldParamsArrayPtr = nullptr;
 
  MoFEMErrorCode getDefaultMaterialParameters() {
    MoFEMFunctionBegin;
 
    double young_modulus = 1;
    double poisson_ratio = 0.25;
    sigmaY = 1;
    H = 0.1;
    K = 0;
    phi = 1;
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-young_modulus", &young_modulus,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-poisson_ratio", &poisson_ratio,
                                 0);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-sigma_y", &sigmaY, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-H", &H, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-K", &K, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-phi", &phi, PETSC_NULLPTR);
 
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "Young modulus: " << young_modulus;
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "Poisson ratio: " << poisson_ratio;
 
    mu = MU(young_modulus, poisson_ratio);
    lambda = LAMBDA(young_modulus, poisson_ratio);
 
    defaultParamsArrayPtr = boost::make_shared<ParamsArray>();
    (*defaultParamsArrayPtr)[LAMBDA] = lambda;
    (*defaultParamsArrayPtr)[MU] = mu;
    (*defaultParamsArrayPtr)[ISOH] = H;
    (*defaultParamsArrayPtr)[KINK] = K;
    (*defaultParamsArrayPtr)[PHI] = phi;
    (*defaultParamsArrayPtr)[SIGMA_Y] = sigmaY;
 
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "LAMBDA: " << (*defaultParamsArrayPtr)[LAMBDA];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "MU: " << (*defaultParamsArrayPtr)[MU];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "H: " << (*defaultParamsArrayPtr)[ISOH];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "K: " << (*defaultParamsArrayPtr)[KINK];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "PHI: " << (*defaultParamsArrayPtr)[PHI];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "SIGMA_Y: " << (*defaultParamsArrayPtr)[SIGMA_Y];
 
    paramsArrayPtr = boost::make_shared<ParamsArray>();
    std::copy(defaultParamsArrayPtr->begin(), defaultParamsArrayPtr->end(),
              paramsArrayPtr->begin());
    oldParamsArrayPtr = boost::make_shared<ParamsArray>();
    std::fill(oldParamsArrayPtr->begin(), oldParamsArrayPtr->end(), 0);
 
    MoFEMFunctionReturn(0);
  };
 
  /**
   * @brief Get material parameters from block
   *
   * @param m_field
   * @param pip
   * @param block_name
   * @param sev
   * @return MoFEMErrorCode
   */
  MoFEMErrorCode
  addMatBlockOps(MoFEM::Interface &m_field,
                 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
                 std::string block_name, Sev sev) {
    MoFEMFunctionBegin;
 
    struct OpMatBlocks : public ForcesAndSourcesCore::UserDataOperator {
 
      using OP = ForcesAndSourcesCore::UserDataOperator;
 
      OpMatBlocks(boost::shared_ptr<ParamsArray> params_array_ptr,
                  boost::shared_ptr<ParamsArray> def_params_array_otr,
                  MoFEM::Interface &m_field, Sev sev,
                  std::vector<const CubitMeshSets *> meshset_vec_ptr)
          : OP(NOSPACE, OP::OPSPACE), paramsArrayPtr(params_array_ptr),
            defParamsArray(def_params_array_otr) {
        CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
                          "Can not get data from block");
      }
 
      MoFEMErrorCode doWork(int side, EntityType type,
                            EntitiesFieldData::EntData &data) {
        MoFEMFunctionBegin;
 
        for (auto &b : blockData) {
 
          if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
            std::copy(b.paramsArray.begin(), b.paramsArray.end(),
                      paramsArrayPtr->begin());
            MoFEMFunctionReturnHot(0);
          }
        }
 
        std::copy(defParamsArray->begin(), defParamsArray->end(),
                  paramsArrayPtr->begin());
 
        MoFEMFunctionReturn(0);
      }
 
    private:
      boost::shared_ptr<ParamsArray> paramsArrayPtr;
      boost::shared_ptr<ParamsArray> defParamsArray;
 
      struct BlockData {
        ParamsArray paramsArray;
        Range blockEnts;
      };
 
      std::vector<BlockData> blockData;
 
      MoFEMErrorCode
      extractBlockData(MoFEM::Interface &m_field,
                       std::vector<const CubitMeshSets *> meshset_vec_ptr,
                       Sev sev) {
        MoFEMFunctionBegin;
 
        for (auto m : meshset_vec_ptr) {
          MOFEM_TAG_AND_LOG("ADOLCPlasticityWold", sev, "JP2 parameters") << *m;
          std::vector<double> block_data;
          CHKERR m->getAttributes(block_data);
          if (block_data.size() != 6) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                    "Expected that block has two attribute");
          }
          auto get_block_ents = [&]() {
            Range ents;
            CHKERR
            m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
            return ents;
          };
 
          blockData.push_back({*defParamsArray, get_block_ents()});
          std::copy(block_data.begin(), block_data.end(),
                    blockData.back().paramsArray.begin());
          const auto young_modulus = blockData.back().paramsArray[LAMBDA];
          const auto poisson_ratio = blockData.back().paramsArray[MU];
 
          // It is assumed that user provide young_modulus and poisson_ratio in
          // first two arguments of the block
 
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "Young modulus: " << (blockData.back().paramsArray)[LAMBDA];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "Poisson ratio: " << (blockData.back().paramsArray)[MU];
 
          blockData.back().paramsArray[LAMBDA] =
              LAMBDA(young_modulus, poisson_ratio);
          blockData.back().paramsArray[MU] = MU(young_modulus, poisson_ratio);
 
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "LAMBDA: " << (blockData.back().paramsArray)[LAMBDA];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "MU: " << (blockData.back().paramsArray)[MU];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "H: " << (blockData.back().paramsArray)[ISOH];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "K: " << (blockData.back().paramsArray)[KINK];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "PHI: " << (blockData.back().paramsArray)[PHI];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "SIGMA_Y: " << (blockData.back().paramsArray)[SIGMA_Y];
        }
        MoFEMFunctionReturn(0);
      }
    };
 
    auto cubit_meshsets_vec_ptr =
        m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
            std::regex((boost::format("%s(.*)") % block_name).str()));
 
    pip.push_back(new OpMatBlocks(paramsArrayPtr, defaultParamsArrayPtr,
                                  m_field, sev, cubit_meshsets_vec_ptr));
 
    MoFEMFunctionReturn(0);
  }
 
  /**
   * @brief Set material parameters at integration point
   *
   * @param tag
   * @param recalculate_elastic_tangent
   * @return MoFEMErrorCode
   */
  MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent) {
    MoFEMFunctionBegin;
    switch (tag) {
    case TypesTags::W: {
      set_param_vec(tapesTags[tag], paramsArrayPtr->size(),
                    paramsArrayPtr->data());
      for (auto p = 0; p != LAST_PARAM; p++) {
        if ((*paramsArrayPtr)[p] != (*oldParamsArrayPtr)[p]) {
          recalculate_elastic_tangent = true;
          std::copy(paramsArrayPtr->begin(), paramsArrayPtr->end(),
                    oldParamsArrayPtr->begin());
          break;
        }
      }
 
    } break;
    case TypesTags::Y:
    case TypesTags::H:
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unknown tag for setParams");
    }
    MoFEMFunctionReturn(0);
  }
 
  //! [Free energy J2]
  MoFEMErrorCode freeHemholtzEnergy() {
    MoFEMFunctionBegin;
 
    // ADOL-C register parameters (those can varied for integration points)
    auto p_lambda = mkparam(lambda);  // 0
    auto p_mu = mkparam(mu);          // 1
    auto p_H = mkparam(H);            // 2
    auto p_K = mkparam(K);            // 3
    auto p_phi = mkparam(phi);        // 4
    auto p_sigma_y = mkparam(sigmaY); // 5
 
    // Operator converting Voight and tensor notation
    auto t_voight_strain_op = voight_to_strain_op();
 
    // Variable in voight notation
    auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
    auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
    auto t_voight_internal_tensor =
        getFTensor1FromPtr<6>(&a_internalVariables[1]);
 
    // Convert variables to tensor notation
    tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
    tPlasticStrain(i, j) =
        t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
    tInternalTensor(i, j) =
        t_voight_strain_op(i, j, Z) * t_voight_internal_tensor(Z);
 
    // Evaluate elastic strain
    tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
    tR = tElasticStrain(i, i);
 
    // Strain energy
    a_w = 0.5 * p_lambda * tR * tR +
          p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
    // Isotropic strain energy
    a_w +=
        a_internalVariables[0] * p_sigma_y +
        0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
    // Kinematic strain energy
    a_w += (0.5 * (1 - p_phi) * p_K) *
           (tInternalTensor(i, j) * tInternalTensor(i, j));
 
    MoFEMFunctionReturn(0);
  }
  //! [Free energy J2]
 
  adouble f;
 
  //! [Yield function J2]
  MoFEMErrorCode evalF() {
    MoFEMFunctionBegin;
    auto t_voight_stress_op = voight_to_stress_op();
    auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
    auto t_back_stress = getFTensor1FromPtr<6>(&a_internalFluxes[1]);
    tStress(i, j) = t_voight_stress_op(i, j, Z) * t_vioght_stress(Z);
    tBackStress(i, j) = t_voight_stress_op(i, j, Z) * t_back_stress(Z);
 
    constexpr double third = boost::math::constants::third<double>();
    constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<double>();
    tR = tStress(i, i) / 3.;
    tDeviator(i, j) = tStress(i, j) - tR * t_kd(i, j) - tBackStress(i, j);
 
    // Fix the constant to mach uniaxial test
    constexpr double c = 3. / 2.;
    f = c * tDeviator(i, j) * tDeviator(i, j);
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode yieldFunction() {
    MoFEMFunctionBegin;
    CHKERR evalF();
    a_y = sqrt(f) - a_internalFluxes[0];
    MoFEMFunctionReturn(0);
  }
  //! [Yield function J2]
 
  MoFEMErrorCode flowPotential() {
    MoFEMFunctionBegin;
    CHKERR evalF();
    a_h = sqrt(f) - a_internalFluxes[0];
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode codedHessianW(vector<double *> hessian)
  {
    MoFEMFunctionBegin;
    lambda = (*paramsArrayPtr)[0];
    mu = (*paramsArrayPtr)[1];
    H = (*paramsArrayPtr)[5];
 
    double lambda_plus_2mu = lambda + 2 * mu;
    // first 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon^p^2}
    // lower symmetric part
    // Diagonal elements
    hessian[0][0] = lambda_plus_2mu;
    hessian[1][1] = lambda_plus_2mu;
    hessian[2][2] = lambda_plus_2mu;
    
    // Off-diagonal elements
    hessian[1][0] = lambda;
    hessian[2][0] = lambda;
    hessian[2][1] = lambda;
    
    // Shear components
    for (int i = 3; i <= 5; ++i) {
        hessian[i][i] = mu;
    }
 
    //next 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon \partial \varepsilon^p}
    double neg_lambda_plus_2mu = -(lambda + 2 * mu);
    double neg_lambda = -lambda;
    double neg_mu = -mu;
 
    for (int i = 6; i <= 8; ++i) {
      for (int j = 0; j <= 2; ++j) {
        hessian[i][j] = (i - 6 == j) ? neg_lambda_plus_2mu : neg_lambda;
      }
    }
 
    hessian[9][3] = neg_mu;
    hessian[10][4] = neg_mu;
    hessian[11][5] = neg_mu;
 
    //next 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon \partial \varepsilon}
    hessian[6][6] = lambda_plus_2mu;
    hessian[7][6] = lambda;
    hessian[7][7] = lambda_plus_2mu;
    hessian[8][6] = lambda;
    hessian[8][7] = lambda;
    hessian[8][8] = lambda_plus_2mu;
    hessian[9][9] = mu;
    hessian[10][10] = mu;
    hessian[11][11] = mu;
 
    // internalVariables is \frac{\partial^2 \psi}{\partial \alpha^2}
    hessian[12][12] = H;
 
    MoFEMFunctionReturn(0);
  }
};
 
//! [J2 2D]
template <> struct J2Plasticity<2> : public J2Plasticity<3> {
  using J2Plasticity<3>::J2Plasticity;
 
  MoFEMErrorCode freeHemholtzEnergy() {
    MoFEMFunctionBegin;
 
    auto p_lambda = mkparam(lambda);  // 0
    auto p_mu = mkparam(mu);          // 1
    auto p_H = mkparam(H);            // 2
    auto p_K = mkparam(K);            // 3
    auto p_phi = mkparam(phi);        // 4
    auto p_sigma_y = mkparam(sigmaY); // 5
 
    auto t_voight_strain_op = voight_to_strain_op();
    auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
    auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
    auto t_voight_internal_tensor =
        getFTensor1FromPtr<6>(&a_internalVariables[1]);
 
    tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
    tPlasticStrain(i, j) =
        t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
    tInternalTensor(i, j) =
        t_voight_strain_op(i, j, Z) * t_voight_internal_tensor(Z);
 
    tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
 
    //! [Plane stress]
    // Calculate st  rain at ezz enforcing plane strain case
    tElasticStrain(2, 2) =
        (-tElasticStrain(0, 0) * p_lambda - tElasticStrain(1, 1) * p_lambda) /
        (p_lambda + 2 * p_mu);
    //! [Plane stress]
 
    tR = tElasticStrain(i, i);
 
    a_w = 0.5 * p_lambda * tR * tR +
          p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
 
    // plastic part
    // isotropic
    a_w +=
        a_internalVariables[0] * p_sigma_y +
        0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
    // kinematic
    a_w += (0.5 * (1 - p_phi) * p_K) *
           (tInternalTensor(i, j) * tInternalTensor(i, j));
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode codedHessianW(vector<double *> hessian) {
    MoFEMFunctionBegin;
    // Not implemented yet
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "2D Energy Hessian not yet implemented");
    MoFEMFunctionReturn(0);
  }
};
//! [J2 2D]
 
/** \brief Paraboloidal yield criterion
 *
 *  See paper for reference \cite ullah2017three and \cite stier2015comparing
 *
 *  Linear hardening
 *  \f[
 *  \psi =
 *  \frac{1}{2} \lambda \textrm{tr}[\varepsilon]^2 + \mu \varepsilon :
 * \varepsilon
 *  +
 *  \sigma_{yt}\alpha_0
 *  +
 *  \frac{1}{2} H_t \alpha_0^2
 *  +
 * \sigma_{yc}\alpha_1
 *  +
 *  \frac{1}{2} H_c \alpha_1^2
 *  \f]
 *
 *  Exponential hardening
 *  \f[
 *  \psi =
 *  \frac{1}{2} \lambda \textrm{tr}[\varepsilon]^2 + \mu \varepsilon :
 * \varepsilon\\
 *  +
 *  (\sigma_{yt}+H_t)\alpha_0
 *  +
 *  \frac{H_t}{n_t} \exp{(-n_t \alpha_0)}\\
 *  +
 *  (\sigma_{yc}+H_c)\alpha_0
 *  +
 *  \frac{H_c}{n_c} \exp{(-n_c \alpha_1)}
 *  \f]
 * 
 *  \f[
 *  I_1 = \textrm{tr} (\boldsymbol{\sigma})
 * \f]
 * 
 * \f[
 * \eta=\textrm{dev}[\boldsymbol{\sigma}]
 * \f]
 * 
 * \f[
 * J_2 = \frac{1}{2} \eta:\eta
 * \f] 
 *
 * \f[
 * y = 6J_2 + 2I_1\left(\overline{\alpha_1}-\overline{\alpha_0}\right) -
 * 2\overline{\alpha_0} \,\overline{\alpha_1} \f] 
 * where
 * \f[
 * \overline{\alpha_0}=\frac{\partial \psi}{\partial \alpha_0}=\sigma_{yt} + H_t
 * \alpha_0 
 * \f]
 * 
 * \f[
 * \overline{\alpha_1}=\frac{\partial \psi}{\partial \alpha_1}=\sigma_{yc} + H_c
 * \alpha_1 
 * \f] 
 *
 * \f[
 * \Psi = 6J_2 + 2\alpha I_1 \left(\overline{\alpha_1}-\overline{\alpha_0}\right)
 * - 2\overline{\alpha_0} \,\overline{\alpha_1} \f] \f[ \alpha=
 * \frac{1-2\nu_p}{1+\nu_p} 
 * \f]
 * 
 * \ingroup adoc_plasticity
 */
struct ParaboloidalPlasticity : public ClosestPointProjection {
 
  ParaboloidalPlasticity() : ClosestPointProjection() {
    nbInternalVariables = 2;
    activeVariablesW.resize(12 + nbInternalVariables, false);
  }
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'Z', 6> Z;
 
  FTensor::Tensor2_symmetric<adouble, 3> tTotalStrain;
  FTensor::Tensor2_symmetric<adouble, 3> tPlasticStrain;
  FTensor::Tensor2_symmetric<adouble, 3> tElasticStrain;
  FTensor::Tensor2_symmetric<adouble, 3> tStress;
  FTensor::Tensor2_symmetric<adouble, 3> tDeviator;
 
  adouble tR;
 
  double lambda;           //< Lamé parameters
  double mu;               //< Lamé parameter
  double nup;              //< Plastic Poisson’s ratio
  double Ht, Hc;           //< Isotropic hardening for tension and compression
  double sigmaYt, sigmaYc; //< Yield stress in tension and compression
  double nt, nc;           //< Hardening parameters for tension and compression
 
  enum ModelParams {
    LAMBDA,
    MU,
    NUP,
    SIGMA_YT,
    SIGMA_YC,
    HT,
    HC,
    NT,
    NC,
    LAST_PARAM
  };
 
  using ParamsArray = std::array<double, LAST_PARAM>;
  boost::shared_ptr<ParamsArray> defaultParamsArrayPtr = nullptr;
  boost::shared_ptr<ParamsArray> paramsArrayPtr = nullptr;
  boost::shared_ptr<ParamsArray> oldParamsArrayPtr = nullptr;
 
  MoFEMErrorCode getDefaultMaterialParameters() {
    MoFEMFunctionBegin;
 
    double young_modulus = 1;
    double poisson_ratio = 0.25;
 
    sigmaYt = 1;
    sigmaYc = 1;
    Ht = 0.1;
    Hc = 0.1;
    nc = 1.;
    nt = 1.;
    nup = 0.;
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-young_modulus", &young_modulus,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-poisson_ratio", &poisson_ratio,
                                 0);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-sigma_yt", &sigmaYt, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-sigma_yc", &sigmaYc, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-Ht", &Ht, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-Hc", &Hc, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-nt", &nt, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-nc", &nc, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-nup", &nup, PETSC_NULLPTR);
 
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "Young modulus: " << young_modulus;
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "Poisson ratio: " << poisson_ratio;
 
    mu = MU(young_modulus, poisson_ratio);
    lambda = LAMBDA(young_modulus, poisson_ratio);
 
    defaultParamsArrayPtr = boost::make_shared<ParamsArray>();
    (*defaultParamsArrayPtr)[LAMBDA] = lambda;
    (*defaultParamsArrayPtr)[MU] = mu;
    (*defaultParamsArrayPtr)[NUP] = nup;
    (*defaultParamsArrayPtr)[SIGMA_YT] = sigmaYt;
    (*defaultParamsArrayPtr)[SIGMA_YC] = sigmaYc;
    (*defaultParamsArrayPtr)[HT] = Ht;
    (*defaultParamsArrayPtr)[HC] = Hc;
    (*defaultParamsArrayPtr)[NT] = nt;
    (*defaultParamsArrayPtr)[NC] = nc;
 
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "LAMBDA: " << (*defaultParamsArrayPtr)[LAMBDA];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "MU: " << (*defaultParamsArrayPtr)[MU];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "NUP: " << (*defaultParamsArrayPtr)[NUP];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "SIGMA_YT: " << (*defaultParamsArrayPtr)[SIGMA_YT];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "SIGMA_YC: " << (*defaultParamsArrayPtr)[SIGMA_YC];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "HT: " << (*defaultParamsArrayPtr)[HT];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "HC: " << (*defaultParamsArrayPtr)[HC];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "NT: " << (*defaultParamsArrayPtr)[NT];
    MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
        << "NC: " << (*defaultParamsArrayPtr)[NC];
 
    paramsArrayPtr = boost::make_shared<ParamsArray>();
    std::copy(defaultParamsArrayPtr->begin(), defaultParamsArrayPtr->end(),
              paramsArrayPtr->begin());
    oldParamsArrayPtr = boost::make_shared<ParamsArray>();
    std::fill(oldParamsArrayPtr->begin(), oldParamsArrayPtr->end(), 0);
 
    MoFEMFunctionReturn(0);
  };
 
  MoFEMErrorCode
  addMatBlockOps(MoFEM::Interface &m_field,
                 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
                 std::string block_name, Sev sev) {
    MoFEMFunctionBegin;
 
    struct OpMatBlocks : public ForcesAndSourcesCore::UserDataOperator {
 
      using OP = ForcesAndSourcesCore::UserDataOperator;
 
      OpMatBlocks(boost::shared_ptr<ParamsArray> params_array_ptr,
                  boost::shared_ptr<ParamsArray> def_params_array_ptr,
                  MoFEM::Interface &m_field, Sev sev,
                  std::vector<const CubitMeshSets *> meshset_vec_ptr)
          : OP(NOSPACE, OP::OPSPACE), paramsArrayPtr(params_array_ptr),
            defParamsArray(def_params_array_ptr) {
        CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
                          "Can not get data from block");
      }
 
      MoFEMErrorCode doWork(int side, EntityType type,
                            EntitiesFieldData::EntData &data) {
        MoFEMFunctionBegin;
 
        for (auto &b : blockData) {
 
          if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
            std::copy(b.paramsArray.begin(), b.paramsArray.end(),
                      paramsArrayPtr->begin());
            MoFEMFunctionReturnHot(0);
          }
        }
 
        std::copy(defParamsArray->begin(), defParamsArray->end(),
                  paramsArrayPtr->begin());
 
        MoFEMFunctionReturn(0);
      }
 
    private:
      boost::shared_ptr<ParamsArray> paramsArrayPtr;
      boost::shared_ptr<ParamsArray> defParamsArray;
 
      struct BlockData {
        ParamsArray paramsArray;
        Range blockEnts;
      };
 
      std::vector<BlockData> blockData;
 
      MoFEMErrorCode
      extractBlockData(MoFEM::Interface &m_field,
                       std::vector<const CubitMeshSets *> meshset_vec_ptr,
                       Sev sev) {
        MoFEMFunctionBegin;
 
        for (auto m : meshset_vec_ptr) {
          MOFEM_TAG_AND_LOG("ADOLCPlasticityWold", sev, "Paraboloidal") << *m;
          std::vector<double> block_data;
          CHKERR m->getAttributes(block_data);
          if (block_data.size() != 9) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                    "Expected that block has nine attributes");
          }
          auto get_block_ents = [&]() {
            Range ents;
            CHKERR
            m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
            return ents;
          };
 
          blockData.push_back({*defParamsArray, get_block_ents()});
          std::copy(block_data.begin(), block_data.end(),
                    blockData.back().paramsArray.begin());
 
          const auto young_modulus = blockData.back().paramsArray[LAMBDA];
          const auto poisson_ratio = blockData.back().paramsArray[MU];
 
          // It is assumed that user provide young_modulus and poisson_ratio in
          // first two arguments of the block
 
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "Young modulus: " << (blockData.back().paramsArray)[LAMBDA];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "Poisson ratio: " << (blockData.back().paramsArray)[MU];
 
          blockData.back().paramsArray[LAMBDA] =
              LAMBDA(young_modulus, poisson_ratio);
          blockData.back().paramsArray[MU] = MU(young_modulus, poisson_ratio);
 
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "LAMBDA: " << (blockData.back().paramsArray)[LAMBDA];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "MU: " << (blockData.back().paramsArray)[MU];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "NUP: " << (blockData.back().paramsArray)[NUP];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "SIGMA_YT: " << (blockData.back().paramsArray)[SIGMA_YT];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "SIGMA_YC: " << (blockData.back().paramsArray)[SIGMA_YC];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "HT: " << (blockData.back().paramsArray)[HT];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "HC: " << (blockData.back().paramsArray)[HC];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "NT: " << (blockData.back().paramsArray)[NT];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "NC: " << (blockData.back().paramsArray)[NC];
        }
        MoFEMFunctionReturn(0);
      }
    };
 
    auto cubit_meshsets_vec_ptr =
        m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
            std::regex((boost::format("%s(.*)") % block_name).str()));
 
    pip.push_back(new OpMatBlocks(paramsArrayPtr, defaultParamsArrayPtr,
                                  m_field, sev, cubit_meshsets_vec_ptr));
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode setParams(short tag, bool &recalculate_elastic_tangent) {
        MoFEMFunctionBegin;
    switch (tag) {
    case TypesTags::W: {
      set_param_vec(tapesTags[tag], paramsArrayPtr->size(),
                    paramsArrayPtr->data());
      for (auto p = 0; p != LAST_PARAM; p++) {
        if ((*paramsArrayPtr)[p] != (*oldParamsArrayPtr)[p]) {
          recalculate_elastic_tangent = true;
          std::copy(paramsArrayPtr->begin(), paramsArrayPtr->end(),
                    oldParamsArrayPtr->begin());
          break;
        }
      }
 
    } break;
    case TypesTags::Y:
    case TypesTags::H:
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Unknown tag for setParams");
    }
    MoFEMFunctionReturn(0);
    return 0;
  }
 
  //! [Paraboloidal free energy]
  MoFEMErrorCode freeHemholtzEnergy() {
    MoFEMFunctionBegin;
 
    auto p_lambda = mkparam(lambda);  // 0
    auto p_mu = mkparam(mu);          // 1
    auto p_nup = mkparam(nup);            // 2
    auto p_sigma_yt = mkparam(sigmaYt); // 3
    auto p_sigma_yc = mkparam(sigmaYc); // 4
    auto p_Ht = mkparam(Ht);            // 5
    auto p_Hc = mkparam(Hc);            // 6
    auto p_nt = mkparam(nt);            // 7
    auto p_nc = mkparam(nc);            // 8
 
    auto t_voight_strain_op = voight_to_strain_op();
    auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
    auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
 
    tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
    tPlasticStrain(i, j) =
        t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
    tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
    tR = tElasticStrain(i, i);
 
    a_w = 0.5 * p_lambda * tR * tR +
          p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
 
    a_w += p_sigma_yt * a_internalVariables[0] +
           0.5 * p_Ht * a_internalVariables[0] * a_internalVariables[0];
    a_w += p_sigma_yc * a_internalVariables[1] +
           0.5 * p_Hc * a_internalVariables[1] * a_internalVariables[1];
 
    MoFEMFunctionReturn(0);
  }
  //! [Paraboloidal free energy]
 
 
  MoFEMErrorCode codedHessianW(vector<double *> hessian)
  {
    MoFEMFunctionBegin;
    // first 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon^p^2}
    // lower symmetric part
 
    lambda = (*paramsArrayPtr)[0];
    mu = (*paramsArrayPtr)[1];
    
    Ht = (*paramsArrayPtr)[5];
    Hc = (*paramsArrayPtr)[6];
 
    double lambda_plus_2mu = lambda + 2 * mu;
    
    // Diagonal elements
    hessian[0][0] = lambda_plus_2mu;
    hessian[1][1] = lambda_plus_2mu;
    hessian[2][2] = lambda_plus_2mu;
    
    // Off-diagonal elements
    hessian[1][0] = lambda;
    hessian[2][0] = lambda;
    hessian[2][1] = lambda;
    
    // Shear components
    for (int i = 3; i <= 5; ++i) {
        hessian[i][i] = mu;
    }
 
    //next 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon \partial \varepsilon^p}
    double neg_lambda_plus_2mu = -(lambda + 2 * mu);
    double neg_lambda = -lambda;
    double neg_mu = -mu;
 
    for (int i = 6; i <= 8; ++i) {
      for (int j = 0; j <= 2; ++j) {
        hessian[i][j] = (i - 6 == j) ? neg_lambda_plus_2mu : neg_lambda;
      }
    }
 
    hessian[9][3] = neg_mu;
    hessian[10][4] = neg_mu;
    hessian[11][5] = neg_mu;
 
    //next 6x6 is \frac{\partial^2 \psi}{\partial \varepsilon \partial \varepsilon}
    hessian[6][6] = lambda_plus_2mu;
    hessian[7][6] = lambda;
    hessian[7][7] = lambda_plus_2mu;
    hessian[8][6] = lambda;
    hessian[8][7] = lambda;
    hessian[8][8] = lambda_plus_2mu;
    hessian[9][9] = mu;
    hessian[10][10] = mu;
    hessian[11][11] = mu;
 
    // internalVariables is \frac{\partial^2 \psi}{\partial \alpha^2}
    hessian[12][12] = Ht;
    hessian[13][13] = Hc;
 
    MoFEMFunctionReturn(0);
  }
 
  adouble I1, J2;
 
//! [Paraboloidal yield function]
  MoFEMErrorCode evalF() {
    MoFEMFunctionBegin;
    auto t_voight_stress_op = voight_to_stress_op();
    auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
    tStress(i, j) = t_voight_stress_op(i, j, Z) * t_vioght_stress(Z);
    tR = tStress(i, i);
    I1 = tR;
    constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<double>();
    tR /= 3.;
    tDeviator(i, j) = tStress(i, j) - tR * t_kd(i, j);
    J2 = 0.5 * (tDeviator(i, j) * tDeviator(i, j));
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode yieldFunction() {
    MoFEMFunctionBegin;
    CHKERR evalF();
    a_y = 6 * J2 + 2 * I1 * (a_internalFluxes[1] - a_internalFluxes[0]) -
          2 * a_internalFluxes[0] * a_internalFluxes[1];
    MoFEMFunctionReturn(0);
  }
//! [Paraboloidal yield function]
 
//! [Paraboloidal flow potential]
  MoFEMErrorCode flowPotential() {
    MoFEMFunctionBegin;
    CHKERR evalF();
    double alpha =
        (1 - 2 * nup) /
        (1 + nup); // relation between alpha and plastic poisson's ratio
    a_h = 6 * J2 +
          2 * alpha * I1 * (a_internalFluxes[1] - a_internalFluxes[0]) -
          2 * a_internalFluxes[0] * a_internalFluxes[1];
    MoFEMFunctionReturn(0);
  }
//! [Paraboloidal flow potential]
};
 
struct HeterogeneousParaboloidalPlasticity : public ParaboloidalPlasticity {
 
  using ParaboloidalPlasticity::ParaboloidalPlasticity;
 
  MoFEMErrorCode
  addMatBlockOps(MoFEM::Interface &m_field,
                 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
                 std::string block_name, Sev sev) override {
    MoFEMFunctionBegin;
 
    struct OpMatBlocks : public ForcesAndSourcesCore::UserDataOperator {
 
      using OP = ForcesAndSourcesCore::UserDataOperator;
 
      double poissonRatio = 0.25;
      Tag thYoungModulus;
      Tag thCompressiveYieldStress;
      Tag thTensionYieldStress;
      double youngsModulus;
      double yieldStrengthC;
      double yieldStrengthT;
 
      OpMatBlocks(boost::shared_ptr<ParamsArray> params_array_ptr,
                  boost::shared_ptr<ParamsArray> def_params_array_ptr,
                  MoFEM::Interface &m_field, Sev sev,
                  std::vector<const CubitMeshSets *> meshset_vec_ptr)
          : OP(NOSPACE, OP::OPSPACE), paramsArrayPtr(params_array_ptr),
            defParamsArray(def_params_array_ptr), mField(m_field) {
        CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-poisson_ratio",
                                     &poissonRatio, 0);
        rval = mField.get_moab().tag_get_handle(
            "YOUNGS_MODULUS", 1, MB_TYPE_DOUBLE, thYoungModulus, MB_TAG_SPARSE);
        if (rval != MB_SUCCESS) {
          MOFEM_LOG("ADOLCPlasticityWold", Sev::warning)
              << "YOUNGS_MODULUS tag does not exist using default value";
          CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-young_modulus",
                                       &youngsModulus, 0);
        }
 
        rval = mField.get_moab().tag_get_handle(
            "Yield_Strength_C", 1, MB_TYPE_DOUBLE, thCompressiveYieldStress,
            MB_TAG_SPARSE);
 
        if (rval != MB_SUCCESS) {
          MOFEM_LOG("ADOLCPlasticityWold", Sev::warning)
              << "Yield_Strength_C tag does not exist using default value";
          yieldStrengthC = defParamsArray->at(SIGMA_YC);
        }
 
        rval = mField.get_moab().tag_get_handle(
            "Yield_Strength_T", 1, MB_TYPE_DOUBLE, thTensionYieldStress,
            MB_TAG_SPARSE);
 
        if (rval != MB_SUCCESS) {
          MOFEM_LOG("ADOLCPlasticityWold", Sev::warning)
              << "Yield_Strength_T tag does not exist using default value";
          yieldStrengthT = defParamsArray->at(SIGMA_YT);
        }
 
        CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
                          "Can not get data from block");
      }
 
      MoFEMErrorCode doWork(int side, EntityType type,
                            EntitiesFieldData::EntData &data) {
        MoFEMFunctionBegin;
 
        for (auto &b : blockData) {
 
          if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
            EntityHandle fe_ent = getFEEntityHandle();
 
            if (thYoungModulus != 0)
              CHKERR mField.get_moab().tag_get_data(thYoungModulus, &fe_ent, 1,
                                                    &youngsModulus);
 
            if (thCompressiveYieldStress != 0)
              CHKERR mField.get_moab().tag_get_data(
                  thCompressiveYieldStress, &fe_ent, 1, &yieldStrengthC);
 
            if (thTensionYieldStress != 0)
              CHKERR mField.get_moab().tag_get_data(
                  thTensionYieldStress, &fe_ent, 1, &yieldStrengthT);
 
            std::copy(b.paramsArray.begin(), b.paramsArray.end(),
                      paramsArrayPtr->begin());
 
            (*paramsArrayPtr)[LAMBDA] = LAMBDA(youngsModulus, poissonRatio);
            (*paramsArrayPtr)[MU] = MU(youngsModulus, poissonRatio);
            (*paramsArrayPtr)[SIGMA_YC] = yieldStrengthC;
            (*paramsArrayPtr)[SIGMA_YT] = yieldStrengthT;
 
            MoFEMFunctionReturnHot(0);
          }
        }
 
        std::copy(defParamsArray->begin(), defParamsArray->end(),
                  paramsArrayPtr->begin());
 
        MoFEMFunctionReturn(0);
      }
 
    private:
      boost::shared_ptr<ParamsArray> paramsArrayPtr;
      boost::shared_ptr<ParamsArray> defParamsArray;
      MoFEM::Interface &mField;
 
      struct BlockData {
        ParamsArray paramsArray;
        Range blockEnts;
      };
 
      std::vector<BlockData> blockData;
 
      MoFEMErrorCode
      extractBlockData(MoFEM::Interface &m_field,
                       std::vector<const CubitMeshSets *> meshset_vec_ptr,
                       Sev sev) {
        MoFEMFunctionBegin;
 
        for (auto m : meshset_vec_ptr) {
          MOFEM_TAG_AND_LOG("ADOLCPlasticityWold", sev,
                            "Heterogeneous Paraboloidal")
              << *m;
          std::vector<double> block_data;
          CHKERR m->getAttributes(block_data);
          if (block_data.size() != 0) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                    "Expected that block has zero attributes");
          }
          auto get_block_ents = [&]() {
            Range ents;
            CHKERR
            m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
            return ents;
          };
 
          blockData.push_back({*defParamsArray, get_block_ents()});
          std::copy(block_data.begin(), block_data.end(),
                    blockData.back().paramsArray.begin());
          const auto young_modulus = blockData.back().paramsArray[LAMBDA];
          const auto poisson_ratio = blockData.back().paramsArray[MU];
 
          // It is assumed that user provide young_modulus and poisson_ratio in
          // first two arguments of the block
 
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "Young modulus: " << (blockData.back().paramsArray)[LAMBDA];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "Poisson ratio: " << (blockData.back().paramsArray)[MU];
 
          blockData.back().paramsArray[LAMBDA] =
              LAMBDA(young_modulus, poisson_ratio);
          blockData.back().paramsArray[MU] = MU(young_modulus, poisson_ratio);
 
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "LAMBDA: " << (blockData.back().paramsArray)[LAMBDA];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "MU: " << (blockData.back().paramsArray)[MU];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "NUP: " << (blockData.back().paramsArray)[NUP];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "SIGMA_YT: " << (blockData.back().paramsArray)[SIGMA_YT];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "SIGMA_YC: " << (blockData.back().paramsArray)[SIGMA_YC];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "HT: " << (blockData.back().paramsArray)[HT];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "HC: " << (blockData.back().paramsArray)[HC];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "NT: " << (blockData.back().paramsArray)[NT];
          MOFEM_LOG("ADOLCPlasticityWold", sev)
              << "NC: " << (blockData.back().paramsArray)[NC];
        }
        MoFEMFunctionReturn(0);
      }
    };
 
    auto cubit_meshsets_vec_ptr =
        m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
            std::regex((boost::format("%s(.*)") % block_name).str()));
 
    pip.push_back(new OpMatBlocks(paramsArrayPtr, defaultParamsArrayPtr,
                                  m_field, sev, cubit_meshsets_vec_ptr));
 
    MoFEMFunctionReturn(0);
  }
};
 
template <typename MODEL>
inline auto createMaterial(
    std::array<int, ClosestPointProjection::LAST_TAPE> tapes_tags = {0, 1, 2}) {
  auto cp_ptr = boost::make_shared<MODEL>();
  cp_ptr->getDefaultMaterialParameters();
  cp_ptr->createMatAVecR();
  cp_ptr->snesCreate();
  cp_ptr->recordTapes();
  cp_ptr->tapesTags = tapes_tags;
  return cp_ptr;
}
 
} // namespace ADOLCPlasticity
 
#endif //__ADOLC_MATERIAL_MODELS_HPP