gradient.cpp

Go to the documentation of this file.

/**
 * @file gradient.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;
 
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();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode setupAdJoint();
  MoFEMErrorCode boundaryCondition();
  MoFEMErrorCode assembleSystem();
 
  // Analysis methods
  MoFEMErrorCode solveElastic(); ///< Solve forward elastic problem
  MoFEMErrorCode postprocessElastic(
      int iter, SmartPetscObj<Vec> gradient_vector = nullptr,
      SmartPetscObj<Vec> adjoint_vector = nullptr,
      SmartPetscObj<Vec> dJ_du = 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, Vec dJ_du);
 
  MoFEMErrorCode testGradient(Vec gradient_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
};
 
//! [Run topology optimization problem]
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
 
  double f;
  auto g = createDMVector(adjointDM, RowColData::ROW);
  auto simple = mField.getInterface<Simple>();
  auto adjoint_vector = createDMVector(simple->getDM());
  auto dJ_du = createDMVector(simple->getDM());
  CHKERR calculateGradient(&f, g, adjoint_vector, dJ_du);
  MOFEM_LOG("WORLD", Sev::inform) << "Objective function value: " << f;
 
  auto grad = createDMVector(simple->getDM());
  CHKERR DMoFEMMeshToLocalVector(adjointDM, g, INSERT_VALUES, SCATTER_REVERSE,
                                 RowColData::ROW);
  CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
      simple->getProblemName(), "U", "ADJOINT_FIELD", RowColData::ROW, grad,
      INSERT_VALUES, SCATTER_FORWARD);
  CHKERR postprocessElastic(0, grad, adjoint_vector,
                            dJ_du); // Post-process initial solution
 
  // Testing gradient
  CHKERR testGradient(g);
 
  // 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]
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");
 
    constexpr int adjoint_field_order = 1;
 
    auto add_field = [&]() {
      MoFEMFunctionBegin;
      CHKERR mField.add_field("ADJOINT_FIELD", H1, base, SPACE_DIM);
      CHKERR mField.add_ents_to_field_by_dim(simple->getBoundaryMeshSet(),
                                             SPACE_DIM - 1, "ADJOINT_FIELD");
      for (auto tt = MBEDGE;
           tt <= moab::CN::TypeDimensionMap[SPACE_DIM - 1].second; ++tt)
        CHKERR mField.set_field_order(simple->getBoundaryMeshSet(), tt,
                                      "ADJOINT_FIELD", adjoint_field_order);
      CHKERR mField.set_field_order(simple->getBoundaryMeshSet(), 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",
                                                        "U");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_DOMAIN_FE",
                                                         "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_DOMAIN_FE",
                                                         "U");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_DOMAIN_FE",
                                                         "GEOMETRY");
 
      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",
                                                        "U");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_BOUNDARY_FE",
                                                         "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_BOUNDARY_FE",
                                                         "U");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_BOUNDARY_FE",
                                                         "GEOMETRY");
 
      CHKERR mField.add_ents_to_finite_element_by_dim(
          simple->getMeshset(), SPACE_DIM, "ADJOINT_DOMAIN_FE");
      CHKERR mField.add_ents_to_finite_element_by_dim(
          simple->getBoundaryMeshSet(), SPACE_DIM - 1, "ADJOINT_BOUNDARY_FE");
      CHKERR mField.build_finite_elements("ADJOINT_DOMAIN_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_FALSE);
      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);
}
//! [Setup adjoint]
 
//! [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");
 
  CHKERR bc_mng->removeBlockDOFsOnEntities("ADJOINT", "REMOVE_X",
                                           "ADJOINT_FIELD", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities("ADJOINT", "REMOVE_Y",
                                           "ADJOINT_FIELD", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities("ADJOINT", "REMOVE_Z",
                                           "ADJOINT_FIELD", 2, 2);
  CHKERR bc_mng->removeBlockDOFsOnEntities("ADJOINT", "REMOVE_ALL",
                                           "ADJOINT_FIELD", 0, 3);
  CHKERR bc_mng->removeBlockDOFsOnEntities("ADJOINT", "FIX_X", "ADJOINT_FIELD",
                                           0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities("ADJOINT", "FIX_Y", "ADJOINT_FIELD",
                                           1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities("ADJOINT", "FIX_Z", "ADJOINT_FIELD",
                                           2, 2);
  CHKERR bc_mng->removeBlockDOFsOnEntities("ADJOINT", "FIX_ALL",
                                           "ADJOINT_FIELD", 0, 3);
 
  MoFEMFunctionReturn(0);
}
//! [Boundary condition]
 
//! [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> gradient_vector,
                                           SmartPetscObj<Vec> adjoint_vector,
                                           SmartPetscObj<Vec> dJ_du) {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  std::vector<std::pair<std::string, SmartPetscObj<Vec>>> additional_vecs;
  if (gradient_vector) {
    additional_vecs.emplace_back("GRADIENT", gradient_vector);
  }
  if (adjoint_vector) {
    additional_vecs.emplace_back("ADJOINT", adjoint_vector);
  }
  if (dJ_du) {
    additional_vecs.emplace_back("dJ_dU", dJ_du);
  }
  CHKERR ElasticPostProc::postProcessElasticResults<SPACE_DIM, DomainEle,
                                                    BoundaryEle, SideEle>(
      mField, simple->getDM(), simple->getDomainFEName(),
      "out_elastic_" + std::to_string(iter) + ".h5m", additional_vecs,
      {}, Sev::noisy);
  MoFEMFunctionReturn(0);
}
//! [Postprocess results]
 
using DomainBaseOp = FormsIntegrators<DomainEleOp>::Assembly<A>::OpBase;
 
 
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;
};
 
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 OpObjective : public ForcesAndSourcesCore::UserDataOperator {
  using OP = ForcesAndSourcesCore::UserDataOperator;
 
  OpObjective(boost::shared_ptr<ObjectiveFunctionData> python_ptr,
              boost::shared_ptr<HookeOps::CommonData> comm_ptr,
              boost::shared_ptr<MatrixDouble> jac_ptr,
              boost::shared_ptr<MatrixDouble> u_ptr,
              boost::shared_ptr<double> glob_objective_ptr)
      : OP(NOSPACE, OP::OPSPACE), pythonPtr(python_ptr), commPtr(comm_ptr),
        jacPtr(jac_ptr), uPtr(u_ptr), globObjectivePtr(glob_objective_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;
 
    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 evaluate_python = [&]() {
      MoFEMFunctionBegin;
      auto &coords = OP::getCoordsAtGaussPts();
      CHKERR pythonPtr->evalInteriorObjectiveFunction(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_ptr);
 
      auto t_obj = getFTensor0FromMat(*objective_ptr);
 
      auto vol = OP::getMeasure();
      auto t_w = getFTensor0IntegrationWeight();
      for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
 
        auto alpha = t_w * vol;
        (*globObjectivePtr) += alpha * t_obj;
 
        ++t_w;
 
        ++t_obj;
      }
      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> uPtr;
  boost::shared_ptr<double> globObjectivePtr;
};
 
struct OpAdJointObjective : public DomainBaseOp {
  using OP = DomainBaseOp;
 
  OpAdJointObjective(boost::shared_ptr<ObjectiveFunctionData> python_ptr,
                     boost::shared_ptr<HookeOps::CommonData> comm_ptr,
                     boost::shared_ptr<MatrixDouble> jac_ptr,
                     boost::shared_ptr<MatrixDouble> u_ptr,
                     boost::shared_ptr<MatrixDouble> grad_lambda_u_ptr,
                     boost::shared_ptr<double> glob_objective_ptr)
      : OP("ADJOINT_FIELD", "ADJOINT_FIELD", OP::OPROW), pythonPtr(python_ptr),
        commPtr(comm_ptr), jacPtr(jac_ptr), uPtr(u_ptr),
        gradLambdaUPtr(grad_lambda_u_ptr),
        globObjectivePtr(glob_objective_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 iNtegrate(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    const auto nb_gauss_pts =
        getGaussPts().size2(); // number of gauss points in the element
    const auto nb_dofs = data.getIndices().size();
    const auto nb_base_funcs = data.getN().size2() / SPACE_DIM;
 
    // 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 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
    // The dJ/du contribution is assembled separately by OpStateSensitivity.
 
    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);
 
      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_grad_lambda_u =
          getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradLambdaUPtr);
 
      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 vol = OP::getMeasure();
      auto t_w = getFTensor0IntegrationWeight();
      auto t_base_diff = data.getFTensor1DiffN<SPACE_DIM>();
      for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
        auto alpha = t_w * vol;
        (*globObjectivePtr) += alpha * t_obj;
 
        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::Tensor4<double, SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM>
            t_diff_inv_jac;
        t_diff_inv_jac(i, j, I, J) =
            -t_inv_jac(i, I) * t_inv_jac(J, j);
        FTensor::Tensor4<double, SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM>
            t_diff_grad;
        t_diff_grad(i, j, I, J) = t_grad_u(i, k) * t_diff_inv_jac(k, j, I, J);
        FTensor::Tensor4<double, SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM>
            t_d_strain;
        t_d_strain(i, j, I, J) =
            t_diff_symm(i, j, k, l) * t_diff_grad(k, l, I, J);
 
        auto t_local_grad_vector = OP::template getNf<SPACE_DIM>();
        int bb = 0;
        for (; bb != nb_dofs/ SPACE_DIM; bb++) {
 
          FTensor::Tensor1<double, SPACE_DIM> t_coef;
          t_coef(I) = (t_inv_jac(J, I) * t_base_diff(J)) * t_det;
 
          t_local_grad_vector(I) +=
              alpha * (
 
                          ((t_obj_dstress(k, l) * t_D(k, l, i, j)) +
                           t_obj_dstrain(i, j)) *
                              (t_d_strain(i, j, I, J) * t_base_diff(J))
 
                          +
 
                          t_obj * t_coef(I)
 
                      );
 
          t_local_grad_vector(I) -=
              alpha * (t_grad_lambda_u(i, j) * t_D(i, j, k, l)) *
              (t_d_strain(k, l, I, J) * t_base_diff(J));
 
          ++t_local_grad_vector;
          ++t_base_diff;
        }
        for (; bb < nb_base_funcs; ++bb) {
          ++t_base_diff;
        }
 
        ++t_w;
        ++t_jac;
 
        ++t_obj;
        ++t_obj_dstress;
        ++t_obj_dstrain;
 
        ++t_grad_u;
        ++t_grad_lambda_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> uPtr;
  boost::shared_ptr<MatrixDouble> gradLambdaUPtr;
 
  boost::shared_ptr<double> globObjectivePtr;
};
 
MoFEMErrorCode Example::calculateGradient(PetscReal *objective_function_value,
                                          Vec objective_function_gradient,
                                          Vec lambda, Vec dJ_du) {
  MOFEM_LOG_CHANNEL("WORLD");
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto dm = simple->getDM();
 
  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_objective_fe = [&](auto lambda_vec, auto glob_objective_ptr,
                              auto fe_rule) {
    auto fe_obj = boost::make_shared<DomainEle>(mField);
    fe_obj->getRuleHook = fe_rule;
    auto &pip = fe_obj->getOpPtrVector();
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
 
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto u_ptr = boost::make_shared<MatrixDouble>();
    auto grad_lambda_ptr = boost::make_shared<MatrixDouble>();
 
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "GEOMETRY", jac_ptr));
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
    auto lambda_smart_vec = SmartPetscObj<Vec>(lambda_vec, true);
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "U", grad_lambda_ptr, lambda_smart_vec));
 
    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, u_ptr,
                                         grad_lambda_ptr,
                                         glob_objective_ptr));
    return fe_obj;
  };
 
  auto evaluate_objective_terms =
      [&](auto objective_fe, auto objective_ptr) {
        MoFEMFunctionBegin;
        *objective_ptr = 0.0;
        CHKERR VecZeroEntries(objective_function_gradient);
        CHKERR VecGhostUpdateBegin(lambda, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(lambda, INSERT_VALUES, SCATTER_FORWARD);
        objective_fe->f = objective_function_gradient;
        CHKERR DMoFEMLoopFiniteElements(adjointDM, "ADJOINT_DOMAIN_FE",
                                        objective_fe);
        MPI_Allreduce(MPI_IN_PLACE, objective_ptr.get(), 1, MPI_DOUBLE, MPI_SUM,
                      mField.get_comm());
        CHKERR VecAssemblyBegin(objective_function_gradient);
        CHKERR VecAssemblyEnd(objective_function_gradient);
        CHKERR VecGhostUpdateBegin(objective_function_gradient, ADD_VALUES,
                                   SCATTER_REVERSE);
        CHKERR VecGhostUpdateEnd(objective_function_gradient, ADD_VALUES,
                                 SCATTER_REVERSE);
        MoFEMFunctionReturn(0);
      };
 
  auto calculate_variance_of_objective_function_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 = [&]() {
    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);
  };
 
  auto fe_rule = [](int, int, int p_data) { return 2 * p_data + p_data - 1; };
  auto objective_ptr_value = boost::make_shared<double>(0.0);
  auto objective_fe = get_objective_fe(lambda, objective_ptr_value, fe_rule);
 
  auto adjoint = [&]() {
    MoFEMFunctionBegin;
 
    CHKERR calculate_variance_of_objective_function_dJ_du();
    CHKERR calculate_adjoint_lambda();
    CHKERR evaluate_objective_terms(objective_fe, objective_ptr_value);
    *objective_function_value = *objective_ptr_value;
 
    MOFEM_LOG("WORLD", Sev::verbose)
        << "Objective function: " << *objective_function_value;
 
    CHKERR VecAssemblyBegin(objective_function_gradient);
    CHKERR VecAssemblyEnd(objective_function_gradient);
    CHKERR VecGhostUpdateBegin(lambda, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(lambda, INSERT_VALUES, SCATTER_FORWARD);
 
    MoFEMFunctionReturn(0);
  };
 
  switch (derivative_type) {
 
  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);
 
  MoFEMFunctionReturn(0);
}
//! [calculateGradient]
 
//! [Finite difference check]
MoFEMErrorCode Example::testGradient(Vec gradient_vector) {
  MoFEMFunctionBegin;
 
  PetscBool gradient_fd_check = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, PETSC_NULLPTR, "-gradient_fd_check",
                             &gradient_fd_check, PETSC_NULLPTR);
  if (gradient_fd_check) {
 
    PetscInt nb_modes = 5;
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR,
                              "-gradient_nb_modes", &nb_modes, PETSC_NULLPTR);
 
    auto simple = mField.getInterface<Simple>();
    auto dm = simple->getDM();
    auto *adj_problem_ptr = getProblemPtr(adjointDM);
 
    auto fe_rule = [](int, int, int p_data) { return 2 * p_data + p_data - 1; };
    auto get_objective_fe = [&](auto glob_objective_ptr, auto fe_rule) {
      auto fe_obj_fe = boost::make_shared<DomainEle>(mField);
      fe_obj_fe->getRuleHook = fe_rule;
      auto &pip = fe_obj_fe->getOpPtrVector();
      AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
      auto jac_ptr = boost::make_shared<MatrixDouble>();
      auto u_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
          "GEOMETRY", jac_ptr));
      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 OpObjective(pythonPtr, common_ptr, jac_ptr, u_ptr,
                                    glob_objective_ptr));
      return fe_obj_fe;
    };
 
    auto f = boost::make_shared<double>(0.);
    auto fe_obj = get_objective_fe(f, fe_rule);
    CacheTupleSharedPtr tmp_cache_ptr = boost::make_shared<CacheTuple>();
    CHKERR mField.cache_problem_entities(simple->getProblemName(),
                                         tmp_cache_ptr);
 
    auto evaluate_objective_terms = [&](auto objective_fe, auto objective_ptr,
                                        double &objective_value) {
      MoFEMFunctionBegin;
      *objective_ptr = 0.0;
      CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                      objective_fe);
      MPI_Allreduce(MPI_IN_PLACE, objective_ptr.get(), 1, MPI_DOUBLE, MPI_SUM,
                    mField.get_comm());
      objective_value = *objective_ptr;
      MoFEMFunctionReturn(0);
    };
 
    auto geometry_bit_number = mField.get_field_bit_number("GEOMETRY");
    for (PetscInt mode = 0; mode < nb_modes; ++mode) {
      MOFEM_LOG_CHANNEL("WORLD");
      MOFEM_LOG("WORLD", Sev::verbose) << "Processing mode: " << mode;
 
      auto &adj_dofs =
          adj_problem_ptr->getNumeredRowDofsPtr()->get<PetscGlobalIdx_mi_tag>();
      auto &dofs = mField.get_dofs()->get<Unique_mi_tag>();
 
      auto a_dof = adj_dofs.find(mode);
      if (a_dof != adj_dofs.end()) {
        auto ent = (*a_dof)->getEnt();
        auto dof_idx = (*a_dof)->getEntDofIdx();
        auto uid = DofEntity::getUniqueIdCalculate(
            dof_idx,
            FieldEntity::getLocalUniqueIdCalculate(geometry_bit_number, ent));
        auto dof = dofs.find(uid);
        if (dof != dofs.end()) {
          auto org_geom_val = (*dof)->getFieldData();
          constexpr double eps = 1e-6;
 
          (*dof)->getFieldData() = org_geom_val + eps;
          CHKERR KSPReset(kspElastic);
          CHKERR solveElastic();
          double f_plus;
          CHKERR evaluate_objective_terms(fe_obj, f, f_plus);
 
          (*dof)->getFieldData() = org_geom_val - eps;
          CHKERR KSPReset(kspElastic);
          CHKERR solveElastic();
          double f_minus;
          CHKERR evaluate_objective_terms(fe_obj, f, f_minus);
 
          double *g_array;
          CHKERR VecGetArray(gradient_vector, &g_array);
          if ((*a_dof)->getHasLocalIndex()) {
            auto adjoint_grad = g_array[(*a_dof)->getPetscLocalDofIdx()];
            auto finite_diff_grad = (f_plus - f_minus) / (2 * eps);
            auto err = std::abs(adjoint_grad - finite_diff_grad) /
                       std::max(std::abs(adjoint_grad),
                                std::abs(finite_diff_grad));
            MOFEM_LOG("WORLD", Sev::inform)
                << "Mode: " << mode << ", Adjoint gradient: " << adjoint_grad
                << ", Finite difference gradient: " << finite_diff_grad
                << ", Relative error: "
                << std::abs(adjoint_grad - finite_diff_grad) /
                       std::max(std::abs(adjoint_grad),
                                std::abs(finite_diff_grad));
            constexpr double tol = 1e-4;
            if (err > tol) {
              SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                      "Gradient check failed for mode %d: relative error %e is "
                      "greater than tolerance %e",
                      mode, err, tol);
            } 
          }
 
          CHKERR VecRestoreArray(gradient_vector, &g_array);
 
          constexpr bool restore_solution = false;
          if constexpr (restore_solution) {
            (*dof)->getFieldData() = org_geom_val;
            CHKERR KSPReset(kspElastic);
            CHKERR solveElastic();
          }
        }
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
//! [Finite difference check]
 
static char help[] = "...\n\n";
 
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);
};