ObjectiveFunctionData.cpp

Go to the documentation of this file.

/**
 * @file ObjectiveFunctionData.cpp
 * @brief Implementation of Python-integrated objective function interface
 * @details Implements Python-C++ bridge for objective function evaluation
 *          using Boost.Python and NumPy for topology optimization
 */
 
#ifdef ENABLE_PYTHON_BINDING
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include <ObjectiveFunctionData.hpp>
 
namespace ShapeOptimization {
 
#ifdef ENABLE_PYTHON_BINDING
 
constexpr int SPACE_DIM =
    EXECUTABLE_DIMENSION; ///< Space dimension of problem (2D or 3D), set at compile time
 
/**
 * @brief Implementation of ObjectiveFunctionData interface using Python integration
 * 
 * This class provides a concrete implementation of the ObjectiveFunctionData interface
 * that bridges MoFEM C++ data structures with Python-defined objective functions.
 * It enables flexible definition of optimization objectives through Python scripting
 * while maintaining high-performance computation in the C++ finite element framework.
 * 
 * Key features:
 * - Python function evaluation with automatic data conversion
 * - Support for objective function and gradient computations
 * - Topology optimization mode definition through Python
 * - Efficient NumPy array interfacing for large data sets
 * - Automatic memory management between C++ and Python
 * 
 * The class handles:
 * 1. Loading and executing Python objective function scripts
 * 2. Converting MoFEM data structures to NumPy arrays
 * 3. Calling Python functions for objective evaluation
 * 4. Converting Python results back to MoFEM format
 * 5. Managing Python interpreter state and namespace
 * 
 * Example Python interface functions that must be defined:
 * - objectiveInteriorFunction(coords, u, stress, strain) -> objective_value
 * - objectiveInteriorGradientStress(coords, u, stress, strain) -> gradient_wrt_stress
 * - numberOfModes(block_id) -> number_of_topology_modes
 * - blockModes(block_id, coords, centroid, bbox) -> mode_vectors
 */
struct ObjectiveFunctionDataImpl : public ObjectiveFunctionData {
  ObjectiveFunctionDataImpl() = default;
  virtual ~ObjectiveFunctionDataImpl() = default;
 
  /// Initialize Python interpreter and load objective function script
  MoFEMErrorCode initPython(const std::string py_file);
 
  /**
   * @brief Evaluate objective function at current state
   * 
   * Calls Python-defined objective function with current displacement,
   * stress, and strain fields. Used during optimization to compute
   * the objective value that drives the optimization process.
   * 
   * @param coords Gauss point coordinates
   * @param u_ptr Displacement field values
   * @param stress_ptr Stress tensor values  
   * @param strain_ptr Strain tensor values
   * @param o_ptr Output objective function values
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode evalInteriorObjectiveFunction(
      MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
      boost::shared_ptr<MatrixDouble> stress_ptr,
      boost::shared_ptr<MatrixDouble> strain_ptr,
      boost::shared_ptr<MatrixDouble> o_ptr, bool symmetrize = true);
 
  /**
   * @brief Compute gradient of objective function with respect to stress
   * 
   * Evaluates ∂f/∂σ where f is objective function and σ is stress tensor.
   * This gradient is used in the adjoint method to compute sensitivities
   * efficiently. Essential for gradient-based topology optimization.
   * 
   * @param coords Gauss point coordinates
   * @param u_ptr Displacement field values
   * @param stress_ptr Stress tensor values
   * @param strain_ptr Strain tensor values  
   * @param o_ptr Output gradient values
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode evalInteriorObjectiveGradientStress(
      MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
      boost::shared_ptr<MatrixDouble> stress_ptr,
      boost::shared_ptr<MatrixDouble> strain_ptr,
      boost::shared_ptr<MatrixDouble> o_ptr, bool symmetrize = true);
 
  /**
   * @brief Compute gradient of objective function with respect to strain
   * 
   * Evaluates ∂f/∂ε where f is objective function and ε is strain tensor.
   * Used when objective function depends directly on strain measures,
   * complementing stress-based gradients in adjoint sensitivity analysis.
   * 
   * @param coords Gauss point coordinates
   * @param u_ptr Displacement field values
   * @param stress_ptr Stress tensor values
   * @param strain_ptr Strain tensor values
   * @param o_ptr Output gradient values  
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode evalInteriorObjectiveGradientStrain(
      MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
      boost::shared_ptr<MatrixDouble> stress_ptr,
      boost::shared_ptr<MatrixDouble> strain_ptr,
      boost::shared_ptr<MatrixDouble> o_ptr, bool symmetrize = true);
 
  /**
   * @brief Compute gradient of objective function with respect to displacement
   * 
   * Evaluates ∂f/∂u where f is objective function and u is displacement field.
   * This provides direct sensitivity of objective to displacement changes,
   * used as right-hand side in adjoint equation K^T * λ = ∂f/∂u.
   * 
   * @param coords Gauss point coordinates  
   * @param u_ptr Displacement field values
   * @param stress_ptr Stress tensor values
   * @param strain_ptr Strain tensor values
   * @param o_ptr Output gradient values
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode evalInteriorObjectiveGradientU(
      MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
      boost::shared_ptr<MatrixDouble> stress_ptr,
      boost::shared_ptr<MatrixDouble> strain_ptr,
      boost::shared_ptr<MatrixDouble> o_ptr, bool symmetrize = true);
 
  MoFEMErrorCode
  evalBoundaryObjectiveFunction(MatrixDouble &coords,
                                boost::shared_ptr<MatrixDouble> u_ptr,
                                boost::shared_ptr<MatrixDouble> t_ptr,
                                boost::shared_ptr<VectorDouble> o_ptr,
                                bool symmetrize = true);
 
  MoFEMErrorCode
  evalBoundaryObjectiveGradientTraction(MatrixDouble &coords,
                                        boost::shared_ptr<MatrixDouble> u_ptr,
                                        boost::shared_ptr<MatrixDouble> t_ptr,
                                        boost::shared_ptr<MatrixDouble> o_ptr);
 
  MoFEMErrorCode
  evalBoundaryObjectiveGradientU(MatrixDouble &coords,
                                 boost::shared_ptr<MatrixDouble> u_ptr,
                                 boost::shared_ptr<MatrixDouble> t_ptr,
                                 boost::shared_ptr<MatrixDouble> o_ptr);
 
  /// Return number of topology optimization modes for given material block
  MoFEMErrorCode numberOfModes(int block_id, int &modes);
 
  /**
   * @brief Define spatial topology modes for design optimization
   * 
   * Generates basis functions that define how the geometry can be modified
   * during topology optimization. These modes serve as design variables
   * and define the design space for optimization.
   * 
   * @param block_id Material block identifier
   * @param coords Element coordinates
   * @param centroid Block centroid coordinates
   * @param bbodx Bounding box dimensions [xmin,xmax,ymin,ymax,zmin,zmax]
   * @param o_ptr Output mode vectors
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords,
                            std::array<double, 3> &centroid,
                            std::array<double, 6> &bbodx, MatrixDouble &o_ptr);
 
private:
  // Python interpreter objects for objective function evaluation
  bp::object mainNamespace; ///< Main Python namespace for script execution
 
  /**
   * @brief Internal implementation for objective function evaluation
   * 
   * Handles low-level Python function call with NumPy array conversion.
   * Converts MoFEM matrices to NumPy arrays, calls Python function,
   * and handles return value conversion.
   * 
   * @param coords NumPy array of coordinates
   * @param u NumPy array of displacements
   * @param stress NumPy array of stress tensors
   * @param strain NumPy array of strain tensors  
   * @param o Output NumPy array for objective values
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode interiorObjectiveFunctionImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  /**
   * @brief Internal implementation for stress gradient computation
   * 
   * Calls Python function to compute ∂f/∂σ with automatic array conversion.
   * Essential for adjoint-based sensitivity analysis in topology optimization.
   * 
   * @param coords NumPy array of coordinates
   * @param u NumPy array of displacements  
   * @param stress NumPy array of stress tensors
   * @param strain NumPy array of strain tensors
   * @param o Output NumPy array for gradient values
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode interiorObjectiveGradientStressImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  /**
   * @brief Internal implementation for strain gradient computation
   * 
   * Evaluates ∂f/∂ε through Python interface with NumPy array handling.
   * Provides strain-based sensitivities for comprehensive gradient computation.
   * 
   * @param coords NumPy array of coordinates
   * @param u NumPy array of displacements
   * @param stress NumPy array of stress tensors  
   * @param strain NumPy array of strain tensors
   * @param o Output NumPy array for gradient values
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode interiorObjectiveGradientStrainImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  /**
   * @brief Internal implementation for displacement gradient computation
   * 
   * Computes ∂f/∂u through Python interface for adjoint equation right-hand side.
   * This gradient drives the adjoint solution that enables efficient sensitivity
   * computation independent of design variable count.
   * 
   * @param coords NumPy array of coordinates
   * @param u NumPy array of displacements
   * @param stress NumPy array of stress tensors
   * @param strain NumPy array of strain tensors
   * @param o Output NumPy array for gradient values  
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode interiorObjectiveGradientUImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  MoFEMErrorCode boundaryObjectiveGradientTractionImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray t, np::ndarray &o
 
  );
 
  MoFEMErrorCode boundaryObjectiveFunctionImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray t, np::ndarray &o
 
  );
 
  MoFEMErrorCode boundaryObjectiveGradientUImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray t, np::ndarray &o
 
  );
 
  /**
   * @brief Internal implementation for topology mode generation
   * 
   * Calls Python function to generate spatial basis functions for topology
   * optimization. These modes define the design space and parametrize
   * allowable geometry modifications during optimization.
   * 
   * @param block_id Material block identifier
   * @param coords NumPy array of element coordinates
   * @param centroid NumPy array of block centroid
   * @param bbodx NumPy array of bounding box dimensions
   * @param o_ptr Output NumPy array for mode vectors
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode blockModesImpl(
 
      int block_id, np::ndarray coords, np::ndarray centroid, np::ndarray bbodx,
      np::ndarray &o_ptr
 
  );
 
  /**
   * @brief Convert std::vector to NumPy array for Python interface
   * 
   * Efficient conversion from MoFEM data structures to NumPy arrays
   * for seamless Python function calls without data copying.
   * 
   * @param data Source vector data
   * @param rows Number of rows in resulting array
   * @param nb_gauss_pts Number of Gauss points (affects array structure)
   * @return np::ndarray NumPy array for Python use
   */
  np::ndarray convertToNumPy(std::vector<double> &data, int rows,
                             int nb_gauss_pts);
 
  /**
   * @brief Convert raw pointer to NumPy array for Python interface
   * 
   * Low-level conversion for direct memory access to create NumPy arrays.
   * Provides zero-copy conversion when possible for performance.
   * 
   * @param ptr Raw data pointer
   * @param s Array size
   * @return np::ndarray NumPy array for Python use  
   */
  np::ndarray convertToNumPy(double *ptr, int s);
 
  /// Convert symmetric tensor storage to full matrix format
  MatrixDouble copyToFull(MatrixDouble &s);
 
  /// Convert full matrix to symmetric tensor storage format
  void copyToSymmetric(double *ptr, MatrixDouble &s);
};
 
/**
 * @brief Factory function to create Python-integrated objective function interface
 * 
 * Creates and initializes an ObjectiveFunctionDataImpl instance that bridges
 * MoFEM finite element computations with Python-defined objective functions.
 * This enables flexible objective function definition for topology optimization
 * while maintaining computational efficiency.
 * 
 * The Python file must define specific functions with correct signatures:
 * - objectiveInteriorFunction(coords, u, stress, strain) -> objective_value
 * - objectiveInteriorGradientStress(coords, u, stress, strain) -> gradient_array
 * - objectiveInteriorGradientStrain(coords, u, stress, strain) -> gradient_array  
 * - objectiveInteriorGradientU(coords, u, stress, strain) -> gradient_array
 * - numberOfModes(block_id) -> integer
 * - blockModes(block_id, coords, centroid, bbox) -> mode_array
 * 
 * @param py_file Path to Python file containing objective function definitions
 * @return boost::shared_ptr<ObjectiveFunctionData> Configured objective function interface
 * @throws MoFEM exception if Python initialization fails
 */
boost::shared_ptr<ObjectiveFunctionData>
create_python_objective_function(std::string py_file) {
  auto ptr = boost::make_shared<ObjectiveFunctionDataImpl>();
  CHK_THROW_MESSAGE(ptr->initPython(py_file), "init python");
  return ptr;
}
 
/**
 * @brief Initialize Python interpreter and load objective function script
 * 
 * This method sets up the Python environment for objective function evaluation
 * by loading a Python script that defines the required optimization functions.
 * It establishes the bridge between MoFEM's C++ finite element computations
 * and user-defined Python objective functions.
 * 
 * The Python script must define the following functions:
 * - f(coords, u, stress, strain): Main objective function
 * - f_stress(coords, u, stress, strain): Gradient w.r.t. stress ∂f/∂σ  
 * - f_strain(coords, u, stress, strain): Gradient w.r.t. strain ∂f/∂ε
 * - f_u(coords, u, stress, strain): Gradient w.r.t. displacement ∂f/∂u
 * - number_of_modes(block_id): Return number of topology modes
 * - block_modes(block_id, coords, centroid, bbox): Define topology modes
 * 
 * All functions receive NumPy arrays and must return NumPy arrays of
 * appropriate dimensions for the finite element computation.
 * 
 * @param py_file Path to Python script containing objective function definitions
 * @return MoFEMErrorCode Success or error code
 * @throws MOFEM_OPERATION_UNSUCCESSFUL if Python script has errors
 */
MoFEMErrorCode
ObjectiveFunctionDataImpl::initPython(const std::string py_file) {
  MoFEMFunctionBegin;
  try {
 
    // Create main Python module and namespace for script execution
    auto main_module = bp::import("__main__");
    mainNamespace = main_module.attr("__dict__");
 
    // Execute the Python script in the main namespace
    bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
 
    // Python callbacks are resolved lazily with explicit existence checks
 
  } catch (bp::error_already_set const &) {
    // Handle Python errors by printing to stderr and throwing MoFEM exception
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MatrixDouble ObjectiveFunctionDataImpl::copyToFull(MatrixDouble &s) {
  const auto nb_gauss_pts = s.size1();
  MatrixDouble f(nb_gauss_pts, 9);
  f.clear();
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  auto t_f =
      getFTensor2FromMat<3, 3, -1, DataLayoutTraits<DataLayout::GaussByCoeffs>>(
          f);
  auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
  for (size_t ii = 0; ii != nb_gauss_pts; ++ii) {
    t_f(i, j) = t_s(i, j);
    ++t_f;
    ++t_s;
  }
  return f;
};
 
void ObjectiveFunctionDataImpl::copyToSymmetric(double *ptr, MatrixDouble &s) {
  const auto nb_gauss_pts = s.size1();
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  auto t_f =
      getFTensor2FromPtr<3, 3, -1, DataLayoutTraits<DataLayout::GaussByCoeffs>>(
          ptr);
 
  auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
  for (size_t ii = 0; ii != nb_gauss_pts; ++ii) {
    t_s(i, j) = (t_f(i, j) || t_f(j, i)) / 2.0;
    ++t_f;
    ++t_s;
  }
}
 
/**
 * @brief Evaluate objective function at current finite element state
 * 
 * This method bridges MoFEM finite element data with Python-defined objective
 * functions for topology optimization. It handles the complete data conversion
 * workflow from MoFEM matrices to NumPy arrays, calls the Python objective
 * function, and converts results back to MoFEM format.
 * 
 * Process:
 * 1. Convert coordinate matrix to NumPy format for Python access
 * 2. Convert displacement field data to NumPy arrays
 * 3. Convert symmetric stress/strain tensors to full 3x3 matrix format
 * 4. Call Python objective function: f(coords, u, stress, strain)
 * 5. Extract results and copy back to MoFEM vector format
 * 
 * The objective function typically computes scalar quantities like:
 * - Compliance: ∫ u^T * f dΩ (minimize structural deformation)
 * - Stress constraints: ∫ ||σ - σ_target||² dΩ (control stress distribution)
 * - Volume constraints: ∫ ρ dΩ (material usage limitations)
 * 
 * @param coords Gauss point coordinates for current element
 * @param u_ptr Displacement field values at Gauss points
 * @param stress_ptr Cauchy stress tensor values (symmetric storage)
 * @param strain_ptr Strain tensor values (symmetric storage)  
 * @param o_ptr Output objective function values at each Gauss point
 * @return MoFEMErrorCode Success or error code
 */
MoFEMErrorCode ObjectiveFunctionDataImpl::evalInteriorObjectiveFunction(
    MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
    boost::shared_ptr<MatrixDouble> stress_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    boost::shared_ptr<MatrixDouble> o_ptr, bool symmetrize) {
  MoFEMFunctionBegin;
  try {
 
    // Convert coordinates to NumPy array for Python function
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
    // Convert displacement field to NumPy array
    auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
 
    // Convert symmetric tensor storage to full matrix format for Python
    // MoFEM stores symmetric tensors in like Voigt notation, Python expects full 3x3 matrices
    auto full_stress = symmetrize ? copyToFull(*(stress_ptr)) : *(stress_ptr);
    auto full_strain = symmetrize ? copyToFull(*(strain_ptr)) : *(strain_ptr);
 
    // Create NumPy arrays for stress and strain tensors
    auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
                                    full_stress.size2());
    auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
                                    full_strain.size2());
 
    // Prepare output array for objective function values
    const auto nb_gauss_pts = full_strain.size1();
    np::ndarray np_output =
        np::empty(bp::make_tuple(nb_gauss_pts), np::dtype::get_builtin<double>());
 
    // Call Python objective function implementation
    CHKERR interiorObjectiveFunctionImpl(np_coords, np_u, np_stress, np_strain,
                                 np_output);
 
    //Check the shape of returned array
    if (np_output.get_nd() != 1 || np_output.get_shape()[0] != nb_gauss_pts) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Wrong shape of Objective Function from python expected (" +
              std::to_string(nb_gauss_pts) + "), got (" +
              std::to_string(np_output.get_shape()[0]) + ")");
    }
 
    // Copy Python results back to a 1 x n matrix, matching common-data storage.
    o_ptr->resize(1, nb_gauss_pts, false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    std::copy(val_ptr, val_ptr + nb_gauss_pts, o_ptr->data().begin());
 
  } catch (bp::error_already_set const &) {
    // Handle Python errors with detailed error reporting
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Compute gradient of objective function with respect to stress tensor
 * 
 * This method evaluates ∂f/∂σ, the partial derivative of the objective function
 * with respect to the Cauchy stress tensor. This gradient is fundamental to the
 * adjoint method for topology optimization, as it provides the driving force
 * for the adjoint equation solution.
 * 
 * Mathematical context:
 * The adjoint method requires ∂f/∂σ to compute sensitivities efficiently.
 * For stress-based objectives like von Mises stress constraints:
 * ∂f/∂σ = ∂/∂σ[∫(σ_vm - σ_target)² dΩ] = 2(σ_vm - σ_target) * ∂σ_vm/∂σ
 * 
 * Process:
 * 1. Convert all field data to NumPy format for Python processing  
 * 2. Call Python function f_stress(coords, u, stress, strain)
 * 3. Python returns full 3x3 gradient matrices for each Gauss point
 * 4. Convert back to symmetric tensor storage used by MoFEM
 * 5. Store results for use in adjoint equation assembly
 * 
 * The resulting gradients drive the adjoint solution that enables efficient
 * computation of design sensitivities independent of design variable count.
 * 
 * @param coords Gauss point coordinates  
 * @param u_ptr Displacement field values
 * @param stress_ptr Current stress tensor values (symmetric storage)
 * @param strain_ptr Current strain tensor values (symmetric storage)
 * @param o_ptr Output stress gradients ∂f/∂σ (symmetric storage)
 * @return MoFEMErrorCode Success or error code
 */
MoFEMErrorCode ObjectiveFunctionDataImpl::evalInteriorObjectiveGradientStress(
    MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
    boost::shared_ptr<MatrixDouble> stress_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    boost::shared_ptr<MatrixDouble> o_ptr, bool symmetrize) {
  MoFEMFunctionBegin;
  try {
 
    // Convert coordinates and displacement field to NumPy format
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
    auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
 
    // Convert symmetric tensors to full 3x3 format for Python processing
    auto full_stress = symmetrize ? copyToFull(*(stress_ptr)) : *(stress_ptr);
    auto full_strain = symmetrize ? copyToFull(*(strain_ptr)) : *(strain_ptr);
 
    // Create NumPy arrays for stress and strain tensors
    auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
                                    full_stress.size2());
    auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
                                    full_strain.size2());
    // Prepare output array for stress gradients (full matrix format)
    np::ndarray np_output =
        np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
                  np::dtype::get_builtin<double>());
 
    // Call Python implementation for stress gradient computation
    CHKERR interiorObjectiveGradientStressImpl(np_coords, np_u, np_stress, np_strain,
                                       np_output);
 
    // Check the shape of returned array
    if (np_output.get_shape()[0] != full_strain.size1() ||
        np_output.get_shape()[1] != full_strain.size2()) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Wrong shape of Objective Gradient from python expected (" +
              std::to_string(full_strain.size1()) + ", " +
              std::to_string(full_strain.size2()) + "), got (" +
              std::to_string(np_output.get_shape()[0]) + ", " +
              std::to_string(np_output.get_shape()[1]) + ")");
    }
 
    // Prepare output matrix in the requested tensor storage format
    o_ptr->resize(stress_ptr->size1(), stress_ptr->size2(), false);
    if (symmetrize) {
      double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
      // Convert full matrix results back to symmetric tensor storage
      copyToSymmetric(val_ptr, *(o_ptr));
    } else {
      double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
      std::copy(val_ptr, val_ptr + stress_ptr->size1() * stress_ptr->size2(),
                o_ptr->data().begin());
    }
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Compute gradient of objective function with respect to strain tensor
 * 
 * This method evaluates ∂f/∂ε, the partial derivative of the objective function
 * with respect to the strain tensor. While many structural objectives depend
 * primarily on stress, strain-based gradients are important for certain
 * optimization formulations and provide additional sensitivity information.
 * 
 * Mathematical context:
 * For strain energy-based objectives: f = ½ε:C:ε
 * The gradient is: ∂f/∂ε = C:ε = σ (stress tensor)
 * 
 * For strain-based constraints or objectives like strain concentration:
 * ∂f/∂ε = ∂/∂ε[∫(ε_vm - ε_target)² dΩ] = 2(ε_vm - ε_target) * ∂ε_vm/∂ε
 * 
 * Process:
 * 1. Convert field data to NumPy format for Python compatibility
 * 2. Call Python function f_strain(coords, u, stress, strain)  
 * 3. Python returns gradient matrices ∂f/∂ε for each Gauss point
 * 4. Convert from full 3x3 format back to symmetric storage
 * 5. Results used in adjoint sensitivity analysis
 * 
 * This gradient complements stress-based gradients in comprehensive
 * sensitivity analysis for topology optimization problems.
 * 
 * @param coords Gauss point coordinates
 * @param u_ptr Displacement field values
 * @param stress_ptr Current stress tensor values (symmetric storage)
 * @param strain_ptr Current strain tensor values (symmetric storage)  
 * @param o_ptr Output strain gradients ∂f/∂ε (symmetric storage)
 * @return MoFEMErrorCode Success or error code
 */
MoFEMErrorCode ObjectiveFunctionDataImpl::evalInteriorObjectiveGradientStrain(
    MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
    boost::shared_ptr<MatrixDouble> stress_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    boost::shared_ptr<MatrixDouble> o_ptr, bool symmetrize) {
  MoFEMFunctionBegin;
  try {
 
    // Convert coordinates and displacement data to NumPy format
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
    auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
 
    // Convert symmetric tensor data to full 3x3 matrices for Python
    auto full_stress = symmetrize ? copyToFull(*(stress_ptr)) : *(stress_ptr);
    auto full_strain = symmetrize ? copyToFull(*(strain_ptr)) : *(strain_ptr);
 
    auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
                                    full_stress.size2());
    auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
                                    full_strain.size2());
 
    // Prepare output array for strain gradients
    np::ndarray np_output =
        np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
                  np::dtype::get_builtin<double>());
 
    // Call Python implementation for strain gradient computation
    CHKERR interiorObjectiveGradientStrainImpl(np_coords, np_u, np_stress, np_strain,
                                       np_output);
 
    // Check the shape of returned array
    if (np_output.get_shape()[0] != full_strain.size1() ||
        np_output.get_shape()[1] != full_strain.size2()) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Wrong shape of Objective Gradient from python expected (" +
              std::to_string(full_strain.size1()) + ", " +
              std::to_string(full_strain.size2()) + "), got (" +
              std::to_string(np_output.get_shape()[0]) + ", " +
              std::to_string(np_output.get_shape()[1]) + ")");
    }
 
    o_ptr->resize(strain_ptr->size1(), strain_ptr->size2(), false);
    if (symmetrize) {
      double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
      copyToSymmetric(val_ptr, *(o_ptr));
    } else {
      double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
      std::copy(val_ptr, val_ptr + strain_ptr->size1() * strain_ptr->size2(),
                o_ptr->data().begin());
    }
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Compute gradient of objective function with respect to displacement field
 * 
 * This method evaluates ∂f/∂u, the partial derivative of the objective function  
 * with respect to the displacement field. This gradient is crucial for the adjoint
 * method as it forms the right-hand side of the adjoint equation: K^T * λ = ∂f/∂u
 * 
 * Mathematical context:
 * The adjoint method solves: K^T * λ = ∂f/∂u
 * where λ are the adjoint variables (Lagrange multipliers)
 * 
 * For compliance minimization: f = ½u^T * K * u
 * The gradient is: ∂f/∂u = K * u (applied forces)
 * 
 * For displacement-based constraints: f = ||u - u_target||²
 * The gradient is: ∂f/∂u = 2(u - u_target)
 * 
 * Process:
 * 1. Convert all field data to NumPy format for Python processing
 * 2. Call Python function f_u(coords, u, stress, strain)
 * 3. Python returns displacement gradients for each component
 * 4. Copy results directly (no tensor conversion needed for vectors)
 * 5. Results drive adjoint equation solution for sensitivity analysis
 * 
 * This gradient is fundamental to adjoint-based topology optimization,
 * enabling efficient sensitivity computation for any number of design variables.
 * 
 * @param coords Gauss point coordinates
 * @param u_ptr Displacement field values
 * @param stress_ptr Current stress tensor values  
 * @param strain_ptr Current strain tensor values
 * @param o_ptr Output displacement gradients ∂f/∂u
 * @return MoFEMErrorCode Success or error code
 */
MoFEMErrorCode ObjectiveFunctionDataImpl::evalInteriorObjectiveGradientU(
    MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
    boost::shared_ptr<MatrixDouble> stress_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    boost::shared_ptr<MatrixDouble> o_ptr, bool symmetrize) {
  MoFEMFunctionBegin;
  try {
 
    // Convert coordinates and displacement field to NumPy format
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
    auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
 
    // Convert stress and strain tensors to full matrix format
    auto full_stress = symmetrize ? copyToFull(*(stress_ptr)) : *(stress_ptr);
    auto full_strain = symmetrize ? copyToFull(*(strain_ptr)) : *(strain_ptr);
 
    auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
                                    full_stress.size2());
    auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
                                    full_strain.size2());
 
    // Prepare output array for displacement gradients (same size as displacement field)
    np::ndarray np_output =
        np::empty(bp::make_tuple(u_ptr->size1(), u_ptr->size2()),
                  np::dtype::get_builtin<double>());
 
    // Call Python implementation for displacement gradient computation
    // Note: This should call interiorObjectiveGradientUImpl, not interiorObjectiveGradientStrainImpl
    CHKERR interiorObjectiveGradientUImpl(np_coords, np_u, np_stress, np_strain,
                                  np_output);
 
    // Check the shape of returned array
    if (np_output.get_shape()[0] != u_ptr->size1() ||
        np_output.get_shape()[1] != u_ptr->size2()) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Wrong shape of Objective Gradient from python expected (" +
              std::to_string(u_ptr->size1()) + ", " +
              std::to_string(u_ptr->size2()) + "), got (" +
              std::to_string(np_output.get_shape()[0]) + ", " +
              std::to_string(np_output.get_shape()[1]) + ")");
    }
 
    // Copy results directly to output matrix (no tensor conversion needed for vectors)
    o_ptr->resize(u_ptr->size1(), u_ptr->size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    std::copy(val_ptr, val_ptr + u_ptr->size1() * u_ptr->size2(),
              o_ptr->data().begin());
 
  } catch (bp::error_already_set const &) {
    // Handle Python errors with detailed reporting
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::evalBoundaryObjectiveGradientTraction(
    MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
    boost::shared_ptr<MatrixDouble> t_ptr,
    boost::shared_ptr<MatrixDouble> o_ptr) {
  MoFEMFunctionBegin;
  try {
 
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
    auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
    auto np_t = convertToNumPy(t_ptr->data(), t_ptr->size1(), t_ptr->size2());
 
    np::ndarray np_output =
        np::empty(bp::make_tuple(u_ptr->size1(), u_ptr->size2()),
                  np::dtype::get_builtin<double>());
 
    CHKERR boundaryObjectiveGradientTractionImpl(np_coords, np_u, np_t, np_output);
 
    // Check the shape of returned array
    if (np_output.get_shape()[0] != u_ptr->size1() ||
        np_output.get_shape()[1] != u_ptr->size2()) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Wrong shape of Objective Gradient from python expected (" +
              std::to_string(u_ptr->size1()) + ", " +
              std::to_string(u_ptr->size2()) + "), got (" +
              std::to_string(np_output.get_shape()[0]) + ", " +
              std::to_string(np_output.get_shape()[1]) + ")");
    }
 
    o_ptr->resize(u_ptr->size1(), u_ptr->size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    std::copy(val_ptr, val_ptr + u_ptr->size1() * u_ptr->size2(),
              o_ptr->data().begin());
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::evalBoundaryObjectiveFunction(
    MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
    boost::shared_ptr<MatrixDouble> t_ptr,
    boost::shared_ptr<VectorDouble> o_ptr, bool symmetrize) {
  MoFEMFunctionBegin;
  try {
    (void)symmetrize;
 
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
    auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
    auto np_t = convertToNumPy(t_ptr->data(), t_ptr->size1(), t_ptr->size2());
 
    np::ndarray np_output = np::empty(bp::make_tuple(t_ptr->size2()),
                                      np::dtype::get_builtin<double>());
 
    CHKERR boundaryObjectiveFunctionImpl(np_coords, np_u, np_t, np_output);
 
    // Check the shape of returned array
    if (np_output.get_nd() != 1 || np_output.get_shape()[0] != t_ptr->size2()) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Wrong shape of Objective Function from python expected (" +
              std::to_string(t_ptr->size2()) + "), got (" +
              std::to_string(np_output.get_shape()[0]) + ")");
    }
 
    o_ptr->resize(t_ptr->size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    std::copy(val_ptr, val_ptr + t_ptr->size2(), o_ptr->data().begin());
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::evalBoundaryObjectiveGradientU(
    MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
    boost::shared_ptr<MatrixDouble> t_ptr,
    boost::shared_ptr<MatrixDouble> o_ptr) {
  MoFEMFunctionBegin;
  try {
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
    auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
    auto np_t = convertToNumPy(t_ptr->data(), t_ptr->size1(), t_ptr->size2());
 
    np::ndarray np_output =
        np::empty(bp::make_tuple(u_ptr->size1(), u_ptr->size2()),
                  np::dtype::get_builtin<double>());
 
    CHKERR boundaryObjectiveGradientUImpl(np_coords, np_u, np_t, np_output);
 
    // Check the shape of returned array
    if (np_output.get_shape()[0] != u_ptr->size1() ||
        np_output.get_shape()[1] != u_ptr->size2()) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Wrong shape of Objective Gradient from python expected (" +
              std::to_string(u_ptr->size1()) + ", " +
              std::to_string(u_ptr->size2()) + "), got (" +
              std::to_string(np_output.get_shape()[0]) + ", " +
              std::to_string(np_output.get_shape()[1]) + ")");
    }
 
    o_ptr->resize(u_ptr->size1(), u_ptr->size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    std::copy(val_ptr, val_ptr + u_ptr->size1() * u_ptr->size2(),
              o_ptr->data().begin());
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Generate spatial topology modes for design optimization
 * 
 * This method defines the design parameterization for topology optimization by
 * generating spatial basis functions (modes) that describe how the geometry
 * can be modified during optimization. These modes serve as design variables
 * and define the feasible design space for the optimization problem.
 * 
 * Mathematical context:
 * The geometry modification is parameterized as: x_new = x_original + Σ(αᵢ * φᵢ(x))
 * where αᵢ are design variables and φᵢ(x) are spatial mode functions
 * 
 * Common mode types:
 * - Radial basis functions: φ(x) = exp(-||x-c||²/σ²) for localized changes
 * - Polynomial modes: φ(x) = xⁿyᵐzᵖ for global shape changes  
 * - Sinusoidal modes: φ(x) = sin(kx)cos(ly) for periodic patterns
 * - Principal component modes: Derived from geometric sensitivity analysis
 * 
 * Process:
 * 1. Query Python function for number of modes for this material block
 * 2. Convert coordinate data and geometric information to NumPy format
 * 3. Call Python function block_modes(block_id, coords, centroid, bbox)
 * 4. Python returns mode vectors for each coordinate at each mode
 * 5. Reshape and store modes for use as design variables in optimization
 * 
 * The modes enable efficient design space exploration and gradient-based
 * optimization while maintaining geometric feasibility and smoothness.
 * 
 * @param block_id Material block identifier for mode generation
 * @param coords Element coordinates where modes are evaluated
 * @param centroid Geometric centroid of the material block [x,y,z]
 * @param bbodx Bounding box dimensions [xmin,xmax,ymin,ymax,zmin,zmax]
 * @param o_ptr Output matrix: modes × (coordinates × spatial_dimension)
 * @return MoFEMErrorCode Success or error code
 */
MoFEMErrorCode ObjectiveFunctionDataImpl::blockModes(
    int block_id, MatrixDouble &coords, std::array<double, 3> &centroid,
    std::array<double, 6> &bbodx, MatrixDouble &o_ptr) {
  MoFEMFunctionBegin;
  try {
 
    // Query Python function for number of topology modes for this block
    int nb_modes;
    CHKERR numberOfModes(block_id, nb_modes);
 
    // Convert coordinate matrix to NumPy format for Python processing
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
 
    // Convert geometric information to NumPy arrays
    auto np_centroid =
        convertToNumPy(centroid.data(), 3); // Block centroid [x,y,z]
    auto np_bbodx = convertToNumPy(
        bbodx.data(), 6); // Bounding box [xmin,xmax,ymin,ymax,zmin,zmax]
 
    // Prepare output array: [modes × (coordinates * spatial_dimensions)]
    np::ndarray np_output =
        np::empty(bp::make_tuple(nb_modes, coords.size1(), coords.size2()),
                  np::dtype::get_builtin<double>());
 
    // Call Python implementation to generate topology modes
    CHKERR blockModesImpl(block_id, np_coords, np_centroid, np_bbodx,
                          np_output);
 
    // Check the shape of returned array
    if (np_output.get_shape()[0] != nb_modes ||
        np_output.get_shape()[1] != coords.size1() ||
        np_output.get_shape()[2] != coords.size2()) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Wrong shape of Modes from python expected (" +
                            std::to_string(nb_modes) + ", " +
                            std::to_string(coords.size1()) + ", " +
                            std::to_string(coords.size2()) + "), got (" +
                            std::to_string(np_output.get_shape()[0]) + ", " +
                            std::to_string(np_output.get_shape()[1]) + ", " +
                            std::to_string(np_output.get_shape()[2]) + ")");
    }
 
    // Reshape output matrix for MoFEM format: [modes × (coordinates * spatial_dimensions)]
    o_ptr.resize(nb_modes, coords.size1() * coords.size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    // Copy flattened mode data to output matrix
    std::copy(val_ptr, val_ptr + coords.size1() * coords.size2() * nb_modes,
              o_ptr.data().begin());
 
  } catch (bp::error_already_set const &) {
    // Handle Python errors in mode generation
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::interiorObjectiveFunctionImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  MoFEMFunctionBegin;
  try {
 
    if (bp::extract<bool>(mainNamespace.attr("__contains__")("f"))) {
      // Deprecated: Check for main objective function 'f' first for backward compatibility
      o = bp::extract<np::ndarray>(
          mainNamespace["f"](coords, u, stress, strain));
    } else if (bp::extract<bool>(
                   mainNamespace.attr("__contains__")("f_interior"))) {
      o = bp::extract<np::ndarray>(
          mainNamespace["f_interior"](coords, u, stress, strain));
    } else {
      CHK_THROW_MESSAGE(
          MOFEM_OPERATION_UNSUCCESSFUL,
          "Python function f_interior(coords,u,stress,strain) is not defined");
    }
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::interiorObjectiveGradientStressImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  MoFEMFunctionBegin;
  try {
 
    if (bp::extract<bool>(mainNamespace.attr("__contains__")("f_stress"))) {
      // Deprecated: keep legacy name first for backward compatibility
      o = bp::extract<np::ndarray>(
          mainNamespace["f_stress"](coords, u, stress, strain));
    } else if (bp::extract<bool>(
                   mainNamespace.attr("__contains__")("f_interior_stress"))) {
      o = bp::extract<np::ndarray>(
          mainNamespace["f_interior_stress"](coords, u, stress, strain));
    } else {
      CHK_THROW_MESSAGE(
          MOFEM_OPERATION_UNSUCCESSFUL,
          "Python function f_interior_stress(coords,u,stress,strain) is not defined");
    }
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::interiorObjectiveGradientStrainImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  MoFEMFunctionBegin;
  try {
 
    if (bp::extract<bool>(mainNamespace.attr("__contains__")("f_strain"))) {
      // Deprecated: keep legacy name first for backward compatibility
      o = bp::extract<np::ndarray>(
          mainNamespace["f_strain"](coords, u, stress, strain));
    } else if (bp::extract<bool>(
                   mainNamespace.attr("__contains__")("f_interior_strain"))) {
      o = bp::extract<np::ndarray>(
          mainNamespace["f_interior_strain"](coords, u, stress, strain));
    } else {
      CHK_THROW_MESSAGE(
          MOFEM_OPERATION_UNSUCCESSFUL,
          "Python function f_interior_strain(coords,u,stress,strain) is not defined");
    }
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::interiorObjectiveGradientUImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  MoFEMFunctionBegin;
  try {
 
    if (bp::extract<bool>(mainNamespace.attr("__contains__")("f_u"))) {
      // Deprecated: keep legacy name first for backward compatibility
      o = bp::extract<np::ndarray>(
          mainNamespace["f_u"](coords, u, stress, strain));
    } else if (bp::extract<bool>(
                   mainNamespace.attr("__contains__")("f_interior_u"))) {
      o = bp::extract<np::ndarray>(
          mainNamespace["f_interior_u"](coords, u, stress, strain));
    } else {
      CHK_THROW_MESSAGE(
          MOFEM_OPERATION_UNSUCCESSFUL,
          "Python function f_interior_u(coords,u,stress,strain) is not defined");
    }
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::boundaryObjectiveGradientTractionImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray t, np::ndarray &o
 
) {
  MoFEMFunctionBegin;
  try {
 
    if (bp::extract<bool>(mainNamespace.attr("__contains__")("f_boundary_t"))) {
      o = bp::extract<np::ndarray>(mainNamespace["f_boundary_t"](coords, u, t));
    } else {
      CHK_THROW_MESSAGE(
          MOFEM_OPERATION_UNSUCCESSFUL,
          "Python function f_boundary_t(coords,u,t) is not defined");
    }
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::boundaryObjectiveFunctionImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray t, np::ndarray &o
 
) {
  MoFEMFunctionBegin;
  try {
 
    if (bp::extract<bool>(mainNamespace.attr("__contains__")("f_boundary"))) {
      o = bp::extract<np::ndarray>(mainNamespace["f_boundary"](coords, u, t));
    } else if (bp::extract<bool>(
                   mainNamespace.attr("__contains__")("f_boundary_function"))) {
      o = bp::extract<np::ndarray>(
          mainNamespace["f_boundary_function"](coords, u, t));
    } else {
      CHK_THROW_MESSAGE(
          MOFEM_OPERATION_UNSUCCESSFUL,
          "Python function f_boundary(coords,u,t) is not defined");
    }
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::boundaryObjectiveGradientUImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray t, np::ndarray &o
 
) {
  MoFEMFunctionBegin;
  try {
 
    if (bp::extract<bool>(mainNamespace.attr("__contains__")("f_boundary_u"))) {
      o = bp::extract<np::ndarray>(mainNamespace["f_boundary_u"](coords, u, t));
    } else {
      CHK_THROW_MESSAGE(
          MOFEM_OPERATION_UNSUCCESSFUL,
          "Python function f_boundary_u(coords,u,t) is not defined");
    }
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::numberOfModes(int block_id,
                                                        int &modes) {
  MoFEMFunctionBegin;
  try {
 
    if (bp::extract<bool>(
            mainNamespace.attr("__contains__")("number_of_modes"))) {
      modes = bp::extract<int>(mainNamespace["number_of_modes"](block_id));
    } else {
      CHK_THROW_MESSAGE(
          MOFEM_OPERATION_UNSUCCESSFUL,
          "Python function number_of_modes(block_id) is not defined");
    }
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::blockModesImpl(int block_id,
                                                         np::ndarray coords,
                                                         np::ndarray centroid,
                                                         np::ndarray bbodx,
                                                         np::ndarray &o) {
  MoFEMFunctionBegin;
  try {
    if (bp::extract<bool>(mainNamespace.attr("__contains__")("block_modes"))) {
      o = bp::extract<np::ndarray>(
          mainNamespace["block_modes"](block_id, coords, centroid, bbodx));
    } else {
      CHK_THROW_MESSAGE(
          MOFEM_OPERATION_UNSUCCESSFUL,
          "Python function block_modes(block_id,coords,centroid,bbox) is not "
          "defined");
    }
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Converts a std::vector<double> to a NumPy ndarray.
 *
 * This function wraps the given vector data into a NumPy array with the
 * specified number of rows and Gauss points. The resulting ndarray shares
 * memory with the input vector, so changes to one will affect the other.
 *
 * @param data Reference to the vector containing double values to be converted.
 * @param rows Number of rows in the resulting NumPy array.
 * @param nb_gauss_pts Number of Gauss points (columns) in the resulting NumPy
 * array.
 * @return np::ndarray NumPy array view of the input data.
 *
 * @note
 * - `size` specifies the shape of the resulting ndarray as a tuple (rows,
 * nb_gauss_pts).
 * - `stride` specifies the step size in bytes to move to the next element in
 * memory. Here, it is set to sizeof(double), indicating contiguous storage for
 * each element.
 */
inline np::ndarray
ObjectiveFunctionDataImpl::convertToNumPy(std::vector<double> &data, int rows,
                                          int nb_gauss_pts) {
  auto dtype = np::dtype::get_builtin<double>();
  auto size = bp::make_tuple(rows, nb_gauss_pts);
  auto stride = bp::make_tuple(nb_gauss_pts * sizeof(double), sizeof(double));
  return (np::from_data(data.data(), dtype, size, stride, bp::object()));
}
 
inline np::ndarray ObjectiveFunctionDataImpl::convertToNumPy(double *ptr,
                                                             int s) {
  auto dtype = np::dtype::get_builtin<double>();
  auto size = bp::make_tuple(s);
  auto stride = bp::make_tuple(sizeof(double));
  return (np::from_data(ptr, dtype, size, stride, bp::object()));
}
 
#else
 
boost::shared_ptr<ObjectiveFunctionData>
create_python_objective_function(std::string) {
  CHK_THROW_MESSAGE(MOFEM_NOT_IMPLEMENTED,
                    "Python objective function requires Boost.NumPy");
  return {};
}
 
#endif // ENABLE_PYTHON_BINDING
 
} // namespace ShapeOptimization