v0.15.0
Loading...
Searching...
No Matches
Public Member Functions | Private Member Functions | Private Attributes | List of all members
ObjectiveFunctionDataImpl Struct Reference

Implementation of ObjectiveFunctionData interface using Python integration. More...

Inheritance diagram for ObjectiveFunctionDataImpl:
[legend]
Collaboration diagram for ObjectiveFunctionDataImpl:
[legend]

Public Member Functions

 ObjectiveFunctionDataImpl ()=default
 
virtual ~ObjectiveFunctionDataImpl ()=default
 
MoFEMErrorCode initPython (const std::string py_file)
 Initialize Python interpreter and load objective function script.
 
MoFEMErrorCode evalObjectiveFunction (MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< VectorDouble > o_ptr)
 Evaluate objective function at current state.
 
MoFEMErrorCode evalObjectiveGradientStress (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)
 Compute gradient of objective function with respect to stress.
 
MoFEMErrorCode evalObjectiveGradientStrain (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)
 Compute gradient of objective function with respect to strain.
 
MoFEMErrorCode evalObjectiveGradientU (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)
 Compute gradient of objective function with respect to displacement.
 
MoFEMErrorCode numberOfModes (int block_id, int &modes)
 Return number of topology optimization modes for given material block.
 
MoFEMErrorCode blockModes (int block_id, MatrixDouble &coords, std::array< double, 3 > &centroid, std::array< double, 6 > &bbodx, MatrixDouble &o_ptr)
 Define spatial topology modes for design optimization.
 
- Public Member Functions inherited from ObjectiveFunctionData
virtual ~ObjectiveFunctionData ()=default
 

Private Member Functions

MoFEMErrorCode objectiveFunctionImpl (np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
 Internal implementation for objective function evaluation.
 
MoFEMErrorCode objectiveGradientStressImpl (np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
 Internal implementation for stress gradient computation.
 
MoFEMErrorCode objectiveGradientStrainImpl (np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
 Internal implementation for strain gradient computation.
 
MoFEMErrorCode objectiveGradientUImpl (np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
 Internal implementation for displacement gradient computation.
 
MoFEMErrorCode blockModesImpl (int block_id, np::ndarray coords, np::ndarray centroid, np::ndarray bbodx, np::ndarray &o_ptr)
 Internal implementation for topology mode generation.
 
np::ndarray convertToNumPy (std::vector< double > &data, int rows, int nb_gauss_pts)
 Convert std::vector to NumPy array for Python interface.
 
np::ndarray convertToNumPy (double *ptr, int s)
 Convert raw pointer to NumPy array for Python interface.
 
MatrixDouble copyToFull (MatrixDouble &s)
 Convert symmetric tensor storage to full matrix format.
 
void copyToSymmetric (double *ptr, MatrixDouble &s)
 Convert full matrix to symmetric tensor storage format

 

Private Attributes

bp::object mainNamespace
 Main Python namespace for script execution.
 
bp::object objectiveFunction
 Python function: f(coords,u,stress,strain) -> objective.
 
bp::object objectiveGradientStress
 Python function: ∂f/∂σ(coords,u,stress,strain) -> gradient

 
bp::object objectiveGradientStrain
 Python function: ∂f/∂ε(coords,u,stress,strain) -> gradient.
 
bp::object objectiveGradientU
 Python function: ∂f/∂u(coords,u,stress,strain) -> gradient.
 
bp::object objectNumberOfModes
 Python function: numberOfModes(block_id) -> int.
 
bp::object objectBlockModes
 Python function: blockModes(block_id,coords,centroid,bbox) -> modes.
 

Detailed Description

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:

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:

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2384 of file adjoint.cpp.

Constructor & Destructor Documentation

◆ ObjectiveFunctionDataImpl()

ObjectiveFunctionDataImpl::ObjectiveFunctionDataImpl ( )
default

◆ ~ObjectiveFunctionDataImpl()

virtual ObjectiveFunctionDataImpl::~ObjectiveFunctionDataImpl ( )
virtualdefault

Member Function Documentation

◆ blockModes()

MoFEMErrorCode ObjectiveFunctionDataImpl::blockModes ( int  block_id,
MatrixDouble coords,
std::array< double, 3 > &  centroid,
std::array< double, 6 > &  bbodx,
MatrixDouble o_ptr 
)
virtual

Define spatial topology modes for design optimization.

Generate 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.

Parameters
block_idMaterial block identifier
coordsElement coordinates
centroidBlock centroid coordinates
bbodxBounding box dimensions [xmin,xmax,ymin,ymax,zmin,zmax]
o_ptrOutput mode vectors
Returns
MoFEMErrorCode Success or error code

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.

Parameters
block_idMaterial block identifier for mode generation
coordsElement coordinates where modes are evaluated
centroidGeometric centroid of the material block [x,y,z]
bbodxBounding box dimensions [xmin,xmax,ymin,ymax,zmin,zmax]
o_ptrOutput matrix: modes × (coordinates × spatial_dimension)
Returns
MoFEMErrorCode Success or error code

Implements ObjectiveFunctionData.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 3093 of file adjoint.cpp.

3095 {
3097 try {
3098
3099 // Query Python function for number of topology modes for this block
3100 int nb_modes = bp::extract<int>(objectNumberOfModes(block_id));
3101
3102 // Convert coordinate matrix to NumPy format for Python processing
3103 auto np_coords =
3104 convertToNumPy(coords.data(), coords.size1(), coords.size2());
3105
3106 // Convert geometric information to NumPy arrays
3107 auto np_centroid = convertToNumPy(centroid.data(), 3); // Block centroid [x,y,z]
3108 auto np_bbodx = convertToNumPy(bbodx.data(), 6); // Bounding box [xmin,xmax,ymin,ymax,zmin,zmax]
3109
3110 // Prepare output array: [coordinates × modes × spatial_dimensions]
3111 np::ndarray np_output =
3112 np::empty(bp::make_tuple(coords.size1(), nb_modes, 3),
3113 np::dtype::get_builtin<double>());
3114
3115 // Call Python implementation to generate topology modes
3116 CHKERR blockModesImpl(block_id, np_coords, np_centroid, np_bbodx,
3117 np_output);
3118
3119 // Reshape output matrix for MoFEM format: [modes × (coordinates * spatial_dimensions)]
3120 o_ptr.resize(nb_modes, coords.size1() * coords.size2(), false);
3121 double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
3122 // Copy flattened mode data to output matrix
3123 std::copy(val_ptr, val_ptr + coords.size1() * coords.size2() * nb_modes,
3124 o_ptr.data().begin());
3125
3126 } catch (bp::error_already_set const &) {
3127 // Handle Python errors in mode generation
3128 PyErr_Print();
3130 }
3132}
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
MoFEMErrorCode blockModesImpl(int block_id, np::ndarray coords, np::ndarray centroid, np::ndarray bbodx, np::ndarray &o_ptr)
Internal implementation for topology mode generation.
Definition adjoint.cpp:3235
bp::object objectNumberOfModes
Python function: numberOfModes(block_id) -> int.
Definition adjoint.cpp:2503
np::ndarray convertToNumPy(std::vector< double > &data, int rows, int nb_gauss_pts)
Convert std::vector to NumPy array for Python interface.
Definition adjoint.cpp:3274

◆ blockModesImpl()

MoFEMErrorCode ObjectiveFunctionDataImpl::blockModesImpl ( int  block_id,
np::ndarray  coords,
np::ndarray  centroid,
np::ndarray  bbodx,
np::ndarray &  o_ptr 
)
private

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.

Parameters
block_idMaterial block identifier
coordsNumPy array of element coordinates
centroidNumPy array of block centroid
bbodxNumPy array of bounding box dimensions
o_ptrOutput NumPy array for mode vectors
Returns
MoFEMErrorCode Success or error code
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 3235 of file adjoint.cpp.

3239 {
3241 try {
3242 // call python function
3243 o = bp::extract<np::ndarray>(
3244 objectBlockModes(block_id, coords, centroid, bbodx));
3245 } catch (bp::error_already_set const &) {
3246 // print all other errors to stderr
3247 PyErr_Print();
3249 }
3251}
bp::object objectBlockModes
Python function: blockModes(block_id,coords,centroid,bbox) -> modes.
Definition adjoint.cpp:2504

◆ convertToNumPy() [1/2]

np::ndarray ObjectiveFunctionDataImpl::convertToNumPy ( double ptr,
int  s 
)
inlineprivate

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.

Parameters
ptrRaw data pointer
sArray size
Returns
np::ndarray NumPy array for Python use

Definition at line 3282 of file adjoint.cpp.

3283 {
3284 auto dtype = np::dtype::get_builtin<double>();
3285 auto size = bp::make_tuple(s);
3286 auto stride = bp::make_tuple(sizeof(double));
3287 return (np::from_data(ptr, dtype, size, stride, bp::object()));
3288}

◆ convertToNumPy() [2/2]

np::ndarray ObjectiveFunctionDataImpl::convertToNumPy ( std::vector< double > &  data,
int  rows,
int  nb_gauss_pts 
)
inlineprivate

Convert std::vector to NumPy array for Python interface.

Converts a std::vector<double> to a NumPy ndarray.

Efficient conversion from MoFEM data structures to NumPy arrays for seamless Python function calls without data copying.

Parameters
dataSource vector data
rowsNumber of rows in resulting array
nb_gauss_ptsNumber of Gauss points (affects array structure)
Returns
np::ndarray NumPy array for Python use

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.

Parameters
dataReference to the vector containing double values to be converted.
rowsNumber of rows in the resulting NumPy array.
nb_gauss_ptsNumber of Gauss points (columns) in the resulting NumPy array.
Returns
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.
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 3274 of file adjoint.cpp.

3275 {
3276 auto dtype = np::dtype::get_builtin<double>();
3277 auto size = bp::make_tuple(rows, nb_gauss_pts);
3278 auto stride = bp::make_tuple(nb_gauss_pts * sizeof(double), sizeof(double));
3279 return (np::from_data(data.data(), dtype, size, stride, bp::object()));
3280}

◆ copyToFull()

MatrixDouble ObjectiveFunctionDataImpl::copyToFull ( MatrixDouble s)
private

Convert symmetric tensor storage to full matrix format.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2724 of file adjoint.cpp.

2724 {
2725 MatrixDouble f(9, s.size2());
2726 f.clear();
2729 auto t_f = getFTensor2FromMat<3, 3>(f);
2730 auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
2731 for (int ii = 0; ii != s.size2(); ++ii) {
2732 t_f(i, j) = t_s(i, j);
2733 ++t_f;
2734 ++t_s;
2735 }
2736 return f;
2737};
#define FTENSOR_INDEX(DIM, I)
constexpr int SPACE_DIM
[Define dimension]
Definition adjoint.cpp:28
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77

◆ copyToSymmetric()

void ObjectiveFunctionDataImpl::copyToSymmetric ( double ptr,
MatrixDouble s 
)
private

Convert full matrix to symmetric tensor storage format

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2739 of file adjoint.cpp.

2739 {
2743
2744 ptr + 0 * s.size2(), ptr + 1 * s.size2(),
2745 ptr + 2 * s.size2(), ptr + 3 * s.size2(),
2746 ptr + 4 * s.size2(), ptr + 5 * s.size2(),
2747 ptr + 6 * s.size2(), ptr + 7 * s.size2(),
2748 ptr + 8 * s.size2()
2749
2750 };
2751 auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
2752 for (int ii = 0; ii != s.size2(); ++ii) {
2753 t_s(i, j) = (t_f(i, j) || t_f(j, i)) / 2.0;
2754 ++t_f;
2755 ++t_s;
2756 }
2757}

◆ evalObjectiveFunction()

MoFEMErrorCode ObjectiveFunctionDataImpl::evalObjectiveFunction ( MatrixDouble coords,
boost::shared_ptr< MatrixDouble u_ptr,
boost::shared_ptr< MatrixDouble stress_ptr,
boost::shared_ptr< MatrixDouble strain_ptr,
boost::shared_ptr< VectorDouble o_ptr 
)
virtual

Evaluate objective function at current state.

Evaluate objective function at current finite element 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.

Parameters
coordsGauss point coordinates
u_ptrDisplacement field values
stress_ptrStress tensor values
strain_ptrStrain tensor values
o_ptrOutput objective function values
Returns
MoFEMErrorCode Success or error code

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)
Parameters
coordsGauss point coordinates for current element
u_ptrDisplacement field values at Gauss points
stress_ptrCauchy stress tensor values (symmetric storage)
strain_ptrStrain tensor values (symmetric storage)
o_ptrOutput objective function values at each Gauss point
Returns
MoFEMErrorCode Success or error code

Implements ObjectiveFunctionData.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2786 of file adjoint.cpp.

2790 {
2792 try {
2793
2794 // Convert coordinates to NumPy array for Python function
2795 auto np_coords =
2796 convertToNumPy(coords.data(), coords.size1(), coords.size2());
2797 // Convert displacement field to NumPy array
2798 auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
2799
2800 // Convert symmetric tensor storage to full matrix format for Python
2801 // MoFEM stores symmetric tensors in Voigt notation, Python expects full 3x3 matrices
2802 auto full_stress = copyToFull(*(stress_ptr));
2803 auto full_strain = copyToFull(*(strain_ptr));
2804
2805 // Create NumPy arrays for stress and strain tensors
2806 auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
2807 full_stress.size2());
2808 auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
2809 full_strain.size2());
2810
2811 // Prepare output array for objective function values
2812 np::ndarray np_output = np::empty(bp::make_tuple(strain_ptr->size2()),
2813 np::dtype::get_builtin<double>());
2814
2815 // Call Python objective function implementation
2816 CHKERR objectiveFunctionImpl(np_coords, np_u, np_stress, np_strain,
2817 np_output);
2818
2819 // Copy Python results back to MoFEM vector format
2820 o_ptr->resize(stress_ptr->size2(), false);
2821 double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
2822 std::copy(val_ptr, val_ptr + strain_ptr->size2(), o_ptr->data().begin());
2823
2824 } catch (bp::error_already_set const &) {
2825 // Handle Python errors with detailed error reporting
2826 PyErr_Print();
2828 }
2830}
MatrixDouble copyToFull(MatrixDouble &s)
Convert symmetric tensor storage to full matrix format.
Definition adjoint.cpp:2724
MoFEMErrorCode objectiveFunctionImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for objective function evaluation.
Definition adjoint.cpp:3134

◆ evalObjectiveGradientStrain()

MoFEMErrorCode ObjectiveFunctionDataImpl::evalObjectiveGradientStrain ( 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 
)
virtual

Compute gradient of objective function with respect to strain.

Compute gradient of objective function with respect to strain tensor.

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.

Parameters
coordsGauss point coordinates
u_ptrDisplacement field values
stress_ptrStress tensor values
strain_ptrStrain tensor values
o_ptrOutput gradient values
Returns
MoFEMErrorCode Success or error code

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.

Parameters
coordsGauss point coordinates
u_ptrDisplacement field values
stress_ptrCurrent stress tensor values (symmetric storage)
strain_ptrCurrent strain tensor values (symmetric storage)
o_ptrOutput strain gradients ∂f/∂ε (symmetric storage)
Returns
MoFEMErrorCode Success or error code

Implements ObjectiveFunctionData.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2939 of file adjoint.cpp.

2943 {
2945 try {
2946
2947 // Convert coordinates and displacement data to NumPy format
2948 auto np_coords =
2949 convertToNumPy(coords.data(), coords.size1(), coords.size2());
2950 auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
2951
2952 // Convert symmetric tensor data to full 3x3 matrices for Python
2953 auto full_stress = copyToFull(*(stress_ptr));
2954 auto full_strain = copyToFull(*(strain_ptr));
2955 auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
2956 full_stress.size2());
2957 auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
2958 full_strain.size2());
2959
2960 // Prepare output array for strain gradients
2961 np::ndarray np_output =
2962 np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
2963 np::dtype::get_builtin<double>());
2964
2965 // Call Python implementation for strain gradient computation
2966 CHKERR objectiveGradientStrainImpl(np_coords, np_u, np_stress, np_strain,
2967 np_output);
2968 o_ptr->resize(strain_ptr->size1(), strain_ptr->size2(), false);
2969 double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
2970 copyToSymmetric(val_ptr, *(o_ptr));
2971
2972 } catch (bp::error_already_set const &) {
2973 PyErr_Print();
2975 }
2977}
MoFEMErrorCode objectiveGradientStrainImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for strain gradient computation.
Definition adjoint.cpp:3177
void copyToSymmetric(double *ptr, MatrixDouble &s)
Convert full matrix to symmetric tensor storage format
Definition adjoint.cpp:2739

◆ evalObjectiveGradientStress()

MoFEMErrorCode ObjectiveFunctionDataImpl::evalObjectiveGradientStress ( 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 
)
virtual

Compute gradient of objective function with respect to stress.

Compute gradient of objective function with respect to stress tensor.

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.

Parameters
coordsGauss point coordinates
u_ptrDisplacement field values
stress_ptrStress tensor values
strain_ptrStrain tensor values
o_ptrOutput gradient values
Returns
MoFEMErrorCode Success or error code

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.

Parameters
coordsGauss point coordinates
u_ptrDisplacement field values
stress_ptrCurrent stress tensor values (symmetric storage)
strain_ptrCurrent strain tensor values (symmetric storage)
o_ptrOutput stress gradients ∂f/∂σ (symmetric storage)
Returns
MoFEMErrorCode Success or error code

Implements ObjectiveFunctionData.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2862 of file adjoint.cpp.

2866 {
2868 try {
2869
2870 // Convert coordinates and displacement field to NumPy format
2871 auto np_coords =
2872 convertToNumPy(coords.data(), coords.size1(), coords.size2());
2873 auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
2874
2875 // Convert symmetric tensors to full 3x3 format for Python processing
2876 auto full_stress = copyToFull(*(stress_ptr));
2877 auto full_strain = copyToFull(*(strain_ptr));
2878
2879 // Create NumPy arrays for stress and strain tensors
2880 auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
2881 full_stress.size2());
2882 auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
2883 full_strain.size2());
2884
2885 // Prepare output array for stress gradients (full matrix format)
2886 np::ndarray np_output =
2887 np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
2888 np::dtype::get_builtin<double>());
2889
2890 // Call Python implementation for stress gradient computation
2891 CHKERR objectiveGradientStressImpl(np_coords, np_u, np_stress, np_strain,
2892 np_output);
2893
2894 // Prepare output matrix in symmetric format
2895 o_ptr->resize(stress_ptr->size1(), stress_ptr->size2(), false);
2896 double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
2897 // Convert full matrix results back to symmetric tensor storage
2898 copyToSymmetric(val_ptr, *(o_ptr));
2899
2900 } catch (bp::error_already_set const &) {
2901 PyErr_Print();
2903 }
2905}
MoFEMErrorCode objectiveGradientStressImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for stress gradient computation.
Definition adjoint.cpp:3155

◆ evalObjectiveGradientU()

MoFEMErrorCode ObjectiveFunctionDataImpl::evalObjectiveGradientU ( 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 
)

Compute gradient of objective function with respect to displacement.

Compute gradient of objective function with respect to displacement field.

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.

Parameters
coordsGauss point coordinates
u_ptrDisplacement field values
stress_ptrStress tensor values
strain_ptrStrain tensor values
o_ptrOutput gradient values
Returns
MoFEMErrorCode Success or error code

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.

Parameters
coordsGauss point coordinates
u_ptrDisplacement field values
stress_ptrCurrent stress tensor values
strain_ptrCurrent strain tensor values
o_ptrOutput displacement gradients ∂f/∂u
Returns
MoFEMErrorCode Success or error code
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 3013 of file adjoint.cpp.

3017 {
3019 try {
3020
3021 // Convert coordinates and displacement field to NumPy format
3022 auto np_coords =
3023 convertToNumPy(coords.data(), coords.size1(), coords.size2());
3024 auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
3025
3026 // Convert stress and strain tensors to full matrix format
3027 auto full_stress = copyToFull(*(stress_ptr));
3028 auto full_strain = copyToFull(*(strain_ptr));
3029 auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
3030 full_stress.size2());
3031 auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
3032 full_strain.size2());
3033
3034 // Prepare output array for displacement gradients (same size as displacement field)
3035 np::ndarray np_output =
3036 np::empty(bp::make_tuple(u_ptr->size1(), u_ptr->size2()),
3037 np::dtype::get_builtin<double>());
3038
3039 // Call Python implementation for displacement gradient computation
3040 // Note: This should call objectiveGradientUImpl, not objectiveGradientStrainImpl
3041 CHKERR objectiveGradientUImpl(np_coords, np_u, np_stress, np_strain,
3042 np_output);
3043
3044 // Copy results directly to output matrix (no tensor conversion needed for vectors)
3045 o_ptr->resize(u_ptr->size1(), u_ptr->size2(), false);
3046 double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
3047 std::copy(val_ptr, val_ptr + u_ptr->size1() * u_ptr->size2(),
3048 o_ptr->data().begin());
3049
3050 } catch (bp::error_already_set const &) {
3051 // Handle Python errors with detailed reporting
3052 PyErr_Print();
3054 }
3056}
MoFEMErrorCode objectiveGradientUImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for displacement gradient computation.
Definition adjoint.cpp:3199

◆ initPython()

MoFEMErrorCode ObjectiveFunctionDataImpl::initPython ( const std::string  py_file)

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.

Parameters
py_filePath to Python script containing objective function definitions
Returns
MoFEMErrorCode Success or error code
Exceptions
MOFEM_OPERATION_UNSUCCESSFULif Python script has errors
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2697 of file adjoint.cpp.

2697 {
2699 try {
2700
2701 // Create main Python module and namespace for script execution
2702 auto main_module = bp::import("__main__");
2703 mainNamespace = main_module.attr("__dict__");
2704
2705 // Execute the Python script in the main namespace
2706 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
2707
2708 // Create references to required Python functions for later calling
2709 objectiveFunction = mainNamespace["f"]; // Main objective function
2710 objectiveGradientStress = mainNamespace["f_stress"]; // ∂f/∂σ gradient function
2711 objectiveGradientStrain = mainNamespace["f_strain"]; // ∂f/∂ε gradient function
2712 objectiveGradientU = mainNamespace["f_u"]; // ∂f/∂u gradient function
2713 objectNumberOfModes = mainNamespace["number_of_modes"]; // Topology mode count function
2714 objectBlockModes = mainNamespace["block_modes"]; // Topology mode definition function
2715
2716 } catch (bp::error_already_set const &) {
2717 // Handle Python errors by printing to stderr and throwing MoFEM exception
2718 PyErr_Print();
2720 }
2722}
bp::object objectiveGradientStress
Python function: ∂f/∂σ(coords,u,stress,strain) -> gradient
Definition adjoint.cpp:2500
bp::object objectiveGradientStrain
Python function: ∂f/∂ε(coords,u,stress,strain) -> gradient.
Definition adjoint.cpp:2501
bp::object objectiveFunction
Python function: f(coords,u,stress,strain) -> objective.
Definition adjoint.cpp:2499
bp::object objectiveGradientU
Python function: ∂f/∂u(coords,u,stress,strain) -> gradient.
Definition adjoint.cpp:2502
bp::object mainNamespace
Main Python namespace for script execution.
Definition adjoint.cpp:2498

◆ numberOfModes()

MoFEMErrorCode ObjectiveFunctionDataImpl::numberOfModes ( int  block_id,
int &  modes 
)
virtual

Return number of topology optimization modes for given material block.

Implements ObjectiveFunctionData.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 3220 of file adjoint.cpp.

3221 {
3223 try {
3224
3225 modes = bp::extract<int>(objectNumberOfModes(block_id));
3226
3227 } catch (bp::error_already_set const &) {
3228 // print all other errors to stderr
3229 PyErr_Print();
3231 }
3233}

◆ objectiveFunctionImpl()

MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveFunctionImpl ( np::ndarray  coords,
np::ndarray  u,
np::ndarray  stress,
np::ndarray  strain,
np::ndarray &  o 
)
private

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.

Parameters
coordsNumPy array of coordinates
uNumPy array of displacements
stressNumPy array of stress tensors
strainNumPy array of strain tensors
oOutput NumPy array for objective values
Returns
MoFEMErrorCode Success or error code
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 3134 of file adjoint.cpp.

3140 {
3142 try {
3143
3144 // call python function
3145 o = bp::extract<np::ndarray>(objectiveFunction(coords, u, stress, strain));
3146
3147 } catch (bp::error_already_set const &) {
3148 // print all other errors to stderr
3149 PyErr_Print();
3151 }
3153}

◆ objectiveGradientStrainImpl()

MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveGradientStrainImpl ( np::ndarray  coords,
np::ndarray  u,
np::ndarray  stress,
np::ndarray  strain,
np::ndarray &  o 
)
private

Internal implementation for strain gradient computation.

Evaluates ∂f/∂ε through Python interface with NumPy array handling. Provides strain-based sensitivities for comprehensive gradient computation.

Parameters
coordsNumPy array of coordinates
uNumPy array of displacements
stressNumPy array of stress tensors
strainNumPy array of strain tensors
oOutput NumPy array for gradient values
Returns
MoFEMErrorCode Success or error code
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 3177 of file adjoint.cpp.

3183 {
3185 try {
3186
3187 // call python function
3188 o = bp::extract<np::ndarray>(
3189 objectiveGradientStrain(coords, u, stress, strain));
3190
3191 } catch (bp::error_already_set const &) {
3192 // print all other errors to stderr
3193 PyErr_Print();
3195 }
3197}

◆ objectiveGradientStressImpl()

MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveGradientStressImpl ( np::ndarray  coords,
np::ndarray  u,
np::ndarray  stress,
np::ndarray  strain,
np::ndarray &  o 
)
private

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.

Parameters
coordsNumPy array of coordinates
uNumPy array of displacements
stressNumPy array of stress tensors
strainNumPy array of strain tensors
oOutput NumPy array for gradient values
Returns
MoFEMErrorCode Success or error code
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 3155 of file adjoint.cpp.

3161 {
3163 try {
3164
3165 // call python function
3166 o = bp::extract<np::ndarray>(
3167 objectiveGradientStress(coords, u, stress, strain));
3168
3169 } catch (bp::error_already_set const &) {
3170 // print all other errors to stderr
3171 PyErr_Print();
3173 }
3175}

◆ objectiveGradientUImpl()

MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveGradientUImpl ( np::ndarray  coords,
np::ndarray  u,
np::ndarray  stress,
np::ndarray  strain,
np::ndarray &  o 
)
private

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.

Parameters
coordsNumPy array of coordinates
uNumPy array of displacements
stressNumPy array of stress tensors
strainNumPy array of strain tensors
oOutput NumPy array for gradient values
Returns
MoFEMErrorCode Success or error code
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 3199 of file adjoint.cpp.

3205 {
3207 try {
3208
3209 // call python function
3210 o = bp::extract<np::ndarray>(objectiveGradientU(coords, u, stress, strain));
3211
3212 } catch (bp::error_already_set const &) {
3213 // print all other errors to stderr
3214 PyErr_Print();
3216 }
3218}

Member Data Documentation

◆ mainNamespace

bp::object ObjectiveFunctionDataImpl::mainNamespace
private

Main Python namespace for script execution.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2498 of file adjoint.cpp.

◆ objectBlockModes

bp::object ObjectiveFunctionDataImpl::objectBlockModes
private

Python function: blockModes(block_id,coords,centroid,bbox) -> modes.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2504 of file adjoint.cpp.

◆ objectiveFunction

bp::object ObjectiveFunctionDataImpl::objectiveFunction
private

Python function: f(coords,u,stress,strain) -> objective.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2499 of file adjoint.cpp.

◆ objectiveGradientStrain

bp::object ObjectiveFunctionDataImpl::objectiveGradientStrain
private

Python function: ∂f/∂ε(coords,u,stress,strain) -> gradient.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2501 of file adjoint.cpp.

◆ objectiveGradientStress

bp::object ObjectiveFunctionDataImpl::objectiveGradientStress
private

Python function: ∂f/∂σ(coords,u,stress,strain) -> gradient

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2500 of file adjoint.cpp.

◆ objectiveGradientU

bp::object ObjectiveFunctionDataImpl::objectiveGradientU
private

Python function: ∂f/∂u(coords,u,stress,strain) -> gradient.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2502 of file adjoint.cpp.

◆ objectNumberOfModes

bp::object ObjectiveFunctionDataImpl::objectNumberOfModes
private

Python function: numberOfModes(block_id) -> int.

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2503 of file adjoint.cpp.


The documentation for this struct was generated from the following file: