EshelbianTopologicalDerivativePush.cpp

Go to the documentation of this file.

/**
 * @file EshelbianTopologicalPush.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>()) {
 
  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()));
 
  if (ep.noStretch) {
    fe->getOpPtrVector().push_back(
        ep.physicalEquations->returnOpCalculateStretchFromStress(
            data_at_pts_ptr, ep.physicalEquations));
  } else {
    fe->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValues<3>(
            ep.stretchTensor, data_at_pts_ptr->getLogStretchTensorAtPts(),
            MBTET));
  }
 
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      ep.rotAxis, data_at_pts_ptr->getRotAxisAtPts(), 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()));
 
  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 (ep.noStretch) {
      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) {
 
  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 skeleton interation 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 skelton, particularly one
  // 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>();
  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) {
 
  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::noStretch) {
    fe->getOpPtrVector().push_back(new OpJ_dP_no_streachImpl(
        ep.piolaStress, data_at_pts_ptr, topo_ptr, nullptr, dJ_dx_vec,
        Tag()));
    fe->getOpPtrVector().push_back(new OpJ_dBubble_no_streachImpl(
        ep.bubbleField, data_at_pts_ptr, topo_ptr, nullptr, dJ_dx_vec,
        Tag()));
  } else {
    fe->getOpPtrVector().push_back(new OpJ_dUImpl(
        ep.stretchTensor, data_at_pts_ptr, topo_ptr, nullptr, dJ_dx_vec,
        Tag()));
  }
  fe->getOpPtrVector().push_back(
      new OpJ_dwImpl(ep.spatialL2Disp, 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;
  // fe->getOpPtrVector().push_back(new dJ_duGammaImpl(
  //     ep.hybridSpatialDisp, data_at_pts_ptr, topo_ptr, nullptr, dJ_dx_vec));
  // fe->getOpPtrVector().push_back(
  //     new dJ_dTractionImpl(broken_data_ptr, topo_ptr, nullptr, dJ_dx_vec));
  MoFEMFunctionReturn(0);
}
 
SmartPetscObj<Vec> pushTopologicalSpatialOps(
    EshelbianCore &ep, boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe,
    ForcesAndSourcesCore::GaussHookFun interior_integration_hook,
    ForcesAndSourcesCore::GaussHookFun boundary_integration_hook) {
 
  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>();
  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));
 
  auto dJ_dx_vec =
      pushInteriorTopological_dJ_dx_Impl(ep, fe, data_at_pts_ptr, topo_ptr);
 
  // Assume for now 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, Tag th_grad_tag) {
  MoFEMFunctionBegin;
 
  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>();
  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));
  fe->getOpPtrVector().push_back(
      new OpInteriorJImpl(ep.materialH1Positions, data_at_pts_ptr, topo_ptr,
                          J_ptr, SmartPetscObj<Vec>(), th_grad_tag));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode pushTopologicalOps_dJ_adjont_gradient(
    EshelbianCore &ep, boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe,
    ForcesAndSourcesCore::GaussHookFun interior_integration_hook,
    ForcesAndSourcesCore::GaussHookFun boundary_integration_hook,
    SmartPetscObj<Vec> lambda_vec, Tag th_grad_tag,
    const double alpha, const double rho, const double alpha_omega) {
 
  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);
    domain_side_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
            ep.piolaStress, data_at_pts_ptr->getVarPiolaPts(), nullptr,
            lambda_vec));
    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()));
 
    using OpTopoC_dHybrid = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
          GAUSS>::OpTopoDerivativeBrokenSpaceConstrainDHybrid<SPACE_DIM>;
    using OpTopoC_dBroken = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
          GAUSS>::OpTopoDerivativeBrokenSpaceConstrainDFlux<SPACE_DIM>;
    op_skeleton_side->getOpPtrVector().push_back(new OpTopoC_dHybrid(
        ep.hybridSpatialDisp, broken_data_ptr,
        data_at_pts_ptr->getVarHybridDispAtPts(), topo_ptr->getTangent1AtPts(),
        topo_ptr->getTangent2AtPts(), boost::make_shared<double>(1.0),
        SmartPetscObj<Vec>(), th_grad_tag));
    op_skeleton_side->getOpPtrVector().push_back(new OpTopoC_dBroken(
        broken_data_ptr, data_at_pts_ptr->getHybridDispAtPts(),
        data_at_pts_ptr->getVarHybridDispAtPts(), topo_ptr->getTangent1AtPts(),
        topo_ptr->getTangent2AtPts(), boost::make_shared<double>(1.0),
        SmartPetscObj<Vec>(), th_grad_tag));
 
    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>();
  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, th_grad_tag, alpha,
      rho, alpha_omega));
 
  MoFEMFunctionReturn(0);
}
    
}