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
	[Adjoint operator declarations] 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< ObjectiveFunctionData >	create_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:

• Adjoint-based sensitivity analysis for structural topology optimization
• Integration with Python for flexible objective function definition
• Higher-order finite element analysis with geometry fields
• Automatic differentiation using MoFEM's adjoint capabilities

  Author

  Anonymous author(s) committing under MIT license

Definition in file adjoint.cpp.

Typedef Documentation

BoundaryEle

using BoundaryEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle

Boundary finite elements.

Definition at line 48 of file adjoint.cpp.

BoundaryEleOp

using BoundaryEleOp = BoundaryEle::UserDataOperator

Boundary element operators

Definition at line 50 of file adjoint.cpp.

BoundaryLhsBCs

using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::BiLinearForm<I>

Boundary LHS natural BCs.

Definition at line 62 of file adjoint.cpp.

BoundaryRhsBCs

using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::LinearForm<I>

Boundary RHS natural BCs

Definition at line 60 of file adjoint.cpp.

DomainBaseOp

using DomainBaseOp = FormsIntegrators<DomainEleOp>::Assembly<A>::OpBase

[Postprocess results]

[calculateGradient]

Examples

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

Definition at line 1537 of file adjoint.cpp.

DomainEle

using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle

Domain finite elements.

Definition at line 47 of file adjoint.cpp.

DomainEleOp

using DomainEleOp = DomainEle::UserDataOperator

Domain element operators.

Definition at line 49 of file adjoint.cpp.

DomainRhsBCs

using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<A>::LinearForm<I>

Domain RHS natural BCs.

Definition at line 58 of file adjoint.cpp.

EntData

using EntData = EntitiesFieldData::EntData

Entity data for field operations.

Definition at line 46 of file adjoint.cpp.

OpBoundaryLhsBCs

using OpBoundaryLhsBCs = BoundaryLhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>

Boundary LHS flux operator.

Definition at line 63 of file adjoint.cpp.

OpBoundaryRhsBCs

using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>

Boundary flux operator.

Definition at line 61 of file adjoint.cpp.

OpDomainRhsBCs

using OpDomainRhsBCs = DomainRhsBCs::OpFlux<DomainBCs, 1, SPACE_DIM>

Domain flux operator.

Definition at line 59 of file adjoint.cpp.

PostProcEleBdy

using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy

Definition at line 81 of file adjoint.cpp.

PostProcEleDomain

using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain

Definition at line 79 of file adjoint.cpp.

SideEle

using SideEle = PostProcEleByDim<SPACE_DIM>::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_file

Path to Python file containing objective function definitions

Returns

boost::shared_ptr<ObjectiveFunctionData> Configured objective function interface

Exceptions

MoFEM

exception if Python initialization fails

Examples

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

Definition at line 2660 of file adjoint.cpp.

                                                    {
  auto ptr = boost::make_shared<ObjectiveFunctionDataImpl>();
  CHK_THROW_MESSAGE(ptr->initPython(py_file), "init python");
  return ptr;
}

diff_symmetrize()

template<int DIM>

auto diff_symmetrize ( FTensor::Number< DIM >  )

inline

Examples

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

Definition at line 1649 of file adjoint.cpp.

                                                                   {
 
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  FTensor::Index<'k', DIM> k;
  FTensor::Index<'l', DIM> l;
 
  FTensor::Tensor4<double, DIM, DIM, DIM, DIM> t_diff;
 
  t_diff(i, j, k, l) = 0;
  t_diff(0, 0, 0, 0) = 1;
  t_diff(1, 1, 1, 1) = 1;
 
  t_diff(1, 0, 1, 0) = 0.5;
  t_diff(1, 0, 0, 1) = 0.5;
 
  t_diff(0, 1, 0, 1) = 0.5;
  t_diff(0, 1, 1, 0) = 0.5;
 
  if constexpr (DIM == 3) {
    t_diff(2, 2, 2, 2) = 1;
 
    t_diff(2, 0, 2, 0) = 0.5;
    t_diff(2, 0, 0, 2) = 0.5;
    t_diff(0, 2, 0, 2) = 0.5;
    t_diff(0, 2, 2, 0) = 0.5;
 
    t_diff(2, 1, 2, 1) = 0.5;
    t_diff(2, 1, 1, 2) = 0.5;
    t_diff(1, 2, 1, 2) = 0.5;
    t_diff(1, 2, 2, 1) = 0.5;
  }
 
  return t_diff;
};

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 3283 of file adjoint.cpp.

                                                                {
  Range r;
 
  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
  auto bcs = mesh_mng->getCubitMeshsetPtr(
 
      std::regex((boost::format("%s(.*)") % block_name).str())
 
  );
 
  for (auto bc : bcs) {
    Range faces;
    CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
                                                       faces, true),
                   "get meshset ents");
    r.merge(faces);
  }
 
  for (auto dd = dim - 1; dd >= 0; --dd) {
    if (dd >= 0) {
      Range ents;
      CHK_MOAB_THROW(m_field.get_moab().get_adjacencies(r, dd, false, ents,
                                                        moab::Interface::UNION),
                     "get adjs");
      r.merge(ents);
    } else {
      Range verts;
      CHK_MOAB_THROW(m_field.get_moab().get_connectivity(r, verts),
                     "get verts");
      r.merge(verts);
    }
    CHK_MOAB_THROW(
        m_field.getInterface<CommInterface>()->synchroniseEntities(r), "comm");
  }
 
  return r;
};

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

argc	Command line argument count
argv	Command 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 2297 of file adjoint.cpp.

                                 {
 
  // Initialize Python environment for objective function interface
  Py_Initialize();
  np::initialize();
 
  // Initialize MoFEM/PETSc and MOAB data structures
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  auto core_log = logging::core::get();
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
 
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSync(), "FieldEvaluator"));
  LogManager::setLog("FieldEvaluator");
  MOFEM_LOG_TAG("FieldEvaluator", "field_eval");
 
  try {
 
    //! [Register MoFEM discrete manager in PETSc]
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    DMType dm_name_mg = "DMMOFEM_MG";
    CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
    //! [Register MoFEM discrete manager in PETSc
 
    //! [Create MoAB]
    moab::Core mb_instance;              ///< mesh database
    moab::Interface &moab = mb_instance; ///< mesh database interface
    //! [Create MoAB]
 
    //! [Create MoFEM]
    MoFEM::Core core(moab);           ///< finite element database
    MoFEM::Interface &m_field = core; ///< finite element database interface
    //! [Create MoFEM]
 
    //! [Example]
    Example ex(m_field);
    CHKERR ex.runProblem();
    //! [Example]
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
 
  if (Py_FinalizeEx() < 0) {
    exit(120);
  }
}

save_range()

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

Definition at line 3322 of file adjoint.cpp.

                                         {
  MoFEMFunctionBegin;
  auto out_meshset = get_temp_meshset_ptr(moab);
  CHKERR moab.add_entities(*out_meshset, r);
  CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
  MoFEMFunctionReturn(0);
}

Variable Documentation

A

constexpr AssemblyType A = AssemblyType::PETSC

constexpr

[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 = young_modulus / (3 * (1 - 2 * poisson_ratio))

constexpr

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

Definition at line 39 of file adjoint.cpp.

help

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

static

[calculateGradient]

Definition at line 2265 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

HenckyOps.hpp, 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 ν

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+ν))

Definition at line 40 of file adjoint.cpp.

SPACE_DIM

constexpr int SPACE_DIM

constexpr

Initial value:

=
    EXECUTABLE_DIMENSION

[Define dimension]

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

Examples

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

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 37 of file adjoint.cpp.