v0.15.0
Loading...
Searching...
No Matches
Classes | Typedefs | Functions | Variables
adjoint.cpp File Reference

Topology optimization using adjoint method with Python objective functions. More...

#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
#include <MoFEM.hpp>
#include <ElasticSpring.hpp>
#include <FluidLevel.hpp>
#include <CalculateTraction.hpp>
#include <NaturalDomainBC.hpp>
#include <NaturalBoundaryBC.hpp>
#include <HookeOps.hpp>

Go to the source code of this file.

Classes

struct  DomainBCs
 [Define entities] More...
 
struct  BoundaryBCs
 Boundary conditions marker. More...
 
struct  PostProcEleByDim< 2 >
 
struct  PostProcEleByDim< 3 >
 
struct  ObjectiveFunctionData
 Abstract interface for Python-defined objective functions. More...
 
struct  Example
 [Example] More...
 
struct  OpAdJointTestOp< SPACE_DIM, IntegrationType::GAUSS, DomainBaseOp >
 Adjoint test operator for elastic problems. More...
 
struct  OpAdJointGradTimesSymTensor< SPACE_DIM, IntegrationType::GAUSS, DomainBaseOp >
 
struct  OpAdJointObjective
 Operator for computing objective function and gradients in topology optimization. More...
 
struct  ObjectiveFunctionDataImpl
 Implementation of ObjectiveFunctionData interface using Python integration. More...
 

Typedefs

using EntData = EntitiesFieldData::EntData
 Entity data for field operations.
 
using DomainEle = PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomainEle
 Domain finite elements.
 
using BoundaryEle = PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::BoundaryEle
 Boundary finite elements.
 
using DomainEleOp = DomainEle::UserDataOperator
 Domain element operators.
 
using BoundaryEleOp = BoundaryEle::UserDataOperator
 
using DomainRhsBCs = NaturalBC< DomainEleOp >::Assembly< A >::LinearForm< I >
 Domain RHS natural BCs.
 
using OpDomainRhsBCs = DomainRhsBCs::OpFlux< DomainBCs, 1, SPACE_DIM >
 Domain flux operator.
 
using BoundaryRhsBCs = NaturalBC< BoundaryEleOp >::Assembly< A >::LinearForm< I >
 Boundary RHS natural BCs

 
using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux< BoundaryBCs, 1, SPACE_DIM >
 Boundary flux operator.
 
using BoundaryLhsBCs = NaturalBC< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >
 Boundary LHS natural BCs.
 
using OpBoundaryLhsBCs = BoundaryLhsBCs::OpFlux< BoundaryBCs, 1, SPACE_DIM >
 Boundary LHS flux operator.
 
using PostProcEleDomain = PostProcEleByDim< SPACE_DIM >::PostProcEleDomain
 
using SideEle = PostProcEleByDim< SPACE_DIM >::SideEle
 
using PostProcEleBdy = PostProcEleByDim< SPACE_DIM >::PostProcEleBdy
 
using DomainBaseOp = FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase
 [Postprocess results]
 

Functions

boost::shared_ptr< ObjectiveFunctionDatacreate_python_objective_function (std::string py_file)
 Factory function to create Python-integrated objective function interface.
 
Range get_range_from_block (MoFEM::Interface &m_field, const std::string block_name, int dim)
 
MoFEMErrorCode save_range (moab::Interface &moab, const std::string name, const Range r)
 
template<int DIM>
auto diff_symmetrize (FTensor::Number< DIM >)
 
int main (int argc, char *argv[])
 Main function for topology optimization tutorial using adjoint method.
 

Variables

constexpr int BASE_DIM = 1
 [Constants and material properties]
 
constexpr int SPACE_DIM
 [Define dimension]
 
constexpr AssemblyType A = AssemblyType::PETSC
 [Define dimension]
 
constexpr IntegrationType I
 Use Gauss quadrature for integration.
 
constexpr double young_modulus = 1
 [Material properties for linear elasticity]
 
constexpr double poisson_ratio = 0.3
 Poisson's ratio ν
 
constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio))
 Bulk modulus K = E/(3(1-2ν))
 
constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio))
 Shear modulus G = E/(2(1+ν))
 
PetscBool is_plane_strain = PETSC_FALSE
 
static char help [] = "...\n\n"
 [calculateGradient]
 

Detailed Description

Topology optimization using adjoint method with Python objective functions.


This tutorial demonstrates:

Definition in file adjoint.cpp.

Typedef Documentation

◆ BoundaryEle

Boundary finite elements.

Definition at line 48 of file adjoint.cpp.

◆ BoundaryEleOp

Boundary element operators

Definition at line 50 of file adjoint.cpp.

◆ BoundaryLhsBCs

Boundary LHS natural BCs.

Definition at line 62 of file adjoint.cpp.

◆ BoundaryRhsBCs

Boundary RHS natural BCs

Definition at line 60 of file adjoint.cpp.

◆ DomainBaseOp

[Postprocess results]

[calculateGradient]

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

Definition at line 1544 of file adjoint.cpp.

◆ DomainEle

Domain finite elements.

Definition at line 47 of file adjoint.cpp.

◆ DomainEleOp

Domain element operators.

Definition at line 49 of file adjoint.cpp.

◆ DomainRhsBCs

Domain RHS natural BCs.

Definition at line 58 of file adjoint.cpp.

◆ EntData

Entity data for field operations.

Definition at line 46 of file adjoint.cpp.

◆ OpBoundaryLhsBCs

Boundary LHS flux operator.

Definition at line 63 of file adjoint.cpp.

◆ OpBoundaryRhsBCs

Boundary flux operator.

Definition at line 61 of file adjoint.cpp.

◆ OpDomainRhsBCs

Domain flux operator.

Definition at line 59 of file adjoint.cpp.

◆ PostProcEleBdy

Definition at line 81 of file adjoint.cpp.

◆ PostProcEleDomain

Definition at line 79 of file adjoint.cpp.

◆ SideEle

Definition at line 80 of file adjoint.cpp.

Function Documentation

◆ create_python_objective_function()

boost::shared_ptr< ObjectiveFunctionData > create_python_objective_function ( std::string  py_file)

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:

  • objectiveFunction(coords, u, stress, strain) -> objective_value
  • objectiveGradientStress(coords, u, stress, strain) -> gradient_array
  • objectiveGradientStrain(coords, u, stress, strain) -> gradient_array
  • objectiveGradientU(coords, u, stress, strain) -> gradient_array
  • numberOfModes(block_id) -> integer
  • blockModes(block_id, coords, centroid, bbox) -> mode_array
Parameters
py_filePath to Python file containing objective function definitions
Returns
boost::shared_ptr<ObjectiveFunctionData> Configured objective function interface
Exceptions
MoFEMexception if Python initialization fails
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 2667 of file adjoint.cpp.

2667 {
2668 auto ptr = boost::make_shared<ObjectiveFunctionDataImpl>();
2669 CHK_THROW_MESSAGE(ptr->initPython(py_file), "init python");
2670 return ptr;
2671}
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.

◆ diff_symmetrize()

template<int DIM>
auto diff_symmetrize ( FTensor::Number< DIM >  )
inline
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 1656 of file adjoint.cpp.

1656 {
1657
1658 FTensor::Index<'i', DIM> i;
1659 FTensor::Index<'j', DIM> j;
1660 FTensor::Index<'k', DIM> k;
1661 FTensor::Index<'l', DIM> l;
1662
1664
1665 t_diff(i, j, k, l) = 0;
1666 t_diff(0, 0, 0, 0) = 1;
1667 t_diff(1, 1, 1, 1) = 1;
1668
1669 t_diff(1, 0, 1, 0) = 0.5;
1670 t_diff(1, 0, 0, 1) = 0.5;
1671
1672 t_diff(0, 1, 0, 1) = 0.5;
1673 t_diff(0, 1, 1, 0) = 0.5;
1674
1675 if constexpr (DIM == 3) {
1676 t_diff(2, 2, 2, 2) = 1;
1677
1678 t_diff(2, 0, 2, 0) = 0.5;
1679 t_diff(2, 0, 0, 2) = 0.5;
1680 t_diff(0, 2, 0, 2) = 0.5;
1681 t_diff(0, 2, 2, 0) = 0.5;
1682
1683 t_diff(2, 1, 2, 1) = 0.5;
1684 t_diff(2, 1, 1, 2) = 0.5;
1685 t_diff(1, 2, 1, 2) = 0.5;
1686 t_diff(1, 2, 2, 1) = 0.5;
1687 }
1688
1689 return t_diff;
1690};
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k

◆ get_range_from_block()

Range get_range_from_block ( MoFEM::Interface m_field,
const std::string  block_name,
int  dim 
)
Examples
mofem/tutorials/vec-7/adjoint.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 3290 of file adjoint.cpp.

3291 {
3292 Range r;
3293
3294 auto mesh_mng = m_field.getInterface<MeshsetsManager>();
3295 auto bcs = mesh_mng->getCubitMeshsetPtr(
3296
3297 std::regex((boost::format("%s(.*)") % block_name).str())
3298
3299 );
3300
3301 for (auto bc : bcs) {
3302 Range faces;
3303 CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
3304 faces, true),
3305 "get meshset ents");
3306 r.merge(faces);
3307 }
3308
3309 for (auto dd = dim - 1; dd >= 0; --dd) {
3310 if (dd >= 0) {
3311 Range ents;
3312 CHK_MOAB_THROW(m_field.get_moab().get_adjacencies(r, dd, false, ents,
3313 moab::Interface::UNION),
3314 "get adjs");
3315 r.merge(ents);
3316 } else {
3317 Range verts;
3318 CHK_MOAB_THROW(m_field.get_moab().get_connectivity(r, verts),
3319 "get verts");
3320 r.merge(verts);
3321 }
3323 m_field.getInterface<CommInterface>()->synchroniseEntities(r), "comm");
3324 }
3325
3326 return r;
3327};
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition ddTensor0.hpp:33
int r
Definition sdf.py:8
Managing BitRefLevels.
MoFEMErrorCode synchroniseEntities(Range &ent, std::map< int, Range > *received_ents, int verb=DEFAULT_VERBOSITY)
synchronize entity range on processors (collective)
virtual moab::Interface & get_moab()=0
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.

◆ main()

int main ( int  argc,
char *  argv[] 
)

Main function for topology optimization tutorial using adjoint method.

This tutorial demonstrates structural topology optimization using:

  • MoFEM finite element library for elasticity analysis
  • Adjoint method for efficient gradient computation
  • TAO optimization library for design optimization
  • Python interface for flexible objective function definition

Workflow:

  1. Initialize MoFEM, PETSc and Python environments
  2. Read mesh and setup finite element problem
  3. Define objective function via Python interface
  4. Compute topology optimization modes as design variables
  5. Run gradient-based optimization using adjoint sensitivities
  6. Post-process optimized design

The adjoint method enables efficient gradient computation with cost independent of the number of design variables, making it suitable for large-scale topology optimization problems.

Required input files:

  • Mesh file (.h5m format from CUBIT)
  • Parameter file (param_file.petsc)
  • Objective function (objective_function.py)
Parameters
argcCommand line argument count
argvCommand line argument values
Returns
int Exit code (0 for success)

[Register MoFEM discrete manager in PETSc]

[Register MoFEM discrete manager in PETSc

[Create MoAB]

< mesh database

< mesh database interface

[Create MoAB]

[Create MoFEM]

< finite element database

< finite element database interface

[Create MoFEM]

[Example]

[Example]

Definition at line 2304 of file adjoint.cpp.

2304 {
2305
2306 // Initialize Python environment for objective function interface
2307 Py_Initialize();
2308 np::initialize();
2309
2310 // Initialize MoFEM/PETSc and MOAB data structures
2311 const char param_file[] = "param_file.petsc";
2312 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
2313
2314 auto core_log = logging::core::get();
2315 core_log->add_sink(
2317
2318 core_log->add_sink(
2319 LogManager::createSink(LogManager::getStrmSync(), "FieldEvaluator"));
2320 LogManager::setLog("FieldEvaluator");
2321 MOFEM_LOG_TAG("FieldEvaluator", "field_eval");
2322
2323 try {
2324
2325 //! [Register MoFEM discrete manager in PETSc]
2326 DMType dm_name = "DMMOFEM";
2327 CHKERR DMRegister_MoFEM(dm_name);
2328 DMType dm_name_mg = "DMMOFEM_MG";
2330 //! [Register MoFEM discrete manager in PETSc
2331
2332 //! [Create MoAB]
2333 moab::Core mb_instance; ///< mesh database
2334 moab::Interface &moab = mb_instance; ///< mesh database interface
2335 //! [Create MoAB]
2336
2337 //! [Create MoFEM]
2338 MoFEM::Core core(moab); ///< finite element database
2339 MoFEM::Interface &m_field = core; ///< finite element database interface
2340 //! [Create MoFEM]
2341
2342 //! [Example]
2343 Example ex(m_field);
2344 CHKERR ex.runProblem();
2345 //! [Example]
2346 }
2348
2350
2351 if (Py_FinalizeEx() < 0) {
2352 exit(120);
2353 }
2354}
static char help[]
[calculateGradient]
Definition adjoint.cpp:2272
#define CATCH_ERRORS
Catch errors.
#define CHKERR
Inline error check.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
[Example]
Definition plastic.cpp:216
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.

◆ save_range()

MoFEMErrorCode save_range ( moab::Interface &  moab,
const std::string  name,
const Range  r 
)

Definition at line 3329 of file adjoint.cpp.

3330 {
3332 auto out_meshset = get_temp_meshset_ptr(moab);
3333 CHKERR moab.add_entities(*out_meshset, r);
3334 CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
3336}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.

Variable Documentation

◆ A

constexpr AssemblyType A = AssemblyType::PETSC
constexpr

[Define dimension]

Use PETSc for matrix/vector assembly

Examples
mofem/tutorials/vec-1/eigen_elastic.cpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 32 of file adjoint.cpp.

◆ BASE_DIM

constexpr int BASE_DIM = 1
constexpr

[Constants and material properties]

Dimension of the base functions

Definition at line 25 of file adjoint.cpp.

◆ bulk_modulus_K

constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio))
constexpr

Bulk modulus K = E/(3(1-2ν))

Examples
mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, and mofem/tutorials/vec-1/eigen_elastic.cpp.

Definition at line 39 of file adjoint.cpp.

◆ help

char help[] = "...\n\n"
static

[calculateGradient]

Definition at line 2272 of file adjoint.cpp.

◆ I

constexpr IntegrationType I
constexpr
Initial value:
=
IntegrationType::GAUSS

Use Gauss quadrature for integration.

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

Definition at line 33 of file adjoint.cpp.

◆ is_plane_strain

PetscBool is_plane_strain = PETSC_FALSE

Flag for plane strain vs plane stress in 2D

Examples
mofem/tutorials/vec-2/src/HenckyOps.hpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 42 of file adjoint.cpp.

◆ poisson_ratio

constexpr double poisson_ratio = 0.3
constexpr

Poisson's ratio ν

Examples
mofem/tutorials/adv-0/plastic.cpp, and mofem/tutorials/vec-1/eigen_elastic.cpp.

Definition at line 38 of file adjoint.cpp.

◆ shear_modulus_G

constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio))
constexpr

Shear modulus G = E/(2(1+ν))

Examples
mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, and mofem/tutorials/vec-1/eigen_elastic.cpp.

Definition at line 40 of file adjoint.cpp.

◆ SPACE_DIM

constexpr int SPACE_DIM
constexpr

◆ young_modulus

constexpr double young_modulus = 1
constexpr

[Material properties for linear elasticity]

[Adjoint operator declarations]

Young's modulus E

Examples
mofem/tutorials/adv-0/plastic.cpp, and mofem/tutorials/vec-1/eigen_elastic.cpp.

Definition at line 37 of file adjoint.cpp.