ElasticOperators.hpp

Go to the documentation of this file.

/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
namespace OpElasticTools {
 
// constexpr auto VecSetValues = MoFEM::VecSetValues<EssentialBcStorage>;
// constexpr auto MatSetValues = MoFEM::MatSetValues<EssentialBcStorage>;
// constexpr auto MatSetValues = MoFEM::MatSetValues<MpEntityStorage>;
 
using OpMassPrecHTensorBound = FormsIntegrators<BoundaryEleOp>::Assembly<
    USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<3, 9>;
 
using OpMassPrecHTensorVol = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<3, 9>;
 
using OpMassPrecScalar = FormsIntegrators<DomainEleOp>::Assembly<
    USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<1, 1>;
 
//! [Common data]
struct CommonData {
  // Ddg<double, 3, 3> tD;
  boost::shared_ptr<MatrixDouble> mtD;
  boost::shared_ptr<MatrixDouble> mtD_Axiator;
  boost::shared_ptr<MatrixDouble> mtD_Deviator;
 
  boost::shared_ptr<MatrixDouble> mGradPtr;
  boost::shared_ptr<MatrixDouble> mStrainPtr;
  boost::shared_ptr<MatrixDouble> mStressPtr;
  boost::shared_ptr<MatrixDouble> mPiolaStressPtr;
 
  boost::shared_ptr<MatrixDouble> mDeformationGradient;
  boost::shared_ptr<MatrixDouble> matC;
  boost::shared_ptr<VectorDouble> detF;
 
  boost::shared_ptr<MatrixDouble> matEigenVal;
  boost::shared_ptr<MatrixDouble> matEigenVector;
 
  boost::shared_ptr<MatrixDouble> materialTangent;
  boost::shared_ptr<MatrixDouble> dE_dF_mat;
 
  array<int, 2> stateArrayPlast;
  array<int, 2> stateArrayCont;
  array<double, 2> stateArrayArea;
  vector<double> bodyReactions;
 
  SmartPetscObj<Vec> reactionsVec;
  Range reactionEnts;
 
  EntData *faceRowData;
  double mMeasure;
  double time;
 
  VectorInt rowIndices;
 
  boost::shared_ptr<MatrixDouble> mDispPtr;
  bool isNeoHook;
 
  Tag reactionTag;
 
  CommonData()
      : stateArrayCont({0, 0}), stateArrayPlast({0, 0}),
        stateArrayArea({0., 0.}), isNeoHook(false) {
  }
 
  ~CommonData() {}
 
  double timeStepSize;
 
  SmartPetscObj<Mat> SEpEp;
  SmartPetscObj<AO> aoSEpEp;
  SmartPetscObj<Mat> STauTau;
  SmartPetscObj<AO> aoSTauTau;
 
  boost::shared_ptr<BlockMatContainer> blockMatContainerPtr;
  MatrixDouble massMatTensor;
  MatrixDouble massMatScalar;
};
 
template <typename T> inline auto to_non_symm(T &symm) {
  Tensor2<double, 3, 3> non_symm;
  Number<0> N0;
  Number<1> N1;
  Number<2> N2;
  non_symm(N0, N0) = symm(N0, N0);
  non_symm(N1, N1) = symm(N1, N1);
  non_symm(N2, N2) = symm(N2, N2);
  non_symm(N0, N1) = non_symm(N1, N0) = symm(N0, N1);
  non_symm(N0, N2) = non_symm(N2, N0) = symm(N0, N2);
  non_symm(N2, N1) = non_symm(N1, N2) = symm(N2, N1);
  return non_symm;
};
 
template <typename T> inline auto to_symm(T &non_symm) {
  Tensor2_symmetric<double, 3> symm;
  Number<0> N0;
  Number<1> N1;
  Number<2> N2;
  symm(N0, N0) = non_symm(N0, N0);
  symm(N1, N1) = non_symm(N1, N1);
  symm(N2, N2) = non_symm(N2, N2);
  symm(N0, N1) = non_symm(N0, N1);
  symm(N0, N2) = non_symm(N0, N2);
  symm(N1, N2) = non_symm(N1, N2);
  return symm;
};
 
inline bool is_eq(const double &a, const double &b) {
  constexpr double eps = 1e-12;
  return abs(a - b) < eps;
};
inline bool is_diff(const double &a, const double &b) {
  constexpr double eps = 1e-12;
  return abs(a - b) > eps;
};
 
template <typename T>
inline Ddg<double, 3, 3> tensor4_to_ddg(Tensor4<T, 3, 3, 3, 3> &t) {
  Ddg<double, 3, 3> loc_tens;
 
  Number<0> N0;
  Number<1> N1;
  Number<2> N2;
 
  loc_tens(N0, N0, N0, N0) = t(N0, N0, N0, N0); // 0
  loc_tens(N0, N0, N0, N1) = t(N0, N0, N0, N1); // 1
  loc_tens(N0, N0, N0, N2) = t(N0, N0, N0, N2); // 2
  loc_tens(N0, N0, N1, N1) = t(N0, N0, N1, N1); // 3
  loc_tens(N0, N0, N1, N2) = t(N0, N0, N1, N2); // 4
  loc_tens(N0, N0, N2, N2) = t(N0, N0, N2, N2); // 5
  loc_tens(N0, N1, N0, N0) = t(N0, N1, N0, N0); // 6
  loc_tens(N0, N1, N0, N1) = t(N0, N1, N0, N1); // 7
  loc_tens(N0, N1, N0, N2) = t(N0, N1, N0, N2); // 8
  loc_tens(N0, N1, N1, N1) = t(N0, N1, N1, N1); // 9
  loc_tens(N0, N1, N1, N2) = t(N0, N1, N1, N2); // 10
  loc_tens(N0, N1, N2, N2) = t(N0, N1, N2, N2); // 11
  loc_tens(N0, N2, N0, N0) = t(N0, N2, N0, N0); // 12
  loc_tens(N0, N2, N0, N1) = t(N0, N2, N0, N1); // 13
  loc_tens(N0, N2, N0, N2) = t(N0, N2, N0, N2); // 14
  loc_tens(N0, N2, N1, N1) = t(N0, N2, N1, N1); // 15
  loc_tens(N0, N2, N1, N2) = t(N0, N2, N1, N2); // 16
  loc_tens(N0, N2, N2, N2) = t(N0, N2, N2, N2); // 17
  loc_tens(N1, N1, N0, N0) = t(N1, N1, N0, N0); // 18
  loc_tens(N1, N1, N0, N1) = t(N1, N1, N0, N1); // 19
  loc_tens(N1, N1, N0, N2) = t(N1, N1, N0, N2); // 20
  loc_tens(N1, N1, N1, N1) = t(N1, N1, N1, N1); // 21
  loc_tens(N1, N1, N1, N2) = t(N1, N1, N1, N2); // 22
  loc_tens(N1, N1, N2, N2) = t(N1, N1, N2, N2); // 23
  loc_tens(N1, N2, N0, N0) = t(N1, N2, N0, N0); // 24
  loc_tens(N1, N2, N0, N1) = t(N1, N2, N0, N1); // 25
  loc_tens(N1, N2, N0, N2) = t(N1, N2, N0, N2); // 26
  loc_tens(N1, N2, N1, N1) = t(N1, N2, N1, N1); // 27
  loc_tens(N1, N2, N1, N2) = t(N1, N2, N1, N2); // 28
  loc_tens(N1, N2, N2, N2) = t(N1, N2, N2, N2); // 29
  loc_tens(N2, N2, N0, N0) = t(N2, N2, N0, N0); // 30
  loc_tens(N2, N2, N0, N1) = t(N2, N2, N0, N1); // 31
  loc_tens(N2, N2, N0, N2) = t(N2, N2, N0, N2); // 32
  loc_tens(N2, N2, N1, N1) = t(N2, N2, N1, N1); // 33
  loc_tens(N2, N2, N1, N2) = t(N2, N2, N1, N2); // 34
  loc_tens(N2, N2, N2, N2) = t(N2, N2, N2, N2); // 35
 
  return loc_tens;
};
 
template <typename T> inline auto get_tensor4_from_ddg(Ddg<T, 3, 3> &ddg) {
  Tensor4<double, 3, 3, 3, 3> tens4;
  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)
          tens4(ii, jj, kk, ll) = ddg(ii, jj, kk, ll);
  return tens4;
};
 
template <typename T> inline bool is_ddg(Tensor4<T, 3, 3, 3, 3> tensor) {
  auto test_ddg = tensor4_to_ddg(tensor);
  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, 1, 2}) {
          const double diff =
              abs(test_ddg(ii, jj, kk, ll) - tensor(ii, jj, kk, ll));
          if (diff > 1e-6)
            return false;
        }
 
  return true;
};
 
template <typename T> inline auto exp_map(const Tensor2_symmetric<T, 3> &X) {
  Tensor2<double, 3, 3> exp_tensor(1., 0., 0., 0., 1., 0., 0., 0., 1.);
  Tensor2<double, 3, 3> X_nt;
 
  constexpr int max_it = 1000;
  constexpr double eps = 1e-12;
 
  int fac_n = 1;
 
  auto Xcopy = to_non_symm(X);
  auto X_n = to_non_symm(X);
  exp_tensor(i, j) += Xcopy(i, j);
 
  for (int n = 2; n != max_it; ++n) {
    fac_n *= n;
    X_nt(i, j) = X_n(i, j);
    X_n(i, k) = X_nt(i, j) * Xcopy(j, k);
    const double inv_nf = 1. / fac_n;
    exp_tensor(i, j) += inv_nf * X_n(i, j);
    const double e = inv_nf * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
    if (e < eps)
      return to_symm(exp_tensor);
  }
 
  return to_symm(exp_tensor);
};
 
template <typename T>
inline auto get_rotation_R(Tensor1<T, 3> t1_omega, double tt) {
 
  auto get_rotation = [&](Tensor1<double, 3> &t_omega) {
    FTensor::Tensor2<double, 3, 3> t_R;
 
    constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
    t_R(i, j) = t_kd(i, j);
 
    const double angle = sqrt(t_omega(i) * t_omega(i));
    if (std::abs(angle) < 1e-18)
      return t_R;
 
    FTensor::Tensor2<double, 3, 3> t_Omega;
    t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
    const double a = sin(angle) / angle;
    const double ss_2 = sin(angle / 2.);
    const double 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;
  };
 
  Tensor1<double, 3> my_c_omega;
  my_c_omega(i) = t1_omega(i) * tt;
  auto t_R = get_rotation(my_c_omega);
 
  return get_rotation(my_c_omega);
};
 
template <typename T>
inline auto get_log_map(const Tensor2_symmetric<T, 3> &X) {
 
  constexpr int max_it = 1000;
  constexpr double eps = 1e-12;
 
  auto log_tensor = to_non_symm(X);
  log_tensor(i, j) -= kronecker_delta(i, j);
 
  auto Xcopy = to_non_symm(log_tensor);
  auto X_n = to_non_symm(log_tensor);
  Tensor2<double, 3, 3> X_nt_tmp;
 
  for (int n = 2; n != max_it; ++n) {
 
    X_nt_tmp(i, j) = X_n(i, j);
    X_n(i, k) = X_nt_tmp(i, j) * Xcopy(j, k);
    const double inv_n = 1. / n;
    log_tensor(i, j) += inv_n * pow(-1, n + 1) * X_n(i, j);
    const double e = inv_n * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
    if (e < eps)
      break;
  }
 
  return to_symm(log_tensor);
};
 
template <typename T> inline auto get_diff_log_map(Tensor2_symmetric<T, 3> &X) {
  Tensor2<double, 3, 3> X_nt_tmp;
  Tensor4<double, 3, 3, 3, 3> diff_log_map;
  Tensor4<double, 3, 3, 3, 3> diff_log_map_tmp;
 
  constexpr int max_it = 1000;
  constexpr double eps = 1e-12;
 
  auto get_Xn = [&](auto &log_tensor) {
    vector<Tensor2<double, 3, 3>> Xn;
    Xn.emplace_back();
    auto &cu_Xn = Xn.back();
    cu_Xn(i, j) = kronecker_delta(i, j);
    for (int n = 2; n != max_it; ++n) {
      const double cof = 1. / n;
      auto X_nt_tmp = Xn.back();
      Xn.emplace_back();
      auto &c_Xn = Xn.back();
      c_Xn(i, k) = X_nt_tmp(i, j) * log_tensor(j, k);
      const double e = cof * std::abs(sqrt(c_Xn(i, j) * c_Xn(i, j)));
      if (e < eps)
        return Xn;
    }
    return Xn;
  };
 
  auto log_tensor = to_non_symm(X);
  log_tensor(i, j) -= kronecker_delta(i, j);
  diff_log_map(i, j, k, l) = 0.;
 
  auto Xn_vec = get_Xn(log_tensor);
 
  for (int n = 1; n != Xn_vec.size() + 1; ++n) {
    const double cof = pow(-1, n + 1) / n;
    for (int m = 1; m != n + 1; ++m) {
      auto &Xn_1 = Xn_vec[m - 1];
      auto &Xn_2 = Xn_vec[n - m];
      diff_log_map(i, j, k, l) += cof * Xn_1(i, k) * Xn_2(l, j);
    }
  }
 
  auto isddg = tensor4_to_ddg((*cache).Is);
  diff_log_map_tmp(i, j, k, l) = isddg(i, j, m, n) * diff_log_map(m, n, k, l);
  return tensor4_to_ddg(diff_log_map_tmp);
};
 
template <typename T>
inline MoFEMErrorCode
get_eigen_val_and_proj_lapack(const Tensor2_symmetric<T, 3> &X,
                              Tensor1<T, 3> &lambda,
                              Tensor2<T, 3, 3> &eig_vec) {
 
  MoFEMFunctionBeginHot;
  Tensor1<double, 3> eig;
  auto eigen_vec_lap = to_non_symm(X);
  constexpr int n = 3;
  constexpr int lda = 3;
  constexpr int lwork = 15;
  std::array<double, 15> work;
  if (lapack_dsyev('V', 'U', n, &(eigen_vec_lap(0, 0)), lda, &eig(0),
                   work.data(), lwork) > 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "The algorithm failed to compute eigenvalues.");
 
  lambda(i) = eig(i);
  eig_vec(i, j) = eigen_vec_lap(i, j);
 
  // perturb
  // const double h = 1e-12;
  // if (is_eq(lambda[0], lambda[1]))
  //   lambda[0] *= (1 + h);
  // if (is_eq(lambda[1], lambda[2]))
  //   lambda[2] *= (1 - h);
 
  MoFEMFunctionReturnHot(0);
}
 
template <typename T1, typename T2, typename T3>
inline auto get_ln_X_lapack(const Tensor2_symmetric<T1, 3> &X,
                            Tensor1<T2, 3> &lambda,
                            Tensor2<T3, 3, 3> &eig_vec) {
  Tensor2_symmetric<double, 3> lnX;
  Tensor2_symmetric<double, 3> diag(0., 0., 0., 0., 0., 0.);
  for (int ii = 0; ii != 3; ++ii)
    diag(ii, ii) = log(lambda(ii));
 
  lnX(i, j) = eig_vec(l, i) ^ (diag(l, k) * eig_vec(k, j));
 
  return lnX;
};
 
template <typename T> inline auto get_ln_X(const Tensor2_symmetric<T, 3> &X) {
  return get_log_map(X);
};
 
template <typename T, typename TP, typename T3>
inline auto
calculate_nominal_moduli_TL(Tensor2_symmetric<T, 3> &X, Tensor2<TP, 3, 3> &F,
                            Tensor2_symmetric<TP, 3> &stress,
                            Tensor1<T, 3> &lambda, Tensor2<T, 3, 3> &eig_vec,
                            Tensor4<T3, 3, 3, 3, 3> &TLs3,
                            bool compute_tangent) {
 
  // Tensor4<double, 3, 3, 3, 3> TLs3;
  Tensor4<double, 3, 3, 3, 3> P3;
  Tensor2<double, 3, 3> invF;
 
  double det_f;
  CHKERR determinantTensor3by3(F, det_f);
  CHKERR invertTensor3by3(F, det_f, invF);
 
  // from https://www.sciencedirect.com/science/article/pii/S0045782502004383
 
  // VectorDouble lambda(3);
  VectorDouble e(3);
  VectorDouble f(3);
  VectorDouble d(3);
  MatrixDouble Vi(3, 3);
  MatrixDouble Ksi(3, 3);
  MatrixDouble Zeta(3, 3);
  Vi.clear();
  Ksi.clear();
  Zeta.clear();
  double eta = 0;
 
  Tensor1<double, 3> _N[3];
  Tensor1<double, 3> _n[3];
 
  for (int dd : {0, 1, 2}) {
    for (int cc : {0, 1, 2})
      _N[dd](cc) = eig_vec(dd, cc);
  }
 
  for (int dd : {0, 1, 2})
    _n[dd](i) = F(i, j) * _N[dd](j);
 
  for (int ii = 0; ii != 3; ++ii) {
    e(ii) = 0.5 * log(lambda(ii));
    d(ii) = 1. / lambda(ii);
    f(ii) = -2. / (lambda(ii) * lambda(ii));
  }
 
  const double &lam0 = lambda(0);
  const double &lam1 = lambda(1);
  const double &lam2 = lambda(2);
 
  auto compute_ksi_eta_2eq = [&](int i_, int j_, int k_) {
    Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
    Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
 
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj)
        if (ii != jj) {
          if (jj == k_ || ii == k_) {
            Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
            Ksi(ii, jj) =
                (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
          }
        }
    eta = Ksi(k_, i_);
  };
 
  if (is_diff(lam0, lam1) && is_diff(lam0, lam2) && is_diff(lam1, lam2)) {
    // all different
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj) {
        if (ii != jj) {
 
          Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
          Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
          for (int kk = 0; kk != 3; ++kk) {
 
            if (kk != ii && kk != jj)
              eta += e(ii) / (2. * (lambda(ii) - lambda(jj)) *
                              (lambda(ii) - lambda(kk)));
          }
        }
      }
 
  } else if (is_eq(lam0, lam1) && is_eq(lam0, lam2) && is_eq(lam1, lam2)) {
    // all the same
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj)
        if (ii != jj) {
 
          Vi(ii, jj) = 0.5 * d(ii);
          Ksi(ii, jj) = 0.125 * f(ii);
        }
    eta = 0.125 * f(0);
  } else if (is_eq(lam0, lam1)) {
 
    compute_ksi_eta_2eq(0, 1, 2);
 
  } else if (is_eq(lam0, lam2)) {
 
    compute_ksi_eta_2eq(0, 2, 1);
 
  } else if (is_eq(lam1, lam2)) {
 
    compute_ksi_eta_2eq(1, 2, 0);
  }
 
  P3(i, j, k, l) = 0.;
 
  Tensor2<double, 3, 3> Mii;
  Tensor2<double, 3, 3> Mij;
  Tensor2<double, 3, 3> Mik;
  Tensor2<double, 3, 3> Mjk;
  Tensor2<double, 3, 3> Mjj;
 
  for (int ii = 0; ii != 3; ++ii) {
    Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
    P3(i, j, k, l) += 0.5 * d(ii) * _N[ii](i) * _N[ii](j) * Mii(k, l);
    for (int jj = 0; jj != 3; ++jj)
      if (ii != jj) {
        Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
        P3(i, j, k, l) += Vi(ii, jj) * _N[ii](i) * _N[jj](j) * Mij(k, l);
      }
  }
 
  if (!compute_tangent)
    return P3;
 
  TLs3(i, j, k, l) = 0.;
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      Zeta(ii, jj) = (stress(i, j) * (_N[ii](i) * _N[jj](j)));
 
  for (int ii = 0; ii != 3; ++ii) {
    Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
    TLs3(i, j, k, l) += 0.25 * f(ii) * Zeta(ii, ii) * Mii(i, j) * Mii(k, l);
  }
 
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      if (jj != ii)
        for (int kk = 0; kk != 3; ++kk)
          if (kk != ii && kk != jj) {
            Mik(i, j) = _n[ii](i) * _N[kk](j) + _n[kk](i) * _N[ii](j);
            Mjk(i, j) = _n[jj](i) * _N[kk](j) + _n[kk](i) * _N[jj](j);
            TLs3(i, j, k, l) += 2. * eta * Zeta(ii, jj) * Mik(i, j) * Mjk(k, l);
          }
 
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      if (jj != ii) {
        Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
        Mjj(i, j) = _n[jj](i) * _N[jj](j) + _n[jj](i) * _N[jj](j);
        const double k2 = 2. * Ksi(ii, jj);
        const double zp = k2 * Zeta(ii, jj);
        const double zp2 = k2 * Zeta(jj, jj);
        TLs3(i, j, k, l) += zp * Mij(i, j) * Mjj(k, l) +
                            zp * Mjj(i, j) * Mij(k, l) +
                            zp2 * Mij(i, j) * Mij(k, l);
      }
 
  Tensor2<double, 3, 3> Piola;
  Tensor2<double, 3, 3> invFPiola;
  Tensor4<double, 3, 3, 3, 3> Ttmp3;
  auto nstress = to_non_symm(stress);
  Piola(k, l) = nstress(m, n) * P3(m, n, k, l);
  invFPiola(i, j) = invF(i, k) * Piola(k, j);
  Ttmp3(i, j, k, l) = 0.5 * invFPiola(i, k) * kronecker_delta(j, l) +
                      0.5 * invFPiola(i, l) * kronecker_delta(j, k);
  TLs3(i, j, k, l) += Ttmp3(i, j, k, l);
 
  return P3;
};
 
template <typename T, typename TP, typename T3>
inline auto
test_projection_lukasz(Tensor2_symmetric<T, 3> &X,
                       Tensor2_symmetric<TP, 3> &stress, Tensor1<T, 3> &lambda,
                       Tensor2<T, 3, 3> &eig_vec, Ddg<T3, 3, 3> &TLs3) {
 
  Tensor1<double, 3> eig{lambda(0), lambda(1), lambda(2)};
  Tensor2<double, 3, 3> eigen_vec;
  eigen_vec(i, j) = eig_vec(i, j);
 
  auto get_uniq_nb = [&](auto eig) {
    return distance(&eig(0), unique(&eig(0), &eig(0) + 3, is_eq));
  };
 
  auto sort_eigen_vals2 = [&](auto &eig, auto &eigen_vec) {
    if (is_eq(eig(0), eig(1))) {
      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)};
      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 {
      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)};
      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);
    }
  };
 
  Tensor2_symmetric<double, 3> Ts;
  Ts(i, j) = stress(i, j);
  auto f = [](double v) { return 0.5 * log(v); };
  auto d_f = [](double v) { return 0.5 / v; };
  auto dd_f = [](double v) { return -0.5 / (v * v); };
 
  auto nb_uniq = get_uniq_nb(eig);
  if (nb_uniq == 2)
    sort_eigen_vals2(eig, eigen_vec);
  auto t_d =
 
      EigenMatrix::getDiffMat(eig, eigen_vec, f, d_f, nb_uniq);
 
  auto t_dd =
 
      EigenMatrix::getDiffDiffMat(eig, eigen_vec, f, d_f, dd_f, Ts, nb_uniq);
 
  TLs3(i, j, k, l) = t_dd(i, j, k, l);
  return t_d;
};
 
template <bool TANGENT, typename T, typename TP, typename T3>
inline auto calculate_lagrangian_moduli_TL(Tensor2_symmetric<TP, 3> &stress,
                                           Tensor1<T, 3> &lambda,
                                           Tensor2<T, 3, 3> &eig_vec,
                                           Ddg<T3, 3, 3> &TLs3) {
 
  Ddg<double, 3, 3> P3;
  // from https://www.sciencedirect.com/science/article/pii/S0045782502004383
 
  VectorDouble e(3);
  VectorDouble f(3);
  VectorDouble d(3);
  MatrixDouble Vi(3, 3);
  MatrixDouble Ksi(3, 3);
  MatrixDouble Zeta(3, 3);
  Vi.clear();
  Ksi.clear();
  Zeta.clear();
  double eta = 0;
 
  Tensor1<double, 3> _N[3];
  for (int dd = 0; dd != 3; ++dd) {
    for (int cc = 0; cc != 3; ++cc)
      _N[dd](cc) = eig_vec(dd, cc);
  }
 
  for (int ii = 0; ii != 3; ++ii) {
    e(ii) = 0.5 * log(lambda(ii));
    d(ii) = 1. / lambda(ii);
    f(ii) = -2. / (lambda(ii) * lambda(ii));
  }
 
  const double &lam0 = lambda(0);
  const double &lam1 = lambda(1);
  const double &lam2 = lambda(2);
 
  auto compute_ksi_eta_2eq = [&](int i_, int j_, int k_) {
    Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
    Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
 
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj)
        if (ii != jj) {
          if (jj == k_ || ii == k_) {
            Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
            Ksi(ii, jj) =
                (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
          }
        }
    eta = Ksi(k_, i_);
  };
 
  if (is_diff(lam0, lam1) && is_diff(lam0, lam2) && is_diff(lam1, lam2)) {
    // all different
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj) {
        if (ii != jj) {
 
          Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
          Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
          for (int kk = 0; kk != 3; ++kk) {
 
            if (kk != ii && kk != jj)
              eta += e(ii) / (2. * (lambda(ii) - lambda(jj)) *
                              (lambda(ii) - lambda(kk)));
          }
        }
      }
 
  } else if (is_eq(lam0, lam1) && is_eq(lam0, lam2) && is_eq(lam1, lam2)) {
    // all the same
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj)
        if (ii != jj) {
 
          Vi(ii, jj) = 0.5 * d(ii);
          Ksi(ii, jj) = 0.125 * f(ii);
        }
    eta = 0.125 * f(0);
  } else if (is_eq(lam0, lam1)) {
 
    compute_ksi_eta_2eq(0, 1, 2);
 
  } else if (is_eq(lam0, lam2)) {
 
    compute_ksi_eta_2eq(0, 2, 1);
 
  } else if (is_eq(lam1, lam2)) {
 
    compute_ksi_eta_2eq(1, 2, 0);
  }
 
  P3(i, j, k, l) = 0.;
  Tensor2_symmetric<double, 3> Mii;
  Tensor2_symmetric<double, 3> Mij;
  Tensor2_symmetric<double, 3> Mik;
  Tensor2_symmetric<double, 3> Mjk;
  Tensor2_symmetric<double, 3> Mjj;
 
  for (int ii = 0; ii != 3; ++ii) {
    Mii(i, j) = (_N[ii](i) ^ _N[ii](j)) + (_N[ii](i) ^ _N[ii](j));
    P3(i, j, k, l) += ((0.5 * d(ii) * _N[ii](i)) ^ _N[ii](j)) * Mii(k, l);
    for (int jj = 0; jj != 3; ++jj)
      if (ii != jj) {
        Mij(i, j) = (_N[ii](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[ii](j));
        P3(i, j, k, l) += ((Vi(ii, jj) * _N[ii](i)) ^ _N[jj](j)) * Mij(k, l);
      }
  }
 
  if (!TANGENT)
    return P3;
 
  TLs3(i, j, k, l) = 0.;
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      Zeta(ii, jj) = stress(i, j) * (_N[ii](i) ^ _N[jj](j));
 
  for (int ii = 0; ii != 3; ++ii) {
    Mii(i, j) = (_N[ii](i) ^ _N[ii](j)) + (_N[ii](i) ^ _N[ii](j));
    TLs3(i, j, k, l) += (0.25 * f(ii) * Zeta(ii, ii) * Mii(i, j)) * Mii(k, l);
  }
 
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      if (jj != ii) {
        Mij(i, j) = (_N[ii](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[ii](j));
        Mjj(i, j) = (_N[jj](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[jj](j));
        const double k2 = 2. * Ksi(ii, jj);
        const double zp = k2 * Zeta(ii, jj);
        const double zp2 = k2 * Zeta(jj, jj);
        TLs3(i, j, k, l) += (zp * Mij(i, j)) * Mjj(k, l) +
                            (zp * Mjj(i, j)) * Mij(k, l) +
                            (zp2 * Mij(i, j)) * Mij(k, l);
        for (int kk = 0; kk != 3; ++kk) {
          if (kk != ii && kk != jj) {
            Mik(i, j) = (_N[ii](i) ^ _N[kk](j)) + (_N[kk](i) ^ _N[ii](j));
            Mjk(i, j) = (_N[jj](i) ^ _N[kk](j)) + (_N[kk](i) ^ _N[jj](j));
            TLs3(i, j, k, l) +=
                (2. * eta * Zeta(ii, jj) * Mik(i, j)) * Mjk(k, l);
          }
        }
      }
 
  return P3;
};
 
template <typename T>
inline Tensor2<double, 3, 3>
getFiniteDiff(const int &kk, const int &ll, CommonData &common_data,
              Tensor2<T, 3, 3> &F, Tensor2<double, 3, 3> &t_stress) {
  MoFEMFunctionBeginHot;
  Tensor2<double, 3, 3> out;
 
  Tensor2<double, 3, 3> strainPl;
  Tensor2<double, 3, 3> strainMin;
  strainPl(i, j) = F(i, j);
  strainMin(i, j) = F(i, j);
  const double h = 1e-8;
  strainPl(kk, ll) += h;
  strainMin(kk, ll) -= h;
  Tensor4<PackPtr<double *, 1>, 3, 3, 3, 3> dummy;
 
  auto calculate_stress = [&](auto &F) {
    Tensor2<double, 3, 3> Piola;
    Tensor2_symmetric<double, 3> stress_sym;
    Tensor2_symmetric<double, 3> C;
    Tensor2_symmetric<double, 3> strain;
 
    C(i, j) = F(m, i) ^ F(m, j);
 
    Tensor1<double, 3> lambda;
    Tensor2<double, 3, 3> eig_vec;
    CHKERR get_eigen_val_and_proj_lapack(C, lambda, eig_vec);
 
    auto lnC = get_ln_X_lapack(C, lambda, eig_vec);
    // auto lnC = get_log_map(C);
    strain(i, j) = 0.5 * lnC(i, j);
    auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*common_data.mtD);
    stress_sym(i, j) = t_D(i, j, k, l) * strain(k, l);
    auto stress_tmp = to_non_symm(stress_sym);
 
    // Tensor4<double, 3, 3, 3, 3> D2;
    // constexpr auto t_kd = Kronecker_Delta<int>();
    // Tensor4<double, 3, 3, 3, 3> dC_dF;
    // dC_dF(i, j, k, l) = (t_kd(i, l) * F(k, j)) + (t_kd(j, l) * F(k, i));
    // auto t_L = get_diff_log_map(C);
    // t_L(i, j, k, l) *= 0.5;
    // D2(i, j, k, l) = t_L(i, j, m, n) * dC_dF(m, n, k, l);
 
    auto D2 = calculate_nominal_moduli_TL(C, F, stress_sym, lambda, eig_vec,
                                          dummy, false);
 
    Piola(k, l) = stress_tmp(i, j) * D2(i, j, k, l);
 
    return Piola;
  };
 
  auto StressPlusH = calculate_stress(strainPl);
  auto StressMinusH = calculate_stress(strainMin);
  // out(i, j) = (StressPlusH(i, j) - t_stress(i, j)) / (h);
 
  out(i, j) = (StressPlusH(i, j) - StressMinusH(i, j)) / (2 * h);
 
  return out;
}
 
template <typename T>
inline Tensor4<double, 3, 3, 3, 3>
getDTensorFunction(CommonData &common_data, Tensor2<T, 3, 3> &F,
                   Tensor2<double, 3, 3> &t_stress) {
  Tensor4<double, 3, 3, 3, 3> my_D;
  for (int k = 0; k != 3; ++k) {
    for (int l = 0; l != 3; ++l) {
      auto d_stress = getFiniteDiff(k, l, common_data, F, t_stress);
      for (int i = 0; i != 3; ++i) {
        for (int j = 0; j != 3; ++j) {
          my_D(i, j, k, l) = d_stress(i, j);
        }
      }
    }
  }
  return my_D;
};
 
//! [Operators definitions]
typedef boost::function<Tensor1<double, 3>(const double, const double,
                                           const double)>
    VectorFun;
 
struct OpStrain : public DomainEleOp {
  OpStrain(const std::string field_name,
           boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpLogStrain : public DomainEleOp {
  OpLogStrain(const std::string field_name,
              boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpStress : public DomainEleOp {
  OpStress(const std::string field_name,
           boost::shared_ptr<CommonData> common_data_ptr, boost::shared_ptr<MatrixDouble> mat_D_Ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  boost::shared_ptr<MatrixDouble> matDPtr;
};
 
template <bool TANGENT>
struct OpCalculateLogStressTangent : public DomainEleOp {
  OpCalculateLogStressTangent(const std::string field_name,
                              boost::shared_ptr<CommonData> common_data_ptr,
                              boost::shared_ptr<MatrixDouble> mat_D_ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  boost::shared_ptr<MatrixDouble> matDPtr;
};
 
struct OpCalculateNeoHookeStressTangent : public DomainEleOp {
  OpCalculateNeoHookeStressTangent(
      const std::string field_name,
      boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
template <bool LOGSTRAIN, bool REACTIONS>
struct OpInternalForceRhs : public DomainEleOpAssembly {
  OpInternalForceRhs(const std::string field_name,
                     boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntData &data);
  MoFEMErrorCode aSsemble(EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
template <bool LOGSTRAIN> struct OpPostProcElastic : public DomainEleOp {
  OpPostProcElastic(const std::string field_name,
                    moab::Interface &post_proc_mesh,
                    std::vector<EntityHandle> &map_gauss_pts,
                    boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  moab::Interface &postProcMesh;
  std::vector<EntityHandle> &mapGaussPts;
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpSetMaterialBlock : public DomainEleOp {
  OpSetMaterialBlock(const std::string field_name,
                     boost::shared_ptr<CommonData> common_data_ptr,
                     std::map<EntityHandle, int> &map_block_tets);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  std::map<EntityHandle, int> mapBlockTets;
};
 
struct OpSetMaterialBlockBoundary : public BoundaryEleOp {
 
  OpSetMaterialBlockBoundary(const std::string field_name,
                             boost::shared_ptr<CommonData> common_data_ptr,
                             std::map<EntityHandle, int> &map_block_tets);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  std::map<EntityHandle, int> mapBlockTets;
};
 
struct OpSaveReactionForces : public DomainEleOp {
  OpSaveReactionForces(MoFEM::Interface &m_field, const std::string field_name,
                       moab::Interface &post_proc_mesh,
                       std::vector<EntityHandle> &map_gauss_pts,
                       boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  MoFEM::Interface &mField;
  moab::Interface &postProcMesh;
  std::vector<EntityHandle> &mapGaussPts;
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpEssentialRHS_U : public BoundaryEleOpAssembly {
 
  OpEssentialRHS_U(std::string disp_name,
                   boost::shared_ptr<CommonData> common_data_ptr,
                   VectorDouble &disp, const Range &essential_ents)
      : BoundaryEleOpAssembly(disp_name, disp_name, BoundaryEleOp::OPROW),
        essentialEnts(essential_ents), disP(disp),
        commonDataPtr(common_data_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
    MoFEMFunctionBegin;
    const int nb_row_dofs = row_data.getIndices().size();
    const int nb_row_base_functions = row_data.getN().size2();
 
    if (nb_row_dofs) {
 
      EntityHandle ent = getFEEntityHandle();
      if (essentialEnts.find(ent) == essentialEnts.end()) {
        MoFEMFunctionReturnHot(0);
      }
 
      auto t_disp = getFTensor1FromMat<3>(*commonDataPtr->mDispPtr);
 
      const int nb_integration_pts = getGaussPts().size2();
      auto t_row_val = row_data.getFTensor0N();
      auto t_w = getFTensor0IntegrationWeight();
      const double vol = getMeasure();
      auto t_coords = getFTensor1CoordsAtGaussPts();
 
      Tensor1<double, 3> bc_disp(disP(0), disP(1), disP(2));
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        const double a = vol * t_w;
        
        // // for testing radial displacement
        // auto theta = atan2(t_coords(1), t_coords(0));
        // bc_disp(0) = disP(0) * cos(theta);
        // bc_disp(1) = disP(1) * sin(theta);
 
        auto t_nf = getFTensor1FromArray<3, 3>(locF);
        int rr = 0; 
        for (;rr != nb_row_dofs / 3; ++rr) {
          t_nf(i) += a * t_row_val * (t_disp(i) - bc_disp(i));
          ++t_row_val;
          ++t_nf;
        }
        for (; rr < nb_row_base_functions; ++rr)
          ++t_row_val;
        ++t_disp;
        ++t_coords;
        ++t_w;
      }
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  VectorDouble &disP;
  Range essentialEnts;
};
 
struct OpRotate : public BoundaryEleOpAssembly {
 
  OpRotate(std::string disp_name, VectorDouble &omega,
           boost::shared_ptr<VectorDouble> c_coords,
           boost::shared_ptr<Range> ents_ptr)
      : BoundaryEleOpAssembly(disp_name, disp_name, OPROW), entsPtr(ents_ptr),
        sourceVec(omega), centerCoords(c_coords) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
    MoFEMFunctionBegin;
    const int nb_row_dofs = row_data.getIndices().size();
    const int nb_row_base_functions = row_data.getN().size2();
 
    if (nb_row_dofs) {
 
      EntityHandle ent = getFEEntityHandle();
      if (entsPtr->find(ent) == entsPtr->end()) {
        MoFEMFunctionReturnHot(0);
      }
 
      auto t_coords = getFTensor1CoordsAtGaussPts();
 
      const int nb_integration_pts = getGaussPts().size2();
      auto t_row_val = row_data.getFTensor0N();
      auto t_w = getFTensor0IntegrationWeight();
      const double vol = getMeasure();
      FTensor::Tensor1<double, 3> rot_disp;
      auto t_omega = getFTensor1FromArray<3, 3>(sourceVec);
      auto t_cen_coords = getFTensor1FromArray<3, 3>(*centerCoords);
      FTensor::Tensor1<double, 3> c_dist;
 
      auto rot = get_rotation_R(t_omega, 1);
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        const double a = vol * t_w;
        c_dist(i) = t_cen_coords(i) - t_coords(i);
        rot_disp(i) = c_dist(i) - rot(j, i) * c_dist(j);
 
        auto t_nf = getFTensor1FromArray<3, 3>(locF);
        int rr = 0;
        for (; rr != nb_row_dofs / 3; ++rr) {
          t_nf(i) -= a * t_row_val * rot_disp(i);
          ++t_row_val;
          ++t_nf;
        }
        for (; rr < nb_row_base_functions; ++rr)
          ++t_row_val;
 
        ++t_coords;
        ++t_w;
      }
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  VectorDouble locRes;
  ScalarFun betaCoeff;
  boost::shared_ptr<Range> entsPtr;
  VectorDouble &sourceVec;
  boost::shared_ptr<VectorDouble> centerCoords;
};
 
struct OpSetUnSetRowData : public ForcesAndSourcesCore::UserDataOperator {
  OpSetUnSetRowData(std::string field_name, bool yes_set,
                    boost::shared_ptr<CommonData> common_data_ptr)
      : ForcesAndSourcesCore::UserDataOperator(
            field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
        yesSet(yes_set), commonDataPtr(common_data_ptr) {}
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBeginHot;
    MoFEMFunctionReturnHot(0);
    if (!data.getIndices().empty())
      if (auto e_ptr = data.getFieldEntities()[0])
        if (auto stored_data_ptr =
                e_ptr->getSharedStoragePtr<EssentialBcStorage>()) {
          if (yesSet) {
            commonDataPtr->rowIndices = data.getIndices();
            data.getIndices() = stored_data_ptr->entityIndices;
          } else {
            data.getIndices() = commonDataPtr->rowIndices;
          }
        }
 
    MoFEMFunctionReturnHot(0);
  }
 
public:
  bool yesSet;
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpSchurBegin : public DomainEleOp {
  OpSchurBegin(const std::string row_field,
               boost::shared_ptr<CommonData> common_data_ptr)
      : DomainEleOp(row_field, row_field, OPROW),
        commonDataPtr(common_data_ptr) {
    sYmm = false;
  }
  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &data);
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpSchurEnd : public DomainEleOp {
  OpSchurEnd(const std::string row_field,
             boost::shared_ptr<CommonData> common_data_ptr,
             const double eps = 1e-12)
      : DomainEleOp(row_field, row_field, OPROW),
        commonDataPtr(common_data_ptr), eps(eps) {
    sYmm = false;
  }
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  const double eps;
 
  MatrixDouble invMat;
  MatrixDouble K;
  VectorInt iPIV;
  VectorDouble lapackWork;
 
  map<std::string, MatrixDouble> invBlockMat;
  map<std::string, BlockMatContainer> blockMat;
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpSchurBeginBoundary : public BoundaryEleOp {
  OpSchurBeginBoundary(const std::string row_field,
               boost::shared_ptr<CommonData> common_data_ptr)
      : BoundaryEleOp(row_field, row_field, OPROW),
        commonDataPtr(common_data_ptr) {
    sYmm = false;
  }
  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &data);
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpSchurEndBoundary : public BoundaryEleOp {
  OpSchurEndBoundary(const std::string row_field,
             boost::shared_ptr<CommonData> common_data_ptr,
             const double eps = 1e-12)
      : BoundaryEleOp(row_field, row_field, OPROW),
        commonDataPtr(common_data_ptr), eps(eps) {
    sYmm = false;
  }
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  const double eps;
 
  MatrixDouble invMat;
  MatrixDouble K;
  VectorInt iPIV;
  VectorDouble lapackWork;
 
  map<std::string, MatrixDouble> invBlockMat;
  map<std::string, BlockMatContainer> blockMat;
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
template <typename TYPE>
struct OpSchurPreconditionMassContact : public TYPE {
  OpSchurPreconditionMassContact(const std::string row_field_name,
                          const std::string col_field_name, ScalarFun beta,
                          boost::shared_ptr<CommonData> common_data_ptr)
      : TYPE(row_field_name, col_field_name, beta),
        commonDataPtr(common_data_ptr) {}
 
  MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data,
                          const bool trans) {
    MoFEMFunctionBeginHot;
    commonDataPtr->massMatTensor = this->locMat;
    MoFEMFunctionReturnHot(0);
  }
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
// template struct OpSchurPreconditionMassContact<OpMassPrecHTensorBound>;
// template struct OpSchurPreconditionMassContact<OpMassPrecHTensorVol>;
 
typedef struct OpSchurPreconditionMassContact<OpMassPrecHTensorBound> OpSchurPreconditionMassContactBound;
typedef struct OpSchurPreconditionMassContact<OpMassPrecHTensorVol> OpSchurPreconditionMassContactVol;
 
struct OpSchurPreconditionMassTau : public OpMassPrecScalar {
  OpSchurPreconditionMassTau(const std::string row_field_name,
                             const std::string col_field_name, ScalarFun beta,
                             boost::shared_ptr<CommonData> common_data_ptr)
      : OpMassPrecScalar(row_field_name, col_field_name, beta),
        commonDataPtr(common_data_ptr) {}
 
  MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data,
                          const bool trans) {
    MoFEMFunctionBeginHot;
    commonDataPtr->massMatScalar = this->locMat;
    MoFEMFunctionReturnHot(0);
  }
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
 
/**
 * \brief Not used at this stage. Could be used to do some calculations,
 * before assembly of local elements.
 */
struct FePrePostProcess : public FEMethod {
 
  MoFEM::Interface &mField;
  boost::ptr_vector<DataFromMove> &bcData;
 
  boost::shared_ptr<MethodForForceScaling> methodsOpForMove;
  boost::shared_ptr<MethodForForceScaling> methodsOpForRollersDisp;
 
  boost::shared_ptr<CommonData> commonDataPtr;
 
  // SmartPetscObj<Mat> SEpEp;
  // SmartPetscObj<AO> aoSEpEp;
 
  bool runOnceFlag;
 
  FePrePostProcess(MoFEM::Interface &m_field,
                   boost::ptr_vector<DataFromMove> &bc_data,
                   boost::shared_ptr<CommonData> common_data_ptr)
      : mField(m_field), bcData(bc_data), commonDataPtr(common_data_ptr) {}
 
  MoFEMErrorCode preProcess() {
    MoFEMFunctionBegin;
    runOnceFlag = true;
    CHKERR VecSetOption(ts_F, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
 
    switch (ts_ctx) {
    case CTX_TSSETIJACOBIAN: {
      snes_ctx = CTX_SNESSETJACOBIAN;
      snes_B = ts_B;
      // CHKERR MatSetOption(ts_B, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
      // CHKERR MatSetOption(ts_B, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
 
      if (commonDataPtr->SEpEp)
        CHKERR MatZeroEntries(commonDataPtr->SEpEp);
      if (commonDataPtr->STauTau)
        CHKERR MatZeroEntries(commonDataPtr->STauTau);
 
      break;
    }
    case CTX_TSSETIFUNCTION: {
      snes_ctx = CTX_SNESSETFUNCTION;
      snes_f = ts_F;
      break;
    }
    default:
      break;
    }
 
    auto time = ts_t;
 
    // forceDataScaled = forceData;
    // pressureDataScaled = pressureData;
    if (methodsOpForMove)
      for (auto &bdata : bcData) {
        bdata.scaledValues = bdata.values;
        CHKERR MethodForForceScaling::applyScale(this, methodsOpForMove,
                                                 bdata.scaledValues);
      }
 
    //FIXME: obsolete, for debugging get rid of this flag
    if (time <= (*cache).rollers_stop_time) {
 
      if (methodsOpForRollersDisp)
        for (auto &roller : (*cache).rollerDataVec) {
          roller.BodyDispScaled = roller.rollerDisp;
          CHKERR MethodForForceScaling::applyScale(
              this, methodsOpForRollersDisp, roller.BodyDispScaled);
        }
 
      for (auto &roller : (*cache).rollerDataVec) {
        if (roller.methodOpForRollerPosition) {
          roller.BodyDispScaled.clear();
          CHKERR MethodForForceScaling::applyScale(
              this, roller.methodOpForRollerPosition, roller.BodyDispScaled);
        }
        if (roller.methodOpForRollerDirection) {
          roller.BodyDirectionScaled.clear();
          CHKERR MethodForForceScaling::applyScale(
              this, roller.methodOpForRollerDirection,
              roller.BodyDirectionScaled);
        }
      }
    }
 
    // if (snes_ctx == CTX_SNESNONE) {
    //   if (!dofsIndices.empty()) {
    //     CHKERR VecSetValues(snes_x, dofsIndices.size(),
    //     &*dofsIndices.begin(),
    //                         &*dofsValues.begin(), INSERT_VALUES);
    //   }
    //   CHKERR VecAssemblyBegin(snes_x);
    //   CHKERR VecAssemblyEnd(snes_x);
    // }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode postProcess() {
    MoFEMFunctionBeginHot;
    // if (ts_ctx == CTX_TSSETIJACOBIAN)
    //   CHKERR MatView(ts_B, PETSC_VIEWER_STDOUT_SELF);
    auto assemble = [](Mat a) {
      MoFEMFunctionBeginHot;
      if (a) {
        CHKERR MatAssemblyBegin(a, MAT_FINAL_ASSEMBLY);
        CHKERR MatAssemblyEnd(a, MAT_FINAL_ASSEMBLY);
      }
      MoFEMFunctionReturnHot(0);
    };
 
    if (ts_ctx == CTX_TSSETIJACOBIAN) {
 
      if (commonDataPtr->SEpEp)
        CHKERR assemble(commonDataPtr->SEpEp);
      if (commonDataPtr->STauTau)
        CHKERR assemble(commonDataPtr->STauTau);
 
      // if (true)
      //   CHKERR MatView(ts_B, PETSC_VIEWER_DRAW_WORLD);
      // if (commonDataPtr->SEpEp)
      //   CHKERR MatView(commonDataPtr->SEpEp, PETSC_VIEWER_DRAW_WORLD);
      // if (commonDataPtr->STauTau)
      //   CHKERR MatView(commonDataPtr->STauTau, PETSC_VIEWER_DRAW_WORLD);
 
    }
 
    MoFEMFunctionReturnHot(0);
  }
  };
}; // namespace OpElasticTools