EshelbianOperators.cpp

Go to the documentation of this file.

/**
 * \file EshelbianOperators.cpp
 * \example
 * mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianOperators.cpp
 *
 * \brief Implementation of operators
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
 
#include <boost/math/constants/constants.hpp>
 
#include <EshelbianAux.hpp>
 
#include <lapack_wrap.h>
 
#include <Lie.hpp>
#include <MatrixFunction.hpp>
#include <optional>
 
namespace EshelbianPlasticity {
 
OpCalculateExternalPressure::OpCalculateExternalPressure(
    VectorPtr external_pressure_ptr,
    boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
    std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
    : VolUserDataOperator(H1, OPLAST),
      externalPressurePtr(std::move(external_pressure_ptr)),
      externalStrainVecPtr(std::move(external_strain_vec_ptr)),
      scalingMethodsMap(std::move(smv)) {
  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
  doEntities[MBVERTEX] = true;
}
 
MoFEMErrorCode OpCalculateExternalPressure::doWork(int, EntityType,
                                                   EntData &) {
  MoFEMFunctionBegin;
 
  if (!externalPressurePtr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "External-pressure integration-point vector is null");
  }
 
  const int nb_integration_pts = getGaussPts().size2();
  externalPressurePtr->resize(nb_integration_pts, false);
  externalPressurePtr->clear();
  if (!externalStrainVecPtr) {
    MoFEMFunctionReturnHot(0);
  }
 
  double time = 0;
  double time_step = 0;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
    time_step = EshelbianCore::physicalDt;
  } else if ((getFEMethod()->data_ctx & PetscData::CtxSetTime).any()) {
    time = getTStime();
    time_step = getTStimeStep();
  }
  const EntityHandle fe_ent = getFEEntityHandle();
  const std::regex analytical_pattern("(.*)ANALYTICAL_EXTERNALSTRAIN(.*)");
 
  for (const auto &block : *externalStrainVecPtr) {
    if (block.ents.find(fe_ent) == block.ents.end()) {
      continue;
    }
    if (!std::isfinite(block.bulkModulusK)) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "External-strain bulk modulus is not finite in block %s",
              block.blockName.c_str());
    }
 
    VectorDouble external_strain(nb_integration_pts);
    if (std::regex_match(block.blockName, analytical_pattern)) {
      auto reference_coordinates = getCoordsAtGaussPts();
      external_strain = analytical_externalstrain_function(
          time_step, time, nb_integration_pts, reference_coordinates,
          block.blockName);
    } else {
      double scale = 1;
      if (const auto it = scalingMethodsMap.find(block.blockName);
          it != scalingMethodsMap.end()) {
        if (!it->second) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Scaling method is null for external-strain block %s",
                  block.blockName.c_str());
        }
        scale = it->second->getScale(time);
      } else {
        MOFEM_LOG("EP", Sev::warning)
            << "No scaling method found for " << block.blockName;
      }
      std::fill(external_strain.begin(), external_strain.end(),
                scale * block.val);
    }
 
    if (external_strain.size() != nb_integration_pts) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Wrong number of analytical external-strain integration points");
    }
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      const double q = 3 * block.bulkModulusK * external_strain[gg];
      if (!std::isfinite(q)) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "External pressure q is not finite in block %s",
                block.blockName.c_str());
      }
      (*externalPressurePtr)[gg] += q;
      if (!std::isfinite((*externalPressurePtr)[gg])) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "Accumulated external pressure q is not finite in block %s",
                block.blockName.c_str());
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <typename TInvD, typename TRotation, typename TDiffRotation,
          typename TStress>
auto getDiffSpatialGradientDR(TInvD &t_d_u_d_b, TRotation &t_R,
                              TDiffRotation &t_diff_R, TStress &t_P) {
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l, m, n);
  FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM>
      t_d_rotated_p_d_omega;
  t_d_rotated_p_d_omega(i, j, k) = t_diff_R(n, i, k) * t_P(n, j);
 
  FTensor::Dg<double, SPACE_DIM, SPACE_DIM> t_d_b_d_omega;
  t_d_b_d_omega(i, j, k) =
      (t_d_rotated_p_d_omega(i, j, k) || t_d_rotated_p_d_omega(j, i, k)) / 2.;
 
  FTensor::Dg<double, SPACE_DIM, SPACE_DIM> t_d_u_d_omega;
  t_d_u_d_omega(i, j, k) = t_d_u_d_b(i, j, l, m) * t_d_b_d_omega(l, m, k);
 
  FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM> t_d_h_d_omega;
  t_d_h_d_omega(i, j, k) = t_R(i, l) * t_d_u_d_omega(j, l, k);
 
  return std::make_tuple(t_d_b_d_omega, t_d_u_d_omega, t_d_h_d_omega);
}
 
MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork(int side,
                                                             EntityType type,
                                                             EntData &data) {
  MoFEMFunctionBegin;
 
  auto ts_ctx = getTSCtx();
  int nb_integration_pts = getGaussPts().size2();
 
  // space size indices
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, m);
  FTENSOR_INDEX(SPACE_DIM, n);
  FTENSOR_INDEX(SPACE_DIM, o);
  FTENSOR_INDEX(SPACE_DIM, I);
  FTENSOR_INDEX(SPACE_DIM, J);
 
  // sym size indices
  FTENSOR_INDEX(size_symm, L);
 
  auto t_L = FTensor::SymmLTensor<double, 3>();
 
  MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::size(
      *dataAtPts->getStretchTensorAtPts(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor4DdgFromMatType<3, 3, -1, DL>, DL>::size(
      *dataAtPts->getDiffStretchTensorAtPts(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
      *dataAtPts->getStretchH1AtPts(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor4FromMatType<3, 3, 3, 3, -1, DL>, DL>::size(
      *dataAtPts->getDiffStretchH1AtPts(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
      *dataAtPts->getAdjointPdstretchAtPts(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor1FromMatType<size_symm, -1, DL>, DL>::size(
      *dataAtPts->getAdjointPdUAtPts(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor3FromMatType<3, 3, size_symm, -1, DL>, DL>::size(
      *dataAtPts->getAdjointPdUdPAtPts(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor2FromMatType<3, size_symm, -1, DL>, DL>::size(
      *dataAtPts->getAdjointPdUdOmegaAtPts(), nb_integration_pts);
 
  MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
      *dataAtPts->getDeformationGradient(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::get(
      *dataAtPts->getPlasticH(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::size(
      *dataAtPts->getPlasticF(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
      *dataAtPts->getInvPlasticF(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor3FromMatType<3, 3, 3, -1, DL>, DL>::size(
      dataAtPts->hdOmegaAtPts, nb_integration_pts);
  MatrixSizeHelper<GetFTensor3FromMatType<3, 3, size_symm, -1, DL>, DL>::size(
      dataAtPts->hdLogStretchAtPts, nb_integration_pts);
 
  MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
      dataAtPts->leviKirchhoffAtPts, nb_integration_pts);
  MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
      dataAtPts->leviKirchhoff0AtPts, nb_integration_pts);
  MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
      dataAtPts->leviKirchhoffdOmegaAtPts, nb_integration_pts);
  MatrixSizeHelper<GetFTensor2FromMatType<3, size_symm, -1, DL>, DL>::size(
      dataAtPts->leviKirchhoffdLogStreatchAtPts, nb_integration_pts);
  MatrixSizeHelper<GetFTensor3FromMatType<3, 3, 3, -1, DL>, DL>::size(
      dataAtPts->leviKirchhoffPAtPts, nb_integration_pts);
 
  MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
      dataAtPts->rotMatAtPts, nb_integration_pts);
  MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
      *dataAtPts->getEigenVals(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
      *dataAtPts->getEigenVecs(), nb_integration_pts);
  dataAtPts->nbUniq.resize(nb_integration_pts, false);
  MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
      dataAtPts->eigenValsC, nb_integration_pts);
  MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
      dataAtPts->eigenVecsC, nb_integration_pts);
  dataAtPts->nbUniqC.resize(nb_integration_pts, false);
 
  MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::size(
      dataAtPts->logStretch2H1AtPts, nb_integration_pts);
  MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::size(
      dataAtPts->logStretchTotalTensorAtPts, nb_integration_pts);
 
  MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
      dataAtPts->internalStressAtPts, nb_integration_pts);
  dataAtPts->internalStressAtPts.clear();
 
  // Reconstruct the fixed plastic deformation before evaluating volume
  // equations.  Their work-conjugate stress is the Piola transform
  // P_bar = P F_p^T / J_p, whereas reference equilibrium continues to use P.
  auto t_log_plasticH = dataAtPts->getFTensorPlasticH(nb_integration_pts);
  auto t_plasticF_reconstruct =
      dataAtPts->getFTensorPlasticF(nb_integration_pts);
  auto t_invPlasticF_reconstruct =
      dataAtPts->getFTensorInvPlasticF(nb_integration_pts);
  const EigenMatrix::Fun<double> exp_fun = [](const double v) {
    return std::exp(v);
  };
  const EigenMatrix::Fun<double> inv_exp_fun = [](const double v) {
    return std::exp(-v);
  };
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    FTensor::Tensor1<double, SPACE_DIM> t_eigen_vals;
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_eigen_vecs;
    t_eigen_vecs(i, j) = t_log_plasticH(i, j);
    if (computeEigenValuesSymmetric(t_eigen_vecs, t_eigen_vals) != MB_SUCCESS)
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "Failed to diagonalise logarithmic plastic deformation");
    const auto t_exp_plasticH =
        EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, exp_fun);
    const auto t_inv_exp_plasticH =
        EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, inv_exp_fun);
    t_plasticF_reconstruct(i, j) = t_exp_plasticH(i, j);
    t_invPlasticF_reconstruct(i, j) = t_inv_exp_plasticH(i, j);
 
#ifndef NDEBUG
    const double det_plasticF =
        determinantTensor3by3(t_plasticF_reconstruct);
    if (!std::isfinite(det_plasticF) ||
        det_plasticF <= std::numeric_limits<double>::epsilon())
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "Plastic deformation gradient must have a positive determinant; "
              "got %g",
              det_plasticF);
#endif
 
    ++t_log_plasticH;
    ++t_plasticF_reconstruct;
    ++t_invPlasticF_reconstruct;
  }
 
  MatrixDouble intermediate_p_at_pts;
  MatrixDouble intermediate_p0_at_pts;
  auto get_intermediate_p =
      MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                       DL>::size(intermediate_p_at_pts, nb_integration_pts);
  auto get_intermediate_p0 =
      MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                       DL>::size(intermediate_p0_at_pts, nb_integration_pts);
  auto t_reference_P = dataAtPts->getFTensorApproxP(nb_integration_pts);
  auto t_reference_P0 = dataAtPts->getFTensorApproxP0(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
  auto t_intermediate_P = get_intermediate_p();
  auto t_intermediate_P0 = get_intermediate_p0();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    t_intermediate_P(i, j) =
        t_reference_P(i, k) * t_plasticF(j, k) / det_plasticF;
    t_intermediate_P0(i, j) =
        t_reference_P0(i, k) * t_plasticF(j, k) / det_plasticF;
    ++t_reference_P;
    ++t_reference_P0;
    ++t_plasticF;
    ++t_intermediate_P;
    ++t_intermediate_P0;
  }
 
  // Calculated values
  auto t_h = dataAtPts->getFTensorSmallH(getGaussPts().size2());
  auto t_h_domega = dataAtPts->getFTensorSmallHdOmega(getGaussPts().size2());
  auto t_h_dlog_u =
      dataAtPts->getFTensorSmallHdLogStretch(getGaussPts().size2());
  auto t_levi_kirchhoff =
      dataAtPts->getFTensorLeviKirchhoff(getGaussPts().size2());
  auto t_levi_kirchhoff0 =
      dataAtPts->getFTensorLeviKirchhoff0(getGaussPts().size2());
  auto t_levi_kirchhoff_domega =
      dataAtPts->getFTensorLeviKirchhoffdOmega(getGaussPts().size2());
  auto t_levi_kirchhoff_dstreach =
      dataAtPts->getFTensorLeviKirchhoffdLogStretch(getGaussPts().size2());
  auto t_levi_kirchhoff_dP =
      dataAtPts->getFTensorLeviKirchhoffP(getGaussPts().size2());
  auto t_approx_P_adjoint_dstretch =
      dataAtPts->getFTensorAdjointPdstretch(getGaussPts().size2());
  auto t_approx_P_adjoint_log_du =
      dataAtPts->getFTensorAdjointPdU(getGaussPts().size2());
  auto t_approx_P_adjoint_log_du_dP =
      dataAtPts->getFTensorAdjointPdUdP(getGaussPts().size2());
  auto t_approx_P_adjoint_log_du_domega =
      dataAtPts->getFTensorAdjointPdUdOmega(getGaussPts().size2());
  auto t_R = dataAtPts->getFTensorRotMat(getGaussPts().size2());
  auto t_u = dataAtPts->getFTensorStretch(getGaussPts().size2());
  auto t_diff_u = dataAtPts->getFTensorDiffStretch(getGaussPts().size2());
  auto t_eigen_vals = dataAtPts->getFTensorEigenVals(getGaussPts().size2());
  auto t_eigen_vecs = dataAtPts->getFTensorEigenVecs(getGaussPts().size2());
  auto &nbUniq = dataAtPts->nbUniq;
  auto t_nb_uniq =
      FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniq.data().data());
  auto t_eigen_vals_C = dataAtPts->getFTensorEigenValsC(nb_integration_pts);
  auto t_eigen_vecs_C = dataAtPts->getFTensorEigenVecsC(nb_integration_pts);
  auto &nbUniqC = dataAtPts->nbUniqC;
  auto t_nb_uniq_C =
      FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniqC.data().data());
 
  auto t_u_h1 = dataAtPts->getFTensorStretchH1(getGaussPts().size2());
  auto t_diff_u_h1 = dataAtPts->getFTensorDiffStretchH1(getGaussPts().size2());
  auto t_log_stretch_total =
      dataAtPts->getFTensorLogStretchTotal(getGaussPts().size2());
  auto t_log_u2_h1 = dataAtPts->getFTensorLogStretch2H1(getGaussPts().size2());
 
  // Field values
  auto t_grad_h1 = dataAtPts->getFTensorSmallWGradH1(getGaussPts().size2());
  auto t_omega = dataAtPts->getFTensorRotAxis(getGaussPts().size2());
  auto t_approx_P = get_intermediate_p();
  auto t_approx_P0 = get_intermediate_p0();
  auto t_log_u = dataAtPts->getFTensorLogStretch(getGaussPts().size2());
 
  // Rot axis 0
  auto t_omega0 = dataAtPts->getFTensorRotAxis0(getGaussPts().size2());
  auto t_log_u0 = dataAtPts->getFTensorLogStretch0(getGaussPts().size2());
 
  auto next = [&]() {
    // calculated values
    ++t_h;
    ++t_h_domega;
    ++t_h_dlog_u;
    ++t_levi_kirchhoff;
    ++t_levi_kirchhoff0;
    ++t_levi_kirchhoff_domega;
    ++t_levi_kirchhoff_dstreach;
    ++t_levi_kirchhoff_dP;
    ++t_approx_P_adjoint_dstretch;
    ++t_approx_P_adjoint_log_du;
    ++t_approx_P_adjoint_log_du_dP;
    ++t_approx_P_adjoint_log_du_domega;
    ++t_R;
    ++t_u;
    ++t_diff_u;
    ++t_eigen_vals;
    ++t_eigen_vecs;
    ++t_nb_uniq;
    ++t_eigen_vals_C;
    ++t_eigen_vecs_C;
    ++t_nb_uniq_C;
    ++t_u_h1;
    ++t_diff_u_h1;
    ++t_log_u2_h1;
    ++t_log_stretch_total;
    // field values
    ++t_omega;
    ++t_omega0;
    ++t_grad_h1;
    ++t_approx_P;
    ++t_approx_P0;
    ++t_log_u;
    ++t_log_u0;
  };
 
  auto t_kd = FTensor::Kronecker_Delta<double>();
  auto t_symm_kd = FTensor::Kronecker_Delta_symmetric<double>();
  constexpr auto t_diff_sym = FTensor::DiffSymmetrize<double>();
 
  auto calculate_stretch_from_log = [&](auto &t_log_u_src, auto &t_u_dst,
                                        auto &t_eigen_vals_dst,
                                        auto &t_eigen_vecs_dst,
                                        int &nb_uniq_dst) {
    MoFEMFunctionBegin;
    FTensor::Tensor1<double, 3> eig;
    FTensor::Tensor2<double, 3, 3> eigen_vec;
    eigen_vec(i, j) = t_log_u_src(i, j);
    if (computeEigenValuesSymmetric(eigen_vec, eig) != MB_SUCCESS) {
      MOFEM_LOG("SELF", Sev::error) << "Failed to compute eigen values";
    }
    // CHKERR bound_eig(eig);
    // rare case when two eigen values are equal
    nb_uniq_dst = getUniqNb<3>(eig);
    if (nb_uniq_dst < 3) {
      CHKERR sortEigenVals<3>(eig, eigen_vec);
    }
    t_eigen_vals_dst(i) = eig(i);
    t_eigen_vecs_dst(i, j) = eigen_vec(i, j);
    t_u_dst(i, j) = EigenMatrix::getMat(t_eigen_vals_dst, t_eigen_vecs_dst,
                                        EshelbianCore::f)(i, j);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_log_stretch = [&]() {
    MoFEMFunctionBegin;
    int nb_uniq_val = 0;
    CHKERR calculate_stretch_from_log(t_log_u, t_u, t_eigen_vals, t_eigen_vecs,
                                      nb_uniq_val);
    t_nb_uniq = nb_uniq_val;
    auto get_t_diff_u = [&]() {
      return EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs,
                                     EshelbianCore::f, EshelbianCore::d_f,
                                     t_nb_uniq);
    };
    t_diff_u(i, j, k, l) = get_t_diff_u()(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);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_total_stretch = [&](auto &t_h1) {
    MoFEMFunctionBegin;
    if (EshelbianCore::gradApproximator == NO_H1_CONFIGURATION) {
 
      t_log_u2_h1(i, j) = 0;
      t_log_stretch_total(i, j) = t_log_u(i, j);
 
    } else {
 
      FTensor::Tensor1<double, 3> t_eig_C;
      FTensor::Tensor1<double, 3> t_coordinate_stretch;
      FTensor::Tensor2<double, 3, 3> t_eigen_vec;
 
      FTensor::Tensor2<double, 3, 3> t_C_h1;
      t_C_h1(i, j) = t_h1(k, i) * t_h1(k, j);
      t_eigen_vec(i, j) = t_C_h1(i, j);
      if (computeEigenValuesSymmetric(t_eigen_vec, t_eig_C) != MB_SUCCESS) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                "Failed to compute eigenvalues of F_H1^T F_H1");
      }
      // rare case when two eigen values are equal
      t_nb_uniq_C = getUniqNb<3>(t_eig_C);
      if (t_nb_uniq_C < 3) {
        CHKERR sortEigenVals<3>(t_eig_C, t_eigen_vec);
      }
      for (int aa = 0; aa != 3; ++aa) {
        if (!std::isfinite(t_eig_C(aa)) || t_eig_C(aa) <= 0.) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                  "F_H1^T F_H1 must be positive definite; eigenvalue %d is "
                  "%g",
                  aa, t_eig_C(aa));
        }
        const double principal_stretch = std::sqrt(t_eig_C(aa));
        const double coordinate_stretch =
            EshelbianCore::inv_f(principal_stretch);
        if (!std::isfinite(coordinate_stretch)) {
          SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
                  "Non-finite H1 coordinate stretch for principal stretch %g",
                  principal_stretch);
        }
        t_coordinate_stretch(aa) = coordinate_stretch;
      }
      t_eigen_vals_C(i) = t_eig_C(i);
      t_eigen_vecs_C(i, j) = t_eigen_vec(i, j);
 
      t_log_u2_h1(i, j) =
          EigenMatrix::getMat(t_coordinate_stretch, t_eigen_vec,
                              [](const double v) { return v; })(i, j);
      // The hand-coded Hencky formulation uses additive stretch coordinates.
      // For logarithmic coordinates this is log(U_H1) + log(U_increment).
      t_log_stretch_total(i, j) = t_log_u2_h1(i, j) + t_log_u(i, j);
    }
    MoFEMFunctionReturn(0);
  };
 
  auto no_h1_loop = [&]() {
    MoFEMFunctionBegin;
 
    switch (EshelbianCore::rotSelector) {
    case LARGE_ROT:
      break;
    case SMALL_ROT:
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "no_h1_loop is only implemented for LARGE_ROT");
    };
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      auto t_kd = FTensor::Kronecker_Delta<double>();
 
      FTensor::Tensor2<double, 3, 3> t_h1;
      t_h1(i, j) = t_kd(i, j);
 
      // calculate streach
      CHKERR calculate_log_stretch();
      FTensor::Tensor2<double, 3, 3> t_u0;
      if ((dataAtPts->physicsPtr->getFeatures() &
           PhysicalEquations::noStretchMask)
              .any()) {
        t_u0(i, j) = t_u(i, j);
      } else {
        FTensor::Tensor1<double, 3> t_eigen_vals_0;
        FTensor::Tensor2<double, 3, 3> t_eigen_vecs_0;
        int nb_uniq_0 = 0;
        CHKERR calculate_stretch_from_log(t_log_u0, t_u0, t_eigen_vals_0,
                                          t_eigen_vecs_0, nb_uniq_0);
      }
      // calculate total stretch
      CHKERR calculate_total_stretch(t_h1);
 
      t_u_h1(i, j) = t_u(i, j);
      t_diff_u_h1(i, j, k, l) = t_diff_u(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::Tensor3<double, 3, 3, 3> t_diff_R;
      FTensor::Tensor4<double, 3, 3, 3, 3> t_diff_diff_R;
 
      auto large_rot = [&]() {
        t_R(i, j) = LieGroups::SO3::exp(t_omega, t_omega.l2())(i, j);
        t_diff_R(i, j, k) =
            LieGroups::SO3::diffExp(t_omega, t_omega.l2())(i, j, k);
        t_diff_diff_R(i, j, k, l) =
            LieGroups::SO3::diffDiffExp(t_omega, t_omega.l2())(i, j, k, l);
 
        FTensor::Tensor3<double, 3, 3, 3> t_diff_R0;
        t_diff_R0(i, j, k) =
            LieGroups::SO3::diffExp(t_omega0, t_omega0.l2())(i, j, k);
 
        t_h(i, k) = t_R(i, l) * t_u(l, k);
 
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_rotated_P;
        t_rotated_P(l, k) = t_R(i, l) * t_approx_P(i, k);
        t_approx_P_adjoint_dstretch(l, k) =
            t_diff_sym(l, k, i, j) * t_rotated_P(i, j);
        t_approx_P_adjoint_log_du(L) =
            t_approx_P_adjoint_dstretch(l, k) * t_Ldiff_u(l, k, L);
 
        t_levi_kirchhoff(m) =
            t_diff_R(i, l, m) * (t_u(l, k) * t_approx_P(i, k));
        t_levi_kirchhoff0(m) =
            t_diff_R0(i, l, m) * (t_u0(l, k) * t_approx_P0(i, k));
 
        if (ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
          t_h_domega(i, k, m) = t_diff_R(i, l, m) * t_u(l, k);
          t_h_dlog_u(i, k, L) = t_R(i, l) * t_Ldiff_u(l, k, L);
 
          t_approx_P_adjoint_log_du_dP(i, k, L) =
              t_R(i, l) * t_Ldiff_u(l, k, L);
 
          FTensor::Tensor3<double, 3, 3, 3> t_A;
          t_A(k, l, m) = t_diff_R(i, l, m) * t_approx_P(i, k);
          t_approx_P_adjoint_log_du_domega(m, L) =
              t_A(k, l, m) * t_Ldiff_u(k, l, L);
 
          t_levi_kirchhoff_dstreach(m, L) =
              t_diff_R(i, l, m) * (t_Ldiff_u(l, k, L) * t_approx_P(i, k));
          t_levi_kirchhoff_dP(m, i, k) = t_diff_R(i, l, m) * t_u(l, k);
          t_levi_kirchhoff_domega(m, n) =
              t_diff_diff_R(i, l, m, n) * (t_u(l, k) * t_approx_P(i, k));
 
          if (dataAtPts->physicsPtr->getFeatures().test(
                  PhysicalEquations::NO_STRETCH_NONLINEAR)) {
            auto t_d_u_d_b = GetFTensor4DdgFromMatImpl<
                SPACE_DIM, SPACE_DIM, -1, DL, MatrixDouble>::get(
                dataAtPts->matInvD, gg, 0);
            auto [t_d_b_d_omega, t_d_u_d_omega, t_d_h_d_omega] =
                getDiffSpatialGradientDR(t_d_u_d_b, t_R, t_diff_R,
                                         t_approx_P);
 
            t_h_domega(i, k, m) += t_d_h_d_omega(i, k, m);
            FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM>
                t_d_u_contract_p;
            t_d_u_contract_p(i, l, n) =
                t_d_u_d_omega(l, k, n) * t_approx_P(i, k);
            t_levi_kirchhoff_domega(m, n) +=
                t_diff_R(i, l, m) * t_d_u_contract_p(i, l, n);
 
            if constexpr (EshelbianCore::symmetrySelector > SYMMETRIC) {
              FTensor::Tensor4<double, SPACE_DIM, SPACE_DIM, SPACE_DIM,
                               SPACE_DIM>
                  t_d_b_d_p;
              t_d_b_d_p(i, j, k, l) =
                  t_diff_sym(i, j, m, l) * t_R(k, m);
              FTensor::Tensor4<double, SPACE_DIM, SPACE_DIM, SPACE_DIM,
                               SPACE_DIM>
                  t_d_u_d_p;
              t_d_u_d_p(i, j, k, l) =
                  t_d_u_d_b(i, j, m, n) * t_d_b_d_p(m, n, k, l);
              t_levi_kirchhoff_dP(m, k, l) +=
                  t_d_u_d_p(i, j, k, l) * t_d_b_d_omega(i, j, m);
            }
          }
        }
      };
 
      auto moderate_rot = [&](auto &t_omega0) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "moderate_rot is not implemented yet");
      };
 
      auto small_rot = [&]() {
        t_u_h1(i, j) = t_u(i, j);
        t_diff_u_h1(i, j, k, l) = t_diff_u(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);
 
        t_R(i, j) = t_kd(i, j) + levi_civita(i, j, k) * t_omega(k);
        t_h(i, j) = levi_civita(i, j, k) * t_omega(k) + t_u(i, j);
 
        t_h_domega(i, j, k) = levi_civita(i, j, k);
        t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
 
        // Adjoint stress
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_rotated_P;
        t_rotated_P(i, j) = t_R(k, i) * t_approx_P(k, j);
        t_approx_P_adjoint_dstretch(i, j) =
            t_diff_sym(i, j, k, l) * t_rotated_P(k, l);
        t_approx_P_adjoint_log_du(L) =
            t_approx_P_adjoint_dstretch(i, j) * t_Ldiff_u(i, j, L);
        t_approx_P_adjoint_log_du_dP(i, j, L) =
            t_R(i, k) * t_Ldiff_u(k, j, L);
        t_approx_P_adjoint_log_du_domega(m, L) =
            levi_civita(k, i, m) * t_approx_P(k, j) *
            t_Ldiff_u(i, j, L);
 
        // Kirchhoff stress
        t_levi_kirchhoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
        t_levi_kirchhoff0(k) = levi_civita(i, j, k) * t_approx_P0(i, j);
        t_levi_kirchhoff_dstreach(m, L) = 0;
        t_levi_kirchhoff_dP(k, i, j) = levi_civita(i, j, k);
        t_levi_kirchhoff_domega(m, n) = 0;
      };
 
      // rotation
      switch (EshelbianCore::rotSelector) {
      case LARGE_ROT:
        large_rot();
        break;
      case MODERATE_ROT:
        moderate_rot(t_omega0);
        break;
      case SMALL_ROT:
        small_rot();
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "rotationSelector not handled");
      }
 
      next();
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto large_loop = [&]() {
    MoFEMFunctionBegin;
 
    switch (EshelbianCore::rotSelector) {
    case LARGE_ROT:
      break;
    case SMALL_ROT:
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "rotSelector should be large or small");
    };
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      auto t_kd = FTensor::Kronecker_Delta<double>();
 
      FTensor::Tensor2<double, 3, 3> t_h1;
      switch (EshelbianCore::gradApproximator) {
      case LARGE_ROT:
        t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Selected grad approximator not handled");
      };
 
      // calculate streach
      CHKERR calculate_log_stretch();
      FTensor::Tensor2<double, 3, 3> t_u0;
      if ((dataAtPts->physicsPtr->getFeatures() &
           PhysicalEquations::noStretchMask)
              .any()) {
        t_u0(i, j) = t_u(i, j);
      } else {
        FTensor::Tensor1<double, 3> t_eigen_vals_0;
        FTensor::Tensor2<double, 3, 3> t_eigen_vecs_0;
        int nb_uniq_0 = 0;
        CHKERR calculate_stretch_from_log(t_log_u0, t_u0, t_eigen_vals_0,
                                          t_eigen_vecs_0, nb_uniq_0);
      }
      // calculate total stretch
      CHKERR calculate_total_stretch(t_h1);
 
      t_u_h1(l, k) = t_u(l, o) * t_h1(o, k);
      FTensor::Tensor2<double, 3, 3> t_u_h10;
      t_u_h10(l, k) = t_u0(l, o) * t_h1(o, k);
      t_diff_u_h1(i, j, k, l) = t_diff_u(i, o, k, l) * t_h1(o, j);
      FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u_h1;
      t_Ldiff_u_h1(l, k, L) = t_diff_u_h1(l, k, i, j) * t_L(i, j, L);
 
      FTensor::Tensor3<double, 3, 3, 3> t_diff_R;
      FTensor::Tensor3<double, 3, 3, 3> t_diff_R0;
      FTensor::Tensor4<double, 3, 3, 3, 3> t_diff_diff_R;
 
      // rotation
      switch (EshelbianCore::rotSelector) {
      case SMALL_ROT:
        t_R(i, k) = t_kd(i, k) + levi_civita(i, k, l) * t_omega(l);
        t_diff_R(i, j, k) = levi_civita(i, j, k);
        t_diff_R0(i, j, k) = levi_civita(i, j, k);
        t_diff_diff_R(i, j, l, m) = 0;
        break;
      case LARGE_ROT:
        t_R(i, j) = LieGroups::SO3::exp(t_omega, t_omega.l2())(i, j);
        t_diff_R(i, j, k) =
            LieGroups::SO3::diffExp(t_omega, t_omega.l2())(i, j, k);
        t_diff_R0(i, j, k) =
            LieGroups::SO3::diffExp(t_omega0, t_omega0.l2())(i, j, k);
        t_diff_diff_R(i, j, k, l) =
            LieGroups::SO3::diffDiffExp(t_omega, t_omega.l2())(i, j, k, l);
        break;
 
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "rotationSelector not handled");
      }
 
      // calculate gradient
      t_h(i, k) = t_R(i, l) * t_u_h1(l, k);
 
      // Adjoint stress
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_rotated_P;
      t_rotated_P(l, o) =
          (t_R(i, l) * t_approx_P(i, k)) * t_h1(o, k);
      t_approx_P_adjoint_dstretch(l, o) =
          t_diff_sym(l, o, i, j) * t_rotated_P(i, j);
      t_approx_P_adjoint_log_du(L) =
          t_R(i, l) * t_approx_P(i, k) * t_Ldiff_u_h1(l, k, L);
 
      // Kirchhoff stress
      t_levi_kirchhoff(m) = t_diff_R(i, l, m) * t_u_h1(l, k) * t_approx_P(i, k);
      t_levi_kirchhoff0(m) =
          t_diff_R0(i, l, m) * t_u_h10(l, k) * t_approx_P0(i, k);
 
      if (ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
 
        t_h_domega(i, k, m) = t_diff_R(i, l, m) * t_u_h1(l, k);
        t_h_dlog_u(i, k, L) = t_R(i, l) * t_Ldiff_u_h1(l, k, L);
 
        t_approx_P_adjoint_log_du_dP(i, k, L) =
            t_R(i, l) * t_Ldiff_u_h1(l, k, L);
 
        FTensor::Tensor4<double, 3, size_symm, 3, 3> t_A;
        t_A(m, L, i, k) = t_diff_R(i, l, m) * t_Ldiff_u_h1(l, k, L);
        t_approx_P_adjoint_log_du_domega(m, L) =
            t_A(m, L, i, k) * t_approx_P(i, k);
 
        t_levi_kirchhoff_dstreach(m, L) =
            t_diff_R(i, l, m) * (t_Ldiff_u_h1(l, k, L) * t_approx_P(i, k));
 
        t_levi_kirchhoff_dP(m, i, k) = t_diff_R(i, l, m) * t_u_h1(l, k);
        t_levi_kirchhoff_domega(m, n) =
            t_diff_diff_R(i, l, m, n) * (t_u_h1(l, k) * t_approx_P(i, k));
      }
 
      next();
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto moderate_loop = [&]() {
    MoFEMFunctionBegin;
 
    switch (EshelbianCore::rotSelector) {
    case LARGE_ROT:
      break;
    case SMALL_ROT:
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "rotSelector should be large or small");
    };
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      auto t_kd = FTensor::Kronecker_Delta<double>();
 
      FTensor::Tensor2<double, 3, 3> t_h1;
      switch (EshelbianCore::gradApproximator) {
      case MODERATE_ROT:
        t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Selected grad approximator not handled");
      };
 
      // calculate streach
      CHKERR calculate_log_stretch();
      // calculate total stretch
      CHKERR calculate_total_stretch(t_h1);
 
      auto t_diff = FTensor::DiffTensor<double>();
 
      t_u_h1(l, k) = (t_kd(l, o) + t_log_u(l, o)) * t_h1(o, k);
      FTensor::Tensor2<double, 3, 3> t_u_h10;
      t_u_h10(l, k) = (t_kd(l, o) + t_log_u0(l, o)) * t_h1(o, k);
      t_diff_u_h1(i, j, k, l) = t_diff(i, o, k, l) * t_h1(o, j);
      FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u_h1;
      t_Ldiff_u_h1(l, k, L) = t_diff_u_h1(l, k, i, j) * t_L(i, j, L);
 
      FTensor::Tensor3<double, 3, 3, 3> t_diff_R;
      FTensor::Tensor3<double, 3, 3, 3> t_diff_R0;
      FTensor::Tensor4<double, 3, 3, 3, 3> t_diff_diff_R;
 
      // rotation
      switch (EshelbianCore::rotSelector) {
      case SMALL_ROT:
        t_R(i, k) = t_kd(i, k) + levi_civita(i, k, l) * t_omega(l);
        t_diff_R(i, j, k) = levi_civita(i, j, k);
        t_diff_R0(i, j, k) = levi_civita(i, j, k);
        t_diff_diff_R(i, j, l, m) = 0;
        break;
      case LARGE_ROT:
        t_R(i, j) = LieGroups::SO3::exp(t_omega, t_omega.l2())(i, j);
        t_diff_R(i, j, k) =
            LieGroups::SO3::diffExp(t_omega, t_omega.l2())(i, j, k);
        t_diff_R0(i, j, k) =
            LieGroups::SO3::diffExp(t_omega0, t_omega0.l2())(i, j, k);
        t_diff_diff_R(i, j, k, l) =
            LieGroups::SO3::diffDiffExp(t_omega, t_omega.l2())(i, j, k, l);
        break;
 
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "rotationSelector not handled");
      }
 
      // calculate gradient
      t_h(i, k) = t_R(i, l) * t_u_h1(l, k);
 
      // Adjoint stress
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_rotated_P;
      t_rotated_P(l, o) =
          (t_R(i, l) * t_approx_P(i, k)) * t_h1(o, k);
      t_approx_P_adjoint_dstretch(l, o) =
          t_diff_sym(l, o, i, j) * t_rotated_P(i, j);
      t_approx_P_adjoint_log_du(L) =
          t_R(i, l) * t_approx_P(i, k) * t_Ldiff_u_h1(l, k, L);
 
      // Kirchhoff stress
      t_levi_kirchhoff(m) = t_diff_R(i, l, m) * t_u_h1(l, k) * t_approx_P(i, k);
      t_levi_kirchhoff0(m) =
          t_diff_R0(i, l, m) * t_u_h10(l, k) * t_approx_P0(i, k);
 
      if (ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
 
        t_h_domega(i, k, m) = t_diff_R(i, l, m) * t_u_h1(l, k);
        t_h_dlog_u(i, k, L) = t_R(i, l) * t_Ldiff_u_h1(l, k, L);
 
        t_approx_P_adjoint_log_du_dP(i, k, L) =
            t_R(i, l) * t_Ldiff_u_h1(l, k, L);
 
        FTensor::Tensor4<double, 3, size_symm, 3, 3> t_A;
        t_A(m, L, i, k) = t_diff_R(i, l, m) * t_Ldiff_u_h1(l, k, L);
        t_approx_P_adjoint_log_du_domega(m, L) =
            t_A(m, L, i, k) * t_approx_P(i, k);
 
        t_levi_kirchhoff_dstreach(m, L) =
            t_diff_R(i, l, m) * (t_Ldiff_u_h1(l, k, L) * t_approx_P(i, k));
 
        t_levi_kirchhoff_dP(m, i, k) = t_diff_R(i, l, m) * t_u_h1(l, k);
        t_levi_kirchhoff_domega(m, n) =
            t_diff_diff_R(i, l, m, n) * (t_u_h1(l, k) * t_approx_P(i, k));
      }
 
      next();
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto small_loop = [&]() {
    MoFEMFunctionBegin;
    switch (EshelbianCore::rotSelector) {
    case SMALL_ROT:
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "rotSelector should be small");
    };
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      FTensor::Tensor2<double, 3, 3> t_h1;
      switch (EshelbianCore::gradApproximator) {
      case SMALL_ROT:
        t_h1(i, j) = t_kd(i, j);
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "gradApproximator not handled");
      };
 
      FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
      if (EshelbianCore::stretchSelector > LINEAR) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "stretchSelector should be linear for small loop");
      } else {
        t_u(i, j) = t_symm_kd(i, j) + t_log_u(i, j);
        t_u_h1(i, j) = t_u(i, j);
        t_diff_u_h1(i, j, k, l) =
            (t_kd(i, k) * t_kd(j, l) + t_kd(i, l) * t_kd(j, k));
        t_diff_u_h1(i, j, k, l) /= 2.;
        t_Ldiff_u(i, j, L) = t_L(i, j, L);
      }
      t_log_u2_h1(i, j) = 0;
      t_log_stretch_total(i, j) = t_log_u(i, j);
 
      t_R(i, j) = t_kd(i, j) + levi_civita(i, j, k) * t_omega(k);
      t_h(i, j) = levi_civita(i, j, k) * t_omega(k) + t_u(i, j);
 
      t_h_domega(i, j, k) = levi_civita(i, j, k);
      t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
 
      // Adjoint stress
      t_approx_P_adjoint_dstretch(i, j) =
          t_diff_sym(i, j, k, l) * t_approx_P(k, l);
      t_approx_P_adjoint_log_du(L) =
          t_approx_P_adjoint_dstretch(i, j) * t_Ldiff_u(i, j, L);
      t_approx_P_adjoint_log_du_dP(i, j, L) = t_Ldiff_u(i, j, L);
      t_approx_P_adjoint_log_du_domega(m, L) = 0;
 
      // Kirchhoff stress
      t_levi_kirchhoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
      t_levi_kirchhoff0(k) = levi_civita(i, j, k) * t_approx_P0(i, j);
      t_levi_kirchhoff_dstreach(m, L) = 0;
      t_levi_kirchhoff_dP(k, i, j) = levi_civita(i, j, k);
      t_levi_kirchhoff_domega(m, n) = 0;
 
      next();
    }
 
    MoFEMFunctionReturn(0);
  };
 
  switch (EshelbianCore::gradApproximator) {
  case NO_H1_CONFIGURATION:
    CHKERR no_h1_loop();
    break;
  case LARGE_ROT:
    CHKERR large_loop();
    break;
  case MODERATE_ROT:
    CHKERR moderate_loop();
    break;
  case SMALL_ROT:
    CHKERR small_loop();
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "gradApproximator not handled");
    break;
  };
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpCalculateTractionFromSideEle::doWork(int side, EntityType type,
                                                     EntData &data) {
  MoFEMFunctionBegin;
  FTENSOR_INDEXES(SPACE_DIM, i, j);
 
  auto n_in_the_loop = getNinTheLoop();
  auto loop_size = getLoopSize();
  auto sense = getSkeletonSense();
  auto nb_gauss_pts = getGaussPts().size2();
  auto t_normal = getFTensor1NormalsAtGaussPts();
 
  auto t_sigma = dataAtPts->getFTensorApproxP(getGaussPts().size2());
  auto get_tracion =
      MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
          dataAtPts->tractionAtPts, nb_gauss_pts);
  if (!n_in_the_loop) {
    dataAtPts->tractionAtPts.clear();
  }
 
  auto t_traction = get_tracion();
  for (int gg = 0; gg != nb_gauss_pts; gg++) {
    t_traction(i) +=
        t_sigma(i, j) * sense * (t_normal(j) / t_normal.l2()) / loop_size;
    ++t_traction;
    ++t_sigma;
    ++t_normal;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpCalculateReactionForces::doWork(int side, EntityType type,
                                                 EntData &data) {
  MoFEMFunctionBegin;
  if (blockEntities.find(getFEEntityHandle()) == blockEntities.end()) {
    MoFEMFunctionReturnHot(0);
  };
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  int nb_integration_pts = getGaussPts().size2();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_traction = dataAtPts->getFTensorTraction(nb_integration_pts);
  auto t_coords = getFTensor1CoordsAtGaussPts();
  auto t_spatial_disp = dataAtPts->getFTensorSmallWL2(nb_integration_pts);
 
  FTensor::Tensor1<double, 3> t_coords_spatial{0., 0., 0.};
  FTensor::Tensor1<double, 3> loc_reaction_forces{0., 0., 0.};
  FTensor::Tensor1<double, 3> loc_moment_forces{0., 0., 0.};
 
  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
    double a = t_w * getMeasure();
    loc_reaction_forces(i) += a*t_traction(i);
    t_coords_spatial(i) = t_coords(i) + t_spatial_disp(i);
    loc_moment_forces(i) +=
        (a * (FTensor::levi_civita<double>(i, j, k) * t_coords_spatial(j))) *
        t_traction(k);
    ++t_coords;
    ++t_spatial_disp;
    ++t_w;
    ++t_traction;
  }
 
  reactionVec[0] += loc_reaction_forces(0);
  reactionVec[1] += loc_reaction_forces(1);
  reactionVec[2] += loc_reaction_forces(2);
  reactionVec[3] += loc_moment_forces(0);
  reactionVec[4] += loc_moment_forces(1);
  reactionVec[5] += loc_moment_forces(2);
 
  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 = dataAtPts->getFTensorDivP(nb_integration_pts);
  auto t_s_dot_w = dataAtPts->getFTensorSmallWL2Dot(nb_integration_pts);
  auto w_l2_dot_dot_at_pts = dataAtPts->getSmallWL2DotDotAtPts();
  const bool reset_w_l2_dot_dot =
      w_l2_dot_dot_at_pts->size1() != nb_integration_pts ||
      w_l2_dot_dot_at_pts->size2() != 3;
  MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
      *w_l2_dot_dot_at_pts, nb_integration_pts);
  if (reset_w_l2_dot_dot) {
    w_l2_dot_dot_at_pts->clear();
  }
  auto t_s_dot_dot_w = dataAtPts->getFTensorSmallWL2DotDot(nb_integration_pts);
 
  auto piola_scale = dataAtPts->piolaScale;
  auto alpha_w = alphaW / piola_scale;
  auto alpha_rho = alphaRho / piola_scale;
 
  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]);
  };
 
  auto next = [&]() {
    ++t_w;
    ++t_div_P;
    ++t_s_dot_w;
    ++t_s_dot_dot_w;
  };
 
  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 * alpha_w * t_s_dot_w(i);
      t_nf(i) += a * t_row_base_fun * alpha_rho * t_s_dot_dot_w(i);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
    next();
  }
 
  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_kirchhoff =
      dataAtPts->getFTensorLeviKirchhoff(nb_integration_pts);
  auto t_omega = dataAtPts->getFTensorRotAxis(nb_integration_pts);
  auto t_omega_grad = dataAtPts->getFTensorRotAxisGrad(nb_integration_pts);
  auto t_omega_dot = dataAtPts->getFTensorRotAxisDot(nb_integration_pts);
  auto t_omega_grad_dot =
      dataAtPts->getFTensorRotAxisGradDot(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
  auto t_invPlasticF = dataAtPts->getFTensorInvPlasticF(nb_integration_pts);
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  auto t_row_grad_fun = data.getFTensor1DiffN<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]);
  };
  // auto time_step = getTStimeStep();
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
    const double det_plasticF = determinantTensor3by3(t_plasticF);
 
    FTensor::Tensor2<double, 3, 3> t_omega_grad_intermediate;
    t_omega_grad_intermediate(k, j) =
        t_omega_grad(k, i) * t_invPlasticF(i, j);
    FTensor::Tensor2<double, 3, 3> t_omega_grad_dot_intermediate;
    t_omega_grad_dot_intermediate(k, j) =
        t_omega_grad_dot(k, i) * t_invPlasticF(i, j);
 
    double a = v * t_w * det_plasticF;
    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_kirchhoff(k);
      t_nf(k) += (a * alphaR) * (t_row_base_fun * t_omega(k));
      FTensor::Tensor1<double, 3> t_row_grad_intermediate;
      t_row_grad_intermediate(j) =
          t_row_grad_fun(i) * t_invPlasticF(i, j);
      t_nf(k) += (a * alphaOmega) *
                 (t_row_grad_intermediate(j) *
                  t_omega_grad_intermediate(k, j));
      t_nf(k) += (a * alphaViscousR /*/ time_step*/) *
                 (t_row_base_fun * t_omega_dot(k));
      t_nf(k) += (a * alphaViscousOmega /*/ time_step*/) *
                 (t_row_grad_intermediate(j) *
                  t_omega_grad_dot_intermediate(k, j));
      ++t_nf;
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
    for (; bb != nb_base_functions; ++bb) {
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
    ++t_w;
    ++t_levi_kirchhoff;
    ++t_omega;
    ++t_omega_grad;
    ++t_omega_dot;
    ++t_omega_grad_dot;
    ++t_plasticF;
    ++t_invPlasticF;
  }
  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();
 
  int nb_base_functions = data.getN().size2() / 3;
  auto t_row_base_fun = data.getFTensor1N<3>();
  FTENSOR_INDEX(3, i);
  FTENSOR_INDEX(3, j);
  FTENSOR_INDEX(3, k);
  FTENSOR_INDEX(3, m);
  FTENSOR_INDEX(3, l);
 
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v[0], &v[1],
                                                              &v[2]);
  };
 
  auto t_h = dataAtPts->getFTensorSmallH(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
  auto t_invPlasticF = dataAtPts->getFTensorInvPlasticF(nb_integration_pts);
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
    auto t_nf = get_ftensor1(nF);
 
    FTensor::Tensor2<double, 3, 3> t_residuum;
    t_residuum(i, j) = t_h(i, j) - t_invPlasticF(i, j);
 
    int bb = 0;
    for (; bb != nb_dofs / 3; ++bb) {
      FTensor::Tensor1<double, 3> t_row_base_piola;
      t_row_base_piola(j) =
          t_plasticF(j, k) * t_row_base_fun(k) / det_plasticF;
      t_nf(i) -= a * t_row_base_piola(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_plasticF;
    ++t_invPlasticF;
  }
 
  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();
 
  int nb_base_functions = data.getN().size2() / 9;
  auto t_row_base_fun = data.getFTensor2N<3, 3>();
  FTENSOR_INDEX(3, i);
  FTENSOR_INDEX(3, j);
  FTENSOR_INDEX(3, k);
  FTENSOR_INDEX(3, m);
  FTENSOR_INDEX(3, l);
 
  auto get_ftensor0 = [](auto &v) {
    return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(&v[0]);
  };
 
  auto t_h = dataAtPts->getFTensorSmallH(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
  auto t_invPlasticF = dataAtPts->getFTensorInvPlasticF(nb_integration_pts);
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
    auto t_nf = get_ftensor0(nF);
 
    FTensor::Tensor2<double, 3, 3> t_residuum;
    t_residuum(i, j) = t_h(i, j) - t_invPlasticF(i, j);
 
    int bb = 0;
    for (; bb != nb_dofs; ++bb) {
      FTensor::Tensor2<double, 3, 3> t_row_base_piola;
      t_row_base_piola(i, j) =
          t_row_base_fun(i, k) * t_plasticF(j, k) / det_plasticF;
      t_nf -= a * t_row_base_piola(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_plasticF;
    ++t_invPlasticF;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialConsistencyDivTerm::integrate(EntData &data) {
  MoFEMFunctionBegin;
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_w_l2 = dataAtPts->getFTensorSmallWL2(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
  auto t_invPlasticF = dataAtPts->getFTensorInvPlasticF(nb_integration_pts);
  int nb_base_functions = data.getN().size2() / 3;
  auto t_row_diff_base_fun = data.getFTensor2DiffN<3, 3>();
  FTENSOR_INDEX(3, i);
  FTENSOR_INDEX(3, j);
  FTENSOR_INDEX(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) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
    auto t_nf = get_ftensor1(nF);
    int bb = 0;
    for (; bb != nb_dofs / 3; ++bb) {
      const double div_row_base =
          (t_plasticF(i, j) * t_row_diff_base_fun(j, k) *
           t_invPlasticF(k, i)) /
          det_plasticF;
      t_nf(i) -= a * div_row_base * t_w_l2(i);
      ++t_nf;
      ++t_row_diff_base_fun;
    }
    for (; bb != nb_base_functions; ++bb) {
      ++t_row_diff_base_fun;
    }
    ++t_w;
    ++t_w_l2;
    ++t_plasticF;
    ++t_invPlasticF;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <>
MoFEMErrorCode OpGetInternalStress<0>::doWork(int side, EntityType type,
                                              EntData &data) {
  MoFEMFunctionBegin;
 
  int nb_integration_pts = getGaussPts().size2();
 
  Tag tag;
  CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
  int tag_length;
  CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
  if (tag_length != 9) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Number of internal stress components should be 9 but is %d",
            tag_length);
  }
 
  VectorDouble const_stress_vec(9);
  auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
  CHKERR getPtrFE() -> mField.get_moab().tag_get_data(
                        tag, &fe_ent, 1, &*const_stress_vec.data().begin());
  auto t_const_stress = getFTensor1FromArray<9, 9>(const_stress_vec);
 
  auto get_internal_stress =
      MatrixSizeHelper<GetFTensor1FromMatType<9, -1, DL>, DL>::size(
          dataAtPts->internalStressAtPts, nb_integration_pts);
  dataAtPts->internalStressAtPts.clear();
  auto t_internal_stress = get_internal_stress();
 
  FTensor::Index<'L', 9> L;
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    t_internal_stress(L) = t_const_stress(L);
    ++t_internal_stress;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <>
MoFEMErrorCode OpGetInternalStress<1>::doWork(int side, EntityType type,
                                              EntData &data) {
  MoFEMFunctionBegin;
 
  int nb_integration_pts = getGaussPts().size2();
 
  Tag tag;
  CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
  int tag_length;
  CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
  if (tag_length != 9) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Number of internal stress components should be 9 but is %d",
            tag_length);
  }
 
  auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
  const EntityHandle *vert_conn;
  int vert_num;
  CHKERR getPtrFE() -> mField.get_moab().get_connectivity(fe_ent, vert_conn,
                                                          vert_num, true);
  VectorDouble vert_data(vert_num * tag_length);
  CHKERR getPtrFE() -> mField.get_moab().tag_get_data(tag, vert_conn, vert_num,
                                                      &vert_data[0]);
 
  auto get_internal_stress =
      MatrixSizeHelper<GetFTensor1FromMatType<9, -1, DL>, DL>::size(
          dataAtPts->internalStressAtPts, nb_integration_pts);
  dataAtPts->internalStressAtPts.clear();
  auto t_internal_stress = get_internal_stress();
 
  auto t_shape_n = data.getFTensor0N();
  int nb_shape_fn = data.getN(NOBASE).size2();
  FTensor::Index<'L', 9> L;
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    auto t_vert_data = getFTensor1FromArray<9, 9>(vert_data);
    for (int bb = 0; bb != nb_shape_fn; ++bb) {
      t_internal_stress(L) += t_vert_data(L) * t_shape_n;
      ++t_vert_data;
      ++t_shape_n;
    }
    ++t_internal_stress;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <>
MoFEMErrorCode
OpSpatialPhysicalInternalStress<false>::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
 
  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]);
  };
 
  auto t_internal_stress =
      dataAtPts->getFTensorInternalStress(nb_integration_pts);
 
  const double time = EshelbianCore::physicalTimeFlg
                          ? EshelbianCore::currentPhysicalTime
                          : getFEMethod()->ts_t;
 
  // default scaling is constant
  double scale = scalingMethodPtr->getScale(time);
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = FTensor::SymmLTensor<double, 3>();
 
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor2(nF);
 
    FTensor::Tensor2<double, 3, 3> t_symm_stress;
    t_symm_stress(i, j) =
        (t_internal_stress(i, j) + t_internal_stress(j, i)) / 2;
 
    FTensor::Tensor1<double, size_symm> t_residual;
    t_residual(L) = t_L(i, j, L) * (scale * t_symm_stress(i, j));
 
    int bb = 0;
    for (; bb != nb_dofs / 6; ++bb) {
      t_nf(L) += a * 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_internal_stress;
  }
  MoFEMFunctionReturn(0);
}
 
template <>
MoFEMErrorCode OpSpatialPhysicalInternalStress<true>::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 get_ftensor2 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 6>, 6>(
        &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
  };
 
  auto t_internal_stress =
      dataAtPts->getFTensorInternalStressVec(nb_integration_pts);
 
  FTensor::Index<'L', size_symm> L;
  FTensor::Index<'M', size_symm> M;
  FTensor::Tensor2<double, size_symm, size_symm> t_L;
  t_L = voigt_to_symm();
 
  const double time = EshelbianCore::physicalTimeFlg
                          ? EshelbianCore::currentPhysicalTime
                          : getFEMethod()->ts_t;
 
  // default is constant
  double scale = scalingMethodPtr->getScale(time);
 
  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::Tensor1<double, size_symm> t_residual;
    t_residual(L) = t_L(M, L) * (scale * t_internal_stress(M));
 
    int bb = 0;
    for (; bb != nb_dofs / 6; ++bb) {
      t_nf(L) += a * 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_internal_stress;
  }
  MoFEMFunctionReturn(0);
}
 
template <AssemblyType A>
MoFEMErrorCode OpDispBcImpl<A, GAUSS>::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcDispPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
      int nb_dofs = data.getIndices().size();
 
      int nb_integration_pts = OP::getGaussPts().size2();
      auto t_normal = OP::getFTensor1NormalsAtGaussPts();
      auto t_w = OP::getFTensor0IntegrationWeight();
      int nb_base_functions = data.getN().size2() / 3;
      auto t_row_base_fun = data.getFTensor1N<3>();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
 
      double scale = 1;
      if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
        if (EshelbianCore::physicalTimeFlg) {
          scale *= scalingMethodsMap.at(bc.blockName)
                       ->getScale(EshelbianCore::currentPhysicalTime);
        } else {
          scale *= scalingMethodsMap.at(bc.blockName)
                       ->getScale(OP::getFEMethod()->ts_t);
        }
      } else {
        MOFEM_LOG("SELF", Sev::warning)
            << "No scaling method found for " << bc.blockName;
      }
 
      // get bc data
      FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
      t_bc_disp(i) *= scale;
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
        int bb = 0;
        for (; bb != nb_dofs / SPACE_DIM; ++bb) {
          t_nf(i) +=
              t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_disp(i) * 0.5;
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
 
        ++t_w;
        ++t_normal;
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpDispBc::iNtegrate(EntData &data) {
  return OP::iNtegrate(data);
}
 
MoFEMErrorCode OpTauStabilizationDispRhsBc::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcDispPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
      int nb_dofs = data.getIndices().size();
 
      int nb_integration_pts = OP::getGaussPts().size2();
      auto t_w = OP::getFTensor0IntegrationWeight();
      int nb_base_functions = data.getN().size2();
      auto t_row_base_fun = data.getFTensor0N();
#ifndef NDEBUG
      if (!this->sourceVec) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Source vector for OpTauStabilizationDispRhsBc is not set");
      }
      if (data.getN().size1() != nb_integration_pts) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Number of integration points in data should be %d but is %d",
                nb_integration_pts, (int)data.getN().size1());
      }
      if (nb_base_functions < nb_dofs / SPACE_DIM) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Number of base functions in data should be %d but is %d",
                nb_base_functions, (int)data.getN().size2() / SPACE_DIM);
      }
 
#endif
 
      auto t_disp_val =
          MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::get(
              *this->sourceVec, nb_integration_pts)();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
 
      double scale = 1;
      if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
        if (EshelbianCore::physicalTimeFlg) {
          scale *= scalingMethodsMap.at(bc.blockName)
                       ->getScale(EshelbianCore::currentPhysicalTime);
        } else {
          scale *= scalingMethodsMap.at(bc.blockName)
                       ->getScale(OP::getFEMethod()->ts_t);
        }
      } else {
        MOFEM_LOG("SELF", Sev::warning)
            << "No scaling method found for " << bc.blockName;
      }
 
      // get bc data
      FTensor::Tensor1<double, 3> t_bc_disp(bc.vals[0], bc.vals[1], bc.vals[2]);
      t_bc_disp(i) *= scale;
 
      auto area = getMeasure();
      auto t_coords = OP::getFTensor1CoordsAtGaussPts();
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        auto tau_scale =
            area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
        FTensor::Tensor1<double, SPACE_DIM> t_bc_residual;
        for (int ii = 0; ii != SPACE_DIM; ++ii)
          t_bc_residual(ii) =
              bc.flags[ii] ? t_disp_val(ii) - t_bc_disp(ii) : 0.;
        auto t_nf = getFTensor1FromPtr<3, 3>(OP::locF.data().data());
        int bb = 0;
        for (; bb != nb_dofs / SPACE_DIM; ++bb) {
          t_nf(i) +=
              (tau_scale * t_row_base_fun) * t_bc_residual(i);
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
 
        ++t_w;
        ++t_coords;
        ++t_disp_val;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpTauStabilizationDispLhsBc::iNtegrate(EntData &row_data,
                                                      EntData &col_data) {
  MoFEMFunctionBegin;
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcDispPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
      int nb_dofs = row_data.getIndices().size();
 
      int nb_integration_pts = OP::getGaussPts().size2();
      auto t_w = OP::getFTensor0IntegrationWeight();
      int nb_base_functions = row_data.getN().size2();
      auto t_row_base_fun = row_data.getFTensor0N();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
 
      FTensor::Tensor1<double, SPACE_DIM> t_bc_mask;
      for (int ii = 0; ii != SPACE_DIM; ++ii)
        t_bc_mask(ii) = bc.flags[ii] ? 1. : 0.;
 
      auto get_t_vec = [&](const int rr) {
        std::array<double *, SPACE_DIM> ptrs;
        for (auto i = 0; i != SPACE_DIM; ++i)
          ptrs[i] = &OP::locMat(rr + i, i);
        return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>,
                                SPACE_DIM>(ptrs);
      };
 
      auto area = getMeasure();
      auto t_coords = OP::getFTensor1CoordsAtGaussPts();
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        auto tau_scale =
            area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
        int rr = 0;
        for (; rr != nb_dofs / SPACE_DIM; ++rr) {
          auto t_mat = get_t_vec(SPACE_DIM * rr);
          auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
          for (int cc = 0; cc != nb_dofs / SPACE_DIM; ++cc) {
            t_mat(i) +=
                tau_scale * (t_row_base_fun * t_col_base_fun) * t_bc_mask(i);
            ++t_col_base_fun;
            ++t_mat;
          }
          ++t_row_base_fun;
        }
        for (; rr != nb_base_functions; ++rr)
          ++t_row_base_fun;
 
        ++t_w;
        ++t_coords;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpTauStabilizationOpAnalyticalDispBc::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  EntityHandle fe_ent = OP::getFEEntityHandle();
  for (auto &bc : (*bcDispPtr)) {
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
      MatrixDouble analytical_expr;
      auto analytical_data = getAnalyticalExpr(this, analytical_expr, bc.blockName);
      auto &v_analytical_expr = std::get<1>(analytical_data);
 
      int nb_dofs = data.getIndices().size();
      int nb_integration_pts = OP::getGaussPts().size2();
      auto t_w = OP::getFTensor0IntegrationWeight();
      int nb_base_functions = data.getN().size2();
      auto t_row_base_fun = data.getFTensor0N();
 
#ifndef NDEBUG
      if (!this->sourceVec) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Source vector for OpTauStabilizationOpAnalyticalDispBc is not "
                "set");
      }
      if (data.getN().size1() != nb_integration_pts) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Number of integration points in data should be %d but is %d",
                nb_integration_pts, (int)data.getN().size1());
      }
      if (nb_base_functions < nb_dofs / SPACE_DIM) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Number of base functions in data should be at least %d but is "
                "%d",
                nb_dofs / SPACE_DIM, nb_base_functions);
      }
#endif
 
      auto t_disp_val =
          MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::get(
              *this->sourceVec, nb_integration_pts)();
      auto t_bc_disp = getFTensor1FromMat<3, -1, DL>(v_analytical_expr);
 
      FTENSOR_INDEX(SPACE_DIM, i);
 
      auto area = getMeasure();
      auto t_coords = OP::getFTensor1CoordsAtGaussPts();
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        auto tau_scale =
            area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
        FTensor::Tensor1<double, SPACE_DIM> t_bc_residual;
        for (int ii = 0; ii != SPACE_DIM; ++ii)
          t_bc_residual(ii) =
              bc.flags[ii] ? t_disp_val(ii) - t_bc_disp(ii) : 0.;
        auto t_nf = getFTensor1FromPtr<3, 3>(OP::locF.data().data());
        int bb = 0;
        for (; bb != nb_dofs / SPACE_DIM; ++bb) {
          t_nf(i) +=
              (tau_scale * t_row_base_fun) * t_bc_residual(i);
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
 
        ++t_w;
        ++t_coords;
        ++t_disp_val;
        ++t_bc_disp;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpTauStabilizationOpAnalyticalDispLhsBc::iNtegrate(
    EntData &row_data, EntData &col_data) {
  MoFEMFunctionBegin;
 
  EntityHandle fe_ent = OP::getFEEntityHandle();
  for (auto &bc : (*bcDispPtr)) {
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
      int nb_dofs = row_data.getIndices().size();
      int nb_integration_pts = OP::getGaussPts().size2();
      auto t_w = OP::getFTensor0IntegrationWeight();
      int nb_base_functions = row_data.getN().size2();
      auto t_row_base_fun = row_data.getFTensor0N();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTensor::Tensor1<double, SPACE_DIM> t_bc_mask;
      for (int ii = 0; ii != SPACE_DIM; ++ii)
        t_bc_mask(ii) = bc.flags[ii] ? 1. : 0.;
 
      auto get_t_vec = [&](const int rr) {
        std::array<double *, SPACE_DIM> ptrs;
        for (auto i = 0; i != SPACE_DIM; ++i)
          ptrs[i] = &OP::locMat(rr + i, i);
        return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>,
                                SPACE_DIM>(ptrs);
      };
 
      auto area = getMeasure();
      auto t_coords = OP::getFTensor1CoordsAtGaussPts();
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        auto tau_scale =
            area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
        int rr = 0;
        for (; rr != nb_dofs / SPACE_DIM; ++rr) {
          auto t_mat = get_t_vec(SPACE_DIM * rr);
          auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
          for (int cc = 0; cc != nb_dofs / SPACE_DIM; ++cc) {
            t_mat(i) +=
                tau_scale * (t_row_base_fun * t_col_base_fun) * t_bc_mask(i);
            ++t_col_base_fun;
            ++t_mat;
          }
          ++t_row_base_fun;
        }
        for (; rr != nb_base_functions; ++rr)
          ++t_row_base_fun;
 
        ++t_w;
        ++t_coords;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <AssemblyType A>
MoFEMErrorCode OpRotationBcImpl<A, GAUSS>::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(3, i);
  FTENSOR_INDEX(3, j);
  FTENSOR_INDEX(3, k);
 
  double time = OP::getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
 
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // interate over all boundary data
  for (auto &bc : (*bcRotPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
      int nb_dofs = data.getIndices().size();
      int nb_integration_pts = OP::getGaussPts().size2();
      auto t_normal = OP::getFTensor1NormalsAtGaussPts();
      auto t_w = OP::getFTensor0IntegrationWeight();
 
      int nb_base_functions = data.getN().size2() / 3;
      auto t_row_base_fun = data.getFTensor1N<3>();
 
      // Note: First three values of bc.vals are the center of rotation
      // 4th is rotation angle in radians, and remaining values are axis of
      // rotation. Also, if rotation axis is not provided, it defaults to the
      // normal vector of the face.
 
      // get bc data
      FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
 
      auto get_rotation_angle = [&]() {
        double theta = bc.theta;
        if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
          theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
        }
        return theta;
      };
 
      auto get_rotation = [&](auto theta) {
        FTensor::Tensor1<double, 3> t_omega;
        if (bc.vals.size() == 7) {
          t_omega(0) = bc.vals[4];
          t_omega(1) = bc.vals[5];
          t_omega(2) = bc.vals[6];
        } else {
          // Use gemetric face normal as rotation axis
          t_omega(i) = OP::getFTensor1Normal()(i);
        }
        if (t_omega.l2() > std::numeric_limits<double>::epsilon()) {
          t_omega.normalize();
        } else {
          MOFEM_LOG("SELF", Sev::warning)
              << "Rotation axis is zero vector for block " << bc.blockName
              << ". This may lead to unexpected results.";
        }
        t_omega(i) *= theta;
        return LieGroups::SO3::exp(t_omega, EshelbianCore::rotSelector ==
                                                    RotSelector::SMALL_ROT
                                                ? 0.
                                                : t_omega.l2());
      };
 
      auto t_R = get_rotation(get_rotation_angle());
      auto t_coords = OP::getFTensor1CoordsAtGaussPts();
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        FTensor::Tensor1<double, 3> t_delta;
        t_delta(i) = t_center(i) - t_coords(i);
        FTensor::Tensor1<double, 3> t_disp;
        t_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
 
        auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
        int bb = 0;
        for (; bb != nb_dofs / SPACE_DIM; ++bb) {
          t_nf(i) += t_w * (t_row_base_fun(j) * t_normal(j)) * t_disp(i) * 0.5;
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
 
        ++t_w;
        ++t_normal;
        ++t_coords;
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpRotationBc::iNtegrate(EntData &data) {
  return OP::iNtegrate(data);
}
 
MoFEMErrorCode OpTauStabilizationRotationRhsBc::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(3, i);
  FTENSOR_INDEX(3, j);
 
  double time = OP::getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
 
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcRotPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
      int nb_dofs = data.getIndices().size();
      int nb_integration_pts = OP::getGaussPts().size2();
      auto t_w = OP::getFTensor0IntegrationWeight();
 
      int nb_base_functions = data.getN().size2();
      auto t_row_base_fun = data.getFTensor0N();
 
      // Note: First three values of bc.vals are the center of rotation
      // 4th is rotation angle in radians, and remaining values are axis of
      // rotation. Also, if rotation axis is not provided, it defaults to the
      // normal vector of the face.
 
      // get bc data
      FTensor::Tensor1<double, 3> t_center(bc.vals[0], bc.vals[1], bc.vals[2]);
 
      auto get_rotation_angle = [&]() {
        double theta = bc.theta;
        if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
          theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
        }
        return theta;
      };
 
      auto get_rotation = [&](auto theta) {
        FTensor::Tensor1<double, 3> t_omega;
        if (bc.vals.size() == 7) {
          t_omega(0) = bc.vals[4];
          t_omega(1) = bc.vals[5];
          t_omega(2) = bc.vals[6];
        } else {
          // Use gemetric face normal as rotation axis
          t_omega(i) = OP::getFTensor1Normal()(i);
        }
        if (t_omega.l2() > std::numeric_limits<double>::epsilon()) {
          t_omega.normalize();
        }
        t_omega(i) *= theta;
        return LieGroups::SO3::exp(t_omega, EshelbianCore::rotSelector ==
                                                    RotSelector::SMALL_ROT
                                                ? 0.
                                                : t_omega.l2());
      };
 
      auto area = getMeasure();
      auto t_R = get_rotation(get_rotation_angle());
      auto t_coords = OP::getFTensor1CoordsAtGaussPts();
      auto t_disp_val =
          MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::get(
              *this->sourceVec, nb_integration_pts)();
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        auto tau_scale =
            area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
 
        FTensor::Tensor1<double, 3> t_delta;
        t_delta(i) = t_center(i) - t_coords(i);
        FTensor::Tensor1<double, 3> t_bc_disp;
        t_bc_disp(i) = t_delta(i) - t_R(i, j) * t_delta(j);
 
        auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
        int bb = 0;
        for (; bb != nb_dofs / SPACE_DIM; ++bb) {
          t_nf(i) +=
              (tau_scale * t_row_base_fun) * (t_disp_val(i) - t_bc_disp(i));
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
 
        ++t_w;
        ++t_coords;
        ++t_disp_val;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpTauStabilizationRotationLhsBc::iNtegrate(EntData &row_data,
                                                          EntData &col_data) {
  MoFEMFunctionBegin;
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate 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 = row_data.getIndices().size();
 
      int nb_integration_pts = OP::getGaussPts().size2();
      auto t_w = OP::getFTensor0IntegrationWeight();
      int nb_base_functions = row_data.getN().size2();
      auto t_row_base_fun = row_data.getFTensor0N();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
 
      auto get_t_vec = [&](const int rr) {
        std::array<double *, SPACE_DIM> ptrs;
        for (auto i = 0; i != SPACE_DIM; ++i)
          ptrs[i] = &OP::locMat(rr + i, i);
        return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>,
                                SPACE_DIM>(ptrs);
      };
 
      auto area = getMeasure();
      auto t_coords = OP::getFTensor1CoordsAtGaussPts();
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        auto tau_scale =
            area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
        int rr = 0;
        for (; rr != nb_dofs / SPACE_DIM; ++rr) {
          auto t_mat = get_t_vec(SPACE_DIM * rr);
          auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
          for (int cc = 0; cc != nb_dofs / SPACE_DIM; ++cc) {
            for (int ii = 0; ii != SPACE_DIM; ++ii) {
              t_mat(ii) += tau_scale * (t_row_base_fun * t_col_base_fun);
            }
            ++t_col_base_fun;
            ++t_mat;
          }
          ++t_row_base_fun;
        }
        for (; rr != nb_base_functions; ++rr)
          ++t_row_base_fun;
 
        ++t_w;
        ++t_coords;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <AssemblyType A>
MoFEMErrorCode OpNormalDispBcRhsImpl<A, GAUSS>::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  double time = OP::getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
 
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcDispPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      for (auto &bd : (*brokenBaseSideDataPtr)) {
 
        auto t_approx_P = getFTensor2FromMat<3, 3, -1, DL>(bd.getFlux());
        auto t_u = getFTensor1FromMat<3, -1, DL>(*hybridDispPtr);
        auto t_normal = OP::getFTensor1NormalsAtGaussPts();
        auto t_w = OP::getFTensor0IntegrationWeight();
 
        FTENSOR_INDEX(SPACE_DIM, i);
        FTENSOR_INDEX(SPACE_DIM, j);
 
        auto t_kd = FTensor::Kronecker_Delta<double>();
 
        double scale = 1;
        if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
          scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
        } else {
          MOFEM_LOG("SELF", Sev::warning)
              << "No scaling method found for " << bc.blockName;
        }
 
        // get bc data
        double val = scale * bc.val;
 
        int nb_dofs = data.getIndices().size();
        int nb_integration_pts = OP::getGaussPts().size2();
        int nb_base_functions = data.getN().size2();
        auto t_row_base = data.getFTensor0N();
        for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
          FTensor::Tensor1<double, 3> t_N;
          t_N(i) = t_normal(i);
          t_N.normalize();
 
          FTensor::Tensor2<double, 3, 3> t_P;
          t_P(i, j) = t_N(i) * t_N(j);
          FTensor::Tensor2<double, 3, 3> t_Q;
          t_Q(i, j) = t_kd(i, j) - t_P(i, j);
 
          FTensor::Tensor1<double, 3> t_traction;
          t_traction(i) = t_approx_P(i, j) * t_N(j);
 
          FTensor::Tensor1<double, 3> t_res;
          t_res(i) =
              t_Q(i, j) * t_traction(j) + t_P(i, j) * 2 * t_u(j) - t_N(i) * val;
 
          auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
          int bb = 0;
          for (; bb != nb_dofs / SPACE_DIM; ++bb) {
            t_nf(i) += (t_w * t_row_base * OP::getMeasure()) * t_res(i);
            ++t_nf;
            ++t_row_base;
          }
          for (; bb != nb_base_functions; ++bb)
            ++t_row_base;
 
          ++t_w;
          ++t_normal;
          ++t_u;
          ++t_approx_P;
        }
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
template <AssemblyType A>
MoFEMErrorCode
OpNormalDispBcLhsImpl_dU<A, GAUSS>::iNtegrate(EntData &row_data,
                                              EntData &col_data) {
  MoFEMFunctionBegin;
 
  double time = OP::getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
 
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
  auto &locMat = OP::locMat;
  locMat.resize(row_nb_dofs, col_nb_dofs, false);
  locMat.clear();
 
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcDispPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      auto t_normal = OP::getFTensor1NormalsAtGaussPts();
      auto t_w = OP::getFTensor0IntegrationWeight();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
 
      double scale = 1;
      if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
        scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
      } else {
        MOFEM_LOG("SELF", Sev::warning)
            << "No scaling method found for " << bc.blockName;
      }
 
      int nb_integration_pts = OP::getGaussPts().size2();
      int row_nb_dofs = row_data.getIndices().size();
      int col_nb_dofs = col_data.getIndices().size();
      int nb_base_functions = row_data.getN().size2();
      auto t_row_base = row_data.getFTensor0N();
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
        FTensor::Tensor1<double, SPACE_DIM> t_N;
        t_N(i) = t_normal(i);
        t_N.normalize();
 
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_P;
        t_P(i, j) = t_N(i) * t_N(j);
 
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_d_res;
        t_d_res(i, j) = 2.0 * t_P(i, j);
 
        int rr = 0;
        for (; rr != row_nb_dofs / SPACE_DIM; ++rr) {
          auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
              locMat, SPACE_DIM * rr);
          auto t_col_base = col_data.getFTensor0N(gg, 0);
          for (auto cc = 0; cc != col_nb_dofs / SPACE_DIM; ++cc) {
            t_mat(i, j) += (t_w * t_row_base * t_col_base) * t_d_res(i, j);
            ++t_mat;
            ++t_col_base;
          }
          ++t_row_base;
        }
 
        for (; rr != nb_base_functions; ++rr)
          ++t_row_base;
 
        ++t_w;
        ++t_normal;
      }
 
      locMat *= OP::getMeasure();
    }
  }
  MoFEMFunctionReturn(0);
}
 
template <AssemblyType A>
MoFEMErrorCode
OpNormalDispBcLhsImpl_dP<A, GAUSS>::iNtegrate(EntData &row_data,
                                              EntData &col_data) {
  MoFEMFunctionBegin;
 
  double time = OP::getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
 
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
  auto &locMat = OP::locMat;
  locMat.resize(row_nb_dofs, col_nb_dofs, false);
  locMat.clear();
 
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcDispPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      auto t_normal = OP::getFTensor1NormalsAtGaussPts();
      auto t_w = OP::getFTensor0IntegrationWeight();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
      FTENSOR_INDEX(SPACE_DIM, k);
 
      auto t_kd = FTensor::Kronecker_Delta<double>();
 
      double scale = 1;
      if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
        scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
      } else {
        MOFEM_LOG("SELF", Sev::warning)
            << "No scaling method found for " << bc.blockName;
      }
 
      int nb_integration_pts = OP::getGaussPts().size2();
      int nb_base_functions = row_data.getN().size2();
      auto t_row_base = row_data.getFTensor0N();
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
        FTensor::Tensor1<double, SPACE_DIM> t_N;
        t_N(i) = t_normal(i);
        t_N.normalize();
 
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_P;
        t_P(i, j) = t_N(i) * t_N(j);
        FTensor::Tensor2<double, 3, 3> t_Q;
        t_Q(i, j) = t_kd(i, j) - t_P(i, j);
 
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_d_res;
        t_d_res(i, j) = t_Q(i, j);
 
        int rr = 0;
        for (; rr != row_nb_dofs / SPACE_DIM; ++rr) {
          auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
              OP::locMat, SPACE_DIM * rr);
          auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
          for (auto cc = 0; cc != col_nb_dofs / SPACE_DIM; ++cc) {
            t_mat(i, j) +=
                ((t_w * t_row_base) * (t_N(k) * t_col_base(k))) * t_d_res(i, j);
            ++t_mat;
            ++t_col_base;
          }
          ++t_row_base;
        }
 
        for (; rr != nb_base_functions; ++rr)
          ++t_row_base;
 
        ++t_w;
        ++t_normal;
      }
 
      locMat *= OP::getMeasure();
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpNormalDispRhsBc::iNtegrate(EntData &data) {
  return OP::iNtegrate(data);
}
 
MoFEMErrorCode OpNormalDispLhsBc_dU::iNtegrate(EntData &row_data,
                                               EntData &col_data) {
  return OP::iNtegrate(row_data, col_data);
}
 
MoFEMErrorCode OpNormalDispLhsBc_dP::iNtegrate(EntData &row_data,
                                               EntData &col_data) {
  return OP::iNtegrate(row_data, col_data);
}
 
MoFEMErrorCode OpSpringRhsBc::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcSpringPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
      
      for (auto &bd : (*brokenBaseSideDataPtr)) {
 
        auto t_approx_P = getFTensor2FromMat<3, 3, -1, DL>(bd.getFlux());
        auto t_u = getFTensor1FromMat<3, -1, DL>(*hybridDispPtr);
        auto t_normal = OP::getFTensor1NormalsAtGaussPts();
        auto t_w = OP::getFTensor0IntegrationWeight();
 
        FTENSOR_INDEX(SPACE_DIM, i);
        FTENSOR_INDEX(SPACE_DIM, j);
 
        auto t_kd = FTensor::Kronecker_Delta<double>();
 
        int nb_dofs = data.getIndices().size();
        int nb_integration_pts = OP::getGaussPts().size2();
        int nb_base_functions = data.getN().size2();
        auto t_row_base = data.getFTensor0N();
        for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
          FTensor::Tensor1<double, 3> t_N;
          t_N(i) = t_normal(i);
          t_N.normalize();
 
          FTensor::Tensor2<double, 3, 3> t_P;
          t_P(i, j) = t_N(i) * t_N(j);
          FTensor::Tensor2<double, 3, 3> t_Q;
          t_Q(i, j) = t_kd(i, j) - t_P(i, j);
 
          FTensor::Tensor1<double, 3> t_traction;
          t_traction(i) = t_approx_P(i, j) * t_N(j);
 
          FTensor::Tensor1<double, 3> t_res;
          t_res(i) = 0.5 *(t_traction(i)) - bc.normalStiffness * t_P(i, j) * t_u(j) -
                     bc.tangentialStiffness * t_Q(i, j) * t_u(j);
 
          auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
          int bb = 0;
          for (; bb != nb_dofs / SPACE_DIM; ++bb) {
            t_nf(i) += (t_w * t_row_base * OP::getMeasure()) * t_res(i);
            ++t_nf;
            ++t_row_base;
          }
          for (; bb != nb_base_functions; ++bb)
            ++t_row_base;
 
          ++t_w;
          ++t_normal;
          ++t_u;
          ++t_approx_P;
        }
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpringLhsBc_dU::iNtegrate(EntData &row_data,
                                           EntData &col_data) {
  MoFEMFunctionBegin;
 
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
  auto &locMat = OP::locMat;
  locMat.resize(row_nb_dofs, col_nb_dofs, false);
  locMat.clear();
 
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcSpringPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      auto t_normal = OP::getFTensor1NormalsAtGaussPts();
      auto t_w = OP::getFTensor0IntegrationWeight();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
 
      int nb_integration_pts = OP::getGaussPts().size2();
      int nb_base_functions = row_data.getN().size2();
      auto t_row_base = row_data.getFTensor0N();
 
      auto t_kd = FTensor::Kronecker_Delta<double>();
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
        FTensor::Tensor1<double, SPACE_DIM> t_N;
        t_N(i) = t_normal(i);
        t_N.normalize();
 
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_P;
        t_P(i, j) = t_N(i) * t_N(j);
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_Q;
        t_Q(i, j) = t_kd(i, j) - t_P(i, j);
 
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_d_res;
        t_d_res(i, j) = -(bc.normalStiffness * t_P(i, j) +
                          bc.tangentialStiffness * t_Q(i, j));
 
        int rr = 0;
        for (; rr != row_nb_dofs / SPACE_DIM; ++rr) {
          auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
              locMat, SPACE_DIM * rr);
          auto t_col_base = col_data.getFTensor0N(gg, 0);
          for (auto cc = 0; cc != col_nb_dofs / SPACE_DIM; ++cc) {
            t_mat(i, j) += (t_w * t_row_base * t_col_base) * t_d_res(i, j);
            ++t_mat;
            ++t_col_base;
          }
          ++t_row_base;
        }
 
        for (; rr != nb_base_functions; ++rr)
          ++t_row_base;
 
        ++t_w;
        ++t_normal;
      }
 
      locMat *= OP::getMeasure();
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpringLhsBc_dP::iNtegrate(EntData &row_data,
                                           EntData &col_data) {
  MoFEMFunctionBegin;
 
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
  auto &locMat = OP::locMat;
  locMat.resize(row_nb_dofs, col_nb_dofs, false);
  locMat.clear();
 
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcSpringPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      auto t_normal = OP::getFTensor1NormalsAtGaussPts();
      auto t_w = OP::getFTensor0IntegrationWeight();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
      FTENSOR_INDEX(SPACE_DIM, k);
 
      int nb_integration_pts = OP::getGaussPts().size2();
      int nb_base_functions = row_data.getN().size2();
      auto t_row_base = row_data.getFTensor0N();
 
      auto t_kd = FTensor::Kronecker_Delta<double>();
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
        FTensor::Tensor1<double, SPACE_DIM> t_N;
        t_N(i) = t_normal(i);
        t_N.normalize();
 
        int rr = 0;
        for (; rr != row_nb_dofs / SPACE_DIM; ++rr) {
          auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
              OP::locMat, SPACE_DIM * rr);
          auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
          for (auto cc = 0; cc != col_nb_dofs / SPACE_DIM; ++cc) {
            t_mat(i, j) +=
                ((t_w * t_row_base) * (t_N(k) * t_col_base(k))) * 0.5 * t_kd(i, j);
            ++t_mat;
            ++t_col_base;
          }
          ++t_row_base;
        }
 
        for (; rr != nb_base_functions; ++rr)
          ++t_row_base;
 
        ++t_w;
        ++t_normal;
      }
 
      locMat *= OP::getMeasure();
    }
  }
  MoFEMFunctionReturn(0);
}
 
template <AssemblyType A>
MoFEMErrorCode OpAnalyticalDispBcImpl<A, GAUSS>::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  double time = OP::getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
 
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  // iterate over all boundary data
  for (auto &bc : (*bcDispPtr)) {
    // check if finite element entity is part of boundary condition
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      MatrixDouble analytical_expr;
      // placeholder to pass boundary block id to python
 
      auto [block_name, v_analytical_expr] =
          getAnalyticalExpr(this, analytical_expr, bc.blockName);
 
      int nb_dofs = data.getIndices().size();
      if (!nb_dofs) {
        continue;
      }
 
      int nb_integration_pts = OP::getGaussPts().size2();
      auto t_normal = OP::getFTensor1NormalsAtGaussPts();
      auto t_w = OP::getFTensor0IntegrationWeight();
      int nb_base_functions = data.getN().size2() / 3;
      auto t_row_base_fun = data.getFTensor1N<3>();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
 
      // get bc data
      auto t_bc_disp = getFTensor1FromMat<3, -1, DL>(v_analytical_expr);
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
 
        int bb = 0;
        for (; bb != nb_dofs / SPACE_DIM; ++bb) {
          t_nf(i) +=
              t_w * (t_row_base_fun(j) * t_normal(j)) * t_bc_disp(i) * 0.5;
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
 
        ++t_bc_disp;
        ++t_w;
        ++t_normal;
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpAnalyticalDispBc::iNtegrate(EntData &data) {
  return OP::iNtegrate(data);
}
 
MoFEMErrorCode OpBrokenTractionBc::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(3, i);
 
  int nb_dofs = data.getFieldData().size();
  int nb_integration_pts = getGaussPts().size2();
  int nb_base_functions = data.getN().size2();
 
  double time = getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
 
#ifndef NDEBUG
  if (this->locF.size() != nb_dofs)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
#endif // NDEBUG
 
  auto integrate_rhs = [&](auto &bc, auto calc_tau, double time_scale) {
    MoFEMFunctionBegin;
 
    auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
    auto t_row_base = data.getFTensor0N();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_normal = getFTensor1NormalsAtGaussPts();
 
    double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      double a = sqrt(t_normal(i) * t_normal(i));
      a /= 2.;
      const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
      auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
      int rr = 0;
      for (; rr != nb_dofs / SPACE_DIM; ++rr) {
        t_f(i) -=
            (time_scale * a * t_w * t_row_base * tau) * (t_val(i) * scale);
        ++t_row_base;
        ++t_f;
      }
 
      for (; rr != nb_base_functions; ++rr)
        ++t_row_base;
      ++t_w;
      ++t_coords;
      ++t_normal;
    }
    MoFEMFunctionReturn(0);
  };
 
  // get entity of face
  EntityHandle fe_ent = getFEEntityHandle();
  for (auto &bc : *(bcData)) {
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      double time_scale = 1;
      if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
        time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
      }
 
      int nb_dofs = data.getFieldData().size();
      if (nb_dofs) {
 
        if (std::regex_match(bc.blockName, std::regex(".*COOK.*"))) {
          auto calc_tau = [](double, double y, double) {
            y -= 44;
            y /= (60 - 44);
            return -y * (y - 1) / 0.25;
          };
          CHKERR integrate_rhs(bc, calc_tau, time_scale);
        } else {
          CHKERR integrate_rhs(
              bc, [](double, double, double) { return 1; }, time_scale);
        }
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpBrokenPressureBc::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(3, i);
 
  int nb_dofs = data.getFieldData().size();
  int nb_integration_pts = getGaussPts().size2();
  int nb_base_functions = data.getN().size2();
 
  double time = getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
 
#ifndef NDEBUG
  if (this->locF.size() != nb_dofs)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
#endif // NDEBUG
 
  auto integrate_rhs = [&](auto &bc, auto calc_tau, double time_scale) {
    MoFEMFunctionBegin;
 
    auto val = bc.val;
    auto t_row_base = data.getFTensor0N();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
    auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
 
    auto t_grad_gamma_u = getFTensor2FromMat<3, 2, -1, DL>(*hybridGradDispPtr);
 
    double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
      FTENSOR_INDEX(SPACE_DIM, k);
      FTensor::Number<0> N0;
      FTensor::Number<1> N1;
 
      FTensor::Tensor1<double, 3> t_normal;
      if (EshelbianCore::stretchSelector == LINEAR &&
          EshelbianCore::gradApproximator < MODERATE_ROT) {
 
        t_normal(i) = (FTensor::levi_civita<double>(i, j, k) * t_tangent1(j)) *
                      t_tangent2(k);
      } else {
        t_normal(i) = (FTensor::levi_civita<double>(i, j, k) *
                       (t_tangent1(j) + t_grad_gamma_u(j, N0))) *
                      (t_tangent2(k) + t_grad_gamma_u(k, N1));
      }
      auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
      auto t_val = FTensor::Tensor1<double, 3>();
      t_val(i) = (time_scale * t_w * tau * scale * val) * t_normal(i);
 
      auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
      int rr = 0;
      for (; rr != nb_dofs / SPACE_DIM; ++rr) {
        t_f(i) += t_row_base * t_val(i);
        ++t_row_base;
        ++t_f;
      }
 
      for (; rr != nb_base_functions; ++rr)
        ++t_row_base;
      ++t_w;
      ++t_coords;
      ++t_tangent1;
      ++t_tangent2;
      ++t_grad_gamma_u;
    }
    this->locF /= 2.;
 
    MoFEMFunctionReturn(0);
  };
 
  // get entity of face
  EntityHandle fe_ent = getFEEntityHandle();
  for (auto &bc : *(bcData)) {
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      double time_scale = 1;
      if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
        time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
      }
 
      int nb_dofs = data.getFieldData().size();
      if (nb_dofs) {
        CHKERR integrate_rhs(
            bc, [](double, double, double) { return 1; }, time_scale);
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
template <AssemblyType A>
MoFEMErrorCode
OpBrokenPressureBcLhsImpl_dU<A, GAUSS>::iNtegrate(EntData &row_data,
                                                  EntData &col_data) {
  MoFEMFunctionBegin;
 
  if (EshelbianCore::stretchSelector == LINEAR &&
      EshelbianCore::gradApproximator < MODERATE_ROT) {
    MoFEMFunctionReturnHot(0);
  }
 
  double time = OP::getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
 
  int nb_base_functions = row_data.getN().size2();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
  int nb_integration_pts = OP::getGaussPts().size2();
  auto &locMat = OP::locMat;
  locMat.resize(row_nb_dofs, col_nb_dofs, false);
  locMat.clear();
 
  auto integrate_lhs = [&](auto &bc, auto calc_tau, double time_scale) {
    MoFEMFunctionBegin;
 
    auto val = bc.val;
    auto t_row_base = row_data.getFTensor0N();
    auto t_w = OP::getFTensor0IntegrationWeight();
    auto t_coords = OP::getFTensor1CoordsAtGaussPts();
    auto t_tangent1 = OP::getFTensor1Tangent1AtGaussPts();
    auto t_tangent2 = OP::getFTensor1Tangent2AtGaussPts();
 
    auto t_grad_gamma_u = getFTensor2FromMat<3, 2, -1, DL>(*hybridGradDispPtr);
    auto t_kd = FTensor::Kronecker_Delta<double>();
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
      FTENSOR_INDEX(SPACE_DIM, k);
      FTENSOR_INDEX(SPACE_DIM, l);
 
      FTensor::Number<0> N0;
      FTensor::Number<1> N1;
 
      auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
      auto t_val = time_scale * t_w * tau * val;
 
      int rr = 0;
      for (; rr != row_nb_dofs / SPACE_DIM; ++rr) {
        auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
            locMat, SPACE_DIM * rr);
        auto t_diff_col_base = col_data.getFTensor1DiffN<2>(gg, 0);
        for (auto cc = 0; cc != col_nb_dofs / SPACE_DIM; ++cc) {
          FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_normal_du;
          t_normal_du(i, l) = (FTensor::levi_civita<double>(i, j, k) *
                               (t_tangent2(k) + t_grad_gamma_u(k, N1))) *
                                  t_kd(j, l) * t_diff_col_base(N0)
 
                              +
 
                              (FTensor::levi_civita<double>(i, j, k) *
                               (t_tangent1(j) + t_grad_gamma_u(j, N0))) *
                                  t_kd(k, l) * t_diff_col_base(N1);
 
          t_mat(i, j) += t_row_base * t_val * t_normal_du(i, j);
          ++t_mat;
          ++t_diff_col_base;
        }
        ++t_row_base;
      }
 
      for (; rr != nb_base_functions; ++rr)
        ++t_row_base;
      ++t_w;
      ++t_coords;
      ++t_tangent1;
      ++t_tangent2;
      ++t_grad_gamma_u;
    }
 
    OP::locMat /= 2.;
 
    MoFEMFunctionReturn(0);
  };
 
  // get entity of face
  EntityHandle fe_ent = OP::getFEEntityHandle();
  for (auto &bc : *(bcData)) {
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      double time_scale = 1;
      if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
        time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
      }
 
      int nb_dofs = row_data.getFieldData().size();
      if (nb_dofs) {
        CHKERR integrate_lhs(
            bc, [](double, double, double) { return 1; }, time_scale);
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpBrokenPressureBcLhs_dU::iNtegrate(EntData &row_data,
                                                   EntData &col_data) {
  return OP::iNtegrate(row_data, col_data);
}
 
MoFEMErrorCode OpBrokenAnalyticalTractionBc::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(3, i);
 
  int nb_dofs = data.getFieldData().size();
  int nb_integration_pts = getGaussPts().size2();
  int nb_base_functions = data.getN().size2();
 
  double time = getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
 
#ifndef NDEBUG
  if (this->locF.size() != nb_dofs)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
#endif // NDEBUG
 
  // get entity of face
  EntityHandle fe_ent = getFEEntityHandle();
  for (auto &bc : *(bcData)) {
    if (bc.faces.find(fe_ent) != bc.faces.end()) {
 
      MatrixDouble analytical_expr;
      // placeholder to pass boundary block id to python
      auto [block_name, v_analytical_expr] =
          getAnalyticalExpr(this, analytical_expr, bc.blockName);
      auto t_val = getFTensor1FromMat<3, -1, DL>(v_analytical_expr);
      auto t_row_base = data.getFTensor0N();
      auto t_w = getFTensor0IntegrationWeight();
      auto t_coords = getFTensor1CoordsAtGaussPts();
 
      double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
        auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
        int rr = 0;
        for (; rr != nb_dofs / SPACE_DIM; ++rr) {
          t_f(i) -= t_w * t_row_base * (t_val(i) * scale);
          ++t_row_base;
          ++t_f;
        }
 
        for (; rr != nb_base_functions; ++rr)
          ++t_row_base;
        ++t_w;
        ++t_coords;
        ++t_val;
      }
      this->locF *= getMeasure();
    }
  }
  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();
 
  auto piola_scale = dataAtPts->piolaScale;
  auto alpha_w = alphaW / piola_scale;
  auto alpha_rho = alphaRho / piola_scale;
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  double ts_scale = alpha_w * getTSa();
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
    ts_scale += alpha_rho * 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(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(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));
  };
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  auto t_approx_P_adjoint_log_du_dP =
      dataAtPts->getFTensorAdjointPdUdP(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
    int rr = 0;
    for (; rr != row_nb_dofs / 6; ++rr) {
 
      auto t_col_base_fun = col_data.getFTensor1N<3>(gg, 0);
      auto t_m = get_ftensor3(K, 6 * rr, 0);
 
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        FTensor::Tensor1<double, 3> t_col_base_piola;
        t_col_base_piola(j) =
            t_plasticF(j, k) * t_col_base_fun(k) / det_plasticF;
        t_m(L, i) -=
            a * (t_approx_P_adjoint_log_du_dP(i, j, L) *
                 t_col_base_piola(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_adjoint_log_du_dP;
    ++t_plasticF;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical_du_dBubble::integrate(EntData &row_data,
                                                       EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(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));
  };
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  int row_nb_base_functions = row_data.getN().size2();
 
  auto t_approx_P_adjoint_log_du_dP =
      dataAtPts->getFTensorAdjointPdUdP(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
    int rr = 0;
    for (; rr != row_nb_dofs / 6; ++rr) {
      auto t_m = get_ftensor2(K, 6 * rr, 0);
      auto t_col_base_fun = col_data.getFTensor2N<3, 3>(gg, 0);
      for (int cc = 0; cc != col_nb_dofs; ++cc) {
        FTensor::Tensor2<double, 3, 3> t_col_base_piola;
        t_col_base_piola(i, j) =
            t_col_base_fun(i, k) * t_plasticF(j, k) / det_plasticF;
        t_m(L) -=
            a * (t_approx_P_adjoint_log_du_dP(i, j, L) *
                 t_col_base_piola(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_adjoint_log_du_dP;
    ++t_plasticF;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialPhysical_du_domega::integrate(EntData &row_data,
                                                      EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
 
  int nb_integration_pts = getGaussPts().size2();
  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_adjoint_log_du_domega =
      dataAtPts->getFTensorAdjointPdUdOmega(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
  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) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
 
    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_adjoint_log_du_domega(k, L);
        ++t_m;
        ++t_col_base_fun;
      }
      ++t_row_base_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_approx_P_adjoint_log_du_domega;
    ++t_plasticF;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialRotation_domega_du::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 *, size_symm>, 3,
                            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)
 
    };
  };
 
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(size_symm, L);
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_levi_kirchhoff_du =
      dataAtPts->getFTensorLeviKirchhoffdLogStretch(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
  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) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
    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 / size_symm; ++cc) {
        t_m(k, L) -= (b * t_col_base_fun) * t_levi_kirchhoff_du(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_kirchhoff_du;
    ++t_plasticF;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialRotation_domega_dP::integrate(EntData &row_data,
                                                      EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, m);
  FTENSOR_INDEX(size_symm, L);
 
  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)
 
    );
  };
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
  auto t_levi_kirchhoff_dP =
      dataAtPts->getFTensorLeviKirchhoffP(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
    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) {
        FTensor::Tensor1<double, 3> t_col_base_piola;
        t_col_base_piola(k) =
            t_plasticF(k, l) * t_col_base_fun(l) / det_plasticF;
        t_m(m, i) -=
            b * (t_levi_kirchhoff_dP(m, i, k) * t_col_base_piola(k));
        ++t_m;
        ++t_col_base_fun;
      }
      ++t_row_base_fun;
    }
    for (; rr != row_nb_base_functions; ++rr) {
      ++t_row_base_fun;
    }
 
    ++t_w;
    ++t_levi_kirchhoff_dP;
    ++t_plasticF;
  }
  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(3, i);
  FTENSOR_INDEX(3, k);
  FTENSOR_INDEX(3, l);
  FTENSOR_INDEX(3, m);
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_levi_kirchoff_dP =
      dataAtPts->getFTensorLeviKirchhoffP(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
  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) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
    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) {
        FTensor::Tensor2<double, 3, 3> t_col_base_piola;
        t_col_base_piola(i, k) =
            t_col_base_fun(i, l) * t_plasticF(k, l) / det_plasticF;
        t_m(m) -=
            b * (t_levi_kirchoff_dP(m, i, k) * t_col_base_piola(i, k));
        ++t_m;
        ++t_col_base_fun;
      }
      ++t_row_base_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr) {
      ++t_row_base_fun;
    }
    ++t_w;
    ++t_levi_kirchoff_dP;
    ++t_plasticF;
  }
  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 t_kd = FTensor::Kronecker_Delta<double>();
 
  auto v = getVolume();
  auto ts_a =
      (std::abs(alphaViscousR) > std::numeric_limits<double>::epsilon() ||
       std::abs(alphaViscousOmega) > std::numeric_limits<double>::epsilon())
                  ? getTSa()
                  : 0.0;
  auto t_w = getFTensor0IntegrationWeight();
  auto t_levi_kirchhoff_domega =
      dataAtPts->getFTensorLeviKirchhoffdOmega(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
  auto t_invPlasticF = dataAtPts->getFTensorInvPlasticF(nb_integration_pts);
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
  auto t_row_grad_fun = row_data.getFTensor1DiffN<3>();
 
  // auto time_step = getTStimeStep();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
    double mass_coeff =
        a * alphaR + (a * alphaViscousR) * (ts_a /*/ time_step*/);
    double grad_coeff =
        a * alphaOmega + (a * alphaViscousOmega) * (ts_a /*/ time_step*/);
 
    int rr = 0;
    for (; rr != row_nb_dofs / 3; ++rr) {
      auto t_m = get_ftensor2(K, 3 * rr, 0);
      const double row_mass = a * t_row_base_fun;
      FTensor::Tensor1<double, 3> t_row_grad_intermediate;
      t_row_grad_intermediate(j) =
          t_row_grad_fun(i) * t_invPlasticF(i, j);
      auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
      auto t_col_grad_fun = col_data.getFTensor1DiffN<3>(gg, 0);
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        FTensor::Tensor1<double, 3> t_col_grad_intermediate;
        t_col_grad_intermediate(j) =
            t_col_grad_fun(i) * t_invPlasticF(i, j);
        t_m(k, l) -=
            (row_mass * t_col_base_fun) * t_levi_kirchhoff_domega(k, l);
        t_m(k, l) +=
            t_kd(k, l) * (mass_coeff * t_row_base_fun * t_col_base_fun);
        t_m(k, l) +=
            t_kd(k, l) *
            (grad_coeff *
             (t_row_grad_intermediate(j) * t_col_grad_intermediate(j)));
        ++t_m;
        ++t_col_base_fun;
        ++t_col_grad_fun;
      }
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
    for (; rr != row_nb_base_functions; ++rr) {
      ++t_row_base_fun;
      ++t_row_grad_fun;
    }
    ++t_w;
    ++t_levi_kirchhoff_domega;
    ++t_plasticF;
    ++t_invPlasticF;
  }
  MoFEMFunctionReturn(0);
}
 
template <typename TInvD, typename TRotation>
auto getDiffSpatialGradientDP(TInvD &t_d_u_d_b, TRotation &t_R) {
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l, m, n);
  constexpr auto t_diff_sym = FTensor::DiffSymmetrize<double>();
 
  FTensor::Tensor4<double, SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM>
      t_d_b_d_p;
  t_d_b_d_p(i, j, k, l) = t_diff_sym(i, j, m, l) * t_R(k, m);
 
  FTensor::Tensor4<double, SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM>
      t_d_u_d_p;
  t_d_u_d_p(i, j, k, l) =
      t_d_u_d_b(i, j, m, n) * t_d_b_d_p(m, n, k, l);
 
  FTensor::Tensor4<double, SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM>
      t_d_h_d_p;
  t_d_h_d_p(i, j, k, l) = t_R(i, m) * t_d_u_d_p(m, j, k, l);
  return t_d_h_d_p;
}
 
MoFEMErrorCode OpSpatialConsistency_dP_dP::integrate(EntData &row_data,
                                                     EntData &col_data) {
  MoFEMFunctionBeginHot;
  if (hasNonhomogeneousMatBlock ||
      dataAtPts->physicsPtr->getFeatures().test(
          PhysicalEquations::NO_STRETCH_NONLINEAR)) {
    MoFEMFunctionReturnHot(
        integrateImpl<size_symm * size_symm>(row_data, col_data));
  } else {
    MoFEMFunctionReturnHot(integrateImpl<0>(row_data, col_data));
  }
  MoFEMFunctionReturnHot(0);
};
 
template <int S>
MoFEMErrorCode OpSpatialConsistency_dP_dP::integrateImpl(EntData &row_data,
                                                         EntData &col_data) {
  MoFEMFunctionBegin;
 
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2)
 
    );
  };
 
  int nb_integration_pts = getGaussPts().size2();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  int row_nb_base_functions = row_data.getN().size2() / 3;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, m);
  FTENSOR_INDEX(SPACE_DIM, n);
 
  auto t_inv_D =
      getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(dataAtPts->matInvD);
  auto t_R = dataAtPts->getFTensorRotMat(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
  auto t_row_base = row_data.getFTensor1N<3>();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
 
    auto assemble = [&](auto &t_diff_h_p) {
      int rr = 0;
      for (; rr != row_nb_dofs / 3; ++rr) {
        FTensor::Tensor1<double, 3> t_row_base_piola;
        t_row_base_piola(j) =
            t_plasticF(j, m) * t_row_base(m) / det_plasticF;
        auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
        auto t_m = get_ftensor2(K, 3 * rr, 0);
        for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
          FTensor::Tensor1<double, 3> t_col_base_piola;
          t_col_base_piola(l) =
              t_plasticF(l, n) * t_col_base(n) / det_plasticF;
          t_m(i, k) -= a * t_row_base_piola(j) *
                       (t_diff_h_p(i, j, k, l) * t_col_base_piola(l));
          ++t_m;
          ++t_col_base;
        }
 
        ++t_row_base;
      }
 
      for (; rr != row_nb_base_functions; ++rr)
        ++t_row_base;
    };
 
    if (dataAtPts->physicsPtr->getFeatures().test(
            PhysicalEquations::NO_STRETCH_NONLINEAR)) {
      auto t_diff_h_p = getDiffSpatialGradientDP(t_inv_D, t_R);
      assemble(t_diff_h_p);
    } else {
      assemble(t_inv_D);
    }
 
    ++t_w;
    ++t_inv_D;
    ++t_R;
    ++t_plasticF;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
OpSpatialConsistency_dBubble_dBubble::integrate(EntData &row_data,
                                                EntData &col_data) {
  MoFEMFunctionBeginHot;
  if (hasNonhomogeneousMatBlock ||
      dataAtPts->physicsPtr->getFeatures().test(
          PhysicalEquations::NO_STRETCH_NONLINEAR)) {
    MoFEMFunctionReturnHot(
        integrateImpl<size_symm * size_symm>(row_data, col_data));
  } else {
    MoFEMFunctionReturnHot(integrateImpl<0>(row_data, col_data));
  }
  MoFEMFunctionReturnHot(0);
};
 
template <int S>
MoFEMErrorCode
OpSpatialConsistency_dBubble_dBubble::integrateImpl(EntData &row_data,
                                                    EntData &col_data) {
  MoFEMFunctionBegin;
 
  int nb_integration_pts = getGaussPts().size2();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  int row_nb_base_functions = row_data.getN().size2() / 9;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, m);
  FTENSOR_INDEX(SPACE_DIM, n);
 
  auto t_inv_D =
      getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(dataAtPts->matInvD);
  auto t_R = dataAtPts->getFTensorRotMat(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
  auto t_row_base = row_data.getFTensor2N<3, 3>();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
 
    auto assemble = [&](auto &t_diff_h_p) {
      int rr = 0;
      for (; rr != row_nb_dofs; ++rr) {
        FTensor::Tensor2<double, 3, 3> t_row_base_piola;
        t_row_base_piola(i, j) =
            t_row_base(i, m) * t_plasticF(j, m) / det_plasticF;
        auto t_col_base = col_data.getFTensor2N<3, 3>(gg, 0);
        for (int cc = 0; cc != col_nb_dofs; ++cc) {
          FTensor::Tensor2<double, 3, 3> t_col_base_piola;
          t_col_base_piola(k, l) =
              t_col_base(k, n) * t_plasticF(l, n) / det_plasticF;
          K(rr, cc) -=
              a * (t_row_base_piola(i, j) *
                   (t_diff_h_p(i, j, k, l) * t_col_base_piola(k, l)));
          ++t_col_base;
        }
 
        ++t_row_base;
      }
 
      for (; rr != row_nb_base_functions; ++rr)
        ++t_row_base;
    };
 
    if (dataAtPts->physicsPtr->getFeatures().test(
            PhysicalEquations::NO_STRETCH_NONLINEAR)) {
      auto t_diff_h_p = getDiffSpatialGradientDP(t_inv_D, t_R);
      assemble(t_diff_h_p);
    } else {
      assemble(t_inv_D);
    }
 
    ++t_w;
    ++t_inv_D;
    ++t_R;
    ++t_plasticF;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialConsistency_dBubble_dP::integrate(EntData &row_data,
                                                          EntData &col_data) {
  MoFEMFunctionBeginHot;
  if (hasNonhomogeneousMatBlock ||
      dataAtPts->physicsPtr->getFeatures().test(
          PhysicalEquations::NO_STRETCH_NONLINEAR)) {
    MoFEMFunctionReturnHot(
        integrateImpl<size_symm * size_symm>(row_data, col_data));
  } else {
    MoFEMFunctionReturnHot(integrateImpl<0>(row_data, col_data));
  }
  MoFEMFunctionReturnHot(0);
};
 
template <int S>
MoFEMErrorCode
OpSpatialConsistency_dBubble_dP::integrateImpl(EntData &row_data,
                                               EntData &col_data) {
  MoFEMFunctionBegin;
 
  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2)
 
    );
  };
 
  int nb_integration_pts = getGaussPts().size2();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  int row_nb_base_functions = row_data.getN().size2() / 9;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, m);
  FTENSOR_INDEX(SPACE_DIM, n);
 
  auto t_inv_D =
      getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(dataAtPts->matInvD);
  auto t_R = dataAtPts->getFTensorRotMat(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
  auto t_row_base = row_data.getFTensor2N<3, 3>();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
 
    auto assemble = [&](auto &t_diff_h_p) {
      auto t_m = get_ftensor1(K, 0, 0);
      int rr = 0;
      for (; rr != row_nb_dofs; ++rr) {
        FTensor::Tensor2<double, 3, 3> t_row_base_piola;
        t_row_base_piola(i, j) =
            t_row_base(i, m) * t_plasticF(j, m) / det_plasticF;
        auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
        for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
          FTensor::Tensor1<double, 3> t_col_base_piola;
          t_col_base_piola(l) =
              t_plasticF(l, n) * t_col_base(n) / det_plasticF;
          t_m(k) -=
              a * (t_row_base_piola(i, j) * t_diff_h_p(i, j, k, l)) *
              t_col_base_piola(l);
          ++t_col_base;
          ++t_m;
        }
 
        ++t_row_base;
      }
 
      for (; rr != row_nb_base_functions; ++rr)
        ++t_row_base;
    };
 
    if (dataAtPts->physicsPtr->getFeatures().test(
            PhysicalEquations::NO_STRETCH_NONLINEAR)) {
      auto t_diff_h_p = getDiffSpatialGradientDP(t_inv_D, t_R);
      assemble(t_diff_h_p);
    } else {
      assemble(t_inv_D);
    }
 
    ++t_w;
    ++t_inv_D;
    ++t_R;
    ++t_plasticF;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpStressGram_dBubble_dP::integrate(EntData &row_data,
                                                  EntData &col_data) {
  MoFEMFunctionBegin;
 
  auto get_ftensor1 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2)
 
    );
  };
 
  int nb_integration_pts = getGaussPts().size2();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  int row_nb_base_functions = row_data.getN().size2() / 9;
 
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
 
  auto t_row_base = row_data.getFTensor2N<3, 3>();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    auto t_m = get_ftensor1(K, 0, 0);
 
    int rr = 0;
    for (; rr != row_nb_dofs; ++rr) {
      auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
      for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
        t_m(k) += a * t_row_base(k, l) * t_col_base(l);
        ++t_col_base;
        ++t_m;
      }
 
      ++t_row_base;
    }
 
    for (; rr != row_nb_base_functions; ++rr)
      ++t_row_base;
    ++t_w;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpSpatialConsistency_dP_domega::integrate(EntData &row_data,
                                                         EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, m);
  FTENSOR_INDEX(SPACE_DIM, n);
 
  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)
 
    );
  };
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  int row_nb_base_functions = row_data.getN().size2() / 3;
  auto t_row_base_fun = row_data.getFTensor1N<3>();
 
  auto t_h_domega = dataAtPts->getFTensorSmallHdOmega(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
 
    int rr = 0;
    for (; rr != row_nb_dofs / 3; ++rr) {
 
      FTensor::Tensor1<double, 3> t_row_base_piola;
      t_row_base_piola(j) =
          t_plasticF(j, m) * t_row_base_fun(m) / det_plasticF;
      FTensor::Tensor2<double, 3, 3> t_PRT;
      t_PRT(i, k) = t_row_base_piola(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;
    ++t_plasticF;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
OpSpatialConsistency_dBubble_domega::integrate(EntData &row_data,
                                               EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, m);
  FTENSOR_INDEX(SPACE_DIM, n);
 
  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));
  };
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  int row_nb_base_functions = row_data.getN().size2() / 9;
  auto t_row_base_fun = row_data.getFTensor2N<3, 3>();
 
  auto t_h_domega = dataAtPts->getFTensorSmallHdOmega(nb_integration_pts);
  auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double det_plasticF = determinantTensor3by3(t_plasticF);
    double a = v * t_w * det_plasticF;
 
    int rr = 0;
    for (; rr != row_nb_dofs; ++rr) {
 
      FTensor::Tensor2<double, 3, 3> t_row_base_piola;
      t_row_base_piola(i, j) =
          t_row_base_fun(i, m) * t_plasticF(j, m) / det_plasticF;
      FTensor::Tensor1<double, 3> t_PRT;
      t_PRT(k) = t_row_base_piola(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;
    ++t_plasticF;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpPostProcDataStructure::doWork(int side, EntityType type,
                                               EntData &data) {
  MoFEMFunctionBegin;
 
  if (tagSense != getSkeletonSense())
    MoFEMFunctionReturnHot(0);
 
  auto create_tag = [this](const std::string tag_name, const int size) {
    double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
    Tag th;
    CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
                                       th, MB_TAG_CREAT | MB_TAG_SPARSE,
                                       def_VAL);
    return th;
  };
 
  Tag th_cauchy_streess = create_tag("CauchyStress", 9);
  Tag th_detF = create_tag("detF", 1);
  Tag th_traction = create_tag("traction", 3);
  Tag th_disp_error = create_tag("DisplacementError", 1);
 
  Tag th_energy = create_tag("Energy", 1);
  Tag th_young_modulus = create_tag("YoungModulus", 1);
 
  const auto nb_gauss_pts = getGaussPts().size2();
  auto t_w = dataAtPts->getFTensorSmallWL2(nb_gauss_pts);
  auto t_h = dataAtPts->getFTensorSmallH(nb_gauss_pts);
  auto t_approx_P = dataAtPts->getFTensorApproxP(nb_gauss_pts);
 
  auto t_normal = getFTensor1NormalsAtGaussPts();
  auto t_disp = dataAtPts->getFTensorSmallWH1(nb_gauss_pts);
 
  // auto sense = getSkeletonSense();
 
  if (dataAtPts->energyAtPts.size() == 0) {
    // that is for case that energy is not calculated
    dataAtPts->energyAtPts.resize(nb_gauss_pts);
    dataAtPts->energyAtPts.clear();
  }
  auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
  std::optional<decltype(getFTensor0FromVec(dataAtPts->youngModulusAtPts))>
      t_youngs_modulus;
  if (dataAtPts->physicsPtr->getFeatures().test(
          PhysicalEquations::NON_HOMOGENEOUS_MATERIAL)) {
    if (dataAtPts->youngModulusAtPts.size() != nb_gauss_pts)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Young's modulus postprocessing requires current material data");
    t_youngs_modulus.emplace(
        getFTensor0FromVec(dataAtPts->youngModulusAtPts));
  }
 
  auto next = [&]() {
    ++t_w;
    ++t_h;
    ++t_approx_P;
    ++t_normal;
    ++t_disp;
    if (t_youngs_modulus)
      ++*t_youngs_modulus;
    ++t_energy;
  };
 
  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);
  };
 
  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
    t_traction(i) *= tagSense;
    CHKERR save_vec_tag(th_traction, t_traction, gg);
 
    double u_error = sqrt((t_disp(i) - t_w(i)) * (t_disp(i) - t_w(i)));
    if (!std::isfinite(u_error))
      u_error = -1.;
    CHKERR save_scal_tag(th_disp_error, u_error, gg);
    CHKERR save_scal_tag(th_energy, t_energy, gg);
    if (t_youngs_modulus)
      CHKERR save_scal_tag(th_young_modulus, *t_youngs_modulus, 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_cauchy_streess, t_cauchy, gg);
    CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
 
    next();
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode AddHOOps<2, 3, 3>::add(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
    std::vector<FieldSpace> spaces, std::string geom_field_name,
    boost::shared_ptr<Range> crack_front_edges_ptr) {
  MoFEMFunctionBegin;
 
  constexpr bool scale_l2 = false;
 
  if (scale_l2) {
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
            "Scale L2 Ainsworth Legendre base is not implemented");
  }
 
  CHKERR MoFEM::AddHOOps<2, 3, 3>::add(pipeline, spaces, geom_field_name);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode AddHOOps<2, 2, 3>::add(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
    std::vector<FieldSpace> spaces, std::string geom_field_name,
    boost::shared_ptr<Range> crack_front_edges_ptr) {
  MoFEMFunctionBegin;
 
  constexpr bool scale_l2 = false;
 
  if (scale_l2) {
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
            "Scale L2 Ainsworth Legendre base is not implemented");
  }
 
  CHKERR MoFEM::AddHOOps<2, 2, 3>::add(pipeline, spaces, geom_field_name);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode AddHOOps<3, 3, 3>::add(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
    std::vector<FieldSpace> spaces, std::string geom_field_name,
    boost::shared_ptr<Range> crack_front_edges_ptr,
    boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
    boost::shared_ptr<MatrixDouble> inv_jac) {
  MoFEMFunctionBegin;
 
  if (!geom_field_name.empty()) {
    if (EshelbianCore::l2UserBaseScale) {
      auto jac = boost::make_shared<MatrixDouble>();
      auto det = boost::make_shared<VectorDouble>();
      pipeline.push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
              geom_field_name, jac));
      pipeline.push_back(new OpInvertMatrix<3>(jac, det, nullptr));
      pipeline.push_back(
          new OpScaleBaseBySpaceInverseOfMeasure(L2, USER_BASE, det));
    }
  }
 
  constexpr bool scale_l2_ainsworth_legendre_base = false;
 
  if (scale_l2_ainsworth_legendre_base) {
 
    struct OpCalculateVectorFieldGradient
        : public MoFEM::OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM> {
 
      using OP = MoFEM::OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>;
 
      OpCalculateVectorFieldGradient(const std::string &field_name,
                                     boost::shared_ptr<MatrixDouble> jac,
                                     boost::shared_ptr<Range> edges_ptr)
          : OP(field_name, jac), edgesPtr(edges_ptr) {}
 
      MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
 
        auto ent = data.getFieldEntities().size()
                       ? data.getFieldEntities()[0]->getEnt()
                       : 0;
 
        if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
          return 0;
        } else {
          return OP::doWork(side, type, data);
        }
      };
 
    private:
      boost::shared_ptr<Range> edgesPtr;
    };
 
    if (!geom_field_name.empty()) {
      auto jac = boost::make_shared<MatrixDouble>();
      auto det = boost::make_shared<VectorDouble>();
      pipeline.push_back(new OpCalculateVectorFieldGradient(
          geom_field_name, jac,
          EshelbianCore::setSingularity ? crack_front_edges_ptr
                                        : boost::make_shared<Range>()));
      pipeline.push_back(new OpInvertMatrix<3>(jac, det, nullptr));
      pipeline.push_back(new OpScaleBaseBySpaceInverseOfMeasure(
          L2, AINSWORTH_LEGENDRE_BASE, det));
    }
  }
 
  CHKERR MoFEM::AddHOOps<3, 3, 3>::add(pipeline, spaces, geom_field_name, jac,
                                       det, inv_jac);
 
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Caluclate face material force and normal pressure at gauss points
 *
 * @param side
 * @param type
 * @param data
 * @return MoFEMErrorCode
 *
 * Reconstruct the full gradient \f$U=\nabla u\f$ on a surface from the
 * symmetric part and the surface gradient.
 *
 * @details
 * Inputs:
 *  - \c t_strain : \f$\varepsilon=\tfrac12(U+U^\top)\f$ (symmetric strain on
 * S),
 *  - \c t_grad_u_gamma : \f$u^\Gamma = U P\f$ (right-projected/surface
 * gradient), with \f$P=I-\mathbf N\otimes\mathbf N\f$,
 *  - \c t_normal : (possibly non‑unit) surface normal.
 *
 * Procedure (pointwise on S):
 *  1) Normalize the normal \f$\mathbf n=\mathbf N/\|\mathbf N\|\f$.
 *  2) Form the residual \f$R=\varepsilon-\operatorname{sym}(u^\Gamma)\f$, where
 *     \f$\operatorname{sym}(A)=\tfrac12(A+A^\top)\f$.
 *  3) Recover the normal directional derivative (a vector)
 *     \f$\mathbf v=\partial_{\mathbf n}u=2R\mathbf n-(\mathbf n^\top R\,\mathbf
 * n)\,\mathbf n\f$. 4) Assemble the full gradient \f$U = u^\Gamma + \mathbf
 * v\otimes \mathbf n\f$.
 *
 * Properties (sanity checks):
 *  - \f$\tfrac12(U+U^\top)=\varepsilon\f$ (matches the given symmetric part),
 *  - \f$U P = u^\Gamma\f$ (tangential/right-projected columns unchanged),
 *  - Only the **normal column** is updated via \f$\mathbf v\otimes\mathbf n\f$.
 *
 * Mapping to variables in this snippet:
 *  - \f$\varepsilon \leftrightarrow\f$ \c t_strain,
 *  - \f$u^\Gamma \leftrightarrow\f$ \c t_grad_u_gamma,
 *  - \f$\mathbf N \leftrightarrow\f$ \c t_normal (normalized into \c t_N),
 *  - \f$R \leftrightarrow\f$ \c t_R,
 *  - \f$U \leftrightarrow\f$ \c t_grad_u.
 *
 * @pre \c t_normal is nonzero; \c t_strain is symmetric.
 * @note All indices use Einstein summation; computation is local to the surface
 * point.
 *
 */
MoFEMErrorCode OpFaceSideMaterialForce::doWork(int side, EntityType type,
                                               EntData &data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, m);
  FTENSOR_INDEX(SPACE_DIM, n);
  FTENSOR_INDEX(SPACE_DIM, J);
  FTENSOR_INDEX(SPACE_DIM, I);
  FTENSOR_INDEX(SPACE_DIM, K);
  FTENSOR_INDEX(SPACE_DIM, M);
  FTENSOR_INDEX(SPACE_DIM, N);
  FTENSOR_INDEX(SPACE_DIM, L);
 
  const auto nb_gauss_pts = getGaussPts().size2();
  MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
      dataAtPts->faceMaterialForceAtPts, nb_gauss_pts);
  dataAtPts->normalPressureAtPts.resize(nb_gauss_pts, false);
  if (getNinTheLoop() == 0) {
    dataAtPts->faceMaterialForceAtPts.clear();
    dataAtPts->normalPressureAtPts.clear();
  }
  auto loop_size = getLoopSize();
  if (loop_size == 1) {
    auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
    auto pstatus = numebered_fe_ptr->getPStatus();
    if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
      loop_size = 2;
    }
  }
 
  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<double>();
 
  auto t_normal = getFTensor1NormalsAtGaussPts();
  auto t_T = dataAtPts->getFTensorFaceMaterialForce(
      nb_gauss_pts); //< face material force
  auto t_p =
      getFTensor0FromVec(dataAtPts->normalPressureAtPts); //< normal pressure
  auto t_P = dataAtPts->getFTensorApproxP(nb_gauss_pts);
  auto t_u_gamma = dataAtPts->getFTensorSmallHybridDisp(nb_gauss_pts);
  auto t_grad_u_gamma = dataAtPts->getFTensorGradHybridDisp(nb_gauss_pts);
  auto t_strain = dataAtPts->getFTensorLogStretch(nb_gauss_pts);
  auto t_omega = dataAtPts->getFTensorRotAxis(nb_gauss_pts);
 
  FTensor::Tensor1<double, SPACE_DIM> t_N;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_A;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_R;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_grad_u;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_Proj;
 
  auto next = [&]() {
    ++t_normal;
    ++t_P;
    // ++t_grad_P;
    ++t_omega;
    ++t_u_gamma;
    ++t_grad_u_gamma;
    ++t_strain;
    ++t_T;
    ++t_p;
  };
 
  switch (EshelbianCore::energyReleaseSelector) {
  case GRIFFITH_FORCE:
    for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
      t_N(I) = t_normal(I);
      t_N.normalize();
 
      t_A(i, j) = levi_civita(i, j, k) * t_omega(k);
      t_R(i, k) = t_kd(i, k) + t_A(i, k);
      t_grad_u(i, j) = t_R(i, j) + t_strain(i, j);
 
      t_T(I) += t_N(J) * (t_grad_u(i, I) * t_P(i, J)) / loop_size;
      // note that works only for Hooke material, for nonlinear material we need
      // strain energy expressed by stress
      t_T(I) -= t_N(I) * ((t_strain(i, K) * t_P(i, K)) / 2.) / loop_size;
 
      t_p += t_N(I) *
             (t_N(J) * ((t_kd(i, I) + t_grad_u_gamma(i, I)) * t_P(i, J))) /
             loop_size;
 
      next();
    }
    break;
  case GRIFFITH_SKELETON:
    for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
 
      // Normalize the normal
      t_N(I) = t_normal(I);
      t_N.normalize();
 
      // R = ε − sym(u^Γ)
      t_R(i, j) =
          t_strain(i, j) - 0.5 * (t_grad_u_gamma(i, j) + t_grad_u_gamma(j, i));
 
      // U = u^Γ + [2 R N − (Nᵀ R N) N] ⊗ N
      t_grad_u(i, J) =
          t_grad_u_gamma(i, J) +
          (2 * t_R(i, K) * t_N(K) - (t_R(k, L) * t_N(k) * t_N(L)) * t_N(i)) *
              t_N(J);
 
      t_T(I) += t_N(J) * (t_grad_u(i, I) * t_P(i, J)) / loop_size;
      // note that works only for Hooke material, for nonlinear material we need
      // strain energy expressed by stress
      t_T(I) -= t_N(I) * ((t_strain(i, K) * t_P(i, K)) / 2.) / loop_size;
 
      // calculate nominal face pressure
      t_p += t_N(I) *
             (t_N(J) * ((t_kd(i, I) + t_grad_u_gamma(i, I)) * t_P(i, J))) /
             loop_size;
 
      next();
    }
    break;
 
  default:
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
            "Grffith energy release "
            "selector not implemented");
  };
 
#ifndef NDEBUG
  auto side_fe_ptr = getSidePtrFE();
  auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
  auto pstatus = side_fe_mi_ptr->getPStatus();
  if (pstatus) {
    auto owner = side_fe_mi_ptr->getOwnerProc();
    MOFEM_LOG("SELF", Sev::noisy)
        << "OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
        << " " << getPtrFE()->mField.get_comm_rank() << " n in the loop "
        << getNinTheLoop() << " loop size " << getLoopSize();
  }
#endif // NDEBUG
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpFaceMaterialForce::doWork(int side, EntityType type,
                                           EntData &data) {
  MoFEMFunctionBegin;
 
#ifndef NDEBUG
  auto fe_mi_ptr = getFEMethod()->numeredEntFiniteElementPtr;
  auto pstatus = fe_mi_ptr->getPStatus();
  if (pstatus) {
    auto owner = fe_mi_ptr->getOwnerProc();
    MOFEM_LOG("SELF", Sev::noisy)
        << "OpFaceMaterialForce: owner proc is not 0, owner proc: " << owner
        << " " << getPtrFE()->mField.get_comm_rank();
  }
#endif // NDEBUG
 
  FTENSOR_INDEX(SPACE_DIM, I);
 
  FTensor::Tensor1<double, SPACE_DIM> t_face_T;
  t_face_T(I) = 0.;
  double face_pressure = 0.;
  auto t_T = dataAtPts->getFTensorFaceMaterialForce(
      getGaussPts().size2()); //< face material force
  auto t_p =
      getFTensor0FromVec(dataAtPts->normalPressureAtPts); //< normal pressure
  auto t_w = getFTensor0IntegrationWeight();
  for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
    t_face_T(I) += t_w * t_T(I);
    face_pressure += t_w * t_p;
    ++t_w;
    ++t_T;
    ++t_p;
  }
  t_face_T(I) *= getMeasure();
  face_pressure *= getMeasure();
 
  auto get_tag = [&](auto name, auto dim) {
    auto &moab = getPtrFE()->mField.get_moab();
    Tag tag;
    double def_val[] = {0., 0., 0.};
    CHK_MOAB_THROW(moab.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
                                       MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
                   "create tag");
    return tag;
  };
 
  auto set_tag = [&](auto &&tag, auto ptr) {
    auto &moab = getPtrFE()->mField.get_moab();
    auto face = getPtrFE()->getFEEntityHandle();
    CHK_MOAB_THROW(moab.tag_set_data(tag, &face, 1, ptr), "set tag");
  };
 
  set_tag(get_tag("MaterialForce", 3), &t_face_T(0));
  set_tag(get_tag("FacePressure", 1), &face_pressure);
 
  MoFEMFunctionReturn(0);
}
 
template <typename OP_PTR>
std::tuple<std::string, MatrixDouble>
getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr,
                  const std::string block_name) {
 
  auto nb_gauss_pts = op_ptr->getGaussPts().size2();
 
  auto ts_time = op_ptr->getTStime();
  auto ts_time_step = op_ptr->getTStimeStep();
 
  if (EshelbianCore::physicalTimeFlg) {
    ts_time = EshelbianCore::currentPhysicalTime;
    ts_time_step = EshelbianCore::physicalDt;
  }
 
  MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
  MatrixDouble m_ref_normals = op_ptr->getNormalsAtGaussPts();
 
  auto v_analytical_expr =
      analytical_expr_function(ts_time_step, ts_time, nb_gauss_pts,
                               m_ref_coords, m_ref_normals, block_name);
 
  if (PetscUnlikely(!v_analytical_expr.size2())) {
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Analytical expression is empty or does not exist, "
                      "check python file");
  }
 
  return std::make_tuple(block_name, v_analytical_expr);
}
 
OpBrokenBaseTimesHybridDisp::OpBrokenBaseTimesHybridDisp(
    boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
    boost::shared_ptr<MatrixDouble> vec, ScalarFun beta_coeff,
    boost::shared_ptr<Range> ents_ptr)
    : OP(broken_base_side_data, ents_ptr) {
  this->sourceVec = vec;
  this->betaCoeff = beta_coeff;
}
 
OpBrokenBaseTimesBrokenDisp::OpBrokenBaseTimesBrokenDisp(
    boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
    ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
    : OP(broken_base_side_data, ents_ptr) {
  this->sourceVec = boost::shared_ptr<MatrixDouble>();
  this->betaCoeff = beta_coeff;
}
 
MoFEMErrorCode
OpBrokenBaseTimesBrokenDisp::doWork(int row_side, EntityType row_type,
                                    EntitiesFieldData::EntData &row_data) {
  MoFEMFunctionBegin;
 
  if (OP::entsPtr) {
    if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
      MoFEMFunctionReturnHot(0);
  }
 
#ifndef NDEBUG
  if (!brokenBaseSideData) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not set");
  }
#endif // NDEBUG
 
  auto do_work_rhs = [this](int row_side, EntityType row_type,
                            EntitiesFieldData::EntData &row_data) {
    MoFEMFunctionBegin;
    // get number of dofs on row
    OP::nbRows = row_data.getIndices().size();
    if (!OP::nbRows)
      MoFEMFunctionReturnHot(0);
    // get number of integration points
    OP::nbIntegrationPts = OP::getGaussPts().size2();
    // get row base functions
    OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
    // resize and clear the right hand side vector
    OP::locF.resize(OP::nbRows, false);
    OP::locF.clear();
    // integrate local vector
    CHKERR this->iNtegrate(row_data);
    // assemble local vector
    CHKERR this->aSsemble(row_data);
    MoFEMFunctionReturn(0);
  };
 
  switch (OP::opType) {
  case OP::OPSPACE:
    for (auto &bd : *brokenBaseSideData) {
      this->sourceVec =
          boost::shared_ptr<MatrixDouble>(brokenBaseSideData, &bd.getFlux());
      CHKERR do_work_rhs(bd.getSide(), bd.getType(), bd.getData());
      this->sourceVec.reset();
    }
    break;
  default:
    CHK_MOAB_THROW(MOFEM_IMPOSSIBLE_CASE,
                   (std::string("wrong op type ") +
                    OpBaseDerivativesBase::OpTypeNames[OP::opType])
                       .c_str());
  }
 
  MoFEMFunctionReturn(0);
}
 
OpHybridBaseTimesBrokenDisp::OpHybridBaseTimesBrokenDisp(
    const std::string row_field,
    boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
    ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
    : OP(row_field, boost::shared_ptr<MatrixDouble>(), beta_coeff, ents_ptr),
      brokenBaseSideDataPtr(broken_base_side_data) {
  this->betaCoeff = beta_coeff;
}
 
MoFEMErrorCode
OpHybridBaseTimesBrokenDisp::iNtegrate(EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
  for (auto &bd : (*brokenBaseSideDataPtr)) {
    this->sourceVec =
        boost::shared_ptr<MatrixDouble>(brokenBaseSideDataPtr, &bd.getFlux());
 
#ifndef NDEBUG
    if (this->sourceVec->size2() != SPACE_DIM) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Inconsistent size of the source vector");
    }
    if (this->sourceVec->size1() != OP::getGaussPts().size2()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Inconsistent size of the source vector");
    }
#endif // NDEBUG
 
    CHKERR OP::iNtegrate(data);
 
    this->sourceVec.reset();
  }
  MoFEMFunctionReturn(0);
}
 
OpHyrbridBaseBrokenBase::OpHyrbridBaseBrokenBase(
    std::string row_field,
    boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
    ScalarFun beta, const bool assmb_transpose, const bool only_transpose,
    boost::shared_ptr<Range> ents_ptr)
    : OP(row_field, broken_base_side_data, assmb_transpose, only_transpose,
         ents_ptr) {
  this->betaCoeff = beta;
  this->sYmm = false;
}
 
OpBrokenBaseBrokenBase::OpBrokenBaseBrokenBase(
    boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
    ScalarFun beta, boost::shared_ptr<Range> ents_ptr)
    : OP(broken_base_side_data, ents_ptr) {
  this->sYmm = false;
  this->betaCoeff = beta;
  OP::assembleTranspose = false;
  OP::onlyTranspose = false;
}
 
MoFEMErrorCode
OpBrokenBaseBrokenBase::doWork(int row_side, EntityType row_type,
                               EntitiesFieldData::EntData &row_data) {
  MoFEMFunctionBegin;
 
  if (OP::entsPtr) {
    if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
      MoFEMFunctionReturnHot(0);
  }
 
#ifndef NDEBUG
  if (!brokenBaseSideData) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not set");
  }
#endif // NDEBUG
 
  auto do_work_lhs = [this](int row_side, int col_side, EntityType row_type,
                            EntityType col_type,
                            EntitiesFieldData::EntData &row_data,
                            EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
 
    auto check_if_assemble_transpose = [&] {
      if (this->sYmm) {
        if (OP::rowSide != OP::colSide || OP::rowType != OP::colType)
          return true;
        else
          return false;
      } else if (OP::assembleTranspose) {
        return true;
      }
      return false;
    };
 
    OP::rowSide = row_side;
    OP::rowType = row_type;
    OP::colSide = col_side;
    OP::colType = col_type;
    OP::nbCols = col_data.getIndices().size();
    OP::locMat.resize(OP::nbRows, OP::nbCols, false);
    OP::locMat.clear();
    CHKERR this->iNtegrate(row_data, col_data);
    CHKERR this->aSsemble(row_data, col_data, check_if_assemble_transpose());
    MoFEMFunctionReturn(0);
  };
 
  switch (OP::opType) {
  case OP::OPSPACE:
 
    for (auto &bd : *brokenBaseSideData) {
 
#ifndef NDEBUG
      if (!bd.getData().getNSharedPtr(bd.getData().getBase())) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "base functions not set");
      }
#endif
 
      OP::nbRows = bd.getData().getIndices().size();
      if (!OP::nbRows)
        MoFEMFunctionReturnHot(0);
      OP::nbIntegrationPts = OP::getGaussPts().size2();
      OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(bd.getData());
 
      if (!OP::nbRows)
        MoFEMFunctionReturnHot(0);
 
      CHKERR do_work_lhs(
 
          // side
          bd.getSide(), bd.getSide(),
 
          // type
          bd.getType(), bd.getType(),
 
          // row_data
          bd.getData(), bd.getData()
 
      );
    }
 
    break;
 
  default:
    CHK_MOAB_THROW(MOFEM_IMPOSSIBLE_CASE,
                   (std::string("wrong op type ") +
                    OpBaseDerivativesBase::OpTypeNames[OP::opType])
                       .c_str());
  }
 
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity