PlasticOperators.cpp

Go to the documentation of this file.

/** \file PlasticOperators.cpp
 */
 
/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
#include <MoFEM.hpp>
using namespace MoFEM;
using namespace FTensor;
 
#include <BasicFiniteElements.hpp>
#include <AnalyticalSurfaces.hpp>
#include <RigidBodies.hpp>
#include <BlockMatData.hpp>
#include <MultifieldPlasticity.hpp>
#include <BasicFeTools.hpp>
#include <MatrixFunction.hpp>
 
#include <ElasticOperators.hpp>
#include <PlasticOperators.hpp>
 
namespace OpPlasticTools {
 
MoFEMErrorCode CommonData::calculateDiffDeviator() {
  MoFEMFunctionBeginHot;
  auto diff_dev = diff_deviator(diff_tensor());
  this->diffDeviator(I, J, k, l) = diff_dev(I, J, k, l);
  MoFEMFunctionReturnHot(0);
}
 
OpCalculatePlasticSurfaceAndFlow::OpCalculatePlasticSurfaceAndFlow(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
 
MoFEMErrorCode OpCalculatePlasticSurfaceAndFlow::doWork(int side,
                                                        EntityType type,
                                                        EntData &data) {
  MoFEMFunctionBegin;
 
  const size_t nb_gauss_pts = commonDataPtr->mStressPtr->size2();
  auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
  auto &diff_dev = commonDataPtr->diffDeviator;
  commonDataPtr->plasticSurfacePtr->resize(nb_gauss_pts, false);
  commonDataPtr->plasticFlowPtr->resize(6, nb_gauss_pts, false);
  auto t_flow = getFTensor2SymmetricFromMat<3>(*commonDataPtr->plasticFlowPtr);
  auto t_plastic_strain =
      getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticStrainPtr));
 
  for (auto &f : *(commonDataPtr->plasticSurfacePtr)) {
 
    auto stress_eff = get_effective_stress(t_stress, t_plastic_strain);
    f = plastic_surface(deviator(stress_eff));
 
    auto t_flow_tmp = plastic_flow(f, deviator(stress_eff), diff_dev);
    t_flow(i, j) = t_flow_tmp(i, j);
 
    ++t_plastic_strain;
    ++t_flow;
    ++t_stress;
  }
 
  MoFEMFunctionReturn(0);
}
 
OpPlasticStress::OpPlasticStress(const std::string field_name,
                                 boost::shared_ptr<CommonData> common_data_ptr,
                                 boost::shared_ptr<MatrixDouble> mat_d_ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr), matDPtr(mat_d_ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
 
//! [Calculate stress]
MoFEMErrorCode OpPlasticStress::doWork(int side, EntityType type,
                                       EntData &data) {
  MoFEMFunctionBegin;
  const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
  commonDataPtr->mStressPtr->resize(6, nb_gauss_pts);
 
  auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*matDPtr);
  auto t_strain = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStrainPtr));
  auto t_plastic_strain =
      getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticStrainPtr));
 
  auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
    t_stress(i, j) =
        t_D(i, j, k, l) * (t_strain(k, l) - t_plastic_strain(k, l));
 
    ++t_strain;
    ++t_plastic_strain;
    ++t_stress;
  }
 
  MoFEMFunctionReturn(0);
}
//! [Calculate stress]
 
OpCalculatePlasticFlowRhs::OpCalculatePlasticFlowRhs(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(field_name, field_name, DomainEleOpAssembly::OPROW),
      commonDataPtr(common_data_ptr) {}
 
MoFEMErrorCode OpCalculatePlasticFlowRhs::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
  const size_t nb_dofs = data.getIndices().size();
  if (nb_dofs) {
 
    auto get_dt = [&]() {
      double dt;
      CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
      dt = 1;
      return dt;
    };
    const auto dt = get_dt();
 
    auto t_flow =
        getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
 
    auto t_plastic_strain_dot =
        getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticStrainDotPtr));
    auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
 
    auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
 
    const size_t nb_integration_pts = data.getN().size1();
    const size_t nb_base_functions = data.getN().size2();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_base = data.getFTensor0N();
 
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
      double alpha = dt * getMeasure() * t_w;
 
      Tensor2<PackPtr<double *, 6>, 3, 3> t_nf{&locF(0), &locF(1), &locF(2),
                                               &locF(1), &locF(3), &locF(4),
                                               &locF(2), &locF(4), &locF(5)};
 
      Tensor2_symmetric<double, 3> t_flow_stress_diff;
      t_flow_stress_diff(i, j) = t_D(i, j, k, l) * (t_plastic_strain_dot(k, l) -
                                                    t_tau_dot * t_flow(k, l));
 
      size_t bb = 0;
      for (; bb != nb_dofs / 6; ++bb) {
        t_nf(i, j) += alpha * t_base * t_flow_stress_diff(i, j);
        ++t_base;
        ++t_nf;
      }
      for (; bb < nb_base_functions; ++bb)
        ++t_base;
 
      ++t_flow;
      ++t_plastic_strain_dot;
      ++t_tau_dot;
      ++t_w;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpCalculateConstraintRhs::OpCalculateConstraintRhs(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(field_name, field_name, DomainEleOpAssembly::OPROW),
      commonDataPtr(common_data_ptr) {}
 
MoFEMErrorCode OpCalculateConstraintRhs::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  const size_t nb_dofs = data.getIndices().size();
  if (nb_dofs) {
 
    auto get_dt = [&]() {
      double dt;
      CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
      dt = 1;
      return dt;
    };
    const auto dt = get_dt();
 
    auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
    auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
    auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
    auto t_w = getFTensor0IntegrationWeight();
 
    auto t_base = data.getFTensor0N();
    const size_t nb_integration_pts = data.getN().size1();
    const size_t nb_base_functions = data.getN().size2();
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
      const double alpha = dt * getMeasure() * t_w * (*cache).scale_constraint;
      const double beta = alpha * constraint(t_tau_dot, t_f, hardening(t_tau));
 
      size_t bb = 0;
      for (; bb != nb_dofs; ++bb) {
        locF(bb) += beta * t_base;
        ++t_base;
      }
      for (; bb < nb_base_functions; ++bb)
        ++t_base;
 
      ++t_tau;
      ++t_tau_dot;
      ++t_f;
      ++t_w;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <bool LOGSTRAIN>
OpCalculatePlasticInternalForceLhs_dEP<LOGSTRAIN>::
    OpCalculatePlasticInternalForceLhs_dEP(
        const std::string row_field_name, const std::string col_field_name,
        boost::shared_ptr<CommonData> common_data_ptr,
        boost::shared_ptr<MatrixDouble> mat_D_ptr)
    : DomainEleOpAssembly(row_field_name, col_field_name,
                          DomainEleOpAssembly::OPROWCOL),
      commonDataPtr(common_data_ptr), matDPtr(mat_D_ptr) {
  sYmm = false;
}
template <bool LOGSTRAIN>
MoFEMErrorCode OpCalculatePlasticInternalForceLhs_dEP<LOGSTRAIN>::iNtegrate(
    EntData &row_data, EntData &col_data) {
  MoFEMFunctionBegin;
 
  const size_t nb_row_dofs = row_data.getIndices().size();
  const size_t nb_col_dofs = col_data.getIndices().size();
  if (nb_row_dofs && nb_col_dofs) {
 
    const size_t nb_integration_pts = row_data.getN().size1();
    const size_t nb_row_base_functions = row_data.getN().size2();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
    auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*matDPtr);
 
    MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
    auto t_dE_dF = getFTensor4FromMat<3, 3, 3, 3>(dE_dF);
    Index<'L', 6> L;
    auto t_L = symm_L_tensor();
    auto t_diff_symmetrize = diff_symmetrize();
 
    Tensor3<double, 3, 3, 6> t_DLs;
    Dg<double, 3, 6> t_DL;
    if (!LOGSTRAIN) {
      t_DL(i, j, L) = 0;
      for (int ii = 0; ii != 3; ++ii)
        for (int jj = ii; jj != 3; ++jj)
          for (int kk = 0; kk != 3; ++kk)
            for (int ll = 0; ll != 3; ++ll)
              for (int LL = 0; LL != 6; ++LL)
                t_DL(ii, jj, LL) += t_D(ii, jj, kk, ll) * t_L(kk, ll, LL);
    }
 
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
      double alpha = getMeasure() * t_w;
 
      if (LOGSTRAIN) {
        Tensor4<double, 3, 3, 3, 3> D2_symm;
        // D2_symm(i, j, k, l) = 0.;
 
        // FIXME: that transpose is not necessary
        // for (int i = 0; i != 3; ++i)
        //   for (int jj = 0; jj != 3; ++jj)
        //     for (int k = 0; k != 3; ++k)
        //       for (int l = 0; l != 3; ++l)
        //         for (int m = 0; m != 3; ++m)
        //           for (int n = 0; n != 3; ++n)
        //             D2_symm(k, l, m, n) +=
        //                 t_D(i, jj, m, n) * t_dE_dF(i, jj, k, l);
        D2_symm(k, l, m, n) = t_D(m, n, i, j) * t_dE_dF(i, j, k, l);
 
        t_DLs(i, j, L) = 0;
        for (int ii = 0; ii != 3; ++ii)
          for (int jj = 0; jj != 3; ++jj)
            for (int kk = 0; kk != 3; ++kk)
              for (int ll = 0; ll != 3; ++ll)
                for (int LL = 0; LL != 6; ++LL)
                  t_DLs(ii, jj, LL) +=
                      D2_symm(ii, jj, kk, ll) * t_L(kk, ll, LL);
      }
 
      size_t rr = 0;
      for (; rr != nb_row_dofs / 3; ++rr) {
 
        Tensor2<PackPtr<double *, 6>, 3, 6> t_mat{
            &locMat(3 * rr + 0, 0), &locMat(3 * rr + 0, 1),
            &locMat(3 * rr + 0, 2), &locMat(3 * rr + 0, 3),
            &locMat(3 * rr + 0, 4), &locMat(3 * rr + 0, 5),
            &locMat(3 * rr + 1, 0), &locMat(3 * rr + 1, 1),
            &locMat(3 * rr + 1, 2), &locMat(3 * rr + 1, 3),
            &locMat(3 * rr + 1, 4), &locMat(3 * rr + 1, 5),
            &locMat(3 * rr + 2, 0), &locMat(3 * rr + 2, 1),
            &locMat(3 * rr + 2, 2), &locMat(3 * rr + 2, 3),
            &locMat(3 * rr + 2, 4), &locMat(3 * rr + 2, 5)};
 
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        Tensor2<double, 3, 6> t_tmp;
 
        if (LOGSTRAIN) {
          t_tmp(i, L) = (t_DLs(i, j, L) * (alpha * t_row_diff_base(j)));
 
        } else
          t_tmp(i, L) = (t_DL(i, j, L) * (alpha * t_row_diff_base(j)));
 
        for (size_t cc = 0; cc != nb_col_dofs / 6; ++cc) {
 
          t_mat(i, L) -= (t_col_base * t_tmp(i, L));
 
          ++t_mat;
          ++t_col_base;
        }
 
        ++t_row_diff_base;
      }
 
      if (LOGSTRAIN)
        ++t_dE_dF;
 
      for (; rr < nb_row_base_functions; ++rr)
        ++t_row_diff_base;
      ++t_w;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <bool LOGSTRAIN>
OpCalculatePlasticFlowLhs_dU<LOGSTRAIN>::OpCalculatePlasticFlowLhs_dU(
    const std::string row_field_name, const std::string col_field_name,
    boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(row_field_name, col_field_name,
                          DomainEleOpAssembly::OPROWCOL),
      commonDataPtr(common_data_ptr) {
  sYmm = false;
}
template <bool LOGSTRAIN>
MoFEMErrorCode OpCalculatePlasticFlowLhs_dU<LOGSTRAIN>::iNtegrate(
 
    EntData &row_data, EntData &col_data) {
  MoFEMFunctionBegin;
 
  const size_t nb_row_dofs = row_data.getIndices().size();
  const size_t nb_col_dofs = col_data.getIndices().size();
  if (nb_row_dofs && nb_col_dofs) {
 
    auto get_dt = [&]() {
      double dt;
      CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
      dt = 1;
      return dt;
    };
    const auto dt = get_dt();
 
    const size_t nb_integration_pts = row_data.getN().size1();
    const size_t nb_row_base_functions = row_data.getN().size2();
    auto t_w = getFTensor0IntegrationWeight();
 
    auto get_row_base = [&]() {
      if (commonDataPtr->isDualBase) {
        double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
        return Tensor0<PackPtr<double *, 1>>(base_ptr);
      } else {
        return row_data.getFTensor0N();
      }
    };
    auto t_row_base = get_row_base();
 
    auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
    auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
    auto &diff_dev = commonDataPtr->diffDeviator;
    auto t_flow =
        getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
    auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
    auto t_D_Deviator =
        getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
 
    MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
    auto t_dE_dF = getFTensor4FromMat<3, 3, 3, 3>(dE_dF);
 
    auto t_diff_symmetrize = diff_symmetrize();
 
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
 
      double alpha = dt * getMeasure() * t_w;
      auto t_diff_plastic_flow_dstrain = diff_plastic_flow_dstrain(
          t_D_Deviator, diff_plastic_flow_dstress(t_f, t_flow, diff_dev));
      Ddg<double, 3, 3> t_flow_stress_dstrain;
      t_flow_stress_dstrain(i, j, k, l) =
          t_D(i, j, m, n) * t_diff_plastic_flow_dstrain(m, n, k, l);
      Tensor4<double, 3, 3, 3, 3> t_diff_plastic_flow_stress_dgrad;
 
      if (!LOGSTRAIN)
        t_diff_plastic_flow_stress_dgrad(i, j, k, l) =
            t_flow_stress_dstrain(i, j, m, n) * t_diff_symmetrize(m, n, k, l);
      else {
 
        Tensor4<double, 3, 3, 3, 3> t_diff_plastic_flow_stress_dgrads;
        t_diff_plastic_flow_stress_dgrads(i, j, k, l) =
            t_flow_stress_dstrain(i, j, m, n) * t_diff_symmetrize(m, n, k, l);
        t_diff_plastic_flow_stress_dgrad(i, j, k, l) = 0.;
        for (int i = 0; i != 3; ++i)
          for (int j = 0; j != 3; ++j)
            for (int k = 0; k != 3; ++k)
              for (int l = 0; l != 3; ++l)
                for (int m = 0; m != 3; ++m)
                  for (int n = 0; n != 3; ++n)
                    t_diff_plastic_flow_stress_dgrad(i, j, k, l) +=
                        t_diff_plastic_flow_stress_dgrads(i, j, m, n) *
                        t_dE_dF(m, n, k, l);
      }
 
      size_t rr = 0;
      for (; rr != nb_row_dofs / 6; ++rr) {
 
        // t3(T d000, T d001, T d002, T d010, T d011, T d012, T d020, T d021,
        //   T d022, T d100, T d101, T d102, T d110, T d111, T d112, T d120,
        //   T d121, T d122, T d200, T d201, T d202, T d210, T d211, T d212,
        //   T d220, T d221, T d222);
 
        Tensor3<PackPtr<double *, 3>, 3, 3, 3> t_mat{
            &locMat(6 * rr + 0, 0), // 000
            &locMat(6 * rr + 0, 1), // 001
            &locMat(6 * rr + 0, 2), // 002
 
            &locMat(6 * rr + 1, 0), // 010
            &locMat(6 * rr + 1, 1), // 011
            &locMat(6 * rr + 1, 2), // 012
 
            &locMat(6 * rr + 2, 0), // 020
            &locMat(6 * rr + 2, 1), // 021
            &locMat(6 * rr + 2, 2), // 022
 
            &locMat(6 * rr + 1, 0), // 100
            &locMat(6 * rr + 1, 1), // 101
            &locMat(6 * rr + 1, 2), // 102
 
            &locMat(6 * rr + 3, 0), // 110
            &locMat(6 * rr + 3, 1), // 111
            &locMat(6 * rr + 3, 2), // 112
 
            &locMat(6 * rr + 4, 0), // 120
            &locMat(6 * rr + 4, 1), // 121
            &locMat(6 * rr + 4, 2), // 122
 
            &locMat(6 * rr + 2, 0), // 200
            &locMat(6 * rr + 2, 1), // 201
            &locMat(6 * rr + 2, 2), // 202
 
            &locMat(6 * rr + 4, 0), // 210
            &locMat(6 * rr + 4, 1), // 211
            &locMat(6 * rr + 4, 2), // 212
 
            &locMat(6 * rr + 5, 0), // 220
            &locMat(6 * rr + 5, 1), // 221
            &locMat(6 * rr + 5, 2)  // 222
        };
 
        const double c0 = alpha * t_row_base * t_tau_dot;
 
        auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
        for (size_t cc = 0; cc != nb_col_dofs / 3; ++cc) {
 
          t_mat(i, j, l) -= c0 * (t_diff_plastic_flow_stress_dgrad(i, j, l, k) *
                                  t_col_diff_base(k));
 
          ++t_mat;
          ++t_col_diff_base;
        }
 
        ++t_row_base;
      }
      for (; rr < nb_row_base_functions; ++rr)
        ++t_row_base;
 
      if (LOGSTRAIN)
        ++t_dE_dF;
 
      ++t_w;
      ++t_f;
      ++t_flow;
      ++t_tau_dot;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpCalculatePlasticFlowLhs_dEP::OpCalculatePlasticFlowLhs_dEP(
    const std::string row_field_name, const std::string col_field_name,
    boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(row_field_name, col_field_name,
                          DomainEleOpAssembly::OPROWCOL),
      commonDataPtr(common_data_ptr) {
  sYmm = false;
}
 
MoFEMErrorCode OpCalculatePlasticFlowLhs_dEP::iNtegrate(EntData &row_data,
                                                        EntData &col_data) {
  MoFEMFunctionBegin;
 
  const size_t nb_row_dofs = row_data.getIndices().size();
  const size_t nb_col_dofs = col_data.getIndices().size();
  if (nb_row_dofs && nb_col_dofs) {
 
    auto get_dt = [&]() {
      double dt;
      CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
      dt = 1;
      return dt;
    };
    const auto dt = get_dt();
 
    const size_t nb_integration_pts = row_data.getN().size1();
    const size_t nb_row_base_functions = row_data.getN().size2();
    auto t_w = getFTensor0IntegrationWeight();
 
    auto get_row_base = [&]() {
      if (commonDataPtr->isDualBase) {
        double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
        return Tensor0<PackPtr<double *, 1>>(base_ptr);
      } else {
        return row_data.getFTensor0N();
      }
    };
    auto t_row_base = get_row_base();
 
    auto &diff_dev = commonDataPtr->diffDeviator;
    auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
    auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
    auto t_flow =
        getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
 
    auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
    auto t_D_Deviator =
        getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
 
    auto t_omega = getFTensor1FromMat<3>(*commonDataPtr->guidingVelocityPtr);
    bool is_rotating = commonDataPtr->guidingVelocityPtr->size2() > 1;
 
    auto t_diff_plastic_strain = diff_tensor();
 
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
      double alpha = dt * getMeasure() * t_w;
      double c1 = alpha * getTSa();
      const double c0 = alpha * t_tau_dot;
 
 
      auto t_diff_plastic_flow_dstrain = diff_plastic_flow_kin_hard_dstrain(
          t_D_Deviator, diff_plastic_flow_dstress(t_f, t_flow, diff_dev));
      Ddg<double, 3, 3> my_tensor;
      my_tensor(i, j, k, l) =
          (t_D(i, j, m, n) * (c1 * t_diff_plastic_strain(m, n, k, l) +
                              c0 * t_diff_plastic_flow_dstrain(m, n, k, l)));
 
      size_t rr = 0;
      for (; rr != nb_row_dofs / 6; ++rr) {
 
        // t4(T d0000, T d0001, T d0002, T d0010, T d0011, T d0012, T d0020,
        //    T d0021, T d0022, T d0100, T d0101, T d0102, T d0110, T d0111,
        //    T d0112, T d0120, T d0121, T d0122, T d0200, T d0201, T d0202,
        //    T d0210, T d0211, T d0212, T d0220, T d0221, T d0222, T d1000,
        //    T d1001, T d1002, T d1010, T d1011, T d1012, T d1020, T d1021,
        //    T d1022, T d1100, T d1101, T d1102, T d1110, T d1111, T d1112,
        //    T d1120, T d1121, T d1122, T d1200, T d1201, T d1202, T d1210,
        //    T d1211, T d1212, T d1220, T d1221, T d1222, T d2000, T d2001,
        //    T d2002, T d2010, T d2011, T d2012, T d2020, T d2021, T d2022,
        //    T d2100, T d2101, T d2102, T d2110, T d2111, T d2112, T d2120,
        //    T d2121, T d2122, T d2200, T d2201, T d2202, T d2210, T d2211,
        //    T d2212, T d2220, T d2221, T d2222);
 
        Tensor4<PackPtr<double *, 6>, 3, 3, 3, 3> t_mat{
 
            &locMat(6 * rr + 0, 0), // 0000
            &locMat(6 * rr + 0, 1), // 0001
            &locMat(6 * rr + 0, 2), // 0002
            &locMat(6 * rr + 0, 1), // 0010
            &locMat(6 * rr + 0, 3), // 0011
            &locMat(6 * rr + 0, 4), // 0012
            &locMat(6 * rr + 0, 2), // 0020
            &locMat(6 * rr + 0, 4), // 0021
            &locMat(6 * rr + 0, 5), // 0022
 
            &locMat(6 * rr + 1, 0), // 0100
            &locMat(6 * rr + 1, 1), // 0101
            &locMat(6 * rr + 1, 2), // 0102
            &locMat(6 * rr + 1, 1), // 0110
            &locMat(6 * rr + 1, 3), // 0111
            &locMat(6 * rr + 1, 4), // 0112
            &locMat(6 * rr + 1, 2), // 0120
            &locMat(6 * rr + 1, 4), // 0121
            &locMat(6 * rr + 1, 5), // 0122
 
            &locMat(6 * rr + 2, 0), // 0200
            &locMat(6 * rr + 2, 1), // 0201
            &locMat(6 * rr + 2, 2), // 0202
            &locMat(6 * rr + 2, 1), // 0210
            &locMat(6 * rr + 2, 3), // 0211
            &locMat(6 * rr + 2, 4), // 0212
            &locMat(6 * rr + 2, 2), // 0220
            &locMat(6 * rr + 2, 4), // 0221
            &locMat(6 * rr + 2, 5), // 0222
 
            &locMat(6 * rr + 1, 0), // 1000
            &locMat(6 * rr + 1, 1), // 1001
            &locMat(6 * rr + 1, 2), // 1002
            &locMat(6 * rr + 1, 1), // 1010
            &locMat(6 * rr + 1, 3), // 1011
            &locMat(6 * rr + 1, 4), // 1012
            &locMat(6 * rr + 1, 2), // 1020
            &locMat(6 * rr + 1, 4), // 1021
            &locMat(6 * rr + 1, 5), // 1022
 
            &locMat(6 * rr + 3, 0), // 1100
            &locMat(6 * rr + 3, 1), // 1101
            &locMat(6 * rr + 3, 2), // 1102
            &locMat(6 * rr + 3, 1), // 1110
            &locMat(6 * rr + 3, 3), // 1111
            &locMat(6 * rr + 3, 4), // 1112
            &locMat(6 * rr + 3, 2), // 1120
            &locMat(6 * rr + 3, 4), // 1121
            &locMat(6 * rr + 3, 5), // 1122
 
            &locMat(6 * rr + 4, 0), // 1200
            &locMat(6 * rr + 4, 1), // 1201
            &locMat(6 * rr + 4, 2), // 1202
            &locMat(6 * rr + 4, 1), // 1210
            &locMat(6 * rr + 4, 3), // 1211
            &locMat(6 * rr + 4, 4), // 1212
            &locMat(6 * rr + 4, 2), // 1220
            &locMat(6 * rr + 4, 4), // 1221
            &locMat(6 * rr + 4, 5), // 1222
 
            &locMat(6 * rr + 2, 0), // 2000
            &locMat(6 * rr + 2, 1), // 2001
            &locMat(6 * rr + 2, 2), // 2002
            &locMat(6 * rr + 2, 1), // 2010
            &locMat(6 * rr + 2, 3), // 2011
            &locMat(6 * rr + 2, 4), // 2012
            &locMat(6 * rr + 2, 2), // 2020
            &locMat(6 * rr + 2, 4), // 2021
            &locMat(6 * rr + 2, 5), // 2022
 
            &locMat(6 * rr + 4, 0), // 2100
            &locMat(6 * rr + 4, 1), // 2101
            &locMat(6 * rr + 4, 2), // 2102
            &locMat(6 * rr + 4, 1), // 2110
            &locMat(6 * rr + 4, 3), // 2111
            &locMat(6 * rr + 4, 4), // 2112
            &locMat(6 * rr + 4, 2), // 2120
            &locMat(6 * rr + 4, 4), // 2121
            &locMat(6 * rr + 4, 5), // 2122
 
            &locMat(6 * rr + 5, 0), // 2200
            &locMat(6 * rr + 5, 1), // 2201
            &locMat(6 * rr + 5, 2), // 2202
            &locMat(6 * rr + 5, 1), // 2210
            &locMat(6 * rr + 5, 3), // 2211
            &locMat(6 * rr + 5, 4), // 2212
            &locMat(6 * rr + 5, 2), // 2220
            &locMat(6 * rr + 5, 4), // 2221
            &locMat(6 * rr + 5, 5)  // 2222
 
        };
 
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
        for (size_t cc = 0; cc != nb_col_dofs / 6; ++cc) {
          t_mat(i, j, k, l) += my_tensor(i, j, k, l) * t_col_base * t_row_base;
 
          if (is_rotating)
            t_mat(i, j, k, l) +=
                t_row_base *
                (t_D(i, j, m, n) * (alpha * t_diff_plastic_strain(m, n, k, l) *
                                    (t_col_diff_base(i) * t_omega(i))));
 
          ++t_mat;
          ++t_col_diff_base;
          ++t_col_base;
        }
 
        ++t_row_base;
      }
      for (; rr < nb_row_base_functions; ++rr)
        ++t_row_base;
 
      if (is_rotating)
        ++t_omega;
 
      ++t_w;
      ++t_f;
      ++t_flow;
      ++t_tau_dot;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpCalculatePlasticFlowLhs_dTAU::OpCalculatePlasticFlowLhs_dTAU(
    const std::string row_field_name, const std::string col_field_name,
    boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(row_field_name, col_field_name,
                          DomainEleOp::OPROWCOL),
      commonDataPtr(common_data_ptr) {
  sYmm = false;
}
 
MoFEMErrorCode OpCalculatePlasticFlowLhs_dTAU::iNtegrate(EntData &row_data,
                                                         EntData &col_data) {
  MoFEMFunctionBegin;
 
  const size_t nb_row_dofs = row_data.getIndices().size();
  const size_t nb_col_dofs = col_data.getIndices().size();
  if (nb_row_dofs && nb_col_dofs) {
 
    auto get_dt = [&]() {
      double dt;
      CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
      dt = 1;
      return dt;
    };
    const auto dt = get_dt();
 
    const size_t nb_integration_pts = row_data.getN().size1();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_flow =
        getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
 
    auto get_row_base = [&]() {
      if (commonDataPtr->isDualBase) {
        double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
        return Tensor0<PackPtr<double *, 1>>(base_ptr);
      } else {
        return row_data.getFTensor0N();
      }
    };
    auto t_row_base = get_row_base();
 
    auto t_omega = getFTensor1FromMat<3>(*commonDataPtr->guidingVelocityPtr);
    bool is_rotating = commonDataPtr->guidingVelocityPtr->size2() > 1;
 
    auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
 
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
      double alpha = dt * t_w * getMeasure();
      double beta = alpha * getTSa();
 
      Tensor2_symmetric<double, 3> t_flow_stress;
      t_flow_stress(i, j) = t_D(i, j, m, n) * t_flow(m, n);
 
      for (size_t rr = 0; rr != nb_row_dofs / 6; ++rr) {
 
        Tensor2<PackPtr<double *, 1>, 3, 3> t_mat{
            &locMat(6 * rr + 0, 0), &locMat(6 * rr + 1, 0),
            &locMat(6 * rr + 2, 0), &locMat(6 * rr + 1, 0),
            &locMat(6 * rr + 3, 0), &locMat(6 * rr + 4, 0),
            &locMat(6 * rr + 2, 0), &locMat(6 * rr + 4, 0),
            &locMat(6 * rr + 5, 0)};
 
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
        for (size_t cc = 0; cc != nb_col_dofs; cc++) {
          t_mat(i, j) -= beta * t_row_base * t_col_base * t_flow_stress(i, j);
          if (is_rotating)
            t_mat(i, j) -= alpha * t_row_base *
                           t_flow_stress(i, j) *
                           (t_col_diff_base(k) * t_omega(k));
          ++t_mat;
          ++t_col_base;
          ++t_col_diff_base;
        }
 
        ++t_row_base;
      }
      if (is_rotating)
        ++t_omega;
      ++t_w;
      ++t_flow;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpCalculatePlasticInternalForceLhs_dTAU::
    OpCalculatePlasticInternalForceLhs_dTAU(
        const std::string row_field_name, const std::string col_field_name,
        boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(row_field_name, col_field_name,
                          DomainEleOp::OPROWCOL),
      commonDataPtr(common_data_ptr) {
  sYmm = false;
}
 
MoFEMErrorCode
OpCalculatePlasticInternalForceLhs_dTAU::iNtegrate(EntData &row_data,
                                                   EntData &col_data) {
  return 0;
}
 
template <bool LOGSTRAIN>
OpCalculateConstraintLhs_dU<LOGSTRAIN>::OpCalculateConstraintLhs_dU(
    const std::string row_field_name, const std::string col_field_name,
    boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(row_field_name, col_field_name,
                          DomainEleOp::OPROWCOL),
      commonDataPtr(common_data_ptr) {
  sYmm = false;
}
template <bool LOGSTRAIN>
MoFEMErrorCode OpCalculateConstraintLhs_dU<LOGSTRAIN>::iNtegrate(
 
    EntData &row_data, EntData &col_data) {
  MoFEMFunctionBegin;
 
  const size_t nb_row_dofs = row_data.getIndices().size();
  const size_t nb_col_dofs = col_data.getIndices().size();
  if (nb_row_dofs && nb_col_dofs) {
 
    auto get_dt = [&]() {
      double dt;
      CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
      dt = 1;
      return dt;
    };
    const auto dt = get_dt();
 
    const size_t nb_integration_pts = row_data.getN().size1();
    const size_t nb_row_base_functions = row_data.getN().size2();
    auto t_w = getFTensor0IntegrationWeight();
 
    auto get_row_base = [&]() {
      if (commonDataPtr->isDualBase) {
        double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
        return Tensor0<PackPtr<double *, 1>>(base_ptr);
      } else {
        return row_data.getFTensor0N();
      }
    };
    auto t_row_base = get_row_base();
 
    auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
    auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
    auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
    auto t_flow =
        getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
    auto t_stress =
        getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
    auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
    auto t_D_Deviator =
        getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
    auto t_diff_symmetrize = diff_symmetrize();
 
    MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
    auto t_dE_dF = getFTensor4FromMat<3, 3, 3, 3>(dE_dF);
 
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
      double alpha = dt * getMeasure() * t_w * (*cache).scale_constraint;
 
      auto t_diff_constrain_dstrain = diff_constrain_dstrain(
          t_D_Deviator,
          diff_constrain_dstress(
              diff_constrain_df(t_tau_dot, t_f, hardening(t_tau)), t_flow));
 
      Tensor2<double, 3, 3> t_diff_constrain_dgrad;
      if (!LOGSTRAIN) {
 
        t_diff_constrain_dgrad(k, l) =
            t_diff_constrain_dstrain(i, j) * t_diff_symmetrize(i, j, k, l);
      } else {
 
        Tensor2<double, 3, 3> t_diff_constrain_dgrads;
        t_diff_constrain_dgrad(k, l) =
            t_diff_constrain_dstrain(i, j) * t_dE_dF(i, j, k, l);
      }
 
      Tensor1<PackPtr<double *, 3>, 3> t_mat{&locMat(0, 0), &locMat(0, 1),
                                             &locMat(0, 2)};
 
      size_t rr = 0;
      for (; rr != nb_row_dofs; ++rr) {
 
        auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
        for (size_t cc = 0; cc != nb_col_dofs / 3; cc++) {
 
          t_mat(i) += alpha * t_row_base * t_diff_constrain_dgrad(i, j) *
                      t_col_diff_base(j);
 
          ++t_mat;
          ++t_col_diff_base;
        }
 
        ++t_row_base;
      }
      for (; rr != nb_row_base_functions; ++rr)
        ++t_row_base;
 
      if (LOGSTRAIN)
        ++t_dE_dF;
 
      ++t_f;
      ++t_tau;
      ++t_tau_dot;
      ++t_flow;
      ++t_stress;
      ++t_w;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpCalculateConstraintLhs_dEP::OpCalculateConstraintLhs_dEP(
    const std::string row_field_name, const std::string col_field_name,
    boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(row_field_name, col_field_name,
                          DomainEleOp::OPROWCOL),
      commonDataPtr(common_data_ptr) {
  sYmm = false;
}
 
MoFEMErrorCode OpCalculateConstraintLhs_dEP::iNtegrate(EntData &row_data,
                                                       EntData &col_data) {
  MoFEMFunctionBegin;
 
  const size_t nb_row_dofs = row_data.getIndices().size();
  const size_t nb_col_dofs = col_data.getIndices().size();
  if (nb_row_dofs && nb_col_dofs) {
 
    auto get_dt = [&]() {
      double dt;
      CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
      dt = 1;
      return dt;
    };
    const auto dt = get_dt();
 
    const size_t nb_integration_pts = row_data.getN().size1();
    const size_t nb_row_base_functions = row_data.getN().size2();
    auto t_w = getFTensor0IntegrationWeight();
 
    auto get_row_base = [&]() {
      if (commonDataPtr->isDualBase) {
        double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
        return Tensor0<PackPtr<double *, 1>>(base_ptr);
      } else {
        return row_data.getFTensor0N();
      }
    };
    auto t_row_base = get_row_base();
 
    auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
    auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
    auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
    auto t_flow =
        getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
 
    auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD);
    auto t_D_Deviator =
        getFTensor4DdgFromMat<3, 3, 0>(*commonDataPtr->mtD_Deviator);
 
    auto t_diff_symmetrize = diff_symmetrize();
    Index<'L', 6> L;
    auto t_L = symm_L_tensor();
 
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
      double alpha = dt * getMeasure() * t_w * (*cache).scale_constraint;
 
      auto mat_ptr = locMat.data().begin();
      auto t_diff_constrain_dstrain = diff_constrain_kin_hard_dstrain(
          t_D_Deviator,
          diff_constrain_dstress(
              diff_constrain_df(t_tau_dot, t_f, hardening(t_tau)), t_flow));
 
      Tensor1<PackPtr<double *, 6>, 6> t_mat(&locMat(0, 0), &locMat(0, 1),
                                             &locMat(0, 2), &locMat(0, 3),
                                             &locMat(0, 4), &locMat(0, 5));
 
      size_t rr = 0;
      for (; rr != nb_row_dofs; ++rr) {
 
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        for (size_t cc = 0; cc != nb_col_dofs / 6; cc++) {
 
          t_mat(L) -= (alpha * t_row_base * t_col_base) *
                      (t_diff_constrain_dstrain(i, j) * t_L(i, j, L));
 
          ++t_mat;
          ++t_col_base;
        }
 
        ++t_row_base;
      }
      for (; rr != nb_row_base_functions; ++rr)
        ++t_row_base;
 
      ++t_f;
      ++t_tau;
      ++t_tau_dot;
      ++t_flow;
      ++t_w;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpCalculateConstraintLhs_dTAU::OpCalculateConstraintLhs_dTAU(
    const std::string row_field_name, const std::string col_field_name,
    boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOpAssembly(row_field_name, col_field_name,
                          DomainEleOp::OPROWCOL),
      commonDataPtr(common_data_ptr) {
  sYmm = false;
}
 
MoFEMErrorCode OpCalculateConstraintLhs_dTAU::iNtegrate(EntData &row_data,
                                                        EntData &col_data) {
  MoFEMFunctionBegin;
 
  const size_t nb_row_dofs = row_data.getIndices().size();
  const size_t nb_col_dofs = col_data.getIndices().size();
  if (nb_row_dofs && nb_col_dofs) {
 
    auto get_dt = [&]() {
      double dt;
      CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
      dt = 1;
      return dt;
    };
    const auto dt = get_dt();
 
    const size_t nb_integration_pts = row_data.getN().size1();
    const size_t nb_row_base_functions = row_data.getN().size2();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
    auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
    auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
 
    auto t_omega = getFTensor1FromMat<3>(*commonDataPtr->guidingVelocityPtr);
    bool is_rotating = commonDataPtr->guidingVelocityPtr->size2() > 1;
 
    auto &cn = (*cache).cn_pl;
 
    auto get_row_base = [&]() {
      if (commonDataPtr->isDualBase) {
        double *base_ptr = &*commonDataPtr->dualBaseMat.data().begin();
        return Tensor0<PackPtr<double *, 1>>(base_ptr);
      } else {
        return row_data.getFTensor0N();
      }
    };
    auto t_row_base = get_row_base();
 
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
      const double alpha = dt * getMeasure() * t_w * (*cache).scale_constraint;
      const double t_a = getTSa();
      const double c0 =
          alpha * /* t_a * */
          diff_constrain_dtau_dot(t_tau_dot, t_f, hardening(t_tau));
      const double c1 =
          alpha * diff_constrain_dsigma_y(t_tau_dot, t_f, hardening(t_tau)) *
          hardening_dtau(t_tau);
 
      auto mat_ptr = locMat.data().begin();
 
      size_t rr = 0;
      for (; rr != nb_row_dofs; ++rr) {
 
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
        for (size_t cc = 0; cc != nb_col_dofs; ++cc) {
          if (!is_rotating) {
            *mat_ptr += (c0 * t_a + c1) * t_row_base * t_col_base;
 
          } else {
            const double c0_p =
                c0 * (t_col_base * t_a + (t_col_diff_base(i) * t_omega(i)));
            *mat_ptr += (c0_p + c1 * t_col_base) * t_row_base;
          }
 
          ++mat_ptr;
          ++t_col_base;
          ++t_col_diff_base;
        }
        ++t_row_base;
      }
      for (; rr < nb_row_base_functions; ++rr)
        ++t_row_base;
 
      if (is_rotating)
        ++t_omega;
 
      ++t_w;
      ++t_f;
      ++t_tau;
      ++t_tau_dot;
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
OpPostProcPlastic::OpPostProcPlastic(
    const std::string field_name, moab::Interface &post_proc_mesh,
    std::vector<EntityHandle> &map_gauss_pts,
    boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
      commonDataPtr(common_data_ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
 
//! [Postprocessing]
MoFEMErrorCode OpPostProcPlastic::doWork(int side, EntityType type,
                                         EntData &data) {
  MoFEMFunctionBegin;
 
  std::array<double, 9> def;
  std::fill(def.begin(), def.end(), 0);
 
  auto get_tag = [&](const std::string name, size_t size) {
    Tag th;
    CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
                                       MB_TAG_CREAT | MB_TAG_SPARSE,
                                       def.data());
    return th;
  };
 
  MatrixDouble3by3 mat(3, 3);
 
  auto set_matrix_3d = [&](auto &t) -> MatrixDouble3by3 & {
    mat.clear();
    for (size_t r = 0; r != 3; ++r)
      for (size_t c = 0; c != 3; ++c)
        mat(r, c) = t(r, c);
    return mat;
  };
 
  auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
    mat.clear();
    mat(0, 0) = t;
    return mat;
  };
 
  auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
    return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
                                     &*mat.data().begin());
  };
 
  auto th_plastic_surface = get_tag("PLASTIC_SURFACE", 1);
  auto th_tau = get_tag("PLASTIC_MULTIPLIER", 1);
 
  auto th_eqv_ep = get_tag("EQUIV_EP", 1);
  auto th_hardening = get_tag("HARDENING", 1);
 
  auto th_plastic_flow = get_tag("PLASTIC_FLOW", 9);
  string strain_name = "PLASTIC_STRAIN";
 
  auto th_plastic_strain = get_tag(strain_name, 9);
 
  auto t_flow =
      getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticFlowPtr));
  auto t_plastic_strain =
      getFTensor2SymmetricFromMat<3>(*(commonDataPtr->plasticStrainPtr));
  auto t_stress = getFTensor2SymmetricFromMat<3>(*(commonDataPtr->mStressPtr));
 
  size_t gg = 0;
  for (; gg != commonDataPtr->plasticSurfacePtr->size(); ++gg) {
    const double f = (*(commonDataPtr->plasticSurfacePtr))[gg];
    const double tau = (*(commonDataPtr->plasticTauPtr))[gg];
 
    auto plast_surface = f / (*cache).young_modulus_inv;
    CHKERR set_tag(th_plastic_surface, gg, set_scalar(plast_surface));
    CHKERR set_tag(th_tau, gg, set_scalar(tau));
    // FIXME: add option for more postprocessing
    // CHKERR set_tag(th_plastic_flow, gg, set_matrix_3d(t_flow));
 
    // t_plastic_strain(i, j) *= 2.;
    // auto exp_plast_strain = exp_map(t_plastic_strain);
    // double det_g;
    // CHKERR determinantTensor3by3(exp_plast_strain, det_g);
    // CHKERR set_tag(th_det_g, gg, set_scalar(det_g));
 
    CHKERR set_tag(th_plastic_strain, gg, set_matrix_3d(t_plastic_strain));
    const double eqv_ep =
        std::sqrt(t_plastic_strain(i, j) * t_plastic_strain(i, j));
    // CHKERR set_tag(th_eqv_ep, gg, set_scalar(eqv_ep));
    const double hard = hardening(tau) / (*cache).young_modulus_inv;
    CHKERR set_tag(th_hardening, gg, set_scalar(hard));
 
    ++t_flow;
    ++t_stress;
    ++t_plastic_strain;
  }
 
  MoFEMFunctionReturn(0);
}
 
template struct OpCalculatePlasticInternalForceLhs_dEP<true>;
template struct OpCalculatePlasticInternalForceLhs_dEP<false>;
 
template struct OpCalculatePlasticFlowLhs_dU<true>;
template struct OpCalculatePlasticFlowLhs_dU<false>;
 
template struct OpCalculateConstraintLhs_dU<true>;
template struct OpCalculateConstraintLhs_dU<false>;
 
} // namespace OpPlasticTools