v0.16.0
Loading...
Searching...
No Matches
Classes | Typedefs | Enumerations | 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>
#include <ElasticPostProc.hpp>
#include <ObjectiveFunctionData.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  Example
 [Example] More...
 
struct  OpAdJointGradTimesSymTensor< SPACE_DIM, IntegrationType::GAUSS, DomainBaseOp >
 
struct  OpStateSensitivity
 
struct  OpAdJointObjective
 

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]
 

Enumerations

enum  SensitivityMethod { DIRECT , ADJOINT }
 

Functions

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)
 
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
 [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
 Bulk modulus K = E/(3(1-2ν))
 
constexpr double shear_modulus_G
 Shear modulus G = E/(2(1+ν))
 
PetscBool is_plane_strain
 
SensitivityMethod derivative_type = ADJOINT
 
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 55 of file adjoint.cpp.

◆ BoundaryEleOp

Boundary element operators

Definition at line 58 of file adjoint.cpp.

◆ BoundaryLhsBCs

Boundary LHS natural BCs.

Definition at line 76 of file adjoint.cpp.

◆ BoundaryRhsBCs

Boundary RHS natural BCs.

Definition at line 71 of file adjoint.cpp.

◆ DomainBaseOp

[Postprocess results]

[calculateGradient]

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

Definition at line 1341 of file adjoint.cpp.

◆ DomainEle

Domain finite elements.

Definition at line 53 of file adjoint.cpp.

◆ DomainEleOp

Domain element operators.

Definition at line 57 of file adjoint.cpp.

◆ DomainRhsBCs

Domain RHS natural BCs.

Definition at line 67 of file adjoint.cpp.

◆ EntData

Entity data for field operations.

Definition at line 51 of file adjoint.cpp.

◆ OpBoundaryLhsBCs

Boundary LHS flux operator.

Definition at line 78 of file adjoint.cpp.

◆ OpBoundaryRhsBCs

Boundary flux operator.

Definition at line 73 of file adjoint.cpp.

◆ OpDomainRhsBCs

Domain flux operator.

Definition at line 69 of file adjoint.cpp.

◆ PostProcEleBdy

Definition at line 98 of file adjoint.cpp.

◆ PostProcEleDomain

Definition at line 96 of file adjoint.cpp.

◆ SideEle

Definition at line 97 of file adjoint.cpp.

Enumeration Type Documentation

◆ SensitivityMethod

Enumerator
DIRECT 
ADJOINT 

Definition at line 107 of file adjoint.cpp.

107{ DIRECT, ADJOINT };
@ DIRECT
Definition adjoint.cpp:107
@ ADJOINT
Definition adjoint.cpp:107

Function Documentation

◆ get_range_from_block()

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

Definition at line 2254 of file adjoint.cpp.

2255 {
2256 Range r;
2257
2258 auto mesh_mng = m_field.getInterface<MeshsetsManager>();
2259 auto bcs = mesh_mng->getCubitMeshsetPtr(
2260
2261 std::regex((boost::format("%s(.*)") % block_name).str())
2262
2263 );
2264
2265 for (auto bc : bcs) {
2266 Range faces;
2267 CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
2268 faces, true),
2269 "get meshset ents");
2270 r.merge(faces);
2271 }
2272
2273 for (auto dd = dim - 1; dd >= 0; --dd) {
2274 if (dd >= 0) {
2275 Range ents;
2276 CHK_MOAB_THROW(m_field.get_moab().get_adjacencies(r, dd, false, ents,
2277 moab::Interface::UNION),
2278 "get adjs");
2279 r.merge(ents);
2280 } else {
2281 Range verts;
2282 CHK_MOAB_THROW(m_field.get_moab().get_connectivity(r, verts),
2283 "get verts");
2284 r.merge(verts);
2285 }
2287 m_field.getInterface<CommInterface>()->synchroniseEntities(r), "comm");
2288 }
2289
2290 return r;
2291};
#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:205
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 2201 of file adjoint.cpp.

2201 {
2202
2203 // Initialize Python environment for objective function interface
2204 Py_Initialize();
2205 np::initialize();
2206
2207 // Initialize MoFEM/PETSc and MOAB data structures
2208 const char param_file[] = "param_file.petsc";
2209 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
2210
2211 auto core_log = logging::core::get();
2212 core_log->add_sink(
2214
2215 core_log->add_sink(
2216 LogManager::createSink(LogManager::getStrmSync(), "FieldEvaluator"));
2217 LogManager::setLog("FieldEvaluator");
2218 MOFEM_LOG_TAG("FieldEvaluator", "field_eval");
2219
2220 try {
2221
2222 //! [Register MoFEM discrete manager in PETSc]
2223 DMType dm_name = "DMMOFEM";
2224 CHKERR DMRegister_MoFEM(dm_name);
2225 DMType dm_name_mg = "DMMOFEM_MG";
2227 //! [Register MoFEM discrete manager in PETSc
2228
2229 //! [Create MoAB]
2230 moab::Core mb_instance; ///< mesh database
2231 moab::Interface &moab = mb_instance; ///< mesh database interface
2232 //! [Create MoAB]
2233
2234 //! [Create MoFEM]
2235 MoFEM::Core core(moab); ///< finite element database
2236 MoFEM::Interface &m_field = core; ///< finite element database interface
2237 //! [Create MoFEM]
2238
2239 //! [Example]
2240
2241 Example ex(m_field);
2242 CHKERR ex.runProblem();
2243 //! [Example]
2244 }
2246
2248
2249 if (Py_FinalizeEx() < 0) {
2250 exit(120);
2251 }
2252}
static char help[]
[calculateGradient]
Definition adjoint.cpp:2169
#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:217
Core (interface) class.
Definition Core.hpp:83
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:68
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:123
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 2293 of file adjoint.cpp.

2294 {
2296 auto out_meshset = get_temp_meshset_ptr(moab);
2297 CHKERR moab.add_entities(*out_meshset, r);
2298 CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
2300};
#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
constexpr
Initial value:
=
AssemblyType::PETSC

[Define dimension]

Use PETSc for matrix/vector assembly

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
constexpr
Initial value:
=
(3 * (1 - 2 * poisson_ratio))
constexpr double poisson_ratio
Poisson's ratio ν
Definition adjoint.cpp:39
constexpr double young_modulus
[Material properties for linear elasticity]
Definition adjoint.cpp:38

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

Definition at line 40 of file adjoint.cpp.

◆ derivative_type

SensitivityMethod derivative_type = ADJOINT

◆ help

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

[calculateGradient]

Definition at line 2169 of file adjoint.cpp.

◆ I

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

Use Gauss quadrature for integration.

Definition at line 34 of file adjoint.cpp.

◆ is_plane_strain

PetscBool is_plane_strain
Initial value:
=
PETSC_FALSE

Flag for plane strain vs plane stress in 2D

Definition at line 46 of file adjoint.cpp.

◆ poisson_ratio

constexpr double poisson_ratio = 0.3
constexpr

Poisson's ratio ν

Definition at line 39 of file adjoint.cpp.

◆ shear_modulus_G

constexpr double shear_modulus_G
constexpr
Initial value:
=

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

Definition at line 43 of file adjoint.cpp.

◆ SPACE_DIM

constexpr int SPACE_DIM
constexpr
Initial value:
=
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13

[Define dimension]

Space dimension of problem (2D or 3D), set at compile time

Definition at line 28 of file adjoint.cpp.

◆ young_modulus

constexpr double young_modulus = 1
constexpr

[Material properties for linear elasticity]

Young's modulus E

Definition at line 38 of file adjoint.cpp.