ObjectiveFunctionDataImpl Struct Reference

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

Inheritance diagram for ObjectiveFunctionDataImpl:

Collaboration diagram for ObjectiveFunctionDataImpl:

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:

• 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:

• objectiveFunction(coords, u, stress, strain) -> objective_value
• objectiveGradientStress(coords, u, stress, strain) -> gradient_wrt_stress
• numberOfModes(block_id) -> number_of_topology_modes
• blockModes(block_id, coords, centroid, bbox) -> mode_vectors

Examples

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

Definition at line 2384 of file adjoint.cpp.

Constructor & Destructor Documentation

ObjectiveFunctionDataImpl()

ObjectiveFunctionDataImpl::ObjectiveFunctionDataImpl ( )

default

Examples

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

~ObjectiveFunctionDataImpl()

virtual ObjectiveFunctionDataImpl::~ObjectiveFunctionDataImpl ( )

virtualdefault

Examples

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

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_id	Material block identifier
coords	Element coordinates
centroid	Block centroid coordinates
bbodx	Bounding box dimensions [xmin,xmax,ymin,ymax,zmin,zmax]
o_ptr	Output 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_id	Material block identifier for mode generation
coords	Element coordinates where modes are evaluated
centroid	Geometric centroid of the material block [x,y,z]
bbodx	Bounding box dimensions [xmin,xmax,ymin,ymax,zmin,zmax]
o_ptr	Output 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.

                                                     {
  MoFEMFunctionBegin;
  try {
 
    // Query Python function for number of topology modes for this block
    int nb_modes = bp::extract<int>(objectNumberOfModes(block_id));
 
    // 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: [coordinates × modes × spatial_dimensions]
    np::ndarray np_output =
        np::empty(bp::make_tuple(coords.size1(), nb_modes, 3),
                  np::dtype::get_builtin<double>());
    
    // Call Python implementation to generate topology modes  
    CHKERR blockModesImpl(block_id, np_coords, np_centroid, np_bbodx,
                          np_output);
 
    // 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);
}

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_id	Material block identifier
coords	NumPy array of element coordinates
centroid	NumPy array of block centroid
bbodx	NumPy array of bounding box dimensions
o_ptr	Output 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.

                                                                       {
  MoFEMFunctionBegin;
  try {
    // call python function
    o = bp::extract<np::ndarray>(
        objectBlockModes(block_id, coords, centroid, bbodx));
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}

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

ptr	Raw data pointer
s	Array size

Returns

np::ndarray NumPy array for Python use

Definition at line 3282 of file adjoint.cpp.

                                                                    {
  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()));
}

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

data	Source vector data
rows	Number of rows in resulting array
nb_gauss_pts	Number 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

data	Reference to the vector containing double values to be converted.
rows	Number of rows in the resulting NumPy array.
nb_gauss_pts	Number 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.

                                                            {
  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()));
}

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.

                                                                  {
  MatrixDouble f(9, s.size2());
  f.clear();
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  auto t_f = getFTensor2FromMat<3, 3>(f);
  auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
  for (int ii = 0; ii != s.size2(); ++ii) {
    t_f(i, j) = t_s(i, j);
    ++t_f;
    ++t_s;
  }
  return f;
};

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.

                                                                            {
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  auto t_f = FTensor::Tensor2<FTensor::PackPtr<double *, 1>, 3, 3>{
 
      ptr + 0 * s.size2(), ptr + 1 * s.size2(),
      ptr + 2 * s.size2(), ptr + 3 * s.size2(),
      ptr + 4 * s.size2(), ptr + 5 * s.size2(),
      ptr + 6 * s.size2(), ptr + 7 * s.size2(),
      ptr + 8 * s.size2()
 
  };
  auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
  for (int ii = 0; ii != s.size2(); ++ii) {
    t_s(i, j) = (t_f(i, j) || t_f(j, i)) / 2.0;
    ++t_f;
    ++t_s;
  }
}

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

coords	Gauss point coordinates
u_ptr	Displacement field values
stress_ptr	Stress tensor values
strain_ptr	Strain tensor values
o_ptr	Output 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

coords	Gauss point coordinates for current element
u_ptr	Displacement field values at Gauss points
stress_ptr	Cauchy stress tensor values (symmetric storage)
strain_ptr	Strain tensor values (symmetric storage)
o_ptr	Output 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.

                                         {
  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 Voigt notation, Python expects full 3x3 matrices
    auto full_stress = copyToFull(*(stress_ptr));
    auto full_strain = copyToFull(*(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
    np::ndarray np_output = np::empty(bp::make_tuple(strain_ptr->size2()),
                                      np::dtype::get_builtin<double>());
    
    // Call Python objective function implementation
    CHKERR objectiveFunctionImpl(np_coords, np_u, np_stress, np_strain,
                                 np_output);
    
    // Copy Python results back to MoFEM vector format
    o_ptr->resize(stress_ptr->size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    std::copy(val_ptr, val_ptr + strain_ptr->size2(), 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);
}

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

coords	Gauss point coordinates
u_ptr	Displacement field values
stress_ptr	Stress tensor values
strain_ptr	Strain tensor values
o_ptr	Output 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

coords	Gauss point coordinates
u_ptr	Displacement field values
stress_ptr	Current stress tensor values (symmetric storage)
strain_ptr	Current strain tensor values (symmetric storage)
o_ptr	Output 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.

                                         {
  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 = copyToFull(*(stress_ptr));
    auto full_strain = copyToFull(*(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 objectiveGradientStrainImpl(np_coords, np_u, np_stress, np_strain,
                                       np_output);
    o_ptr->resize(strain_ptr->size1(), strain_ptr->size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    copyToSymmetric(val_ptr, *(o_ptr));
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}

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

coords	Gauss point coordinates
u_ptr	Displacement field values
stress_ptr	Stress tensor values
strain_ptr	Strain tensor values
o_ptr	Output 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

coords	Gauss point coordinates
u_ptr	Displacement field values
stress_ptr	Current stress tensor values (symmetric storage)
strain_ptr	Current strain tensor values (symmetric storage)
o_ptr	Output 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.

                                         {
  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 = copyToFull(*(stress_ptr));
    auto full_strain = copyToFull(*(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 objectiveGradientStressImpl(np_coords, np_u, np_stress, np_strain,
                                       np_output);
    
    // Prepare output matrix in symmetric format
    o_ptr->resize(stress_ptr->size1(), stress_ptr->size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    // Convert full matrix results back to symmetric tensor storage
    copyToSymmetric(val_ptr, *(o_ptr));
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}

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

coords	Gauss point coordinates
u_ptr	Displacement field values
stress_ptr	Stress tensor values
strain_ptr	Strain tensor values
o_ptr	Output 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

coords	Gauss point coordinates
u_ptr	Displacement field values
stress_ptr	Current stress tensor values
strain_ptr	Current strain tensor values
o_ptr	Output 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.

                                         {
  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 = copyToFull(*(stress_ptr));
    auto full_strain = copyToFull(*(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 objectiveGradientUImpl, not objectiveGradientStrainImpl
    CHKERR objectiveGradientUImpl(np_coords, np_u, np_stress, np_strain,
                                  np_output);
    
    // 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);
}

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_file

Path to Python script containing objective function definitions

Returns

MoFEMErrorCode Success or error code

Exceptions

MOFEM_OPERATION_UNSUCCESSFUL

if Python script has errors

Examples

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

Definition at line 2697 of file adjoint.cpp.

                                                             {
  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);
    
    // Create references to required Python functions for later calling
    objectiveFunction = mainNamespace["f"];                    // Main objective function
    objectiveGradientStress = mainNamespace["f_stress"];       // ∂f/∂σ gradient function  
    objectiveGradientStrain = mainNamespace["f_strain"];       // ∂f/∂ε gradient function
    objectiveGradientU = mainNamespace["f_u"];                 // ∂f/∂u gradient function
    objectNumberOfModes = mainNamespace["number_of_modes"];    // Topology mode count function
    objectBlockModes = mainNamespace["block_modes"];           // Topology mode definition function
 
  } 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);
}

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.

                                                                    {
  MoFEMFunctionBegin;
  try {
 
    modes = bp::extract<int>(objectNumberOfModes(block_id));
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}

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

coords	NumPy array of coordinates
u	NumPy array of displacements
stress	NumPy array of stress tensors
strain	NumPy array of strain tensors
o	Output 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.

  {
  MoFEMFunctionBegin;
  try {
 
    // call python function
    o = bp::extract<np::ndarray>(objectiveFunction(coords, u, stress, strain));
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}

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

coords	NumPy array of coordinates
u	NumPy array of displacements
stress	NumPy array of stress tensors
strain	NumPy array of strain tensors
o	Output 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.

  {
  MoFEMFunctionBegin;
  try {
 
    // call python function
    o = bp::extract<np::ndarray>(
        objectiveGradientStrain(coords, u, stress, strain));
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}

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

coords	NumPy array of coordinates
u	NumPy array of displacements
stress	NumPy array of stress tensors
strain	NumPy array of strain tensors
o	Output 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.

  {
  MoFEMFunctionBegin;
  try {
 
    // call python function
    o = bp::extract<np::ndarray>(
        objectiveGradientStress(coords, u, stress, strain));
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}

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

coords	NumPy array of coordinates
u	NumPy array of displacements
stress	NumPy array of stress tensors
strain	NumPy array of strain tensors
o	Output 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.

  {
  MoFEMFunctionBegin;
  try {
 
    // call python function
    o = bp::extract<np::ndarray>(objectiveGradientU(coords, u, stress, strain));
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}

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:

• tutorials/vec-7/adjoint.cpp