mofem/tutorials/vec-7/adjoint.cpp

/**
 * @file adjoint.cpp  
 * @brief Topology optimization using adjoint method with Python objective functions
 * @details 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
 * @example mofem/tutorials/vec-7/adjoint.cpp
 *
 */
 
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
 
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
//! [Constants and material properties]
constexpr int BASE_DIM = 1; ///< Dimension of the base functions
 
//! [Define dimension]
constexpr int SPACE_DIM =
    EXECUTABLE_DIMENSION; ///< Space dimension of problem (2D or 3D), set at compile time
 
//! [Define dimension]
constexpr AssemblyType A = AssemblyType::PETSC; ///< Use PETSc for matrix/vector assembly
constexpr IntegrationType I =
    IntegrationType::GAUSS; ///< Use Gauss quadrature for integration
 
//! [Material properties for linear elasticity]
constexpr double young_modulus = 1;      ///< Young's modulus E
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; ///< Flag for plane strain vs plane stress in 2D
//! [Constants and material properties]
 
//! [Define finite element types and operators]
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; ///< Boundary element operators
//! [Define finite element types and operators]
 
//! [Boundary condition types]
struct DomainBCs {};   ///< Domain boundary conditions marker
struct BoundaryBCs {}; ///< Boundary conditions marker
 
//! [Natural boundary condition operators]
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
 
template <int DIM> struct PostProcEleByDim;
 
template <> struct PostProcEleByDim<2> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
};
 
template <> struct PostProcEleByDim<3> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
};
 
using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain;
using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle;
using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy;
 
//! [Adjoint operator declarations]
/// Forward declaration of adjoint test operator for geometric sensitivity
template <int SPACE_DIM, IntegrationType I, typename OpBase>
struct OpAdJointTestOp;
 
/// Forward declaration of operator for gradient times symmetric tensor operations
template <int SPACE_DIM, IntegrationType I, typename OpBase>  
struct OpAdJointGradTimesSymTensor;
//! [Adjoint operator declarations]
 
#include <ElasticSpring.hpp>
#include <FluidLevel.hpp>
#include <CalculateTraction.hpp>
#include <NaturalDomainBC.hpp>
#include <NaturalBoundaryBC.hpp>
#include <HookeOps.hpp>
 
/**
 * @brief Abstract interface for Python-defined objective functions
 * 
 * This interface allows users to define custom objective functions in Python
 * that can be called from C++ during the optimization process. The objective
 * function and its gradients are evaluated at integration points using:
 * - Coordinates of integration points
 * - Displacement field values
 * - Stress and strain tensors
 * 
 * Key methods:
 * - evalObjectiveFunction: Compute f(x,u,σ,ε) 
 * - evalObjectiveGradientStress: Compute ∂f/∂σ
 * - evalObjectiveGradientStrain: Compute ∂f/∂ε
 * - numberOfModes: Return number of topology modes for optimization
 * - blockModes: Define spatial topology modes for shape optimization
 */
struct ObjectiveFunctionData {
 
  /// Evaluate objective function f(coords, u, stress, strain)
  virtual 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) = 0;
 
  /// Evaluate gradient of objective function w.r.t. stress: ∂f/∂σ
  virtual 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) = 0;
 
  /// Evaluate gradient of objective function w.r.t. strain: ∂f/∂ε  
  virtual 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) = 0;
 
  /// Return number of topology optimization modes for given block
  virtual MoFEMErrorCode numberOfModes(int block_id, int &modes) = 0;
 
  /// Define spatial modes for topology optimization
  virtual MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords,
                                    std::array<double, 3> &centroid,
                                    std::array<double, 6> &bbox,
                                    MatrixDouble &o_ptr) = 0;
 
  virtual ~ObjectiveFunctionData() = default;
};
 
boost::shared_ptr<ObjectiveFunctionData>
create_python_objective_function(std::string py_file);
 
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);
/**
 * @brief Main class for topology optimization using adjoint method
 * 
 * This class implements a complete topology optimization workflow:
 * 1. Setup finite element problems for elasticity and adjoint analysis
 * 2. Define topology modes for design parametrization  
 * 3. Use TAO optimization library to minimize Python-defined objective function
 * 4. Compute sensitivities using adjoint method for efficient gradient calculation
 * 
 * The optimization process alternates between:
 * - Forward analysis: solve elastic problem for current geometry
 * - Adjoint analysis: solve adjoint problem for objective function gradients
 * - Geometry update: modify design using optimization algorithm (L-BFGS)
 */
struct Example {
 
  Example(MoFEM::Interface &m_field) : mField(m_field) {}
 
  /// Main driver function for the optimization process
  MoFEMErrorCode runProblem();
 
private:
  MoFEM::Interface &mField; ///< Reference to MoFEM interface
 
  boost::shared_ptr<MatrixDouble> vectorFieldPtr = nullptr; ///< Field values at evaluation points
 
  // Problem setup methods
  MoFEMErrorCode readMesh();          ///< Read mesh from file and setup meshsets
  MoFEMErrorCode setupProblem();      ///< Setup fields, approximation spaces and DOFs
  MoFEMErrorCode setupAdJoint();      ///< Setup adjoint fields and finite elements  
  MoFEMErrorCode boundaryCondition(); ///< Apply essential boundary conditions
  MoFEMErrorCode topologyModes();     ///< Compute topology optimization modes
  MoFEMErrorCode assembleSystem();    ///< Setup operators in finite element pipeline
  
  // Analysis methods  
  MoFEMErrorCode solveElastic();      ///< Solve forward elastic problem
  MoFEMErrorCode postprocessElastic(int iter, SmartPetscObj<Vec> adjoint_vector = nullptr); ///< Post-process and output results
  
  /// Calculate objective function gradient using adjoint method
  MoFEMErrorCode calculateGradient(PetscReal *objective_function_value,
                                   Vec objective_function_gradient,
                                   Vec adjoint_vector);
 
  // Problem configuration
  FieldApproximationBase base; ///< Choice of finite element basis functions  
  int fieldOrder = 2;          ///< Polynomial order for approximation
 
  // Solver and data management
  SmartPetscObj<KSP> kspElastic;  ///< Linear solver for elastic problem
  SmartPetscObj<DM> adjointDM;    ///< Data manager for adjoint problem
  boost::shared_ptr<ObjectiveFunctionData> pythonPtr; ///< Interface to Python objective function
  
  // Topology optimization data
  std::vector<SmartPetscObj<Vec>> modeVecs;           ///< Topology mode vectors (design variables)
  std::vector<std::array<double, 3>> modeCentroids;   ///< Centroids of optimization blocks  
  std::vector<std::array<double, 6>> modeBBoxes;      ///< Bounding boxes of optimization blocks
  SmartPetscObj<Vec> initialGeometry;                 ///< Initial geometry field
};
 
//! [Run topology optimization problem]
/**
 * @brief Main driver for topology optimization using adjoint sensitivity analysis
 * 
 * This function orchestrates the complete topology optimization workflow:
 * 1. Initialize Python objective function interface
 * 2. Setup finite element problems (forward and adjoint)
 * 3. Compute initial elastic solution
 * 4. Generate topology optimization modes 
 * 5. Run TAO optimization loop with adjoint-based gradients
 * 
 * The optimization uses TAO (Toolkit for Advanced Optimization) with L-BFGS
 * algorithm. At each iteration:
 * - Update geometry based on current design variables
 * - Solve forward elastic problem  
 * - Compute objective function and gradients using adjoint method
 * - Post-process results
 * 
 * @return MoFEMErrorCode Success or error code
 */
MoFEMErrorCode Example::runProblem() {
  MoFEMFunctionBegin;
 
  // Read objective function from Python file
  char objective_function_file_name[255] = "objective_function.py";
  CHKERR PetscOptionsGetString(
      PETSC_NULLPTR, PETSC_NULLPTR, "-objective_function",
      objective_function_file_name, 255, PETSC_NULLPTR);
      
  // Verify that the Python objective function file exists
  auto file_exists = [](std::string myfile) {
    std::ifstream file(myfile.c_str());
    if (file) {
      return true;
    }
    return false;
  };
  if (!file_exists(objective_function_file_name)) {
    MOFEM_LOG("WORLD", Sev::error) << "Objective function file NOT found: "
                                   << objective_function_file_name;
    CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "file NOT found");
  }
 
  // Create Python interface for objective function
  pythonPtr = create_python_objective_function(objective_function_file_name);
 
  // Setup finite element problems
  CHKERR readMesh();           // Read mesh and meshsets
  CHKERR setupProblem();       // Setup displacement field and geometry field  
  CHKERR setupAdJoint();       // Setup adjoint field for sensitivity analysis
  CHKERR boundaryCondition();  // Apply essential boundary conditions
  CHKERR assembleSystem();     // Setup finite element operators
 
  // Create linear solver for elastic problem
  auto pip = mField.getInterface<PipelineManager>();
  kspElastic = pip->createKSP();
  CHKERR KSPSetFromOptions(kspElastic);
 
  // Solve initial elastic problem
  CHKERR solveElastic();
  CHKERR postprocessElastic(-1); // Post-process initial solution
 
  auto create_vec_modes = [&](auto block_name) {
    MoFEMFunctionBegin;
    auto mesh_mng = mField.getInterface<MeshsetsManager>();
    auto bcs = mesh_mng->getCubitMeshsetPtr(
 
        std::regex((boost::format("%s(.*)") % block_name).str())
 
    );
 
    int nb_total_modes = 0;
    for (auto &bc : bcs) {
      auto id = bc->getMeshsetId();
      int nb_modes;
      CHKERR pythonPtr->numberOfModes(id, nb_modes);
      nb_total_modes += nb_modes;
    }
 
    MOFEM_LOG("WORLD", Sev::inform)
        << "Total number of modes to apply: " << nb_total_modes;
 
    modeVecs.resize(nb_total_modes);
 
    MoFEMFunctionReturn(0);
  };
 
  auto get_modes_bounding_boxes = [&](auto block_name) {
    MoFEMFunctionBegin;
    auto mesh_mng = mField.getInterface<MeshsetsManager>();
    auto bcs = mesh_mng->getCubitMeshsetPtr(
 
        std::regex((boost::format("%s(.*)") % block_name).str())
 
    );
 
    for (auto &bc : bcs) {
      auto meshset = bc->getMeshset();
      Range ents;
      CHKERR mField.get_moab().get_entities_by_handle(meshset, ents, true);
      Range verts;
      CHKERR mField.get_moab().get_connectivity(ents, verts, false);
      std::vector<double> x(verts.size());
      std::vector<double> y(verts.size());
      std::vector<double> z(verts.size());
      CHKERR mField.get_moab().get_coords(verts, x.data(), y.data(), z.data());
      std::array<double, 3> centroid = {0, 0, 0};
      for (int i = 0; i != verts.size(); ++i) {
        centroid[0] += x[i];
        centroid[1] += y[i];
        centroid[2] += z[i];
      }
      MPI_Allreduce(MPI_IN_PLACE, centroid.data(), 3, MPI_DOUBLE, MPI_SUM,
                    mField.get_comm());
      int nb_vertex = verts.size();
      MPI_Allreduce(MPI_IN_PLACE, &nb_vertex, 1, MPI_INT, MPI_SUM,
                    mField.get_comm());
      if (nb_vertex) {
        centroid[0] /= nb_vertex;
        centroid[1] /= nb_vertex;
        centroid[2] /= nb_vertex;
      }
      std::array<double, 6> bbox = {centroid[0], centroid[1], centroid[2],
                                    centroid[0], centroid[1], centroid[2]};
      for (int i = 0; i != verts.size(); ++i) {
        bbox[0] = std::min(bbox[0], x[i]);
        bbox[1] = std::min(bbox[1], y[i]);
        bbox[2] = std::min(bbox[2], z[i]);
        bbox[3] = std::max(bbox[3], x[i]);
        bbox[4] = std::max(bbox[4], y[i]);
        bbox[5] = std::max(bbox[5], z[i]);
      }
      MPI_Allreduce(MPI_IN_PLACE, &bbox[0], 3, MPI_DOUBLE, MPI_MIN,
                    mField.get_comm());
      MPI_Allreduce(MPI_IN_PLACE, &bbox[3], 3, MPI_DOUBLE, MPI_MAX,
                    mField.get_comm());
 
      MOFEM_LOG("WORLD", Sev::inform)
          << "Block: " << bc->getName() << " centroid: " << centroid[0] << " "
          << centroid[1] << " " << centroid[2];
      MOFEM_LOG("WORLD", Sev::inform)
          << "Block: " << bc->getName() << " bbox: " << bbox[0] << " "
          << bbox[1] << " " << bbox[2] << " " << bbox[3] << " " << bbox[4]
          << " " << bbox[5];
 
      modeCentroids.push_back(centroid);
      modeBBoxes.push_back(bbox);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto eval_objective_and_gradient = [](Tao tao, Vec x, PetscReal *f, Vec g,
                                        void *ctx) -> PetscErrorCode {
    MoFEMFunctionBegin;
    auto *ex_ptr = (Example *)ctx;
 
    int iter;
    CHKERR TaoGetIterationNumber(tao, &iter);
 
    auto set_geometry = [&](Vec x) {
      MoFEMFunctionBegin;
 
      VecScatter ctx;
      Vec x_local;
      CHKERR VecScatterCreateToAll(x, &ctx, &x_local);
      // scatter as many times as you need
      CHKERR VecScatterBegin(ctx, x, x_local, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecScatterEnd(ctx, x, x_local, INSERT_VALUES, SCATTER_FORWARD);
      // destroy scatter context and local vector when no longer needed
      CHKERR VecScatterDestroy(&ctx);
 
      auto current_geometry = vectorDuplicate(ex_ptr->initialGeometry);
      CHKERR VecCopy(ex_ptr->initialGeometry, current_geometry);
      const double *a;
      CHKERR VecGetArrayRead(x_local, &a);
      const double *coeff = a;
      for (auto &mode_vec : ex_ptr->modeVecs) {
        MOFEM_LOG("WORLD", Sev::verbose)
            << "Adding mode with coeff: " << *coeff;
        CHKERR VecAXPY(current_geometry, (*coeff), mode_vec);
        ++coeff;
      }
      CHKERR VecRestoreArrayRead(x_local, &a);
      CHKERR VecGhostUpdateBegin(current_geometry, INSERT_VALUES,
                                 SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(current_geometry, INSERT_VALUES,
                               SCATTER_FORWARD);
      CHKERR ex_ptr->mField.getInterface<VecManager>()
          ->setOtherLocalGhostVector("ADJOINT", "ADJOINT_FIELD", "GEOMETRY",
                                     RowColData::ROW, current_geometry,
                                     INSERT_VALUES, SCATTER_REVERSE);
 
      CHKERR VecDestroy(&x_local);
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_geometry(x);
    CHKERR KSPReset(ex_ptr->kspElastic);
    CHKERR ex_ptr->solveElastic();
    auto simple = ex_ptr->mField.getInterface<Simple>();
    auto adjoint_vector = createDMVector(simple->getDM());
    CHKERR ex_ptr->calculateGradient(f, g, adjoint_vector);
    CHKERR ex_ptr->postprocessElastic(iter, adjoint_vector);
 
    MoFEMFunctionReturn(0);
  };
 
  // Create topology optimization modes and setup TAO solver
  CHKERR create_vec_modes("OPTIMISE");
  CHKERR get_modes_bounding_boxes("OPTIMISE");
  
  /**
   * Setup TAO (Toolkit for Advanced Optimization) solver for topology optimization
   * 
   * TAO provides various optimization algorithms. Here we use TAOLMVM (Limited Memory
   * Variable Metric) which is a quasi-Newton method suitable for unconstrained
   * optimization with gradient information.
   * 
   * The optimization variables are coefficients for the topology modes, and
   * gradients are computed using the adjoint method for efficiency.
   */
  auto tao = createTao(mField.get_comm());
  CHKERR TaoSetType(tao, TAOLMVM);
 
  auto rank = mField.get_comm_rank();
  auto g = createVectorMPI(mField.get_comm(), (!rank) ? modeVecs.size() : 0,
                           PETSC_DECIDE);
  CHKERR TaoSetObjectiveAndGradient(tao, g, eval_objective_and_gradient,
                                    (void *)this);
 
  // Store initial geometry for reference during optimization
  initialGeometry = createDMVector(adjointDM);
  CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
      "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW, initialGeometry,
      INSERT_VALUES, SCATTER_FORWARD);
 
  // Setup topology optimization modes
  /**
   * Generate modes for topology optimization design parameterization
   * 
   * These modes represent perturbations to the geometry that can be used
   * as design variables in topology optimization. The modes are defined
   * through the Python interface and provide spatial basis functions
   * for shape modifications.
   */
  CHKERR topologyModes();
  
  // Initialize optimization variables and run solver
  /**
   * Start optimization with zero initial guess for design variables
   * 
   * The TAO solver will iteratively:
   * 1. Evaluate objective function at current design point
   * 2. Compute gradients using adjoint sensitivity analysis  
   * 3. Update design variables using L-BFGS algorithm
   * 4. Check convergence criteria
   * 
   * Each function evaluation involves solving the forward elastic problem
   * and the adjoint problem to compute sensitivities efficiently.
   */
  auto x0 = vectorDuplicate(g);
 
  CHKERR VecSet(x0, 0.0);
  CHKERR TaoSetSolution(tao, x0);
  CHKERR TaoSetFromOptions(tao);
  CHKERR TaoSolve(tao);
 
  // Optimization complete - results available in solution vectors
  MOFEM_LOG("WORLD", Sev::inform) << "Topology optimization completed";
 
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Read mesh]
/**
 * @brief Read mesh from file and setup material/boundary condition meshsets
 * 
 * This function loads the finite element mesh from file and processes
 * associated meshsets that define material properties and boundary conditions.
 * The mesh is typically generated using CUBIT and exported in .h5m format.
 * 
 * Meshsets are used to group elements/faces by:
 * - Material properties (for different material blocks)
 * - Boundary conditions (for applying loads and constraints)
 * - Optimization regions (for topology optimization)
 * 
 * @return MoFEMErrorCode Success or error code
 */
MoFEMErrorCode Example::readMesh() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  CHKERR simple->getOptions();        // Read command line options
  CHKERR simple->loadFile();          // Load mesh from file
  // Add meshsets if config file provided
  CHKERR mField.getInterface<MeshsetsManager>()->setMeshsetFromFile();
  MoFEMFunctionReturn(0);
}
//! [Read mesh]
 
//! [Set up problem]
/**
 * @brief Setup finite element fields, approximation spaces and degrees of freedom
 * 
 * This function configures the finite element problem by:
 * 1. Setting up the displacement field "U" with vector approximation
 * 2. Setting up the geometry field "GEOMETRY" for mesh deformation  
 * 3. Defining polynomial approximation order and basis functions
 * 4. Creating degrees of freedom on mesh entities
 * 
 * The displacement field uses H1 vector space for standard elasticity.
 * The geometry field allows mesh modification during topology optimization.
 * Different basis functions (Ainsworth-Legendre vs Demkowicz) can be selected.
 * 
 * @return MoFEMErrorCode Success or error code
 */
MoFEMErrorCode Example::setupProblem() {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  // Select basis functions for finite element approximation
  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
  PetscInt choice_base_value = AINSWORTH;
  CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
                              LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
 
  switch (choice_base_value) {
  case AINSWORTH:
    base = AINSWORTH_LEGENDRE_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set AINSWORTH_LEGENDRE_BASE for displacements";
    break;
  case DEMKOWICZ:
    base = DEMKOWICZ_JACOBI_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set DEMKOWICZ_JACOBI_BASE for displacements";
    break;
  default:
    base = LASTBASE;
    break;
  }
 
  // Add finite element fields
  /**
   * Setup displacement field "U" - the primary unknown in elasticity
   * This field represents displacement vector at each node/DOF
   */
  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &fieldOrder,
                            PETSC_NULLPTR);
 
  /**
   * Setup geometry field "GEOMETRY" - used for mesh deformation in optimization
   * This field stores current nodal coordinates and can be modified
   * during topology optimization to represent design changes
   */
  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
 
  // Set polynomial approximation order for both fields
  CHKERR simple->setFieldOrder("U", fieldOrder);
  CHKERR simple->setFieldOrder("GEOMETRY", fieldOrder);
  CHKERR simple->setUp();
 
  // Project higher-order geometry representation onto geometry field
  /**
   * For higher-order elements, this projects the exact geometry
   * onto the geometry field to maintain curved boundaries accurately
   */
  auto project_ho_geometry = [&]() {
    Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
    return mField.loop_dofs("GEOMETRY", ent_method);
  };
  CHKERR project_ho_geometry();
 
  // Check if plane strain assumption should be used
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-plane_strain",
                             &is_plane_strain, PETSC_NULLPTR);
 
  MoFEMFunctionReturn(0);
}
//! [Set up problem]
 
//! [Setup adjoint]
/**
 * @brief Setup adjoint fields and finite elements for sensitivity analysis
 * 
 * The adjoint method is used to efficiently compute gradients of the objective
 * function with respect to design variables. This function sets up:
 * 
 * 1. ADJOINT_FIELD - stores adjoint variables (Lagrange multipliers)
 * 2. ADJOINT_DM - data manager for adjoint problem
 * 3. Adjoint finite elements for domain and boundary
 * 
 * The adjoint equation is: K^T * λ = ∂f/∂u
 * where λ are adjoint variables, K is stiffness matrix, f is objective
 * 
 * @return MoFEMErrorCode Success or error code  
 */
MoFEMErrorCode Example::setupAdJoint() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  // Create adjoint data manager and field
  auto create_adjoint_dm = [&]() {
    auto adjoint_dm = createDM(mField.get_comm(), "DMMOFEM");
 
    auto add_field = [&]() {
      MoFEMFunctionBegin;
      CHKERR mField.add_field("ADJOINT_FIELD", H1, base, SPACE_DIM);
      CHKERR mField.add_ents_to_field_by_dim(simple->getMeshset(), SPACE_DIM,
                                             "ADJOINT_FIELD");
      for (auto tt = MBEDGE; tt <= moab::CN::TypeDimensionMap[SPACE_DIM].second;
           ++tt)
        CHKERR mField.set_field_order(simple->getMeshset(), tt, "ADJOINT_FIELD",
                                      fieldOrder);
      CHKERR mField.set_field_order(simple->getMeshset(), MBVERTEX,
                                    "ADJOINT_FIELD", 1);
      CHKERR mField.build_fields();
      MoFEMFunctionReturn(0);
    };
 
    auto add_adjoint_fe_impl = [&]() {
      MoFEMFunctionBegin;
      CHKERR mField.add_finite_element("ADJOINT_DOMAIN_FE");
      CHKERR mField.modify_finite_element_add_field_row("ADJOINT_DOMAIN_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_col("ADJOINT_DOMAIN_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_DOMAIN_FE",
                                                         "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_DOMAIN_FE",
                                                         "GEOMETRY");
      CHKERR mField.add_ents_to_finite_element_by_dim(
          simple->getMeshset(), SPACE_DIM, "ADJOINT_DOMAIN_FE");
      CHKERR mField.build_finite_elements("ADJOINT_DOMAIN_FE");
 
      CHKERR mField.add_finite_element("ADJOINT_BOUNDARY_FE");
      CHKERR mField.modify_finite_element_add_field_row("ADJOINT_BOUNDARY_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_col("ADJOINT_BOUNDARY_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_BOUNDARY_FE",
                                                         "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_BOUNDARY_FE",
                                                         "GEOMETRY");
 
      auto block_name = "OPTIMISE";
      auto mesh_mng = mField.getInterface<MeshsetsManager>();
      auto bcs = mesh_mng->getCubitMeshsetPtr(
 
          std::regex((boost::format("%s(.*)") % block_name).str())
 
      );
 
      for (auto bc : bcs) {
        CHKERR mField.add_ents_to_finite_element_by_dim(
            bc->getMeshset(), SPACE_DIM - 1, "ADJOINT_BOUNDARY_FE");
      }
 
      CHKERR mField.build_finite_elements("ADJOINT_BOUNDARY_FE");
 
      CHKERR mField.build_adjacencies(simple->getBitRefLevel(),
                                      simple->getBitRefLevelMask());
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_adjoint_dm_imp = [&]() {
      MoFEMFunctionBegin;
      CHKERR DMMoFEMCreateMoFEM(adjoint_dm, &mField, "ADJOINT",
                                simple->getBitRefLevel(),
                                simple->getBitRefLevelMask());
      CHKERR DMMoFEMSetDestroyProblem(adjoint_dm, PETSC_TRUE);
      CHKERR DMSetFromOptions(adjoint_dm);
      CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_DOMAIN_FE");
      CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_BOUNDARY_FE");
      CHKERR DMMoFEMSetSquareProblem(adjoint_dm, PETSC_TRUE);
      CHKERR DMMoFEMSetIsPartitioned(adjoint_dm, PETSC_TRUE);
      mField.getInterface<ProblemsManager>()->buildProblemFromFields =
          PETSC_TRUE;
      CHKERR DMSetUp(adjoint_dm);
      mField.getInterface<ProblemsManager>()->buildProblemFromFields =
          PETSC_FALSE;
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(add_field(), "add adjoint field");
    CHK_THROW_MESSAGE(add_adjoint_fe_impl(), "add adjoint fe");
    CHK_THROW_MESSAGE(set_adjoint_dm_imp(), "set adjoint dm");
 
    return adjoint_dm;
  };
 
  adjointDM = create_adjoint_dm();
 
  MoFEMFunctionReturn(0);
}
//
 
//! [Boundary condition]
MoFEMErrorCode Example::boundaryCondition() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
 
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
                                           "U", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
                                           "U", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
                                           "U", 2, 2);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                           "REMOVE_ALL", "U", 0, 3);
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "U");
 
  MoFEMFunctionReturn(0);
}
//! [Boundary condition]
 
//! [Adjoint modes]
MoFEMErrorCode Example::topologyModes() {
  MoFEMFunctionBegin;
 
  auto opt_ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
  auto subset_dm_bdy = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMSetSquareProblem(subset_dm_bdy, PETSC_TRUE);
  CHKERR DMMoFEMCreateSubDM(subset_dm_bdy, adjointDM, "SUBSET_BDY");
  CHKERR DMMoFEMAddElement(subset_dm_bdy, "ADJOINT_BOUNDARY_FE");
  CHKERR DMMoFEMAddSubFieldRow(subset_dm_bdy, "ADJOINT_FIELD",
                               boost::make_shared<Range>(opt_ents));
  CHKERR DMMoFEMAddSubFieldCol(subset_dm_bdy, "ADJOINT_FIELD",
                               boost::make_shared<Range>(opt_ents));
  CHKERR DMSetUp(subset_dm_bdy);
 
  auto subset_dm_domain = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMSetSquareProblem(subset_dm_domain, PETSC_TRUE);
  CHKERR DMMoFEMCreateSubDM(subset_dm_domain, adjointDM, "SUBSET_DOMAIN");
  CHKERR DMMoFEMAddElement(subset_dm_domain, "ADJOINT_DOMAIN_FE");
  CHKERR DMMoFEMAddSubFieldRow(subset_dm_domain, "ADJOINT_FIELD");
  CHKERR DMMoFEMAddSubFieldCol(subset_dm_domain, "ADJOINT_FIELD");
  CHKERR DMSetUp(subset_dm_domain);
 
  // remove dofs on boundary of the domain
  auto remove_dofs = [&]() {
    MoFEMFunctionBegin;
 
    std::array<Range, 3> remove_dim_ents;
    remove_dim_ents[0] =
        get_range_from_block(mField, "OPT_REMOVE_X", SPACE_DIM - 1);
    remove_dim_ents[1] =
        get_range_from_block(mField, "OPT_REMOVE_Y", SPACE_DIM - 1);
    remove_dim_ents[2] =
        get_range_from_block(mField, "OPT_REMOVE_Z", SPACE_DIM - 1);
 
    for (int d = 0; d != 3; ++d) {
      MOFEM_LOG("WORLD", Sev::inform)
          << "Removing topology modes on block OPT_REMOVE_" << (char)('X' + d)
          << " with " << remove_dim_ents[d].size() << " entities";
    }
 
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents,
                                                       true);
    auto skin = moab::Skinner(&mField.get_moab());
    Range boundary_ents;
    CHKERR skin.find_skin(0, body_ents, false, boundary_ents);
    for (int d = 0; d != 3; ++d) {
      boundary_ents = subtract(boundary_ents, remove_dim_ents[d]);
    }
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(boundary_ents,
                                 PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, nullptr);
    for (auto d = SPACE_DIM - 2; d >= 0; --d) {
      if (d >= 0) {
        Range ents;
        CHKERR mField.get_moab().get_adjacencies(boundary_ents, d, false, ents,
                                                 moab::Interface::UNION);
        boundary_ents.merge(ents);
      } else {
        Range verts;
        CHKERR mField.get_moab().get_connectivity(boundary_ents, verts);
        boundary_ents.merge(verts);
      }
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          boundary_ents);
    }
    boundary_ents.merge(opt_ents);
    CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
        "SUBSET_DOMAIN", "ADJOINT_FIELD", boundary_ents);
    for (int d = 0; d != 3; ++d) {
      CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
          "SUBSET_DOMAIN", "ADJOINT_FIELD", remove_dim_ents[d], d, d);
    }
 
    // #ifndef NDEBUG
    if (mField.get_comm_rank() == 0) {
      CHKERR save_range(mField.get_moab(), "topoMode_boundary_ents.vtk",
                        boundary_ents);
    }
    // #endif
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR remove_dofs();
 
  auto get_lhs_fe = [&]() {
    auto fe_lhs = boost::make_shared<BoundaryEle>(mField);
    fe_lhs->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe_lhs->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
                                                              "GEOMETRY");
    using OpMass = FormsIntegrators<BoundaryEleOp>::Assembly<A>::BiLinearForm<
        I>::OpMass<1, SPACE_DIM>;
    pip.push_back(new OpMass("ADJOINT_FIELD", "ADJOINT_FIELD",
                             [](double, double, double) { return 1.; }));
    return fe_lhs;
  };
 
  auto get_rhs_fe = [&]() {
    auto fe_rhs = boost::make_shared<BoundaryEle>(mField);
    fe_rhs->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe_rhs->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
                                                              "GEOMETRY");
 
    return fe_rhs;
  };
 
  auto block_name = "OPTIMISE";
  auto mesh_mng = mField.getInterface<MeshsetsManager>();
  auto bcs = mesh_mng->getCubitMeshsetPtr(
 
      std::regex((boost::format("%s(.*)") % block_name).str())
 
  );
 
  for (auto &v : modeVecs) {
    v = createDMVector(subset_dm_bdy);
  }
 
  using OP = FormsIntegrators<BoundaryEleOp>::Assembly<A>::OpBase;
  struct OpMode : public OP {
    OpMode(const std::string name,
           boost::shared_ptr<ObjectiveFunctionData> python_ptr, int id,
           std::vector<SmartPetscObj<Vec>> mode_vecs,
           std::vector<std::array<double, 3>> mode_centroids,
           std::vector<std::array<double, 6>> mode_bboxes, int block_counter,
           int mode_counter, boost::shared_ptr<Range> range = nullptr)
        : OP(name, name, OP::OPROW, range), pythonPtr(python_ptr), iD(id),
          modeVecs(mode_vecs), modeCentroids(mode_centroids),
          modeBboxes(mode_bboxes), blockCounter(block_counter),
          modeCounter(mode_counter) {}
 
    MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
      MoFEMFunctionBegin;
 
      if (OP::entsPtr) {
        if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
          MoFEMFunctionReturnHot(0);
      }
 
      auto nb_rows = data.getIndices().size();
      if (!nb_rows) {
        MoFEMFunctionReturnHot(0);
      }
      auto nb_base_functions = data.getN().size2();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      CHKERR pythonPtr->blockModes(iD, OP::getCoordsAtGaussPts(),
                                   modeCentroids[blockCounter],
                                   modeBboxes[blockCounter], blockModes);
 
      auto nb_integration_pts = getGaussPts().size2();
      if (blockModes.size2() != 3 * nb_integration_pts) {
        MOFEM_LOG("WORLD", Sev::error)
            << "Number of modes does not match number of integration points: "
            << blockModes.size2() << "!=" << 3 * nb_integration_pts;
        CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "modes/integration points");
      }
 
      VectorDouble nf(nb_rows);
 
      int nb_modes = blockModes.size1();
      for (auto mode = 0; mode != nb_modes; ++mode) {
        nf.clear();
        // get mode
        auto t_mode = getFTensor1FromPtr<3>(&blockModes(mode, 0));
        // get element volume
        const double vol = OP::getMeasure();
        // get integration weights
        auto t_w = OP::getFTensor0IntegrationWeight();
        // get base function gradient on rows
        auto t_base = data.getFTensor0N();
        // loop over integration points
        for (int gg = 0; gg != nb_integration_pts; gg++) {
 
          // take into account Jacobian
          const double alpha = t_w * vol;
          // loop over rows base functions
          auto t_nf = getFTensor1FromPtr<SPACE_DIM>(nf.data().data());
          int rr = 0;
          for (; rr != nb_rows / SPACE_DIM; ++rr) {
            t_nf(i) += alpha * t_base * t_mode(i);
            ++t_base;
            ++t_nf;
          }
          for (; rr < nb_base_functions; ++rr)
            ++t_base;
          ++t_w;    // move to another integration weight
          ++t_mode; // move to another mode
        }
        Vec vec = modeVecs[modeCounter + mode];
        auto size = data.getIndices().size();
        auto *indices = data.getIndices().data().data();
        auto *nf_data = nf.data().data();
        CHKERR VecSetValues(vec, size, indices, nf_data, ADD_VALUES);
      }
 
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
    MatrixDouble blockModes;
    std::vector<std::array<double, 3>> modeCentroids;
    std::vector<std::array<double, 6>> modeBboxes;
    int iD;
    std::vector<SmartPetscObj<Vec>> modeVecs;
    int blockCounter;
    int modeCounter;
  };
 
  auto solve_bdy = [&]() {
    MoFEMFunctionBegin;
 
    auto fe_lhs = get_lhs_fe();
    auto fe_rhs = get_rhs_fe();
    int block_counter = 0;
    int mode_counter = 0;
    for (auto &bc : bcs) {
      auto id = bc->getMeshsetId();
      Range ents;
      CHKERR mField.get_moab().get_entities_by_handle(bc->getMeshset(), ents,
                                                      true);
      auto range = boost::make_shared<Range>(ents);
      auto &pip_rhs = fe_rhs->getOpPtrVector();
      pip_rhs.push_back(new OpMode("ADJOINT_FIELD", pythonPtr, id, modeVecs,
                                   modeCentroids, modeBBoxes, block_counter,
                                   mode_counter, range));
      CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
                                      fe_rhs);
      pip_rhs.pop_back();
      int nb_modes;
      CHKERR pythonPtr->numberOfModes(id, nb_modes);
      ++block_counter;
      mode_counter += nb_modes;
      MOFEM_LOG("WORLD", Sev::inform)
          << "Setting mode block block: " << bc->getName()
          << " with ID: " << bc->getMeshsetId()
          << " total modes: " << mode_counter;
    }
 
    for (auto &v : modeVecs) {
      CHKERR VecAssemblyBegin(v);
      CHKERR VecAssemblyEnd(v);
      CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
    }
 
    auto M = createDMMatrix(subset_dm_bdy);
    fe_lhs->B = M;
    CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
                                    fe_lhs);
    CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
 
    auto solver = createKSP(mField.get_comm());
    CHKERR KSPSetOperators(solver, M, M);
    CHKERR KSPSetFromOptions(solver);
    CHKERR KSPSetUp(solver);
    auto v = createDMVector(subset_dm_bdy);
    for (auto &f : modeVecs) {
      CHKERR KSPSolve(solver, f, v);
      CHKERR VecSwap(f, v);
    }
 
    for (auto &v : modeVecs) {
      CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR solve_bdy();
 
  auto get_elastic_fe_lhs = [&]() {
    auto fe = boost::make_shared<DomainEle>(mField);
    fe->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
    CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
        mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
    return fe;
  };
 
  auto get_elastic_fe_rhs = [&]() {
    auto fe = boost::make_shared<DomainEle>(mField);
    fe->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
    CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
        mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
    return fe;
  };
 
  auto adjoint_gradient_postprocess = [&](auto mode) {
    MoFEMFunctionBegin;
    auto post_proc_mesh = boost::make_shared<moab::Core>();
    auto post_proc_begin =
        boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(mField,
                                                              post_proc_mesh);
    auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
        mField, post_proc_mesh);
 
    auto geom_vec = boost::make_shared<MatrixDouble>();
 
    auto post_proc_fe =
        boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        post_proc_fe->getOpPtrVector(), {H1}, "GEOMETRY");
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("ADJOINT_FIELD", geom_vec,
                                                    modeVecs[mode]));
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
            {},
 
            {{"MODE", geom_vec}},
 
            {},
 
            {}
 
            )
 
    );
 
    CHKERR DMoFEMPreProcessFiniteElements(adjointDM,
                                          post_proc_begin->getFEMethod());
    CHKERR DMoFEMLoopFiniteElements(adjointDM, "ADJOINT_DOMAIN_FE",
                                    post_proc_fe);
    CHKERR DMoFEMPostProcessFiniteElements(adjointDM,
                                           post_proc_begin->getFEMethod());
 
    CHKERR post_proc_end->writeFile("mode_" + std::to_string(mode) + ".h5m");
 
    MoFEMFunctionReturn(0);
  };
 
  auto solve_domain = [&]() {
    MoFEMFunctionBegin;
    auto fe_lhs = get_elastic_fe_lhs();
    auto fe_rhs = get_elastic_fe_rhs();
    auto v = createDMVector(subset_dm_domain);
    auto F = vectorDuplicate(v);
    fe_rhs->f = F;
 
    auto M = createDMMatrix(subset_dm_domain);
    fe_lhs->B = M;
    CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
                                    fe_lhs);
    CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
 
    auto solver = createKSP(mField.get_comm());
    CHKERR KSPSetOperators(solver, M, M);
    CHKERR KSPSetFromOptions(solver);
    CHKERR KSPSetUp(solver);
 
    int mode_counter = 0;
    for (auto &f : modeVecs) {
      CHKERR mField.getInterface<FieldBlas>()->setField(0, "ADJOINT_FIELD");
      CHKERR DMoFEMMeshToLocalVector(subset_dm_bdy, f, INSERT_VALUES,
                                     SCATTER_REVERSE);
      CHKERR VecZeroEntries(F);
      CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
                                      fe_rhs);
      CHKERR VecAssemblyBegin(F);
      CHKERR VecAssemblyEnd(F);
      CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
      CHKERR KSPSolve(solver, F, v);
      CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR DMoFEMMeshToLocalVector(subset_dm_domain, v, INSERT_VALUES,
                                     SCATTER_REVERSE);
      auto m = createDMVector(adjointDM);
      CHKERR DMoFEMMeshToLocalVector(adjointDM, m, INSERT_VALUES,
                                     SCATTER_FORWARD);
      f = m;
      ++mode_counter;
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR solve_domain();
 
  for (int i = 0; i < modeVecs.size(); ++i) {
    CHKERR adjoint_gradient_postprocess(i);
  }
 
  MoFEMFunctionReturn(0);
}
//! [Adjoint modes]
 
//! [Push operators to pipeline]
MoFEMErrorCode Example::assembleSystem() {
  MoFEMFunctionBegin;
  auto pip = mField.getInterface<PipelineManager>();
 
  //! [Integration rule]
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
 
  auto integration_rule_bc = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
 
  CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
  CHKERR pip->setDomainLhsIntegrationRule(integration_rule);
  CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
  CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
  //! [Integration rule]
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
 
  //! [Push domain stiffness matrix to pipeline]
  // Add LHS operator for elasticity (stiffness matrix)
  CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
      mField, pip->getOpDomainLhsPipeline(), "U", "MAT_ELASTIC", Sev::noisy);
  //! [Push domain stiffness matrix to pipeline]
 
  // Add RHS operator for internal forces
  CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
      mField, pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC", Sev::noisy);
 
  //! [Push Body forces]
  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
      pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
  //! [Push Body forces]
 
  //! [Push natural boundary conditions]
  // Add force boundary condition
  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
      pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
  // Add case for mix type of BCs
  CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
      pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::noisy);
  //! [Push natural boundary conditions]
  MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
 
//! [Solve]
MoFEMErrorCode Example::solveElastic() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto dm = simple->getDM();
  auto f = createDMVector(dm);
  auto d = vectorDuplicate(f);
  CHKERR VecZeroEntries(d);
  CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
 
  auto set_essential_bc = [&]() {
    MoFEMFunctionBegin;
    // This is low level pushing finite elements (pipelines) to solver
 
    auto ksp_ctx_ptr = getDMKspCtx(dm);
    auto pre_proc_rhs = boost::make_shared<FEMethod>();
    auto post_proc_rhs = boost::make_shared<FEMethod>();
    auto post_proc_lhs = boost::make_shared<FEMethod>();
 
    auto get_pre_proc_hook = [&]() {
      return EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_rhs,
                                                       {});
    };
    pre_proc_rhs->preProcessHook = get_pre_proc_hook();
 
    auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
                                                           post_proc_rhs, 1.)();
      MoFEMFunctionReturnHot(0);
    };
 
    auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(mField,
                                                           post_proc_lhs, 1.)();
      MoFEMFunctionReturnHot(0);
    };
 
    post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
    post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
 
    ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
    ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
    ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
    MoFEMFunctionReturn(0);
  };
 
  auto setup_and_solve = [&](auto solver) {
    MoFEMFunctionBegin;
    BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp";
    CHKERR KSPSetUp(solver);
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp <= Done";
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve";
    CHKERR KSPSolve(solver, f, d);
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve <= Done";
    MoFEMFunctionReturn(0);
  };
 
  MOFEM_LOG_CHANNEL("TIMER");
  MOFEM_LOG_TAG("TIMER", "timer");
 
  CHKERR set_essential_bc();
  CHKERR setup_and_solve(kspElastic);
 
  CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
 
  auto evaluate_field_at_the_point = [&]() {
    MoFEMFunctionBegin;
 
    int coords_dim = 3;
    std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
    PetscBool do_eval_field = PETSC_FALSE;
    CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
                                    field_eval_coords.data(), &coords_dim,
                                    &do_eval_field);
 
    if (do_eval_field) {
 
      vectorFieldPtr = boost::make_shared<MatrixDouble>();
      auto field_eval_data =
          mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
 
      CHKERR mField.getInterface<FieldEvaluatorInterface>()
          ->buildTree<SPACE_DIM>(field_eval_data, simple->getDomainFEName());
 
      field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
      auto no_rule = [](int, int, int) { return -1; };
      auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
      field_eval_fe_ptr->getRuleHook = no_rule;
 
      field_eval_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", vectorFieldPtr));
 
      CHKERR mField.getInterface<FieldEvaluatorInterface>()
          ->evalFEAtThePoint<SPACE_DIM>(
              field_eval_coords.data(), 1e-12, simple->getProblemName(),
              simple->getDomainFEName(), field_eval_data,
              mField.get_comm_rank(), mField.get_comm_rank(), nullptr, MF_EXIST,
              QUIET);
 
      if (vectorFieldPtr->size1()) {
        auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
        if constexpr (SPACE_DIM == 2)
          MOFEM_LOG("FieldEvaluator", Sev::inform)
              << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
        else
          MOFEM_LOG("FieldEvaluator", Sev::inform)
              << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
              << " U_Z: " << t_disp(2);
      }
 
      MOFEM_LOG_SYNCHRONISE(mField.get_comm());
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR evaluate_field_at_the_point();
 
  MoFEMFunctionReturn(0);
}
//! [Solve]
 
//! [Postprocess results]
MoFEMErrorCode Example::postprocessElastic(int iter,
                                           SmartPetscObj<Vec> adjoint_vector) {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>();
  auto det_ptr = boost::make_shared<VectorDouble>();
  auto jac_ptr = boost::make_shared<MatrixDouble>();
  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
  //! [Postprocess clean]
  pip->getDomainRhsFE().reset();
  pip->getDomainLhsFE().reset();
  pip->getBoundaryRhsFE().reset();
  pip->getBoundaryLhsFE().reset();
  //! [Postprocess clean]
 
  //! [Postprocess initialise]
  auto post_proc_mesh = boost::make_shared<moab::Core>();
  auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
      mField, post_proc_mesh);
  auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
      mField, post_proc_mesh);
  //! [Postprocess initialise]
 
  // lambada function to calculate displacement, strain and stress
  auto calculate_stress_ops = [&](auto &pip) {
    // Add H1 geometry ops
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
 
    // Use HookeOps commonDataFactory to get matStrainPtr and matCauchyStressPtr
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEle>(
        mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
 
    // Store U and GEOMETRY values if needed
    auto u_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
    auto x_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
 
    using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
 
    DataMapMat vec_map = {{"U", u_ptr}, {"GEOMETRY", x_ptr}};
    DataMapMat sym_map = {{"STRAIN", common_ptr->getMatStrain()},
                          {"STRESS", common_ptr->getMatCauchyStress()}};
 
    if (adjoint_vector) {
      auto adjoint_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
          "U", adjoint_ptr, adjoint_vector));
      vec_map["ADJOINT"] = adjoint_ptr;
    }
 
    // Return what you need: displacements, coordinates, strain, stress
    return boost::make_tuple(u_ptr, x_ptr, common_ptr->getMatStrain(),
                             common_ptr->getMatCauchyStress(), vec_map,
                             sym_map);
  };
 
  auto get_tag_id_on_pmesh = [&](bool post_proc_skin) {
    int def_val_int = 0;
    Tag tag_mat;
    CHKERR mField.get_moab().tag_get_handle(
        "MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
        MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
    auto meshset_vec_ptr =
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
            std::regex((boost::format("%s(.*)") % "MAT_ELASTIC").str()));
 
    Range skin_ents;
    std::unique_ptr<Skinner> skin_ptr;
    if (post_proc_skin) {
      skin_ptr = std::make_unique<Skinner>(&mField.get_moab());
      auto boundary_meshset = simple->getBoundaryMeshSet();
      CHKERR mField.get_moab().get_entities_by_handle(boundary_meshset,
                                                      skin_ents, true);
    }
 
    for (auto m : meshset_vec_ptr) {
      Range ents_3d;
      CHKERR mField.get_moab().get_entities_by_handle(m->getMeshset(), ents_3d,
                                                      true);
      int const id = m->getMeshsetId();
      ents_3d = ents_3d.subset_by_dimension(SPACE_DIM);
      if (post_proc_skin) {
        Range skin_faces;
        CHKERR skin_ptr->find_skin(0, ents_3d, false, skin_faces);
        ents_3d = intersect(skin_ents, skin_faces);
      }
      CHKERR mField.get_moab().tag_clear_data(tag_mat, ents_3d, &id);
    }
 
    return tag_mat;
  };
 
  auto post_proc_domain = [&](auto post_proc_mesh) {
    auto post_proc_fe =
        boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
        calculate_stress_ops(post_proc_fe->getOpPtrVector());
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
            {},
 
            vec_map,
 
            {},
 
            sym_map
 
            )
 
    );
 
    post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(false)});
    return post_proc_fe;
  };
 
  auto post_proc_boundary = [&](auto post_proc_mesh) {
    auto post_proc_fe =
        boost::make_shared<PostProcEleBdy>(mField, post_proc_mesh);
    AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
    auto op_loop_side =
        new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
    // push ops to side element, through op_loop_side operator
    auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
        calculate_stress_ops(op_loop_side->getOpPtrVector());
    post_proc_fe->getOpPtrVector().push_back(op_loop_side);
    auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
    post_proc_fe->getOpPtrVector().push_back(
        new ElasticOps::OpCalculateTraction(mat_stress_ptr, mat_traction_ptr));
    vec_map["T"] = mat_traction_ptr;
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
            {},
 
            vec_map,
 
            {},
 
            sym_map
 
            )
 
    );
 
    post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(true)});
    return post_proc_fe;
  };
 
  PetscBool post_proc_skin_only = PETSC_FALSE;
  if (SPACE_DIM == 3) {
    post_proc_skin_only = PETSC_TRUE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin_only",
                               &post_proc_skin_only, PETSC_NULLPTR);
  }
  if (post_proc_skin_only == PETSC_FALSE) {
    pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
  } else {
    pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
  }
 
  CHKERR DMoFEMPreProcessFiniteElements(simple->getDM(),
                                        post_proc_begin->getFEMethod());
  CHKERR pip->loopFiniteElements();
  CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
                                         post_proc_end->getFEMethod());
 
  CHKERR post_proc_end->writeFile("out_elastic_" + std::to_string(iter) +
                                  ".h5m");
  MoFEMFunctionReturn(0);
}
//! [Postprocess results]
 
//! [calculateGradient]
 
using DomainBaseOp = FormsIntegrators<DomainEleOp>::Assembly<A>::OpBase;
 
/**
 * @brief Adjoint test operator for elastic problems
 * 
 * This operator implements the adjoint (transpose) of the elastic stiffness operator.
 * It's used in the adjoint equation K^T * λ = ∂f/∂u to compute sensitivities.
 * 
 * The operator computes: ∫_Ω ∇v : σ dΩ
 * where v are test functions and σ is the Cauchy stress tensor.
 * 
 * This is the transpose of the standard elastic operator ∫_Ω ∇u : C : ∇v dΩ
 * 
 * @tparam SPACE_DIM Spatial dimension (2 or 3)
 */
template <int SPACE_DIM>
struct OpAdJointTestOp<SPACE_DIM, IntegrationType::GAUSS, DomainBaseOp>
    : public DomainBaseOp {
 
  using OP = DomainBaseOp;
 
  OpAdJointTestOp(const std::string field_name,
                  boost::shared_ptr<HookeOps::CommonData> comm_ptr)
      : OP(field_name, field_name, DomainEleOp::OPROW), commPtr(comm_ptr) {}
 
protected:
  boost::shared_ptr<HookeOps::CommonData> commPtr;
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
};
 
/**
 * @brief Integration implementation for adjoint test operator
 * 
 * Computes the integral: ∫_Ω ∇v : σ dΩ
 * where v are test functions and σ is Cauchy stress tensor
 * 
 * This forms the right-hand side of the adjoint equation K^T * λ = ∂f/∂u
 * 
 * @param row_data Test function data for current element
 * @return MoFEMErrorCode Success or error code
 */
template <int SPACE_DIM>
MoFEMErrorCode
OpAdJointTestOp<SPACE_DIM, IntegrationType::GAUSS, DomainBaseOp>::iNtegrate(
    EntitiesFieldData::EntData &row_data) {
  MoFEMFunctionBegin;
  // Define tensor indices for Einstein notation
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, I);
  FTENSOR_INDEX(SPACE_DIM, J);
 
  // Get element geometric information
  const double vol = OP::getMeasure();
  // Get integration weights for Gauss quadrature
  auto t_w = OP::getFTensor0IntegrationWeight();
  // Get test function gradients ∇v
  auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
  // Get Cauchy stress tensor σ from forward solution
  auto t_cauchy_stress =
      getFTensor2SymmetricFromMat<SPACE_DIM>(*(commPtr->getMatCauchyStress()));
  
  // Loop over integration points for numerical integration
  for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
    // take into account Jacobian
    const double alpha = t_w * vol;
 
    // get rhs vector
    auto t_nf = OP::template getNf<SPACE_DIM>();
    // loop over rows base functions
    int rr = 0;
    for (; rr != OP::nbRows / SPACE_DIM; rr++) {
      // Variation due to change in geometry (diff base)
      t_nf(j) += alpha * t_row_grad(i) * t_cauchy_stress(i, j);
 
      ++t_row_grad;
      ++t_nf;
    }
    for (; rr < OP::nbRowBaseFunctions; ++rr) {
      ++t_row_grad;
    }
 
    ++t_cauchy_stress;
    ++t_w; // move to another integration weight
  }
  MoFEMFunctionReturn(0);
}
 
template <int SPACE_DIM>
struct OpAdJointGradTimesSymTensor<SPACE_DIM, IntegrationType::GAUSS,
                                   DomainBaseOp> : public DomainBaseOp {
 
  using OP = DomainBaseOp;
 
  OpAdJointGradTimesSymTensor(const std::string field_name,
                              boost::shared_ptr<HookeOps::CommonData> comm_ptr,
                              boost::shared_ptr<MatrixDouble> jac,
                              boost::shared_ptr<MatrixDouble> diff_jac,
                              boost::shared_ptr<VectorDouble> cof_vals)
      : OP(field_name, field_name, DomainEleOp::OPROW), commPtr(comm_ptr),
        jac(jac), diffJac(diff_jac), cofVals(cof_vals) {}
 
protected:
  boost::shared_ptr<HookeOps::CommonData> commPtr;
  boost::shared_ptr<MatrixDouble> jac;
  boost::shared_ptr<MatrixDouble> diffJac;
  boost::shared_ptr<VectorDouble> cofVals;
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
};
 
template <int DIM> inline auto diff_symmetrize(FTensor::Number<DIM>) {
 
  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;
};
 
template <int SPACE_DIM>
MoFEMErrorCode
OpAdJointGradTimesSymTensor<SPACE_DIM, IntegrationType::GAUSS, DomainBaseOp>::
    iNtegrate(EntitiesFieldData::EntData &row_data) {
  MoFEMFunctionBegin;
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, I);
  FTENSOR_INDEX(SPACE_DIM, J);
 
  auto t_diff_symm = diff_symmetrize(FTensor::Number<SPACE_DIM>());
 
  // get element volume
  const double vol = OP::getMeasure();
  // get integration weights
  auto t_w = OP::getFTensor0IntegrationWeight();
  // get Jacobian values
  auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(jac));
  // get diff Jacobian values
  auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(diffJac));
  // get cofactor values
  auto t_cof = getFTensor0FromVec(*(cofVals));
  // get base function gradient on rows
  auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
  // get fradient of the field
  auto t_grad_u =
      getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(commPtr->matGradPtr));
  // get field gradient values
  auto t_cauchy_stress =
      getFTensor2SymmetricFromMat<SPACE_DIM>(*(commPtr->getMatCauchyStress()));
  // material stiffness tensor
  auto t_D =
      getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
  // loop over integration points
  for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
    // take into account Jacobian
    const double alpha = t_w * vol;
 
    auto t_det = determinantTensor(t_jac);
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_inv_jac;
    CHKERR invertTensor(t_jac, t_det, t_inv_jac);
 
    // Calculate the variation of the gradient due to geometry change
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_inv_jac;
    t_diff_inv_jac(i, j) =
        -(t_inv_jac(i, I) * t_diff_jac(I, J)) * t_inv_jac(J, j);
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_grad;
    t_diff_grad(i, j) = t_grad_u(i, k) * t_diff_inv_jac(k, j);
 
    // Calculate the variation of the strain tensor
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_strain;
    t_diff_strain(i, j) = t_diff_symm(i, j, k, l) * t_diff_grad(k, l);
 
    // Calculate the variation of the stress tensor
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_stress;
    t_diff_stress(i, j) = t_D(i, j, k, l) * t_diff_strain(k, l);
 
    // get rhs vector
    auto t_nf = OP::template getNf<SPACE_DIM>();
    // loop over rows base functions
    int rr = 0;
    for (; rr != OP::nbRows / SPACE_DIM; rr++) {
 
      FTensor::Tensor1<double, SPACE_DIM> t_diff_row_grad;
      t_diff_row_grad(k) = t_row_grad(j) * t_diff_inv_jac(j, k);
 
      // Variation due to change in geometry (diff base)
      t_nf(j) += alpha * t_diff_row_grad(i) * t_cauchy_stress(i, j);
 
      // Variation due to change in domain (cofactor)
      t_nf(j) += (alpha * t_cof) * t_row_grad(i) * t_cauchy_stress(i, j);
 
      // Variation due to change in stress (diff stress)
      t_nf(j) += alpha * t_row_grad(i) * t_diff_stress(i, j);
 
      ++t_row_grad;
      ++t_nf;
    }
    for (; rr < OP::nbRowBaseFunctions; ++rr) {
      ++t_row_grad;
    }
 
    ++t_grad_u;
    ++t_cauchy_stress;
    ++t_jac;
    ++t_diff_jac;
    ++t_cof;
    ++t_w; // move to another integration weight
  }
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Operator for computing objective function and gradients in topology optimization
 * 
 * This operator interfaces with Python-defined objective functions to:
 * 1. Evaluate objective function f(u, x) at current design point
 * 2. Compute gradients ∂f/∂u and ∂f/∂x for adjoint sensitivity analysis
 * 
 * The objective function typically depends on:
 * - u: displacement field from forward elastic solution
 * - x: design variables (topology parameters)
 * - σ: stress field computed from displacement
 * 
 * This enables flexible objective function definition (compliance, stress
 * constraints, volume constraints, etc.) through Python scripting.
 */
struct OpAdJointObjective : public ForcesAndSourcesCore::UserDataOperator {
  using OP = ForcesAndSourcesCore::UserDataOperator;
  
  OpAdJointObjective(boost::shared_ptr<ObjectiveFunctionData> python_ptr,
                     boost::shared_ptr<HookeOps::CommonData> comm_ptr,
                     boost::shared_ptr<MatrixDouble> jac_ptr,
                     boost::shared_ptr<MatrixDouble> diff_jac,
                     boost::shared_ptr<VectorDouble> cof_vals,
                     boost::shared_ptr<MatrixDouble> d_grad_ptr,
                     boost::shared_ptr<MatrixDouble> u_ptr,
                     boost::shared_ptr<double> glob_objective_ptr,
                     boost::shared_ptr<double> glob_objective_grad_ptr)
      : OP(NOSPACE, OP::OPSPACE), pythonPtr(python_ptr), commPtr(comm_ptr),
        jacPtr(jac_ptr), diffJacPtr(diff_jac), cofVals(cof_vals),
        dGradPtr(d_grad_ptr), uPtr(u_ptr), globObjectivePtr(glob_objective_ptr),
        globObjectiveGradPtr(glob_objective_grad_ptr) {}
 
  /**
   * @brief Compute objective function contributions at element level
   * 
   * Evaluates Python objective function with current displacement and stress
   * state, and accumulates global objective value and gradients.
   */
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    // Define tensor indices for calculations
    FTENSOR_INDEX(SPACE_DIM, i);
    FTENSOR_INDEX(SPACE_DIM, j);
    FTENSOR_INDEX(SPACE_DIM, k);
    FTENSOR_INDEX(SPACE_DIM, l);
 
    FTENSOR_INDEX(SPACE_DIM, I);
    FTENSOR_INDEX(SPACE_DIM, J);
 
    constexpr auto symm_size = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
 
    auto t_diff_symm = diff_symmetrize(FTensor::Number<SPACE_DIM>());
 
    auto nb_gauss_pts = getGaussPts().size2();
    auto objective_ptr = boost::make_shared<VectorDouble>(nb_gauss_pts);
    auto objective_dstress =
        boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
    auto objective_dstrain =
        boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
    auto obj_grad = boost::make_shared<MatrixDouble>(SPACE_DIM, nb_gauss_pts);
 
    auto evaluate_python = [&]() {
      MoFEMFunctionBegin;
      auto &coords = OP::getCoordsAtGaussPts();
      CHKERR pythonPtr->evalObjectiveFunction(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_ptr);
      CHKERR pythonPtr->evalObjectiveGradientStress(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_dstress);
      CHKERR pythonPtr->evalObjectiveGradientStrain(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_dstrain);
 
      auto t_grad_u =
          getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(commPtr->matGradPtr));
      auto t_D =
          getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
      auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(jacPtr));
      auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(diffJacPtr));
      auto t_cof = getFTensor0FromVec(*(cofVals));
      auto t_d_grad = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(dGradPtr));
 
      auto t_obj = getFTensor0FromVec(*objective_ptr);
      auto t_obj_dstress =
          getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstress);
      auto t_obj_dstrain =
          getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstrain);
 
      auto vol = OP::getMeasure();
      auto t_w = getFTensor0IntegrationWeight();
      for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
 
        auto t_det = determinantTensor(t_jac);
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_inv_jac;
        CHKERR invertTensor(t_jac, t_det, t_inv_jac);
 
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_inv_jac;
        t_diff_inv_jac(i, j) =
            -(t_inv_jac(i, I) * t_diff_jac(I, J)) * t_inv_jac(J, j);
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_grad;
        t_diff_grad(i, j) = t_grad_u(i, k) * t_diff_inv_jac(k, j);
 
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_d_strain;
        t_d_strain(i, j) = t_diff_symm(i, j, k, l) * (
 
                                                         t_d_grad(k, l)
 
                                                         +
 
                                                         t_diff_grad(k, l)
 
                                                     );
 
        auto alpha = t_w * vol;
 
        (*globObjectivePtr) += alpha * t_obj;
        (*globObjectiveGradPtr) +=
            alpha *
            (
 
                t_obj_dstress(i, j) * (t_D(i, j, k, l) * t_d_strain(k, l))
 
                +
 
                t_obj_dstrain(i, j) * t_d_strain(i, j)
 
                +
 
                t_obj * t_cof
 
            );
 
        ++t_w;
        ++t_jac;
        ++t_diff_jac;
        ++t_cof;
 
        ++t_obj;
        ++t_obj_dstress;
        ++t_obj_dstrain;
 
        ++t_grad_u;
        ++t_d_grad;
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR evaluate_python();
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
  boost::shared_ptr<HookeOps::CommonData> commPtr;
  boost::shared_ptr<MatrixDouble> jacPtr;
  boost::shared_ptr<MatrixDouble> diffJacPtr;
  boost::shared_ptr<VectorDouble> cofVals;
  boost::shared_ptr<MatrixDouble> dGradPtr;
  boost::shared_ptr<MatrixDouble> uPtr;
 
  boost::shared_ptr<double> globObjectivePtr;
  boost::shared_ptr<double> globObjectiveGradPtr;
};
 
MoFEMErrorCode Example::calculateGradient(PetscReal *objective_function_value,
                                          Vec objective_function_gradient,
                                          Vec adjoint_vector) {
  MOFEM_LOG_CHANNEL("WORLD");
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  auto ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
 
  auto get_essential_fe = [this]() {
    auto post_proc_rhs = boost::make_shared<FEMethod>();
    auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
                                                           post_proc_rhs, 0)();
      MoFEMFunctionReturnHot(0);
    };
    post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
    return post_proc_rhs;
  };
 
  auto get_fd_direvative_fe = [&]() {
    auto fe = boost::make_shared<DomainEle>(mField);
    fe->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe->getOpPtrVector();
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
    // Add RHS operators for internal forces
    constexpr bool debug = false;
    if constexpr (debug) {
      auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
          mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
      pip.push_back(
          new OpAdJointTestOp<SPACE_DIM, I, DomainBaseOp>("U", common_ptr));
    } else {
      HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
          mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
    }
    return fe;
  };
 
  auto calulate_fd_residual = [&](auto eps, auto diff_vec, auto fd_vec) {
    MoFEMFunctionBegin;
 
    constexpr bool debug = false;
 
    auto geom_norm = [](MoFEM::Interface &mField) {
      MoFEMFunctionBegin;
      auto field_blas = mField.getInterface<FieldBlas>();
      double nrm2 = 0.0;
      auto norm2_field = [&](const double val) {
        nrm2 += val * val;
        return val;
      };
      CHKERR field_blas->fieldLambdaOnValues(norm2_field, "GEOMETRY");
      MPI_Allreduce(MPI_IN_PLACE, &nrm2, 1, MPI_DOUBLE, MPI_SUM,
                    mField.get_comm());
      MOFEM_LOG("WORLD", Sev::inform) << "Geometry norm: " << sqrt(nrm2);
      MoFEMFunctionReturn(0);
    };
 
    if constexpr (debug)
      CHKERR geom_norm(mField);
 
    auto initial_current_geometry = createDMVector(adjointDM);
    CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
        "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
        initial_current_geometry, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecAssemblyBegin(initial_current_geometry);
    CHKERR VecAssemblyEnd(initial_current_geometry);
 
    if constexpr (debug)
      CHKERR geom_norm(mField);
 
    auto perturb_geometry = [&](auto eps, auto diff_vec) {
      MoFEMFunctionBegin;
      auto current_geometry = vectorDuplicate(initial_current_geometry);
      CHKERR VecCopy(initial_current_geometry, current_geometry);
      CHKERR VecAXPY(current_geometry, eps, diff_vec);
      CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
          "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
          current_geometry, INSERT_VALUES, SCATTER_REVERSE);
      MoFEMFunctionReturn(0);
    };
 
    auto fe = get_fd_direvative_fe();
    auto fp = vectorDuplicate(diff_vec);
    auto fm = vectorDuplicate(diff_vec);
    auto calc_impl = [&](auto f, auto eps) {
      MoFEMFunctionBegin;
      CHKERR VecZeroEntries(f);
      fe->f = f;
      CHKERR perturb_geometry(eps, diff_vec);
      CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
                                      simple->getDomainFEName(), fe);
      CHKERR VecAssemblyBegin(f);
      CHKERR VecAssemblyEnd(f);
      CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
      auto post_proc_rhs = get_essential_fe();
      post_proc_rhs->f = f;
      CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
                                             post_proc_rhs.get());
      MoFEMFunctionReturn(0);
    };
    CHKERR calc_impl(fp, eps);
    CHKERR calc_impl(fm, -eps);
    CHKERR VecWAXPY(fd_vec, -1.0, fm, fp);
    CHKERR VecScale(fd_vec, 1.0 / (2.0 * eps));
 
    CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
        "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
        initial_current_geometry, INSERT_VALUES, SCATTER_REVERSE);
 
    if constexpr (debug)
      CHKERR geom_norm(mField);
 
    MoFEMFunctionReturn(0);
  };
 
  auto get_direvative_fe = [&](auto diff_vec) {
    auto fe_adjoint = boost::make_shared<DomainEle>(mField);
    fe_adjoint->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe_adjoint->getOpPtrVector();
 
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
    auto cof_ptr = boost::make_shared<VectorDouble>();
 
    using OpCoFactor =
        AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
 
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
 
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "GEOMETRY", jac_ptr));
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "U", diff_jac_ptr,
        diff_vec)); // Note: that vector is stored on displacemnt vector, that
                    // why is used here
 
    pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
 
    // Add RHS operators for internal forces
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
        mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
    pip.push_back(new OpAdJointGradTimesSymTensor<SPACE_DIM, I, DomainBaseOp>(
        "U", common_ptr, jac_ptr, diff_jac_ptr, cof_ptr));
 
    return fe_adjoint;
  };
 
  auto get_objective_fe = [&](auto diff_vec, auto grad_vec,
                              auto glob_objective_ptr,
                              auto glob_objective_grad_ptr) {
    auto fe_adjoint = boost::make_shared<DomainEle>(mField);
    fe_adjoint->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data +  p_data - 1;
    };
    auto &pip = fe_adjoint->getOpPtrVector();
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
 
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
    auto cof_ptr = boost::make_shared<VectorDouble>();
    auto d_grad_ptr = boost::make_shared<MatrixDouble>();
    auto u_ptr = boost::make_shared<MatrixDouble>();
 
    using OpCoFactor =
        AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "GEOMETRY", jac_ptr));
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "U", diff_jac_ptr,
        diff_vec)); // Note: that vector is stored on displacemnt vector, that
                    // why is used here
    pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "U", d_grad_ptr, grad_vec));
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
 
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
        mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
    pip.push_back(new OpAdJointObjective(
        pythonPtr, common_ptr, jac_ptr, diff_jac_ptr, cof_ptr, d_grad_ptr,
        u_ptr, glob_objective_ptr, glob_objective_grad_ptr));
 
    return fe_adjoint;
  };
 
  auto dm = simple->getDM();
  auto f = createDMVector(dm);
  auto d = vectorDuplicate(f);
  auto dm_diff_vec = vectorDuplicate(d);
 
  auto adjoint_fe = get_direvative_fe(dm_diff_vec);
  auto glob_objective_ptr = boost::make_shared<double>(0.0);
  auto glob_objective_grad_ptr = boost::make_shared<double>(0.0);
  auto objective_fe = get_objective_fe(dm_diff_vec, d, glob_objective_ptr,
                                       glob_objective_grad_ptr);
 
  auto set_variance_of_geometry = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    CHKERR DMoFEMMeshToLocalVector(adjointDM, mod_vec, INSERT_VALUES,
                                   SCATTER_REVERSE);
    CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
        simple->getProblemName(), "U", "ADJOINT_FIELD", RowColData::ROW,
        dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateBegin(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_variance_internal_forces = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    CHKERR VecZeroEntries(f);
    CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
    adjoint_fe->f = f;
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), adjoint_fe);
    CHKERR VecAssemblyBegin(f);
    CHKERR VecAssemblyEnd(f);
    CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
    auto post_proc_rhs = get_essential_fe();
    post_proc_rhs->f = f;
    CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
                                           post_proc_rhs.get());
    CHKERR VecScale(f, -1.0);
 
#ifndef NDEBUG
    constexpr bool debug = false;
    if constexpr (debug) {
      double norm0;
      CHKERR VecNorm(f, NORM_2, &norm0);
      auto fd_check = vectorDuplicate(f);
      double eps = 1e-5;
      CHKERR calulate_fd_residual(eps, dm_diff_vec, fd_check);
      double nrm;
      CHKERR VecAXPY(fd_check, -1.0, f);
      CHKERR VecNorm(fd_check, NORM_2, &nrm);
      MOFEM_LOG("WORLD", Sev::inform)
          << " FD check for internal forces [ " << mode << " ]: " << nrm
          << " / " << norm0 << " ( " << (nrm / norm0) << " )";
    }
#endif
 
    MoFEMFunctionReturn(0);
  };
 
  auto calulate_variance_of_displacement = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    CHKERR KSPSolve(kspElastic, f, d);
    CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_variance_of_objective_function = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    *glob_objective_ptr = 0.0;
    *glob_objective_grad_ptr = 0.0;
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                    objective_fe);
    CHKERR VecSetValue(objective_function_gradient, mode,
                       *glob_objective_grad_ptr, ADD_VALUES);
    std::array<double, 2> array = {*glob_objective_ptr,
                                   *glob_objective_grad_ptr};
    MPI_Allreduce(MPI_IN_PLACE, array.data(), 2, MPI_DOUBLE, MPI_SUM,
                  mField.get_comm());
    *glob_objective_ptr = array[0];
    *glob_objective_grad_ptr = array[1];
    (*objective_function_value) += *glob_objective_ptr;
    MOFEM_LOG_CHANNEL("WORLD");
    MOFEM_LOG("WORLD", Sev::verbose) << " Objective gradient [ " << mode
                                     << " ]: " << *glob_objective_grad_ptr;
    MoFEMFunctionReturn(0);
  };
 
  CHKERR VecZeroEntries(objective_function_gradient);
  CHKERR VecZeroEntries(adjoint_vector);
  *objective_function_value = 0.0;
  int mode = 0;
  for (auto mod_vec : modeVecs) {
 
    CHKERR set_variance_of_geometry(mode, mod_vec);
    CHKERR calculate_variance_internal_forces(mode, mod_vec);
    CHKERR calulate_variance_of_displacement(mode, mod_vec);
    CHKERR calculate_variance_of_objective_function(mode, mod_vec);
 
    CHKERR VecAXPY(adjoint_vector, *glob_objective_grad_ptr, dm_diff_vec);
    ++mode;
  }
 
  (*objective_function_value) /= modeVecs.size();
  MOFEM_LOG("WORLD", Sev::verbose)
      << "Objective function: " << *glob_objective_ptr;
 
  CHKERR VecAssemblyBegin(objective_function_gradient);
  CHKERR VecAssemblyEnd(objective_function_gradient);
 
  CHKERR VecAssemblyBegin(adjoint_vector);
  CHKERR VecAssemblyEnd(adjoint_vector);
  CHKERR VecGhostUpdateBegin(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
 
  MoFEMFunctionReturn(0);
}
//! [calculateGradient]
 
static char help[] = "...\n\n";
 
/**
 * @brief 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)
 * 
 * @param argc Command line argument count
 * @param argv Command line argument values  
 * @return int Exit code (0 for success)
 */
int main(int argc, char *argv[]) {
 
  // 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);
  }
}
 
/**
 * @brief 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
 */
struct ObjectiveFunctionDataImpl : public ObjectiveFunctionData {
  ObjectiveFunctionDataImpl() = default;
  virtual ~ObjectiveFunctionDataImpl() = default;
 
  /// Initialize Python interpreter and load objective function script
  MoFEMErrorCode initPython(const std::string py_file);
 
  /**
   * @brief Evaluate objective function at current 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.
   * 
   * @param coords Gauss point coordinates
   * @param u_ptr Displacement field values
   * @param stress_ptr Stress tensor values  
   * @param strain_ptr Strain tensor values
   * @param o_ptr Output objective function values
   * @return MoFEMErrorCode Success or error code
   */
  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);
 
  /**
   * @brief Compute gradient of objective function with respect to stress
   * 
   * 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.
   * 
   * @param coords Gauss point coordinates
   * @param u_ptr Displacement field values
   * @param stress_ptr Stress tensor values
   * @param strain_ptr Strain tensor values  
   * @param o_ptr Output gradient values
   * @return MoFEMErrorCode Success or error code
   */
  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);
 
  /**
   * @brief Compute gradient of objective function with respect to strain
   * 
   * 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.
   * 
   * @param coords Gauss point coordinates
   * @param u_ptr Displacement field values
   * @param stress_ptr Stress tensor values
   * @param strain_ptr Strain tensor values
   * @param o_ptr Output gradient values  
   * @return MoFEMErrorCode Success or error code
   */
  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);
 
  /**
   * @brief Compute gradient of objective function with respect to displacement
   * 
   * 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.
   * 
   * @param coords Gauss point coordinates  
   * @param u_ptr Displacement field values
   * @param stress_ptr Stress tensor values
   * @param strain_ptr Strain tensor values
   * @param o_ptr Output gradient values
   * @return MoFEMErrorCode Success or error code
   */
  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);
 
  /// Return number of topology optimization modes for given material block
  MoFEMErrorCode numberOfModes(int block_id, int &modes);
 
  /**
   * @brief Define 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.
   * 
   * @param block_id Material block identifier
   * @param coords Element coordinates
   * @param centroid Block centroid coordinates
   * @param bbodx Bounding box dimensions [xmin,xmax,ymin,ymax,zmin,zmax]
   * @param o_ptr Output mode vectors
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords,
                            std::array<double, 3> &centroid,
                            std::array<double, 6> &bbodx, MatrixDouble &o_ptr);
 
private:
  // Python interpreter objects for objective function evaluation
  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
 
  /**
   * @brief 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.
   * 
   * @param coords NumPy array of coordinates
   * @param u NumPy array of displacements
   * @param stress NumPy array of stress tensors
   * @param strain NumPy array of strain tensors  
   * @param o Output NumPy array for objective values
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode objectiveFunctionImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  /**
   * @brief 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.
   * 
   * @param coords NumPy array of coordinates
   * @param u NumPy array of displacements  
   * @param stress NumPy array of stress tensors
   * @param strain NumPy array of strain tensors
   * @param o Output NumPy array for gradient values
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode objectiveGradientStressImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  /**
   * @brief Internal implementation for strain gradient computation
   * 
   * Evaluates ∂f/∂ε through Python interface with NumPy array handling.
   * Provides strain-based sensitivities for comprehensive gradient computation.
   * 
   * @param coords NumPy array of coordinates
   * @param u NumPy array of displacements
   * @param stress NumPy array of stress tensors  
   * @param strain NumPy array of strain tensors
   * @param o Output NumPy array for gradient values
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode objectiveGradientStrainImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  /**
   * @brief 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.
   * 
   * @param coords NumPy array of coordinates
   * @param u NumPy array of displacements
   * @param stress NumPy array of stress tensors
   * @param strain NumPy array of strain tensors
   * @param o Output NumPy array for gradient values  
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode objectiveGradientUImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  /**
   * @brief 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.
   * 
   * @param block_id Material block identifier
   * @param coords NumPy array of element coordinates
   * @param centroid NumPy array of block centroid
   * @param bbodx NumPy array of bounding box dimensions
   * @param o_ptr Output NumPy array for mode vectors
   * @return MoFEMErrorCode Success or error code
   */
  MoFEMErrorCode blockModesImpl(
 
      int block_id, np::ndarray coords, np::ndarray centroid, np::ndarray bbodx,
      np::ndarray &o_ptr
 
  );
 
  /**
   * @brief Convert std::vector to NumPy array for Python interface
   * 
   * Efficient conversion from MoFEM data structures to NumPy arrays
   * for seamless Python function calls without data copying.
   * 
   * @param data Source vector data
   * @param rows Number of rows in resulting array
   * @param nb_gauss_pts Number of Gauss points (affects array structure)
   * @return np::ndarray NumPy array for Python use
   */
  np::ndarray convertToNumPy(std::vector<double> &data, int rows,
                             int nb_gauss_pts);
 
  /**
   * @brief 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.
   * 
   * @param ptr Raw data pointer
   * @param s Array size
   * @return np::ndarray NumPy array for Python use  
   */
  np::ndarray convertToNumPy(double *ptr, int s);
 
  /// Convert symmetric tensor storage to full matrix format
  MatrixDouble copyToFull(MatrixDouble &s);
  
  /// Convert full matrix to symmetric tensor storage format  
  void copyToSymmetric(double *ptr, MatrixDouble &s);
};
 
/**
 * @brief 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
 * 
 * @param py_file Path to Python file containing objective function definitions
 * @return boost::shared_ptr<ObjectiveFunctionData> Configured objective function interface
 * @throws MoFEM exception if Python initialization fails
 */
boost::shared_ptr<ObjectiveFunctionData>
create_python_objective_function(std::string py_file) {
  auto ptr = boost::make_shared<ObjectiveFunctionDataImpl>();
  CHK_THROW_MESSAGE(ptr->initPython(py_file), "init python");
  return ptr;
}
 
/**
 * @brief 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.
 * 
 * @param py_file Path to Python script containing objective function definitions
 * @return MoFEMErrorCode Success or error code
 * @throws MOFEM_OPERATION_UNSUCCESSFUL if Python script has errors
 */
MoFEMErrorCode
ObjectiveFunctionDataImpl::initPython(const std::string py_file) {
  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);
}
 
MatrixDouble ObjectiveFunctionDataImpl::copyToFull(MatrixDouble &s) {
  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;
};
 
void ObjectiveFunctionDataImpl::copyToSymmetric(double *ptr, MatrixDouble &s) {
  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;
  }
}
 
/**
 * @brief Evaluate objective function at current finite element state
 * 
 * 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)
 * 
 * @param coords Gauss point coordinates for current element
 * @param u_ptr Displacement field values at Gauss points
 * @param stress_ptr Cauchy stress tensor values (symmetric storage)
 * @param strain_ptr Strain tensor values (symmetric storage)  
 * @param o_ptr Output objective function values at each Gauss point
 * @return MoFEMErrorCode Success or error code
 */
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) {
  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);
}
 
/**
 * @brief Compute gradient of objective function with respect to stress tensor
 * 
 * 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.
 * 
 * @param coords Gauss point coordinates  
 * @param u_ptr Displacement field values
 * @param stress_ptr Current stress tensor values (symmetric storage)
 * @param strain_ptr Current strain tensor values (symmetric storage)
 * @param o_ptr Output stress gradients ∂f/∂σ (symmetric storage)
 * @return MoFEMErrorCode Success or error code
 */
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) {
  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);
}
 
/**
 * @brief Compute gradient of objective function with respect to strain tensor
 * 
 * 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.
 * 
 * @param coords Gauss point coordinates
 * @param u_ptr Displacement field values
 * @param stress_ptr Current stress tensor values (symmetric storage)
 * @param strain_ptr Current strain tensor values (symmetric storage)  
 * @param o_ptr Output strain gradients ∂f/∂ε (symmetric storage)
 * @return MoFEMErrorCode Success or error code
 */
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) {
  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);
}
 
/**
 * @brief Compute gradient of objective function with respect to displacement field
 * 
 * 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.
 * 
 * @param coords Gauss point coordinates
 * @param u_ptr Displacement field values
 * @param stress_ptr Current stress tensor values  
 * @param strain_ptr Current strain tensor values
 * @param o_ptr Output displacement gradients ∂f/∂u
 * @return MoFEMErrorCode Success or error code
 */
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) {
  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);
}
 
/**
 * @brief Generate spatial topology modes for design optimization
 * 
 * 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.
 * 
 * @param block_id Material block identifier for mode generation
 * @param coords Element coordinates where modes are evaluated
 * @param centroid Geometric centroid of the material block [x,y,z]
 * @param bbodx Bounding box dimensions [xmin,xmax,ymin,ymax,zmin,zmax]
 * @param o_ptr Output matrix: modes × (coordinates × spatial_dimension)
 * @return MoFEMErrorCode Success or error code
 */
MoFEMErrorCode ObjectiveFunctionDataImpl::blockModes(
    int block_id, MatrixDouble &coords, std::array<double, 3> &centroid,
    std::array<double, 6> &bbodx, MatrixDouble &o_ptr) {
  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);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveFunctionImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  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);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveGradientStressImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  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);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveGradientStrainImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  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);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveGradientUImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  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);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::numberOfModes(int block_id,
                                                        int &modes) {
  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);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::blockModesImpl(int block_id,
                                                         np::ndarray coords,
                                                         np::ndarray centroid,
                                                         np::ndarray bbodx,
                                                         np::ndarray &o) {
  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);
}
 
/**
 * @brief Converts a std::vector<double> to a NumPy ndarray.
 *
 * 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.
 *
 * @param data Reference to the vector containing double values to be converted.
 * @param rows Number of rows in the resulting NumPy array.
 * @param nb_gauss_pts Number of Gauss points (columns) in the resulting NumPy
 * array.
 * @return 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.
 */
inline np::ndarray
ObjectiveFunctionDataImpl::convertToNumPy(std::vector<double> &data, int rows,
                                          int nb_gauss_pts) {
  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()));
}
 
inline np::ndarray ObjectiveFunctionDataImpl::convertToNumPy(double *ptr,
                                                             int s) {
  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()));
}
 
Range get_range_from_block(MoFEM::Interface &m_field,
                           const std::string block_name, int dim) {
  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;
};
 
MoFEMErrorCode save_range(moab::Interface &moab, const std::string name,
                          const Range r) {
  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);
};