EshelbianTopologicalDerivative.hpp

Go to the documentation of this file.

/**
 * @file EshelbianTopologicalDerivative.hpp
 * @brief 
 * @version 0.1
 * @date 2026-02-11
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
#ifndef ESHELBIAN_TOPOLOGICAL_DERIVATIVE_HPP
#define ESHELBIAN_TOPOLOGICAL_DERIVATIVE_HPP
 
namespace EshelbianPlasticity {
 
struct MatOpsTopoData
    : public boost::enable_shared_from_this<MatOpsTopoData> {
 
  MatOpsTopoData() {
    matOpsData = MatOps::createMatOpsDataPtr();
    matOpsData->insertCommonData("jacobianAtPts");
    matOpsData->insertCommonData("invJacobianAtPtr");
    matOpsData->insertCommonData("detJacobianAtPts");
    matOpsData->insertCommonData("tangent1AtPts");
    matOpsData->insertCommonData("tangent2AtPts");
    matOpsData->insertCommonData("normalAtPts");
    matOpsData->insertCommonData("objAtPts");
    matOpsData->insertCommonData("objDStrainAtPts");
    matOpsData->insertCommonData("objDStressAtPts");
    matOpsData->insertCommonData("objDRotationAtPts");
    matOpsData->insertCommonData("objDDisplacementAtPts");
  }
 
  boost::shared_ptr<MatOps::MatOpsData> matOpsData;
};
 
enum ObjectiveModelType { PYTHON_MODEL, HENCKY_MODEL, LAST_MODEL };
 
struct TopologicalData : public MatOpsTopoData {
 
  using MatOpsTopoData::MatOpsTopoData;
 
  template <class Tensor>
  static inline auto getFTensor(MatrixDouble &m, const int nb_gauss_pts) {
    return MatrixSizeHelper<Tensor, DL>::get(m, nb_gauss_pts)();
  }
 
  template <class Tensor>
  static inline auto getFTensor(MatrixPtr m, const int nb_gauss_pts) {
    return getFTensor<Tensor>(*m, nb_gauss_pts);
  }
 
  static inline auto getFTensor0(MatrixPtr m, const int nb_gauss_pts) {
#ifndef NDEBUG
    if (!m)
      THROW_MESSAGE("TopologicalData matrix pointer is not available");
 
    const auto nb = static_cast<size_t>(nb_gauss_pts);
    if (!((m->size1() == nb && m->size2() == 1) ||
          (m->size1() == 1 && m->size2() == nb)))
      THROW_MESSAGE("Scalar matrix size is not compatible with the number of "
                    "integration points");
#endif
    return getFTensor0FromMat(*m);
  }
 
  inline auto getFTensorJacobian(const int nb_gauss_pts) {
    return getFTensor<
        GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("jacobianAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorInvJacobian(const int nb_gauss_pts) {
    return getFTensor<
        GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("invJacobianAtPtr"), nb_gauss_pts);
  }
 
  inline auto getFTensorDetJacobian(const int nb_gauss_pts) {
    return getFTensor0(matOpsData->getCommonDataPtr("detJacobianAtPts"),
                       nb_gauss_pts);
  }
 
  inline auto getFTensorTangent1(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("tangent1AtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorTangent2(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("tangent2AtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorNormal(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("normalAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorObj(const int nb_gauss_pts) {
    return getFTensor0(matOpsData->getCommonDataPtr("objAtPts"),
                       nb_gauss_pts);
  }
 
  inline auto getFTensorObjDStrain(const int nb_gauss_pts) {
    return getFTensor<
        GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("objDStrainAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorObjDStress(const int nb_gauss_pts) {
    return getFTensor<
        GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("objDStressAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorObjDDisplacement(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("objDDisplacementAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorObjDRotation(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("objDRotationAtPts"), nb_gauss_pts);
  }
 
  inline auto getJacobianAtPts() {
    return matOpsData->getCommonDataPtr("jacobianAtPts");
  }
 
  inline auto getInvJacobianAtPts() {
    return matOpsData->getCommonDataPtr("invJacobianAtPtr");
  }
 
  inline auto getDetJacobianAtPts() {
    return matOpsData->getCommonDataPtr("detJacobianAtPts");
  }
 
  inline auto getTangent1AtPts() {
    return matOpsData->getCommonDataPtr("tangent1AtPts");
  }
 
  inline auto getTangent2AtPts() {
    return matOpsData->getCommonDataPtr("tangent2AtPts");
  }
 
  inline auto getNormalAtPts() {
    return matOpsData->getCommonDataPtr("normalAtPts");
  }
 
  inline auto getObjAtPts() {
    return matOpsData->getCommonDataPtr("objAtPts");
  }
 
  inline auto getObjDStrainAtPts() {
    return matOpsData->getCommonDataPtr("objDStrainAtPts");
  }
 
  inline auto getObjDStressAtPts() {
    return matOpsData->getCommonDataPtr("objDStressAtPts");
  }
 
  inline auto getObjDDisplacementAtPts() {
    return matOpsData->getCommonDataPtr("objDDisplacementAtPts");
  }
 
  inline auto getObjDRotationAtPts() {
    return matOpsData->getCommonDataPtr("objDRotationAtPts");
  }
 
};
 
struct TopologicalTAOCtx {
  TopologicalTAOCtx() = default;
  virtual ~TopologicalTAOCtx() = default;
};
 
 
boost::shared_ptr<TopologicalTAOCtx> createTopologicalTAOCtx(
    EshelbianCore *ep,
    ForcesAndSourcesCore::GaussHookFun set_integration_at_interior,
    ForcesAndSourcesCore::GaussHookFun set_integration_at_face,
    SmartPetscObj<TS> time_solver);
 
MoFEMErrorCode topologicalEvaluateObjectiveAndGradient(Tao tao, Vec sol,
                                                       PetscReal *f, Vec g,
                                                       void *ctx);
 
MoFEMErrorCode evaluateGradient(TopologicalTAOCtx *ctx_ptr, double *f, Vec g,
                                ObjectiveModelType eval_energy_model);
 
MoFEMErrorCode testTopologicalDerivative(TopologicalTAOCtx *ctx_ptr, Vec sol,
                                         PetscReal *f, Vec g,
                                         ObjectiveModelType eval_energy_model);
}
 
#endif // ESHELBIAN_TOPOLOGICAL_DERIVATIVE_HPP