ADOLCPlasticity.cpp

Go to the documentation of this file.

/** 
 * \file ADOLCPlasticity.cpp
 * \ingroup adoc_plasticity
 * \brief Operators and data structures for ADOL-C plasticity
 *
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <ADOLCPlasticity.hpp>
 
namespace ADOLCPlasticity {
 
static FTensor::Dg<double, 3, 6> tStrainToVoightOp = strain_to_voight_op();
 
template <typename T1, typename T2, int DIM1, int DIM2>
double calcAveStrain(bool b_bar, const int nb_gauss_pts,
                            FTensor::Tensor2<T1, DIM1, DIM2> &&t_grad,
                            FTensor::Tensor0<T2> &&t_w) {
  if constexpr (DIM1 != DIM2)
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Case of mixed dimension by gradient not implemented");
 
 
  FTensor::Index<'i', DIM1> i;
  double ave_tr_strain = 0;
  if (b_bar) {
    double v = 0;
    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
      v += t_w;
      ave_tr_strain += t_w * t_grad(i, i);
      ++t_grad;
      ++t_w;
    }
    ave_tr_strain /= (DIM1 * v);
  }
  return ave_tr_strain;
}
 
template <typename T1, typename T2, int DIM1, int DIM2, int DIM3>
FTensor::Tensor1<T2, DIM3>
calcStrain(bool b_bar, double ave_tr_strain,
           FTensor::Tensor2<T1, DIM1, DIM2> &t_grad,
           FTensor::Tensor1<T2, DIM3> &&t_voight_strain) {
 
  if constexpr (DIM1 != DIM2)
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Case of mixed dimension by gradient not implemented");
 
  FTensor::Index<'i', DIM1> i;
  FTensor::Index<'j', DIM2> j;
  FTensor::Index<'Z', DIM3> Z;
  FTensor::Tensor2_symmetric<double, 3> t_strain{0., 0., 0., 0., 0., 0.};
  t_strain(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
  double tr_strain = t_grad(i, i) / DIM1;
  if (b_bar) {
    FTensor::Index<'i', DIM1> i;
    FTensor::Index<'j', DIM2> j;
    auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
    t_strain(i, j) += (ave_tr_strain - tr_strain) * t_kd(i, j);
  } 
  {
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    t_voight_strain(Z) = tStrainToVoightOp(i, j, Z) * t_strain(i, j);
  }
  return t_voight_strain;
};
 
template <typename T1, typename T2, int DIM1, int DIM3>
FTensor::Tensor1<T2, DIM3>
calcStrain(bool b_bar, double ave_tr_strain,
           FTensor::Tensor2_symmetric<T1, DIM1> &t_strain,
           FTensor::Tensor1<T2, DIM3> &&t_voight_strain) {
 
// TODO: implement this function to calculate b_bar strain
 
  FTensor::Index<'Z', DIM3> Z;
  {
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    t_voight_strain(Z) = tStrainToVoightOp(i, j, Z) * t_strain(i, j);
  }
  return t_voight_strain;
};
 
//! [BBar method]
/**
 * @brief Calculate base for B-bar method
 * 
 * @tparam DIM Dimension of problem
 * @param data access to base functions
 * @param storage storage for base functions at integration points
 * @param b_bar flag to sith on B-bar method 
 * @param nb_integration_pts number of integration points 
 * @param w_ptr integration weights
 * @return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3> 
 */
template <int DIM>
FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3>
getFTensor2SymmetricDiffBaseImpl(DataForcesAndSourcesCore::EntData &data,
                                 MatrixDouble &storage, const bool b_bar,
                                 const int nb_integration_pts, double *w_ptr,
                                 FTensor::Number<DIM>) {
  const auto nb_dofs = data.getFieldData().size();
  const auto nb_bases = data.getN().size2();
 
  storage.resize(nb_integration_pts, 6 * nb_dofs, false);
 
  FTensor::Index<'I', DIM> I;
  FTensor::Index<'J', DIM> J;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Number<0> N0;
  FTensor::Number<1> N1;
  FTensor::Number<2> N2;
 
  // get tensorial base
  auto get_ftensor2_symmetric = [&](const int gg, const int rr) {
    return getFTensor2SymmetricFromPtr<3>(&storage(gg, 6 * rr));
  };
 
  // calculate base 
  auto calc_base = [&]() {
    auto t_t2_diff = get_ftensor2_symmetric(0, 0);
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      auto t_diff = data.getFTensor1DiffN<DIM>(gg, 0);
      for (auto b = 0; b != nb_dofs / DIM; ++b) {
        FTensor::Tensor2<double, 3, 3> t_grad;
 
        t_grad(i, j) = 0;
 
        t_grad(N0, J) = t_diff(J);
        t_t2_diff(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
        ++t_t2_diff;
 
        t_grad(i, j) = 0;
        t_grad(N1, J) = t_diff(J);
        t_t2_diff(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
        ++t_t2_diff;
 
        t_grad(i, j) = 0;
        if constexpr (DIM == 3) {
          t_grad(N2, J) = t_diff(J);
          t_t2_diff(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
          ++t_t2_diff;
        }
 
        ++t_diff;
      }
    }
    return get_ftensor2_symmetric(0, 0);
  };
 
  // calculate volume average of trace
  auto calc_vol = [&](auto &&t_t2_dff) {
    std::vector<double> vol(nb_dofs, 0);
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      for (auto b = 0; b != nb_dofs; ++b) {
        vol[b] += w_ptr[gg] * t_t2_dff(i, i) / DIM;
        ++t_t2_dff;
      }
    }
    double sum = 0;
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      sum += w_ptr[gg];
    }
    for (auto &v : vol) {
      v /= sum;
    }
    return vol;
  };
 
  // modify base for B-bar method
  auto make_b_bar = [&](auto &&vol) {
    auto t_t2_diff = get_ftensor2_symmetric(0, 0);
    constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      for (auto b = 0; b != nb_dofs; ++b) {
        const auto trace = t_t2_diff(J, J) / DIM;
        t_t2_diff(I, J) += (vol[b] - trace) * t_kd(I, J);
        ++t_t2_diff;
      }
    }
    return get_ftensor2_symmetric(0, 0);
  };
 
  return b_bar ? make_b_bar(calc_vol(calc_base())) : calc_base();
};
//! [BBar method]
 
FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3>
MakeB::getFTensor2SymmetricDiffBase(DataForcesAndSourcesCore::EntData &data,
                                    MatrixDouble &storage, const bool b_bar,
                                    const int nb_integration_pts, double *w_ptr,
                                    FTensor::Number<3>) {
  return getFTensor2SymmetricDiffBaseImpl(
      data, storage, b_bar, nb_integration_pts, w_ptr, FTensor::Number<3>());
}
 
FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3>
MakeB::getFTensor2SymmetricDiffBase(DataForcesAndSourcesCore::EntData &data,
                                    MatrixDouble &storage, const bool b_bar,
                                    const int nb_integration_pts, double *w_ptr,
                                    FTensor::Number<2>) {
  return getFTensor2SymmetricDiffBaseImpl(
      data, storage, b_bar, nb_integration_pts, w_ptr, FTensor::Number<2>());
}
 
struct OpUpdate : public ForcesAndSourcesCore::UserDataOperator {
 
  OpUpdate(boost::shared_ptr<CommonData> common_data_ptr);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
 
protected:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
OpUpdate::OpUpdate(boost::shared_ptr<CommonData> common_data_ptr)
    : ForcesAndSourcesCore::UserDataOperator(NOSPACE,
                                             UserDataOperator::OPSPACE),
      commonDataPtr(common_data_ptr) {}
 
MoFEMErrorCode OpUpdate::doWork(int side, EntityType type,
                                DataForcesAndSourcesCore::EntData &data) {
  MoFEMFunctionBegin;
 
  int nb_gauss_pts = getGaussPts().size2();
 
  for (int gg = 0; gg < nb_gauss_pts; gg++) {
    if (commonDataPtr->deltaGamma[gg] > 0) {
      auto cp_plastic_strain = commonDataPtr->getPlasticStrain(gg);
      auto cp_internal_variables = commonDataPtr->getInternalVariables(gg);
      const auto nb_internal_variables = cp_internal_variables.size();
      auto plastic_strain =
          getVectorAdaptor(&commonDataPtr->plasticStrainPtr[gg * 6], 6);
      auto internal_variables = getVectorAdaptor(
          &commonDataPtr->internalVariablesPtr[gg * nb_internal_variables],
          nb_internal_variables);
      plastic_strain = cp_plastic_strain;
      internal_variables = cp_internal_variables;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, StrainType STRAIN>
struct OpCalculateStressImpl : public ForcesAndSourcesCore::UserDataOperator {
 
  OpCalculateStressImpl(MoFEM::Interface &m_field,
                        boost::shared_ptr<CommonData> common_data_ptr,
                        boost::shared_ptr<ClosestPointProjection> cp_ptr,
                        bool calc_lhs);
 
protected:
  MoFEM::Interface &mField;
  boost::shared_ptr<CommonData> commonDataPtr;
  boost::shared_ptr<ClosestPointProjection> cpPtr;
  bool calcLhs;
 
  Tag thPlasticStrain;
  Tag thInternalVariables;
 
  MoFEMErrorCode getTags();
 
  MoFEMErrorCode setTagsData(const EntityHandle tet, const int nb_gauss_pts,
                             const int nb_internal_variables);
};
 
template <int DIM, StrainType STRAIN>
struct OpCalculateStress : public OpCalculateStressImpl<DIM, STRAIN> {
 
  using OpCalculateStressImpl<DIM, STRAIN>::OpCalculateStressImpl;
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data) override;
};
template <int DIM>
struct OpCalculateStress<DIM, LARGE_STRAIN> : public OpCalculateStressImpl<DIM, LARGE_STRAIN> {
 
  using OpCalculateStressImpl<DIM, LARGE_STRAIN>::OpCalculateStressImpl;
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data) override;
};
 
template <int DIM>
struct OpCalculateStress<DIM, SMALL_STRAIN> : public OpCalculateStressImpl<DIM, SMALL_STRAIN> {
 
  using OpCalculateStressImpl<DIM, SMALL_STRAIN>::OpCalculateStressImpl;
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data) override;
};
 
template <int DIM, StrainType STRAIN>
OpCalculateStressImpl<DIM, STRAIN>::OpCalculateStressImpl(
    MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs)
    : ForcesAndSourcesCore::UserDataOperator::UserDataOperator(
          NOSPACE, UserDataOperator::OPSPACE),
      mField(m_field), commonDataPtr(common_data_ptr), cpPtr(cp_ptr),
      calcLhs(calc_lhs) {
  CHK_THROW_MESSAGE(getTags(), "get tags");
}
 
template <int DIM, StrainType STRAIN>
MoFEMErrorCode OpCalculateStressImpl<DIM, STRAIN>::getTags() {
  MoFEMFunctionBegin;
  int def_length = 0;
  CHKERR mField.get_moab().tag_get_handle(
      "PLASTIC_STRAIN", def_length, MB_TYPE_DOUBLE, thPlasticStrain,
      MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN, NULL);
  CHKERR mField.get_moab().tag_get_handle(
      "INTERNAL_VARIABLES", def_length, MB_TYPE_DOUBLE, thInternalVariables,
      MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN, NULL);
  MoFEMFunctionReturn(0);
}
 
template <int DIM, StrainType STRAIN>
MoFEMErrorCode OpCalculateStressImpl<DIM, STRAIN>::setTagsData(
    const EntityHandle tet, const int nb_gauss_pts,
    const int nb_internal_variables) {
  MoFEMFunctionBegin;
 
  VectorDouble v;
  {
    rval = mField.get_moab().tag_get_by_ptr(
        thPlasticStrain, &tet, 1,
        (const void **)&commonDataPtr->plasticStrainPtr,
        &commonDataPtr->plasticStrainSize);
    if (rval != MB_SUCCESS ||
        commonDataPtr->plasticStrainSize != 6 * nb_gauss_pts) {
      v.resize(6 * nb_gauss_pts, false);
      v.clear();
      int tag_size[1];
      tag_size[0] = v.size();
      void const *tag_data[] = {&v[0]};
      CHKERR mField.get_moab().tag_set_by_ptr(thPlasticStrain, &tet, 1,
                                              tag_data, tag_size);
      CHKERR mField.get_moab().tag_get_by_ptr(
          thPlasticStrain, &tet, 1,
          (const void **)&commonDataPtr->plasticStrainPtr,
          &commonDataPtr->plasticStrainSize);
    }
  }
  if (nb_internal_variables > 0) {
    rval = mField.get_moab().tag_get_by_ptr(
        thInternalVariables, &tet, 1,
        (const void **)&commonDataPtr->internalVariablesPtr,
        &commonDataPtr->internalVariablesSize);
    if (rval != MB_SUCCESS || commonDataPtr->internalVariablesSize !=
                                  nb_internal_variables * nb_gauss_pts) {
      v.resize(nb_internal_variables * nb_gauss_pts, false);
      v.clear();
      int tag_size[1];
      tag_size[0] = v.size();
      void const *tag_data[] = {&v[0]};
      CHKERR mField.get_moab().tag_set_by_ptr(thInternalVariables, &tet, 1,
                                              tag_data, tag_size);
      CHKERR mField.get_moab().tag_get_by_ptr(
          thInternalVariables, &tet, 1,
          (const void **)&commonDataPtr->internalVariablesPtr,
          &commonDataPtr->internalVariablesSize);
    }
  }
  MoFEMFunctionReturn(0);
}
 
template <int DIM>
MoFEMErrorCode
OpCalculateStress<DIM, LARGE_STRAIN>::doWork(int side, EntityType type,
                                  DataForcesAndSourcesCore::EntData &data) {
  MoFEMFunctionBegin;
 
  auto do_work_impl = [this](auto commonDataPtr, auto cpPtr) {
    MoFEMFunctionBegin;
 
    int nb_gauss_pts = this->getGaussPts().size2();
 
    int nb_internal_variables = cpPtr->nbInternalVariables;
    auto tet = this->getNumeredEntFiniteElementPtr()->getEnt();
    CHKERR this->setTagsData(tet, nb_gauss_pts, nb_internal_variables);
 
    commonDataPtr->activeVariablesW.resize(
        nb_gauss_pts, 12 + cpPtr->nbInternalVariables, false);
    commonDataPtr->gradientW.resize(nb_gauss_pts,
                                    12 + cpPtr->nbInternalVariables, false);
    commonDataPtr->materialTangentOperator.resize(36, nb_gauss_pts, false);
    commonDataPtr->deltaGamma.resize(nb_gauss_pts);
 
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    FTensor::Index<'k', 3> k;
    FTensor::Index<'l', 3> l;
    FTensor::Index<'Z', 6> Z;
    FTensor::Index<'Y', 6> Y;
 
    // Code uses trial stress, until first solve of linear system of equations
    int iter = 1;
    if (this->getFEMethod()->snes_ctx != SnesMethod::CTX_SNESNONE) {
      CHKERR SNESGetLinearSolveIterations(this->getFEMethod()->snes, &iter);
    }
 
    double ave_tr_strain = 0;
 
    auto t_voight_stress_op = voight_to_stress_op();
 
    auto t_grad =
        getFTensor2SymmetricFromMat<DIM>(commonDataPtr->gradAtGaussPts);
 
    auto t_Cep =
        getFTensor4DdgFromMat<3, 3, 1>(commonDataPtr->materialTangentOperator);
    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
 
      auto tmp_active =
          getVectorAdaptor(&(commonDataPtr->activeVariablesW(gg, 0)),
                           12 + cpPtr->nbInternalVariables);
      cpPtr->activeVariablesW.swap(tmp_active);
      auto tmp_gradient = getVectorAdaptor(&(commonDataPtr->gradientW(gg, 0)),
                                           12 + cpPtr->nbInternalVariables);
      this->cpPtr->gradientW.swap(tmp_gradient);
 
      auto t_voigt_strain =
          calcStrain(commonDataPtr->bBar, ave_tr_strain, t_grad,
                     getFTensor1FromPtr<6>(&(cpPtr->activeVariablesW[6])));
      ++t_grad;
 
      // Copy plastic strains and internal variables from tags
      auto copy_plastic_strain_and_internal_variables = [&]() {
        MoFEMFunctionBegin;
 
        // Get data stored on Tags
        VectorAdaptor plastic_strain_tags =
            VectorAdaptor(6, ublas::shallow_array_adaptor<double>(
                                 6, &commonDataPtr->plasticStrainPtr[gg * 6]));
        VectorAdaptor internal_variables_tags = VectorAdaptor(
            nb_internal_variables,
            ublas::shallow_array_adaptor<double>(
                nb_internal_variables,
                &commonDataPtr
                     ->internalVariablesPtr[gg * nb_internal_variables]));
 
        // Set values to closest point projection
        cpPtr->plasticStrain0.resize(6, false);
        noalias(cpPtr->plasticStrain0) = plastic_strain_tags;
        cpPtr->internalVariables0.resize(nb_internal_variables, false);
        noalias(cpPtr->internalVariables0) = internal_variables_tags;
 
        auto cp_plastic_strain = cpPtr->getPlasticStrain();
        auto cp_internal_variables = cpPtr->getInternalVariables();
        noalias(cp_plastic_strain) = plastic_strain_tags;
        noalias(cp_internal_variables) = internal_variables_tags;
        MoFEMFunctionReturn(0);
      };
 
      CHKERR copy_plastic_strain_and_internal_variables();
 
      auto closest_point_projection = [&](auto &recalculate_elastic_tangent) {
        MoFEMFunctionBegin;
        for (auto t : {ClosestPointProjection::H, ClosestPointProjection::Y,
                       ClosestPointProjection::W}) {
          CHKERR cpPtr->setParams(t, recalculate_elastic_tangent);
        }
        CHKERR cpPtr->playW_NoHessian();
        CHKERR cpPtr->playY_NoGradient();
        double y = cpPtr->y;
        cpPtr->deltaGamma = 0; // zero lagrange multiplier
        if (iter > 0 && y > std::numeric_limits<double>::epsilon()) {
          CHKERR cpPtr->solveClosestProjection();
        }
        commonDataPtr->deltaGamma[gg] = cpPtr->deltaGamma;
        MoFEMFunctionReturn(0);
      };
 
      auto evaluate_consistent_tangent_matrix =
          [&](auto &recalculate_elastic_tangent) {
            MoFEMFunctionBegin;
            if (recalculate_elastic_tangent)
              CHKERR cpPtr->playW();
 
            if (iter > 0 &&
                cpPtr->deltaGamma > std::numeric_limits<double>::epsilon()) {
 
              CHKERR cpPtr->consistentTangent();
 
              auto &m = cpPtr->Cep;
              // for generic case of non-associated plasticity consistent
              // tangent matrix is non-symmetric
              auto t_voight_cep =
                  FTensor::Tensor2<FTensor::PackPtr<double *, 0>, 6, 6>(
                      &m(0, 0), &m(0, 1), &m(0, 2), &m(0, 3), &m(0, 4),
                      &m(0, 5),
 
                      &m(1, 0), &m(1, 1), &m(1, 2), &m(1, 3), &m(1, 4),
                      &m(1, 5),
 
                      &m(2, 0), &m(2, 1), &m(2, 2), &m(2, 3), &m(2, 4),
                      &m(2, 5),
 
                      &m(3, 0), &m(3, 1), &m(3, 2), &m(3, 3), &m(3, 4),
                      &m(3, 5),
 
                      &m(4, 0), &m(4, 1), &m(4, 2), &m(4, 3), &m(4, 4),
                      &m(4, 5),
 
                      &m(5, 0), &m(5, 1), &m(5, 2), &m(5, 3), &m(5, 4), &m(5, 5)
 
                  );
 
              t_Cep(i, j, k, l) =
                  t_voight_stress_op(i, j, Z) *
                  (t_voight_cep(Z, Y) * t_voight_stress_op(k, l, Y));
 
            } else {
 
              auto &m = cpPtr->C;
              auto t_voight_cep =
                  FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 0>, 6>(
                      &m(0, 0), &m(0, 1), &m(0, 2), &m(0, 3), &m(0, 4),
                      &m(0, 5),
 
                      &m(1, 1), &m(1, 2), &m(1, 3), &m(1, 4), &m(1, 5),
 
                      &m(2, 2), &m(2, 3), &m(2, 4), &m(2, 5),
 
                      &m(3, 3), &m(3, 4), &m(3, 5),
 
                      &m(4, 4), &m(4, 5),
 
                      &m(5, 5));
 
              t_Cep(i, j, k, l) =
                  t_voight_stress_op(i, j, Z) *
                  (t_voight_cep(Z, Y) * t_voight_stress_op(k, l, Y));
            }
 
            ++t_Cep;
            MoFEMFunctionReturn(0);
          };
 
      // Always recalculate elastic tangent if first element
      bool recalculate_elastic_tangent =
          (this->getNinTheLoop() == 0) ? true : false;
      CHKERR closest_point_projection(recalculate_elastic_tangent);
      // always calculate consistent tangent matrix for large strain
      CHKERR evaluate_consistent_tangent_matrix(recalculate_elastic_tangent);
    }
 
    auto calcs_stress_matrix = [&]() {
      MoFEMFunctionBegin;
      auto t_voight_stress_op = voight_to_stress_op();
      commonDataPtr->stressMatrix.resize(
          6, this->commonDataPtr->gradientW.size1(), false);
      auto t_stess =
          getFTensor2SymmetricFromMat<3>(commonDataPtr->stressMatrix);
      auto t_voight_stress = commonDataPtr->getFTensor1StressAtGaussPts();
      for (auto gg = 0; gg != commonDataPtr->gradientW.size1(); ++gg) {
        t_stess(i, j) = t_voight_stress_op(i, j, Z) * t_voight_stress(Z);
        ++t_voight_stress;
        ++t_stess;
      }
      MoFEMFunctionReturn(0);
    };
 
    auto calcs_plastic_strain_matrix = [&]() {
      MoFEMFunctionBegin;
      auto t_voight_strain_op = voight_to_stress_op();
      commonDataPtr->plasticStrainMatrix.resize(
          6, commonDataPtr->activeVariablesW.size1(), false);
      auto t_plastic_strain =
          getFTensor2SymmetricFromMat<3>(commonDataPtr->plasticStrainMatrix);
      auto t_voight_plastic_strain =
          commonDataPtr->getFTensor1PlasticStrainAtGaussPts();
      for (auto gg = 0; gg != commonDataPtr->activeVariablesW.size1(); ++gg) {
        t_plastic_strain(i, j) =
            t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
        ++t_voight_plastic_strain;
        ++t_plastic_strain;
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR calcs_stress_matrix();
    CHKERR calcs_plastic_strain_matrix();
    MoFEMFunctionReturn(0);
  };
 
  CHKERR do_work_impl(this->commonDataPtr, this->cpPtr);
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM>
MoFEMErrorCode
OpCalculateStress<DIM, SMALL_STRAIN>::doWork(int side, EntityType type,
                                  DataForcesAndSourcesCore::EntData &data) {
  MoFEMFunctionBegin;
 
  auto do_work_impl = [this](auto commonDataPtr, auto cpPtr) {
    MoFEMFunctionBegin;
 
    int nb_gauss_pts = this->getGaussPts().size2();
 
    int nb_internal_variables = cpPtr->nbInternalVariables;
    auto tet = this->getNumeredEntFiniteElementPtr()->getEnt();
    CHKERR this->setTagsData(tet, nb_gauss_pts, nb_internal_variables);
 
    commonDataPtr->activeVariablesW.resize(
        nb_gauss_pts, 12 + cpPtr->nbInternalVariables, false);
    commonDataPtr->gradientW.resize(nb_gauss_pts,
                                    12 + cpPtr->nbInternalVariables, false);
    commonDataPtr->materialTangentOperator.resize(36, nb_gauss_pts, false);
    commonDataPtr->deltaGamma.resize(nb_gauss_pts);
 
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    FTensor::Index<'k', 3> k;
    FTensor::Index<'l', 3> l;
    FTensor::Index<'Z', 6> Z;
    FTensor::Index<'Y', 6> Y;
 
    // Code uses trial stress, until first solve of linear system of equations
    int iter = 1;
    if (this->getFEMethod()->snes_ctx != SnesMethod::CTX_SNESNONE) {
      CHKERR SNESGetLinearSolveIterations(this->getFEMethod()->snes, &iter);
    }
 
    // b-bar
    auto ave_tr_strain = calcAveStrain(
        this->commonDataPtr->bBar, nb_gauss_pts,
        getFTensor2FromMat<DIM, DIM>(commonDataPtr->gradAtGaussPts),
        this->getFTensor0IntegrationWeight());
 
    auto t_voight_stress_op = voight_to_stress_op();
 
    auto t_grad = getFTensor2FromMat<DIM, DIM>(commonDataPtr->gradAtGaussPts);
 
    auto t_Cep = getFTensor4DdgFromMat<3, 3, 1>(
        this->commonDataPtr->materialTangentOperator);
    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
 
      auto tmp_active =
          getVectorAdaptor(&(commonDataPtr->activeVariablesW(gg, 0)),
                           12 + cpPtr->nbInternalVariables);
      cpPtr->activeVariablesW.swap(tmp_active);
      auto tmp_gradient = getVectorAdaptor(&(commonDataPtr->gradientW(gg, 0)),
                                           12 + cpPtr->nbInternalVariables);
      cpPtr->gradientW.swap(tmp_gradient);
 
      auto t_voigt_strain =
          calcStrain(commonDataPtr->bBar, ave_tr_strain, t_grad,
                     getFTensor1FromPtr<6>(&(cpPtr->activeVariablesW[6])));
      ++t_grad;
 
      // Copy plastic strains and internal variables from tags
      auto copy_plastic_strain_and_internal_variables = [&]() {
        MoFEMFunctionBegin;
 
        // Get data stored on Tags
        VectorAdaptor plastic_strain_tags =
            VectorAdaptor(6, ublas::shallow_array_adaptor<double>(
                                 6, &commonDataPtr->plasticStrainPtr[gg * 6]));
        VectorAdaptor internal_variables_tags = VectorAdaptor(
            nb_internal_variables,
            ublas::shallow_array_adaptor<double>(
                nb_internal_variables,
                &commonDataPtr
                     ->internalVariablesPtr[gg * nb_internal_variables]));
 
        // Set values to closest point projection
        cpPtr->plasticStrain0.resize(6, false);
        noalias(cpPtr->plasticStrain0) = plastic_strain_tags;
        cpPtr->internalVariables0.resize(nb_internal_variables, false);
        noalias(cpPtr->internalVariables0) = internal_variables_tags;
 
        auto cp_plastic_strain = cpPtr->getPlasticStrain();
        auto cp_internal_variables = cpPtr->getInternalVariables();
        noalias(cp_plastic_strain) = plastic_strain_tags;
        noalias(cp_internal_variables) = internal_variables_tags;
        MoFEMFunctionReturn(0);
      };
 
      CHKERR copy_plastic_strain_and_internal_variables();
 
      auto closest_point_projection = [&](auto &recalculate_elastic_tangent) {
        MoFEMFunctionBegin;
        for (auto t : {ClosestPointProjection::H, ClosestPointProjection::Y,
                       ClosestPointProjection::W}) {
          CHKERR cpPtr->setParams(t, recalculate_elastic_tangent);
        }
        CHKERR cpPtr->playW_NoHessian();
        CHKERR cpPtr->playY_NoGradient();
        double y = cpPtr->y;
        cpPtr->deltaGamma = 0; // zero lagrange multiplier
        if (iter > 0 && y > std::numeric_limits<double>::epsilon()) {
          CHKERR cpPtr->solveClosestProjection();
        }
        commonDataPtr->deltaGamma[gg] = cpPtr->deltaGamma;
        MoFEMFunctionReturn(0);
      };
 
      auto evaluate_consistent_tangent_matrix =
          [&](auto &recalculate_elastic_tangent) {
            MoFEMFunctionBegin;
            if (recalculate_elastic_tangent)
              CHKERR cpPtr->playW();
 
            if (iter > 0 &&
                cpPtr->deltaGamma > std::numeric_limits<double>::epsilon()) {
 
              CHKERR cpPtr->consistentTangent();
 
              auto &m = cpPtr->Cep;
              // for generic case of non-associated plasticity consistent
              // tangent matrix is non-symmetric
              auto t_voight_cep =
                  FTensor::Tensor2<FTensor::PackPtr<double *, 0>, 6, 6>(
                      &m(0, 0), &m(0, 1), &m(0, 2), &m(0, 3), &m(0, 4),
                      &m(0, 5),
 
                      &m(1, 0), &m(1, 1), &m(1, 2), &m(1, 3), &m(1, 4),
                      &m(1, 5),
 
                      &m(2, 0), &m(2, 1), &m(2, 2), &m(2, 3), &m(2, 4),
                      &m(2, 5),
 
                      &m(3, 0), &m(3, 1), &m(3, 2), &m(3, 3), &m(3, 4),
                      &m(3, 5),
 
                      &m(4, 0), &m(4, 1), &m(4, 2), &m(4, 3), &m(4, 4),
                      &m(4, 5),
 
                      &m(5, 0), &m(5, 1), &m(5, 2), &m(5, 3), &m(5, 4), &m(5, 5)
 
                  );
 
              t_Cep(i, j, k, l) =
                  t_voight_stress_op(i, j, Z) *
                  (t_voight_cep(Z, Y) * t_voight_stress_op(k, l, Y));
 
            } else {
 
              auto &m = cpPtr->C;
              auto t_voight_cep =
                  FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 0>, 6>(
                      &m(0, 0), &m(0, 1), &m(0, 2), &m(0, 3), &m(0, 4),
                      &m(0, 5),
 
                      &m(1, 1), &m(1, 2), &m(1, 3), &m(1, 4), &m(1, 5),
 
                      &m(2, 2), &m(2, 3), &m(2, 4), &m(2, 5),
 
                      &m(3, 3), &m(3, 4), &m(3, 5),
 
                      &m(4, 4), &m(4, 5),
 
                      &m(5, 5));
 
              t_Cep(i, j, k, l) =
                  t_voight_stress_op(i, j, Z) *
                  (t_voight_cep(Z, Y) * t_voight_stress_op(k, l, Y));
            }
 
            ++t_Cep;
            MoFEMFunctionReturn(0);
          };
 
      // Always recalculate elastic tangent if first element
      bool recalculate_elastic_tangent =
          (this->getNinTheLoop() == 0) ? true : false;
      CHKERR closest_point_projection(recalculate_elastic_tangent);
      if (this->calcLhs)
        CHKERR evaluate_consistent_tangent_matrix(recalculate_elastic_tangent);
    }
 
    auto calcs_stress_matrix = [&]() {
      MoFEMFunctionBegin;
      auto t_voight_stress_op = voight_to_stress_op();
      commonDataPtr->stressMatrix.resize(6, commonDataPtr->gradientW.size1(),
                                         false);
      auto t_stess =
          getFTensor2SymmetricFromMat<3>(commonDataPtr->stressMatrix);
      auto t_voight_stress = commonDataPtr->getFTensor1StressAtGaussPts();
      for (auto gg = 0; gg != commonDataPtr->gradientW.size1(); ++gg) {
        t_stess(i, j) = t_voight_stress_op(i, j, Z) * t_voight_stress(Z);
        ++t_voight_stress;
        ++t_stess;
      }
      MoFEMFunctionReturn(0);
    };
 
    auto calcs_plastic_strain_matrix = [&]() {
      MoFEMFunctionBegin;
      auto t_voight_strain_op = voight_to_stress_op();
      commonDataPtr->plasticStrainMatrix.resize(
          6, commonDataPtr->activeVariablesW.size1(), false);
      auto t_plastic_strain =
          getFTensor2SymmetricFromMat<3>(commonDataPtr->plasticStrainMatrix);
      auto t_voight_plastic_strain =
          commonDataPtr->getFTensor1PlasticStrainAtGaussPts();
      for (auto gg = 0; gg != commonDataPtr->activeVariablesW.size1(); ++gg) {
        t_plastic_strain(i, j) =
            t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
        ++t_voight_plastic_strain;
        ++t_plastic_strain;
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR calcs_stress_matrix();
    CHKERR calcs_plastic_strain_matrix();
    MoFEMFunctionReturn(0);
  };
 
  CHKERR do_work_impl(this->commonDataPtr, this->cpPtr);
 
  MoFEMFunctionReturn(0);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<3, LARGE_STRAIN>(
    MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs) {
  return new OpCalculateStress<3, LARGE_STRAIN>(m_field, common_data_ptr, cp_ptr,
                                     calc_lhs);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<3, SMALL_STRAIN>(
    MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs) {
  return new OpCalculateStress<3, SMALL_STRAIN>(m_field, common_data_ptr, cp_ptr,
                                     calc_lhs);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<2, LARGE_STRAIN>(
    MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs) {
  return new OpCalculateStress<2, LARGE_STRAIN>(m_field, common_data_ptr, cp_ptr,
                                     calc_lhs);
}
 
template <>
ForcesAndSourcesCore::UserDataOperator *getRawPtrOpCalculateStress<2, SMALL_STRAIN>(
    MoFEM::Interface &m_field, boost::shared_ptr<CommonData> common_data_ptr,
    boost::shared_ptr<ClosestPointProjection> cp_ptr, bool calc_lhs) {
  return new OpCalculateStress<2, SMALL_STRAIN>(m_field, common_data_ptr, cp_ptr,
                                     calc_lhs);
}
 
template <int DIM>
MoFEMErrorCode
opFactoryDomainUpdateImpl(MoFEM::Interface &m_field, Pip &pip,
                          std::string block_name,
                          boost::shared_ptr<CommonData> common_data_ptr,
                          boost::shared_ptr<ClosestPointProjection> cp_ptr) {
  MoFEMFunctionBegin;
  pip.push_back(new OpUpdate(common_data_ptr));
  MoFEMFunctionReturn(0);
};
 
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) {
  return opFactoryDomainUpdateImpl<3>(m_field, pip, block_name, common_data_ptr,
                                      cp_ptr);
};
 
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) {
  return opFactoryDomainUpdateImpl<2>(m_field, pip, block_name, common_data_ptr,
                                      cp_ptr);
};
 
 
//! [TSUpdateImpl]
/**
 * @brief TS update history variables
 */
struct TSUpdateImpl : public TSUpdate {
 
  TSUpdateImpl(std::string fe_name, boost::shared_ptr<FEMethod> fe_ptr);
  MoFEMErrorCode postProcess(TS ts);
 
private:
  boost::shared_ptr<FEMethod> fePtr;
  const std::string feName;
};
 
TSUpdateImpl::TSUpdateImpl(std::string fe_name,
                           boost::shared_ptr<FEMethod> fe_ptr)
    : feName(fe_name), fePtr(fe_ptr) {}
 
MoFEMErrorCode TSUpdateImpl::postProcess(TS ts) {
  MoFEMFunctionBegin;
  DM dm;
  CHKERR TSGetDM(ts, &dm);
  // This is where update is preformed. Element list is iterated and on each
  // element pipeline accessed throng fePtr is executed.
  CHKERR DMoFEMLoopFiniteElements(dm, feName, fePtr);
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Create test update object
 * 
 * @param fe_name 
 * @param fe_ptr 
 * @return boost::shared_ptr<TSUpdate> 
 */
boost::shared_ptr<TSUpdate> createTSUpdate(std::string fe_name,
                                           boost::shared_ptr<FEMethod> fe_ptr) {
  return boost::make_shared<TSUpdateImpl>(fe_name, fe_ptr);
}
//! [TSUpdateImpl]
 
} // namespace ADOLCPlasticity