EshelbianCohesive.cpp

Go to the documentation of this file.

/**
 * @file EshelbianCohesive.cpp
 * @author your name (you@domain.com)
 * @brief
 * @version 0.1
 * @date 2025-12-21
 *
 * @copyright Copyright (c) 2025
 *
 */
 
#define SINGULARITY
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
#include <TSElasticPostStep.hpp>
 
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/lambert_w.hpp>
 
#include <EshelbianCohesive.hpp>
 
namespace EshelbianPlasticity {
 
struct GriffithCohesiveLaw {
 
  static inline constexpr double eps = 1e-8;
 
  static inline constexpr double pi = boost::math::constants::pi<double>();
  static inline constexpr double root_pi =
      boost::math::constants::root_pi<double>();
 
  template <typename T> static inline T getTau(const T &k, double gc) {
    const T c = static_cast<T>(gc) / root_pi;
    return c * std::exp(-k) / std::sqrt(k);
  }
 
  template <typename T>
  static inline auto getDiffTau(const T &tau, const T &kappa, double gc) {
    return -tau * (1.0 + (1.0 / (2.0 * kappa)));
  }
 
  template <typename T>
  static T getAlpha(const T &kappa, double gc, double min_stiffness) {
    auto tau = getTau(kappa, gc);
    return kappa / (kappa * min_stiffness + tau);
  }
 
  template <typename T>
  static T getDiffAlpha(const T &kappa, double gc, double min_stiffness) {
    T tau = getTau(kappa, gc);
    T diff_tau = getDiffTau(tau, kappa, gc);
    T dnom = kappa * min_stiffness + tau;
    return (tau - kappa * diff_tau) / (dnom * dnom);
  }
 
  template <typename T> static auto invTau(const T &tau, double Gf) {
    const T z = (2.0 * Gf * Gf) / (pi * tau * tau);
    return T(0.5) * boost::math::lambert_w0(z); // k=0 == principal branch
  }
 
  static double getInitialStrength(double kappa, double gc) {
    const double tau = getTau(kappa, gc);
    return tau * std::sqrt(2.0 / (kappa + 1.5));
  }
 
  static double invInitialStrength(double strength, double gc) {
    // For negligible residual stiffness, the Mode-I damage-onset condition
    //
    //   0.5 * strength^2 * alpha'(kappa) = tau(kappa)
    //
    // gives strength = tau(kappa) * sqrt(2 / (kappa + 3/2)). This
    // function is strictly decreasing, so use invTau(strength) as an
    // initial upper bound and enlarge it until the root is bracketed.
    double lower_kappa = 0;
    double upper_kappa = invTau(strength, gc);
    while (getInitialStrength(upper_kappa, gc) > strength)
      upper_kappa *= 2;
 
    constexpr int max_bisection_iterations = 100;
    for (int i = 0; i != max_bisection_iterations; ++i) {
      const double mid_kappa = 0.5 * (lower_kappa + upper_kappa);
      if (getInitialStrength(mid_kappa, gc) > strength)
        lower_kappa = mid_kappa;
      else
        upper_kappa = mid_kappa;
    }
 
    return 0.5 * (lower_kappa + upper_kappa);
  }
 
  template <typename T>
  static auto calculateY(const T t_eff, const T &kappa, double gc,
                         double min_stiffness) {
    T diff_alpha = getDiffAlpha(kappa, gc, min_stiffness);
    return 0.5 * t_eff * t_eff * diff_alpha;
  }
 
  template <typename T>
  static auto calculateDissipation(const T &delta_kappa, const T t_eff,
                                   const T &kappa, double gc,
                                   double min_stiffness) {
    T Y = calculateY(t_eff, kappa + delta_kappa, gc, min_stiffness);
    return Y * delta_kappa;
  }
 
  template <typename T>
  static auto calculateDissipationSurplus(const T &delta_kappa, const T t_eff,
                                          const T &kappa, double gc,
                                          double min_stiffness) {
    T M = calculateY(t_eff, kappa + delta_kappa, gc, min_stiffness) -
          getTau(kappa + delta_kappa, gc);
    return -M * delta_kappa;
  }
 
  template <typename T>
  static auto calculateDissipationSurplusDiffKappa(const T &delta_kappa,
                                                   const T &t_eff,
                                                   const T &kappa, double gc,
                                                   double min_stiffness) {
 
    std::complex<T> cpx_delta = delta_kappa;
    std::complex<T> cpx_t_eff = t_eff;
    std::complex<T> cpx_kappa = kappa;
    cpx_delta += eps * 1i;
    std::complex<T> cpx_M = calculateDissipationSurplus(
        cpx_delta, cpx_t_eff, cpx_kappa, gc, min_stiffness);
    return cpx_M.imag() / eps;
  }
 
  template <typename T>
  static auto calculateDissipationSurplusDiffTraction(
      const FTensor::Tensor1<T, 3> &t_traction, const T &delta_kappa,
      const T &kappa, FTensor::Tensor1<double, 3> &t_n_normalize, double gc,
      double beta, double min_stiffness) {
    FTENSOR_INDEX(3, i);
    FTensor::Tensor1<double, 3> t_diff_traction;
    std::complex<T> cpx_delta = delta_kappa;
    std::complex<T> cpx_kappa = kappa;
    FTensor::Tensor1<std::complex<T>, 3> t_cpx_traction;
    for (auto jj = 0; jj != 3; ++jj) {
      t_cpx_traction(i) = t_traction(i);
      t_cpx_traction(jj) += eps * 1i;
      std::complex<T> cpx_teff =
          calculateEffectiveTraction(t_cpx_traction, t_n_normalize, beta);
      std::complex<T> cpx_M = calculateDissipationSurplus(
          cpx_delta, cpx_teff, cpx_kappa, gc, min_stiffness);
      t_diff_traction(jj) = cpx_M.imag() / eps;
    }
    return t_diff_traction;
  }
 
  template <typename T>
  static auto calculateGap(const FTensor::Tensor1<T, 3> &t_traction,
                           FTensor::Tensor1<double, 3> &t_n_normalize,
                           double alpha, double beta,
                           bool sign_sensitive = true) {
    FTENSOR_INDEX(3, i);
    FTENSOR_INDEX(3, j);
    constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
    FTensor::Tensor2<T, 3, 3> t_P;
    t_P(i, j) = t_n_normalize(i) * t_n_normalize(j);
    FTensor::Tensor2<T, 3, 3> t_Q;
    t_Q(i, j) = t_kd(i, j) - t_P(i, j);
    FTensor::Tensor1<T, 3> t_normal;
    t_normal(i) = t_P(i, j) * t_traction(j);
    FTensor::Tensor1<T, 3> t_tangential;
    t_tangential(i) = t_Q(i, j) * t_traction(j);
 
    if (sign_sensitive) {
      T s = std::sqrt((t_normal(i) * t_normal(i)) / 4. +
                      (1. / beta) * t_tangential(i) * t_tangential(i));
      T teff = (t_n_normalize(i) * t_normal(i) / 2.) + s;
 
      FTensor::Tensor1<T, 3> t_gap{0., 0., 0.};
      if (std::real(s) > std::numeric_limits<double>::epsilon()) {
        t_gap(i) = alpha * (
 
                               t_n_normalize(i) * teff / 2. +
 
                               (1.0 / 4.0) * (teff / s) * t_normal(i) +
 
                               (1.0 / beta) * (teff / s) * t_tangential(i)
 
                           );
      }
      return std::make_pair(teff, t_gap);
    } else {
      T teff = std::sqrt((t_normal(i) * t_normal(i)) +
                         (1. / beta) * t_tangential(i) * t_tangential(i));
      FTensor::Tensor1<T, 3> t_gap;
      t_gap(i) = alpha * (
 
                             t_normal(i) + (1.0 / beta) * t_tangential(i)
 
                         );
      return std::make_pair(teff, t_gap);
    }
  }
 
  template <typename T>
  static auto
  calculateDiffGapDTraction(const FTensor::Tensor1<T, 3> &t_traction,
                            FTensor::Tensor1<double, 3> &t_n_normalize,
                            double alpha, double beta,
                            bool sign_sensitive = true) {
    FTENSOR_INDEX(3, i);
    FTensor::Tensor2<T, 3, 3> t_dgap;
    FTensor::Tensor1<std::complex<T>, 3> t_cpx_delta;
    FTensor::Tensor1<std::complex<T>, 3> t_cpx_traction;
    for (auto jj = 0; jj != 3; ++jj) {
      t_cpx_traction(i) = t_traction(i);
      t_cpx_traction(jj) += eps * 1i;
      auto [teff_cpx, t_cpx_gap] = calculateGap(t_cpx_traction, t_n_normalize,
                                                alpha, beta, sign_sensitive);
      for (auto ii = 0; ii != 3; ++ii) {
        auto v = t_cpx_gap(ii).imag();
        t_dgap(ii, jj) = v / eps;
      }
    }
    return t_dgap;
  }
 
  template <typename T>
  static auto
  calculateEffectiveTraction(const FTensor::Tensor1<T, 3> &t_traction,
                             FTensor::Tensor1<double, 3> &t_n_normalize,
                             double beta) {
    auto [teff, t_gap] = calculateGap(t_traction, t_n_normalize, 1.0, beta);
    return teff;
  }
};  
 
struct OpGetParameters : FaceElementForcesAndSourcesCore::UserDataOperator {
  OpGetParameters(boost::shared_ptr<double> gc_ptr, Sev severity = Sev::inform)
      : FaceElementForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE),
        gcPtr(gc_ptr), logSev(severity) {
    CHK_THROW_MESSAGE(getOptions(), "Failed to get EshelbianCohesive options");
  }
 
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        EntitiesFieldData::EntData &row_data) {
    MoFEMFunctionBegin;
    *gcPtr = defaultGc;
    MoFEMFunctionReturn(0);
  }
 
  static inline double strength = -1;
  static inline double min_kappa = 1e-12;
  static inline double min_stiffness = 1e-2;
  static inline double kappa0 = 1;
  static inline double beta = 1;
 
private:
  MoFEMErrorCode getOptions() {
    MoFEMFunctionBegin;
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "interface_", "", "none");
 
    CHKERR PetscOptionsScalar("-gc", "Griffith energy release rate", "",
                              defaultGc, &defaultGc, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-min_stiffness", "Minimal interface stiffness",
                              "", min_stiffness, &min_stiffness, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-strength", "Strength of interface", "",
                              strength, &strength, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-min_kappa",
                              "Minimal kappa to avoid singularity", "",
                              min_kappa, &min_kappa, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-kappa0", "Characteristic length kappa0", "",
                              kappa0, &kappa0, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-beta", "Cohesive tangential coupling", "",
                              beta, &beta, PETSC_NULLPTR);
 
    PetscOptionsEnd();
 
    MOFEM_LOG("EP", logSev)
        << "Interface Griffith energy release rate Gc -interface_gc = "
        << defaultGc;
    MOFEM_LOG("EP", logSev)
        << "Interface min stiffness -interface_min_stiffness = "
        << min_stiffness;
    MOFEM_LOG("EP", logSev)
        << "Interface strength -interface_strength = " << strength;
    MOFEM_LOG("EP", logSev)
        << "Interface minimal kappa -interface_min_kappa = " << min_kappa;
    MOFEM_LOG("EP", logSev)
        << "Interface characteristic length  kappa0 -interface_kappa0 = "
        << kappa0;
    MOFEM_LOG("EP", logSev)
        << "Interface tangential coupling -interface_beta = " << beta;
 
    if (strength > 0) {
      const double kappa_strength = GriffithCohesiveLaw::invInitialStrength(
          static_cast<double>(strength), defaultGc);
      min_kappa = std::max(min_kappa, kappa_strength);
      MOFEM_LOG("EP", logSev)
          << "Adjusted interface min kappa for initial Mode-I strength = "
          << min_kappa << ", resulting strength = "
          << GriffithCohesiveLaw::getInitialStrength(min_kappa, defaultGc);
    }
 
    MoFEMFunctionReturn(0);
  }
 
  double defaultGc = 1.0;
  boost::shared_ptr<double> gcPtr;
  Sev logSev;
};
 
 
static Tag get_tag(moab::Interface &moab, std::string tag_name, int size) {
  std::vector<double> dummy(size, 0.);
  Tag tag;
  rval = moab.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE, tag,
                             MB_TAG_CREAT | MB_TAG_SPARSE, dummy.data());
  if (rval == MB_ALREADY_ALLOCATED)
    rval = MB_SUCCESS;
  CHK_MOAB_THROW(rval, "Failed to get tag " + tag_name);
  return tag;
};
 
Tag get_kappa_tag(moab::Interface &moab) {
  return get_tag(moab, "KAPPA", 1);
}
 
Tag get_delta_kappa_tag(moab::Interface &moab) {
  return get_tag(moab, "DELTA_KAPPA", 1);
}
 
struct OpBrokenBaseCohesive
    : public FormsIntegrators<FaceElementForcesAndSourcesCore::
                                  UserDataOperator>::Assembly<A>::OpBrokenBase {
 
  using OP =
      FormsIntegrators<FaceElementForcesAndSourcesCore::UserDataOperator>::
          Assembly<A>::OpBrokenBase;
  OpBrokenBaseCohesive(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_flux_data_ptr,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_flux_data_ptr, ents_ptr) {}
 
  MoFEMErrorCode 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;
#ifndef NDEBUG
      auto base = row_data.getBase();
      if (base < 0 && base >= LASTBASE) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
                "row base not set properly");
      }
#endif // NDEBUG
 
      // 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();
      OP::locMat.resize(OP::nbRows, OP::nbRows, false);
      OP::locMat.clear();
      if (OP::nbRows) {
        // 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) {
        fluxMatPtr =
            boost::shared_ptr<MatrixDouble>(brokenBaseSideData, &bd.getFlux());
        faceSense = bd.getSense();
        CHKERR do_work_rhs(bd.getSide(), bd.getType(), bd.getData());
        fluxMatPtr.reset();
        faceSense = 0;
      }
      break;
    default:
      CHK_MOAB_THROW(MOFEM_IMPOSSIBLE_CASE,
                     (std::string("wrong op type ") +
                      OpBaseDerivativesBase::OpTypeNames[OP::opType])
                         .c_str());
    }
 
    MoFEMFunctionReturn(0);
  }
 
protected:
  boost::shared_ptr<MatrixDouble> fluxMatPtr;
  int faceSense = 0;
};
 
struct OpCohesiveRhs : public OpBrokenBaseCohesive {
 
  using OP = OpBrokenBaseCohesive;
  OpCohesiveRhs(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<double> gc_ptr,
      boost::shared_ptr<VectorDouble> kappa_ptr,
      boost::shared_ptr<VectorDouble> kappa_delta_ptr,
      boost::shared_ptr<std::array<MatrixDouble, 2>> lambda_ptr = nullptr,
      Tag dissipation_tags = Tag(), Tag grad_dissipation_tags = Tag(),
      SmartPetscObj<Vec> vec_dJ_dx = SmartPetscObj<Vec>(),
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpBrokenBaseCohesive(broken_base_side_data, ents_ptr), gcPtr(gc_ptr),
        kappaPtr(kappa_ptr), kappaDeltaPtr(kappa_delta_ptr),
        lambdaPtr(lambda_ptr), dissipationTag(dissipation_tags),
        gradDissipationTag(grad_dissipation_tags), vec_dJdu(vec_dJ_dx) {}
 
  MoFEMErrorCode iNtegrate(EntData &data) {
    MoFEMFunctionBegin;
 
    auto get_sense_index = [this]() { return (faceSense == 1) ? 0 : 1; };
 
    auto nb_dofs = data.getIndices().size();
    FTENSOR_INDEX(3, i);
    FTENSOR_INDEX(3, J);
 
    int nb_integration_pts = OP::getGaussPts().size2();
 
    auto t_P = getFTensor2FromMat<3, 3>(*(OP::fluxMatPtr));
    auto t_kappa = getFTensor0FromVec<0>(*kappaPtr);
    auto t_delta_kappa = getFTensor0FromVec<0>(*kappaDeltaPtr);
    auto t_face_normal = getFTensor1NormalsAtGaussPts();
 
    int nb_base_functions = data.getN().size2() / 3;
    auto t_row_base_fun = data.getFTensor1N<3>();
    auto t_w = OP::getFTensor0IntegrationWeight();
 
    auto next = [&]() {
      ++t_P;
      ++t_face_normal;
      ++t_kappa;
      ++t_delta_kappa;
      ++t_w;
    };
 
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      FTensor::Tensor1<double, 3> t_normal;
      t_normal(J) = t_face_normal(J) * (faceSense / 2.);
      FTensor::Tensor1<double, 3> t_normalized_normal;
      t_normalized_normal(J) = t_normal(J);
      t_normalized_normal.normalize();
      FTensor::Tensor1<double, 3> t_traction;
      t_traction(i) = t_P(i, J) * t_normalized_normal(J);
 
      auto fracture = [this](auto &t_traction, auto &t_normalized_normal,
                             auto &t_kappa, auto &t_delta_kappa, auto gc) {
        double kappa = static_cast<double>(t_kappa + t_delta_kappa);
        auto alpha = GriffithCohesiveLaw::getAlpha(
            kappa, gc, OpGetParameters::min_stiffness);
        auto [teff, t_gap] = GriffithCohesiveLaw::calculateGap(
            t_traction, t_normalized_normal, alpha, OpGetParameters::beta,
            true);
        FTENSOR_INDEX(3, i);
        FTensor::Tensor1<double, 3> t_gap_double;
        t_gap_double(i) = -t_gap(i) / 2.;
        return t_gap_double;
      };
 
      auto assemble = [&](const FTensor::Tensor1<double, 3> &&t_gap) {
        auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
        int bb = 0;
        for (; bb != nb_dofs / SPACE_DIM; ++bb) {
          double row_base = t_w * (t_row_base_fun(J) * t_normal(J));
          t_nf(i) += row_base * t_gap(i);
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
      };
 
      assemble(fracture(t_traction, t_normalized_normal, t_kappa, t_delta_kappa,
                        *gcPtr));
 
      next();
    }
 
    MoFEMFunctionReturn(0);
  };
 
protected:
  boost::shared_ptr<MatrixDouble> uGammaPtr;
  boost::shared_ptr<double> gcPtr;
  boost::shared_ptr<VectorDouble> kappaPtr;
  boost::shared_ptr<VectorDouble> kappaDeltaPtr;
  boost::shared_ptr<std::array<MatrixDouble, 2>> lambdaPtr;
  boost::shared_ptr<double> totalDissipation;
  boost::shared_ptr<double> tatalDissipationGrad;
  SmartPetscObj<Vec> vec_dJdu;
  Tag dissipationTag;
  Tag gradDissipationTag;
};
 
struct OpCohesiveLhs_dPdP : public OpCohesiveRhs {
 
  using OP = OpCohesiveRhs;
  using OpCohesiveRhs::OpCohesiveRhs;
 
  MoFEMErrorCode iNtegrate(EntData &data) {
    MoFEMFunctionBegin;
 
    auto get_sense_index = [this]() { return (faceSense == 1) ? 0 : 1; };
 
    auto nb_dofs = data.getIndices().size();
    FTENSOR_INDEX(3, i);
    FTENSOR_INDEX(3, j);
    FTENSOR_INDEX(3, J);
 
    int nb_integration_pts = OP::getGaussPts().size2();
 
    auto t_P = getFTensor2FromMat<3, 3>(*(OP::fluxMatPtr));
    auto t_kappa = getFTensor0FromVec<0>(*kappaPtr);
    auto t_delta_kappa = getFTensor0FromVec<0>(*kappaDeltaPtr);
    auto t_face_normal = getFTensor1NormalsAtGaussPts();
 
    int nb_base_functions = data.getN().size2() / 3;
    auto t_row_base_fun = data.getFTensor1N<3>();
    auto t_w = OP::getFTensor0IntegrationWeight();
 
    auto next = [&]() {
      ++t_P;
      ++t_kappa;
      ++t_delta_kappa;
      ++t_face_normal;
      ++t_w;
    };
 
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      FTensor::Tensor1<double, 3> t_normal;
      t_normal(J) = t_face_normal(J) * (faceSense / 2.);
      FTensor::Tensor1<double, 3> t_normalized_normal;
      t_normalized_normal(J) = t_normal(J);
      t_normalized_normal.normalize();
      FTensor::Tensor1<double, 3> t_traction;
      t_traction(i) = t_P(i, J) * t_normalized_normal(J);
 
      auto fracture = [this](auto &t_traction, auto &t_normalized_normal,
                             auto &t_kappa, auto &t_delta_kappa, auto gc) {
        double kappa = static_cast<double>(t_kappa + t_delta_kappa);
        auto alpha = GriffithCohesiveLaw::getAlpha(
            kappa, gc, OpGetParameters::min_stiffness);
        auto t_dgap = GriffithCohesiveLaw::calculateDiffGapDTraction(
            t_traction, t_normalized_normal, alpha, OpGetParameters::beta,
            true);
        FTENSOR_INDEX(3, i);
        FTENSOR_INDEX(3, j);
        FTensor::Tensor2<double, 3, 3> t_dgap_double;
        t_dgap_double(i, j) = -t_dgap(i, j) / 2.;
        return t_dgap_double;
      };
 
      auto assemble = [&](const FTensor::Tensor2<double, 3, 3> &&t_dgap) {
        int rr = 0;
        for (; rr != nb_dofs / SPACE_DIM; ++rr) {
          auto t_mat = getFTensor2FromArray<3, 3, 3>(OP::locMat, 3 * rr);
          double row_base = t_w * (t_row_base_fun(J) * t_normal(J));
          auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
          for (int cc = 0; cc != nb_dofs / SPACE_DIM; ++cc) {
            double col_base = t_col_base_fun(J) * t_normalized_normal(J);
            t_mat(i, j) += (row_base * col_base) * t_dgap(i, j);
            ++t_mat;
            ++t_col_base_fun;
          }
          ++t_row_base_fun;
        }
        for (; rr != nb_base_functions; ++rr)
          ++t_row_base_fun;
      };
 
      assemble(fracture(t_traction, t_normalized_normal, t_kappa, t_delta_kappa,
                        *gcPtr));
 
      next();
    }
 
    MoFEMFunctionReturn(0);
  };
 
  MoFEMErrorCode aSsemble(EntData &data) {
    MoFEMFunctionBegin;
 
    if (!this->timeScalingFun.empty())
      this->locMat *= this->timeScalingFun(this->getFEMethod()->ts_t);
    if (!this->feScalingFun.empty())
      this->locMat *= this->feScalingFun(this->getFEMethod());
    // assemble local matrix
    CHKERR this->matSetValuesHook(this, data, data, this->locMat);
 
    MoFEMFunctionReturn(0);
  }
 
private:
};
 
struct OpCohesive_dJ_dkappa : public OpCohesiveRhs {
 
  using OP = OpCohesiveRhs;
  using OpCohesiveRhs::OpCohesiveRhs;
 
  MoFEMErrorCode iNtegrate(EntData &data) {
    MoFEMFunctionBegin;
 
    auto get_sense_index = [this]() { return (faceSense == 1) ? 0 : 1; };
 
    FTENSOR_INDEX(3, i);
    FTENSOR_INDEX(3, J);
    FTENSOR_INDEX(3, K);
 
    int nb_integration_pts = OP::getGaussPts().size2();
 
    auto t_P = getFTensor2FromMat<3, 3>(*(OP::fluxMatPtr));
    auto t_lambda = getFTensor2FromMat<3, 3>(lambdaPtr->at(get_sense_index()));
    auto t_kappa = getFTensor0FromVec<0>(*kappaPtr);
    auto t_delta_kappa = getFTensor0FromVec<0>(*kappaDeltaPtr);
    auto t_face_normal = getFTensor1NormalsAtGaussPts();
    auto t_w = OP::getFTensor0IntegrationWeight();
 
    auto next = [&]() {
      ++t_P;
      ++t_face_normal;
      ++t_kappa;
      ++t_delta_kappa;
      ++t_lambda;
      ++t_w;
    };
 
    double face_dissipation = 0.0;
    double face_grad_dissipation = 0.0;
    auto face_handle = OP::getFEEntityHandle();
    CHKERR getMoab().tag_get_data(dissipationTag, &face_handle, 1,
                                  &face_dissipation);
    CHKERR getMoab().tag_get_data(gradDissipationTag, &face_handle, 1,
                                  &face_grad_dissipation);
 
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      FTensor::Tensor1<double, 3> t_normal;
      t_normal(J) = t_face_normal(J) * (faceSense / 2.);
      FTensor::Tensor1<double, 3> t_normalized_normal;
      t_normalized_normal(J) = t_normal(J);
      t_normalized_normal.normalize();
      FTensor::Tensor1<double, 3> t_traction;
      t_traction(i) = t_P(i, J) * t_normalized_normal(J);
 
      auto dJ_dkappa = [](auto &t_delta_kappa, auto &t_traction,
                          auto &t_normalized_normal, auto &t_kappa, auto gc) {
        double kappa = static_cast<double>(t_kappa);
        double delta_kappa = static_cast<double>(t_delta_kappa);
        double teff = GriffithCohesiveLaw::calculateEffectiveTraction(
            t_traction, t_normalized_normal, OpGetParameters::beta);
        double m = GriffithCohesiveLaw::calculateDissipationSurplus(
            delta_kappa, teff, kappa, gc, OpGetParameters::min_stiffness);
        double m_grad =
            GriffithCohesiveLaw::calculateDissipationSurplusDiffKappa(
                delta_kappa, teff, kappa, gc, OpGetParameters::min_stiffness);
        return boost::make_tuple(m, m_grad);
      };
 
      auto dr_kappa = [](auto &t_delta_kappa, auto &t_traction,
                         auto &t_normalized_normal, auto &t_kappa, auto gc) {
        double kappa = static_cast<double>(t_kappa);
        double delta_kappa = static_cast<double>(t_delta_kappa);
        double kappa_plus_delta = kappa + delta_kappa;
        double diff_alpha = GriffithCohesiveLaw::getDiffAlpha(
            kappa_plus_delta, gc, OpGetParameters::min_stiffness);
        auto [teff, t_gap] = GriffithCohesiveLaw::calculateGap(
            t_traction, t_normalized_normal, diff_alpha, OpGetParameters::beta);
        FTENSOR_INDEX(3, i);
        FTensor::Tensor1<double, 3> t_gap_double;
        t_gap_double(i) = -t_gap(i);
        return t_gap_double;
      };
 
      auto [J, dJ] = dJ_dkappa(t_delta_kappa, t_traction, t_normalized_normal,
                               t_kappa, *gcPtr);
      face_dissipation += t_w * J * t_normal.l2();
      face_grad_dissipation += t_w * dJ * t_normal.l2();
 
      auto t_dr = dr_kappa(t_delta_kappa, t_traction, t_normalized_normal,
                           t_kappa, *gcPtr);
      FTensor::Tensor1<double, 3> t_l;
      t_l(i) = t_lambda(i, K) * t_normal(K);
      face_grad_dissipation -= t_w * t_l(i) * t_dr(i);
 
      next();
    }
 
    CHKERR getMoab().tag_set_data(dissipationTag, &face_handle, 1,
                                  &face_dissipation);
    CHKERR getMoab().tag_set_data(gradDissipationTag, &face_handle, 1,
                                  &face_grad_dissipation);
 
 
    MoFEMFunctionReturn(0);
  };
 
  MoFEMErrorCode aSsemble(EntData &data) {
    MoFEMFunctionBegin;
    MoFEMFunctionReturn(0);
  }
 
protected:
};
 
struct OpCohesive_dJ_dP : public OpCohesiveRhs {
 
  using OP = OpCohesiveRhs;
  using OpCohesiveRhs::OpCohesiveRhs;
 
  MoFEMErrorCode iNtegrate(EntData &data) {
    MoFEMFunctionBegin;
 
    FTENSOR_INDEX(3, i);
    FTENSOR_INDEX(3, J);
 
    int nb_integration_pts = OP::getGaussPts().size2();
    auto nb_dofs = data.getIndices().size();
 
    int nb_base_functions = data.getN().size2() / 3;
    auto t_row_base_fun = data.getFTensor1N<3>();
    auto t_w = OP::getFTensor0IntegrationWeight();
 
    auto t_P = getFTensor2FromMat<3, 3>(*(OP::fluxMatPtr));
    auto t_kappa = getFTensor0FromVec<0>(*kappaPtr);
    auto t_delta_kappa = getFTensor0FromVec<0>(*kappaDeltaPtr);
    auto t_face_normal = getFTensor1NormalsAtGaussPts();
 
    auto next = [&]() {
      ++t_P;
      ++t_face_normal;
      ++t_kappa;
      ++t_delta_kappa;
      ++t_w;
    };
 
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
      FTensor::Tensor1<double, 3> t_normal;
      t_normal(J) = t_face_normal(J) * (faceSense / 2.);
      FTensor::Tensor1<double, 3> t_normalized_normal;
      t_normalized_normal(J) = t_normal(J);
      t_normalized_normal.normalize();
      FTensor::Tensor1<double, 3> t_traction;
      t_traction(i) = t_P(i, J) * t_normalized_normal(J);
 
      auto dJ_dtraction = [](auto &t_traction, auto &t_normalized_normal,
                             auto &t_delta_kappa, auto &t_kappa, auto gc) {
        double kappa = static_cast<double>(t_kappa);
        double delta_kappa = static_cast<double>(t_delta_kappa);
        FTENSOR_INDEX(3, i);
        FTensor::Tensor1<double, 3> t_traction_double;
        t_traction_double(i) = t_traction(i);
        auto t_dM =
            GriffithCohesiveLaw::calculateDissipationSurplusDiffTraction(
                t_traction_double, delta_kappa, kappa, t_normalized_normal, gc,
                OpGetParameters::beta, OpGetParameters::min_stiffness);
        FTensor::Tensor1<double, 3> t_dJ_double;
        t_dJ_double(i) = t_dM(i);
        return t_dJ_double;
      };
 
      auto assemble = [&](const FTensor::Tensor1<double, 3> &&t_dJ) {
 
        auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
        int bb = 0;
        for (; bb != nb_dofs / SPACE_DIM; ++bb) {
          double row_base = t_w * (t_row_base_fun(J) * t_normal(J));
          t_nf(i) += row_base * t_dJ(i);
          ++t_nf;
          ++t_row_base_fun;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_row_base_fun;
      };
 
      assemble(dJ_dtraction(t_traction, t_normalized_normal, t_delta_kappa,
                            t_kappa, *gcPtr));
 
      // cerr << "AAAA "  << OP::locF << endl;
      next();
    }
 
    MoFEMFunctionReturn(0);
  };
 
  MoFEMErrorCode aSsemble(EntData &data) {
    MoFEMFunctionBegin;
    return VecSetValues<AssemblyTypeSelector<A>>(vec_dJdu, data, OP::locF,
                                                 ADD_VALUES);
    MoFEMFunctionReturn(0);
  }
 
protected:
};
 
enum TagGetType { TagSET, TagGET };
 
struct OpGetTag : public OpBrokenBaseCohesive {
 
  using OP = OpBrokenBaseCohesive;
 
  OpGetTag(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      Tag tag, TagGetType tag_get_type,
      boost::shared_ptr<VectorDouble> tag_data_ptr,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpBrokenBaseCohesive(broken_base_side_data, ents_ptr), tagHandle(tag),
        tagGetType(tag_get_type), tagDataPtr(tag_data_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntData &data) {
    MoFEMFunctionBegin;
 
    auto get_sense_index = [this]() { return (faceSense == 1) ? 0 : 1; };
    auto &v = *tagDataPtr;
 
    switch (tagGetType) {
    case TagSET:
      break;
    case TagGET:
      v.resize(1);
      v.clear();
      break;
    }
 
    auto get_data = [&](auto &v) {
      MoFEMFunctionBegin;
 
      auto &moab = getMoab();
      auto fe_ent = getFEEntityHandle();
 
      int size;
      double *data;
      rval = moab.tag_get_by_ptr(tagHandle, &fe_ent, 1, (const void **)&data,
                                 &size);
      if (
 
          rval == MB_SUCCESS && size > 0 && size != v.size()
 
      ) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Inconsistent size of tag data");
      } else {
        if (rval != MB_SUCCESS || tagGetType == TagSET) {
          int tag_size[1];
          tag_size[0] = v.size();
          void const *tag_data[] = {&v[0]};
          CHKERR moab.tag_set_by_ptr(tagHandle, &fe_ent, 1, tag_data, tag_size);
        } else {
          CHKERR moab.tag_get_by_ptr(tagHandle, &fe_ent, 1,
                                     (const void **)&data, &size);
          std::copy(data, data + size, v.begin());
        }
      }
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR get_data(v);
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode aSsemble(EntData &data) {
    MoFEMFunctionBegin;
    MoFEMFunctionReturn(0);
  }
 
private:
  Tag tagHandle;
  boost::shared_ptr<VectorDouble> tagDataPtr;
  TagGetType tagGetType;
};
 
static auto pushCohesiveOpsDomainImpl(
    EshelbianCore &ep,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    SmartPetscObj<Vec> lambda_vec = SmartPetscObj<Vec>()) {
 
  using EleOnSide =
      PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
  using SideEleOp = EleOnSide::UserDataOperator;
 
  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pip, {HDIV, H1, L2}, ep.materialH1Positions, ep.frontAdjEdges);
 
  auto domain_side_flux = [&](auto &pip) {
    // flux
    auto broken_data_ptr =
        boost::make_shared<std::vector<BrokenBaseSideData>>();
    pip.push_back(new OpGetBrokenBaseSideData<SideEleOp>(ep.piolaStress,
                                                         broken_data_ptr));
    auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
        ep.piolaStress, flux_mat_ptr));
    pip.push_back(new OpSetFlux<SideEleOp>(broken_data_ptr, flux_mat_ptr));
 
    return broken_data_ptr;
  };
 
  auto get_lambda = [&](auto &pip) {
    boost::shared_ptr<std::array<MatrixDouble, 2>> array_lambda_ptr;
    if (lambda_vec) {
      array_lambda_ptr = boost::make_shared<std::array<MatrixDouble, 2>>();
      auto lambda_mat_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
          ep.piolaStress, lambda_mat_ptr, boost::make_shared<double>(1.0),
          lambda_vec));
      using OP = VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator;
      auto op = new OP(NOSPACE, OP::OPSPACE);
      op->doWorkRhsHook =
          [lambda_mat_ptr, array_lambda_ptr](
              DataOperator *base_op_ptr, int side, EntityType type,
              EntitiesFieldData::EntData &data) -> MoFEMErrorCode {
        MoFEMFunctionBegin;
        auto op_ptr = static_cast<OP *>(base_op_ptr);
        auto get_sense_index = [op_ptr]() {
          return (op_ptr->getSkeletonSense() == 1) ? 0 : 1;
        };
        array_lambda_ptr->at(get_sense_index()) = *lambda_mat_ptr;
        MoFEMFunctionReturn(0);
      };
      pip.push_back(op);
    }
    return array_lambda_ptr;
  };
 
  return std::make_pair(domain_side_flux(pip), get_lambda(pip));
};
 
static auto pushCohesiveOpsImpl(
    EshelbianCore &ep,
    ForcesAndSourcesCore::GaussHookFun set_integration_at_front_face,
    boost::shared_ptr<Range> interface_range_ptr,
    SmartPetscObj<Vec> lambda_vec = SmartPetscObj<Vec>()) {
 
  auto &m_field = ep.mField;
 
  using BoundaryEle =
      PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
  using EleOnSide =
      PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
  using SideEleOp = EleOnSide::UserDataOperator;
  using BdyEleOp = BoundaryEle::UserDataOperator;
 
  auto face_side = [&]() {
    // First: Iterate over skeleton FEs adjacent to Domain FEs
    // Note:  BoundaryEle, i.e. uses skeleton interation rule
    auto op_loop_skeleton_side =
        new OpLoopSide<BoundaryEle>(m_field, ep.skeletonElement, SPACE_DIM - 1,
                                    interface_range_ptr, Sev::noisy);
    op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
      return -1;
    };
    op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
        set_integration_at_front_face;
 
    EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        op_loop_skeleton_side->getOpPtrVector(), {L2}, ep.materialH1Positions,
        ep.frontAdjEdges);
 
    return op_loop_skeleton_side;
  };
 
  auto op_loop_skeleton_side = face_side();
 
  // Second: Iterate over domain FEs adjacent to skelton, particularly one
  // domain element.
  // Note: EleOnSide, i.e. uses on domain projected skeleton rule
  auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
      m_field, ep.elementVolumeName, SPACE_DIM, Sev::noisy);
  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
  auto [broken_data_flux_ptr, lambda_ptr] = pushCohesiveOpsDomainImpl(
      ep, op_loop_domain_side->getOpPtrVector(), lambda_vec);
  auto op_reset_broken_data = new ForcesAndSourcesCore::UserDataOperator(
      NOSPACE, ForcesAndSourcesCore::UserDataOperator::OPSPACE);
  op_reset_broken_data->doWorkRhsHook =
      [broken_data_flux_ptr](DataOperator *base_op_ptr, int side,
                             EntityType type,
                             EntitiesFieldData::EntData &data) {
        broken_data_flux_ptr->resize(0);
        return 0;
      };
  op_loop_skeleton_side->getOpPtrVector().push_back(op_reset_broken_data);
  op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
 
  auto kappa_ptr = boost::make_shared<VectorDouble>();
  auto kappa_delta_ptr = boost::make_shared<VectorDouble>();
  auto gc_ptr = boost::make_shared<double>();
 
  op_loop_skeleton_side->getOpPtrVector().push_back(
      new OpGetParameters(gc_ptr));
  op_loop_skeleton_side->getOpPtrVector().push_back(
      new OpGetTag(broken_data_flux_ptr, get_kappa_tag(m_field.get_moab()),
                   TagGetType::TagGET, kappa_ptr));
  op_loop_skeleton_side->getOpPtrVector().push_back(new OpGetTag(
      broken_data_flux_ptr, get_delta_kappa_tag(m_field.get_moab()),
      TagGetType::TagGET, kappa_delta_ptr));
 
  return std::make_tuple(op_loop_skeleton_side, broken_data_flux_ptr, gc_ptr,
                         kappa_ptr, kappa_delta_ptr);
};
 
MoFEMErrorCode pushCohesiveOpsRhs(
    EshelbianCore &ep,
    ForcesAndSourcesCore::GaussHookFun set_integration_at_front_face,
    boost::shared_ptr<Range> interface_range_ptr,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip) {
  MoFEMFunctionBegin;
 
  auto [op_loop_skeleton_side, broken_data_flux_ptr, gc_ptr, kappa_ptr,
        kappa_delta_ptr] =
      pushCohesiveOpsImpl(ep, set_integration_at_front_face,
                          interface_range_ptr);
 
  auto u_gamma_ptr = boost::make_shared<MatrixDouble>();
  op_loop_skeleton_side->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(ep.hybridSpatialDisp,
                                                  u_gamma_ptr));
  op_loop_skeleton_side->getOpPtrVector().push_back(new OpCohesiveRhs(
      broken_data_flux_ptr, gc_ptr, kappa_ptr, kappa_delta_ptr));
 
  pip.push_back(op_loop_skeleton_side);
 
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode pushCohesiveOpsLhs(
    EshelbianCore &ep,
    ForcesAndSourcesCore::GaussHookFun set_integration_at_front_face,
    boost::shared_ptr<Range> interface_range_ptr,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip) {
  MoFEMFunctionBegin;
 
  auto [op_loop_skeleton_side, broken_data_flux_ptr, gc_ptr, kappa_ptr,
        kappa_delta_ptr] =
      pushCohesiveOpsImpl(ep, set_integration_at_front_face,
                          interface_range_ptr);
  auto u_gamma_ptr = boost::make_shared<MatrixDouble>();
  op_loop_skeleton_side->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(ep.hybridSpatialDisp,
                                                  u_gamma_ptr));
  op_loop_skeleton_side->getOpPtrVector().push_back(new OpCohesiveLhs_dPdP(
      broken_data_flux_ptr, gc_ptr, kappa_ptr, kappa_delta_ptr));
 
  pip.push_back(op_loop_skeleton_side);
 
  MoFEMFunctionReturn(0);
};
 
static auto pushCohesive_dJ_dkappa_Impl(
    EshelbianCore &ep,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    SmartPetscObj<Vec> lambda_vec) {
  auto &m_field = ep.mField;
 
  using EleOnSide =
      PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
  auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
      m_field, ep.elementVolumeName, SPACE_DIM, Sev::noisy);
  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
  auto [broken_data_flux_ptr, lambda_ptr] = pushCohesiveOpsDomainImpl(
      ep, op_loop_domain_side->getOpPtrVector(), lambda_vec);
 
  auto op_reset_broken_data = new ForcesAndSourcesCore::UserDataOperator(
      NOSPACE, ForcesAndSourcesCore::UserDataOperator::OPSPACE);
  op_reset_broken_data->doWorkRhsHook =
      [broken_data_flux_ptr](DataOperator *base_op_ptr, int side,
                             EntityType type,
                             EntitiesFieldData::EntData &data) {
        broken_data_flux_ptr->resize(0);
        return 0;
      };
  pip.push_back(op_reset_broken_data);
  pip.push_back(op_loop_domain_side);
 
  auto tag_dissipation = get_tag(m_field.get_moab(), "COHESIVE_DISSIPATION", 1);
  auto tag_grad_dissipation =
      get_tag(m_field.get_moab(), "COHESIVE_DISSIPATION_GRAD", 1);
 
  auto kappa_ptr = boost::make_shared<VectorDouble>();
  auto kappa_delta_ptr = boost::make_shared<VectorDouble>();
  auto gc_ptr = boost::make_shared<double>();
 
  pip.push_back(new OpGetParameters(gc_ptr, Sev::noisy));
  pip.push_back(new OpGetTag(broken_data_flux_ptr,
                             get_kappa_tag(m_field.get_moab()),
                             TagGetType::TagGET, kappa_ptr));
  pip.push_back(new OpGetTag(broken_data_flux_ptr,
                             get_delta_kappa_tag(m_field.get_moab()),
                             TagGetType::TagGET, kappa_delta_ptr));
 
  pip.push_back(new OpCohesive_dJ_dkappa(
      broken_data_flux_ptr, gc_ptr, kappa_ptr, kappa_delta_ptr, lambda_ptr,
      tag_dissipation, tag_grad_dissipation));
 
  return std::make_pair(tag_dissipation, tag_grad_dissipation);
}
 
static auto pushCohesive_dJ_dx_Impl(
    EshelbianCore &ep,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip) {
  auto &m_field = ep.mField;
 
  using EleOnSide =
      PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
  auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
      m_field, ep.elementVolumeName, SPACE_DIM, Sev::noisy);
  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
  auto [broken_data_flux_ptr, lambda_ptr] =
      pushCohesiveOpsDomainImpl(ep, op_loop_domain_side->getOpPtrVector());
 
  auto op_reset_broken_data = new ForcesAndSourcesCore::UserDataOperator(
      NOSPACE, ForcesAndSourcesCore::UserDataOperator::OPSPACE);
  op_reset_broken_data->doWorkRhsHook =
      [broken_data_flux_ptr](DataOperator *base_op_ptr, int side,
                             EntityType type,
                             EntitiesFieldData::EntData &data) {
        broken_data_flux_ptr->resize(0);
        return 0;
      };
  pip.push_back(op_reset_broken_data);
  pip.push_back(op_loop_domain_side);
 
  auto kappa_ptr = boost::make_shared<VectorDouble>();
  auto kappa_delta_ptr = boost::make_shared<VectorDouble>();
  auto gc_ptr = boost::make_shared<double>();
 
  pip.push_back(new OpGetParameters(gc_ptr, Sev::noisy));
  pip.push_back(new OpGetTag(broken_data_flux_ptr,
                             get_kappa_tag(m_field.get_moab()),
                             TagGetType::TagGET, kappa_ptr));
  pip.push_back(new OpGetTag(broken_data_flux_ptr,
                             get_delta_kappa_tag(m_field.get_moab()),
                             TagGetType::TagGET, kappa_delta_ptr));
  auto dJ_dx_vec = createDMVector(ep.dmElastic);
  pip.push_back(new OpCohesive_dJ_dP(broken_data_flux_ptr, gc_ptr, kappa_ptr,
                                     kappa_delta_ptr, lambda_ptr, Tag(), Tag(),
                                     dJ_dx_vec));
 
  return dJ_dx_vec;
}
 
static MoFEMErrorCode evaluateDissipationAndGradImpl(
    EshelbianCore &ep,
    ForcesAndSourcesCore::GaussHookFun set_integration_at_front_face,
    SmartPetscObj<Vec> lambda_vec,
    CommInterface::EntitiesPetscVector &dissipation_vec,
    CommInterface::EntitiesPetscVector &grad_dissipation_vec) {
  MoFEMFunctionBegin;
 
  using BoundaryEle =
      PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
  using BdyEleOp = BoundaryEle::UserDataOperator;
 
  auto get_face_ele = [&]() {
    auto fe_ptr = boost::make_shared<BoundaryEle>(ep.mField);
    fe_ptr->getRuleHook = [](int, int, int) { return -1; };
    fe_ptr->setRuleHook = set_integration_at_front_face;
    EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        fe_ptr->getOpPtrVector(), {L2}, ep.materialH1Positions,
        ep.frontAdjEdges);
 
    auto interface_face = [&](FEMethod *fe_method_ptr) {
      auto ent = fe_method_ptr->getFEEntityHandle();
      if (
 
          ep.interfaceFaces->find(ent) != ep.interfaceFaces->end()
 
      ) {
        return true;
      };
      return false;
    };
 
    fe_ptr->exeTestHook = interface_face;
 
    return fe_ptr;
  };
 
  auto face_fe = get_face_ele();
  auto [tag_dissipation, tag_grad_dissipation] =
      pushCohesive_dJ_dkappa_Impl(ep, face_fe->getOpPtrVector(), lambda_vec);
 
  constexpr double zero = 0.0;
  CHKERR ep.mField.get_moab().tag_clear_data(tag_dissipation,
                                             *(ep.interfaceFaces), &zero);
  CHKERR ep.mField.get_moab().tag_clear_data(tag_grad_dissipation,
                                             *(ep.interfaceFaces), &zero);
  CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.skeletonElement, face_fe);
  CHKERR CommInterface::updateEntitiesPetscVector(
      ep.mField.get_moab(), dissipation_vec, tag_dissipation);
  CHKERR CommInterface::updateEntitiesPetscVector(
      ep.mField.get_moab(), grad_dissipation_vec, tag_grad_dissipation);
 
  MoFEMFunctionReturn(0);
}
 
static MoFEMErrorCode evaluateCohesiveLambdaImpl(
    EshelbianCore &ep,
    ForcesAndSourcesCore::GaussHookFun set_integration_at_front_face,
    SmartPetscObj<KSP> ksp, SmartPetscObj<Vec> lambda_vec) {
  MoFEMFunctionBegin;
 
  using BoundaryEle =
      PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
  using BdyEleOp = BoundaryEle::UserDataOperator;
 
  auto get_face_ele = [&]() {
    auto fe_ptr = boost::make_shared<BoundaryEle>(ep.mField);
    fe_ptr->getRuleHook = [](int, int, int) { return -1; };
    fe_ptr->setRuleHook = set_integration_at_front_face;
    EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        fe_ptr->getOpPtrVector(), {L2}, ep.materialH1Positions,
        ep.frontAdjEdges);
 
    auto interface_face = [&](FEMethod *fe_method_ptr) {
      auto ent = fe_method_ptr->getFEEntityHandle();
      if (
 
          ep.interfaceFaces->find(ent) != ep.interfaceFaces->end()
 
      ) {
        return true;
      };
      return false;
    };
 
    fe_ptr->exeTestHook = interface_face;
 
    return fe_ptr;
  };
 
  auto face_fe = get_face_ele();
  auto dJ_dx = pushCohesive_dJ_dx_Impl(ep, face_fe->getOpPtrVector());
  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
      ep.dmElastic, ep.skeletonElement, face_fe, 0, ep.mField.get_comm_size());
  CHKERR VecAssemblyBegin(dJ_dx);
  CHKERR VecAssemblyEnd(dJ_dx);
  CHKERR VecGhostUpdateBegin(dJ_dx, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateEnd(dJ_dx, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateBegin(dJ_dx, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(dJ_dx, INSERT_VALUES, SCATTER_FORWARD);
  double dJ_dx_norm2;
  CHKERR VecNorm(dJ_dx, NORM_2, &dJ_dx_norm2);
  MOFEM_LOG("EP", Sev::inform)
      << "evaluateCohesiveLambdaImpl: Norm of dJ/dx vector: " << dJ_dx_norm2;
  constexpr double tol = 1e-16;
  if (dJ_dx_norm2 < tol) {
    CHKERR VecZeroEntries(lambda_vec);
  } else {
    CHKERR KSPSolveTranspose(ksp, dJ_dx, lambda_vec);
  }
  CHKERR VecGhostUpdateBegin(lambda_vec, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(lambda_vec, INSERT_VALUES, SCATTER_FORWARD);
  double lambda_norm2;
  CHKERR VecNorm(lambda_vec, NORM_2, &lambda_norm2);
  MOFEM_LOG("EP", Sev::inform)
      << "evaluateCohesiveLambdaImpl: Norm of lambda vector: " << lambda_norm2;
 
  MoFEMFunctionReturn(0);
}
 
static MoFEMErrorCode evaluatePrimalProblemCohesiveImpl(EshelbianCore &ep,
                                                        SmartPetscObj<TS> ts,
                                                        SmartPetscObj<Vec> x) {
  MoFEMFunctionBegin;
 
  CHKERR TSSetFromOptions(ts);
  CHKERR TSSetStepNumber(ts, 0);
  CHKERR TSSetTime(ts, 0);
  CHKERR TSSetSolution(ts, x);
  double dt;
  CHKERR TSGetTimeStep(ts, &dt);
  MOFEM_LOG("EP", Sev::inform)
      << "evaluatePrimalProblemCohesiveImpl: Time step dt: " << dt;
  CHKERR TSSolve(ts, PETSC_NULLPTR);
  CHKERR TSSetTimeStep(ts, dt);
 
  CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, x, INSERT_VALUES,
                                 SCATTER_FORWARD);
  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
 
  double x_norm2;
  CHKERR VecNorm(x, NORM_2, &x_norm2);
  MOFEM_LOG("EP", Sev::inform)
      << "evaluatePrimalProblemCohesiveImpl: Norm of displacement vector: "
      << x_norm2;
 
  MoFEMFunctionReturn(0);
}
 
struct CohesiveTAOCtxImpl : public CohesiveTAOCtx {
  CohesiveTAOCtxImpl(
      EshelbianCore *ep,
      ForcesAndSourcesCore::GaussHookFun set_integration_at_front_face,
      SmartPetscObj<TS> time_solver)
      : ep_ptr(ep), setIntegrationAtFrontFace(set_integration_at_front_face),
        timeSolver(time_solver),
 
        kspSolVec(createDMVector(ep->dmElastic)),
        lambdaVec(createDMVector(ep->dmElastic)),
        kappaVec(CommInterface::createEntitiesPetscVector(
            ep->mField.get_comm(), ep->mField.get_moab(),
            [&ep](Range r) { return intersect(*(ep->interfaceFaces), r); }, 1,
            Sev::noisy)),
        dissipationVec(CommInterface::createEntitiesPetscVector(
            ep->mField.get_comm(), ep->mField.get_moab(),
            [&ep](Range r) { return intersect(*(ep->interfaceFaces), r); }, 1,
            Sev::noisy)),
        gradDissipationVec(CommInterface::createEntitiesPetscVector(
            ep->mField.get_comm(), ep->mField.get_moab(),
            [&ep](Range r) { return intersect(*(ep->interfaceFaces), r); }, 1,
            Sev::noisy)) {}
 
  SmartPetscObj<Vec> duplicateKappaVec() override {
    return vectorDuplicate(kappaVec.second);
  }
 
  SmartPetscObj<Vec> duplicateDissipationVec() override {
    return vectorDuplicate(dissipationVec.second);
  }
 
  SmartPetscObj<Vec> duplicateGradientVec() override {
    return vectorDuplicate(gradDissipationVec.second);
  }
 
  CommInterface::EntitiesPetscVector &getKappaVec() override {
    return kappaVec;
  }
 
private:
  EshelbianCore *ep_ptr;
  ForcesAndSourcesCore::GaussHookFun setIntegrationAtFrontFace;
  SmartPetscObj<TS> timeSolver;
 
  SmartPetscObj<Vec> kspSolVec;
  SmartPetscObj<Vec> lambdaVec;
  CommInterface::EntitiesPetscVector kappaVec;
  CommInterface::EntitiesPetscVector dissipationVec;
  CommInterface::EntitiesPetscVector gradDissipationVec;
  friend MoFEMErrorCode cohesiveEvaluateObjectiveAndGradient(Tao tao, Vec x,
                                                             PetscReal *f,
                                                             Vec g, void *ctx);
};
 
boost::shared_ptr<CohesiveTAOCtx> createCohesiveTAOCtx(
    EshelbianCore *ep_ptr,
    ForcesAndSourcesCore::GaussHookFun set_integration_at_front_face,
    SmartPetscObj<TS> time_solver) {
  return boost::make_shared<CohesiveTAOCtxImpl>(
      ep_ptr, set_integration_at_front_face, time_solver);
}
 
MoFEMErrorCode cohesiveEvaluateObjectiveAndGradient(Tao tao, Vec delta_kappa,
                                                    PetscReal *f, Vec g,
                                                    void *ctx) {
 
  MoFEMFunctionBegin;
  auto cohesive_ctx = static_cast<CohesiveTAOCtxImpl *>(ctx);
  auto &ep = *(cohesive_ctx->ep_ptr);
 
#ifndef NDEBUG
  CHKERR VecView(delta_kappa, PETSC_VIEWER_STDOUT_WORLD);
#endif
  // Set delta_kappa values to the tag
  CHKERR VecCopy(delta_kappa, cohesive_ctx->kappaVec.second);
  CHKERR VecGhostUpdateBegin(cohesive_ctx->kappaVec.second, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(cohesive_ctx->kappaVec.second, INSERT_VALUES,
                           SCATTER_FORWARD);
  CHKERR CommInterface::setTagFromVector(
      ep.mField.get_moab(), cohesive_ctx->kappaVec,
      get_delta_kappa_tag(ep.mField.get_moab()));
 
  // solve primal problem
  CHKERR evaluatePrimalProblemCohesiveImpl(ep, cohesive_ctx->timeSolver,
                                           cohesive_ctx->kspSolVec);
  // solve adjoint problem to get lambda
  auto ksp = snesGetKSP(tsGetSNES(cohesive_ctx->timeSolver));
  CHKERR evaluateCohesiveLambdaImpl(ep, cohesive_ctx->setIntegrationAtFrontFace,
                                    ksp, cohesive_ctx->lambdaVec);
  // evaluate dissipation and its gradient
  CHKERR evaluateDissipationAndGradImpl(
      ep, cohesive_ctx->setIntegrationAtFrontFace, cohesive_ctx->lambdaVec,
      cohesive_ctx->dissipationVec, cohesive_ctx->gradDissipationVec);
 
  CHKERR VecSum(cohesive_ctx->dissipationVec.second, f);
  CHKERR VecCopy(cohesive_ctx->gradDissipationVec.second, g);
 
  MOFEM_LOG("EP", Sev::inform)
      << "Cohesive objective function (negative total dissipation): " << *f;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode initializeCohesiveKappaField(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  auto th = get_kappa_tag(ep.mField.get_moab());
  CHKERR ep.mField.get_moab().tag_clear_data(th, *(ep.interfaceFaces),
                                             &OpGetParameters::min_kappa);
  ep.listTagsToTransfer.push_back(th);
  MoFEMFunctionReturn(0);
}
 
}; // namespace EshelbianPlasticity