ADOLCPlasticity.hpp

Operators and data structures for small strain plasticity

/** \file ADOLCPlasticity.hpp
 * \ingroup adoc_plasticity
 * \example ADOLCPlasticity.hpp
 *
 * \brief Operators and data structures for small strain plasticity
 *
 * \defgroup adoc_plasticity ADOL-C plasticity
 * \ingroup user_modules
 * \defgroup user_modules User modules
 *
**/
 
#ifndef __ADOLCPLASTICITY_HPP_
#define __ADOLCPLASTICITY_HPP_
 
#ifndef WITH_ADOL_C
#error "MoFEM need to be compiled with ADOL-C"
#endif
 
/**
 *  \ingroup adoc_plasticity
 * 
 */
namespace ADOLCPlasticity {
 
/**
 * \brief Op convert Vight strain vector to strain tensor
*/
inline FTensor::Dg<double, 3, 6> voight_to_strain_op() {
  return FTensor::Dg<double, 3, 6>{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                                   0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0,
                                   0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                                   0.0, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
};
 
/**
 * \brief Op convert strain tensor to Vight strain vector
*/
inline FTensor::Dg<double, 3, 6> strain_to_voight_op() {
  return FTensor::Dg<double, 3, 6>{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                                   1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
                                   0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                                   0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
};
 
/**
 * \brief Op convert Vight stress vector to stress tensor
*/
inline FTensor::Dg<double, 3, 6> voight_to_stress_op() {
  return FTensor::Dg<double, 3, 6>{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                                   1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
                                   0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                                   0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
};
 
/** \brief common data used by volume elements
 * \ingroup nonlinear_elastic_elem
 */
struct CommonData : boost::enable_shared_from_this<CommonData> {
 
  MatrixDouble activeVariablesW;
  MatrixDouble gradientW;
 
  inline auto getFTensor1StressAtGaussPts(int gg = 0) {
    return FTensor::Tensor1<double *, 6>{&gradientW(gg, 6),
                                         &gradientW(gg, 7),
                                         &gradientW(gg, 8),
                                         &gradientW(gg, 9),
                                         &gradientW(gg, 10),
                                         &gradientW(gg, 11),
                                         static_cast<int>(gradientW.size2())};
  };
 
  inline auto getFTensor1PlasticStrainAtGaussPts(int gg = 0) {
    return FTensor::Tensor1<double *, 6>{
        &activeVariablesW(gg, 0),
        &activeVariablesW(gg, 1),
        &activeVariablesW(gg, 2),
        &activeVariablesW(gg, 3),
        &activeVariablesW(gg, 4),
        &activeVariablesW(gg, 5),
        static_cast<int>(activeVariablesW.size2())};
  };
 
  MatrixDouble gradAtGaussPts;
  MatrixDouble materialTangentOperator;
 
  inline auto getPlasticStrain(int gg = 0) {
    return getVectorAdaptor(&(activeVariablesW(gg, 0)), 6);
  }
 
  inline auto getInternalVariables(int gg = 0) {
    return getVectorAdaptor(&(activeVariablesW(gg, 12)),
                            activeVariablesW.size2() - 12);
  }
 
  vector<double> deltaGamma; //< Lagrange plastic multiplier
  double *plasticStrainPtr;
  double *internalVariablesPtr;
  int plasticStrainSize;
  int internalVariablesSize;
 
  inline auto getGardAtGaussPtsPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &gradAtGaussPts);
  }
 
  bool bBar = true;
 
  MoFEMErrorCode getDefaultMaterialParameters() {
    MoFEMFunctionBegin;
    PetscBool b_bar = bBar ? PETSC_TRUE : PETSC_FALSE;
    CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-b_bar", &b_bar,
                               PETSC_NULL);
    bBar = b_bar;
    MoFEMFunctionReturn(0);
  }
 
  boost::shared_ptr<MatrixDouble> getStressMatrixPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &stressMatrix);
  }
  boost::shared_ptr<MatrixDouble> getPlasticStrainMatrixPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &plasticStrainMatrix);
  }
 
  MatrixDouble stressMatrix;
  MatrixDouble plasticStrainMatrix;
 
  CommonData() {
 
    CHK_THROW_MESSAGE(getDefaultMaterialParameters(), "get parameters failed");
  }
};
 
/** \brief Closest point projection algorithm
* \ingroup small_strain_plasticity
 
*/
struct ClosestPointProjection {
 
  int nbInternalVariables;
  boost::function<int(int, int, int)> integrationRule = [](int, int, int p) {
    return 2 * (p - 1);
  };
 
  VectorDouble internalVariables0;
  VectorDouble plasticStrain0;
 
  inline VectorAdaptor getPlasticStrain() {
    return getVectorAdaptor(&(activeVariablesW[0]), 6);
  }
  inline VectorAdaptor getTotalStrain() {
    return getVectorAdaptor(&(activeVariablesW[6]), 6);
  }
  inline VectorAdaptor getInternalVariables() {
    return getVectorAdaptor(&(activeVariablesW[12]), nbInternalVariables);
  }
  inline VectorAdaptor getActiveVariablesYH() {
    return getVectorAdaptor(&(gradientW[6]), 6 + nbInternalVariables);
  }
  inline VectorAdaptor getStress() {
    return getVectorAdaptor(&(gradientW[6]), 6);
  }
  inline VectorAdaptor getInternalFluxes() {
    return getVectorAdaptor(&(gradientW[12]), nbInternalVariables);
  }
 
  enum TypesTags { W = 0, Y, H, LAST_TAPE }; //< W - energy, Y - yield, H - flow
  std::array<int, LAST_TAPE> tapesTags;      //< Tapes nmbers
 
  ClosestPointProjection();
 
  VectorAdaptor activeVariablesW;
  VectorAdaptor gradientW;
 
  double w;
  double y;
  double h;
 
  double deltaGamma; // Increment of plastic multiplier
  MatrixDouble Ep, Cp, Cep;
  ublas::symmetric_matrix<double, ublas::lower> C, D;
 
  /**
   * \brief Record strain energy
   */
  MoFEMErrorCode recordW();
 
  /**
   * \brief Record yield function
   */
  MoFEMErrorCode recordY();
 
  /**
   * \brief Record flow potential
   */
  MoFEMErrorCode recordH();
 
  MatrixDouble hessianW; //< Hessian of energy
  MoFEMErrorCode playW();
  MoFEMErrorCode playW_NoHessian();
 
  VectorDouble gradientY; //< Gradient of yield function
  MoFEMErrorCode playY();
  MoFEMErrorCode playY_NoGradient();
 
  VectorDouble gradientH; //< Gradient of flow potential
  MatrixDouble hessianH;  //< Hessian of flow potential
  MoFEMErrorCode playH();
  MoFEMErrorCode playH_NoHessian();
 
  MoFEMErrorCode createMatAVecR(); //< For integration point SNES solver
 
  MoFEMErrorCode evaluatePotentials(); //< Evaluate potentials
 
  MoFEMErrorCode
  playPotentials(); //< Play potentials from recorded ADOl-C tapes
 
  MoFEMErrorCode
  playPotentials_NoHessian(); //< Play potentials from recorded ADOl-C tapes
 
  MoFEMErrorCode calculateR(Vec R); //< Calculate residual
 
  MoFEMErrorCode calculateA(); //< Calculate tangent matrix
 
  /**
   * \brief Function executed by nested SENES at evaluationg residual
   */
  friend MoFEMErrorCode ADOLCPlasticityRes(SNES, Vec, Vec, void *ctx);
 
  /**
   * \brief Function executed by nested SENES at evaluationg tangent matrix
   */
  friend MoFEMErrorCode ADOLCPlasticityJac(SNES, Vec, Mat, Mat, void *ctx);
 
  /** 
   * \brief Create nested snes
  */
  MoFEMErrorCode snesCreate();
 
  /**
   * \brief Solve nonlinear system of equations to find stress, internal
   * fluxes, and Lagrange plastic multiplier
   */
  MoFEMErrorCode solveClosetProjection();
 
  /**
   * \brief Calculate consistent tangent matrix
   */
  MoFEMErrorCode consistentTangent();
 
  /**
   * \brief Record tapes
   */
  MoFEMErrorCode recordTapes();
 
  /**
   * \brief Get model parameters from blocks
   */
  virtual MoFEMErrorCode
  addMatBlockOps(MoFEM::Interface &m_field,
                 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
                 std::string block_name, Sev sev) {
    return 0;
  }
 
  /**
   * \brief Set parameters for ADOL-C of constitutive relations
   *
   * \param tag [in] - tag of the tape
   * \param recalculate_elastic_tangent [out] - if setParam set it to true,
   * tangent matrix for elastic domain should be recalculated
   */
  virtual MoFEMErrorCode setParams(short tag,
                                   bool &recalculate_elastic_tangent) {
    return 0;
  }
 
  /**
   * \brief Set Hemholtz energy
   */
  virtual MoFEMErrorCode freeHemholtzEnergy() = 0;
 
  /**
   * \brief Set yield function
   */
  virtual MoFEMErrorCode yieldFunction() = 0;
 
  /**
   * \brief Set flow potential
   */
  virtual MoFEMErrorCode flowPotential() = 0;
 
  // protected:
 
  MatrixDouble partialWStrainPlasticStrain; //< Partial derivative of free energy
                                            // with respect to plastic strain
  VectorAdaptor partialYSigma; //< Partial derivative of yield function with
                               // respect to stress
  VectorAdaptor partialYFlux; //< Partial derivative of yield function with
                              // respect to internal fluxes
  VectorAdaptor partialHSigma; //< Partial derivative of flow potential with
                               // respect to stress
  VectorAdaptor partialHFlux; //< Partial derivative of flow potential with
                              // respect to internal fluxes
 
  /// Second partial derivative of flow potential with respect to stresses
  /// and internal
  ublas::symmetric_matrix<double, ublas::lower> partial2HSigma;
  /// Second partial derivative of flow potential with respect to internal
  /// fluxes
  ublas::symmetric_matrix<double, ublas::lower> partial2HFlux;
  /// Mixed decond partial derivative of flow potential with respect to stresses
  /// and internal fluxes
  MatrixDouble partial2HSigmaFlux;
 
  SmartPetscObj<Mat> A; //< Nested SNES tangent matrix 
  SmartPetscObj<Vec> R; //< Nested SNES residual
  SmartPetscObj<Vec> Chi; //< Nested SNES unknown vector
  SmartPetscObj<Vec> dChi; //< Nested SNES increment vector
  ublas::matrix<double, ublas::column_major> dataA;
  SmartPetscObj<SNES> sNes; //< Nested SNES solver
 
  ublas::vector<adouble> a_sTrain;
  ublas::vector<adouble> a_plasticStrain;
  ublas::vector<adouble> a_internalVariables;
  ublas::vector<adouble> a_sTress, a_internalFluxes;
  adouble a_w;
  adouble a_y;
  adouble a_h;
};
 
/**
 * \brief Internal SNES function used at integration points to calulate stress
*/
MoFEMErrorCode ADOLCPlasticityRes(SNES snes, Vec chi, Vec r, void *ctx);
 
/**
 * \brief Internal SNES function used at integration points to calulate
 *  tangent matrix
*/
MoFEMErrorCode ADOLCPlasticityJac(SNES, Vec, Mat, Mat, void *ctx);
 
/**
 * \brief Get opreator to calulate stress
*/
template <int DIM>
ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress(
    MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs);
 
template <>
ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<3>(
    MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs);
 
template <>
ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<2>(
    MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs);
 
/**
 * \brief Assemble right hand side
 *
 * \param field_name [in] - name of field
 * \param common_data_ptr [in] - shared pointer to common data
 *
 * \ingroup small_strain_plasticity
 */
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpRhsImpl;
 
/**
 * \brief Assemble left hand side
 *
 * \param field_name [in] - name of field
 * \param common_data_ptr [in] - shared pointer to common data
 *
 * \ingroup small_strain_plasticity
 */
template <int DIM, IntegrationType I, typename AssemblyDomainEleOp>
struct OpLhsImpl;
 
template <typename DomainEleOp> struct ADOLCPlasticityIntegrators {
 
  template <AssemblyType A> struct Assembly {
 
    using AssemblyDomainEleOp =
        typename FormsIntegrators<DomainEleOp>::template Assembly<A>::OpBase;
 
    template <int DIM, IntegrationType I>
    using OpRhs = OpRhsImpl<DIM, I, AssemblyDomainEleOp>;
 
    template <int DIM, IntegrationType I>
    using OpLhs = OpLhsImpl<DIM, I, AssemblyDomainEleOp>;
  };
};
 
using Pip = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>;
 
/**
 * @brief Assemble the left hand side, i.e. tangent matrix
 * @ingroup adoc_plasticity
 * 
 * @tparam DIM dimension of the problem
 * @tparam A assembly type
 * @tparam I integration type
 * @tparam DomainEleOp operator type
 * @param m_field 
 * @param field_name 
 * @param pip 
 * @param block_name esh block name caring material parameters
 * @param common_data_ptr 
 * @param cp_ptr 
 * @return MoFEMErrorCode 
 */
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode opFactoryDomainRhs(
    MoFEM::Interface &m_field, string field_name, Pip &pip,
    std::string block_name, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<ClosestPointProjection> cp_ptr, Sev sev = Sev::inform) {
  MoFEMFunctionBegin;
  using P = ADOLCPlasticityIntegrators<DomainEleOp>;
  CHKERR cp_ptr->addMatBlockOps(m_field, pip, block_name, sev);
  pip.push_back(
      getRawPtrOpCalculateStress<DIM>(m_field, common_data_ptr, cp_ptr, false));
  pip.push_back(new typename P::template Assembly<A>::template OpRhs<DIM, I>(
      field_name, common_data_ptr));
  MoFEMFunctionReturn(0);
}
/**
 * @brief Assemble the left hand side, i.e. tangent matrix
 * @ingroup adoc_plasticity
 * 
 * @tparam DIM dimension of the problem
 * @tparam A assembly type
 * @tparam I integration type
 * @tparam DomainEleOp operator type
 * @param m_field 
 * @param field_name 
 * @param pip 
 * @param block_name esh block name caring material parameters
 * @param common_data_ptr 
 * @param cp_ptr 
 * @return MoFEMErrorCode 
 */
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode
opFactoryDomainLhs(MoFEM::Interface &m_field, string field_name, Pip &pip,
                   std::string block_name,
                   boost::shared_ptr<CommonData> common_data_ptr,
                   boost::shared_ptr<ClosestPointProjection> cp_ptr) {
  MoFEMFunctionBegin;
  using P = ADOLCPlasticityIntegrators<DomainEleOp>;
  CHKERR cp_ptr->addMatBlockOps(m_field, pip, block_name, Sev::noisy);
  pip.push_back(
      getRawPtrOpCalculateStress<DIM>(m_field, common_data_ptr, cp_ptr, true));
  pip.push_back(new typename P::template Assembly<A>::template OpLhs<DIM, I>(
      field_name, common_data_ptr));
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Push operators to update history variables
 * @ingroup adoc_plasticity
 * 
 * @tparam DIM  dimension of the problem  
 * @param m_field core interface
 * @param pip 
 * @param block_name mesh block name caring material parameters
 * @param common_data_ptr 
 * @param cp_ptr 
 * @return MoFEMErrorCode 
 */
template <int DIM>
MoFEMErrorCode
opFactoryDomainUpdate(MoFEM::Interface &m_field, Pip &pip,
                      std::string block_name,
                      boost::shared_ptr<CommonData> common_data_ptr,
                      boost::shared_ptr<ClosestPointProjection> cp_ptr);
 
/**
 * @copydoc ADOLCPlasticity::opFactoryDomainUpdate
 */
template <>
MoFEMErrorCode
opFactoryDomainUpdate<3>(MoFEM::Interface &m_field, Pip &pip,
                         std::string block_name,
                         boost::shared_ptr<CommonData> common_data_ptr,
                         boost::shared_ptr<ClosestPointProjection> cp_ptr);
 
/**
 * @copydoc ADOLCPlasticity::opFactoryDomainUpdate
 */
template <>
MoFEMErrorCode
opFactoryDomainUpdate<2>(MoFEM::Interface &m_field, Pip &pip,
                         std::string block_name,
                         boost::shared_ptr<CommonData> common_data_ptr,
                         boost::shared_ptr<ClosestPointProjection> cp_ptr);
 
 
 
/**
 * \brief Update internal fluxes (update history variables)
 * \ingroup adoc_plasticity
 */
struct TSUpdate {
  virtual MoFEMErrorCode postProcess(TS ts) = 0;
};
boost::shared_ptr<TSUpdate> createTSUpdate(std::string fe_name,
                                           boost::shared_ptr<FEMethod> fe_ptr);
 
/**
 * \brief Calculate tensorial base functions. Apply bBar method when needed
 */
struct MakeB {
 
  static FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3>
  getFTensor2SymmetricDiffBase(DataForcesAndSourcesCore::EntData &data,
                               MatrixDouble &storage, const bool b_bar,
                               const int nb_integration_pts, double *w_ptr,
                               FTensor::Number<2>);
 
  static FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3>
  getFTensor2SymmetricDiffBase(DataForcesAndSourcesCore::EntData &data,
                               MatrixDouble &storage, const bool b_bar,
                               const int nb_integration_pts, double *w_ptr,
                               FTensor::Number<3>);
};
                             
 
template <int DIM, typename AssemblyDomainEleOp>
struct OpRhsImpl<DIM, GAUSS, AssemblyDomainEleOp> : public AssemblyDomainEleOp {
  OpRhsImpl(std::string field_name,
            boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  MatrixDouble baseStorage; ///< Store tensorial base functions
};
 
template <int DIM, typename AssemblyDomainEleOp>
struct OpLhsImpl<DIM, GAUSS, AssemblyDomainEleOp> : public AssemblyDomainEleOp {
  OpLhsImpl(std::string field_name,
            boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data,
                           DataForcesAndSourcesCore::EntData &col_data);
 
protected:
  boost::shared_ptr<CommonData> commonDataPtr;
  MatrixDouble baseRowStorage; ///< Store tensorial base functions
  MatrixDouble baseColStorage; ///< Store tensorial base functions
};
 
template <int DIM, typename AssemblyDomainEleOp>
OpRhsImpl<DIM, GAUSS, AssemblyDomainEleOp>::OpRhsImpl(
    string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
      commonDataPtr(common_data_ptr) {}
 
//! [OpRhsImpl integrate]
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode OpRhsImpl<DIM, GAUSS, AssemblyDomainEleOp>::iNtegrate(
    EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
  using OP = AssemblyDomainEleOp;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  double *w_ptr = &(OP::getGaussPts()(DIM, 0));
  // Calulate tensorial base functions. Apply bBar method when needed
  auto t_diff = MakeB::getFTensor2SymmetricDiffBase(
      data, baseStorage, commonDataPtr->bBar, OP::getGaussPts().size2(), w_ptr,
      FTensor::Number<DIM>());
  auto t_stress = getFTensor2SymmetricFromMat<3>(commonDataPtr->stressMatrix);
  const auto vol = OP::getMeasure();
  auto t_w = OP::getFTensor0IntegrationWeight();
  for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
    double alpha = vol * t_w;
    for (int bb = 0; bb != OP::nbRows; ++bb) {
      OP::locF[bb] += alpha * (t_stress(i, j) * t_diff(i, j));
      ++t_diff;
    }
    ++t_w;
    ++t_stress;
  }
  MoFEMFunctionReturn(0);
}
//! [OpRhsImpl integrate]
 
template <int DIM, typename AssemblyDomainEleOp>
OpLhsImpl<DIM, GAUSS, AssemblyDomainEleOp>::OpLhsImpl(
    string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : AssemblyDomainEleOp(field_name, field_name,
                          AssemblyDomainEleOp::OPROWCOL),
      commonDataPtr(common_data_ptr) {
  this->sYmm = false; // It has to be false for not-associtive plasticity
}
 
template <int DIM, typename AssemblyDomainEleOp>
MoFEMErrorCode OpLhsImpl<DIM, GAUSS, AssemblyDomainEleOp>::iNtegrate(
    DataForcesAndSourcesCore::EntData &row_data,
    DataForcesAndSourcesCore::EntData &col_data) {
  MoFEMFunctionBegin;
  using OP = AssemblyDomainEleOp;
 
  double *w_ptr = &(OP::getGaussPts()(DIM, 0));
  auto t_diff_row = MakeB::getFTensor2SymmetricDiffBase(
      row_data, baseRowStorage, commonDataPtr->bBar, OP::getGaussPts().size2(),
      w_ptr, FTensor::Number<DIM>());
  auto t_diff_col = MakeB::getFTensor2SymmetricDiffBase(
      col_data, baseColStorage, commonDataPtr->bBar, OP::getGaussPts().size2(),
      w_ptr, FTensor::Number<DIM>());
 
  auto get_ftensor2_symmetric = [&](auto &storage, const auto gg,
                                    const auto rr) {
    return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3>{
        &storage(gg, 6 * rr + 0), &storage(gg, 6 * rr + 1),
        &storage(gg, 6 * rr + 2), &storage(gg, 6 * rr + 3),
        &storage(gg, 6 * rr + 4), &storage(gg, 6 * rr + 5)};
  };
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  auto t_Cep =
      getFTensor4DdgFromMat<3, 3, 1>(commonDataPtr->materialTangentOperator);
  const auto vol = OP::getMeasure();
  auto t_w = OP::getFTensor0IntegrationWeight();
  for (int gg = 0; gg != OP::nbIntegrationPts; ++gg) {
    const double alpha = vol * t_w;
    ++t_w;
 
    for (auto rr = 0; rr != OP::nbRows; ++rr) {
      FTensor::Tensor2_symmetric<double, 3> t_rowCep;
      t_rowCep(k, l) = t_Cep(i, j, k, l) * t_diff_row(i, j);
      auto t_diff_col = get_ftensor2_symmetric(baseColStorage, gg, 0);
      for (auto cc = 0; cc != OP::nbCols; ++cc) {
        OP::locMat(rr, cc) += alpha * (t_rowCep(k, l) * t_diff_col(k, l));
        ++t_diff_col;
      }
      ++t_diff_row;
    }
 
    ++t_Cep;
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace ADOLCPlasticity
 
#endif //__ADOLCPLASTICITY_HPP_