EshelbianOperators.cpp

Implementation of operators

/**
 * \file EshelbianOperators.cpp
 * \example EshelbianOperators.cpp
 *
 * \brief Implementation of operators
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
 
#include <boost/math/constants/constants.hpp>
 
#include <EshelbianAux.hpp>
 
#include <lapack_wrap.h>
 
#include <MatrixFunction.hpp>
 
namespace EshelbianPlasticity {
 
template <typename T>
auto get_rotation_form_vector(FTensor::Tensor1<T, 3> &t_omega,
                              RotSelector rotSelector = LARGE_ROT) {
 
  using D = typename TensorTypeExtractor<T>::Type;
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Tensor2<D, 3, 3> t_R;
 
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
  t_R(i, j) = t_kd(i, j);
 
  FTensor::Tensor2<D, 3, 3> t_Omega;
  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
 
  if (rotSelector == SMALL_ROT) {
    t_R(i, j) += t_Omega(i, j);
    return t_R;
  }
 
  const auto angle =
      sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
  const auto a = sin(angle) / angle;
  t_R(i, j) += a * t_Omega(i, j);
  if (rotSelector == MODERATE_ROT)
    return t_R;
 
  const auto ss_2 = sin(angle / 2.) / angle;
  const auto b = 2. * ss_2 * ss_2;
  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
 
  return t_R;
}
 
template <typename T>
auto get_diff_rotation_form_vector(FTensor::Tensor1<T, 3> &t_omega,
                                   RotSelector rotSelector = LARGE_ROT) {
 
  using D = typename TensorTypeExtractor<T>::Type;
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  FTensor::Tensor3<D, 3, 3, 3> t_diff_R;
 
  const auto angle =
      sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
  if (rotSelector == SMALL_ROT) {
    t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, j, k);
    return t_diff_R;
  }
 
  const auto ss = sin(angle);
  const auto a = ss / angle;
  t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
 
  FTensor::Tensor2<D, 3, 3> t_Omega;
  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
 
  const auto angle2 = angle * angle;
  const auto cc = cos(angle);
  const auto diff_a = (angle * cc - ss) / angle2;
  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
  if (rotSelector == MODERATE_ROT)
    return t_diff_R;
 
  const auto ss_2 = sin(angle / 2.);
  const auto cc_2 = cos(angle / 2.);
  const auto b = 2. * ss_2 * ss_2 / angle2;
  t_diff_R(i, j, k) +=
      b * (t_Omega(i, l) * FTensor::levi_civita<double>(l, j, k) +
           FTensor::levi_civita<double>(i, l, k) * t_Omega(l, j));
  const auto diff_b =
      (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
  t_diff_R(i, j, k) +=
      diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
 
  return t_diff_R;
}
 
template <typename T>
auto get_diff2_rotation_form_vector(FTensor::Tensor1<T, 3> &t_omega,
                                    RotSelector rotSelector = LARGE_ROT) {
 
  using D = typename TensorTypeExtractor<T>::Type;
 
  FTensor::Tensor4<D, 3, 3, 3, 3> t_diff2_R;
  FTensor::Index<'i', 3> i;
 
  constexpr double eps = 1e-10;
  for (int l = 0; l != 3; ++l) {
    FTensor::Tensor1<std::complex<double>, 3> t_omega_c;
    t_omega_c(i) = t_omega(i);
    t_omega_c(l) += std::complex<double>(0, eps);
    auto t_diff_R_c = get_diff_rotation_form_vector(t_omega_c, rotSelector);
    for (int i = 0; i != 3; ++i) {
      for (int j = 0; j != 3; ++j) {
        for (int k = 0; k != 3; ++k) {
          t_diff2_R(i, j, k, l) = t_diff_R_c(i, j, k).imag() / eps;
        }
      }
    }
  }
 
  return t_diff2_R;
}
 
struct isEq {
  static inline auto check(const double &a, const double &b) {
    constexpr double eps = std::numeric_limits<float>::epsilon();
    return std::abs(a - b) / absMax < eps;
  }
  static double absMax;
};
 
double isEq::absMax = 1;
 
inline auto is_eq(const double &a, const double &b) {
  return isEq::check(a, b);
};
 
template <int DIM> inline auto get_uniq_nb(double *ptr) {
  std::array<double, DIM> tmp;
  std::copy(ptr, &ptr[DIM], tmp.begin());
  std::sort(tmp.begin(), tmp.end());
  isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
  constexpr double eps = std::numeric_limits<float>::epsilon();
  isEq::absMax = std::max(isEq::absMax, static_cast<double>(eps));
  return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
}
 
template <int DIM>
inline auto sort_eigen_vals(FTensor::Tensor1<double, DIM> &eig,
                            FTensor::Tensor2<double, DIM, DIM> &eigen_vec) {
 
  isEq::absMax =
      std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
 
  int i = 0, j = 1, k = 2;
 
  if (is_eq(eig(0), eig(1))) {
    i = 0;
    j = 2;
    k = 1;
  } else if (is_eq(eig(0), eig(2))) {
    i = 0;
    j = 1;
    k = 2;
  } else if (is_eq(eig(1), eig(2))) {
    i = 1;
    j = 0;
    k = 2;
  }
 
  FTensor::Tensor2<double, 3, 3> eigen_vec_c{
      eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
 
      eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
 
      eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
 
  FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
 
  {
    FTensor::Index<'i', DIM> i;
    FTensor::Index<'j', DIM> j;
    eig(i) = eig_c(i);
    eigen_vec(i, j) = eigen_vec_c(i, j);
  }
}
 
 
 
MoFEMErrorCode OpCalculateEshelbyStress::doWork(int side, EntityType type,
                                                EntData &data) {
  MoFEMFunctionBegin;
  FTENSOR_INDEX(SPACE_DIM, i);  
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, m);
 
  int nb_integration_pts = getGaussPts().size2();
 
  auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
  auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->hAtPts);
  auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
 
  dataAtPts->SigmaAtPts.resize(SPACE_DIM * SPACE_DIM, nb_integration_pts,
                                 false);
  auto t_eshelby_stress =
      getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->SigmaAtPts);
 
  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
 
  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
    t_eshelby_stress(i, j) = t_energy * t_kd(i, j) - t_F(m, i) * t_P(m, j);
    ++t_energy;
    ++t_P;
    ++t_F;
    ++t_eshelby_stress;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork(int side,
                                                             EntityType type,
                                                             EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = symm_L_tensor();
 
  int nb_integration_pts = getGaussPts().size2();
  i_FTIndex<SPACE_DIM> i;
  j_FTIndex<SPACE_DIM> j;
  k_FTIndex<SPACE_DIM> k;
  l_FTIndex<SPACE_DIM> l;
  m_FTIndex<SPACE_DIM> m;
  n_FTIndex<SPACE_DIM> n;
 
  dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
  dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts, false);
  dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts, false);
  dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts, false);
  dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts, false);
  dataAtPts->leviKirchhoffOmegaAtPts.resize(9, nb_integration_pts, false);
 
  dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts, false);
  dataAtPts->adjointPdUAtPts.resize(size_symm, nb_integration_pts, false);
  dataAtPts->adjointPdUdPAtPts.resize(9 * size_symm, nb_integration_pts, false);
  dataAtPts->adjointPdUdOmegaAtPts.resize(3 * size_symm, nb_integration_pts,
                                          false);
  dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
  dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts, false);
  dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts, false);
  dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts, false);
  dataAtPts->eigenVals.resize(3, nb_integration_pts, false);
  dataAtPts->eigenVecs.resize(9, nb_integration_pts, false);
  dataAtPts->nbUniq.resize(nb_integration_pts, false);
 
  dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts, false);
  dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts, false);
 
  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
  auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
  auto t_h_dlog_u =
      getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
  auto t_levi_kirchoff_dP =
      getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
  auto t_levi_kirchoff_domega =
      getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
  auto t_approx_P_adjont_dstretch =
      getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
  auto t_approx_P_adjont_log_du =
      getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
  auto t_approx_P_adjont_log_du_dP =
      getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
  auto t_approx_P_adjont_log_du_domega =
      getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_diff_R = getFTensor3FromMat<3, 3, 3>(dataAtPts->diffRotMatAtPts);
  auto t_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
  auto t_diff_u =
      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
 
  auto &nbUniq = dataAtPts->nbUniq;
  auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
  auto t_log_stretch_total =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
  auto t_log_u2_h1 =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
  auto t_kd = FTensor::Kronecker_Delta<int>();
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
    FTensor::Tensor1<double, 3> eig;
    FTensor::Tensor2<double, 3, 3> eigen_vec;
    FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
    FTensor::Tensor2<double, 3, 3> t_h1;
    FTensor::Tensor2<double, 3, 3> t_Ru;
    FTensor::Tensor2<double, 3, 3> t_approx_P_intermidiate;
 
    switch (EshelbianCore::gradApproximator) {
    case NO_H1_CONFIGURATION:
      t_h1(i, j) = t_kd(i, j);
      break;
    case LARGE_ROT:
    case MODERATE_ROT:
      t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
      break;
    case SMALL_ROT:
      t_h1(i, j) = t_kd(i, j);
      break;
    };
 
    auto calculate_rotation = [&]() {
      auto t0_diff =
          get_diff_rotation_form_vector(t_omega, EshelbianCore::rotSelector);
      auto t0 = get_rotation_form_vector(t_omega, EshelbianCore::rotSelector);
      t_diff_R(i, j, k) = t0_diff(i, j, k);
      t_R(i, j) = t0(i, j);
    };
 
    auto calculate_stretch = [&]() {
      if (EshelbianCore::stretchSelector == StretchSelector::LINEAR) {
        constexpr auto t_kd_sym = FTensor::Kronecker_Delta_symmetric<int>();
        t_u(i, j) = t_log_u(i, j) + t_kd_sym(i, j);
        t_diff_u(i, j, k, l) = (t_kd_sym(i, k) ^ t_kd_sym(j, l)) / 4.;
        t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
      } else {
        eigen_vec(i, j) = t_log_u(i, j);
        CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
        // rare case when two eigen values are equal
        nbUniq[gg] = get_uniq_nb<3>(&eig(0));
        if (nbUniq[gg] < 3) {
          sort_eigen_vals<3>(eig, eigen_vec);
        }
        t_eigen_vals(i) = eig(i);
        t_eigen_vecs(i, j) = eigen_vec(i, j);
        auto t_u_tmp =
            EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, EshelbianCore::f);
        t_u(i, j) = t_u_tmp(i, j);
        auto t_diff_u_tmp = EigenMatrix::getDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            nbUniq[gg]);
        t_diff_u(i, j, k, l) = t_diff_u_tmp(i, j, k, l);
        t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
      }
    };
 
    calculate_rotation();
 
    if (!EshelbianCore::noStretch) {
 
      calculate_stretch();
 
      switch (EshelbianCore::gradApproximator) {
      case NO_H1_CONFIGURATION:
      case LARGE_ROT:
 
        // rotation can be large, moderate or small
        t_Ru(i, m) = t_R(i, l) * t_u(l, m);
        t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
        t_h_domega(i, j, k) = (t_diff_R(i, l, k) * t_u(l, m)) * t_h1(m, j);
        t_h_dlog_u(i, j, L) = (t_R(i, l) * t_Ldiff_u(l, m, L)) * t_h1(m, j);
 
        // conservation of angular momentum at current configuration
        t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
        t_approx_P_adjont_dstretch(l, m) =
            t_approx_P_intermidiate(i, m) * t_R(i, l);
 
        t_levi_kirchoff(k) =
            levi_civita(l, m, k) * (t_approx_P_adjont_dstretch(l, m));
        t_levi_kirchoff_dP(i, j, k) =
            (levi_civita(l, m, k) * t_R(i, l)) * t_h1(m, j);
        t_levi_kirchoff_domega(k, n) =
            levi_civita(l, m, k) *
            (t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
 
        t_approx_P_adjont_log_du(L) =
            t_Ldiff_u(l, m, L) * t_approx_P_adjont_dstretch(l, m);
        t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
        t_approx_P_adjont_log_du_domega(n, L) =
            t_Ldiff_u(l, m, L) *
            (t_approx_P_intermidiate(i, m) * t_diff_R(i, l, n));
 
        break;
      case MODERATE_ROT:
 
      {
 
        // assume small current rotation
        FTensor::Tensor2<double, 3, 3> t_Omega;
        t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
        t_Ru(i, m) = t_Omega(i, m) + t_u(i, m);
        t_h(i, j) = t_Ru(i, m) * t_h1(m, j);
        t_h_domega(i, j, k) =
            FTensor::levi_civita<double>(i, m, k) * t_h1(m, j);
        t_h_dlog_u(i, j, L) = t_Ldiff_u(i, m, L) * t_h1(m, j);
 
        // conservation of angular momentum at intermediate configuration
        t_approx_P_intermidiate(i, m) = t_approx_P(i, j) * t_h1(m, j);
        t_approx_P_adjont_dstretch(i, m) = t_approx_P_intermidiate(i, m);
 
        t_levi_kirchoff(k) =
            levi_civita(i, m, k) * (t_approx_P_adjont_dstretch(i, m));
        t_levi_kirchoff_dP(i, j, k) = levi_civita(i, m, k) * t_h1(m, j);
        t_levi_kirchoff_domega(k, n) = 0;
 
        t_approx_P_adjont_log_du(L) =
            t_Ldiff_u(i, m, L) * t_approx_P_adjont_dstretch(i, m);
        t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
        t_approx_P_adjont_log_du_domega(n, L) = 0;
 
      }
 
      break;
 
      case SMALL_ROT:
 
      {
 
        // assume small  current rotation
        FTensor::Tensor2<double, 3, 3> t_Omega;
        t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
        t_h(i, j) = t_Omega(i, j) + t_u(i, j);
        t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, j, k);
        t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
 
        // conservation of angular momentum at reference condition
        t_levi_kirchoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
        t_levi_kirchoff_dP(i, j, k) = levi_civita(i, j, k);
        t_levi_kirchoff_domega(k, l) = 0;
 
        t_approx_P_adjont_dstretch(i, j) = t_approx_P(i, j);
        t_approx_P_adjont_log_du(L) =
            t_Ldiff_u(i, j, L) * t_approx_P_adjont_dstretch(i, j);
        t_approx_P_adjont_log_du_dP(i, j, L) = t_h_dlog_u(i, j, L);
        t_approx_P_adjont_log_du_domega(m, L) = 0;
 
      }
 
      break;
      }
 
      FTensor::Tensor2<double, 3, 3> t_C_h1;
      t_C_h1(i, j) = t_h1(k, i) * t_h1(k, j);
      eigen_vec(i, j) = t_C_h1(i, j);
      CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
      switch (EshelbianCore::stretchSelector) {
      case StretchSelector::LINEAR:
        break;
      case StretchSelector::LOG:
        // for (int ii = 0; ii != 3; ++ii) {
        //   eig(ii) = std::max(
        //       eig(ii), 2 * EshelbianCore::inv_f(EshelbianCore::f(
        //                        std::numeric_limits<double>::min_exponent)));
        // }
        break;
      }
 
      switch (EshelbianCore::gradApproximator) {
      case NO_H1_CONFIGURATION:
        t_log_u2_h1(i, j) = 0;
        t_log_stretch_total(i, j) = t_log_u(i, j);
        break;
      case LARGE_ROT:
      case MODERATE_ROT: {
        auto t_log_u2_h1_tmp =
            EigenMatrix::getMat(eig, eigen_vec, EshelbianCore::inv_f);
        t_log_u2_h1(i, j) = t_log_u2_h1_tmp(i, j);
        t_log_stretch_total(i, j) = t_log_u2_h1_tmp(i, j) / 2 + t_log_u(i, j);
      } break;
      case SMALL_ROT:
        t_log_u2_h1(i, j) = 0;
        t_log_stretch_total(i, j) = t_log_u(i, j);
        break;
      };
 
    } else {
 
      auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
 
      FTensor::Tensor2<double, 3, 3> t_Omega;
      t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
      t_u(i, j) = t_kd(i, j) + t_log_u(i, j);
      t_h(i, j) = t_Omega(i, j) + t_u(i, j);
      t_h_domega(i, j, k) = FTensor::levi_civita<double>(i, j, k);
 
      // conservation of angular momentum at reference condition
      t_levi_kirchoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
      t_levi_kirchoff_dP(i, j, k) = levi_civita(i, j, k);
      t_levi_kirchoff_domega(k, l) = 0;
 
      t_log_stretch_total(i, j) = t_log_u(i, j);
    }
 
    ++t_h;
    ++t_h_domega;
    ++t_h_dlog_u;
    ++t_levi_kirchoff;
    ++t_levi_kirchoff_dP;
    ++t_levi_kirchoff_domega;
    ++t_approx_P_adjont_dstretch;
    ++t_approx_P_adjont_log_du;
    ++t_approx_P_adjont_log_du_dP;
    ++t_approx_P_adjont_log_du_domega;
    ++t_approx_P;
    ++t_R;
    ++t_diff_R;
    ++t_log_u;
    ++t_u;
    ++t_diff_u;
    ++t_eigen_vals;
    ++t_eigen_vecs;
    ++t_omega;
    ++t_grad_h1;
    ++t_log_u2_h1;
    ++t_log_stretch_total;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialEquilibrium::integrate(EntData &data) {
  MoFEMFunctionBegin;
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
  auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
  if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
      dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
    dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
    dataAtPts->wL2DotDotAtPts.clear();
  }
  auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
 
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
 
  FTensor::Index<'i', 3> i;
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
                                                              &v[2]);
  };
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor1(nF);
    int bb = 0;
    for (; bb != nb_dofs / 3; ++bb) {
      t_nf(i) += a * t_row_base_fun * t_div_P(i);
      t_nf(i) -= a * t_row_base_fun * alphaW * t_s_dot_w(i);
      t_nf(i) -= a * t_row_base_fun * alphaRho * t_s_dot_dot_w(i);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
    ++t_w;
    ++t_div_P;
    ++t_s_dot_w;
    ++t_s_dot_dot_w;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialRotation::integrate(EntData &data) {
  MoFEMFunctionBegin;
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = getGaussPts().size2();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
                                                              &v[2]);
  };
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor1(nF);
    int bb = 0;
    for (; bb != nb_dofs / 3; ++bb) {
      t_nf(k) += (a * t_row_base_fun) * t_levi_kirchoff(k);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb) {
      ++t_row_base_fun;
    }
    ++t_w;
    ++t_levi_kirchoff;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialConsistencyP::integrate(EntData &data) {
  MoFEMFunctionBegin;
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
 
  int nb_base_functions = data.getN().size2() / 3;
  auto t_row_base_fun = data.getFTensor1N<3>();
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'m', 3> m;
 
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
                                                              &v[2]);
  };
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor1(nF);
 
    constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
    FTensor::Tensor2<double, 3, 3> t_residuum;
    t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
 
    int bb = 0;
    for (; bb != nb_dofs / 3; ++bb) {
      t_nf(i) += a * t_row_base_fun(j) * t_residuum(i, j);
      ++t_nf;
      ++t_row_base_fun;
    }
 
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_h;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialConsistencyBubble::integrate(EntData &data) {
  MoFEMFunctionBegin;
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
 
  int nb_base_functions = data.getN().size2() / 9;
  auto t_row_base_fun = data.getFTensor2N<3, 3>();
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'m', 3> m;
 
  auto get_ftensor0 = [](auto &v) {
    return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(&v[0]);
  };
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor0(nF);
 
    constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
    FTensor::Tensor2<double, 3, 3> t_residuum;
    t_residuum(i, j) = t_h(i, j) - t_kd(i, j);
 
    int bb = 0;
    for (; bb != nb_dofs; ++bb) {
      t_nf += a * t_row_base_fun(i, j) * t_residuum(i, j);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb) {
      ++t_row_base_fun;
    }
    ++t_w;
    ++t_h;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialConsistencyDivTerm::integrate(EntData &data) {
  MoFEMFunctionBegin;
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
  int nb_base_functions = data.getN().size2() / 3;
  auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
  FTensor::Index<'i', 3> i;
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
                                                              &v[2]);
  };
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor1(nF);
    int bb = 0;
    for (; bb != nb_dofs / 3; ++bb) {
      double div_row_base = t_row_diff_base_fun(i, i);
      t_nf(i) += a * div_row_base * t_w_l2(i);
      ++t_nf;
      ++t_row_diff_base_fun;
    }
    for (; bb != nb_base_functions; ++bb) {
      ++t_row_diff_base_fun;
    }
    ++t_w;
    ++t_w_l2;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <AssemblyType A>
MoFEMErrorCode OpDispBcImpl<A, GAUSS>::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcDispPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
      int nb_dofs = data.getIndices().size();
      int nb_integration_pts = OP::getGaussPts().size2();
      auto t_normal = OP::getFTensor1NormalsAtGaussPts();
      auto t_w = OP::getFTensor0IntegrationWeight();
      int nb_base_functions = data.getN().size2() / 3;
      auto t_row_base_fun = data.getFTensor1N<3>();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
 
      double scale = 1;
      for (auto &sm : scalingMethodsVec) {
        scale *= sm->getScale(OP::getFEMethod()->ts_t);
      }
 
      // get bc data
      FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
      t_bc_disp(i) *= scale;
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
        int bb = 0;
        for (; bb != nb_dofs / SPACE_DIM; ++bb) {
          t_nf(i) -=
              t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_disp(i) * 0.5;
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
 
        ++t_w;
        ++t_normal;
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpDispBc::iNtegrate(EntData &data) {
  return OP::iNtegrate(data);
}
 
template <AssemblyType A>
MoFEMErrorCode OpRotationBcImpl<A, GAUSS>::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(3, i);
  FTENSOR_INDEX(3, j);
  FTENSOR_INDEX(3, k);
 
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // interate over all boundary data
  for (auto &bc : (*bcRotPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
      int nb_dofs = data.getIndices().size();
      int nb_integration_pts = OP::getGaussPts().size2();
      auto t_normal = OP::getFTensor1NormalsAtGaussPts();
      auto t_w = OP::getFTensor0IntegrationWeight();
 
      int nb_base_functions = data.getN().size2() / 3;
      auto t_row_base_fun = data.getFTensor1N<3>();
 
      auto get_ftensor1 = [](auto &v) {
        return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
                                                                  &v[2]);
      };
 
      // get bc data
      FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
 
      double theta = bc.theta;
      theta *= OP::getFEMethod()->ts_t;
 
      FTensor::Tensor1<double, 3> t_omega;
      const double a = sqrt(t_normal(i) * t_normal(i));
      t_omega(i) = t_normal(i) * (theta / a);
      auto t_R = get_rotation_form_vector(t_omega, LARGE_ROT);
      auto t_coords = OP::getFTensor1CoordsAtGaussPts();
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        FTensor::Tensor1<double, 3> t_delta;
        t_delta(i) = t_center(i) - t_coords(i);
        FTensor::Tensor1<double, 3> t_disp;
        t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
 
        auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
        int bb = 0;
        for (; bb != nb_dofs / 3; ++bb) {
          t_nf(i) -= t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
 
        ++t_w;
        ++t_normal;
        ++t_coords;
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpRotationBc::iNtegrate(EntData &data) {
  return OP::iNtegrate(data);
}
 
MoFEMErrorCode OpBrokenTractionBc::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(3, i);
 
  int nb_dofs = data.getFieldData().size();
  int nb_integration_pts = getGaussPts().size2();
  int nb_base_functions = data.getN().size2();
  double ts_t = getFEMethod()->ts_t;
 
#ifndef NDEBUG
  if (this->locF.size() != nb_dofs)
    SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
             "Size of locF %d != nb_dofs %d", this->locF.size(), nb_dofs);
#endif // NDEBUG
 
  auto integrate_rhs = [&](auto &bc) {
    MoFEMFunctionBegin;
 
    auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
    auto t_row_base = data.getFTensor0N();
    auto t_w = getFTensor0IntegrationWeight();
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
      int rr = 0;
      for (; rr != nb_dofs / SPACE_DIM; ++rr) {
        t_f(i) -= ts_t * t_w * t_row_base * t_val(i);
        ++t_row_base;
        ++t_f;
      }
      for (; rr != nb_base_functions; ++rr)
        ++t_row_base;
      ++t_w;
    }
 
    this->locF *= getMeasure();
    MoFEMFunctionReturn(0);
  };
 
  auto integrate_rhs_cook = [&](auto &bc) {
    MoFEMFunctionBegin;
 
    auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
    auto t_row_base = data.getFTensor0N();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      auto calc_tau = [](double y) {
        y -= 44;
        y /= (60 - 44);
        return -y * (y - 1) / 0.25;
      };
 
      const auto tau = calc_tau(t_coords(1));
      auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
      int rr = 0;
      for (; rr != nb_dofs / SPACE_DIM; ++rr) {
        t_f(i) -= ts_t * t_w * t_row_base * tau * t_val(i);
        ++t_row_base;
        ++t_f;
      }
 
      for (; rr != nb_base_functions; ++rr)
        ++t_row_base;
      ++t_w;
      ++t_coords;
    }
 
    // Multiply by 2 since we integrate on triangle, and hybrid constrain is
    // multiplied by norm.
    this->locF *= 2. * getMeasure();
 
    MoFEMFunctionReturn(0);
  };
 
  // get entity of face
  EntityHandle fe_ent = getFEEntityHandle();
  for (auto &bc : *(bcData)) {
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      int nb_dofs = data.getFieldData().size();
      if (nb_dofs) {
 
        if (std::regex_match(bc.blockName, std::regex(".*COOK.*")))
          CHKERR integrate_rhs_cook(bc);
        else
          CHKERR integrate_rhs(bc);
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialEquilibrium_dw_dP::integrate(EntData &row_data,
                                                     EntData &col_data) {
  MoFEMFunctionBegin;
  int nb_integration_pts = row_data.getN().size1();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
        &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2));
  };
  FTensor::Index<'i', 3> i;
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    int rr = 0;
    for (; rr != row_nb_dofs / 3; ++rr) {
      auto t_col_diff_base_fun = col_data.getFTensor2DiffN<3, 3>(gg, 0);
      auto t_m = get_ftensor1(K, 3 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        double div_col_base = t_col_diff_base_fun(i, i);
        t_m(i) += a * t_row_base_fun * div_col_base;
        ++t_m;
        ++t_col_diff_base_fun;
      }
      ++t_row_base_fun;
    }
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base_fun;
    ++t_w;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialEquilibrium_dw_dw::integrate(EntData &row_data,
                                                     EntData &col_data) {
  MoFEMFunctionBegin;
 
  if (alphaW < std::numeric_limits<double>::epsilon() &&
      alphaRho < std::numeric_limits<double>::epsilon())
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = row_data.getN().size1();
  const int row_nb_dofs = row_data.getIndices().size();
  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
        &m(r + 0, c + 0), &m(r + 1, c + 1), &m(r + 2, c + 2)
 
    );
  };
  FTensor::Index<'i', 3> i;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  double ts_scale = alphaW * getTSa();
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
    ts_scale += alphaRho * getTSaa();
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w * ts_scale;
 
    int rr = 0;
    for (; rr != row_nb_dofs / 3; ++rr) {
 
      auto t_col_base_fun = row_data.getFTensor0N(gg, 0);
      auto t_m = get_ftensor1(K, 3 * rr, 0);
      for (int cc = 0; cc != row_nb_dofs / 3; ++cc) {
        const double b = a * t_row_base_fun * t_col_base_fun;
        t_m(i) -= b;
        ++t_m;
        ++t_col_base_fun;
      }
 
      ++t_row_base_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base_fun;
 
    ++t_w;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical_du_dP::integrate(EntData &row_data,
                                                  EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
 
  int nb_integration_pts = row_data.getN().size1();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, size_symm, 3>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
 
        &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
 
        &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
 
        &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2));
  };
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P_adjont_log_du_dP =
      getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    int rr = 0;
    for (; rr != row_nb_dofs / 6; ++rr) {
 
      auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
      auto t_m = get_ftensor3(K, 6 * rr, 0);
 
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        t_m(L, i) +=
            a * (t_approx_P_adjont_log_du_dP(i, j, L) * t_col_base_fun(j)) *
            t_row_base_fun;
        ++t_col_base_fun;
        ++t_m;
      }
 
      ++t_row_base_fun;
    }
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base_fun;
    ++t_w;
    ++t_approx_P_adjont_log_du_dP;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical_du_dBubble::integrate(EntData &row_data,
                                                       EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
 
  int nb_integration_pts = row_data.getN().size1();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, size_symm>(
        &m(r + 0, c), &m(r + 1, c), &m(r + 2, c), &m(r + 3, c), &m(r + 4, c),
        &m(r + 5, c));
  };
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P_adjont_log_du_dP =
      getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    int rr = 0;
    for (; rr != row_nb_dofs / 6; ++rr) {
      auto t_m = get_ftensor2(K, 6 * rr, 0);
      auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
      for (int cc = 0; cc != col_nb_dofs; ++cc) {
        t_m(L) +=
            a * (t_approx_P_adjont_log_du_dP(i, j, L) * t_col_base_fun(i, j)) *
            t_row_base_fun;
        ++t_m;
        ++t_col_base_fun;
      }
      ++t_row_base_fun;
    }
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base_fun;
    ++t_w;
    ++t_approx_P_adjont_log_du_dP;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical_du_domega::integrate(EntData &row_data,
                                                      EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = symm_L_tensor();
 
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
  auto get_ftensor3 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, size_symm, 3>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2),
 
        &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2),
 
        &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2),
 
        &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2)
 
    );
  };
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P_adjont_log_du_domega =
      getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  int nb_integration_pts = row_data.getN().size1();
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    int rr = 0;
    for (; rr != row_nb_dofs / 6; ++rr) {
      auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
      auto t_m = get_ftensor3(K, 6 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        double v = a * t_row_base_fun * t_col_base_fun;
        t_m(L, k) += v * t_approx_P_adjont_log_du_domega(k, L);
        ++t_m;
        ++t_col_base_fun;
      }
      ++t_row_base_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_approx_P_adjont_log_du_domega;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialRotation_domega_dP::integrate(EntData &row_data,
                                                      EntData &col_data) {
  MoFEMFunctionBegin;
  int nb_integration_pts = getGaussPts().size2();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
 
    );
  };
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_levi_kirchoff_dP =
      getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    int rr = 0;
    for (; rr != row_nb_dofs / 3; ++rr) {
      double b = a * t_row_base_fun;
      auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
      auto t_m = get_ftensor2(K, 3 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        t_m(k, i) += b * (t_levi_kirchoff_dP(i, l, k) * t_col_base_fun(l));
        ++t_m;
        ++t_col_base_fun;
      }
      ++t_row_base_fun;
    }
    for (; rr != row_nb_base_functions; ++rr) {
      ++t_row_base_fun;
    }
 
    ++t_w;
    ++t_levi_kirchoff_dP;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialRotation_domega_dBubble::integrate(EntData &row_data,
                                                           EntData &col_data) {
  MoFEMFunctionBegin;
  int nb_integration_pts = getGaussPts().size2();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
        &m(r + 0, c), &m(r + 1, c), &m(r + 2, c));
  };
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Index<'m', 3> m;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_levi_kirchoff_dP =
      getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    int rr = 0;
    for (; rr != row_nb_dofs / 3; ++rr) {
      double b = a * t_row_base_fun;
      auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
      auto t_m = get_ftensor1(K, 3 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs; ++cc) {
        t_m(k) += b * (t_levi_kirchoff_dP(i, j, k) * t_col_base_fun(i, j));
        ++t_m;
        ++t_col_base_fun;
      }
      ++t_row_base_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr) {
      ++t_row_base_fun;
    }
    ++t_w;
    ++t_levi_kirchoff_dP;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialRotation_domega_domega::integrate(EntData &row_data,
                                                          EntData &col_data) {
  MoFEMFunctionBegin;
  int nb_integration_pts = getGaussPts().size2();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
 
    );
  };
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_levi_kirchoff_domega =
      getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    int rr = 0;
    for (; rr != row_nb_dofs / 3; ++rr) {
      auto t_m = get_ftensor2(K, 3 * rr, 0);
      const double b = a * t_row_base_fun;
      auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        t_m(k, l) += (b * t_col_base_fun) * t_levi_kirchoff_domega(k, l);
        ++t_m;
        ++t_col_base_fun;
      }
      ++t_row_base_fun;
    }
    for (; rr != row_nb_base_functions; ++rr) {
      ++t_row_base_fun;
    }
    ++t_w;
    ++t_levi_kirchoff_domega;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialConsistency_dP_dP::integrate(EntData &row_data,
                                                     EntData &col_data) {
  MoFEMFunctionBeginHot;
  if (dataAtPts->matInvD.size1() == size_symm &&
      dataAtPts->matInvD.size2() == size_symm) {
    return integrateImpl<0>(row_data, col_data);
  } else {
    return integrateImpl<size_symm>(row_data, col_data);
  }
  MoFEMFunctionReturnHot(0);
};
 
template <int S>
MoFEMErrorCode OpSpatialConsistency_dP_dP::integrateImpl(EntData &row_data,
                                                         EntData &col_data) {
  MoFEMFunctionBegin;
 
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
 
    );
  };
 
  int nb_integration_pts = getGaussPts().size2();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  int row_nb_base_functions = row_data.getN().size2() / 3;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
 
  auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
      &*dataAtPts->matInvD.data().begin());
 
  auto t_row_base = row_data.getFTensor1N<3>();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    int rr = 0;
    for (; rr != row_nb_dofs / 3; ++rr) {
      auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
      auto t_m = get_ftensor2(K, 3 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        t_m(i, k) += a * t_row_base(j) * (t_inv_D(i, j, k, l) * t_col_base(l));
        ++t_m;
        ++t_col_base;
      }
 
      ++t_row_base;
    }
 
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base;
 
    ++t_w;
    ++t_inv_D;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
OpSpatialConsistency_dBubble_dBubble::integrate(EntData &row_data,
                                                EntData &col_data) {
  MoFEMFunctionBeginHot;
  if (dataAtPts->matInvD.size1() == size_symm &&
      dataAtPts->matInvD.size2() == size_symm) {
    return integrateImpl<0>(row_data, col_data);
  } else {
    return integrateImpl<size_symm>(row_data, col_data);
  }
  MoFEMFunctionReturnHot(0);
};
 
template <int S>
MoFEMErrorCode
OpSpatialConsistency_dBubble_dBubble::integrateImpl(EntData &row_data,
                                                    EntData &col_data) {
  MoFEMFunctionBegin;
 
  int nb_integration_pts = getGaussPts().size2();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  int row_nb_base_functions = row_data.getN().size2() / 9;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
 
  auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
      &*dataAtPts->matInvD.data().begin());
 
  auto t_row_base = row_data.getFTensor2N<3, 3>();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    int rr = 0;
    for (; rr != row_nb_dofs; ++rr) {
      auto t_col_base = col_data.getFTensor2N<3, 3>(gg, 0);
      for (int cc = 0; cc != col_nb_dofs; ++cc) {
        K(rr, cc) +=
            a * (t_row_base(i, j) * (t_inv_D(i, j, k, l) * t_col_base(k, l)));
        ++t_col_base;
      }
 
      ++t_row_base;
    }
 
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base;
    ++t_w;
    ++t_inv_D;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialConsistency_dBubble_dP::integrate(EntData &row_data,
                                                          EntData &col_data) {
  MoFEMFunctionBeginHot;
  if (dataAtPts->matInvD.size1() == size_symm &&
      dataAtPts->matInvD.size2() == size_symm) {
    return integrateImpl<0>(row_data, col_data);
  } else {
    return integrateImpl<size_symm>(row_data, col_data);
  }
  MoFEMFunctionReturnHot(0);
};
 
template <int S>
MoFEMErrorCode OpSpatialConsistency_dBubble_dP::integrateImpl(EntData &row_data,
                                                          EntData &col_data) {
  MoFEMFunctionBegin;
 
  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2)
 
    );
  };
 
  int nb_integration_pts = getGaussPts().size2();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  int row_nb_base_functions = row_data.getN().size2() / 9;
 
  i_FTIndex<SPACE_DIM> i;
  j_FTIndex<SPACE_DIM> j;
  k_FTIndex<SPACE_DIM> k;
  l_FTIndex<SPACE_DIM> l;
  m_FTIndex<SPACE_DIM> m;
  n_FTIndex<SPACE_DIM> n;
 
  auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
      &*dataAtPts->matInvD.data().begin());
 
  auto t_row_base = row_data.getFTensor2N<3, 3>();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    auto t_m = get_ftensor1(K, 0, 0);
 
    int rr = 0;
    for (; rr != row_nb_dofs; ++rr) {
      auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        t_m(k) += a * (t_row_base(i, j) * t_inv_D(i, j, k, l)) * t_col_base(l);
        ++t_col_base;
        ++t_m;
      }
 
      ++t_row_base;
    }
 
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base;
    ++t_w;
    ++t_inv_D;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialConsistency_dP_domega::integrate(EntData &row_data,
                                                         EntData &col_data) {
  MoFEMFunctionBegin;
 
  int nb_integration_pts = row_data.getN().size1();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
 
    );
  };
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
  int row_nb_base_functions = row_data.getN().size2() / 3;
  auto t_row_base_fun = row_data.getFTensor1N<3>();
 
  const double ts_a = getTSa();
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
 
    int rr = 0;
    for (; rr != row_nb_dofs / 3; ++rr) {
 
      FTensor::Tensor2<double, 3, 3> t_PRT;
      t_PRT(i, k) = t_row_base_fun(j) * t_h_domega(i, j, k);
 
      auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
      auto t_m = get_ftensor2(K, 3 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        t_m(i, j) += (a * t_col_base_fun) * t_PRT(i, j);
        ++t_m;
        ++t_col_base_fun;
      }
 
      ++t_row_base_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base_fun;
    ++t_w;
    ++t_h_domega;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
OpSpatialConsistency_dBubble_domega::integrate(EntData &row_data,
                                               EntData &col_data) {
  MoFEMFunctionBegin;
 
  int nb_integration_pts = row_data.getN().size1();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
        &m(r, c + 0), &m(r, c + 1), &m(r, c + 2));
  };
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
  int row_nb_base_functions = row_data.getN().size2() / 9;
  auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
 
  const double ts_a = getTSa();
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    int rr = 0;
    for (; rr != row_nb_dofs; ++rr) {
 
      FTensor::Tensor1<double, 3> t_PRT;
      t_PRT(k) = t_row_base_fun(i, j) * t_h_domega(i, j, k);
 
      auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
      auto t_m = get_ftensor2(K, rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        t_m(j) += (a * t_col_base_fun) * t_PRT(j);
        ++t_m;
        ++t_col_base_fun;
      }
 
      ++t_row_base_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_h_domega;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpPostProcDataStructure::doWork(int side, EntityType type,
                                               EntData &data) {
  MoFEMFunctionBegin;
 
  if(tagSense != getSkeletonSense())
    MoFEMFunctionReturnHot(0);
 
  auto create_tag = [this](const std::string tag_name, const int size) {
    double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
    Tag th;
    CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
                                       th, MB_TAG_CREAT | MB_TAG_SPARSE,
                                       def_VAL);
    return th;
  };
 
  Tag th_w = create_tag("SpatialDisplacement", 3);
  Tag th_omega = create_tag("Omega", 3);
  Tag th_approxP = create_tag("Piola1Stress", 9);
  Tag th_calcSigma = create_tag("EshelbyStress", 9);
  Tag th_sigma = create_tag("CauchyStress", 9);
  Tag th_log_u = create_tag("LogSpatialStretch", 9);
  Tag th_u = create_tag("SpatialStretch", 9);
  Tag th_h = create_tag("h", 9);
  Tag th_X = create_tag("X", 3);
  Tag th_detF = create_tag("detF", 1);
  Tag th_angular_momentum = create_tag("AngularMomentum", 3);
 
  // Tag th_u_eig_vec = create_tag("SpatialStretchEigenVec", 9);
  // Tag th_u_eig_vals = create_tag("SpatialStretchEigenVals", 3);
  Tag th_traction = create_tag("traction", 3);
 
  Tag th_disp = create_tag("H1Displacement", 3);
  Tag th_disp_error = create_tag("DisplacementError", 1);
  Tag th_lambda_disp = create_tag("ContactDisplacement", 3);
 
  Tag th_energy = create_tag("Energy", 1);
 
  auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
  auto t_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
  if (dataAtPts->SigmaAtPts.size2() != dataAtPts->approxPAtPts.size2()) {
    // that is for case that Eshelby stress is not calculated
    dataAtPts->SigmaAtPts.resize(dataAtPts->approxPAtPts.size1(),
                                  dataAtPts->approxPAtPts.size2(), false);
    dataAtPts->SigmaAtPts.clear();
  }
 
  auto t_calcSigma_P = getFTensor2FromMat<3, 3>(dataAtPts->SigmaAtPts);
  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
  auto t_coords = getFTensor1FromMat<3>(dataAtPts->XH1AtPts);
  auto t_normal = getFTensor1NormalsAtGaussPts();
  auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
  auto t_lambda_disp = getFTensor1FromMat<3>(dataAtPts->contactL2AtPts);
 
  if (dataAtPts->energyAtPts.size() == 0) {
    // that is for case that energy is not calculated
    dataAtPts->energyAtPts.resize(getGaussPts().size2());
    dataAtPts->energyAtPts.clear();
  }
  auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  auto set_float_precision = [](const double x) {
    if (std::abs(x) < std::numeric_limits<float>::epsilon())
      return 0.;
    else
      return x;
  };
 
  // scalars
  auto save_scal_tag = [&](auto &th, auto v, const int gg) {
    MoFEMFunctionBegin;
    v = set_float_precision(v);
    CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, &v);
    MoFEMFunctionReturn(0);
  };
 
  // vectors
  VectorDouble3 v(3);
  FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
  auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
    MoFEMFunctionBegin;
    t_v(i) = t_d(i);
    for (auto &a : v.data())
      a = set_float_precision(a);
    CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
                                     &*v.data().begin());
    MoFEMFunctionReturn(0);
  };
 
  // tensors
 
  MatrixDouble3by3 m(3, 3);
  FTensor::Tensor2<FTensor::PackPtr<double *, 0>, 3, 3> t_m(
      &m(0, 0), &m(0, 1), &m(0, 2),
 
      &m(1, 0), &m(1, 1), &m(1, 2),
 
      &m(2, 0), &m(2, 1), &m(2, 2));
 
  auto save_mat_tag = [&](auto &th, auto &t_d, const int gg) {
    MoFEMFunctionBegin;
    t_m(i, j) = t_d(i, j);
    for (auto &v : m.data())
      v = set_float_precision(v);
    CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
                                     &*m.data().begin());
    MoFEMFunctionReturn(0);
  };
 
  const auto nb_gauss_pts = getGaussPts().size2();
  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
 
    FTensor::Tensor1<double, 3> t_traction;
    t_traction(i) = t_approx_P(i, j) * t_normal(j) / t_normal.l2();
 
    // vectors
    CHKERR save_vec_tag(th_w, t_w, gg);
    CHKERR save_vec_tag(th_X, t_coords, gg);
    CHKERR save_vec_tag(th_omega, t_omega, gg);
    CHKERR save_vec_tag(th_traction, t_traction, gg);
 
    // tensors
    CHKERR save_mat_tag(th_h, t_h, gg);
 
    FTensor::Tensor2<double, 3, 3> t_log_u_tmp;
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj)
        t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
 
    CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
 
    FTensor::Tensor2<double, 3, 3> t_u_tmp;
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj)
        t_u_tmp(ii, jj) = t_u(ii, jj);
 
    CHKERR save_mat_tag(th_u, t_u_tmp, gg);
    CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
    CHKERR save_mat_tag(th_calcSigma, t_calcSigma_P, gg);
    CHKERR save_vec_tag(th_disp, t_disp, gg);
    CHKERR save_vec_tag(th_lambda_disp, t_lambda_disp, gg);
 
    double u_error = sqrt((t_disp(i) - t_w(i)) * (t_disp(i) - t_w(i)));
    CHKERR save_scal_tag(th_disp_error, u_error, gg);
    CHKERR save_scal_tag(th_energy, t_energy, gg);
 
    const double jac = determinantTensor3by3(t_h);
    FTensor::Tensor2<double, 3, 3> t_cauchy;
    t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
    CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
    CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
 
    FTensor::Tensor1<double, 3> t_levi;
    t_levi(k) = t_levi_kirchoff(k);
    CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
                                     &t_levi(0));
 
    auto get_eiegn_vector_symmetric = [&](auto &t_u) {
      MoFEMFunctionBegin;
 
      // for (int ii = 0; ii != 3; ++ii)
      //   for (int jj = 0; jj != 3; ++jj)
      //     t_m(ii, jj) = t_u(ii, jj);
 
      // VectorDouble3 eigen_values(3);
      // auto t_eigen_values = getFTensor1FromArray<3>(eigen_values);
      // if (computeEigenValuesSymmetric(t_m, t_eigen_values) != 0) {
      // };
 
      // CHKERR postProcMesh.tag_set_data(th_u_eig_vec, &mapGaussPts[gg], 1,
      //                                  &*m.data().begin());
      // CHKERR postProcMesh.tag_set_data(th_u_eig_vals, &mapGaussPts[gg], 1,
      //                                  &*eigen_values.data().begin());
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR get_eiegn_vector_symmetric(t_u);
 
    ++t_w;
    ++t_h;
    ++t_log_u;
    ++t_u;
    ++t_omega;
    ++t_R;
    ++t_approx_P;
    ++t_calcSigma_P;
    ++t_levi_kirchoff;
    ++t_coords;
    ++t_normal;
    ++t_disp;
    ++t_lambda_disp;
    ++t_energy;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode AddHOOps<2, 3, 3>::add(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
    std::vector<FieldSpace> spaces, std::string geom_field_name,
    boost::shared_ptr<Range> crack_front_edges_ptr) {
  MoFEMFunctionBegin;
 
  {
 
    struct OpGetHONormalsOnFace
        : public MoFEM::OpGetHONormalsOnFace<SPACE_DIM> {
 
      using OP = MoFEM::OpGetHONormalsOnFace<SPACE_DIM>;
 
      OpGetHONormalsOnFace(std::string &field_name,
                           boost::shared_ptr<Range> edges_ptr)
          : OP(field_name), edgesPtr(edges_ptr) {}
 
      MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
 
        auto ent = data.getFieldEntities().size()
                       ? data.getFieldEntities()[0]->getEnt()
                       : 0;
 
        if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
          return 0;
        } else {
          return OP::doWork(side, type, data);
        }
      };
 
    private:
      boost::shared_ptr<Range> edgesPtr;
    };
 
    auto jac = boost::make_shared<MatrixDouble>();
    auto det = boost::make_shared<VectorDouble>();
    pipeline.push_back(new OpGetHONormalsOnFace(
        geom_field_name, EshelbianCore::setSingularity
                             ? crack_front_edges_ptr
                             : boost::make_shared<Range>()));
    pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac));
    pipeline.push_back(new OpInvertMatrix<3>(jac, det, nullptr));
    pipeline.push_back(new OpScaleBaseBySpaceInverseOfMeasure(
        L2, AINSWORTH_LEGENDRE_BASE, det));
  }
 
  CHKERR MoFEM::AddHOOps<2, 3, 3>::add(pipeline, spaces, geom_field_name);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode AddHOOps<3, 3, 3>::add(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
    std::vector<FieldSpace> spaces, std::string geom_field_name,
    boost::shared_ptr<Range> crack_front_edges_ptr) {
  MoFEMFunctionBegin;
 
  {
    auto jac = boost::make_shared<MatrixDouble>();
    auto det = boost::make_shared<VectorDouble>();
    pipeline.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        geom_field_name, jac));
    pipeline.push_back(new OpInvertMatrix<3>(jac, det, nullptr));
    pipeline.push_back(
        new OpScaleBaseBySpaceInverseOfMeasure(L2, USER_BASE, det));
  }
 
  {
 
    struct OpCalculateVectorFieldGradient
        : public MoFEM::OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM> {
 
      using OP = MoFEM::OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>;
 
      OpCalculateVectorFieldGradient(const std::string &field_name,
                                     boost::shared_ptr<MatrixDouble> jac,
                                     boost::shared_ptr<Range> edges_ptr)
          : OP(field_name, jac), edgesPtr(edges_ptr) {}
 
      MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
 
        auto ent = data.getFieldEntities().size()
                       ? data.getFieldEntities()[0]->getEnt()
                       : 0;
 
        if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
          return 0;
        } else {
          return OP::doWork(side, type, data);
        }
      };
 
    private:
      boost::shared_ptr<Range> edgesPtr;
    };
 
    auto jac = boost::make_shared<MatrixDouble>();
    auto det = boost::make_shared<VectorDouble>();
    pipeline.push_back(new OpCalculateVectorFieldGradient(
        geom_field_name, jac,
        EshelbianCore::setSingularity ? crack_front_edges_ptr
                                      : boost::make_shared<Range>()));
    pipeline.push_back(new OpInvertMatrix<3>(jac, det, nullptr));
    pipeline.push_back(new OpScaleBaseBySpaceInverseOfMeasure(
        L2, AINSWORTH_LEGENDRE_BASE, det));
  }
 
  CHKERR MoFEM::AddHOOps<3, 3, 3>::add(pipeline, spaces, geom_field_name,
                                       nullptr, nullptr, nullptr);
 
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity