EshelbianOperators.cpp

Implementation of operators.

Implementation of operators

/**
 * \file EshelbianOperators.cpp
 * \example EshelbianOperators.cpp
 *
 * \brief Implementation of operators
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
#include <boost/math/constants/constants.hpp>
 
#include <lapack_wrap.h>
 
#include <MatrixFunction.hpp>
 
namespace EshelbianPlasticity {
 
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);
 
  const auto angle =
      sqrt(t_omega(i) * t_omega(i)) + std::numeric_limits<double>::epsilon();
  if (rotSelector == SMALL_ROT) {
    t_R(i, j) += t_Omega(i, j);
    return t_R;
  }
 
  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.);
  const auto b = 2. * ss_2 * ss_2 / (angle * angle);
  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->hdstretchAtPts.resize(9 * 6, nb_integration_pts, false);
  dataAtPts->leviKirchoffAtPts.resize(3, nb_integration_pts, false);
  dataAtPts->leviKirchoffdPAtPts.resize(9 * 3, nb_integration_pts, false);
  dataAtPts->leviKirchoffdOmegaAtPts.resize(9, nb_integration_pts, false);
 
  dataAtPts->adjontPdstretchAtPts.resize(9, nb_integration_pts, false);
  dataAtPts->adjontPdUAtPts.resize(size_symm, nb_integration_pts, false);
  dataAtPts->adjontPdUdPAtPts.resize(9 * size_symm, nb_integration_pts, false);
  dataAtPts->adjontPdUdOmegaAtPts.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);
 
  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
  auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
  auto t_h_du = getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdstretchAtPts);
  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchoffAtPts);
  auto t_levi_kirchoff_dP =
      getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchoffdPAtPts);
  auto t_levi_kirchoff_domega =
      getFTensor2FromMat<3, 3>(dataAtPts->leviKirchoffdOmegaAtPts);
  auto t_approx_P_adjont_dstretch =
      getFTensor2FromMat<3, 3>(dataAtPts->adjontPdstretchAtPts);
  auto t_approx_P_adjont_du =
      getFTensor1FromMat<size_symm>(dataAtPts->adjontPdUAtPts);
  auto t_approx_P_adjont_du_dP =
      getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjontPdUdPAtPts);
  auto t_approx_P_adjont_du_domega =
      getFTensor2FromMat<3, size_symm>(dataAtPts->adjontPdUdOmegaAtPts);
  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_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::Tensor2<double, 3, 3> t_h1;
    t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
    FTensor::Tensor2<double, 3, 3> t_C1;
    t_C1(i, j) = t_h1(k, i) * t_h1(k, j);
 
    eigen_vec(i, j) = t_C1(i, j);
    CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
    auto nb_uniq = get_uniq_nb<3>(&eig(0));
    if (nb_uniq < 3) {
      sort_eigen_vals<3>(eig, eigen_vec);
    }
 
    auto t_inv_u_h1 = EigenMatrix::getMat(eig, eigen_vec, [](const double v) {
      return 1. / (std::sqrt(v) + std::numeric_limits<double>::epsilon());
    });
 
    FTensor::Tensor2<double, 3, 3> t_R_h1;
    t_R_h1(i, j) = t_h1(i, k) * t_inv_u_h1(k, j);
 
    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);
 
    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);
    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);
 
    FTensor::Tensor2<double, 3, 3> t_RR;
    FTensor::Tensor3<double, 3, 3, 3> t_diff_RR;
 
    switch (EshelbianCore::gradApperoximator) {
    case LARGE_ROT:
 
      t_RR(l, j) = t_R_h1(l, i) * t_R(i, j);
      t_diff_RR(l, j, k) = t_R_h1(l, i) * t_diff_R(i, j, k);
 
      t_h(i, j) = t_RR(i, k) * t_u(k, j);
      t_h_domega(i, j, k) = t_diff_RR(i, l, k) * t_u(l, j);
      t_h_du(i, j, L) = t_RR(i, k) * t_Ldiff_u(k, j, L);
 
      t_levi_kirchoff(k) =
          levi_civita(i, j, k) * (t_approx_P(k, j) * t_RR(k, i));
      t_levi_kirchoff_dP(l, j, k) = levi_civita(i, j, k) * t_RR(l, i);
      t_levi_kirchoff_domega(k, l) =
          levi_civita(i, j, k) * (t_approx_P(m, j) * t_diff_RR(m, i, l));
 
      t_approx_P_adjont_dstretch(i, j) = t_approx_P(k, j) * t_RR(k, i);
      t_approx_P_adjont_du(L) =
          t_Ldiff_u(i, j, L) * t_approx_P_adjont_dstretch(i, j);
      t_approx_P_adjont_du_dP(i, j, L) = t_h_du(i, j, L);
      t_approx_P_adjont_du_domega(m, L) =
          t_Ldiff_u(i, j, L) * (t_diff_RR(k, i, m) * t_approx_P(k, j));
 
      break;
    case MODERATE_ROT:
 
      t_h(l, j) = t_R_h1(l, i) * (t_kd(i, j) + (t_R(i, j) - t_kd(i, j)) +
                                  (t_u(i, j) - t_kd(i, j)));
      t_h_domega(l, j, k) = t_R_h1(l, i) * t_diff_R(i, j, k);
      t_h_du(l, j, L) = t_R_h1(l, i) *  t_Ldiff_u(i, j, L);
 
      t_levi_kirchoff(k) =
          levi_civita(i, j, k) * (t_approx_P(l, j) * t_R_h1(l, i));
      t_levi_kirchoff_dP(l, j, k) = levi_civita(i, j, k) * t_R_h1(l, i);
      t_levi_kirchoff_domega(k, l) = 0;
 
      t_approx_P_adjont_dstretch(i, j) = t_approx_P(k, j) * t_R_h1(k, i);
      t_approx_P_adjont_du(L) =
          t_Ldiff_u(i, j, L) * t_approx_P_adjont_dstretch(i, j);
      t_approx_P_adjont_du_dP(i, j, L) = t_h_du(i, j, L);
      t_approx_P_adjont_du_domega(m, L) = 0;
 
      break;
 
    case SMALL_ROT:
 
      t_h(i, j) =
          t_kd(i, j) + (t_R(i, j) - t_kd(i, j)) + (t_u(i, j) - t_kd(i, j));
      t_h_domega(i, j, k) = t_diff_R(i, j, k);
      t_h_du(i, j, L) = t_Ldiff_u(i, j, L);
 
      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_du(L) =
          t_Ldiff_u(i, j, L) * t_approx_P_adjont_dstretch(i, j);
      t_approx_P_adjont_du_dP(i, j, L) = t_h_du(i, j, L);
      t_approx_P_adjont_du_domega(m, L) = 0;
 
      break;
    }
 
    ++t_h;
    ++t_h_domega;
    ++t_h_du;
    ++t_levi_kirchoff;
    ++t_levi_kirchoff_dP;
    ++t_levi_kirchoff_domega;
    ++t_approx_P_adjont_dstretch;
    ++t_approx_P_adjont_du;
    ++t_approx_P_adjont_du_dP;
    ++t_approx_P_adjont_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;
  }
 
  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->leviKirchoffAtPts);
  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 {
    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_du =
      getFTensor1FromMat<size_symm>(dataAtPts->adjontPdUAtPts);
  auto t_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
  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_u(k, l) + alphaU * t_dot_log_u(k, l));
    FTensor::Tensor1<double, size_symm> t_residual;
    t_residual(L) = a * (t_approx_P_adjont_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_du;
    ++t_log_u;
    ++t_dot_log_u;
  }
  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_du =
      getFTensor1FromMat<size_symm>(dataAtPts->adjontPdUAtPts);
  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_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_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_s_w = 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_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, EshelbianCore::rotSelector);
      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_du_dP =
      getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjontPdUdPAtPts);
  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_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_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_du_dP =
      getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjontPdUdPAtPts);
 
  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_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_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 {
    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;
 
  FTensor::Number<0> N0;
  FTensor::Number<1> N1;
  FTensor::Number<2> N2;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
 
  auto t_approx_P_adjont_dstretch =
      getFTensor2FromMat<3, 3>(dataAtPts->adjontPdstretchAtPts);
  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;
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
 
  const double ts_a = getTSa();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    FTensor::Tensor2<double, size_symm, size_symm> t_dP;
 
    if (EshelbianCore::stretchSelector > LINEAR) {
      // 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.;
    } else {
      t_dP(L, J) = 0;
    }
 
    t_dP(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)));
 
    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_approx_P_adjont_dstretch;
    ++t_dot_log_u;
    ++t_u;
    ++t_diff_u;
    ++t_eigen_vals;
    ++t_eigen_vecs;
  }
  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;
 
  FTensor::Number<0> N0;
  FTensor::Number<1> N1;
  FTensor::Number<2> N2;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_P_dh0 = getFTensor3FromMat<3, 3, 3>(dataAtPts->P_dh0);
  auto t_P_dh1 = getFTensor3FromMat<3, 3, 3>(dataAtPts->P_dh1);
  auto t_P_dh2 = getFTensor3FromMat<3, 3, 3>(dataAtPts->P_dh2);
  auto t_approx_P_adjont_dstretch =
      getFTensor2FromMat<3, 3>(dataAtPts->adjontPdstretchAtPts);
  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
  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;
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  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_deltaP;
    t_deltaP(i, j) = t_approx_P_adjont_dstretch(i, j) - t_P(i, j);
 
    // Work of symmetric tensor on undefined tensor is equal to the work of the
    // symmetric part of it
    FTensor::Tensor2<double, size_symm, size_symm> t_dP;
 
    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));
    } else {
      t_dP(L, J) = 0;
    }
 
    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) * (t_P_dh(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)));
 
    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_P_dh0;
    ++t_P_dh1;
    ++t_P_dh2;
    ++t_P;
    ++t_approx_P_adjont_dstretch;
    ++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;
 
  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_du_domega =
      getFTensor2FromMat<3, size_symm>(dataAtPts->adjontPdUdOmegaAtPts);
 
  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_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_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->leviKirchoffdPAtPts);
 
  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->leviKirchoffdPAtPts);
 
  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->leviKirchoffdOmegaAtPts);
  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("U", 3);
  Tag th_disp_error = create_tag("U_ERROR", 1);
 
  auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
  auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
  auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
  auto t_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->StretchTensorAtPts);
  auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
  auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
  auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchoffAtPts);
  auto t_coords = getFTensor1CoordsAtGaussPts();
  auto t_normal = getFTensor1NormalsAtGaussPts();
  auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
 
  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);
 
    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;
  }
 
  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