PlasticOperators.hpp

Go to the documentation of this file.

/* 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/>. */
 
/** \file OpPlasticOperators.hpp
 
\f[
\left\{
\begin{array}{ll}
\frac{\partial \sigma_{ij}}{\partial x_j} - b_i = 0 & \forall x \in \Omega \\
\varepsilon_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} +
\frac{\partial u_j}{\partial x_i} \right)\\
\sigma_{ij} = D_{ijkl}\left(\varepsilon_{kl}-\varepsilon^p_{kl}\right) \\
\dot{\varepsilon}^p_{kl} - \dot{\tau} \left( \left. \frac{\partial f}{\partial
\sigma_{kl}} \right|_{(\sigma,\tau) } \right) = 0 \\
f(\sigma, \tau) \leq 0,\; \dot{\tau} \geq 0,\;\dot{\tau}f(\sigma, \tau)=0\\
u_i = \overline{u}_i & \forall x \in \partial\Omega_u \\
\sigma_{ij}n_j = \overline{t}_i & \forall x \in \partial\Omega_\sigma \\
\Omega_u \cup \Omega_\sigma = \Omega \\
\Omega_u \cap \Omega_\sigma = \emptyset
\end{array}
\right.
\f]
 
\f[
\left\{
\begin{array}{ll}
\left(\frac{\partial \delta u_i}{\partial x_j},\sigma_{ij}\right)_\Omega-(\delta
u_i,b_i)_\Omega -(\delta u_i,\overline{t}_i)_{\partial\Omega_\sigma}=0 & \forall
\delta u_i \in H^1(\Omega)\\ \left(\delta\varepsilon^p_{kl} ,D_{ijkl}\left(
\dot{\varepsilon}^p_{kl} - \dot{\tau} A_{kl} \right)\right) = 0
& \forall \delta\varepsilon^p_{ij} \in L^2(\Omega) \cap \mathcal{S} \\
\left(\delta\tau,c_n\dot{\tau} - \frac{1}{2}\left\{c_n \dot{\tau} +
(f(\pmb\sigma,\tau) - \sigma_y) +
\| c_n \dot{\tau} + (f(\pmb\sigma,\tau) - \sigma_y) \|\right\}\right) = 0 &
\forall \delta\tau \in L^2(\Omega) \end{array} \right. \f]
 
*/
 
namespace OpPlasticTools {
 
//! [Common data]
struct CommonData : public OpElasticTools::CommonData {
  boost::shared_ptr<VectorDouble> plasticSurfacePtr;
  boost::shared_ptr<MatrixDouble> plasticFlowPtr;
  boost::shared_ptr<VectorDouble> plasticTauPtr;
  boost::shared_ptr<VectorDouble> plasticTauDotPtr;
  boost::shared_ptr<MatrixDouble> plasticStrainPtr;
  boost::shared_ptr<MatrixDouble> plasticStrainDotPtr;
 
  boost::shared_ptr<MatrixDouble> plasticGradTauPtr;
  boost::shared_ptr<MatrixDouble> plasticGradStrainPtr;
 
  boost::shared_ptr<VectorDouble> plasticTauJumpPtr;
  boost::shared_ptr<MatrixDouble> plasticStrainJumpPtr;
  boost::shared_ptr<MatrixDouble> guidingVelocityPtr;
 
  Ddg<double, 3, 3> diffDeviator;
 
  MatrixDouble dualBaseMat;
  bool isDualBase;
 
  // data for skeleton computation
  map<int, MatrixDouble> plasticStrainSideMap;
  map<int, VectorDouble> plasticTauSideMap;
  map<int, MatrixDouble> velocityVecSideMap;
 
  MoFEMErrorCode calculateDiffDeviator();
};
 
struct OpCalculatePlasticSurfaceAndFlow : public DomainEleOp {
  OpCalculatePlasticSurfaceAndFlow(
      const std::string field_name,
      boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpPlasticStress : public DomainEleOp {
  OpPlasticStress(const std::string field_name,
                  boost::shared_ptr<CommonData> common_data_ptr,
                  boost::shared_ptr<MatrixDouble> mat_d_ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  boost::shared_ptr<MatrixDouble> matDPtr;
};
 
struct OpCalculatePlasticFlowRhs : public DomainEleOpAssembly {
  OpCalculatePlasticFlowRhs(const std::string field_name,
                            boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpCalculateConstraintRhs : public DomainEleOpAssembly {
  OpCalculateConstraintRhs(const std::string field_name,
                           boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
template <bool LOGSTRAIN>
struct OpCalculatePlasticInternalForceLhs_dEP : public DomainEleOpAssembly {
  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);
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  boost::shared_ptr<MatrixDouble> matDPtr;
};
 
template <bool LOGSTRAIN>
struct OpCalculatePlasticFlowLhs_dU : public DomainEleOpAssembly {
  OpCalculatePlasticFlowLhs_dU(const std::string row_field_name,
                               const std::string col_field_name,
                               boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpCalculatePlasticFlowLhs_dEP : public DomainEleOpAssembly {
  OpCalculatePlasticFlowLhs_dEP(const std::string row_field_name,
                                const std::string col_field_name,
                                boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpCalculatePlasticFlowLhs_dTAU : public DomainEleOpAssembly {
  OpCalculatePlasticFlowLhs_dTAU(const std::string row_field_name,
                                 const std::string col_field_name,
                                 boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpCalculatePlasticInternalForceLhs_dTAU : public DomainEleOpAssembly {
  OpCalculatePlasticInternalForceLhs_dTAU(const std::string row_field_name,
                                 const std::string col_field_name,
                                 boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
template <bool LOGSTRAIN>
struct OpCalculateConstraintLhs_dU : public DomainEleOpAssembly {
  OpCalculateConstraintLhs_dU(const std::string row_field_name,
                              const std::string col_field_name,
                              boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpCalculateConstraintLhs_dEP : public DomainEleOpAssembly {
  OpCalculateConstraintLhs_dEP(const std::string row_field_name,
                               const std::string col_field_name,
                               boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpCalculateConstraintLhs_dTAU : public DomainEleOpAssembly {
  OpCalculateConstraintLhs_dTAU(const std::string row_field_name,
                                const std::string col_field_name,
                                boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpPostProcPlastic : public DomainEleOp {
  OpPostProcPlastic(const std::string field_name,
                    moab::Interface &post_proc_mesh,
                    std::vector<EntityHandle> &map_gauss_pts,
                    boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  moab::Interface &postProcMesh;
  std::vector<EntityHandle> &mapGaussPts;
  boost::shared_ptr<CommonData> commonDataPtr;
};
//! [Operators definitions]
 
//! [Lambda functions]
inline auto diff_tensor() {
  Ddg<double, 3, 3> t_diff;
  constexpr auto t_kd = Kronecker_Delta_symmetric<int>();
  t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
  return t_diff;
};
 
inline auto symm_L_tensor() {
  Dg<double, 3, 6> L;
  Index<'L', 6> Li;
  // dg(T 000, T 001, T 002, T 010, T 011, T 012, T 020, T 021, T 022,
  //    T 110, T 111, T 112, T 120, T 121, T 122, T 220, T 221, T 222);
 
  // Tensor2_symmetric_constructor(T data[(3 * 4) / 2], T d00, T d01, T d02, T
  // d11,
  //                               T d12, T d22);
  // Tensor2<double, 3, 3> test_tensor{1., 2., 3., 2., 5., 6., 3., 6., 9.};
  // Tensor1<double, 6> test_vec;
 
  L(i, j, Li) = 0;
  L(0, 0, 0) = 1;
  L(0, 1, 1) = 1;
  L(0, 2, 2) = 1;
  L(1, 0, 1) = 1;
  L(1, 1, 3) = 1;
  L(1, 2, 4) = 1;
  L(2, 0, 2) = 1;
  L(2, 1, 4) = 1;
  L(2, 2, 5) = 1;
 
  // test_vec(Li) = test_tensor(i, j) * L(i, j, Li);
 
  //  (00,  01,  02,  11,  12,  22)
  return L;
}
 
inline auto diff_symmetrize() { return (*cache).Is; };
 
template <typename T1>
inline auto deviator(Tensor2_symmetric<T1, 3> &t_stress) {
  Tensor2_symmetric<double, 3> t_dev;
  t_dev(I, J) = 0;
  constexpr double third = boost::math::constants::third<double>();
  const double trace = t_stress(i, i);
 
  t_dev(i, j) = t_stress(i, j);
  t_dev(0, 0) -= trace * third;
  t_dev(1, 1) -= trace * third;
  t_dev(2, 2) -= trace * third;
  return t_dev;
};
 
inline auto diff_deviator(Ddg<double, 3, 3> &&t_diff_stress) {
  Ddg<double, 3, 3> t_diff_deviator;
  t_diff_deviator(I, J, K, 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 = kk; ll != 3; ++ll)
          t_diff_deviator(ii, jj, kk, ll) = t_diff_stress(ii, jj, kk, ll);
 
  constexpr double third = boost::math::constants::third<double>();
 
  t_diff_deviator(0, 0, 0, 0) = 1. - third;
  t_diff_deviator(0, 0, 1, 1) = -third;
  t_diff_deviator(0, 0, 2, 2) = -third;
 
  t_diff_deviator(1, 1, 0, 0) = -third;
  t_diff_deviator(1, 1, 1, 1) = 1. - third;
  t_diff_deviator(1, 1, 2, 2) = -third;
 
  t_diff_deviator(2, 2, 0, 0) = -third;
  t_diff_deviator(2, 2, 1, 1) = -third;
  t_diff_deviator(2, 2, 2, 2) = 1. - third;
 
  return t_diff_deviator;
};
 
template <typename T>
inline auto get_effective_stress(Tensor2_symmetric<T, 3> &t_stress,
                                 Tensor2_symmetric<T, 3> &t_plastic_strain) {
 
  Tensor2_symmetric<double, 3> stress_tmp;
 
  stress_tmp(i, j) =
      t_stress(i, j) - 1.5 * (*cache).C1_k * t_plastic_strain(i, j);
 
  return stress_tmp;
};
 
constexpr double sqrt2by3 = 0.816496580927726;
inline auto hardening(long double tau) {
  return (*cache).H * tau + (*cache).Q_inf * (1. - exp(-(*cache).b_iso * tau)) +
         (*cache).sigmaY;
}
 
inline auto hardening_dtau(long double tau) {
  return (*cache).H +
         (*cache).Q_inf * (*cache).b_iso * exp(-(*cache).b_iso * tau);
}
 
/**
 *
 
\f[
\begin{split}
f&=\sqrt{s_{ij}s_{ij}}\\
A_{ij}&=\frac{\partial f}{\partial \sigma_{ij}}=
\frac{1}{f} s_{kl} \frac{\partial s_{kl}}{\partial \sigma_{ij}}\\
\frac{\partial A_{ij}}{\partial \sigma_{kl}}&= \frac{\partial^2 f}{\partial
\sigma_{ij}\partial\sigma_{mn}}= \frac{1}{f} \left( \frac{\partial
s_{kl}}{\partial \sigma_{mn}}\frac{\partial s_{kl}}{\partial \sigma_{ij}}
-A_{mn}A_{ij}
\right)\\
\frac{\partial f}{\partial \varepsilon_{ij}}&=A_{mn}D_{mnij}
\\
\frac{\partial A_{ij}}{\partial \varepsilon_{kl}}&=
\frac{\partial A_{ij}}{\partial \sigma_{mn}} \frac{\partial
\sigma_{mn}}{\partial \varepsilon_{kl}}= \frac{\partial A_{ij}}{\partial
\sigma_{mn}} D_{mnkl}
\end{split}
\f]
 
 */
inline auto plastic_surface(Tensor2_symmetric<double, 3> &&t_stress_deviator) {
  return std::sqrt(1.5 * t_stress_deviator(I, J) * t_stress_deviator(I, J)) +
         std::numeric_limits<double>::epsilon();
};
 
inline auto plastic_flow(long double f,
                         Tensor2_symmetric<double, 3> &&t_dev_stress,
                         Ddg<double, 3, 3> &t_diff_deviator) {
  Tensor2_symmetric<double, 3> t_diff_f;
  t_diff_f(k, l) =
      ((1.5) * (t_dev_stress(I, J)) * t_diff_deviator(I, J, k, l)) / f;
 
  return t_diff_f;
};
 
// FIXME: get rid of second projection for efficiency
// possibly we can also remove D : from the flow rule for the same reason
 
template <typename T>
inline auto diff_plastic_flow_dstress(long double f,
                                      Tensor2_symmetric<T, 3> &t_flow,
                                      Ddg<double, 3, 3> &t_diff_deviator) {
  Ddg<double, 3, 3> t_diff_flow;
  t_diff_flow(i, j, k, l) =
      ((1.5) * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l) -
                2. / 3. * t_flow(i, j) * t_flow(k, l))) /
      f;
 
  return t_diff_flow;
};
 
template <typename T>
inline auto
diff_plastic_flow_dstrain(Ddg<T, 3, 3> &t_D,
                          Ddg<double, 3, 3> &&t_diff_plastic_flow_dstress) {
  Ddg<double, 3, 3> t_diff_flow;
  t_diff_flow(i, j, k, l) =
      t_diff_plastic_flow_dstress(i, j, m, n) * t_D(m, n, k, l);
 
  return t_diff_flow;
};
 
template <typename T>
inline auto diff_plastic_flow_kin_hard_dstrain(
    Ddg<T, 3, 3> &t_D, Ddg<double, 3, 3> &&t_diff_plastic_flow_dstress) {
  Ddg<double, 3, 3> t_diff_flow_kin_hard;
  t_diff_flow_kin_hard(i, j, k, l) =
      t_diff_plastic_flow_dstress(i, j, m, n) * t_D(m, n, k, l) +
      1.5 * (*cache).C1_k * t_diff_plastic_flow_dstress(i, j, k, l);
 
  return t_diff_flow_kin_hard;
};
 
inline double constrain_abs(long double x) { return std::abs(x); };
 
inline double constraint_sign(long double x) {
  if (x > 0)
    return 1;
  else if (x < 0)
    return -1;
  else
    return 0;
};
 
inline double w(long double tau_dot, long double f, long double hardening) {
  auto &cn = (*cache).cn_pl;
  auto &sigmaY = (*cache).sigmaY;
  return (f - hardening) / sigmaY + cn * tau_dot;
};
 
inline double constraint(long double tau_dot, long double f,
                         long double hardening) {
  auto &cn = (*cache).cn_pl;
  auto &visH = (*cache).visH;
  auto &sigmaY = (*cache).sigmaY;
  return visH * tau_dot +
         (sigmaY / 2) * ((cn * tau_dot - (f - hardening) / sigmaY) -
                         constrain_abs(w(tau_dot, f, hardening)));
};
 
inline double diff_constrain_dtau_dot(long double tau_dot, long double f,
                                      long double hardening) {
  auto &cn = (*cache).cn_pl;
  auto &visH = (*cache).visH;
  auto &sigmaY = (*cache).sigmaY;
  return visH +
         (sigmaY / 2) * (cn - cn * constraint_sign(w(tau_dot, f, hardening)));
};
 
inline auto diff_constrain_df(long double tau_dot, long double f,
                              long double hardening) {
  return (-1 - constraint_sign(w(tau_dot, f, hardening))) / 2;
};
 
inline auto diff_constrain_dsigma_y(long double tau_dot, long double f,
                                    long double hardening) {
  return (1 + constraint_sign(w(tau_dot, f, hardening))) / 2;
}
 
template <typename T>
inline auto diff_constrain_dstress(long double &&diff_constrain_df,
                                   Tensor2_symmetric<T, 3> &t_plastic_flow) {
  Tensor2_symmetric<double, 3> t_diff_constrain_dstress;
  t_diff_constrain_dstress(i, j) = diff_constrain_df * t_plastic_flow(i, j);
  return t_diff_constrain_dstress;
};
 
template <typename T1, typename T2>
inline auto
diff_constrain_dstrain(Ddg<T1, 3, 3> &t_D,
                       Tensor2_symmetric<T2, 3> &&t_diff_constrain_dstress) {
  Tensor2_symmetric<double, 3> t_diff_constrain_dstrain;
  t_diff_constrain_dstrain(k, l) =
      t_diff_constrain_dstress(i, j) * t_D(i, j, k, l);
  return t_diff_constrain_dstrain;
};
 
template <typename T1, typename T2>
inline auto diff_constrain_kin_hard_dstrain(
    Ddg<T1, 3, 3> &t_D, Tensor2_symmetric<T2, 3> &&t_diff_constrain_dstress) {
  Tensor2_symmetric<double, 3> t_diff_constrain_kin_hard_dstrain;
  t_diff_constrain_kin_hard_dstrain(k, l) =
      t_diff_constrain_dstress(i, j) * t_D(i, j, k, l) +
      1.5 * (*cache).C1_k * t_diff_constrain_dstress(k, l);
  return t_diff_constrain_kin_hard_dstrain;
};
 
struct OpGetGaussPtsPlasticState : public DomainEleOp {
 
  OpGetGaussPtsPlasticState(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);
  }
 
  //! [Calculate stress]
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
    // const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
    auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
    auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
    auto &data_vec = commonDataPtr->stateArrayPlast;
    for (auto &f : *(commonDataPtr->plasticSurfacePtr)) {
 
      if (w(t_tau_dot, f, hardening(t_tau)) > 0)
        data_vec[0] += 1;
      data_vec[1] += 1;
 
      ++t_tau_dot;
      ++t_tau;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct BubbleTauPolynomialBase : public BaseFunction {
 
  BubbleTauPolynomialBase() = default;
  ~BubbleTauPolynomialBase() = default;
 
  /**
   * @brief Return interface to this class when one ask for for tetrahedron,
   * otherisw return interface class for generic class.
   *
   * @param iface interface class
   * @return MoFEMErrorCode
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const {
    MoFEMFunctionBegin;
    *iface = const_cast<BubbleTauPolynomialBase *>(this);
    MoFEMFunctionReturn(0);
  }
 
  /**
   * @brief Calculate base functions at intergeneration points
   *
   * @param pts
   * @param ctx_ptr
   * @return MoFEMErrorCode
   */
  MoFEMErrorCode getValue(MatrixDouble &pts,
                          boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
    MoFEMFunctionBeginHot;
 
    cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
 
    int nb_gauss_pts = pts.size2();
    if (!nb_gauss_pts) {
      MoFEMFunctionReturnHot(0);
    }
 
    if (pts.size1() < 3) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Wrong dimension of pts, should be at least 3 rows with "
              "coordinates");
    }
 
    switch (cTx->sPace) {
    case L2:
      // CHKERR getValueL2ForTauBubble(pts);
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
    }
 
    MoFEMFunctionReturnHot(0);
  }
 
private:
  EntPolynomialBaseCtx *cTx;
 
  // MatrixDouble shapeFun;
 
  // MoFEMErrorCode getValueL2ForTauBubble(MatrixDouble &pts) {
  //   MoFEMFunctionBegin;
 
  //   const FieldApproximationBase base = cTx->bAse;
  //   // This should be used only in case USER_BASE is selected
  //   if (cTx->bAse != USER_BASE) {
  //     SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
  //             "Wrong base, should be USER_BASE");
  //   }
 
  //   EntitiesFieldData &data = cTx->dAta;
  //   int nb_gauss_pts = pts.size2();
 
  //   data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4, false);
  //   data.dataOnEntities[MBTET][0].getN(base).resize(
  //       nb_gauss_pts,
  //       NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getDataOrder()), false);
  //   data.dataOnEntities[MBTET][0].getDiffN(base).resize(
  //       nb_gauss_pts,
  //       3 * NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getDataOrder()),
  //       false);
 
  //   CHKERR L2_Ainsworth_ShapeFunctions_MBTET(
  //       data.dataOnEntities[MBTET][0].getDataOrder(),
  //       &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
  //       &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
  //       &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
  //       &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
  //       nb_gauss_pts, cTx->basePolynomialsType0);
 
  //   auto tetN = data.dataOnEntities[MBTET][0].getN(base);
  //   data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 2);
  //   for (int gg = 0; gg != nb_gauss_pts; ++gg) {
  //     data.dataOnEntities[MBTET][0].getN(base)(gg, 0) = 1;
  //     auto prod_n0123 = tetN(gg, 0) * tetN(gg, 1) * tetN(gg, 2) * tetN(gg, 3);
  //     data.dataOnEntities[MBTET][0].getN(base)(gg, 1) = prod_n0123;
  //   }
 
  //   MoFEMFunctionReturn(0);
  // }
};
 
//! [Postprocessing]
}; // namespace OpPlasticTools