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 <BasicFiniteElements.hpp>
 
#include <EshelbianPlasticity.hpp>
#include <boost/math/constants/constants.hpp>
 
#include <lapack_wrap.h>
 
#include <MatrixFunction.hpp>
 
namespace EshelbianPlasticity {
 
double f(const double v) { return exp(v); }
double d_f(const double v) { return exp(v); }
double dd_f(const double v) { return exp(v); }
 
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>
FTensor::Tensor2<typename TensorTypeExtractor<T>::Type, 3,
                 3> static get_rotation_form_vector(FTensor::Tensor1<T, 3>
                                                        &t_omega) {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Tensor2<typename TensorTypeExtractor<T>::Type, 3, 3> t_R;
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
  t_R(i, j) = t_kd(i, j);
 
  const typename TensorTypeExtractor<T>::Type angle =
      sqrt(t_omega(i) * t_omega(i));
 
  constexpr double tol = 1e-18;
  if (std::abs(angle) < tol) {
    return t_R;
  }
 
  FTensor::Tensor2<typename TensorTypeExtractor<T>::Type, 3, 3> t_Omega;
  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
  const typename TensorTypeExtractor<T>::Type a = sin(angle) / angle;
  const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
  const typename TensorTypeExtractor<T>::Type b =
      2. * ss_2 * ss_2 / (angle * angle);
  t_R(i, j) += a * t_Omega(i, j);
  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
 
  return t_R;
}
 
template <typename T>
FTensor::Tensor3<typename TensorTypeExtractor<T>::Type, 3, 3,
                 3> static get_diff_rotation_form_vector(FTensor::Tensor1<T, 3>
                                                             &t_omega) {
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  const typename TensorTypeExtractor<T>::Type angle =
      sqrt(t_omega(i) * t_omega(i));
 
  constexpr double tol = 1e-18;
  if (std::abs(angle) < tol) {
    FTensor::Tensor3<typename TensorTypeExtractor<T>::Type, 3, 3, 3> t_diff_R;
    auto t_R = get_rotation_form_vector(t_omega);
    t_diff_R(i, j, k) = FTensor::levi_civita<double>(i, l, k) * t_R(l, j);
    return t_diff_R;
  }
 
  FTensor::Tensor3<typename TensorTypeExtractor<T>::Type, 3, 3, 3> t_diff_R;
  FTensor::Tensor2<typename TensorTypeExtractor<T>::Type, 3, 3> t_Omega;
  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
 
  const typename TensorTypeExtractor<T>::Type angle2 = angle * angle;
  const typename TensorTypeExtractor<T>::Type ss = sin(angle);
  const typename TensorTypeExtractor<T>::Type a = ss / angle;
  const typename TensorTypeExtractor<T>::Type ss_2 = sin(angle / 2.);
  const typename TensorTypeExtractor<T>::Type b = 2. * ss_2 * ss_2 / angle2;
  t_diff_R(i, j, k) = 0;
  t_diff_R(i, j, k) = a * FTensor::levi_civita<double>(i, j, k);
  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 typename TensorTypeExtractor<T>::Type cc = cos(angle);
  const typename TensorTypeExtractor<T>::Type cc_2 = cos(angle / 2.);
  const typename TensorTypeExtractor<T>::Type diff_a =
      (angle * cc - ss) / angle2;
 
  const typename TensorTypeExtractor<T>::Type diff_b =
      (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / 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 T1, typename T2, typename T3>
inline MoFEMErrorCode get_eigen_val_and_proj_lapack(T1 &X, T2 &eig,
                                                    T3 &eig_vec) {
  MoFEMFunctionBeginHot;
  for (int ii = 0; ii != 3; ii++)
    for (int jj = 0; jj != 3; jj++)
      eig_vec(ii, jj) = X(ii, jj);
 
  int n = 3;
  int lda = 3;
  int lwork = (3 + 2) * 3;
  std::array<double, (3 + 2) * 3> work;
 
  if (lapack_dsyev('V', 'U', n, &(eig_vec(0, 0)), lda, &eig(0), work.data(),
                   lwork) > 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "The algorithm failed to compute eigenvalues.");
  MoFEMFunctionReturnHot(0);
}
 
inline bool is_eq(const double &a, const double &b) {
  constexpr double eps = 1e-8;
  return abs(a - b) < eps;
}
 
template <typename T> inline int get_uniq_nb(T &ec) {
  return distance(&ec(0), unique(&ec(0), &ec(0) + 3, is_eq));
}
 
template <typename T1, typename T2>
void sort_eigen_vals(T1 &eig, T2 &eigen_vec) {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  if (is_eq(eig(0), eig(1))) {
    FTensor::Tensor2<double, 3, 3> eigen_vec_c{
        eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
        eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2),
        eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2)};
    FTensor::Tensor1<double, 3> eig_c{eig(0), eig(2), eig(1)};
    eig(i) = eig_c(i);
    eigen_vec(i, j) = eigen_vec_c(i, j);
  } else {
    FTensor::Tensor2<double, 3, 3> eigen_vec_c{
        eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2),
        eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
        eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2)};
    FTensor::Tensor1<double, 3> eig_c{eig(1), eig(0), eig(2)};
    eig(i) = eig_c(i);
    eigen_vec(i, j) = eigen_vec_c(i, j);
  }
}
 
MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork(int side,
                                                             EntityType type,
                                                             EntData &data) {
  MoFEMFunctionBegin;
  if (type != MBTET)
    MoFEMFunctionReturnHot(0);
  int nb_integration_pts = data.getN().size1();
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  dataAtPts->hAtPts.resize(9, nb_integration_pts, false);
  dataAtPts->rotMatAtPts.resize(9, nb_integration_pts, false);
  dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts, false);
  dataAtPts->streachTensorAtPts.resize(6, nb_integration_pts, false);
  dataAtPts->diffStreachTensorAtPts.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);
 
  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
  auto t_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachTensorAtPts);
  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
 
  auto t_diff_u =
      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
  auto &nbUniq = dataAtPts->nbUniq;
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
    auto t0_diff = get_diff_rotation_form_vector(t_omega);
    auto t0 = get_rotation_form_vector(t_omega);
 
    t_diff_R(i, j, k) = t0_diff(i, j, k);
    t_R(i, j) = t0(i, j);
 
    FTensor::Tensor1<double, 3> eig;
    FTensor::Tensor2<double, 3, 3> eigen_vec;
    CHKERR get_eigen_val_and_proj_lapack(t_log_u, eig, eigen_vec);
 
    t_eigen_vals(i) = eig(i);
    t_eigen_vecs(i, j) = eigen_vec(i, j);
 
    // rare case when two eigen values are equal
    nbUniq[gg] = get_uniq_nb(eig);
    if (nbUniq[gg] == 2)
      sort_eigen_vals(eig, eigen_vec);
 
    auto t_u_tmp = EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, f);
    auto t_diff_u_tmp =
        EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs, f, d_f, nbUniq[gg]);
    for (int ii = 0; ii != 3; ++ii) {
      for (int jj = ii; jj != 3; ++jj) {
        for (int kk = 0; kk != 3; ++kk) {
          for (int ll = 0; ll < kk; ++ll) {
            t_diff_u_tmp(ii, jj, kk, ll) *= 2;
          }
        }
      }
    }
 
    t_u(i, j) = t_u_tmp(i, j);
    t_diff_u(i, j, k, l) = t_diff_u_tmp(i, j, k, l);
    t_h(i, j) = t_R(i, k) * t_u(k, j);
 
    ++t_h;
    ++t_R;
    ++t_diff_R;
    ++t_log_u;
    ++t_u;
    ++t_diff_u;
    ++t_eigen_vals;
    ++t_eigen_vecs;
    ++t_omega;
  }
 
  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->wDotAtPts);
  if (dataAtPts->wDotDotAtPts.size1() == 0 &&
      dataAtPts->wDotDotAtPts.size2() != nb_integration_pts) {
    dataAtPts->wDotDotAtPts.resize(3, nb_integration_pts);
    dataAtPts->wDotDotAtPts.clear();
  }
  auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wDotDotAtPts);
 
  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 = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
 
  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;
  FTensor::Index<'l', 3> l;
  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);
    FTensor::Tensor2<double, 3, 3> t_PRT;
    t_PRT(i, m) = t_approx_P(i, j) * t_R(m, j);
    FTensor::Tensor1<double, 3> t_leviPRT;
    t_leviPRT(k) = levi_civita(i, m, k) * t_PRT(i, m);
    int bb = 0;
    for (; bb != nb_dofs / 3; ++bb) {
      t_nf(k) += a * t_row_base_fun * t_leviPRT(k);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
    ++t_w;
    ++t_approx_P;
    ++t_R;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical::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_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_dot_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachDotTensorAtPts);
  auto t_diff_u =
      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
 
  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::Tensor2_symmetric<FTensor::PackPtr<double *, 6>, 3>(
        &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<double, 3, 3> t_RTP;
    t_RTP(i, j) = t_R(k, i) * t_approx_P(k, j);
    FTensor::Tensor2<double, 3, 3> t_deltaP;
    t_deltaP(i, j) = t_RTP(i, j) - t_P(i, j);
 
    FTensor::Tensor2_symmetric<double, 3> t_viscous_stress;
    t_viscous_stress(i, j) = alphaU * t_dot_log_u(i, j);
 
    FTensor::Tensor2_symmetric<double, 3> t1;
    t1(i, j) =
        a * (t_diff_u(k, l, i, j) * (t_deltaP(k, l) - t_viscous_stress(k, l)));
 
    int bb = 0;
    for (; bb != nb_dofs / 6; ++bb) {
 
      t_nf(i, j) += t_row_base_fun * t1(i, j);
 
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_approx_P;
    ++t_P;
    ++t_R;
    ++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_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
  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_R(i, k) * t_u(k, j) - t_kd(i, j);
    FTensor::Tensor2<double, 3, 3> t_residuum_omega;
    t_residuum_omega(i, j) =
        (levi_civita(i, m, k) * t_omega_dot(k)) * t_R(m, 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(i) += a * t_row_base_fun(j) * t_residuum_omega(i, j);
 
      ++t_nf;
      ++t_row_base_fun;
    }
 
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_u;
    ++t_R;
    ++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_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
  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_R(i, k) * t_u(k, j) - t_kd(i, j);
    FTensor::Tensor2<double, 3, 3> t_residuum_omega;
    t_residuum_omega(i, j) =
        (levi_civita(i, m, k) * t_omega_dot(k)) * t_R(m, j);
 
    int bb = 0;
    for (; bb != nb_dofs; ++bb) {
      t_nf += a * t_row_base_fun(i, j) * t_residuum(i, j);
      t_nf += a * t_row_base_fun(i, j) * t_residuum_omega(i, j);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb) {
      ++t_row_base_fun;
    }
    ++t_w;
    ++t_R;
    ++t_u;
    ++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_s_w = getFTensor1FromMat<3>(dataAtPts->wAtPts);
  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_s_w(i);
      ++t_nf;
      ++t_row_diff_base_fun;
    }
    for (; bb != nb_base_functions; ++bb) {
      ++t_row_diff_base_fun;
    }
    ++t_w;
    ++t_s_w;
  }
 
  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]);
      };
 
      // get bc data
      FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
      t_bc_disp(i) *= getFEMethod()->ts_t;
 
      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);
      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 (bc.blockName.compare(20, 4, "COOK") == 0)
          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;
  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::Dg<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),
 
        &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<'l', 3> l;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
 
  auto t_diff_u =
      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
 
  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) {
 
        FTensor::Tensor3<double, 3, 3, 3> t_RTP_dP;
        t_RTP_dP(i, j, k) = t_R(k, i) * t_col_base_fun(j) * t_row_base_fun;
 
        for (int kk = 0; kk != 3; ++kk)
          for (int ii = 0; ii != 3; ++ii)
            for (int jj = ii; jj != 3; ++jj) {
              for (int mm = 0; mm != 3; ++mm)
                for (int nn = 0; nn != 3; ++nn)
                  t_m(ii, jj, kk) +=
                      a * t_diff_u(mm, nn, ii, jj) * t_RTP_dP(mm, nn, kk);
            }
 
        ++t_col_base_fun;
        ++t_m;
      }
 
      ++t_row_base_fun;
    }
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base_fun;
    ++t_w;
    ++t_R;
    ++t_diff_u;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical_du_dBubble::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_symmetric<FTensor::PackPtr<double *, 1>, 3>(
        &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;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_diff_u =
      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
 
  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) {
 
        FTensor::Tensor2<double, 3, 3> t_RTP_dBubble;
        t_RTP_dBubble(i, j) =
            a * t_row_base_fun * t_R(k, i) * t_col_base_fun(k, j);
        t_m(i, j) += t_diff_u(m, n, i, j) * t_RTP_dBubble(m, n);
 
        ++t_m;
        ++t_col_base_fun;
      }
 
      ++t_row_base_fun;
    }
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base_fun;
    ++t_w;
    ++t_R;
    ++t_diff_u;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical_du_du::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_ftensor4 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Ddg<FTensor::PackPtr<double *, 6>, 3, 3>(
 
        &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;
 
  FTensor::Number<0> N0;
  FTensor::Number<1> N1;
  FTensor::Number<2> N2;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_P_dh0 = getFTensor3FromMat(dataAtPts->P_dh0);
  auto t_P_dh1 = getFTensor3FromMat(dataAtPts->P_dh1);
  auto t_P_dh2 = getFTensor3FromMat(dataAtPts->P_dh2);
  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
  auto t_dot_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachDotTensorAtPts);
  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
  auto t_diff_u =
      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
  auto &nbUniq = dataAtPts->nbUniq;
 
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
  FTensor::Tensor4<double, 3, 3, 3, 3> t_one4;
  t_one4(i, j, k, l) = t_kd(j, l) * t_kd(i, k);
 
  FTensor::Tensor4<double, 3, 3, 3, 3> t_one4_symmetric;
  t_one4_symmetric(i, j, k, l) = 0;
  for (auto ii : {0, 1, 2})
    for (auto jj : {0, 1, 2})
      for (auto kk : {0, 1, 2})
        for (auto ll : {0, 1, 2}) {
 
          if (ll != kk)
            t_one4_symmetric(ii, jj, kk, ll) =
                t_one4(ii, jj, kk, ll) + t_one4(ii, jj, ll, kk);
          else
            t_one4_symmetric(ii, jj, kk, ll) = t_one4(ii, jj, kk, ll);
        }
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  FTensor::Tensor2_symmetric<double, 3> t_one;
  t_one(i, j) = 0;
  for (auto ii : {0, 1, 2})
    t_one(ii, ii) = 1;
 
  const double ts_a = getTSa();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    FTensor::Tensor4<double, 3, 3, 3, 3> t_P_dh;
    t_P_dh(i, j, N0, k) = t_P_dh0(i, j, k);
    t_P_dh(i, j, N1, k) = t_P_dh1(i, j, k);
    t_P_dh(i, j, N2, k) = t_P_dh2(i, j, k);
 
    FTensor::Tensor2<double, 3, 3> t_RTP;
    t_RTP(i, j) = t_R(k, i) * t_approx_P(k, j);
    FTensor::Tensor2<double, 3, 3> t_deltaP;
    t_deltaP(i, j) = t_RTP(i, j) - t_P(i, j);
 
    FTensor::Tensor2<double, 3, 3> t_dot_u;
    t_dot_u(i, j) = t_u(i, k) * t_dot_log_u(k, j);
 
    FTensor::Tensor4<double, 3, 3, 3, 3> t_stress_diff;
    t_stress_diff(i, j, k, l) = 0;
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj)
        for (int kk = 0; kk != 3; ++kk)
          for (int ll = 0; ll != 3; ++ll)
            for (int mm = 0; mm != 3; ++mm)
              for (int nn = 0; nn != 3; ++nn)
                t_stress_diff(ii, jj, mm, nn) -=
                    t_P_dh(ii, jj, kk, ll) * t_diff_u(kk, ll, mm, nn);
 
    FTensor::Tensor2_symmetric<double, 3> t_viscous_stress;
    t_viscous_stress(i, j) = alphaU * t_dot_log_u(i, j);
    FTensor::Tensor4<double, 3, 3, 3, 3> t_viscous_stress_diff;
    t_viscous_stress_diff(i, j, k, l) = t_one4_symmetric(i, j, k, l);
    t_viscous_stress_diff(i, j, k, l) *= (ts_a * alphaU);
 
    FTensor::Tensor2<double, 3, 3> t_deltaP2;
    t_deltaP2(i, j) = t_deltaP(i, j) - t_viscous_stress(i, j);
 
    auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
        t_eigen_vals, t_eigen_vecs, f, d_f, dd_f, t_deltaP2, nbUniq[gg]);
    for (int ii = 0; ii != 3; ++ii) {
      for (int jj = 0; jj < ii; ++jj) {
        for (int kk = 0; kk != 3; ++kk) {
          for (int ll = 0; ll != kk; ++ll) {
            t_diff2_uP2(ii, jj, kk, ll) *= 2;
          }
        }
      }
    }
    for (int ii = 0; ii != 3; ++ii) {
      for (int jj = ii; jj != 3; ++jj) {
        for (int kk = 0; kk != 3; ++kk) {
          for (int ll = 0; ll < kk; ++ll) {
            t_diff2_uP2(ii, jj, kk, ll) *= 2;
          }
        }
      }
    }
 
    int rr = 0;
    for (; rr != row_nb_dofs / 6; ++rr) {
 
      auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
      auto t_m = get_ftensor4(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;
 
        for (int ii = 0; ii != 3; ++ii)
          for (int jj = ii; jj != 3; ++jj)
            for (int kk = 0; kk != 3; ++kk)
              for (int ll = kk; ll != 3; ++ll)
                for (int mm = 0; mm != 3; ++mm)
                  for (int nn = 0; nn != 3; ++nn)
                    t_m(ii, jj, kk, ll) +=
                        b * t_diff_u(mm, nn, ii, jj) *
                        (t_stress_diff(mm, nn, kk, ll) -
                         t_viscous_stress_diff(mm, nn, kk, ll));
 
        t_m(i, j, k, l) += b * t_diff2_uP2(i, j, 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_P_dh0;
    ++t_P_dh1;
    ++t_P_dh2;
    ++t_R;
    ++t_approx_P;
    ++t_P;
    ++t_dot_log_u;
    ++t_u;
    ++t_diff_u;
    ++t_eigen_vals;
    ++t_eigen_vecs;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical_du_domega::integrate(EntData &row_data,
                                                      EntData &col_data) {
  MoFEMFunctionBegin;
  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::Dg<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),
 
        &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;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
  auto t_diff_u =
      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStreachTensorAtPts);
 
  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;
 
    FTensor::Tensor3<double, 3, 3, 3> t_RT_domega_P;
    t_RT_domega_P(i, j, m) = t_diff_R(k, i, m) * t_P(k, j);
    FTensor::Dg<double, 3, 3> t;
    t(i, j, k) = 0;
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = ii; jj != 3; ++jj)
        for (int mm = 0; mm != 3; ++mm)
          for (int nn = 0; nn != 3; ++nn)
            for (int kk = 0; kk != 3; ++kk)
              t(ii, jj, kk) +=
                  t_diff_u(mm, nn, ii, jj) * t_RT_domega_P(mm, nn, kk);
 
    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(i, j, k) += v * t(i, j, k);
 
        ++t_m;
        ++t_col_base_fun;
      }
 
      ++t_row_base_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_P;
    ++t_diff_R;
    ++t_diff_u;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialRotation_domega_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_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_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
 
  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_base_fun = col_data.getFTensor1N<3>(gg, 0);
      auto t_m = get_ftensor2(K, 3 * rr, 0);
      FTensor::Tensor3<double, 3, 3, 3> t_levi_R;
      t_levi_R(i, j, k) =
          (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        FTensor::Tensor1<double, 3> t_pushed_base;
        t_m(k, i) += t_levi_R(i, j, k) * t_col_base_fun(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_R;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialRotation_domega_dBubble::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 *, 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_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
 
  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_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
      auto t_m = get_ftensor1(K, 3 * rr, 0);
 
      FTensor::Tensor3<double, 3, 3, 3> t_levi_r;
      t_levi_r(i, j, k) =
          (a * t_row_base_fun) * (levi_civita(i, m, k) * t_R(m, j));
 
      for (int cc = 0; cc != col_nb_dofs; ++cc) {
        FTensor::Tensor2<double, 3, 3> t_pushed_base;
        t_m(k) += t_levi_r(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_R;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialRotation_domega_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_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
 
  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;
 
    FTensor::Tensor2<double, 3, 3> t_PRT;
    t_PRT(i, m) = t_approx_P(i, j) * t_R(m, j);
 
    FTensor::Tensor2<double, 3, 3> t_leviPRT_domega;
    t_leviPRT_domega(n, l) =
        levi_civita(i, j, n) * (t_approx_P(i, k) * t_diff_R(j, k, l));
 
    int rr = 0;
    for (; rr != row_nb_dofs / 3; ++rr) {
 
      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) {
 
        const double b = a * t_row_base_fun * t_col_base_fun;
 
        t_m(k, m) += b * t_leviPRT_domega(k, m);
 
        ++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;
    ++t_R;
    ++t_diff_R;
  }
  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_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
  auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
  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_diff_R(i, l, k) * t_u(l, j));
 
      FTensor::Tensor2<double, 3, 3> t_levi1, t_levi2;
      t_levi1(i, k) = (levi_civita(i, m, n) * t_omega_dot(n)) *
                      (t_row_base_fun(j) * t_diff_R(m, j, k));
      t_levi2(i, k) = levi_civita(i, m, k) * (t_row_base_fun(j) * t_R(m, j));
 
      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(i, j) += (a * t_col_base_fun) * t_levi1(i, j);
        t_m(i, j) += ((a * ts_a) * t_col_base_fun) * t_levi2(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_R;
    ++t_diff_R;
    ++t_u;
    ++t_omega_dot;
  }
  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_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_diff_R = getFTensor3FromMat(dataAtPts->diffRotMatAtPts);
  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
  auto t_omega_dot = getFTensor1FromMat<3>(dataAtPts->rotAxisDotAtPts);
  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_diff_R(i, l, k) * t_u(l, j));
 
      FTensor::Tensor2<double, 3, 3> t_levi_omega;
      t_levi_omega(i, m) = levi_civita(i, m, n) * t_omega_dot(n);
      FTensor::Tensor3<double, 3, 3, 3> t_levi_omegaDiffR;
      t_levi_omegaDiffR(i, j, k) = t_levi_omega(i, m) * t_diff_R(m, j, k);
      FTensor::Tensor3<double, 3, 3, 3> t_levi_omegaR;
      t_levi_omegaR(i, j, k) = levi_civita(i, m, k) * t_R(m, j);
 
      FTensor::Tensor1<double, 3> t_levi1, t_levi2;
      t_levi1(k) = t_row_base_fun(i, j) * t_levi_omegaDiffR(i, j, k);
      t_levi2(k) = t_row_base_fun(i, j) * t_levi_omegaR(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(j) += (a * t_col_base_fun) * t_levi1(j);
        t_m(j) += ((a * ts_a) * t_col_base_fun) * t_levi2(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_R;
    ++t_diff_R;
    ++t_u;
    ++t_omega_dot;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpPostProcDataStructure::doWork(int side, EntityType type,
                                               EntData &data) {
  MoFEMFunctionBegin;
  if (type != MBVERTEX)
    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_sigma = create_tag("CauchyStress", 9);
  Tag th_log_u = create_tag("LogSpatialStreach", 9);
  Tag th_u = create_tag("SpatialStreach", 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("SpatialStreachEigenVec", 9);
  Tag th_u_eig_vals = create_tag("SpatialStreachEigenVals", 3);
 
  int nb_gauss_pts = data.getN().size1();
  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Nb. of integration points is not equal to number points on "
            "post-processing mesh");
  }
 
  auto t_w = getFTensor1FromMat<3>(dataAtPts->wAtPts);
  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
  auto t_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStreachTensorAtPts);
  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
  auto t_coords = getFTensor1CoordsAtGaussPts();
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  // 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);
    CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
                                     &*v.data().begin());
    MoFEMFunctionReturn(0);
  };
 
  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);
    CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
                                     &*m.data().begin());
    MoFEMFunctionReturn(0);
  };
 
  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
 
    // vetors
    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);
 
    // 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);
 
    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::Tensor2<double, 3, 3> t_PhT;
    t_PhT(i, k) = t_approx_P(i, j) * t_R(k, j);
    FTensor::Tensor1<double, 3> t_leviPRT;
    t_leviPRT(k) = levi_civita(i, l, k) * t_PhT(i, l);
 
    CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
                                     &t_leviPRT(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);
 
      // LAPACK - eigenvalues and vectors. Applied twice for initial creates
      // memory space
      int n = 3, lda = 3, info, lwork = -1;
      double wkopt;
      info = lapack_dsyev('V', 'U', n, &(m.data()[0]), lda,
                          &(eigen_values.data()[0]), &wkopt, lwork);
      if (info != 0)
        SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "is something wrong with lapack_dsyev info = %d", info);
      lwork = (int)wkopt;
      double work[lwork];
      info = lapack_dsyev('V', 'U', n, &(m.data()[0]), lda,
                          &(eigen_values.data()[0]), work, lwork);
      if (info != 0)
        SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "is something wrong with lapack_dsyev info = %d", info);
 
      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_coords;
  }
 
  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;
      (*energy) += a * t_P(i, J) * t_h(i, J);
      ++t_w;
      ++t_P;
      ++t_h;
    }
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity