EshelbianOperators.cpp

Go to the documentation of this file.

/**
 * \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 <lapack_wrap.h>
 
#include <MatrixFunction.hpp>
 
namespace EshelbianPlasticity {
 
inline auto diff_deviator(FTensor::Ddg<double, 3, 3> &&t_diff_stress) {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  FTensor::Ddg<double, 3, 3> t_diff_deviator;
  t_diff_deviator(i, j, k, l) = t_diff_stress(i, j, k, l);
 
  constexpr double third = boost::math::constants::third<double>();
 
  t_diff_deviator(0, 0, 0, 0) -= third;
  t_diff_deviator(0, 0, 1, 1) -= third;
 
  t_diff_deviator(1, 1, 0, 0) -= third;
  t_diff_deviator(1, 1, 1, 1) -= third;
 
  t_diff_deviator(2, 2, 0, 0) -= third;
  t_diff_deviator(2, 2, 1, 1) -= third;
 
  t_diff_deviator(0, 0, 2, 2) -= third;
  t_diff_deviator(1, 1, 2, 2) -= third;
  t_diff_deviator(2, 2, 2, 2) -= third;
 
  return t_diff_deviator;
}
 
constexpr static auto size_symm = (3 * (3 + 1)) / 2;
 
inline auto diff_tensor() {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Ddg<double, 3, 3> t_diff;
  constexpr auto t_kd = FTensor::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() {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'L', size_symm> L;
  FTensor::Dg<double, 3, size_symm> t_L;
  t_L(i, j, L) = 0;
  t_L(0, 0, 0) = 1;
  t_L(1, 1, 3) = 1;
  t_L(2, 2, 5) = 1;
  t_L(1, 0, 1) = 1;
  t_L(2, 0, 2) = 1;
  t_L(2, 1, 4) = 1;
  return t_L;
}
 
template <class T> struct TensorTypeExtractor {
  typedef typename std::remove_pointer<T>::type Type;
};
template <class T, int I> struct TensorTypeExtractor<FTensor::PackPtr<T, I>> {
  typedef typename std::remove_pointer<T>::type Type;
};
 
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 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();
  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<'l', 3> 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->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_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;
 
    t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
 
    auto calulate_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 calulate_streach = [&]() {
      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);
    };
 
    calulate_rotation();
    calulate_streach();
 
    switch (EshelbianCore::gradApperoximator) {
    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;
    }
 
    auto t_log_u2_h1_tmp =
        EigenMatrix::getMat(eig, eigen_vec, EshelbianCore::inv_f);
 
    switch (EshelbianCore::gradApperoximator) {
    case LARGE_ROT:
    case MODERATE_ROT:
      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_stretch_total(i, j) = t_log_u(i, j);
      break;
    };
 
    ++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_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 OpSpatialPhysical::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
    CHKERR integratePiola(data);
  } else {
    if (polyConvex) {
      CHKERR integratePolyconvexHencky(data);
    } else {
      CHKERR integrateHencky(data);
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical::integrateHencky(EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = symm_L_tensor();
 
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P_adjont_log_du =
      getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
  auto t_log_streach_h1 =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
  auto t_dot_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
 
  auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  auto get_ftensor2 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 6>, 6>(
        &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
  };
 
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor2(nF);
 
    FTensor::Tensor2_symmetric<double, 3> t_T;
    t_T(i, j) =
        t_D(i, j, k, l) * (t_log_streach_h1(k, l) + alphaU * t_dot_log_u(k, l));
    FTensor::Tensor1<double, size_symm> t_residual;
    t_residual(L) =
        a * (t_approx_P_adjont_log_du(L) - t_L(i, j, L) * t_T(i, j));
 
    int bb = 0;
    for (; bb != nb_dofs / 6; ++bb) {
      t_nf(L) += t_row_base_fun * t_residual(L);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_approx_P_adjont_log_du;
    ++t_dot_log_u;
    ++t_log_streach_h1;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical::integratePolyconvexHencky(EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = symm_L_tensor();
 
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P_adjont_log_du =
      getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
  auto t_log_streach_h1 =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
  auto t_dot_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
 
  auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  auto get_ftensor2 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 6>, 6>(
        &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
  };
 
  constexpr double nohat_k = 1. / 4;
  constexpr double hat_k = 1. / 8;
  double mu = dataAtPts->mu;
  double lambda = dataAtPts->lambda;
 
  constexpr double third = boost::math::constants::third<double>();
  auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
  auto t_diff_deviator = diff_deviator(diff_tensor());
 
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor2(nF);
 
    double log_det = t_log_streach_h1(i, i);
    double log_det2 = log_det * log_det;
    FTensor::Tensor2_symmetric<double, 3> t_dev;
    t_dev(i, j) = t_log_streach_h1(i, j) - t_kd(i, j) * (third * log_det);
    double dev_norm2 = t_dev(i, j) * t_dev(i, j);
 
    FTensor::Tensor2_symmetric<double, 3> t_T;
    auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
    auto B = lambda * std::exp(hat_k * log_det2) * log_det;
    t_T(i, j) =
 
        A * (t_dev(k, l) * t_diff_deviator(k, l, i, j))
 
        +
 
        B * t_kd(i, j)
 
        +
 
        alphaU * t_D(i, j, k, l) * t_log_streach_h1(k, l);
 
    FTensor::Tensor1<double, size_symm> t_residual;
    t_residual(L) =
        a * (t_approx_P_adjont_log_du(L) - t_L(i, j, L) * t_T(i, j));
 
    int bb = 0;
    for (; bb != nb_dofs / 6; ++bb) {
      t_nf(L) += t_row_base_fun * t_residual(L);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_approx_P_adjont_log_du;
    ++t_dot_log_u;
    ++t_log_streach_h1;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical::integratePiola(EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = symm_L_tensor();
 
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P_adjont_log_du =
      getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
  auto t_dot_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
  auto t_diff_u =
      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  auto get_ftensor2 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 6>, 6>(
        &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
  };
 
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor2(nF);
 
    FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
    t_Ldiff_u(i, j, L) = t_diff_u(i, j, k, l) * t_L(k, l, L);
 
    FTensor::Tensor2_symmetric<double, 3> t_viscous_P;
    t_viscous_P(i, j) =
        alphaU *
        t_dot_log_u(i,
                    j); // That is chit, should be split on axiator and deviator
 
    FTensor::Tensor1<double, size_symm> t_residual;
    t_residual(L) =
        a * (t_approx_P_adjont_log_du(L) - t_Ldiff_u(i, j, L) * t_P(i, j) -
             t_L(i, j, L) * t_viscous_P(i, j));
 
    int bb = 0;
    for (; bb != nb_dofs / 6; ++bb) {
      t_nf(L) += t_row_base_fun * t_residual(L);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_approx_P_adjont_log_du;
    ++t_P;
    ++t_dot_log_u;
    ++t_diff_u;
  }
  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);
  auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
 
  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;
    ++t_omega_dot;
  }
 
  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);
  auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
 
  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;
    ++t_omega_dot;
  }
 
  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);
}
 
MoFEMErrorCode OpDispBc::integrate(EntData &data) {
  MoFEMFunctionBegin;
  // get entity of face
  EntityHandle fe_ent = getFEEntityHandle();
  // interate 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 = data.getN().size1();
      auto t_normal = getFTensor1Normal();
      auto t_w = getFTensor0IntegrationWeight();
      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;
      auto get_ftensor1 = [](auto &v) {
        return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
                                                                  &v[2]);
      };
 
      double scale = 1;
      for (auto &sm : scalingMethodsVec) {
        scale *= sm->getScale(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) {
        FTensor::Tensor1<double, 3> t_bc_res;
        t_bc_res(i) = t_bc_disp(i);
        auto t_nf = get_ftensor1(nF);
        int bb = 0;
        for (; bb != nb_dofs / 3; ++bb) {
          t_nf(i) -=
              t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_res(i) * 0.5;
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
 
        ++t_w;
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpRotationBc::integrate(EntData &data) {
  MoFEMFunctionBegin;
  // get entity of face
  EntityHandle fe_ent = 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 = data.getN().size1();
      auto t_normal = getFTensor1Normal();
      auto t_w = getFTensor0IntegrationWeight();
 
      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;
      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 *= 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 = 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 = get_ftensor1(nF);
        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_coords;
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode FeTractionBc::preProcess() {
  MoFEMFunctionBegin;
 
  struct FaceRule {
    int operator()(int p_row, int p_col, int p_data) const {
      return 2 * (p_data + 1);
    }
  };
 
  if (ts_ctx == CTX_TSSETIFUNCTION) {
 
    // Loop boundary elements with traction boundary conditions
    FaceElementForcesAndSourcesCore fe(mField);
    AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(fe.getOpPtrVector(),
                                                       {HDIV});
    fe.getOpPtrVector().push_back(
        new OpTractionBc(std::string("P") /* + "_RT"*/, *this));
    fe.getRuleHook = FaceRule();
    fe.ts_t = ts_t;
    CHKERR mField.loop_finite_elements(problemPtr->getName(), "ESSENTIAL_BC",
                                       fe, nullptr, MF_ZERO,
                                       this->getCacheWeakPtr());
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
 
  auto t_normal = getFTensor1Normal();
  const double nrm2 = sqrt(t_normal(i) * t_normal(i));
  FTensor::Tensor1<double, 3> t_unit_normal;
  t_unit_normal(i) = t_normal(i) / nrm2;
  int nb_dofs = data.getFieldData().size();
  int nb_integration_pts = data.getN().size1();
  int nb_base_functions = data.getN().size2() / 3;
  double ts_t = getFEMethod()->ts_t;
 
  auto integrate_matrix = [&]() {
    MoFEMFunctionBegin;
 
    auto t_w = getFTensor0IntegrationWeight();
    matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
    matPP.clear();
 
    auto t_row_base_fun = data.getFTensor1N<3>();
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      int rr = 0;
      for (; rr != nb_dofs / 3; ++rr) {
        const double a = t_row_base_fun(i) * t_unit_normal(i);
        auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
        for (int cc = 0; cc != nb_dofs / 3; ++cc) {
          const double b = t_col_base_fun(i) * t_unit_normal(i);
          matPP(rr, cc) += t_w * a * b;
          ++t_col_base_fun;
        }
        ++t_row_base_fun;
      }
 
      for (; rr != nb_base_functions; ++rr)
        ++t_row_base_fun;
 
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto integrate_rhs = [&](auto &bc) {
    MoFEMFunctionBegin;
 
    auto t_w = getFTensor0IntegrationWeight();
    vecPv.resize(3, nb_dofs / 3, false);
    vecPv.clear();
 
    auto t_row_base_fun = data.getFTensor1N<3>();
    double ts_t = getFEMethod()->ts_t;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      int rr = 0;
      for (; rr != nb_dofs / 3; ++rr) {
        const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
        for (int dd = 0; dd != 3; ++dd)
          if (bc.flags[dd])
            vecPv(dd, rr) += t * bc.vals[dd];
        ++t_row_base_fun;
      }
      for (; rr != nb_base_functions; ++rr)
        ++t_row_base_fun;
      ++t_w;
    }
    MoFEMFunctionReturn(0);
  };
 
  auto integrate_rhs_cook = [&](auto &bc) {
    MoFEMFunctionBegin;
 
    vecPv.resize(3, nb_dofs / 3, false);
    vecPv.clear();
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_row_base_fun = data.getFTensor1N<3>();
    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 double tau = calc_tau(t_coords(1));
 
      int rr = 0;
      for (; rr != nb_dofs / 3; ++rr) {
        const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
        for (int dd = 0; dd != 3; ++dd)
          if (bc.flags[dd])
            vecPv(dd, rr) += t * tau * bc.vals[dd];
        ++t_row_base_fun;
      }
 
      for (; rr != nb_base_functions; ++rr)
        ++t_row_base_fun;
      ++t_w;
      ++t_coords;
    }
    MoFEMFunctionReturn(0);
  };
 
  // get entity of face
  EntityHandle fe_ent = getFEEntityHandle();
  for (auto &bc : *(bcFe.bcData)) {
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      int nb_dofs = data.getFieldData().size();
      if (nb_dofs) {
 
        CHKERR integrate_matrix();
        if (std::regex_match(bc.blockName, std::regex(".*COOK.*")))
          CHKERR integrate_rhs_cook(bc);
        else
          CHKERR integrate_rhs(bc);
 
        const auto info =
            lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
                         nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
        if (info > 0)
          SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                   "The leading minor of order %i is "
                   "not positive; definite;\nthe "
                   "solution could not be computed",
                   info);
 
        for (int dd = 0; dd != 3; ++dd)
          if (bc.flags[dd])
            for (int rr = 0; rr != nb_dofs / 3; ++rr)
              data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
      }
    }
  }
 
  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_du::integrate(EntData &row_data,
                                                  EntData &col_data) {
  MoFEMFunctionBegin;
  if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
    CHKERR integratePiola(row_data, col_data);
  } else {
    if (polyConvex) {
      CHKERR integratePolyconvexHencky(row_data, col_data);
    } else {
      CHKERR integrateHencky(row_data, col_data);
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical_du_du::integrateHencky(EntData &row_data,
                                                        EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  FTensor::Index<'J', size_symm> J;
  auto t_L = symm_L_tensor();
  auto t_diff = diff_tensor();
 
  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 *, 6>, size_symm,
                            size_symm>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
        &m(r + 0, c + 4), &m(r + 0, c + 5),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
        &m(r + 1, c + 4), &m(r + 1, c + 5),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
        &m(r + 2, c + 4), &m(r + 2, c + 5),
 
        &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
        &m(r + 3, c + 4), &m(r + 3, c + 5),
 
        &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
        &m(r + 4, c + 4), &m(r + 4, c + 5),
 
        &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
        &m(r + 5, c + 4), &m(r + 5, c + 5)
 
    );
  };
  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_approx_P_adjont_dstretch =
      getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
  auto &nbUniq = dataAtPts->nbUniq;
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  auto get_dP = [&]() {
    dP.resize(size_symm * size_symm, nb_integration_pts, false);
    auto ts_a = getTSa();
 
    auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
    FTensor::Tensor2<double, size_symm, size_symm> t_dP_tmp;
    t_dP_tmp(L, J) = -(1 + alphaU * ts_a) *
                     (t_L(i, j, L) *
                      ((t_D(i, j, m, n) * t_diff(m, n, k, l)) * t_L(k, l, J)));
 
 
    if (EshelbianCore::stretchSelector > LINEAR) {
      auto t_approx_P_adjont_dstretch =
          getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
      auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
      auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
      auto &nbUniq = dataAtPts->nbUniq;
 
      auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
      for (auto gg = 0; gg != nb_integration_pts; ++gg) {
 
        // Work of symmetric tensor on undefined tensor is equal to the work
        // of the symmetric part of it
        FTensor::Tensor2_symmetric<double, 3> t_sym;
        t_sym(i, j) = (t_approx_P_adjont_dstretch(i, j) ||
                       t_approx_P_adjont_dstretch(j, i));
        t_sym(i, j) /= 2.0;
        auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_sym, nbUniq[gg]);
        t_dP(L, J) = t_L(i, j, L) *
                         ((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
                          t_L(k, l, J)) /
                         2. +
                     t_dP_tmp(L, J);
 
        ++t_dP;
        ++t_approx_P_adjont_dstretch;
        ++t_eigen_vals;
        ++t_eigen_vecs;
      }
    } else {
      auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
      for (auto gg = 0; gg != nb_integration_pts; ++gg) {
        t_dP(L, J) = t_dP_tmp(L, J);
        ++t_dP;
      }
    }
 
    return getFTensor2FromMat<size_symm, size_symm>(dP);
  };
 
  auto t_dP = get_dP();
  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_ftensor2(K, 6 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
        const double b = a * t_row_base_fun * t_col_base_fun;
        t_m(L, J) += b * t_dP(L, 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_dP;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
OpSpatialPhysical_du_du::integratePolyconvexHencky(EntData &row_data,
                                                   EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  FTensor::Index<'J', size_symm> J;
  auto t_L = symm_L_tensor();
  auto t_diff = diff_tensor();
 
  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 *, 6>, size_symm,
                            size_symm>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
        &m(r + 0, c + 4), &m(r + 0, c + 5),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
        &m(r + 1, c + 4), &m(r + 1, c + 5),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
        &m(r + 2, c + 4), &m(r + 2, c + 5),
 
        &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
        &m(r + 3, c + 4), &m(r + 3, c + 5),
 
        &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
        &m(r + 4, c + 4), &m(r + 4, c + 5),
 
        &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
        &m(r + 5, c + 4), &m(r + 5, c + 5)
 
    );
  };
  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();
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  auto get_dP = [&]() {
    dP.resize(size_symm * size_symm, nb_integration_pts, false);
    auto ts_a = getTSa();
 
    auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
 
    constexpr double nohat_k = 1. / 4;
    constexpr double hat_k = 1. / 8;
    double mu = dataAtPts->mu;
    double lambda = dataAtPts->lambda;
 
    constexpr double third = boost::math::constants::third<double>();
    auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
    auto t_diff_deviator = diff_deviator(diff_tensor());
 
    auto t_approx_P_adjont_dstretch =
        getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
    auto t_log_streach_h1 =
        getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
    auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
    auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
    auto &nbUniq = dataAtPts->nbUniq;
 
    auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
 
      double log_det = t_log_streach_h1(i, i);
      double log_det2 = log_det * log_det;
      FTensor::Tensor2_symmetric<double, 3> t_dev;
      t_dev(i, j) = t_log_streach_h1(i, j) - t_kd(i, j) * (third * log_det);
      double dev_norm2 = t_dev(i, j) * t_dev(i, j);
 
      auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
      auto B = lambda * std::exp(hat_k * log_det2) * log_det;
 
      FTensor::Tensor2_symmetric<double, 3> t_A_diff, t_B_diff;
      t_A_diff(i, j) =
          (A * 2 * nohat_k) * (t_dev(k, l) * t_diff_deviator(k, l, i, j));
      t_B_diff(i, j) = (B * 2 * hat_k) * log_det * t_kd(i, j) +
                       lambda * std::exp(hat_k * log_det2) * t_kd(i, j);
      FTensor::Ddg<double, 3, 3> t_dT;
      t_dT(i, j, k, l) =
          t_A_diff(i, j) * (t_dev(m, n) * t_diff_deviator(m, n, k, l))
 
          +
 
          A * t_diff_deviator(m, n, i, j) * t_diff_deviator(m, n, k, l)
 
          +
 
          t_B_diff(i, j) * t_kd(k, l);
 
      t_dP(L, J) = -t_L(i, j, L) *
                   ((
 
                        t_dT(i, j, k, l)
 
                        +
 
                        (alphaU * ts_a) * (t_D(i, j, m, n) * t_diff(m, n, k, l)
 
                                               )) *
                    t_L(k, l, J));
 
      // Work of symmetric tensor on undefined tensor is equal to the work
      // of the symmetric part of it
      if (EshelbianCore::stretchSelector > LINEAR) {
        FTensor::Tensor2_symmetric<double, 3> t_sym;
        t_sym(i, j) = (t_approx_P_adjont_dstretch(i, j) ||
                       t_approx_P_adjont_dstretch(j, i));
        t_sym(i, j) /= 2.0;
        auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_sym, nbUniq[gg]);
        t_dP(L, J) += t_L(i, j, L) *
                      ((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
                       t_L(k, l, J)) /
                      2.;
      }
 
      ++t_dP;
      ++t_approx_P_adjont_dstretch;
      ++t_log_streach_h1;
      ++t_eigen_vals;
      ++t_eigen_vecs;
    }
 
    return getFTensor2FromMat<size_symm, size_symm>(dP);
  };
 
  auto t_dP = get_dP();
  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_ftensor2(K, 6 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
        const double b = a * t_row_base_fun * t_col_base_fun;
        t_m(L, J) += b * t_dP(L, 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_dP;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical_du_du::integratePiola(EntData &row_data,
                                                       EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  FTensor::Index<'J', size_symm> J;
  auto t_L = symm_L_tensor();
  auto t_diff = diff_tensor();
 
  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 *, 6>, size_symm,
                            size_symm>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
        &m(r + 0, c + 4), &m(r + 0, c + 5),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
        &m(r + 1, c + 4), &m(r + 1, c + 5),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
        &m(r + 2, c + 4), &m(r + 2, c + 5),
 
        &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
        &m(r + 3, c + 4), &m(r + 3, c + 5),
 
        &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
        &m(r + 4, c + 4), &m(r + 4, c + 5),
 
        &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
        &m(r + 5, c + 4), &m(r + 5, c + 5)
 
    );
  };
  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 get_dP = [&]() {
    dP.resize(size_symm * size_symm, nb_integration_pts, false);
    auto ts_a = getTSa();
 
    auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
    auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(dataAtPts->P_du);
    auto t_approx_P_adjont_dstretch =
        getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
    auto t_dot_log_u =
        getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
    auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
    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_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
 
      FTensor::Tensor2<double, 3, 3> t_deltaP;
      t_deltaP(i, j) = t_approx_P_adjont_dstretch(i, j) - t_P(i, j);
 
      FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
      t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
      t_dP(L, J) =
          -t_Ldiff_u(i, j, L) * (r_P_du(i, j, k, l) * t_Ldiff_u(k, l, J));
      // viscous stress derivative
      t_dP(L, J) -= (alphaU * ts_a) *
                    (t_L(i, j, L) * (t_diff(i, j, k, l) * t_L(k, l, J)));
 
      if (EshelbianCore::stretchSelector > LINEAR) {
        FTensor::Tensor2_symmetric<double, 3> t_deltaP_sym;
        t_deltaP_sym(i, j) = (t_deltaP(i, j) || t_deltaP(j, i));
        t_deltaP_sym(i, j) /= 2.0;
        auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_deltaP_sym, nbUniq[gg]);
        t_dP(L, J) += t_L(i, j, L) * (t_diff2_uP2(i, j, k, l) * t_L(k, l, J));
      }
 
      ++t_P;
      ++t_dP;
      ++r_P_du;
      ++t_approx_P_adjont_dstretch;
      ++t_dot_log_u;
      ++t_u;
      ++t_diff_u;
      ++t_eigen_vals;
      ++t_eigen_vecs;
    }
 
    return getFTensor2FromMat<size_symm, size_symm>(dP);
  };
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  auto t_dP = get_dP();
  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_ftensor2(K, 6 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
        const double b = a * t_row_base_fun * t_col_base_fun;
        t_m(L, J) += b * t_dP(L, 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_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_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;
 
  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_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);
 
  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);
  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
  auto t_coords = getFTensor1CoordsAtGaussPts();
  auto t_normal = getFTensor1NormalsAtGaussPts();
  auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
  auto t_lambda_disp = getFTensor1FromMat<3>(dataAtPts->contactL2AtPts);
 
  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_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);
 
    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);
      CHKERR computeEigenValuesSymmetric(t_m, t_eigen_values);
 
      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_levi_kirchoff;
    ++t_coords;
    ++t_normal;
    ++t_disp;
    ++t_lambda_disp;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
OpCalculateStrainEnergy::doWork(int side, EntityType type,
                                EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
  if (type == MBTET) {
    int nb_integration_pts = data.getN().size1();
    auto v = getVolume();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
    auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
 
    FTensor::Index<'i', 3> i;
    FTensor::Index<'J', 3> J;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      const double a = t_w * v;
      // FIXME: this is wrong, energy should be calculated in material
      (*energy) += a * t_P(i, J) * t_h(i, J);
      ++t_w;
      ++t_P;
      ++t_h;
    }
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity