EshelbianTopologicalDerivativePush.cpp

Go to the documentation of this file.

/**
 * @file EshelbianTopologicalDerivativePush.cpp
 * @author your name (you@domain.com)
 * @brief 
 * @version 0.1
 * @date 2026-02-25
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
namespace EshelbianPlasticity {
 
static auto pushInteriorTopologicalOpsImpl(
    EshelbianCore &ep, boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe,
    ForcesAndSourcesCore::GaussHookFun interior_integration_hook,
    SmartPetscObj<Vec> lambda_vec = SmartPetscObj<Vec>()) {
 
  // This calculates field data in the interior, used by other operators.
 
  auto data_at_pts_ptr = boost::make_shared<DataAtIntegrationPts>();
 
  auto bubble_cache =
      boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0, MatrixDouble());
  fe->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      fe->getOpPtrVector(), {HDIV, H1, L2}, ep.materialH1Positions,
      ep.frontAdjEdges);
 
  // set integration rule
  fe->getRuleHook = [](int, int, int) { return -1; };
  fe->setRuleHook = interior_integration_hook;
 
  // calculate fields values
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
      ep.piolaStress, data_at_pts_ptr->getApproxPAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
      ep.bubbleField, data_at_pts_ptr->getApproxPAtPts(), MBMAXTYPE));
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
      ep.piolaStress, data_at_pts_ptr->getDivPAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      ep.rotAxis, data_at_pts_ptr->getRotAxisAtPts(), MBTET));
 
  if (EshelbianCore::isNoStretch()) {
    pushOpCalculateStretchFromStress(
        fe->getOpPtrVector(), ep.physicalEquations, data_at_pts_ptr,
        ep.externalStrainVecPtr, ep.timeScaleMap);
  } else {
    fe->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValues<3>(
            ep.stretchTensor, data_at_pts_ptr->getLogStretchTensorAtPts(),
            MBTET));
  }
  CHK_THROW_MESSAGE(VecSetDM(ep.solTSStep, PETSC_NULLPTR), "VecSetDM failed");
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
      ep.piolaStress, data_at_pts_ptr->getApproxP0AtPts(), nullptr,
      ep.solTSStep));
  if (!EshelbianCore::isNoStretch()) {
    fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
        ep.stretchTensor, data_at_pts_ptr->getLogStretchTensor0AtPts(),
        ep.solTSStep, MBTET));
  }
 
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      ep.rotAxis, data_at_pts_ptr->getRotAxis0AtPts(), ep.solTSStep, MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      ep.spatialL2Disp, data_at_pts_ptr->getSmallWL2AtPts(), MBTET));
 
  // H1 displacements
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      ep.spatialH1Disp, data_at_pts_ptr->getSmallWH1AtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
      ep.spatialH1Disp, data_at_pts_ptr->getSmallWGradH1AtPts()));
  fe->getOpPtrVector().push_back(
      new OpCalculateRotationAndSpatialGradient(data_at_pts_ptr));
 
  if (EshelbianCore::isNoStretch()) {
  } else {
    fe->getOpPtrVector().push_back(ep.physicalEquations->returnOpJacobian(
        false, false, data_at_pts_ptr, ep.physicalEquations));
  }
 
  if (lambda_vec) {
    fe->getOpPtrVector().push_back(
        new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
            ep.piolaStress, data_at_pts_ptr->getVarPiolaPts(), nullptr,
            lambda_vec));
    fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
        ep.bubbleField, data_at_pts_ptr->getVarPiolaPts(), nullptr, lambda_vec,
        MBMAXTYPE));
    fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
        ep.piolaStress, data_at_pts_ptr->getDivVarPiolaPts(), lambda_vec));
 
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        ep.spatialL2Disp, data_at_pts_ptr->getVarWL2Pts(), lambda_vec, MBTET));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        ep.rotAxis, data_at_pts_ptr->getVarRotAxisPts(), lambda_vec, MBTET));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
        ep.rotAxis, data_at_pts_ptr->getVarGradRotAxisPts(), lambda_vec, MBTET));
 
    if (EshelbianCore::isNoStretch()) {
      fe->getOpPtrVector().push_back(
          ep.physicalEquations->returnOpCalculateVarStretchFromStress(
              data_at_pts_ptr, ep.physicalEquations));
    } else {
      fe->getOpPtrVector().push_back(
          new OpCalculateTensor2SymmetricFieldValues<3>(
              ep.stretchTensor, data_at_pts_ptr->getVarLogStreachPts(),
              lambda_vec, MBTET));
    }
  }
 
  return data_at_pts_ptr;
}
 
static auto pushBoundaryTopologicalOpsImpl(
    EshelbianCore &ep, const std::string &fe_name,
    boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe,
    ForcesAndSourcesCore::GaussHookFun boundary_integration_hook) {
 
  // This calculates field data on the boundary, used by other operators.
 
  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;
 
  // First: Iterate over skeleton FEs adjacent to Domain FEs
  // Note: BoundaryEle, i.e. uses the skeleton integration rule
  auto op_loop_skeleton_side =
      new OpLoopSide<BoundaryEle>(m_field, fe_name, SPACE_DIM - 1, Sev::noisy);
  op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
    return -1;
  };
  op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
      boundary_integration_hook;
 
  CHKERR
  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      op_loop_skeleton_side->getOpPtrVector(), {L2}, ep.materialH1Positions,
      ep.frontAdjEdges);
 
  // Second: iterate over domain FEs adjacent to the skeleton, in particular
  // the domain element.
  auto broken_data_ptr = boost::make_shared<std::vector<BrokenBaseSideData>>();
  // 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);
  CHKERR
  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
      ep.materialH1Positions, ep.frontAdjEdges);
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpGetBrokenBaseSideData<SideEleOp>(ep.piolaStress, broken_data_ptr));
  auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(ep.piolaStress,
                                                           flux_mat_ptr));
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpSetFlux<SideEleOp>(broken_data_ptr, flux_mat_ptr));
  op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
 
  return std::make_tuple(op_loop_skeleton_side,
                         op_loop_domain_side->getSideFEPtr(), broken_data_ptr);
}
 
static auto pushInteriorTopological_dJ_dx_Impl(
    EshelbianCore &ep, boost::shared_ptr<ForcesAndSourcesCore> fe,
    boost::shared_ptr<DataAtIntegrationPts> data_at_pts_ptr,
    boost::shared_ptr<TopologicalData> topo_ptr) {
 
  // Calculate the interior objective derivative for the adjoint problem.
 
  auto dJ_dx_vec = createDMVector(ep.dmElastic);
 
  fe->getOpPtrVector().push_back(
      new OpJ_dPImpl(ep.piolaStress, data_at_pts_ptr, topo_ptr, nullptr,
                     dJ_dx_vec, Tag()));
  fe->getOpPtrVector().push_back(new OpJ_dBubbleImpl(
      ep.bubbleField, data_at_pts_ptr, topo_ptr, nullptr, dJ_dx_vec, Tag()));
  if (!EshelbianCore::isNoStretch()) {
    fe->getOpPtrVector().push_back(new OpJ_dUImpl(
        ep.stretchTensor, data_at_pts_ptr, topo_ptr, nullptr, dJ_dx_vec,
        Tag()));
  }
 
  // If we explicitly add rotation to J, add the corresponding rotational
  // derivative operator.
 
  fe->getOpPtrVector().push_back(
      new OpJ_dwImpl(ep.spatialL2Disp, data_at_pts_ptr, topo_ptr, nullptr,
                     dJ_dx_vec, Tag()));
  fe->getOpPtrVector().push_back(new OpJ_dOmegaImpl(
      ep.rotAxis, data_at_pts_ptr, topo_ptr, nullptr, dJ_dx_vec, Tag()));
 
  return dJ_dx_vec;
}
 
static auto pushBoundaryTopological_dJ_dx_Impl(
    EshelbianCore &ep, boost::shared_ptr<ForcesAndSourcesCore> fe,
    boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_data_ptr,
    boost::shared_ptr<DataAtIntegrationPts> data_at_pts_ptr,
    boost::shared_ptr<TopologicalData> topo_ptr, SmartPetscObj<Vec> dJ_dx_vec) {
  MoFEMFunctionBegin;
 
  // Calculate the boundary objective derivative for the adjoint problem.
 
  fe->getOpPtrVector().push_back(new dJ_duGammaImpl(
      ep.hybridSpatialDisp, data_at_pts_ptr, topo_ptr, nullptr, dJ_dx_vec,
      Tag()));
  fe->getOpPtrVector().push_back(
      new dJ_dTractionImpl(broken_data_ptr, topo_ptr, nullptr, dJ_dx_vec,
                           Tag()));
  MoFEMFunctionReturn(0);
}
 
SmartPetscObj<Vec> pushTopologicalSpatialOps(
    EshelbianCore &ep, boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe,
    ForcesAndSourcesCore::GaussHookFun interior_integration_hook,
    ForcesAndSourcesCore::GaussHookFun boundary_integration_hook,
    ObjectiveModelType eval_energy_model) {
 
  auto data_at_pts_ptr =
      pushInteriorTopologicalOpsImpl(ep, fe, interior_integration_hook);
  auto [op_skeleton_side, domain_side_fe_ptr, broken_data_ptr] =
      pushBoundaryTopologicalOpsImpl(ep, ep.skinElement, fe,
                                     boundary_integration_hook);
  fe->getOpPtrVector().push_back(op_skeleton_side);
 
  auto topo_ptr = boost::make_shared<TopologicalData>();
 
  switch (eval_energy_model) {
  case ObjectiveModelType::PYTHON_MODEL: {
    char objective_function_file_name[255] = "objective_function.py";
    CHKERR PetscOptionsGetString(
        PETSC_NULLPTR, PETSC_NULLPTR, "-objective_function",
        objective_function_file_name, 255, PETSC_NULLPTR);
    auto python_ptr =
        create_python_objective_function(objective_function_file_name);
    fe->getOpPtrVector().push_back(new OpTopologicalObjectivePythonImpl(
        data_at_pts_ptr, topo_ptr, python_ptr, eval_energy_model));
  } break;
  case ObjectiveModelType::HENCKY_MODEL: {
    fe->getOpPtrVector().push_back(new OpTopologicalObjectivePythonImpl(
        data_at_pts_ptr, topo_ptr, nullptr, eval_energy_model));
  } break;
  default:
    CHK_THROW_MESSAGE(MOFEM_INVALID_DATA, "Unknown objective model type");
  }
 
  auto dJ_dx_vec =
      pushInteriorTopological_dJ_dx_Impl(ep, fe, data_at_pts_ptr, topo_ptr);
 
  //   Assume for (temporarily) that we do not have boundary contributions to
  //   dJ/dx. We add this once interior contributions are working and verified,
  //   and once we have a test case with non-zero boundary contribution. CHKERR
  //   pushBoundaryTopological_dJ_dx_Impl(
  //       ep, fe, broken_data_ptr, data_at_pts_ptr, topo_ptr, dJ_dx_vec);
 
  return dJ_dx_vec;
}
 
MoFEMErrorCode pushTopologicalMaterialOps(
    EshelbianCore &ep, boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe,
    ForcesAndSourcesCore::GaussHookFun interior_integration_hook,
    ForcesAndSourcesCore::GaussHookFun boundary_integration_hook,
    boost::shared_ptr<double> J_ptr, SmartPetscObj<Vec> dJ_dX_vec,
    ObjectiveModelType eval_energy_model) {
  MoFEMFunctionBegin;
 
  // This is derivative of dJ_dX, i.e. derivative with respect to material
  // position, used by other operators.
 
  auto data_at_pts_ptr =
      pushInteriorTopologicalOpsImpl(ep, fe, interior_integration_hook);
  auto [op_skeleton_side, domain_side_fe_ptr, broken_data_ptr] =
      pushBoundaryTopologicalOpsImpl(ep, ep.skinElement, fe,
                                     boundary_integration_hook);
  fe->getOpPtrVector().push_back(op_skeleton_side);
 
  auto topo_ptr = boost::make_shared<TopologicalData>();
  switch (eval_energy_model) {
  case ObjectiveModelType::PYTHON_MODEL: {
    char objective_function_file_name[255] = "objective_function.py";
    CHKERR PetscOptionsGetString(
        PETSC_NULLPTR, PETSC_NULLPTR, "-objective_function",
        objective_function_file_name, 255, PETSC_NULLPTR);
    auto python_ptr =
        create_python_objective_function(objective_function_file_name);
    fe->getOpPtrVector().push_back(new OpTopologicalObjectivePythonImpl(
        data_at_pts_ptr, topo_ptr, python_ptr, eval_energy_model));
  } break;
  case ObjectiveModelType::HENCKY_MODEL: {
    fe->getOpPtrVector().push_back(new OpTopologicalObjectivePythonImpl(
        data_at_pts_ptr, topo_ptr, nullptr, eval_energy_model));
  } break;
  default:
    CHK_THROW_MESSAGE(MOFEM_INVALID_DATA, "Unknown objective model type");
  }
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
      ForcesAndSourcesCore::UserDataOperator::OPROW, ep.materialH1Positions,
      topo_ptr->getJacobianAtPts(), SmartPetscObj<Vec>()));
  fe->getOpPtrVector().push_back(new OpInvertMatrix<3>(
      topo_ptr->getJacobianAtPts(), topo_ptr->getDetJacobianAtPts(),
      topo_ptr->getInvJacobianAtPts()));
  fe->getOpPtrVector().push_back(new OpInteriorJImpl(ep.materialH1Positions,
                                                     data_at_pts_ptr, topo_ptr,
                                                     J_ptr, dJ_dX_vec, Tag()));
 
  // The boundary derivative is still missing here. The objective function has
  // two integrals, one in the interior and one on the boundary. Only the
  // interior integral is currently implemented, so OpBoundaryJImpl still needs
  // to be added.
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode pushTopologicalInteriorOps_dJ_adjoint_gradient(
    EshelbianCore &ep, boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe,
    ForcesAndSourcesCore::GaussHookFun interior_integration_hook,
    ForcesAndSourcesCore::GaussHookFun boundary_integration_hook,
    SmartPetscObj<Vec> lambda_vec, SmartPetscObj<Vec> dJ_dX_vec,
    const double alpha, const double rho, const double alpha_viscous_omega,
    boost::shared_ptr<double> J_ptr,
    SmartPetscObj<Vec> topo_vec = SmartPetscObj<Vec>()) {
 
  MoFEMFunctionBegin;
 
  auto data_at_pts_ptr = pushInteriorTopologicalOpsImpl(
      ep, fe, interior_integration_hook, lambda_vec);
 
  auto add_hybridised_dJ_gradient = [&](auto fe_ptr) {
    MoFEMFunctionBegin;
    using BoundaryEle =
        PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
    using BdyEleOp = BoundaryEle::UserDataOperator;
 
    auto [op_skeleton_side, domain_side_fe_ptr, broken_data_ptr] =
        pushBoundaryTopologicalOpsImpl(ep, ep.skeletonElement, fe_ptr,
                                       boundary_integration_hook);
 
    using EleOnSide =
        PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
    using SideEleOp = EleOnSide::UserDataOperator;
 
    auto var_piola_mat_ptr = boost::make_shared<MatrixDouble>();
    domain_side_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
            ep.piolaStress, var_piola_mat_ptr, nullptr, lambda_vec));
    domain_side_fe_ptr->getOpPtrVector().push_back(
        new OpSetVarFlux<SideEleOp>(broken_data_ptr, var_piola_mat_ptr));
 
    auto broken_disp_data_ptr =
        boost::make_shared<std::vector<BrokenBaseSideData>>();
    domain_side_fe_ptr->getOpPtrVector().push_back(
        new OpGetBrokenBaseSideData<SideEleOp>(ep.spatialL2Disp,
                                               broken_disp_data_ptr));
    auto disp_mat_ptr = boost::make_shared<MatrixDouble>();
    domain_side_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>(ep.spatialL2Disp,
                                                    disp_mat_ptr));
    domain_side_fe_ptr->getOpPtrVector().push_back(
        new OpSetFlux<SideEleOp>(broken_disp_data_ptr, disp_mat_ptr));
    auto var_disp_mat_ptr = boost::make_shared<MatrixDouble>();
    domain_side_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>(
            ep.spatialL2Disp, var_disp_mat_ptr, lambda_vec));
    domain_side_fe_ptr->getOpPtrVector().push_back(
        new OpSetVarFlux<SideEleOp>(broken_disp_data_ptr, var_disp_mat_ptr));
 
    op_skeleton_side->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>(
            ep.hybridSpatialDisp, data_at_pts_ptr->getVarHybridDispAtPts(),
            lambda_vec));
    op_skeleton_side->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>(
            ep.hybridSpatialDisp, data_at_pts_ptr->getHybridDispAtPts()));
    auto topo_ptr = boost::make_shared<TopologicalData>();
    op_skeleton_side->getOpPtrVector().push_back(new OpGetHONormalsOnFace(
        ep.materialH1Positions, topo_ptr->getTangent1AtPts(),
        topo_ptr->getTangent2AtPts(), topo_ptr->getNormalAtPts()));
 
    if (ep.alphaTau > 0.0) {
      op_skeleton_side->getOpPtrVector().push_back(new OpTauStabilisation_dX(
          ep.materialH1Positions, broken_disp_data_ptr,
          data_at_pts_ptr->getHybridDispAtPts(),
          data_at_pts_ptr->getVarHybridDispAtPts(), topo_ptr, ep.alphaTau,
          dJ_dX_vec, J_ptr));
    }
 
    using OpTopoC_dHybrid =
        FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
            GAUSS>::OpTopoDerivativeBrokenSpaceConstrainDHybrid<SPACE_DIM>;
    using OpTopoC_dFlux =
        FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
            GAUSS>::OpTopoDerivativeBrokenSpaceConstrainDFlux<SPACE_DIM>;
    op_skeleton_side->getOpPtrVector().push_back(new OpTopoC_dHybrid(
        ep.materialH1Positions, broken_data_ptr,
        data_at_pts_ptr->getVarHybridDispAtPts(), topo_ptr->getTangent1AtPts(),
        topo_ptr->getTangent2AtPts(), boost::make_shared<double>(1.0),
        dJ_dX_vec, Tag(), J_ptr));
    op_skeleton_side->getOpPtrVector().push_back(new OpTopoC_dFlux(
        ep.materialH1Positions, broken_data_ptr,
        data_at_pts_ptr->getHybridDispAtPts(), topo_ptr->getTangent1AtPts(),
        topo_ptr->getTangent2AtPts(), boost::make_shared<double>(1.0),
        dJ_dX_vec, Tag(), J_ptr));
 
    fe_ptr->getOpPtrVector().push_back(op_skeleton_side);
    MoFEMFunctionReturn(0);
  };
 
  CHK_THROW_MESSAGE(add_hybridised_dJ_gradient(fe),
                    "Failed to add hybridised dJ/dx gradient operators");
 
  auto topo_ptr = boost::make_shared<TopologicalData>();
  if (topo_vec) {
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
        ForcesAndSourcesCore::UserDataOperator::OPROW, ep.materialH1Positions,
        topo_ptr->getJacobianAtPts(), topo_vec));
  } else {
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
        ep.materialH1Positions, topo_ptr->getJacobianAtPts()));
  }
 
  fe->getOpPtrVector().push_back(new OpInvertMatrix<3>(
      topo_ptr->getJacobianAtPts(), topo_ptr->getDetJacobianAtPts(),
      topo_ptr->getInvJacobianAtPts()));
 
  fe->getOpPtrVector().push_back(new OpSensitivityInteriorGradient(
      ep.materialH1Positions, data_at_pts_ptr, topo_ptr, dJ_dX_vec, alpha, rho,
      alpha_viscous_omega, J_ptr));
  if (auto op = ep.physicalEquations->returnOpTopoSpatialPhysical(
          ep.materialH1Positions, data_at_pts_ptr, dJ_dX_vec, topo_ptr,
          ep.alphaU, J_ptr)) {
    fe->getOpPtrVector().push_back(op);
  }
 
  std::string body_force_history;
  CHKERR ep.getStringArgumentFromJsonBlocksets("BODY_FORCE", "load_history",
                                               body_force_history);
  if (body_force_history.empty()) {
    body_force_history = "body_force.txt";
  } else {
    MOFEM_LOG("EP", Sev::inform)
        << "Body force load history from JSON: " << body_force_history;
  }
  auto body_time_scale =
      boost::make_shared<EshelbianCore::DynamicRelaxationTimeScale>(
          body_force_history);
  auto add_body_force_opv = [&](auto &&meshset_vec_ptr) {
    for (auto m : meshset_vec_ptr) {
      auto op = new OpBodyForce_dX(
          ep.mField, m->getMeshsetId(), ep.materialH1Positions, data_at_pts_ptr,
          dJ_dX_vec, topo_ptr,
          std::vector<boost::shared_ptr<ScalingMethod>>{body_time_scale},
          J_ptr);
      fe->getOpPtrVector().push_back(op);
    }
  };
 
  const std::string body_force_block_name = "BODY_FORCE";
 
  add_body_force_opv(
 
      ep.mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
          (boost::format("%s(.*)") % body_force_block_name).str()
 
              ))
 
  );
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode pushTopologicalBoundaryOps_dJ_adjoint_gradient(
    EshelbianCore &ep, boost::shared_ptr<FaceElementForcesAndSourcesCore> fe,
    ForcesAndSourcesCore::GaussHookFun interior_integration_hook,
    ForcesAndSourcesCore::GaussHookFun boundary_integration_hook,
    SmartPetscObj<Vec> lambda_vec, SmartPetscObj<Vec> dJ_dX_vec,
    boost::shared_ptr<double> J_ptr,
    SmartPetscObj<Vec> topo_vec = SmartPetscObj<Vec>()) {
  MoFEMFunctionBegin;
 
  if (!lambda_vec)
    SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_NULL,
            "Lambda vector is required for boundary dJ/dx gradient");
 
  fe->getRuleHook = [](int, int, int) { return -1; };
  fe->setRuleHook = boundary_integration_hook;
  CHKERR
  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      fe->getOpPtrVector(), {L2}, ep.materialH1Positions, ep.frontAdjEdges);
 
  auto get_broken_op_side = [&ep, lambda_vec](auto &pip) {
    using EleOnSide =
        PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
    using SideEleOp = EleOnSide::UserDataOperator;
    // Iterate over domain FEs adjacent to boundary.
    auto broken_data_ptr =
        boost::make_shared<std::vector<BrokenBaseSideData>>();
    // Note: EleOnSide, i.e. uses on domain projected skeleton rule
    auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
        ep.mField, ep.elementVolumeName, SPACE_DIM, Sev::noisy);
    op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
        boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
    CHKERR
    EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
        ep.materialH1Positions, ep.frontAdjEdges);
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpGetBrokenBaseSideData<SideEleOp>(ep.piolaStress,
                                               broken_data_ptr));
    auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(ep.piolaStress,
                                                             flux_mat_ptr));
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpSetFlux<SideEleOp>(broken_data_ptr, flux_mat_ptr));
    auto var_flux_mat_ptr = boost::make_shared<MatrixDouble>();
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
            ep.piolaStress, var_flux_mat_ptr, nullptr, lambda_vec));
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpSetVarFlux<SideEleOp>(broken_data_ptr, var_flux_mat_ptr));
    pip.push_back(op_loop_domain_side);
    return broken_data_ptr;
  };
 
  auto broken_data_ptr = get_broken_op_side(fe->getOpPtrVector());
  auto topo_ptr = boost::make_shared<TopologicalData>();
 
  auto lambda_hybrid_dip_ptr = boost::make_shared<MatrixDouble>();
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
      ep.hybridSpatialDisp, lambda_hybrid_dip_ptr, lambda_vec));
  fe->getOpPtrVector().push_back(new OpDispBc_dX(
      ep.materialH1Positions, broken_data_ptr, ep.bcSpatialDispVecPtr,
      ep.timeScaleMap, topo_ptr, dJ_dX_vec, J_ptr));
  fe->getOpPtrVector().push_back(new OpAnalyticalDispBc_dX(
      ep.materialH1Positions, broken_data_ptr,
      ep.bcSpatialAnalyticalDisplacementVecPtr, ep.timeScaleMap, topo_ptr,
      dJ_dX_vec, J_ptr));
  fe->getOpPtrVector().push_back(new OpBrokenTractionBc_dX(
      ep.materialH1Positions, ep.bcSpatialTractionVecPtr, lambda_hybrid_dip_ptr,
      topo_ptr, ep.timeScaleMap, dJ_dX_vec, J_ptr));
  fe->getOpPtrVector().push_back(new OpBrokenAnalyticalTractionBc_dX(
      ep.materialH1Positions, ep.bcSpatialAnalyticalTractionVecPtr,
      lambda_hybrid_dip_ptr, topo_ptr, ep.timeScaleMap, dJ_dX_vec, J_ptr));
 
  // Boundary tau stabilisation still needs a material derivative if enabled.
  // The skeleton jump term is handled in
  // pushTopologicalInteriorOps_dJ_adjoint_gradient; boundary displacement and
  // rotation penalties need their own BC-specific derivative operators.
 
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity