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)
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
	[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:

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

BoundaryEleOp

using BoundaryEleOp = BoundaryEle::UserDataOperator

Boundary element operators

Definition at line 58 of file adjoint.cpp.

BoundaryLhsBCs

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

Boundary LHS natural BCs.

Definition at line 76 of file adjoint.cpp.

BoundaryRhsBCs

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

Boundary RHS natural BCs.

Definition at line 71 of file adjoint.cpp.

DomainBaseOp

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

[Postprocess results]

[calculateGradient]

Examples

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

Definition at line 1341 of file adjoint.cpp.

DomainEle

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

Domain finite elements.

Definition at line 53 of file adjoint.cpp.

DomainEleOp

using DomainEleOp = DomainEle::UserDataOperator

Domain element operators.

Definition at line 57 of file adjoint.cpp.

DomainRhsBCs

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

Domain RHS natural BCs.

Definition at line 67 of file adjoint.cpp.

EntData

using EntData = EntitiesFieldData::EntData

Entity data for field operations.

Definition at line 51 of file adjoint.cpp.

OpBoundaryLhsBCs

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

Boundary LHS flux operator.

Definition at line 78 of file adjoint.cpp.

OpBoundaryRhsBCs

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

Boundary flux operator.

Definition at line 73 of file adjoint.cpp.

OpDomainRhsBCs

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

Domain flux operator.

Definition at line 69 of file adjoint.cpp.

PostProcEleBdy

using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy

Definition at line 98 of file adjoint.cpp.

PostProcEleDomain

using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain

Definition at line 96 of file adjoint.cpp.

SideEle

using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle

Definition at line 97 of file adjoint.cpp.

Enumeration Type Documentation

SensitivityMethod

enum SensitivityMethod

Enumerator
DIRECT
ADJOINT

Definition at line 107 of file adjoint.cpp.

{ DIRECT, ADJOINT };

Function Documentation

diff_symmetrize()

template<int DIM>

auto diff_symmetrize ( FTensor::Number< DIM >  )

inline

Examples

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

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

/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 2291 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 2238 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 2330 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

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:

=
    young_modulus /
    (3 * (1 - 2 * poisson_ratio))

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

Definition at line 40 of file adjoint.cpp.

derivative_type

SensitivityMethod derivative_type = ADJOINT

Examples

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

Definition at line 109 of file adjoint.cpp.

help

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

static

[calculateGradient]

Definition at line 2206 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:

=
    young_modulus / (2 * (1 + poisson_ratio))

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

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

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.