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) {
 
    primalProblemVec = createDMVector(ep->dmElastic, RowColData::COL);
    adjointProblemVec = createDMVector(ep->dmElastic, RowColData::COL);
  }
 
private:
  EshelbianCore *ep_ptr;
  ForcesAndSourcesCore::GaussHookFun integrationAtInterior;
  ForcesAndSourcesCore::GaussHookFun integrationAtFace;
  SmartPetscObj<TS> timeSolver;
 
  friend MoFEMErrorCode
  evaluatePrimalProblemTopologicalImpl(TopologicalTAOCtxImpl *ctx_impl_ptr);
 
  friend MoFEMErrorCode
  evaluateAdjointProblemImpl(TopologicalTAOCtxImpl *ctx_impl_ptr);
 
  friend MoFEMErrorCode
  evaluateGradientImpl(TopologicalTAOCtxImpl *ctx_impl_ptr, double *f, Vec g);
 
  friend MoFEMErrorCode
  finiteDifferenceGradientTest(TopologicalTAOCtxImpl *ctx_impl_ptr, Vec sol,
                               double *f, Vec g, double epsilon);
 
  friend MoFEMErrorCode
  topologicalEvaluateObjectiveAndGradient(Tao tao, Vec sol, PetscReal *f, Vec g,
                                          void *ctx);
 
  friend MoFEMErrorCode testTopologicalDerivative(TopologicalTAOCtx *ctx_ptr,
                                                  Vec sol, PetscReal *f,
                                                  Vec g);
 
  SmartPetscObj<Vec> primalProblemVec;
  SmartPetscObj<Vec> adjointProblemVec;
};
 
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
evaluatePrimalProblemTopologicalImpl(TopologicalTAOCtxImpl *ctx_impl_ptr) {
  MoFEMFunctionBegin;
  auto &ep = *ctx_impl_ptr->ep_ptr;
  auto ts = ctx_impl_ptr->timeSolver;
 
  CHKERR VecGhostUpdateBegin(ctx_impl_ptr->primalProblemVec, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(ctx_impl_ptr->primalProblemVec, INSERT_VALUES,
                           SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, ctx_impl_ptr->primalProblemVec,
                                 INSERT_VALUES, SCATTER_FORWARD,
                                 RowColData::ROW);
 
  CHKERR TSSetFromOptions(ts);
  CHKERR TSSetStepNumber(ts, 0);
  CHKERR TSSetTime(ts, 0);
  CHKERR TSSetSolution(ts, ctx_impl_ptr->primalProblemVec);
  double dt;
  CHKERR TSGetTimeStep(ts, &dt);
  MOFEM_LOG("EP", Sev::inform)
      << "evaluatePrimalProblemTopologicalImpl: Time step dt: " << dt;
  CHKERR TSSolve(ts, PETSC_NULLPTR);
  CHKERR TSSetTimeStep(ts, dt);
 
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode evaluateAdjointProblemImpl(TopologicalTAOCtxImpl *ctx_impl_ptr) {
  MoFEMFunctionBegin;
  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 fe_spatial =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
  auto dJ_dx_vec = pushTopologicalSpatialOps(
      ep, fe_spatial, interior_integration_hook, boundary_integration_hook);
  CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName,
                                  fe_spatial);
  CHKERR VecAssemblyBegin(dJ_dx_vec);
  CHKERR VecAssemblyEnd(dJ_dx_vec);
  CHKERR VecGhostUpdateBegin(dJ_dx_vec, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateEnd(dJ_dx_vec, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateBegin(dJ_dx_vec, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(dJ_dx_vec, INSERT_VALUES, SCATTER_FORWARD);
 
  auto ksp = snesGetKSP(tsGetSNES(ctx_impl_ptr->timeSolver));
  CHKERR KSPSolveTranspose(ksp, dJ_dx_vec, ctx_impl_ptr->adjointProblemVec);
  CHKERR VecGhostUpdateBegin(ctx_impl_ptr->adjointProblemVec, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(ctx_impl_ptr->adjointProblemVec, INSERT_VALUES,
                           SCATTER_FORWARD);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode evaluateGradientImpl(TopologicalTAOCtxImpl *ctx_impl_ptr,
                                    double *f, Vec g) {
  MoFEMFunctionBegin;
  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 fe_material =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
  auto fe_adjoint =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
 
  boost::shared_ptr<double> J_ptr(fe_material, f);
  auto dJ_dX_vec = SmartPetscObj<Vec>(g, true);
  CHKERR VecZeroEntries(dJ_dX_vec);
  CHKERR VecGhostUpdateBegin(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
 
  CHKERR pushTopologicalMaterialOps(ep, fe_material, interior_integration_hook,
                                    boundary_integration_hook, J_ptr,
                                    dJ_dX_vec);
  auto alpha = ep.alphaW;
  auto rho = ep.alphaRho;
  auto alpha_omega = ep.alphaOmega;
  CHKERR pushTopologicalOps_dJ_adjoint_gradient(
      ep, fe_adjoint, interior_integration_hook, boundary_integration_hook,
      ctx_impl_ptr->adjointProblemVec, dJ_dX_vec, alpha, rho, alpha_omega);
 
  CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
                                  fe_material);
  CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
                                  fe_adjoint);
  CHKERR VecAssemblyBegin(dJ_dX_vec);
  CHKERR VecAssemblyEnd(dJ_dX_vec);
  CHKERR VecGhostUpdateBegin(dJ_dX_vec, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateEnd(dJ_dX_vec, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateBegin(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode finiteDifferenceGradientTest(TopologicalTAOCtxImpl *ctx_impl_ptr,
                                            Vec sol, double *f, Vec g,
                                            double epsilon) {
  MoFEMFunctionBegin;
 
  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 fe_material =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
  auto J_ptr = boost::make_shared<double>(0);
  CHKERR pushTopologicalMaterialOps(ep, fe_material, interior_integration_hook,
                                    boundary_integration_hook, J_ptr,
                                    SmartPetscObj<Vec>());
 
  auto opt = ep.mField.getInterface<OperatorsTester>();
  constexpr double eps = 1e-6;
  auto random_vec = opt->setRandomFields(
      ep.dmMaterial, {{ep.materialH1Positions, {-eps, eps}}}, nullptr,
      RowColData::ROW);
  auto a_vec = vectorDuplicate(random_vec);
  CHKERR VecCopy(sol, a_vec);
  CHKERR VecAXPY(a_vec, 1, random_vec);
  auto b_vec = vectorDuplicate(random_vec);
  CHKERR VecCopy(sol, b_vec);
  CHKERR VecAXPY(b_vec, -1, random_vec);
 
  *J_ptr = 0;
  CHKERR VecGhostUpdateBegin(a_vec, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(a_vec, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(ep.dmMaterial, a_vec, INSERT_VALUES,
                                 SCATTER_REVERSE, RowColData::ROW);
  CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
                                  fe_material);
  double J_plus = *J_ptr;
  *J_ptr = 0;
  CHKERR VecGhostUpdateBegin(b_vec, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(b_vec, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(ep.dmMaterial, b_vec, INSERT_VALUES,
                                 SCATTER_REVERSE, RowColData::ROW);
  CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
                                  fe_material);
  double J_minus = *J_ptr;
  double dJ_da = (J_plus - J_minus) / (2 * eps);
 
  double exact_dJ;
  CHKERR VecDot(g, random_vec, &exact_dJ);
  double error = std::abs(dJ_da - exact_dJ);
 
  MOFEM_LOG("EP", Sev::inform)
      << "J = " << *f << ", dJ/da = " << dJ_da << ", exact dJ/da = " << exact_dJ
      << ", error = " << error;
 
  CHKERR VecGhostUpdateBegin(sol, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(sol, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(ep.dmMaterial, sol, INSERT_VALUES,
                                 SCATTER_REVERSE, RowColData::ROW);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode testTopologicalDerivative(TopologicalTAOCtx *ctx_ptr, Vec sol,
                                         PetscReal *f, Vec g) {
 
  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;
 
  CHKERR VecGhostUpdateBegin(sol, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(sol, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(ep.dmMaterial, sol, INSERT_VALUES,
                                 SCATTER_REVERSE, RowColData::ROW);
 
  CHKERR VecZeroEntries(g);
  *f = 0;
 
  CHKERR evaluatePrimalProblemTopologicalImpl(ctx_impl_ptr);
  double norm2_x;
  CHKERR VecNorm(ctx_impl_ptr->primalProblemVec, NORM_2, &norm2_x);
  MOFEM_LOG("EP", Sev::inform)
      << "evaluatePrimalProblemTopologicalImpl: Norm of displacement vector: "
      << norm2_x;
  if (norm2_x < 1e-12) {
    MoFEMFunctionReturnHot(0);
  }
 
  CHKERR evaluateAdjointProblemImpl(ctx_impl_ptr);
  CHKERR evaluateGradientImpl(ctx_impl_ptr, f, g);
 
  CHKERR finiteDifferenceGradientTest(ctx_impl_ptr, sol, f, g, 1e-6);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode topologicalEvaluateObjectiveAndGradient(Tao tao, Vec sol,
                                                       PetscReal *f, Vec g,
                                                       void *ctx) {
 
  MoFEMFunctionBegin;
//   auto *ctx_impl_ptr = static_cast<TopologicalTAOCtxImpl *>(ctx);
//   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 alpha = ep.alphaW;
//   auto rho = ep.alphaRho;
//   auto alpha_omega = ep.alphaOmega;
 
 
  MoFEMFunctionReturn(0);
}
}