EshelbianTopologicalDerivativeOperators.cpp

Go to the documentation of this file.

/**
 * @file EshelbianTopologicalDerivative.cpp
 * @brief
 * @version 0.1
 * @date 2026-02-11
 *
 * @copyright Copyright (c) 2026
 *
 */
 
using namespace EshelbianPlasticity;
using namespace ShapeOptimization;
 
namespace EshelbianPlasticity {
 
// dr/dX topological derivative for Eshelbian plasticity model
 
template <typename AssembleOp>
struct OpAssembleTopologicalObjectiveDerivativeImplBase : public AssembleOp {
  using OP = AssembleOp;
 
  OpAssembleTopologicalObjectiveDerivativeImplBase(
      const std::string &field_name,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<TopologicalData> topo_ptr,
      boost::shared_ptr<double> J_ptr, SmartPetscObj<Vec> assemble_vec,
      Tag topo_tag)
      : OP(field_name, data_ptr, OP::OPROW), JPtr(J_ptr),
        assembleVec(assemble_vec), topoTag(topo_tag), topoData(topo_ptr) {}
 
  MoFEMErrorCode assemble(int side, EntityType type, EntData &data) override {
    MoFEMFunctionBegin;
    if (type == MBVERTEX) {
      // Note it is iterated over vertices, since geometry is now in H1.
      if (JPtr) {
        *JPtr += locJ;
      }
    }
    if (assembleVec) {
      double *vec_ptr = OP::nF.data().data();
      const int nb_dofs = data.getIndices().size();
      int *ind_ptr = data.getIndices().data().data();
      CHKERR VecSetValues(assembleVec, nb_dofs, ind_ptr, vec_ptr, ADD_VALUES);
    }
    if (topoTag) {
      const auto field_ents = data.getFieldEntities();
      std::vector<EntityHandle> ents(field_ents.size());
      std::transform(field_ents.begin(), field_ents.end(), ents.begin(),
                     [](const auto *fe) { return fe->getEnt(); });
      if (field_ents.empty())
        MoFEMFunctionReturnHot(0);
      if (type_from_handle(ents[0]) != MBVERTEX)
        MoFEMFunctionReturnHot(0);
      auto &moab = OP::getMoab();
      VectorDouble topo_values(OP::nF.size());
      CHKERR moab.tag_set_data(topoTag, ents.data(), ents.size(),
                               topo_values.data().data());
      noalias(topo_values) += OP::nF;
      CHKERR moab.tag_set_data(topoTag, ents.data(), ents.size(),
                               OP::nF.data().data());
    }
    MoFEMFunctionReturn(0);
  }
 
protected:
  double locJ;
  boost::shared_ptr<double> JPtr;
  SmartPetscObj<Vec> assembleVec;
  Tag topoTag;
  boost::shared_ptr<TopologicalData> topoData;
};
 
struct OpAssembleVolumeTopologicalDerivativeImpl
    : public OpAssembleTopologicalObjectiveDerivativeImplBase<
          OpAssembleVolume> {
  using OpAssembleTopologicalObjectiveDerivativeImplBase<
      OpAssembleVolume>::OpAssembleTopologicalObjectiveDerivativeImplBase;
};
 
struct OpAssembleFaceTopologicalDerivativeImpl
    : public OpAssembleTopologicalObjectiveDerivativeImplBase<OpAssembleFace> {
  using OpAssembleTopologicalObjectiveDerivativeImplBase<
      OpAssembleFace>::OpAssembleTopologicalObjectiveDerivativeImplBase;
};
 
struct OpAssembleBrokenFaceTopologicalDerivativeImplBase
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBrokenBase {
 
  using OP = typename FormsIntegrators<FaceUserDataOperator>::Assembly<
      A>::OpBrokenBase;
 
  OpAssembleBrokenFaceTopologicalDerivativeImplBase(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<TopologicalData> topo_ptr,
      boost::shared_ptr<double> J_ptr, SmartPetscObj<Vec> assemble_vec,
      Tag topo_tag, boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_base_side_data, ents_ptr), JPtr(J_ptr),
        assembleVec(assemble_vec), topoTag(topo_tag), topoData(topo_ptr) {}
 
  MoFEMErrorCode aSsemble(EntData &data) override {
    MoFEMFunctionBegin;
    if (assembleVec) {
      double *vec_ptr = OP::locF.data().data();
      const int nb_dofs = data.getIndices().size();
      int *ind_ptr = data.getIndices().data().data();
      CHKERR VecSetValues(assembleVec, nb_dofs, ind_ptr, vec_ptr, ADD_VALUES);
    }
    if (topoTag) {
      const auto field_ents = data.getFieldEntities();
      std::vector<EntityHandle> ents(field_ents.size());
      std::transform(field_ents.begin(), field_ents.end(), ents.begin(),
                     [](const auto *fe) { return fe->getEnt(); });
      if (field_ents.empty())
        MoFEMFunctionReturnHot(0);
      if (type_from_handle(ents[0]) != MBVERTEX)
        MoFEMFunctionReturnHot(0);
      auto &moab = getMoab();
      VectorDouble topo_values(OP::locF.size());
      CHKERR moab.tag_set_data(topoTag, ents.data(), ents.size(),
                               topo_values.data().data());
      topo_values += OP::locF;
      CHKERR moab.tag_set_data(topoTag, ents.data(), ents.size(),
                               OP::locF.data().data());
    }
    MoFEMFunctionReturn(0);
  }
 
protected:
  boost::shared_ptr<double> JPtr;
  SmartPetscObj<Vec> assembleVec;
  Tag topoTag;
  boost::shared_ptr<TopologicalData> topoData;
};
 
struct OpJ_dPImpl : public OpAssembleVolumeTopologicalDerivativeImpl {
  using OpAssembleVolumeTopologicalDerivativeImpl::
      OpAssembleVolumeTopologicalDerivativeImpl;
 
  MoFEMErrorCode integrate(EntData &data);
};
 
struct OpJ_dBubbleImpl : public OpAssembleVolumeTopologicalDerivativeImpl {
  using OpAssembleVolumeTopologicalDerivativeImpl::
      OpAssembleVolumeTopologicalDerivativeImpl;
 
  MoFEMErrorCode integrate(EntData &data);
};
 
struct OpJ_dwImpl : public OpAssembleVolumeTopologicalDerivativeImpl {
  using OpAssembleVolumeTopologicalDerivativeImpl::
      OpAssembleVolumeTopologicalDerivativeImpl;
 
  MoFEMErrorCode integrate(EntData &data);
};
 
struct OpJ_dOmegaImpl : public OpAssembleVolumeTopologicalDerivativeImpl {
  using OpAssembleVolumeTopologicalDerivativeImpl::
      OpAssembleVolumeTopologicalDerivativeImpl;
 
  MoFEMErrorCode integrate(EntData &data);
};
 
struct dJ_duGammaImpl : public OpAssembleFaceTopologicalDerivativeImpl {
  using OpAssembleFaceTopologicalDerivativeImpl::
      OpAssembleFaceTopologicalDerivativeImpl;
 
  MoFEMErrorCode integrate(EntData &data);
};
 
struct dJ_dTractionImpl
    : public OpAssembleBrokenFaceTopologicalDerivativeImplBase {
  using OpAssembleBrokenFaceTopologicalDerivativeImplBase::
      OpAssembleBrokenFaceTopologicalDerivativeImplBase;
 
  MoFEMErrorCode iNtegrate(EntData &data);
};
 
struct OpTauStabilisation_dX : public OpAssembleFaceTopologicalDerivativeImpl {
  OpTauStabilisation_dX(
      const std::string &field_name,
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_disp_data_ptr,
      boost::shared_ptr<MatrixDouble> hybrid_disp_ptr,
      boost::shared_ptr<MatrixDouble> var_hybrid_disp_ptr,
      boost::shared_ptr<TopologicalData> topo_ptr, const double alpha_tau,
      SmartPetscObj<Vec> vec, boost::shared_ptr<double> J_ptr = nullptr,
      Tag tag = Tag())
      : OpAssembleFaceTopologicalDerivativeImpl(field_name, nullptr, topo_ptr,
                                                J_ptr, vec, tag),
        brokenDispDataPtr(broken_disp_data_ptr), hybridDispPtr(hybrid_disp_ptr),
        varHybridDispPtr(var_hybrid_disp_ptr), alphaTau(alpha_tau) {}
 
  MoFEMErrorCode integrate(int side, EntityType type, EntData &data) override;
 
private:
  boost::shared_ptr<std::vector<BrokenBaseSideData>> brokenDispDataPtr;
  boost::shared_ptr<MatrixDouble> hybridDispPtr;
  boost::shared_ptr<MatrixDouble> varHybridDispPtr;
  double alphaTau;
};
 
struct OpDispBc_dX : public OpAssembleFaceTopologicalDerivativeImpl {
  OpDispBc_dX(
      const std::string &field_name,
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_side_data_ptr,
      boost::shared_ptr<BcDispVec> bc_disp_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv,
      boost::shared_ptr<TopologicalData> topo_ptr, SmartPetscObj<Vec> vec,
      boost::shared_ptr<double> J_ptr = nullptr, Tag tag = Tag())
      : OpAssembleFaceTopologicalDerivativeImpl(field_name, nullptr, topo_ptr,
                                                J_ptr, vec, tag),
        brokenSideDataPtr(broken_side_data_ptr), bcDispPtr(bc_disp_ptr),
        scalingMethodsMap(smv) {}
 
  MoFEMErrorCode integrate(int side, EntityType type, EntData &data) override;
 
private:
  boost::shared_ptr<std::vector<BrokenBaseSideData>> brokenSideDataPtr;
  boost::shared_ptr<BcDispVec> bcDispPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
template <typename OP_PTR>
static MatrixDouble
getTopologicalAnalyticalExpr(OP_PTR op_ptr, 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 v_analytical_expr;
}
 
struct OpAnalyticalDispBc_dX : public OpAssembleFaceTopologicalDerivativeImpl {
  OpAnalyticalDispBc_dX(
      const std::string &field_name,
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_side_data_ptr,
      boost::shared_ptr<AnalyticalDisplacementBcVec> bc_disp_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv,
      boost::shared_ptr<TopologicalData> topo_ptr, SmartPetscObj<Vec> vec,
      boost::shared_ptr<double> J_ptr = nullptr, Tag tag = Tag())
      : OpAssembleFaceTopologicalDerivativeImpl(field_name, nullptr, topo_ptr,
                                                J_ptr, vec, tag),
        brokenSideDataPtr(broken_side_data_ptr), bcDispPtr(bc_disp_ptr),
        scalingMethodsMap(smv) {}
 
  MoFEMErrorCode integrate(int side, EntityType type, EntData &data) override;
 
private:
  boost::shared_ptr<std::vector<BrokenBaseSideData>> brokenSideDataPtr;
  boost::shared_ptr<AnalyticalDisplacementBcVec> bcDispPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
struct OpBrokenTractionBc_dX : public OpAssembleFaceTopologicalDerivativeImpl {
  OpBrokenTractionBc_dX(
      std::string field_name, boost::shared_ptr<TractionBcVec> bc_data,
      boost::shared_ptr<MatrixDouble> lambda_hybrid_ptr,
      boost::shared_ptr<TopologicalData> topo_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv,
      SmartPetscObj<Vec> vec, boost::shared_ptr<double> J_ptr = nullptr,
      Tag tag = Tag())
      : OpAssembleFaceTopologicalDerivativeImpl(field_name, nullptr, topo_ptr,
                                                J_ptr, vec, tag),
        bcData(bc_data), lambdaHybridPtr(lambda_hybrid_ptr), topoData(topo_ptr),
        scalingMethodsMap(smv) {}
 
  MoFEMErrorCode integrate(int side, EntityType type, EntData &data) override;
 
protected:
  boost::shared_ptr<TractionBcVec> bcData;
  boost::shared_ptr<MatrixDouble> lambdaHybridPtr;
  boost::shared_ptr<TopologicalData> topoData;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
struct OpBrokenAnalyticalTractionBc_dX
    : public OpAssembleFaceTopologicalDerivativeImpl {
  OpBrokenAnalyticalTractionBc_dX(
      std::string field_name,
      boost::shared_ptr<AnalyticalTractionBcVec> bc_data,
      boost::shared_ptr<MatrixDouble> lambda_hybrid_ptr,
      boost::shared_ptr<TopologicalData> topo_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv,
      SmartPetscObj<Vec> vec, boost::shared_ptr<double> J_ptr = nullptr,
      Tag tag = Tag())
      : OpAssembleFaceTopologicalDerivativeImpl(field_name, nullptr, topo_ptr,
                                                J_ptr, vec, tag),
        bcData(bc_data), lambdaHybridPtr(lambda_hybrid_ptr), topoData(topo_ptr),
        scalingMethodsMap(smv) {}
 
  MoFEMErrorCode integrate(int side, EntityType type, EntData &data) override;
 
protected:
  boost::shared_ptr<AnalyticalTractionBcVec> bcData;
  boost::shared_ptr<MatrixDouble> lambdaHybridPtr;
  boost::shared_ptr<TopologicalData> topoData;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
struct OpJ_dUImpl : public OpAssembleVolumeTopologicalDerivativeImpl {
  using OpAssembleVolumeTopologicalDerivativeImpl::
      OpAssembleVolumeTopologicalDerivativeImpl;
 
  MoFEMErrorCode integrate(EntData &data);
};
 
struct OpInteriorJImpl : public OpAssembleVolumeTopologicalDerivativeImpl {
  using OpAssembleVolumeTopologicalDerivativeImpl::
      OpAssembleVolumeTopologicalDerivativeImpl;
 
  MoFEMErrorCode integrate(EntData &data);
};
 
struct OpTopologicalObjectivePythonImpl
    : public ForcesAndSourcesCore::UserDataOperator {
  using OP = ForcesAndSourcesCore::UserDataOperator;
 
  OpTopologicalObjectivePythonImpl(
      boost::shared_ptr<DataAtIntegrationPts> data_at_pts_ptr,
      boost::shared_ptr<TopologicalData> topo_p,
      boost::shared_ptr<ObjectiveFunctionData> python_ptr,
      const ObjectiveModelType eval_energy_model = PYTHON_MODEL)
      : OP(NOSPACE, OP::OPSPACE), dataAtPts(data_at_pts_ptr), topoData(topo_p),
        pythonPtr(python_ptr), evalEnergyModel(eval_energy_model) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  ObjectiveModelType evalEnergyModel = PYTHON_MODEL;
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  boost::shared_ptr<TopologicalData> topoData;
  boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
};
 
MoFEMErrorCode OpTopologicalObjectivePythonImpl::doWork(int side,
                                                        EntityType type,
                                                        EntData &data) {
  MoFEMFunctionBegin;
 
#ifndef NDEBUG
  if (!dataAtPts)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "DataAtIntegrationPts pointer is null");
  if (!topoData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Topological data pointer is null");
#endif // NDEBUG
 
  const int nb_gauss_pts = getGaussPts().size2();
  if (!nb_gauss_pts)
    MoFEMFunctionReturnHot(0);
 
  MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
      *topoData->objConstantRotation, nb_gauss_pts);
  topoData->objConstantRotation->clear();
  const auto &features = dataAtPts->physicsPtr->getFeatures();
  const auto variable_compliance_mask =
      PhysicalEquations::Features{}
          .set(PhysicalEquations::NON_HOMOGENEOUS_MATERIAL)
          .set(PhysicalEquations::NO_STRETCH_NONLINEAR);
 
  auto stress_full_ptr = boost::make_shared<MatrixDouble>();
  auto get_stress_full =
      MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                       DL>::size(*stress_full_ptr, nb_gauss_pts);
  stress_full_ptr->clear();
  auto strain_full_ptr = boost::make_shared<MatrixDouble>();
  auto get_strain_full =
      MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                       DL>::size(*strain_full_ptr, nb_gauss_pts);
  strain_full_ptr->clear();
 
  auto t_stress = get_stress_full();
  auto t_strain = get_strain_full();
 
  auto t_biot = dataAtPts->getFTensorAdjointPdstretch(nb_gauss_pts);
  auto t_u = dataAtPts->getFTensorStretch(nb_gauss_pts);
 
  auto next = [&]() {
    ++t_stress;
    ++t_strain;
    ++t_biot;
    ++t_u;
  };
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l);
 
  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
    // we have to handle all variants, that will render how the Jacobian
    // gradient is evaluated.
    t_stress(i, j) = t_biot(i, j);
    t_strain(i, j) = t_u(i, j);
    next();
  }
 
  auto evaluate_python_objective = [&]() {
    MoFEMFunctionBegin;
    if (!pythonPtr)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "ObjectiveFunctionData pointer is null");
 
    auto &coords = OP::getCoordsAtGaussPts();
    CHKERR pythonPtr->evalInteriorObjectiveFunction(
        coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
        topoData->getObjAtPts(), false);
    CHKERR pythonPtr->evalInteriorObjectiveGradientStrain(
        coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
        topoData->getObjDStrainAtPts(), false);
    CHKERR pythonPtr->evalInteriorObjectiveGradientU(
        coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
        topoData->getObjDDisplacementAtPts(), false);
    CHKERR pythonPtr->evalInteriorObjectiveGradientStress(
        coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
        topoData->getObjDStressAtPts(), false);
    MoFEMFunctionReturn(0);
  };
 
  auto evaluate_energy_of_hencky_model = [&]() {
    MoFEMFunctionBegin;
 
    MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
        *topoData->getObjDDisplacementAtPts(), nb_gauss_pts);
    topoData->getObjDDisplacementAtPts()->clear();
    MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                     DL>::size(*topoData->getObjDStressAtPts(), nb_gauss_pts);
    topoData->getObjDStressAtPts()->clear();
    MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
        *topoData->getObjDRotationAtPts(), nb_gauss_pts);
 
    auto eval_evergy = [&](auto &&t_D) {
      auto get_obj =
          MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::size(
              *topoData->getObjAtPts(), nb_gauss_pts);
      auto get_dstrain_obj =
          MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                           DL>::size(*topoData->getObjDStrainAtPts(),
                                     nb_gauss_pts);
      auto t_obj = get_obj();
      auto t_dstrain_obj = get_dstrain_obj();
      auto t_log_u = dataAtPts->getFTensorLogStretchTotal(nb_gauss_pts);
      for (int gg = 0; gg != nb_gauss_pts; ++gg) {
        t_obj(0) = 0.5 * (t_log_u(i, j) * t_D(i, j, k, l) * t_log_u(k, l));
        t_dstrain_obj(i, j) = t_D(i, j, k, l) * t_log_u(k, l);
        ++t_log_u;
        ++t_obj;
        ++t_dstrain_obj;
        ++t_D;
      }
    };
 
    if (!features.test(PhysicalEquations::NON_HOMOGENEOUS_MATERIAL)) {
      eval_evergy(
          getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(dataAtPts->matD));
    } else {
      eval_evergy(getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->matD));
    }
 
    topoData->getObjDRotationAtPts()->clear();
    MoFEMFunctionReturn(0);
  };
 
  auto evaluate_energy_of_hencky_model_nostreach = [&]() {
    MoFEMFunctionBegin;
 
    MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
        *topoData->getObjDDisplacementAtPts(), nb_gauss_pts);
    topoData->getObjDDisplacementAtPts()->clear();
    MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                     DL>::size(*topoData->getObjDStrainAtPts(), nb_gauss_pts);
    topoData->getObjDStrainAtPts()->clear();
    MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
        *topoData->getObjDRotationAtPts(), nb_gauss_pts);
    topoData->getObjDRotationAtPts()->clear();
 
    auto eval_evergy = [&](auto &&t_inv_D) {
      auto get_obj =
          MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::size(
              *topoData->getObjAtPts(), nb_gauss_pts);
      auto get_dstress_obj =
          MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                           DL>::size(*topoData->getObjDStressAtPts(),
                                     nb_gauss_pts);
      auto t_obj = get_obj();
      auto t_dstress_obj = get_dstress_obj();
      auto t_stress = dataAtPts->getFTensorAdjointPdstretch(nb_gauss_pts);
      for (int gg = 0; gg != nb_gauss_pts; ++gg) {
        t_obj(0) =
            0.5 * (t_stress(i, j) * t_inv_D(i, j, k, l) * t_stress(k, l));
        t_dstress_obj(i, j) = t_inv_D(i, j, k, l) * t_stress(k, l);
        ++t_stress;
        ++t_obj;
        ++t_dstress_obj;
        ++t_inv_D;
      }
    };
 
    if ((features & variable_compliance_mask).none()) {
      eval_evergy(
          getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(dataAtPts->matInvD));
    } else {
      eval_evergy(
          getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->matInvD));
    }
    MoFEMFunctionReturn(0);
  };
 
  auto conversion_of_biot_stress = [&]() {
    MoFEMFunctionBegin;
    // Pull Python's Biot-stress covector back to Piola and both rotation
    // factors.
    MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
        *topoData->getObjDRotationAtPts(), nb_gauss_pts);
    topoData->getObjDRotationAtPts()->clear();
    auto objective = topoData->getFTensorObjDStress(nb_gauss_pts);
    auto R = dataAtPts->getFTensorRotMat(nb_gauss_pts);
    auto P = dataAtPts->getFTensorApproxP(nb_gauss_pts);
    auto gradient = dataAtPts->getFTensorSmallWGradH1(nb_gauss_pts);
    FTENSOR_INDEXES(3, i, k, l, m, o);
    constexpr auto identity = FTensor::Kronecker_Delta<double>();
    for (int gg = 0; gg != nb_gauss_pts; ++gg) {
      FTensor::Tensor2<double, 3, 3> covector, h1;
      covector(l, o) = objective(l, o);
      h1(o, k) = identity(o, k);
      if (EshelbianCore::gradApproximator != SMALL_ROT) {
        if (EshelbianCore::gradApproximator != NO_H1_CONFIGURATION)
          h1(o, k) += gradient(o, k);
        objective(i, k) = R(i, l) * (covector(l, o) * h1(o, k));
        for (int mode = 0; mode != 2; ++mode) {
          FTensor::Tensor3<double, 3, 3, 3> derivative;
          if (EshelbianCore::rotSelector == LARGE_ROT) {
            using namespace LogarithmicDirectionalDerivatives;
            const auto read = [&](MatrixPtr values) {
              Vector result;
              auto view = GetFTensor1FromMatImpl<3, -1, DL, MatrixDouble>::get(
                  *values, gg, 0);
              result(i) = view(i);
              return result;
            };
            const SplitRotation rotation(read(dataAtPts->constantRotation),
                                         read(dataAtPts->higherRotation));
            derivative(i, l, m) = rotation.firstTensor(mode)(i, l, m);
          } else {
            derivative(i, l, m) = levi_civita(i, l, m);
          }
          auto rotation_covector =
              topoData->getFTensorObjDRotation(nb_gauss_pts, mode, gg);
          rotation_covector(m) =
              covector(l, o) * (derivative(i, l, m) * P(i, k)) * h1(o, k);
        }
      }
      ++objective;
      ++R;
      ++P;
      ++gradient;
    }
    MoFEMFunctionReturn(0);
  };
 
  auto conversion_of_stretch = [&]() {
    MoFEMFunctionBegin;
    // Python sees the physical stretch tensor. The stretch field stores log
    // stretch, so convert dJ/dU to dJ/dlogU once here.
    auto t_obj_dstretch = topoData->getFTensorObjDStrain(nb_gauss_pts);
    auto t_diff_stretch = dataAtPts->getFTensorDiffStretch(nb_gauss_pts);
 
    FTENSOR_INDEX(SPACE_DIM, i);
    FTENSOR_INDEX(SPACE_DIM, j);
    FTENSOR_INDEX(SPACE_DIM, k);
    FTENSOR_INDEX(SPACE_DIM, l);
 
    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_dJ_dstretch;
      t_dJ_dstretch(i, j) = t_obj_dstretch(i, j);
 
      t_obj_dstretch(k, l) = t_dJ_dstretch(i, j) * t_diff_stretch(i, j, k, l);
 
      ++t_obj_dstretch;
      ++t_diff_stretch;
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto conversion_of_stretch_to_stress_for_no_stretch = [&](auto t_inv_D) {
    MoFEMFunctionBegin;
    // In no-stretch mode, log stretch is computed from the stress field. Fold
    // the objective stretch derivative into the stress derivative buffer.
    auto t_obj_dstress = topoData->getFTensorObjDStress(nb_gauss_pts);
    auto t_obj_dstretch = topoData->getFTensorObjDStrain(nb_gauss_pts);
 
    FTENSOR_INDEX(SPACE_DIM, i);
    FTENSOR_INDEX(SPACE_DIM, j);
    FTENSOR_INDEX(SPACE_DIM, k);
    FTENSOR_INDEX(SPACE_DIM, l);
 
    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_dstretch_dstress;
      t_dstretch_dstress(i, j) =
          ((t_obj_dstretch(k, l) || t_obj_dstretch(l, k)) / 2.) *
          t_inv_D(k, l, i, j);
 
      t_obj_dstress(i, j) += t_dstretch_dstress(i, j);
 
      ++t_obj_dstress;
      ++t_obj_dstretch;
      ++t_inv_D;
    }
 
    MoFEMFunctionReturn(0);
  };
 
  switch (evalEnergyModel) {
  case PYTHON_MODEL:
    CHKERR evaluate_python_objective();
    CHKERR conversion_of_biot_stress();
    if ((features & PhysicalEquations::noStretchMask).any()) {
      if ((features & variable_compliance_mask).none()) {
        CHKERR conversion_of_stretch_to_stress_for_no_stretch(
            getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(dataAtPts->matInvD));
      } else {
        CHKERR conversion_of_stretch_to_stress_for_no_stretch(
            getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->matInvD));
      }
    } else {
      CHKERR conversion_of_stretch();
    }
    break;
  case HENCKY_MODEL:
    if ((features & PhysicalEquations::noStretchMask).any()) {
      CHKERR evaluate_energy_of_hencky_model_nostreach();
      CHKERR conversion_of_biot_stress();
    } else {
      CHKERR evaluate_energy_of_hencky_model();
    }
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Objective model type not handled");
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpInteriorJImpl::integrate(EntData &data) {
  MoFEMFunctionBegin;
  locJ = 0;
 
#ifndef NDEBUG
  if (!topoData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Topological data pointer is null");
#endif // NDEBUG
 
  const int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = getGaussPts().size2();
 
  const auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_obj = topoData->getFTensorObj(nb_integration_pts);
  auto t_obj_dP = topoData->getFTensorObjDStress(nb_integration_pts);
  auto t_obj_dStrain = topoData->getFTensorObjDStrain(nb_integration_pts);
  auto t_obj_dU = topoData->getFTensorObjDDisplacement(nb_integration_pts);
  auto t_P = dataAtPts->getFTensorApproxP(nb_integration_pts);
  auto t_det = topoData->getFTensorDetJacobian(nb_integration_pts);
  auto t_inv_jac = topoData->getFTensorInvJacobian(nb_integration_pts);
  auto t_jac = topoData->getFTensorJacobian(nb_integration_pts);
 
  auto next = [&]() {
    ++t_w;
    ++t_obj;
    ++t_obj_dP;
    ++t_obj_dStrain;
    ++t_obj_dU;
    ++t_P;
    ++t_det;
    ++t_inv_jac;
    ++t_jac;
  };
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, I, J);
  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
 
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
        &v[0], &v[1], &v[2]);
  };
 
  nF.clear();
 
  const int nb_base_functions = data.getN().size2();
  auto t_base_diff = data.getFTensor1DiffN<SPACE_DIM>();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    locJ += (t_w * v * t_det) * t_obj;
 
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_cof;
    t_cof(i, j) = t_det * t_inv_jac(j, i);
 
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_dJ_dX;
    t_dJ_dX(I, J) =
 
        t_obj * t_cof(I, J)
 
        +
 
        t_obj_dP(i, j) * (t_kd(j, I) * (t_kd(k, J) * t_P(i, k)) -
                          t_inv_jac(I, J) * t_jac(j, k) * t_P(i, k));
 
    auto t_nf = get_ftensor1(nF);
    int bb = 0;
    for (; bb != nb_dofs / SPACE_DIM; ++bb) {
      t_nf(i) += (t_w * v) * t_dJ_dX(i, j) * t_base_diff(j);
      ++t_nf;
      ++t_base_diff;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_base_diff;
 
    next();
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpJ_dPImpl::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
#ifndef NDEBUG
  if (!topoData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Topological data pointer is null");
#endif // NDEBUG
 
  const int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = data.getN().size1();
 
  const auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  const int nb_base_functions = data.getN().size2() / SPACE_DIM;
  auto t_row_base_fun = data.getFTensor1N<SPACE_DIM>();
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
 
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
        &v[0], &v[1], &v[2]);
  };
 
  auto t_obj_dP = topoData->getFTensorObjDStress(nb_integration_pts);
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double a = v * t_w;
    auto t_nf = get_ftensor1(nF);
 
    int bb = 0;
    for (; bb != nb_dofs / SPACE_DIM; ++bb) {
      t_nf(i) += a * t_row_base_fun(j) * t_obj_dP(i, j);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_obj_dP;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpJ_dBubbleImpl::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
#ifndef NDEBUG
  if (!topoData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Topological data pointer is null");
#endif // NDEBUG
 
  const int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = data.getN().size1();
 
  const auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  const int nb_base_functions = data.getN().size2() / (SPACE_DIM * SPACE_DIM);
  auto t_row_base_fun = data.getFTensor2N<SPACE_DIM, SPACE_DIM>();
  auto t_obj_dP = topoData->getFTensorObjDStress(nb_integration_pts);
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
 
  auto get_ftensor0 = [](auto &v) {
    return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(&v[0]);
  };
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double a = v * t_w;
    auto t_nf = get_ftensor0(nF);
 
    int bb = 0;
    for (; bb != nb_dofs; ++bb) {
      t_nf += a * t_row_base_fun(i, j) * t_obj_dP(i, j);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_obj_dP;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpJ_dwImpl::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
#ifndef NDEBUG
  if (!topoData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Topological data pointer is null");
#endif // NDEBUG
 
  const int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = data.getN().size1();
 
  const auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  const int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  auto t_obj_dw = topoData->getFTensorObjDDisplacement(nb_integration_pts);
 
  FTENSOR_INDEX(SPACE_DIM, i);
 
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
        &v[0], &v[1], &v[2]);
  };
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double a = v * t_w;
    auto t_nf = get_ftensor1(nF);
 
    int bb = 0;
    for (; bb != nb_dofs / SPACE_DIM; ++bb) {
      t_nf(i) += a * t_row_base_fun * t_obj_dw(i);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_obj_dw;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpJ_dOmegaImpl::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
#ifndef NDEBUG
  if (!topoData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Topological data pointer is null");
#endif // NDEBUG
 
  const int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = data.getN().size1();
 
  const auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  const int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
 
  FTENSOR_INDEX(SPACE_DIM, i);
 
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
        &v[0], &v[1], &v[2]);
  };
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double a = v * t_w;
    auto t_nf = get_ftensor1(nF);
 
    int bb = 0;
    for (; bb != nb_dofs / SPACE_DIM; ++bb) {
      auto objective = topoData->getFTensorObjDRotation(nb_integration_pts,
                                                        bb == 0 ? 0 : 1, gg);
      t_nf(i) += a * t_row_base_fun * objective(i);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode dJ_duGammaImpl::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
#ifndef NDEBUG
  if (!topoData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Topological data pointer is null");
#endif // NDEBUG
 
  const int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = data.getN().size1();
 
  auto t_w = getFTensor0IntegrationWeight();
  const int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  auto t_obj_du_gamma =
      topoData->getFTensorObjDDisplacement(nb_integration_pts);
 
  FTENSOR_INDEX(SPACE_DIM, i);
 
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
        &v[0], &v[1], &v[2]);
  };
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double a = t_w * getMeasure();
    auto t_nf = get_ftensor1(nF);
 
    int bb = 0;
    for (; bb != nb_dofs / SPACE_DIM; ++bb) {
      t_nf(i) += a * t_row_base_fun * t_obj_du_gamma(i);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_obj_du_gamma;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode dJ_dTractionImpl::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
#ifndef NDEBUG
  if (!topoData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Topological data pointer is null");
#endif // NDEBUG
 
  const int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = OP::getGaussPts().size2();
 
  auto t_normal = OP::getFTensor1NormalsAtGaussPts();
  auto t_w = OP::getFTensor0IntegrationWeight();
  const int nb_base_functions = data.getN().size2() / SPACE_DIM;
  auto t_row_base_fun = data.getFTensor1N<SPACE_DIM>();
  auto t_obj_dtraction =
      topoData->getFTensorObjDDisplacement(nb_integration_pts);
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    auto t_nf = getFTensor1FromPtr<SPACE_DIM>(&*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_obj_dtraction(i) * 0.5;
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_normal;
    ++t_obj_dtraction;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpTauStabilisation_dX::integrate(int, EntityType,
                                                EntData &data) {
  MoFEMFunctionBegin;
  locJ = 0;
 
#ifndef NDEBUG
  if (!topoData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Topological data pointer is null");
  if (!brokenDispDataPtr)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Broken displacement data pointer is null");
  if (!hybridDispPtr)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Hybrid displacement pointer is null");
  if (!varHybridDispPtr)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Adjoint hybrid displacement pointer is null");
#endif // NDEBUG
 
  const int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = getGaussPts().size2();
  const int nb_base_functions = data.getN().size2();
 
#ifndef NDEBUG
  if (this->nF.size() != nb_dofs)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Size of nF %ld != nb_dofs %d", this->nF.size(), nb_dofs);
  if (data.getDiffN().size1() != nb_integration_pts)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Differential of base functions should have the same number of "
            "integration points as the data");
  if (data.getDiffN().size2() != nb_base_functions * 2)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Differential of base functions should have the same number of "
            "base functions as the data");
#endif // NDEBUG
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, I);
  FTensor::Number<0> N0;
  FTensor::Number<1> N1;
 
  auto &coords = getCoords();
  // Tau scale is based on the mesh triangle coordinates, not on the perturbed
  // material-position field. Treat h as constant with respect to X.
  const double h =
      std::get<2>(Tools::getTricircumcenter3d(coords.data().data()));
 
  auto t_w = getFTensor0IntegrationWeight();
  auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
  auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
  auto t_u_hybrid = getFTensor1FromMat<SPACE_DIM, -1, DL>(*hybridDispPtr);
  auto t_var_u_hybrid =
      getFTensor1FromMat<SPACE_DIM, -1, DL>(*varHybridDispPtr);
  auto t_diff_base = data.getFTensor1DiffN<2>();
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    FTensor::Tensor1<double, SPACE_DIM> t_normal;
    t_normal(j) =
        (FTensor::levi_civita(i, j, k) * t_tangent1(k)) * t_tangent2(i);
 
    double area = std::sqrt(t_normal(i) * t_normal(i));
    FTensor::Tensor1<double, SPACE_DIM> t_da;
    t_da(i) = t_normal(i) / area;
    area /= 2.;
    t_da(i) /= 2.;
 
    double tau_density = 0;
    for (auto &bd : *brokenDispDataPtr) {
      auto t_u_broken = getFTensor1FromMat<SPACE_DIM, -1, DL>(
          bd.getFlux(), nb_integration_pts);
      auto t_var_u_broken = getFTensor1FromMat<SPACE_DIM, -1, DL>(
          bd.getVarFlux(), nb_integration_pts);
      for (int ss = 0; ss != gg; ++ss) {
        ++t_u_broken;
        ++t_var_u_broken;
      }
 
      // This is the adjoint-weighted material derivative of the four tau
      // stabilisation residual blocks:
      //   u_gamma-u_gamma, L2-L2, u_gamma-L2, and L2-u_gamma.
      const double hybrid_hybrid = t_var_u_hybrid(i) * t_u_hybrid(i);
      const double broken_broken = t_var_u_broken(i) * t_u_broken(i);
      const double hybrid_broken = -t_var_u_hybrid(i) * t_u_broken(i);
      const double broken_hybrid = -t_var_u_broken(i) * t_u_hybrid(i);
      tau_density +=
          hybrid_hybrid + broken_broken + hybrid_broken + broken_hybrid;
    }
 
    const double tau = alphaTau / h;
    locJ += t_w * tau * area * tau_density;
 
    auto t_nf = getFTensor1FromArray<SPACE_DIM, SPACE_DIM>(nF);
    int rr = 0;
    for (; rr != nb_dofs / SPACE_DIM; ++rr) {
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_normal_dX;
      t_normal_dX(j, I) =
          (FTensor::levi_civita(i, j, I) * t_tangent2(i)) * t_diff_base(N0) +
          (FTensor::levi_civita(I, j, k) * t_tangent1(k)) * t_diff_base(N1);
 
      t_nf(I) +=
          t_w * alphaTau * tau_density * (t_da(i) * t_normal_dX(i, I)) / h;
      ++t_diff_base;
      ++t_nf;
    }
    for (; rr != nb_base_functions; ++rr)
      ++t_diff_base;
 
    ++t_w;
    ++t_tangent1;
    ++t_tangent2;
    ++t_u_hybrid;
    ++t_var_u_hybrid;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpDispBc_dX::integrate(int, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
  locJ = 0;
 
#ifndef NDEBUG
  if (!brokenSideDataPtr)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Broken side data pointer is null");
#endif // NDEBUG
 
  const int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = getGaussPts().size2();
  const int nb_base_functions = data.getN().size2();
 
#ifndef NDEBUG
  if (data.getDiffN().size1() != nb_integration_pts)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Differential of base functions should have the same number of "
            "integration points as the data");
  if (data.getDiffN().size2() != nb_base_functions * 2)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Differential of base functions should have the same number of "
            "base functions as the data");
#endif // NDEBUG
 
  double time = getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg)
    time = EshelbianCore::currentPhysicalTime;
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, I);
 
  FTensor::Number<0> N0;
  FTensor::Number<1> N1;
 
  const EntityHandle fe_ent = getFEEntityHandle();
  for (auto &bc : *bcDispPtr) {
    if (bc.faces.find(fe_ent) == bc.faces.end())
      continue;
 
    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;
    }
 
    FTensor::Tensor1<double, SPACE_DIM> t_bc_disp(bc.vals[0], bc.vals[1],
                                                  bc.vals[2]);
    t_bc_disp(i) *= scale;
 
    for (auto &bd : *brokenSideDataPtr) {
      auto t_w = getFTensor0IntegrationWeight();
      auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
      auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
      auto t_var_flux =
          getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(bd.getVarFlux());
      auto t_diff_base = data.getFTensor1DiffN<2>();
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        const double a = 0.5 * bd.getSense() * t_w;
 
        FTensor::Tensor1<double, SPACE_DIM> t_normal;
        t_normal(j) =
            (FTensor::levi_civita(i, j, k) * t_tangent1(k)) * t_tangent2(i);
 
        locJ += a * t_bc_disp(i) * (t_var_flux(i, j) * t_normal(j));
 
        auto t_nf = getFTensor1FromArray<SPACE_DIM, SPACE_DIM>(nF);
        int bb = 0;
        for (; bb != nb_dofs / SPACE_DIM; ++bb) {
          // The H(div) contravariant Piola trace P . n dA is invariant with
          // respect to material-position perturbations. The normal variation is
          // cancelled by the Piola variation of the flux, so this BC
          // contributes to locJ but not to the exact dX vector.
          ++t_nf;
          ++t_diff_base;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_diff_base;
 
        ++t_w;
        ++t_tangent1;
        ++t_tangent2;
        ++t_var_flux;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpAnalyticalDispBc_dX::integrate(int, EntityType type,
                                                EntData &data) {
  MoFEMFunctionBegin;
  locJ = 0;
 
#ifndef NDEBUG
  if (!brokenSideDataPtr)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Broken side data pointer is null");
#endif // NDEBUG
 
  const int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = getGaussPts().size2();
  const int nb_base_functions = data.getN().size2();
 
#ifndef NDEBUG
  if (data.getDiffN().size1() != nb_integration_pts)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Differential of base functions should have the same number of "
            "integration points as the data");
  if (data.getDiffN().size2() != nb_base_functions * 2)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Differential of base functions should have the same number of "
            "base functions as the data");
#endif // NDEBUG
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
 
  const EntityHandle fe_ent = getFEEntityHandle();
  for (auto &bc : *bcDispPtr) {
    if (bc.faces.find(fe_ent) == bc.faces.end())
      continue;
 
    auto v_analytical_expr = getTopologicalAnalyticalExpr(this, bc.blockName);
 
    for (auto &bd : *brokenSideDataPtr) {
      auto t_w = getFTensor0IntegrationWeight();
      auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
      auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
      auto t_var_flux =
          getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(bd.getVarFlux());
      auto t_diff_base = data.getFTensor1DiffN<2>();
      auto t_bc_disp = getFTensor1FromMat<SPACE_DIM, -1, DL>(v_analytical_expr);
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        const double a = 0.5 * bd.getSense() * t_w;
 
        FTensor::Tensor1<double, SPACE_DIM> t_normal;
        t_normal(j) =
            (FTensor::levi_civita(i, j, k) * t_tangent1(k)) * t_tangent2(i);
 
        locJ += a * t_bc_disp(i) * (t_var_flux(i, j) * t_normal(j));
 
        auto t_nf = getFTensor1FromArray<SPACE_DIM, SPACE_DIM>(nF);
        int bb = 0;
        for (; bb != nb_dofs / SPACE_DIM; ++bb) {
          // The H(div) contravariant Piola trace P . n dA is invariant with
          // respect to material-position perturbations. The normal variation is
          // cancelled by the Piola variation of the flux, so this BC
          // contributes to locJ but not to the exact dX vector.
          ++t_nf;
          ++t_diff_base;
        }
        for (; bb != nb_base_functions; ++bb)
          ++t_diff_base;
 
        ++t_w;
        ++t_tangent1;
        ++t_tangent2;
        ++t_var_flux;
        ++t_bc_disp;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpBrokenTractionBc_dX::integrate(int side, EntityType type,
                                                EntData &data) {
  MoFEMFunctionBegin;
  locJ = 0;
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, I);
  FTensor::Number<0> N0;
  FTensor::Number<1> N1;
 
  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->nF.size() != nb_dofs)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Size of nF %ld != nb_dofs %d", this->nF.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_diff_base = data.getFTensor1DiffN<2>();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    auto t_var_u_gamma = getFTensor1FromMat<SPACE_DIM>(lambdaHybridPtr);
    auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
    auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      FTensor::Tensor1<double, SPACE_DIM> t_normal;
      t_normal(j) =
          (FTensor::levi_civita(i, j, k) * t_tangent1(k)) * t_tangent2(i);
 
      double a = sqrt(t_normal(i) * t_normal(i));
      FTensor::Tensor1<double, SPACE_DIM> t_da;
      t_da(i) = t_normal(i) / a;
      a /= 2.;
      t_da(i) /= 2.;
 
      const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
      locJ -= (time_scale * t_w * a * tau) * (t_val(i) * t_var_u_gamma(i));
 
      auto t_nf = getFTensor1FromArray<SPACE_DIM, SPACE_DIM>(nF);
      int rr = 0;
      for (; rr != nb_dofs / SPACE_DIM; ++rr) {
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_normal_dX;
        t_normal_dX(j, I) =
            (FTensor::levi_civita(i, j, I) * t_tangent2(i)) * t_diff_base(N0) +
            (FTensor::levi_civita(I, j, k) * t_tangent1(k)) * t_diff_base(N1);
        t_nf(I) -= (time_scale * t_w * tau) * (t_val(i) * t_var_u_gamma(i)) *
                   (t_da(i) * t_normal_dX(i, I));
        ++t_diff_base;
        ++t_nf;
      }
      for (; rr != nb_base_functions; ++rr)
        ++t_diff_base;
 
      ++t_w;
      ++t_coords;
      ++t_var_u_gamma;
      ++t_tangent1;
      ++t_tangent2;
    }
 
    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);
      }
 
      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 OpBrokenAnalyticalTractionBc_dX::integrate(int side,
                                                          EntityType type,
                                                          EntData &data) {
  MoFEMFunctionBegin;
  locJ = 0;
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, I);
  FTensor::Number<0> N0;
  FTensor::Number<1> N1;
 
  int nb_dofs = data.getFieldData().size();
  int nb_integration_pts = getGaussPts().size2();
  int nb_base_functions = data.getN().size2();
 
#ifndef NDEBUG
  if (this->nF.size() != nb_dofs)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Size of nF %ld != nb_dofs %d", this->nF.size(), nb_dofs);
#endif // NDEBUG
 
  auto integrate_rhs = [&](auto &bc) {
    MoFEMFunctionBegin;
 
    auto v_analytical_expr = getTopologicalAnalyticalExpr(this, bc.blockName);
 
    auto t_val = getFTensor1FromMat<SPACE_DIM, -1, DL>(v_analytical_expr);
    auto t_diff_base = data.getFTensor1DiffN<2>();
    auto t_w = getFTensor0IntegrationWeight();
 
    auto t_var_u_gamma = getFTensor1FromMat<SPACE_DIM>(lambdaHybridPtr);
    auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
    auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      FTensor::Tensor1<double, SPACE_DIM> t_normal;
      t_normal(j) =
          (FTensor::levi_civita(i, j, k) * t_tangent1(k)) * t_tangent2(i);
 
      double a = sqrt(t_normal(i) * t_normal(i));
      FTensor::Tensor1<double, SPACE_DIM> t_da;
      t_da(i) = t_normal(i) / a;
      a /= 2.;
      t_da(i) /= 2.;
 
      locJ -= (t_w * a) * (t_val(i) * t_var_u_gamma(i));
 
      auto t_nf = getFTensor1FromArray<SPACE_DIM, SPACE_DIM>(nF);
      int rr = 0;
      for (; rr != nb_dofs / SPACE_DIM; ++rr) {
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_normal_dX;
        t_normal_dX(j, I) =
            (FTensor::levi_civita(i, j, I) * t_tangent2(i)) * t_diff_base(N0) +
            (FTensor::levi_civita(I, j, k) * t_tangent1(k)) * t_diff_base(N1);
        t_nf(I) -=
            t_w * (t_val(i) * t_var_u_gamma(i)) * (t_da(i) * t_normal_dX(i, I));
        ++t_diff_base;
        ++t_nf;
      }
      for (; rr != nb_base_functions; ++rr)
        ++t_diff_base;
 
      ++t_w;
      ++t_val;
      ++t_var_u_gamma;
      ++t_tangent1;
      ++t_tangent2;
    }
 
    MoFEMFunctionReturn(0);
  };
 
  EntityHandle fe_ent = getFEEntityHandle();
  for (auto &bc : *(bcData)) {
    if (bc.faces.find(fe_ent) != bc.faces.end() && nb_dofs) {
      CHKERR integrate_rhs(bc);
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpJ_dUImpl::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
#ifndef NDEBUG
  if (!topoData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Topological data pointer is null");
#endif // NDEBUG
 
  const int nb_dofs = data.getIndices().size();
  if (!nb_dofs)
    MoFEMFunctionReturnHot(0);
 
  const int nb_integration_pts = data.getN().size1();
 
  const auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  const int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
 
  auto t_obj_dlog_stretch = topoData->getFTensorObjDStrain(nb_integration_pts);
 
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTensor::Index<'L', size_symm> L;
  auto t_L = FTensor::SymmLTensor<double, 3>();
 
  auto get_ftensor1 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, size_symm>, size_symm>(
        &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
  };
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    const double a = v * t_w;
    auto t_nf = get_ftensor1(OP::nF);
 
    FTensor::Tensor1<double, size_symm> t_obj_dU;
    t_obj_dU(L) = t_obj_dlog_stretch(i, j) * t_L(i, j, L);
 
    int bb = 0;
    for (; bb != nb_dofs / size_symm; ++bb) {
      t_nf(L) += a * t_row_base_fun * t_obj_dU(L);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_obj_dlog_stretch;
  }
 
  MoFEMFunctionReturn(0);
}
 
struct OpSensitivityInteriorGradient
    : public OpAssembleVolumeTopologicalDerivativeImpl {
  OpSensitivityInteriorGradient(
      const std::string field_name,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<TopologicalData> topo_ptr,
      SmartPetscObj<Vec> assemble_vec, const double alpha, const double rho,
      const double alpha_viscous_omega = 0,
      boost::shared_ptr<double> J_ptr = nullptr)
      : OpAssembleVolumeTopologicalDerivativeImpl(
            field_name, data_ptr, topo_ptr, J_ptr, assemble_vec, Tag()),
        alphaW(alpha), alphaRho(rho), alphaViscousOmega(alpha_viscous_omega) {}
 
  MoFEMErrorCode integrate(EntData &data) {
    MoFEMFunctionBegin;
#ifndef NDEBUG
    if (EshelbianCore::l2UserBaseScale) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "L2 user base scale is set to %d, current inplementation only "
              "hanlde case for false",
              EshelbianCore::l2UserBaseScale);
    }
#endif // NDEBUG
 
#ifndef NDEBUG
    if (!dataAtPts)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "DataAtIntegrationPts pointer is null");
#endif // NDEBUG
 
    const int nb_dofs = data.getIndices().size();
    if (!nb_dofs)
      MoFEMFunctionReturnHot(0);
 
    const int nb_integration_pts = getGaussPts().size2();
 
    const auto v = getVolume();
    auto t_w = getFTensor0IntegrationWeight();
    // OpSpatialEquilibrium
    auto t_div_P = dataAtPts->getFTensorDivP(nb_integration_pts);
    auto t_var_w = dataAtPts->getFTensorVarWL2(nb_integration_pts);
    // OpSpatialRotation
    auto t_approx_P = dataAtPts->getFTensorApproxP(nb_integration_pts);
    auto t_u_h1 = dataAtPts->getFTensorStretchH1(nb_integration_pts);
    auto t_var_omega = dataAtPts->getFTensorVarRotAxis(nb_integration_pts);
    // OpSpatialConsistencyP
    auto t_h = dataAtPts->getFTensorSmallH(nb_integration_pts);
    auto t_var_P = dataAtPts->getFTensorVarPiola(nb_integration_pts);
    auto t_w_l2 = dataAtPts->getFTensorSmallWL2(nb_integration_pts);
    auto t_var_div_P = dataAtPts->getFTensorDivVarPiola(nb_integration_pts);
 
    auto t_jac = topoData->getFTensorJacobian(nb_integration_pts);
 
    auto get_ftensor1 = [](auto &v) {
      return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
          &v[0], &v[1], &v[2]);
    };
 
    auto next = [&]() {
      ++t_w;
      ++t_div_P;
      ++t_var_w;
      ++t_approx_P;
      ++t_u_h1;
      ++t_var_omega;
      ++t_h;
      ++t_var_P;
      ++t_w_l2;
      ++t_var_div_P;
      ++t_jac;
    };
 
    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, I);
    FTENSOR_INDEX(SPACE_DIM, J);
 
    constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
 
    locJ = 0;
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      double a = v * t_w;
 
      // The adjoint direction varies both element-constant and higher modes.
      FTensor::Tensor2<double, 3, 3> variation_R;
      if (EshelbianCore::rotSelector == LARGE_ROT) {
        using namespace LogarithmicDirectionalDerivatives;
        const auto read = [&](MatrixPtr values) {
          Vector result;
          auto view = GetFTensor1FromMatImpl<3, -1, DL, MatrixDouble>::get(
              *values, gg, 0);
          result(i) = view(i);
          return result;
        };
        const SplitRotation rotation(read(dataAtPts->constantRotation),
                                     read(dataAtPts->higherRotation));
        variation_R(i, j) =
            rotation.first(read(dataAtPts->varConstantRotation),
                           read(dataAtPts->varHigherRotation))(i, j);
      } else {
        variation_R(i, j) = levi_civita(i, j, k) * t_var_omega(k);
      }
 
      FTensor::Tensor4<double, SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM>
          t_diff;
      t_diff(i, j, l, k) = t_kd(i, l) * t_kd(j, k);
 
      // OpSpatialEqulibrium
      {
        const auto beta = a * (t_div_P(i) * t_var_w(i));
        locJ -= beta;
      }
      // OpSpatialEquilibrium
      {
        double beta;
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_beta_dX;
        switch (EshelbianCore::gradApproximator) {
        case NO_H1_CONFIGURATION:
        case LARGE_ROT:
        case MODERATE_ROT:
          beta =
 
              (variation_R(j, k) * t_u_h1(k, l))
 
              * (t_approx_P(j, n) * t_jac(l, n));
 
          t_beta_dX(I, J) =
 
              (variation_R(j, k) * t_u_h1(k, l))
 
              * (t_approx_P(j, n) * t_diff(l, n, I, J));
 
          break;
        case SMALL_ROT:
          beta =
 
              variation_R(i, j)
 
              * (t_approx_P(i, n) * t_jac(j, n));
 
          t_beta_dX(I, J) =
 
              variation_R(i, j)
 
              * (t_approx_P(i, n) * t_diff(j, n, I, J));
 
          break;
        default:
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "gradApproximator not handled");
          break;
        };
 
        locJ -= a * beta;
        auto t_nf = get_ftensor1(nF);
        auto t_base_diff = data.getFTensor1DiffN<3>(gg, 0);
        for (int bb = 0; bb != nb_dofs / SPACE_DIM; ++bb) {
          t_nf(i) -= a * (t_beta_dX(i, j) * t_base_diff(j));
          ++t_nf;
          ++t_base_diff;
        }
      }
      // OpSpatialConsistency
      {
        const auto beta =
            (t_h(i, j) - t_kd(i, j)) * (t_var_P(i, n) * t_jac(j, n));
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_beta_dX;
        t_beta_dX(I, J) =
            (t_h(i, j) - t_kd(i, j)) * (t_var_P(i, n) * t_diff(j, n, I, J));
        locJ -= a * beta;
        auto t_nf = get_ftensor1(nF);
        auto t_base_diff = data.getFTensor1DiffN<3>(gg, 0);
        for (int bb = 0; bb != nb_dofs / SPACE_DIM; ++bb) {
          t_nf(i) -= a * (t_beta_dX(i, j) * t_base_diff(j));
          ++t_nf;
          ++t_base_diff;
        }
      }
      // OpSpatialConsistency, cont
      {
        const auto beta = t_w_l2(i) * t_var_div_P(i);
        locJ -= a * beta;
      }
 
      next();
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  const double alphaW;
  const double alphaRho;
  const double alphaViscousOmega;
};
 
struct OpSensitivityInterior_dX
    : public OpAssembleVolumeTopologicalDerivativeImpl {
  OpSensitivityInterior_dX(const std::string &field_name,
                           boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                           SmartPetscObj<Vec> assemble_vec,
                           boost::shared_ptr<TopologicalData> topo_ptr,
                           const double alpha, const double rho,
                           const double alpha_viscous_omega = 0)
      : OpAssembleVolumeTopologicalDerivativeImpl(
            field_name, data_ptr, topo_ptr, nullptr, assemble_vec, Tag()),
        alphaW(alpha), alphaRho(rho) {
    if (alpha_viscous_omega)
      CHK_THROW_MESSAGE(
          MOFEM_NOT_IMPLEMENTED,
          "OpSensitivityInterior_dX with alpha_viscous_omega != 0 is not "
          "implemented yet");
  }
 
  MoFEMErrorCode integrate(EntData &data) {
    MoFEMFunctionBegin;
 
#ifndef NDEBUG
    if (!dataAtPts)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "DataAtIntegrationPts pointer is null");
#endif // NDEBUG
 
    const int nb_dofs = data.getIndices().size();
    if (!nb_dofs)
      MoFEMFunctionReturnHot(0);
 
    const int nb_integration_pts = getGaussPts().size2();
 
    const auto v = getVolume();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_div_P = dataAtPts->getFTensorDivP(nb_integration_pts);
    auto t_w_l2 = dataAtPts->getFTensorSmallWL2(nb_integration_pts);
    auto t_s_dot_w = dataAtPts->getFTensorSmallWL2Dot(nb_integration_pts);
    auto t_s_dot_dot_w =
        dataAtPts->getFTensorSmallWL2DotDot(nb_integration_pts);
    auto t_h = dataAtPts->getFTensorSmallH(nb_integration_pts);
    auto t_levi_kirchhoff =
        dataAtPts->getFTensorLeviKirchhoff(nb_integration_pts);
    auto t_omega_grad_dot =
        dataAtPts->getFTensorRotAxisGradDot(nb_integration_pts);
    auto t_R = dataAtPts->getFTensorRotMat(nb_integration_pts);
    auto t_u = dataAtPts->getFTensorStretch(nb_integration_pts);
 
    auto t_var_w = dataAtPts->getFTensorVarWL2(nb_integration_pts);
    auto t_var_grad_omega =
        dataAtPts->getFTensorVarGradRotAxis(nb_integration_pts);
    auto t_var_P = dataAtPts->getFTensorVarPiola(nb_integration_pts);
    auto t_var_div_P = dataAtPts->getFTensorDivVarPiola(nb_integration_pts);
    auto t_det = topoData->getFTensorDetJacobian(nb_integration_pts);
    auto t_inv_jac = topoData->getFTensorInvJacobian(nb_integration_pts);
 
    auto w_l2_dot_dot_at_pts = dataAtPts->getSmallWL2DotDotAtPts();
    if (w_l2_dot_dot_at_pts->size1() != nb_integration_pts ||
        w_l2_dot_dot_at_pts->size2() != SPACE_DIM) {
      MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
          *w_l2_dot_dot_at_pts, nb_integration_pts);
      w_l2_dot_dot_at_pts->clear();
    }
 
    const auto piola_scale = dataAtPts->piolaScale;
    const auto alpha_w = alphaW / piola_scale;
    const auto alpha_rho = alphaRho / piola_scale;
 
    const int nb_base_functions = data.getN().size2();
    auto t_base_diff = data.getFTensor1DiffN<3>();
 
    auto get_ftensor1 = [](auto &v) {
      return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
          &v[0], &v[1], &v[2]);
    };
 
    auto next = [&]() {
      ++t_w;
      ++t_div_P;
      ++t_w_l2;
      ++t_s_dot_w;
      ++t_s_dot_dot_w;
      ++t_h;
      ++t_levi_kirchhoff;
      ++t_omega_grad_dot;
      ++t_R;
      ++t_u;
 
      ++t_var_w;
      ++t_var_grad_omega;
      ++t_var_P;
      ++t_var_div_P;
 
      ++t_det;
      ++t_inv_jac;
    };
 
    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);
    constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
      // Calculate the variation of the gradient due to geometry change
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_cof;
      t_cof(i, j) = t_det * t_inv_jac(j, i);
 
      double consistency_residual = 0;
      if (EshelbianCore::gradApproximator == NO_H1_CONFIGURATION) {
        consistency_residual = -0.5 * t_var_P(k, m) * (t_R(k, l) * t_u(l, m)) -
                               0.5 * t_var_P(k, l) * (t_R(k, m) * t_u(l, m)) +
                               t_var_P(k, l) * t_kd(k, l);
      } else {
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_residuum_P;
        t_residuum_P(k, m) = t_h(k, m) - t_kd(k, m);
        consistency_residual = t_var_P(k, m) * (-t_residuum_P(k, m));
      }
 
      auto t_nf = get_ftensor1(nF);
      int bb = 0;
      for (; bb != nb_dofs / SPACE_DIM; ++bb) {
 
        FTensor::Tensor1<double, SPACE_DIM> t_div_base;
        t_div_base(i) = -(1 / t_det) * (t_inv_jac(j, i) * t_base_diff(j));
 
        // OpSpatialEquilibrium
        t_nf(i) += (t_w * v) *
                   (t_var_w(k) * (-t_div_P(k) + alpha_w * t_s_dot_w(k) +
                                  alpha_rho * t_s_dot_dot_w(k))) *
                   t_cof(i, j) * t_base_diff(j);
        t_nf(i) += (t_w * v) * (-(t_var_w(k) * t_div_P(k))) * t_div_base(i);
 
        // OpSpatialRotation
        auto constant_work =
            dataAtPts->getFTensorLeviKirchhoff(nb_integration_pts, 0, gg);
        auto constant_variation =
            GetFTensor1FromMatImpl<3, -1, DL, MatrixDouble>::get(
                *dataAtPts->varConstantRotation, gg, 0);
        auto higher_variation =
            GetFTensor1FromMatImpl<3, -1, DL, MatrixDouble>::get(
                *dataAtPts->varHigherRotation, gg, 0);
        const double rotation_work = constant_variation(k) * constant_work(k) +
                                     higher_variation(k) * t_levi_kirchhoff(k);
        t_nf(i) -= (t_w * v * rotation_work) * t_cof(i, j) * t_base_diff(j);
 
        // OpSpatialConsistencyP
        t_nf(i) +=
            (t_w * v * consistency_residual) * t_cof(i, j) * t_base_diff(j);
 
        // OpSpatialConsistencyDivTerm
        t_nf(i) += (t_w * v) * (t_var_div_P(k) * (-t_w_l2(k))) * t_cof(i, j) *
                   t_base_diff(j);
        t_nf(i) += (t_w * v) * (t_var_div_P(k) * (-t_w_l2(k))) * t_div_base(i);
 
        ++t_nf;
        ++t_base_diff;
      }
      for (; bb != nb_base_functions; ++bb)
        ++t_base_diff;
 
      next();
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  const double alphaW;
  const double alphaRho;
};
 
struct OpBodyForce_dX : public OpAssembleVolumeTopologicalDerivativeImpl {
  OpBodyForce_dX(MoFEM::Interface &m_field, int ms_id,
                 const std::string &field_name,
                 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                 SmartPetscObj<Vec> assemble_vec,
                 boost::shared_ptr<TopologicalData> topo_ptr,
                 std::vector<boost::shared_ptr<ScalingMethod>> smv,
                 boost::shared_ptr<double> J_ptr = nullptr)
      : OpAssembleVolumeTopologicalDerivativeImpl(
            field_name, data_ptr, topo_ptr, J_ptr, assemble_vec, Tag()),
        scalingMethods(smv) {
    CHK_THROW_MESSAGE(getMeshsetData(m_field, ms_id), "Get meshset data");
  }
 
  MoFEMErrorCode integrate(EntData &data) {
    MoFEMFunctionBegin;
 
    locJ = 0;
 
    if (entsPtr) {
      if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
        MoFEMFunctionReturnHot(0);
    }
 
#ifndef NDEBUG
    if (!dataAtPts)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "DataAtIntegrationPts pointer is null");
#endif // NDEBUG
 
    const int nb_dofs = data.getIndices().size();
    if (!nb_dofs)
      MoFEMFunctionReturnHot(0);
 
    const int nb_integration_pts = getGaussPts().size2();
 
    const auto v = getVolume();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_var_w_l2 = dataAtPts->getFTensorVarWL2(nb_integration_pts);
    auto t_det = topoData->getFTensorDetJacobian(nb_integration_pts);
    auto t_inv_jac = topoData->getFTensorInvJacobian(nb_integration_pts);
 
    const int nb_base_functions = data.getN().size2();
    auto t_base_diff = data.getFTensor1DiffN<3>();
 
    auto get_ftensor1 = [](auto &v) {
      return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
          &v[0], &v[1], &v[2]);
    };
 
    auto next = [&]() {
      ++t_w;
      ++t_var_w_l2;
      ++t_inv_jac;
      ++t_det;
    };
 
    FTENSOR_INDEXES(SPACE_DIM, i, j, I, J);
 
    auto get_scale = [&](const double t) {
      double s = 1;
      for (auto &o : scalingMethods) {
 
        s *= o->getScale(t);
      }
      return s;
    };
 
    auto scale = get_scale(getFEMethod()->ts_t);
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      const double alpha = scale * t_w * v;
 
      double adjount = t_var_w_l2(i) * tForce(i);
      locJ += (alpha * t_det) * adjount;
 
      FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_cof;
      t_cof(i, j) = t_det * t_inv_jac(j, i);
 
      auto t_nf = get_ftensor1(nF);
      int bb = 0;
      for (; bb != nb_dofs / SPACE_DIM; ++bb) {
        t_nf(i) += (alpha * adjount) * t_cof(i, J) * t_base_diff(J);
 
        ++t_nf;
        ++t_base_diff;
      }
      for (; bb != nb_base_functions; ++bb)
        ++t_base_diff;
 
      next();
    }
 
    MoFEMFunctionReturn(0);
  }
 
protected:
  MoFEMErrorCode getMeshsetData(MoFEM::Interface &m_field, int ms_id) {
    MoFEMFunctionBegin;
 
    auto cubit_meshset_ptr =
        m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(ms_id,
                                                                    BLOCKSET);
 
    std::vector<double> block_data;
    CHKERR cubit_meshset_ptr->getAttributes(block_data);
 
    if (block_data.size() != SPACE_DIM) {
      MOFEM_LOG("SELF", Sev::warning)
          << "BLOCKSET is expected to have " << SPACE_DIM
          << " attributes but has size " << block_data.size();
      if (block_data.size() < SPACE_DIM) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Size of attribute in BLOCKSET is too small");
      }
    }
 
    for (unsigned int ii = 0; ii != SPACE_DIM; ++ii) {
      tForce(ii) = block_data[ii];
    }
 
    MOFEM_LOG("WORLD", Sev::noisy)
        << "Flux blockset " << cubit_meshset_ptr->getName();
    MOFEM_LOG("WORLD", Sev::noisy)
        << "Number of attributes " << block_data.size();
 
    this->entsPtr = boost::make_shared<Range>();
    CHKERR m_field.get_moab().get_entities_by_handle(cubit_meshset_ptr->meshset,
                                                     *(entsPtr), true);
 
    MOFEM_LOG("WORLD", Sev::noisy) << "tForce vector initialised: " << tForce;
    MOFEM_LOG("WORLD", Sev::noisy) << "Number of elements " << entsPtr->size();
 
    MoFEMFunctionReturn(0);
  }
 
  FTensor::Tensor1<double, SPACE_DIM> tForce;
  boost::shared_ptr<Range> entsPtr;
  std::vector<boost::shared_ptr<ScalingMethod>> scalingMethods;
};
 
} // namespace EshelbianPlasticity