EshelbianTopologicalDerivative.cpp

Go to the documentation of this file.

/**
 * @file EshelbianTopologicalDerivative.cpp
 * @author your name (you@domain.com)
 * @brief 
 * @version 0.1
 * @date 2026-02-21
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
#define SINGULARITY
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
#include <EshelbianAux.hpp>
#include <TSElasticPostStep.hpp>
 
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/lambert_w.hpp>
 
#include <EshelbianTopologicalDerivative.hpp>
#include <ObjectiveFunctionData.hpp>
 
using namespace EshelbianPlasticity;
using namespace ShapeOptimization;
 
namespace EshelbianPlasticity {
 
struct TopologicalTAOCtxImpl : public TopologicalTAOCtx {
  TopologicalTAOCtxImpl(
      EshelbianCore *ep,
      ForcesAndSourcesCore::GaussHookFun set_integration_at_interior,
      ForcesAndSourcesCore::GaussHookFun set_integration_at_face,
      SmartPetscObj<TS> time_solver)
      : ep_ptr(ep), integrationAtInterior(set_integration_at_interior),
        integrationAtFace(set_integration_at_face), timeSolver(time_solver),
 
        geometryVec(CommInterface::createEntitiesPetscVector(
            ep->mField.get_comm(), ep->mField.get_moab(), 0, 3, Sev::noisy)),
        gradJVec(CommInterface::createEntitiesPetscVector(
            ep->mField.get_comm(), ep->mField.get_moab(), 0, 3, Sev::noisy)) {}
 
private:
  EshelbianCore *ep_ptr;
  ForcesAndSourcesCore::GaussHookFun integrationAtInterior;
  ForcesAndSourcesCore::GaussHookFun integrationAtFace;
  SmartPetscObj<TS> timeSolver;
 
  CommInterface::EntitiesPetscVector geometryVec;
  CommInterface::EntitiesPetscVector gradJVec;
 
  auto getMaterialH1PositionsStructureRawPtr() {
    return ep_ptr->mField.get_field_structure(ep_ptr->materialH1Positions);
  }
 
  auto getVertexMaterialH1PositionsTag() {
    auto field_ptr = getMaterialH1PositionsStructureRawPtr();
    return field_ptr->th_FieldDataVerts;
  }
 
  Tag getTag(moab::Interface &moab, std::string tag_name, int size) {
    Tag tag;
    rval = moab.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE, tag,
                               MB_TAG_CREAT | MB_TAG_DENSE, NULL);
    if (rval == MB_ALREADY_ALLOCATED)
      rval = MB_SUCCESS;
    CHK_MOAB_THROW(rval, "Failed to get tag " + tag_name);
    return tag;
  };
 
  Tag getJdXTag() { return getTag(ep_ptr->mField.get_moab(), "DJ_DX", 3); }
 
  friend MoFEMErrorCode cohesiveEvaluateObjectiveAndGradient(Tao tao, Vec x,
                                                             PetscReal *f,
                                                             Vec g, void *ctx);
 
  friend MoFEMErrorCode testTopologicalDerivative(TopologicalTAOCtx *ctx_ptr,
                                                  const double alpha,
                                                  const double rho,
                                                  const double alpha_omega);
};
 
boost::shared_ptr<TopologicalTAOCtx> createTopologicalTAOCtx(
    EshelbianCore *ep,
    ForcesAndSourcesCore::GaussHookFun set_integration_at_interior,
    ForcesAndSourcesCore::GaussHookFun set_integration_at_face,
    SmartPetscObj<TS> time_solver) {
  return boost::make_shared<TopologicalTAOCtxImpl>(
      ep, set_integration_at_interior, set_integration_at_face, time_solver);
}
}
 
#include <impl/EshelbianTopologicalDerivativeOperators.cpp>
#include <impl/EshelbianTopologicalDerivativePush.cpp>
 
namespace EshelbianPlasticity {
 
MoFEMErrorCode testTopologicalDerivative(TopologicalTAOCtx *ctx_ptr,
                                         const double alpha, const double rho,
                                         const double alpha_omega) {
 
  MoFEMFunctionBegin;
  auto *ctx_impl_ptr = dynamic_cast<TopologicalTAOCtxImpl *>(ctx_ptr);
  if (!ctx_impl_ptr)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Invalid context pointer type");
  auto &ep = *ctx_impl_ptr->ep_ptr;
  auto &m_field = ep.mField;
  auto interior_integration_hook = ctx_impl_ptr->integrationAtInterior;
  auto boundary_integration_hook = ctx_impl_ptr->integrationAtFace;
  auto th_top_tag = ctx_impl_ptr->getVertexMaterialH1PositionsTag();
  auto th_grad_tag = ctx_impl_ptr->getJdXTag(m_field.get_moab());
 
  boost::shared_ptr<double> J_ptr = boost::make_shared<double>(0.);
  auto dJ_dX_vec = ctx_impl_ptr->gradJVec.second;
 
  auto fe_spatial =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
  auto dJ_dx_vec = pushTopologicalSpatialOps(
      ep, fe_spatial, interior_integration_hook, boundary_integration_hook);
  auto fe_material =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
  CHKERR pushTopologicalMaterialOps(ep, fe_material, interior_integration_hook,
                                    boundary_integration_hook, J_ptr,
                                    th_grad_tag);
 
  auto lambda_vec = createDMVector(ep.dmElastic);
  auto ksp = snesGetKSP(tsGetSNES(ctx_impl_ptr->timeSolver));
  CHKERR KSPSolveTranspose(ksp, dJ_dx_vec, lambda_vec);
 
  auto fe_adjoint =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
  auto dJ_adjoint_vec = pushTopologicalOps_dJ_adjont_gradient(
      ep, fe_adjoint, interior_integration_hook, boundary_integration_hook,
      lambda_vec, dJ_dX_vec, alpha, rho, alpha_omega);
 
  MoFEMFunctionReturn(0);
}
 
}