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 dimessnion 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;
};
 
//! [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>
struct OpCalculateStress : public ForcesAndSourcesCore::UserDataOperator {
 
  OpCalculateStress(MoFEM::Interface &m_field,
                    boost::shared_ptr<CommonData> common_data_ptr,
                    boost::shared_ptr<ClosestPointProjection> cp_ptr,
                    bool calc_lhs);
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data);
 
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>
OpCalculateStress<DIM>::OpCalculateStress(
    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> MoFEMErrorCode OpCalculateStress<DIM>::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>
MoFEMErrorCode
OpCalculateStress<DIM>::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>::doWork(int side, EntityType type,
                               DataForcesAndSourcesCore::EntData &data) {
  MoFEMFunctionBegin;
 
  int nb_gauss_pts = getGaussPts().size2();
 
  int nb_internal_variables = cpPtr->nbInternalVariables;
  auto tet = getNumeredEntFiniteElementPtr()->getEnt();
  CHKERR 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 use trail stress, until first solve of linear system of equations
  int iter = 1;
  if (getFEMethod()->snes_ctx != SnesMethod::CTX_SNESNONE) {
    CHKERR SNESGetLinearSolveIterations(getFEMethod()->snes, &iter);
  }
 
  // b-bar
  double ave_tr_strain = calcAveStrain(
      commonDataPtr->bBar, nb_gauss_pts,
      getFTensor2FromMat<DIM, DIM>(commonDataPtr->gradAtGaussPts),
      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>(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->solveClosetProjection();
      }
      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 = (getNinTheLoop() == 0) ? true : false;
    CHKERR closest_point_projection(recalculate_elastic_tangent);
    if (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);
}
 
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) {
  return new OpCalculateStress<3>(m_field, common_data_ptr, cp_ptr, 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) {
  return new OpCalculateStress<2>(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;
  CHKERR cp_ptr->addMatBlockOps(m_field, pip, block_name, Sev::noisy);
  pip.push_back(
      getRawPtrOpCalculateStress<DIM>(m_field, common_data_ptr, cp_ptr, false));
  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