adjoint.cpp

Go to the documentation of this file.

/**
 * @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_shape_optimisation/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;
};
// Here you can see how the template is being used
using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain;
using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle;
using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy;
// Forward declaration
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;
 
enum SensitivityMethod { DIRECT, ADJOINT };
 
SensitivityMethod derivative_type = ADJOINT;
 
#include <ElasticSpring.hpp>
#include <FluidLevel.hpp>
#include <CalculateTraction.hpp>
#include <NaturalDomainBC.hpp>
#include <NaturalBoundaryBC.hpp>
#include <HookeOps.hpp>
#include <ElasticPostProc.hpp>
 
#include <ObjectiveFunctionData.hpp>
using namespace ShapeOptimization;
 
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);
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);
 
  char sensitivity_method_name[32] = "adjoint";
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR,
                               "-sensitivity_method", sensitivity_method_name,
                               sizeof(sensitivity_method_name), PETSC_NULLPTR);
  std::string sensitivity_method = sensitivity_method_name;
  std::transform(sensitivity_method.begin(), sensitivity_method.end(),
                 sensitivity_method.begin(),
                 [](unsigned char c) { return std::tolower(c); });
  if (sensitivity_method == "direct") {
    derivative_type = DIRECT;
  } else if (sensitivity_method == "adjoint") {
    derivative_type = ADJOINT;
  } else {
    SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
            "Unknown -sensitivity_method. Use 'direct' or 'adjoint'.");
  }
  MOFEM_LOG("WORLD", Sev::inform)
      << "Sensitivity method: " << sensitivity_method;
 
  // 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;
      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>();
  CHKERR ElasticPostProc::postProcessElasticResults<SPACE_DIM, DomainEle,
                                                    BoundaryEle, SideEle>(
      mField, simple->getDM(), simple->getDomainFEName(),
      "out_elastic_" + std::to_string(iter) + ".h5m",
      {{"ADJOINT", adjoint_vector}}, {}, Sev::noisy);
  MoFEMFunctionReturn(0);
}
//! [Postprocess results]
 
//! [calculateGradient]
 
using DomainBaseOp = FormsIntegrators<DomainEleOp>::Assembly<A>::OpBase;
 
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);
}
 
//J = 1/2 * sigma : epsilon
// calculate dJ/du = (dJ/dsigma * dsigma/depsilon + dJ/depsilon) depsilon / dgrad_u * dgrad_u / du
struct OpStateSensitivity : public DomainBaseOp {
  using OP = DomainBaseOp;
  OpStateSensitivity(const std::string field_name,
                     boost::shared_ptr<ObjectiveFunctionData> python_ptr,
                     boost::shared_ptr<HookeOps::CommonData> comm_ptr,
                     boost::shared_ptr<MatrixDouble> u_ptr)
      : OP(field_name, field_name, DomainEleOp::OPROW), pythonPtr(python_ptr),
        commPtr(comm_ptr), uPtr(u_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
    FTENSOR_INDEX(SPACE_DIM, i);
    FTENSOR_INDEX(SPACE_DIM, j);
    FTENSOR_INDEX(SPACE_DIM, k);
    FTENSOR_INDEX(SPACE_DIM, l);
 
    constexpr auto symm_size = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
    auto nb_gauss_pts = getGaussPts().size2();
 
    auto objective_dstress =
        boost::make_shared<MatrixDouble>(nb_gauss_pts, symm_size);
    auto objective_dstrain =
        boost::make_shared<MatrixDouble>(nb_gauss_pts, symm_size);
    auto objective_du =
        boost::make_shared<MatrixDouble>(nb_gauss_pts, SPACE_DIM);
 
    auto evaluate_python = [&]() {
      MoFEMFunctionBegin;
      auto &coords = OP::getCoordsAtGaussPts();
      CHKERR pythonPtr->evalInteriorObjectiveGradientStress(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_dstress);
      CHKERR pythonPtr->evalInteriorObjectiveGradientStrain(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_dstrain);
      CHKERR pythonPtr->evalInteriorObjectiveGradientU(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_du);
 
      auto vol = OP::getMeasure();
      auto t_w = OP::getFTensor0IntegrationWeight();
 
      auto t_D =
          getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
      auto t_row_grad = data.getFTensor1DiffN<SPACE_DIM>();
      auto t_row_base = data.getFTensor0N();
 
      auto t_obj_dstress =
          getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstress);
      auto t_obj_dstrain =
          getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstrain);
      auto t_obj_du = getFTensor1FromMat<SPACE_DIM>(*objective_du);
 
      for (int gg = 0; gg != nb_gauss_pts; ++gg) {
        const double alpha = t_w * vol;
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_adjoint_stress;
        t_adjoint_stress(i, j) =
            t_D(i, j, k, l) * t_obj_dstress(k, l) + t_obj_dstrain(i, j);
 
        auto t_nf = OP::template getNf<SPACE_DIM>();
        int rr = 0;
        for (; rr != OP::nbRows / SPACE_DIM; rr++) {
          t_nf(j) += alpha * t_row_grad(i) * t_adjoint_stress(i, j);
          t_nf(j) += alpha * t_row_base * t_obj_du(j);
 
          ++t_row_grad;
          ++t_row_base;
          ++t_nf;
        }
 
        for (; rr < OP::nbRowBaseFunctions; ++rr) {
          ++t_row_grad;
          ++t_row_base;
        }
        ++t_obj_dstrain;
        ++t_obj_dstress;
        ++t_obj_du;
        ++t_w;
      }
      MoFEMFunctionReturn(0);
    };
    CHKERR evaluate_python();
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
  boost::shared_ptr<HookeOps::CommonData> commPtr;
  boost::shared_ptr<MatrixDouble> uPtr;
};
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> d_u_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), dUPtr(d_u_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; // size of symmetric tensor in Voigt notation
 
    auto t_diff_symm = diff_symmetrize(
        FTensor::Number<
            SPACE_DIM>()); // fourth order tensor for symmetrization to convert from gradient to strain
 
    auto nb_gauss_pts =
        getGaussPts().size2(); // number of gauss points in the element
    auto objective_ptr = boost::make_shared<MatrixDouble>(
        1, nb_gauss_pts); // objective function values at gauss points
    auto objective_dstress = boost::make_shared<MatrixDouble>(
        nb_gauss_pts,
        symm_size); // objective function gradient w.r.t. stress at gauss points
    auto objective_dstrain = boost::make_shared<MatrixDouble>(
        nb_gauss_pts,
        symm_size); // objective function gradient w.r.t. strain at gauss points
    auto objective_du =
        boost::make_shared<MatrixDouble>(nb_gauss_pts, SPACE_DIM);
    // We have dJ/dsigma and dJ/depsilon at gauss points, we need to convert it to dJ/du
 
    auto evaluate_python = [&]() {
      MoFEMFunctionBegin;
      auto &coords = OP::getCoordsAtGaussPts();
      CHKERR pythonPtr->evalInteriorObjectiveFunction(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_ptr);
      CHKERR pythonPtr->evalInteriorObjectiveGradientStress(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_dstress);
      CHKERR pythonPtr->evalInteriorObjectiveGradientStrain(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_dstrain);
      CHKERR pythonPtr->evalInteriorObjectiveGradientU(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_du);
 
      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 = getFTensor0FromMat(*objective_ptr);
      auto t_obj_dstress =
          getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstress);
      auto t_obj_dstrain =
          getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstrain);
      auto t_obj_du = getFTensor1FromMat<SPACE_DIM>(*objective_du);
      auto t_d_u = getFTensor1FromMat<SPACE_DIM>(*dUPtr);
 
      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_du(i) * t_d_u(i)
 
                +
 
                t_obj * t_cof
 
            );
 
        ++t_w;
        ++t_jac;
        ++t_diff_jac;
        ++t_cof;
 
        ++t_obj;
        ++t_obj_dstress;
        ++t_obj_dstrain;
        ++t_obj_du;
 
        ++t_grad_u;
        ++t_d_grad;
        ++t_d_u;
      }
      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> dUPtr;
  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
    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) { // is this finite difference!
      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);
  };
  // here starts the preperation for the derivative
  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));
 
    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 d_u_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", d_u_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,
        d_u_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 zero_diff_vec = vectorDuplicate(dm_diff_vec);
  auto zero_state_sensitivity = vectorDuplicate(d);
  CHKERR VecZeroEntries(zero_diff_vec);
  CHKERR VecZeroEntries(zero_state_sensitivity);
  CHKERR VecGhostUpdateBegin(zero_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(zero_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateBegin(zero_state_sensitivity, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(zero_state_sensitivity, INSERT_VALUES,
                           SCATTER_FORWARD);
 
  auto adjoint_fe = get_direvative_fe(dm_diff_vec);
  auto objective_ptr_direct = boost::make_shared<double>(0.0);
  auto objective_grad_ptr_direct = boost::make_shared<double>(0.0);
  auto objective_fe_direct = get_objective_fe(
      dm_diff_vec, d, objective_ptr_direct, objective_grad_ptr_direct);
  auto objective_ptr_explicit = boost::make_shared<double>(0.0);
  auto objective_grad_ptr_explicit = boost::make_shared<double>(0.0);
  auto objective_fe_explicit =
      get_objective_fe(dm_diff_vec, zero_state_sensitivity,
                       objective_ptr_explicit, objective_grad_ptr_explicit);
  auto objective_ptr_value = boost::make_shared<double>(0.0);
  auto objective_grad_ptr_value = boost::make_shared<double>(0.0);
  auto objective_fe_value =
      get_objective_fe(zero_diff_vec, zero_state_sensitivity,
                       objective_ptr_value, objective_grad_ptr_value);
 
  auto set_variance_of_geometry =
      [&](auto mode, auto mod_vec) { // think of mod_vec as X(tau) = X + tau*v_h
        // take the mod_vec and write it into the adjoint DM then copy that
        // field to dm_diff_vec using the same layout as U
        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 = true;
    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 calculate_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 evaluate_objective_terms =
      [&](auto objective_fe, auto objective_ptr, auto objective_grad_ptr,
          double &objective_value, double &objective_gradient) {
        MoFEMFunctionBegin;
        *objective_ptr = 0.0;
        *objective_grad_ptr = 0.0;
        CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                        objective_fe);
        std::array<double, 2> array = {*objective_ptr, *objective_grad_ptr};
        MPI_Allreduce(MPI_IN_PLACE, array.data(), 2, MPI_DOUBLE, MPI_SUM,
                      mField.get_comm());
        objective_value = array[0];
        objective_gradient = array[1];
        MoFEMFunctionReturn(0);
      };
 
  auto calculate_objective_value = [&]() {
    MoFEMFunctionBegin;
    double objective_value = 0;
    double objective_gradient = 0;
    CHKERR evaluate_objective_terms(objective_fe_value, objective_ptr_value,
                                    objective_grad_ptr_value, objective_value,
                                    objective_gradient);
    *objective_function_value = objective_value;
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_variance_of_objective_function_dJ_du = [&](Vec dJ_du) {
    MoFEMFunctionBegin;
 
    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");
 
    auto u_ptr = boost::make_shared<MatrixDouble>();
    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 OpStateSensitivity("U", pythonPtr, common_ptr, u_ptr));
 
    CHKERR VecZeroEntries(dJ_du);
    fe->f = dJ_du;
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), fe);
    CHKERR VecAssemblyBegin(dJ_du);
    CHKERR VecAssemblyEnd(dJ_du);
 
    auto post_proc_rhs = get_essential_fe();
    post_proc_rhs->f = dJ_du;
    CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
                                           post_proc_rhs.get());
 
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_adjoint_lambda = [&](auto lambda, auto dJ_du) {
    MoFEMFunctionBegin;
 
    MOFEM_LOG("WORLD", Sev::inform) << "Solving for adjoint variable lambda";
    CHKERR VecZeroEntries(lambda);
    CHKERR KSPSolveTranspose(kspElastic, dJ_du, lambda);
    CHKERR VecGhostUpdateBegin(lambda, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(lambda, INSERT_VALUES, SCATTER_FORWARD);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR VecZeroEntries(objective_function_gradient);
  CHKERR VecZeroEntries(adjoint_vector);
 
  CHKERR calculate_objective_value();
  MOFEM_LOG("WORLD", Sev::verbose)
      << "Objective function: " << *objective_function_value;
 
  auto direct = [&]() {
    MoFEMFunctionBegin;
    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 calculate_variance_of_displacement(mode, mod_vec);
      double objective_value = 0;
      double objective_gradient = 0;
      CHKERR evaluate_objective_terms(objective_fe_direct, objective_ptr_direct,
                                      objective_grad_ptr_direct,
                                      objective_value, objective_gradient);
      CHKERR VecSetValue(objective_function_gradient, mode, objective_gradient,
                         INSERT_VALUES);
      CHKERR VecAXPY(adjoint_vector, objective_gradient, dm_diff_vec);
      ++mode;
    }
    MoFEMFunctionReturn(0);
  };
 
  auto lambda = vectorDuplicate(f);
  auto dJ_du = vectorDuplicate(f);
 
  auto adjoint = [&]() {
    MoFEMFunctionBegin;
 
    CHKERR calculate_variance_of_objective_function_dJ_du(dJ_du);
    CHKERR calculate_adjoint_lambda(lambda, dJ_du);
 
    int mode = 0;
 
    for (auto mod_vec : modeVecs) {
 
      CHKERR set_variance_of_geometry(mode, mod_vec);
      CHKERR calculate_variance_internal_forces(mode, mod_vec);
      double objective_value = 0;
      double objective_gradient_explicit = 0;
      CHKERR evaluate_objective_terms(
          objective_fe_explicit, objective_ptr_explicit,
          objective_grad_ptr_explicit, objective_value,
          objective_gradient_explicit);
      double lambda_dot_residual_variation = 0;
      CHKERR VecDot(lambda, f, &lambda_dot_residual_variation);
      const double dJ_dp =
          objective_gradient_explicit + lambda_dot_residual_variation;
      CHKERR VecSetValue(objective_function_gradient, mode, dJ_dp,
                         INSERT_VALUES);
      CHKERR VecAXPY(adjoint_vector, dJ_dp, dm_diff_vec);
      ++mode;
    }
    CHKERR VecAssemblyBegin(objective_function_gradient);
    CHKERR VecAssemblyEnd(objective_function_gradient);
    MoFEMFunctionReturn(0);
  };
 
  switch (derivative_type) {
  case DIRECT:
    MOFEM_LOG("WORLD", Sev::inform) << "Running Direct Sensitivity...";
    CHKERR direct();
    break;
 
  case ADJOINT:
    MOFEM_LOG("WORLD", Sev::inform) << "Running Adjoint Sensitivity...";
    CHKERR adjoint();
    break;
 
  default:
    SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
            "Wrong sensitivity type selected");
  }
 
  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);
  }
}
 
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);
};