adjoint.cpp

Go to the documentation of this file.

/**
 * @file elastic.cpp
 * @brief elastic example
 * @author Anonymous author(s) committing under MIT license
 * @example mofem/tutorials/vec-7/adjoint.cpp
 *
 */
 
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
 
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
constexpr int BASE_DIM = 1; //< Dimension of the base functions
 
//! [Define dimension]
constexpr int SPACE_DIM =
    EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
 
//! [Define dimension]
constexpr AssemblyType A = AssemblyType::PETSC;
constexpr IntegrationType I =
    IntegrationType::GAUSS; //< selected integration type
 
//! [Define entities]
using EntData = EntitiesFieldData::EntData;
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
using BoundaryEle =
    PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
//! [Define entities]
 
struct DomainBCs {};
struct BoundaryBCs {};
 
using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<A>::LinearForm<I>;
using OpDomainRhsBCs = DomainRhsBCs::OpFlux<DomainBCs, 1, SPACE_DIM>;
using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::LinearForm<I>;
using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::BiLinearForm<I>;
using OpBoundaryLhsBCs = BoundaryLhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
 
template <int DIM> struct PostProcEleByDim;
 
template <> struct PostProcEleByDim<2> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
};
 
template <> struct PostProcEleByDim<3> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
};
 
using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain;
using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle;
using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy;
 
template <int SPACE_DIM, IntegrationType I, typename OpBase>
struct OpAdJointTestOp;
 
template <int SPACE_DIM, IntegrationType I, typename OpBase>
struct OpAdJointGradTimesSymTensor;
 
constexpr double young_modulus = 1;
constexpr double poisson_ratio = 0.3;
constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
 
PetscBool is_plane_strain = PETSC_FALSE;
 
#include <ElasticSpring.hpp>
#include <FluidLevel.hpp>
#include <CalculateTraction.hpp>
#include <NaturalDomainBC.hpp>
#include <NaturalBoundaryBC.hpp>
#include <HookeOps.hpp>
 
struct ObjectiveFunctionData {
 
  virtual MoFEMErrorCode
  evalObjectiveFunction(MatrixDouble &coords,
                        boost::shared_ptr<MatrixDouble> u_ptr,
                        boost::shared_ptr<MatrixDouble> stress_ptr,
                        boost::shared_ptr<MatrixDouble> strain_ptr,
                        boost::shared_ptr<VectorDouble> o_ptr) = 0;
 
  virtual MoFEMErrorCode
  evalObjectiveGradientStress(MatrixDouble &coords,
                              boost::shared_ptr<MatrixDouble> u_ptr,
                              boost::shared_ptr<MatrixDouble> stress_ptr,
                              boost::shared_ptr<MatrixDouble> strain_ptr,
                              boost::shared_ptr<MatrixDouble> o_ptr) = 0;
 
  virtual MoFEMErrorCode
  evalObjectiveGradientStrain(MatrixDouble &coords,
                              boost::shared_ptr<MatrixDouble> u_ptr,
                              boost::shared_ptr<MatrixDouble> stress_ptr,
                              boost::shared_ptr<MatrixDouble> strain_ptr,
                              boost::shared_ptr<MatrixDouble> o_ptr) = 0;
 
  virtual MoFEMErrorCode numberOfModes(int block_id, int &modes) = 0;
 
  virtual MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords,
                                    std::array<double, 3> &centroid,
                                    std::array<double, 6> &bbox,
                                    MatrixDouble &o_ptr) = 0;
 
  virtual ~ObjectiveFunctionData() = default;
};
 
boost::shared_ptr<ObjectiveFunctionData>
create_python_objective_function(std::string py_file);
 
Range get_range_from_block(MoFEM::Interface &m_field,
                           const std::string block_name, int dim);
 
MoFEMErrorCode save_range(moab::Interface &moab, const std::string name,
                          const Range r);
struct Example {
 
  Example(MoFEM::Interface &m_field) : mField(m_field) {}
 
  MoFEMErrorCode runProblem();
 
private:
  MoFEM::Interface &mField;
 
  boost::shared_ptr<MatrixDouble> vectorFieldPtr = nullptr;
 
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode setupAdJoint();
  MoFEMErrorCode boundaryCondition();
  MoFEMErrorCode topologyModes();
  MoFEMErrorCode assembleSystem();
  MoFEMErrorCode solveElastic();
  MoFEMErrorCode
  postprocessElastic(int iter, SmartPetscObj<Vec> adjoint_vector = nullptr);
  MoFEMErrorCode calculateGradient(PetscReal *objective_function_value,
                                   Vec objective_function_gradient,
                                   Vec adjoint_vector);
 
  FieldApproximationBase base;
  int fieldOrder = 2;
 
  SmartPetscObj<KSP> kspElastic;
  SmartPetscObj<DM> adjointDM;
  boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
  std::vector<SmartPetscObj<Vec>> modeVecs;
  std::vector<std::array<double, 3>> modeCentroids;
  std::vector<std::array<double, 6>> modeBBoxes;
  SmartPetscObj<Vec> initialGeometry;
};
 
//! [Run problem]
MoFEMErrorCode Example::runProblem() {
  MoFEMFunctionBegin;
 
  char objective_function_file_name[255] = "objective_function.py";
  CHKERR PetscOptionsGetString(
      PETSC_NULLPTR, PETSC_NULLPTR, "-objective_function",
      objective_function_file_name, 255, PETSC_NULLPTR);
  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");
  }
 
  pythonPtr = create_python_objective_function(objective_function_file_name);
 
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR setupAdJoint();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
 
  auto pip = mField.getInterface<PipelineManager>();
  kspElastic = pip->createKSP();
  CHKERR KSPSetFromOptions(kspElastic);
 
  CHKERR solveElastic();
  CHKERR postprocessElastic(-1);
 
  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);
  };
 
  CHKERR create_vec_modes("OPTIMISE");
  CHKERR get_modes_bounding_boxes("OPTIMISE");
  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);
 
  initialGeometry = createDMVector(adjointDM);
  CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
      "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW, initialGeometry,
      INSERT_VALUES, SCATTER_FORWARD);
 
  CHKERR topologyModes();
  auto x0 = vectorDuplicate(g);
 
  CHKERR VecSet(x0, 0.0);
  CHKERR TaoSetSolution(tao, x0);
  CHKERR TaoSetFromOptions(tao);
  CHKERR TaoSolve(tao);
 
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Read mesh]
MoFEMErrorCode Example::readMesh() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
  // Add meshsets if config file provided
  CHKERR mField.getInterface<MeshsetsManager>()->setMeshsetFromFile();
  MoFEMFunctionReturn(0);
}
//! [Read mesh]
 
//! [Set up problem]
MoFEMErrorCode Example::setupProblem() {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  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 field
  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &fieldOrder,
                            PETSC_NULLPTR);
 
  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
 
  CHKERR simple->setFieldOrder("U", fieldOrder);
  CHKERR simple->setFieldOrder("GEOMETRY", fieldOrder);
  CHKERR simple->setUp();
 
  auto project_ho_geometry = [&]() {
    Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
    return mField.loop_dofs("GEOMETRY", ent_method);
  };
  CHKERR project_ho_geometry();
 
  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>();
 
  // ad adjoint dm
 
  auto create_adjoint_dm = [&]() {
    auto adjoint_dm = createDM(mField.get_comm(), "DMMOFEM");
 
    auto add_field = [&]() {
      MoFEMFunctionBegin;
      CHKERR mField.add_field("ADJOINT_FIELD", H1, base, SPACE_DIM);
      CHKERR mField.add_ents_to_field_by_dim(simple->getMeshset(), SPACE_DIM,
                                             "ADJOINT_FIELD");
      for (auto tt = MBEDGE; tt <= moab::CN::TypeDimensionMap[SPACE_DIM].second;
           ++tt)
        CHKERR mField.set_field_order(simple->getMeshset(), tt, "ADJOINT_FIELD",
                                      fieldOrder);
      CHKERR mField.set_field_order(simple->getMeshset(), MBVERTEX,
                                    "ADJOINT_FIELD", 1);
      CHKERR mField.build_fields();
      MoFEMFunctionReturn(0);
    };
 
    auto add_adjoint_fe_impl = [&]() {
      MoFEMFunctionBegin;
      CHKERR mField.add_finite_element("ADJOINT_DOMAIN_FE");
      CHKERR mField.modify_finite_element_add_field_row("ADJOINT_DOMAIN_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_col("ADJOINT_DOMAIN_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_DOMAIN_FE",
                                                         "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_DOMAIN_FE",
                                                         "GEOMETRY");
      CHKERR mField.add_ents_to_finite_element_by_dim(
          simple->getMeshset(), SPACE_DIM, "ADJOINT_DOMAIN_FE");
      CHKERR mField.build_finite_elements("ADJOINT_DOMAIN_FE");
 
      CHKERR mField.add_finite_element("ADJOINT_BOUNDARY_FE");
      CHKERR mField.modify_finite_element_add_field_row("ADJOINT_BOUNDARY_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_col("ADJOINT_BOUNDARY_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_BOUNDARY_FE",
                                                         "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_BOUNDARY_FE",
                                                         "GEOMETRY");
 
      auto block_name = "OPTIMISE";
      auto mesh_mng = mField.getInterface<MeshsetsManager>();
      auto bcs = mesh_mng->getCubitMeshsetPtr(
 
          std::regex((boost::format("%s(.*)") % block_name).str())
 
      );
 
      for (auto bc : bcs) {
        CHKERR mField.add_ents_to_finite_element_by_dim(
            bc->getMeshset(), SPACE_DIM - 1, "ADJOINT_BOUNDARY_FE");
      }
 
      CHKERR mField.build_finite_elements("ADJOINT_BOUNDARY_FE");
 
      CHKERR mField.build_adjacencies(simple->getBitRefLevel(),
                                      simple->getBitRefLevelMask());
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_adjoint_dm_imp = [&]() {
      MoFEMFunctionBegin;
      CHKERR DMMoFEMCreateMoFEM(adjoint_dm, &mField, "ADJOINT",
                                simple->getBitRefLevel(),
                                simple->getBitRefLevelMask());
      CHKERR DMMoFEMSetDestroyProblem(adjoint_dm, PETSC_TRUE);
      CHKERR DMSetFromOptions(adjoint_dm);
      CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_DOMAIN_FE");
      CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_BOUNDARY_FE");
      CHKERR DMMoFEMSetSquareProblem(adjoint_dm, PETSC_TRUE);
      CHKERR DMMoFEMSetIsPartitioned(adjoint_dm, PETSC_TRUE);
      mField.getInterface<ProblemsManager>()->buildProblemFromFields =
          PETSC_TRUE;
      CHKERR DMSetUp(adjoint_dm);
      mField.getInterface<ProblemsManager>()->buildProblemFromFields =
          PETSC_FALSE;
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(add_field(), "add adjoint field");
    CHK_THROW_MESSAGE(add_adjoint_fe_impl(), "add adjoint fe");
    CHK_THROW_MESSAGE(set_adjoint_dm_imp(), "set adjoint dm");
 
    return adjoint_dm;
  };
 
  adjointDM = create_adjoint_dm();
 
  MoFEMFunctionReturn(0);
}
//
 
//! [Boundary condition]
MoFEMErrorCode Example::boundaryCondition() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
 
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
                                           "U", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
                                           "U", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
                                           "U", 2, 2);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                           "REMOVE_ALL", "U", 0, 3);
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "U");
 
  MoFEMFunctionReturn(0);
}
//! [Boundary condition]
 
//! [Adjoint modes]
MoFEMErrorCode Example::topologyModes() {
  MoFEMFunctionBegin;
 
  auto opt_ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
  auto subset_dm_bdy = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMSetSquareProblem(subset_dm_bdy, PETSC_TRUE);
  CHKERR DMMoFEMCreateSubDM(subset_dm_bdy, adjointDM, "SUBSET_BDY");
  CHKERR DMMoFEMAddElement(subset_dm_bdy, "ADJOINT_BOUNDARY_FE");
  CHKERR DMMoFEMAddSubFieldRow(subset_dm_bdy, "ADJOINT_FIELD",
                               boost::make_shared<Range>(opt_ents));
  CHKERR DMMoFEMAddSubFieldCol(subset_dm_bdy, "ADJOINT_FIELD",
                               boost::make_shared<Range>(opt_ents));
  CHKERR DMSetUp(subset_dm_bdy);
 
  auto subset_dm_domain = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMSetSquareProblem(subset_dm_domain, PETSC_TRUE);
  CHKERR DMMoFEMCreateSubDM(subset_dm_domain, adjointDM, "SUBSET_DOMAIN");
  CHKERR DMMoFEMAddElement(subset_dm_domain, "ADJOINT_DOMAIN_FE");
  CHKERR DMMoFEMAddSubFieldRow(subset_dm_domain, "ADJOINT_FIELD");
  CHKERR DMMoFEMAddSubFieldCol(subset_dm_domain, "ADJOINT_FIELD");
  CHKERR DMSetUp(subset_dm_domain);
 
  // remove dofs on boundary of the domain
  auto remove_dofs = [&]() {
    MoFEMFunctionBegin;
 
    std::array<Range, 3> remove_dim_ents;
    remove_dim_ents[0] =
        get_range_from_block(mField, "OPT_REMOVE_X", SPACE_DIM - 1);
    remove_dim_ents[1] =
        get_range_from_block(mField, "OPT_REMOVE_Y", SPACE_DIM - 1);
    remove_dim_ents[2] =
        get_range_from_block(mField, "OPT_REMOVE_Z", SPACE_DIM - 1);
 
    for (int d = 0; d != 3; ++d) {
      MOFEM_LOG("WORLD", Sev::inform)
          << "Removing topology modes on block OPT_REMOVE_" << (char)('X' + d)
          << " with " << remove_dim_ents[d].size() << " entities";
    }
 
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents,
                                                       true);
    auto skin = moab::Skinner(&mField.get_moab());
    Range boundary_ents;
    CHKERR skin.find_skin(0, body_ents, false, boundary_ents);
    for (int d = 0; d != 3; ++d) {
      boundary_ents = subtract(boundary_ents, remove_dim_ents[d]);
    }
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(boundary_ents,
                                 PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, nullptr);
    for (auto d = SPACE_DIM - 2; d >= 0; --d) {
      if (d >= 0) {
        Range ents;
        CHKERR mField.get_moab().get_adjacencies(boundary_ents, d, false, ents,
                                                 moab::Interface::UNION);
        boundary_ents.merge(ents);
      } else {
        Range verts;
        CHKERR mField.get_moab().get_connectivity(boundary_ents, verts);
        boundary_ents.merge(verts);
      }
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          boundary_ents);
    }
    boundary_ents.merge(opt_ents);
    CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
        "SUBSET_DOMAIN", "ADJOINT_FIELD", boundary_ents);
    for (int d = 0; d != 3; ++d) {
      CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
          "SUBSET_DOMAIN", "ADJOINT_FIELD", remove_dim_ents[d], d, d);
    }
 
    // #ifndef NDEBUG
    if (mField.get_comm_rank() == 0) {
      CHKERR save_range(mField.get_moab(), "topoMode_boundary_ents.vtk",
                        boundary_ents);
    }
    // #endif
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR remove_dofs();
 
  auto get_lhs_fe = [&]() {
    auto fe_lhs = boost::make_shared<BoundaryEle>(mField);
    fe_lhs->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe_lhs->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
                                                              "GEOMETRY");
    using OpMass = FormsIntegrators<BoundaryEleOp>::Assembly<A>::BiLinearForm<
        I>::OpMass<1, SPACE_DIM>;
    pip.push_back(new OpMass("ADJOINT_FIELD", "ADJOINT_FIELD",
                             [](double, double, double) { return 1.; }));
    return fe_lhs;
  };
 
  auto get_rhs_fe = [&]() {
    auto fe_rhs = boost::make_shared<BoundaryEle>(mField);
    fe_rhs->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe_rhs->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
                                                              "GEOMETRY");
 
    return fe_rhs;
  };
 
  auto block_name = "OPTIMISE";
  auto mesh_mng = mField.getInterface<MeshsetsManager>();
  auto bcs = mesh_mng->getCubitMeshsetPtr(
 
      std::regex((boost::format("%s(.*)") % block_name).str())
 
  );
 
  for (auto &v : modeVecs) {
    v = createDMVector(subset_dm_bdy);
  }
 
  using OP = FormsIntegrators<BoundaryEleOp>::Assembly<A>::OpBase;
  struct OpMode : public OP {
    OpMode(const std::string name,
           boost::shared_ptr<ObjectiveFunctionData> python_ptr, int id,
           std::vector<SmartPetscObj<Vec>> mode_vecs,
           std::vector<std::array<double, 3>> mode_centroids,
           std::vector<std::array<double, 6>> mode_bboxes, int block_counter,
           int mode_counter, boost::shared_ptr<Range> range = nullptr)
        : OP(name, name, OP::OPROW, range), pythonPtr(python_ptr), iD(id),
          modeVecs(mode_vecs), modeCentroids(mode_centroids),
          modeBboxes(mode_bboxes), blockCounter(block_counter),
          modeCounter(mode_counter) {}
 
    MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
      MoFEMFunctionBegin;
 
      if (OP::entsPtr) {
        if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
          MoFEMFunctionReturnHot(0);
      }
 
      auto nb_rows = data.getIndices().size();
      if (!nb_rows) {
        MoFEMFunctionReturnHot(0);
      }
      auto nb_base_functions = data.getN().size2();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      CHKERR pythonPtr->blockModes(iD, OP::getCoordsAtGaussPts(),
                                   modeCentroids[blockCounter],
                                   modeBboxes[blockCounter], blockModes);
 
      auto nb_integration_pts = getGaussPts().size2();
      if (blockModes.size2() != 3 * nb_integration_pts) {
        MOFEM_LOG("WORLD", Sev::error)
            << "Number of modes does not match number of integration points: "
            << blockModes.size2() << "!=" << 3 * nb_integration_pts;
        CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "modes/integration points");
      }
 
      VectorDouble nf(nb_rows);
 
      int nb_modes = blockModes.size1();
      for (auto mode = 0; mode != nb_modes; ++mode) {
        nf.clear();
        // get mode
        auto t_mode = getFTensor1FromPtr<3>(&blockModes(mode, 0));
        // get element volume
        const double vol = OP::getMeasure();
        // get integration weights
        auto t_w = OP::getFTensor0IntegrationWeight();
        // get base function gradient on rows
        auto t_base = data.getFTensor0N();
        // loop over integration points
        for (int gg = 0; gg != nb_integration_pts; gg++) {
 
          // take into account Jacobian
          const double alpha = t_w * vol;
          // loop over rows base functions
          auto t_nf = getFTensor1FromPtr<SPACE_DIM>(nf.data().data());
          int rr = 0;
          for (; rr != nb_rows / SPACE_DIM; ++rr) {
            t_nf(i) += alpha * t_base * t_mode(i);
            ++t_base;
            ++t_nf;
          }
          for (; rr < nb_base_functions; ++rr)
            ++t_base;
          ++t_w;    // move to another integration weight
          ++t_mode; // move to another mode
        }
        Vec vec = modeVecs[modeCounter + mode];
        auto size = data.getIndices().size();
        auto *indices = data.getIndices().data().data();
        auto *nf_data = nf.data().data();
        CHKERR VecSetValues(vec, size, indices, nf_data, ADD_VALUES);
      }
 
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
    MatrixDouble blockModes;
    std::vector<std::array<double, 3>> modeCentroids;
    std::vector<std::array<double, 6>> modeBboxes;
    int iD;
    std::vector<SmartPetscObj<Vec>> modeVecs;
    int blockCounter;
    int modeCounter;
  };
 
  auto solve_bdy = [&]() {
    MoFEMFunctionBegin;
 
    auto fe_lhs = get_lhs_fe();
    auto fe_rhs = get_rhs_fe();
    int block_counter = 0;
    int mode_counter = 0;
    for (auto &bc : bcs) {
      auto id = bc->getMeshsetId();
      Range ents;
      CHKERR mField.get_moab().get_entities_by_handle(bc->getMeshset(), ents,
                                                      true);
      auto range = boost::make_shared<Range>(ents);
      auto &pip_rhs = fe_rhs->getOpPtrVector();
      pip_rhs.push_back(new OpMode("ADJOINT_FIELD", pythonPtr, id, modeVecs,
                                   modeCentroids, modeBBoxes, block_counter,
                                   mode_counter, range));
      CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
                                      fe_rhs);
      pip_rhs.pop_back();
      int nb_modes;
      CHKERR pythonPtr->numberOfModes(id, nb_modes);
      ++block_counter;
      mode_counter += nb_modes;
      MOFEM_LOG("WORLD", Sev::inform)
          << "Setting mode block block: " << bc->getName()
          << " with ID: " << bc->getMeshsetId()
          << " total modes: " << mode_counter;
    }
 
    for (auto &v : modeVecs) {
      CHKERR VecAssemblyBegin(v);
      CHKERR VecAssemblyEnd(v);
      CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
    }
 
    auto M = createDMMatrix(subset_dm_bdy);
    fe_lhs->B = M;
    CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
                                    fe_lhs);
    CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
 
    auto solver = createKSP(mField.get_comm());
    CHKERR KSPSetOperators(solver, M, M);
    CHKERR KSPSetFromOptions(solver);
    CHKERR KSPSetUp(solver);
    auto v = createDMVector(subset_dm_bdy);
    for (auto &f : modeVecs) {
      CHKERR KSPSolve(solver, f, v);
      CHKERR VecSwap(f, v);
    }
 
    for (auto &v : modeVecs) {
      CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR solve_bdy();
 
  auto get_elastic_fe_lhs = [&]() {
    auto fe = boost::make_shared<DomainEle>(mField);
    fe->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
    CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
        mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
    return fe;
  };
 
  auto get_elastic_fe_rhs = [&]() {
    auto fe = boost::make_shared<DomainEle>(mField);
    fe->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
    CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
        mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
    return fe;
  };
 
  auto adjoint_gradient_postprocess = [&](auto mode) {
    MoFEMFunctionBegin;
    auto post_proc_mesh = boost::make_shared<moab::Core>();
    auto post_proc_begin =
        boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(mField,
                                                              post_proc_mesh);
    auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
        mField, post_proc_mesh);
 
    auto geom_vec = boost::make_shared<MatrixDouble>();
 
    auto post_proc_fe =
        boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        post_proc_fe->getOpPtrVector(), {H1}, "GEOMETRY");
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("ADJOINT_FIELD", geom_vec,
                                                    modeVecs[mode]));
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
            {},
 
            {{"MODE", geom_vec}},
 
            {},
 
            {}
 
            )
 
    );
 
    CHKERR DMoFEMPreProcessFiniteElements(adjointDM,
                                          post_proc_begin->getFEMethod());
    CHKERR DMoFEMLoopFiniteElements(adjointDM, "ADJOINT_DOMAIN_FE",
                                    post_proc_fe);
    CHKERR DMoFEMPostProcessFiniteElements(adjointDM,
                                           post_proc_begin->getFEMethod());
 
    CHKERR post_proc_end->writeFile("mode_" + std::to_string(mode) + ".h5m");
 
    MoFEMFunctionReturn(0);
  };
 
  auto solve_domain = [&]() {
    MoFEMFunctionBegin;
    auto fe_lhs = get_elastic_fe_lhs();
    auto fe_rhs = get_elastic_fe_rhs();
    auto v = createDMVector(subset_dm_domain);
    auto F = vectorDuplicate(v);
    fe_rhs->f = F;
 
    auto M = createDMMatrix(subset_dm_domain);
    fe_lhs->B = M;
    CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
                                    fe_lhs);
    CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
 
    auto solver = createKSP(mField.get_comm());
    CHKERR KSPSetOperators(solver, M, M);
    CHKERR KSPSetFromOptions(solver);
    CHKERR KSPSetUp(solver);
 
    int mode_counter = 0;
    for (auto &f : modeVecs) {
      CHKERR mField.getInterface<FieldBlas>()->setField(0, "ADJOINT_FIELD");
      CHKERR DMoFEMMeshToLocalVector(subset_dm_bdy, f, INSERT_VALUES,
                                     SCATTER_REVERSE);
      CHKERR VecZeroEntries(F);
      CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
                                      fe_rhs);
      CHKERR VecAssemblyBegin(F);
      CHKERR VecAssemblyEnd(F);
      CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
      CHKERR KSPSolve(solver, F, v);
      CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR DMoFEMMeshToLocalVector(subset_dm_domain, v, INSERT_VALUES,
                                     SCATTER_REVERSE);
      auto m = createDMVector(adjointDM);
      CHKERR DMoFEMMeshToLocalVector(adjointDM, m, INSERT_VALUES,
                                     SCATTER_FORWARD);
      f = m;
      ++mode_counter;
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR solve_domain();
 
  for (int i = 0; i < modeVecs.size(); ++i) {
    CHKERR adjoint_gradient_postprocess(i);
  }
 
  MoFEMFunctionReturn(0);
}
//! [Adjoint modes]
 
//! [Push operators to pipeline]
MoFEMErrorCode Example::assembleSystem() {
  MoFEMFunctionBegin;
  auto pip = mField.getInterface<PipelineManager>();
 
  //! [Integration rule]
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
 
  auto integration_rule_bc = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
 
  CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
  CHKERR pip->setDomainLhsIntegrationRule(integration_rule);
  CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
  CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
  //! [Integration rule]
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
 
  //! [Push domain stiffness matrix to pipeline]
  // Add LHS operator for elasticity (stiffness matrix)
  CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
      mField, pip->getOpDomainLhsPipeline(), "U", "MAT_ELASTIC", Sev::noisy);
  //! [Push domain stiffness matrix to pipeline]
 
  // Add RHS operator for internal forces
  CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
      mField, pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC", Sev::noisy);
 
  //! [Push Body forces]
  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
      pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
  //! [Push Body forces]
 
  //! [Push natural boundary conditions]
  // Add force boundary condition
  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
      pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
  // Add case for mix type of BCs
  CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
      pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::noisy);
  //! [Push natural boundary conditions]
  MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
 
//! [Solve]
MoFEMErrorCode Example::solveElastic() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto dm = simple->getDM();
  auto f = createDMVector(dm);
  auto d = vectorDuplicate(f);
  CHKERR VecZeroEntries(d);
  CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
 
  auto set_essential_bc = [&]() {
    MoFEMFunctionBegin;
    // This is low level pushing finite elements (pipelines) to solver
 
    auto ksp_ctx_ptr = getDMKspCtx(dm);
    auto pre_proc_rhs = boost::make_shared<FEMethod>();
    auto post_proc_rhs = boost::make_shared<FEMethod>();
    auto post_proc_lhs = boost::make_shared<FEMethod>();
 
    auto get_pre_proc_hook = [&]() {
      return EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_rhs,
                                                       {});
    };
    pre_proc_rhs->preProcessHook = get_pre_proc_hook();
 
    auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
                                                           post_proc_rhs, 1.)();
      MoFEMFunctionReturnHot(0);
    };
 
    auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(mField,
                                                           post_proc_lhs, 1.)();
      MoFEMFunctionReturnHot(0);
    };
 
    post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
    post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
 
    ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
    ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
    ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
    MoFEMFunctionReturn(0);
  };
 
  auto setup_and_solve = [&](auto solver) {
    MoFEMFunctionBegin;
    BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp";
    CHKERR KSPSetUp(solver);
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp <= Done";
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve";
    CHKERR KSPSolve(solver, f, d);
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve <= Done";
    MoFEMFunctionReturn(0);
  };
 
  MOFEM_LOG_CHANNEL("TIMER");
  MOFEM_LOG_TAG("TIMER", "timer");
 
  CHKERR set_essential_bc();
  CHKERR setup_and_solve(kspElastic);
 
  CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
 
  auto evaluate_field_at_the_point = [&]() {
    MoFEMFunctionBegin;
 
    int coords_dim = 3;
    std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
    PetscBool do_eval_field = PETSC_FALSE;
    CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
                                    field_eval_coords.data(), &coords_dim,
                                    &do_eval_field);
 
    if (do_eval_field) {
 
      vectorFieldPtr = boost::make_shared<MatrixDouble>();
      auto field_eval_data =
          mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
 
      CHKERR mField.getInterface<FieldEvaluatorInterface>()
          ->buildTree<SPACE_DIM>(field_eval_data, simple->getDomainFEName());
 
      field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
      auto no_rule = [](int, int, int) { return -1; };
      auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
      field_eval_fe_ptr->getRuleHook = no_rule;
 
      field_eval_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", vectorFieldPtr));
 
      CHKERR mField.getInterface<FieldEvaluatorInterface>()
          ->evalFEAtThePoint<SPACE_DIM>(
              field_eval_coords.data(), 1e-12, simple->getProblemName(),
              simple->getDomainFEName(), field_eval_data,
              mField.get_comm_rank(), mField.get_comm_rank(), nullptr, MF_EXIST,
              QUIET);
 
      if (vectorFieldPtr->size1()) {
        auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
        if constexpr (SPACE_DIM == 2)
          MOFEM_LOG("FieldEvaluator", Sev::inform)
              << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
        else
          MOFEM_LOG("FieldEvaluator", Sev::inform)
              << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
              << " U_Z: " << t_disp(2);
      }
 
      MOFEM_LOG_SYNCHRONISE(mField.get_comm());
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR evaluate_field_at_the_point();
 
  MoFEMFunctionReturn(0);
}
//! [Solve]
 
//! [Postprocess results]
MoFEMErrorCode Example::postprocessElastic(int iter,
                                           SmartPetscObj<Vec> adjoint_vector) {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>();
  auto det_ptr = boost::make_shared<VectorDouble>();
  auto jac_ptr = boost::make_shared<MatrixDouble>();
  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
  //! [Postprocess clean]
  pip->getDomainRhsFE().reset();
  pip->getDomainLhsFE().reset();
  pip->getBoundaryRhsFE().reset();
  pip->getBoundaryLhsFE().reset();
  //! [Postprocess clean]
 
  //! [Postprocess initialise]
  auto post_proc_mesh = boost::make_shared<moab::Core>();
  auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
      mField, post_proc_mesh);
  auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
      mField, post_proc_mesh);
  //! [Postprocess initialise]
 
  // lambada function to calculate displacement, strain and stress
  auto calculate_stress_ops = [&](auto &pip) {
    // Add H1 geometry ops
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
 
    // Use HookeOps commonDataFactory to get matStrainPtr and matCauchyStressPtr
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEle>(
        mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
 
    // Store U and GEOMETRY values if needed
    auto u_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
    auto x_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
 
    using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
 
    DataMapMat vec_map = {{"U", u_ptr}, {"GEOMETRY", x_ptr}};
    DataMapMat sym_map = {{"STRAIN", common_ptr->getMatStrain()},
                          {"STRESS", common_ptr->getMatCauchyStress()}};
 
    if (adjoint_vector) {
      auto adjoint_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
          "U", adjoint_ptr, adjoint_vector));
      vec_map["ADJOINT"] = adjoint_ptr;
    }
 
    // Return what you need: displacements, coordinates, strain, stress
    return boost::make_tuple(u_ptr, x_ptr, common_ptr->getMatStrain(),
                             common_ptr->getMatCauchyStress(), vec_map,
                             sym_map);
  };
 
  auto get_tag_id_on_pmesh = [&](bool post_proc_skin) {
    int def_val_int = 0;
    Tag tag_mat;
    CHKERR mField.get_moab().tag_get_handle(
        "MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
        MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
    auto meshset_vec_ptr =
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
            std::regex((boost::format("%s(.*)") % "MAT_ELASTIC").str()));
 
    Range skin_ents;
    std::unique_ptr<Skinner> skin_ptr;
    if (post_proc_skin) {
      skin_ptr = std::make_unique<Skinner>(&mField.get_moab());
      auto boundary_meshset = simple->getBoundaryMeshSet();
      CHKERR mField.get_moab().get_entities_by_handle(boundary_meshset,
                                                      skin_ents, true);
    }
 
    for (auto m : meshset_vec_ptr) {
      Range ents_3d;
      CHKERR mField.get_moab().get_entities_by_handle(m->getMeshset(), ents_3d,
                                                      true);
      int const id = m->getMeshsetId();
      ents_3d = ents_3d.subset_by_dimension(SPACE_DIM);
      if (post_proc_skin) {
        Range skin_faces;
        CHKERR skin_ptr->find_skin(0, ents_3d, false, skin_faces);
        ents_3d = intersect(skin_ents, skin_faces);
      }
      CHKERR mField.get_moab().tag_clear_data(tag_mat, ents_3d, &id);
    }
 
    return tag_mat;
  };
 
  auto post_proc_domain = [&](auto post_proc_mesh) {
    auto post_proc_fe =
        boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
        calculate_stress_ops(post_proc_fe->getOpPtrVector());
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
            {},
 
            vec_map,
 
            {},
 
            sym_map
 
            )
 
    );
 
    post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(false)});
    return post_proc_fe;
  };
 
  auto post_proc_boundary = [&](auto post_proc_mesh) {
    auto post_proc_fe =
        boost::make_shared<PostProcEleBdy>(mField, post_proc_mesh);
    AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
    auto op_loop_side =
        new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
    // push ops to side element, through op_loop_side operator
    auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
        calculate_stress_ops(op_loop_side->getOpPtrVector());
    post_proc_fe->getOpPtrVector().push_back(op_loop_side);
    auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
    post_proc_fe->getOpPtrVector().push_back(
        new ElasticExample::OpCalculateTraction(mat_stress_ptr,
                                                mat_traction_ptr));
    vec_map["T"] = mat_traction_ptr;
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
            {},
 
            vec_map,
 
            {},
 
            sym_map
 
            )
 
    );
 
    post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(true)});
    return post_proc_fe;
  };
 
  PetscBool post_proc_skin_only = PETSC_FALSE;
  if (SPACE_DIM == 3) {
    post_proc_skin_only = PETSC_TRUE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin_only",
                               &post_proc_skin_only, PETSC_NULLPTR);
  }
  if (post_proc_skin_only == PETSC_FALSE) {
    pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
  } else {
    pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
  }
 
  CHKERR DMoFEMPreProcessFiniteElements(simple->getDM(),
                                        post_proc_begin->getFEMethod());
  CHKERR pip->loopFiniteElements();
  CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
                                         post_proc_end->getFEMethod());
 
  CHKERR post_proc_end->writeFile("out_elastic_" + std::to_string(iter) +
                                  ".h5m");
  MoFEMFunctionReturn(0);
}
//! [Postprocess results]
 
//! [calculateGradient]
 
using DomainBaseOp = FormsIntegrators<DomainEleOp>::Assembly<A>::OpBase;
 
template <int SPACE_DIM>
struct OpAdJointTestOp<SPACE_DIM, IntegrationType::GAUSS, DomainBaseOp>
    : public DomainBaseOp {
 
  using OP = DomainBaseOp;
 
  OpAdJointTestOp(const std::string field_name,
                  boost::shared_ptr<HookeOps::CommonData> comm_ptr)
      : OP(field_name, field_name, DomainEleOp::OPROW), commPtr(comm_ptr) {}
 
protected:
  boost::shared_ptr<HookeOps::CommonData> commPtr;
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
};
 
template <int SPACE_DIM>
MoFEMErrorCode
OpAdJointTestOp<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);
 
  // get element volume
  const double vol = OP::getMeasure();
  // get integration weights
  auto t_w = OP::getFTensor0IntegrationWeight();
  // get base function gradient on rows
  auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
  // get stress
  auto t_cauchy_stress =
      getFTensor2SymmetricFromMat<SPACE_DIM>(*(commPtr->getMatCauchyStress()));
  // loop over integration points
  for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
    // take into account Jacobian
    const double alpha = t_w * vol;
 
    // get rhs vector
    auto t_nf = OP::template getNf<SPACE_DIM>();
    // loop over rows base functions
    int rr = 0;
    for (; rr != OP::nbRows / SPACE_DIM; rr++) {
      // Variation due to change in geometry (diff base)
      t_nf(j) += alpha * t_row_grad(i) * t_cauchy_stress(i, j);
 
      ++t_row_grad;
      ++t_nf;
    }
    for (; rr < OP::nbRowBaseFunctions; ++rr) {
      ++t_row_grad;
    }
 
    ++t_cauchy_stress;
    ++t_w; // move to another integration weight
  }
  MoFEMFunctionReturn(0);
}
 
template <int SPACE_DIM>
struct OpAdJointGradTimesSymTensor<SPACE_DIM, IntegrationType::GAUSS,
                                   DomainBaseOp> : public DomainBaseOp {
 
  using OP = DomainBaseOp;
 
  OpAdJointGradTimesSymTensor(const std::string field_name,
                              boost::shared_ptr<HookeOps::CommonData> comm_ptr,
                              boost::shared_ptr<MatrixDouble> jac,
                              boost::shared_ptr<MatrixDouble> diff_jac,
                              boost::shared_ptr<VectorDouble> cof_vals)
      : OP(field_name, field_name, DomainEleOp::OPROW), commPtr(comm_ptr),
        jac(jac), diffJac(diff_jac), cofVals(cof_vals) {}
 
protected:
  boost::shared_ptr<HookeOps::CommonData> commPtr;
  boost::shared_ptr<MatrixDouble> jac;
  boost::shared_ptr<MatrixDouble> diffJac;
  boost::shared_ptr<VectorDouble> cofVals;
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
};
 
template <int DIM> inline auto diff_symmetrize(FTensor::Number<DIM>) {
 
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  FTensor::Index<'k', DIM> k;
  FTensor::Index<'l', DIM> l;
 
  FTensor::Tensor4<double, DIM, DIM, DIM, DIM> t_diff;
 
  t_diff(i, j, k, l) = 0;
  t_diff(0, 0, 0, 0) = 1;
  t_diff(1, 1, 1, 1) = 1;
 
  t_diff(1, 0, 1, 0) = 0.5;
  t_diff(1, 0, 0, 1) = 0.5;
 
  t_diff(0, 1, 0, 1) = 0.5;
  t_diff(0, 1, 1, 0) = 0.5;
 
  if constexpr (DIM == 3) {
    t_diff(2, 2, 2, 2) = 1;
 
    t_diff(2, 0, 2, 0) = 0.5;
    t_diff(2, 0, 0, 2) = 0.5;
    t_diff(0, 2, 0, 2) = 0.5;
    t_diff(0, 2, 2, 0) = 0.5;
 
    t_diff(2, 1, 2, 1) = 0.5;
    t_diff(2, 1, 1, 2) = 0.5;
    t_diff(1, 2, 1, 2) = 0.5;
    t_diff(1, 2, 2, 1) = 0.5;
  }
 
  return t_diff;
};
 
template <int SPACE_DIM>
MoFEMErrorCode
OpAdJointGradTimesSymTensor<SPACE_DIM, IntegrationType::GAUSS, DomainBaseOp>::
    iNtegrate(EntitiesFieldData::EntData &row_data) {
  MoFEMFunctionBegin;
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  FTENSOR_INDEX(SPACE_DIM, k);
  FTENSOR_INDEX(SPACE_DIM, l);
  FTENSOR_INDEX(SPACE_DIM, I);
  FTENSOR_INDEX(SPACE_DIM, J);
 
  auto t_diff_symm = diff_symmetrize(FTensor::Number<SPACE_DIM>());
 
  // get element volume
  const double vol = OP::getMeasure();
  // get integration weights
  auto t_w = OP::getFTensor0IntegrationWeight();
  // get Jacobian values
  auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(jac));
  // get diff Jacobian values
  auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(diffJac));
  // get cofactor values
  auto t_cof = getFTensor0FromVec(*(cofVals));
  // get base function gradient on rows
  auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
  // get fradient of the field
  auto t_grad_u =
      getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(commPtr->matGradPtr));
  // get field gradient values
  auto t_cauchy_stress =
      getFTensor2SymmetricFromMat<SPACE_DIM>(*(commPtr->getMatCauchyStress()));
  // material stiffness tensor
  auto t_D =
      getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
  // loop over integration points
  for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
    // take into account Jacobian
    const double alpha = t_w * vol;
 
    auto t_det = determinantTensor(t_jac);
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_inv_jac;
    CHKERR invertTensor(t_jac, t_det, t_inv_jac);
 
    // Calculate the variation of the gradient due to geometry change
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_inv_jac;
    t_diff_inv_jac(i, j) =
        -(t_inv_jac(i, I) * t_diff_jac(I, J)) * t_inv_jac(J, j);
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_grad;
    t_diff_grad(i, j) = t_grad_u(i, k) * t_diff_inv_jac(k, j);
 
    // Calculate the variation of the strain tensor
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_strain;
    t_diff_strain(i, j) = t_diff_symm(i, j, k, l) * t_diff_grad(k, l);
 
    // Calculate the variation of the stress tensor
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_stress;
    t_diff_stress(i, j) = t_D(i, j, k, l) * t_diff_strain(k, l);
 
    // get rhs vector
    auto t_nf = OP::template getNf<SPACE_DIM>();
    // loop over rows base functions
    int rr = 0;
    for (; rr != OP::nbRows / SPACE_DIM; rr++) {
 
      FTensor::Tensor1<double, SPACE_DIM> t_diff_row_grad;
      t_diff_row_grad(k) = t_row_grad(j) * t_diff_inv_jac(j, k);
 
      // Variation due to change in geometry (diff base)
      t_nf(j) += alpha * t_diff_row_grad(i) * t_cauchy_stress(i, j);
 
      // Variation due to change in domain (cofactor)
      t_nf(j) += (alpha * t_cof) * t_row_grad(i) * t_cauchy_stress(i, j);
 
      // Variation due to change in stress (diff stress)
      t_nf(j) += alpha * t_row_grad(i) * t_diff_stress(i, j);
 
      ++t_row_grad;
      ++t_nf;
    }
    for (; rr < OP::nbRowBaseFunctions; ++rr) {
      ++t_row_grad;
    }
 
    ++t_grad_u;
    ++t_cauchy_stress;
    ++t_jac;
    ++t_diff_jac;
    ++t_cof;
    ++t_w; // move to another integration weight
  }
  MoFEMFunctionReturn(0);
}
 
struct OpAdJointObjective : public ForcesAndSourcesCore::UserDataOperator {
  using OP = ForcesAndSourcesCore::UserDataOperator;
  OpAdJointObjective(boost::shared_ptr<ObjectiveFunctionData> python_ptr,
                     boost::shared_ptr<HookeOps::CommonData> comm_ptr,
                     boost::shared_ptr<MatrixDouble> jac_ptr,
                     boost::shared_ptr<MatrixDouble> diff_jac,
                     boost::shared_ptr<VectorDouble> cof_vals,
                     boost::shared_ptr<MatrixDouble> d_grad_ptr,
                     boost::shared_ptr<MatrixDouble> u_ptr,
                     boost::shared_ptr<double> glob_objective_ptr,
                     boost::shared_ptr<double> glob_objective_grad_ptr)
      : OP(NOSPACE, OP::OPSPACE), pythonPtr(python_ptr), commPtr(comm_ptr),
        jacPtr(jac_ptr), diffJacPtr(diff_jac), cofVals(cof_vals),
        dGradPtr(d_grad_ptr), uPtr(u_ptr), globObjectivePtr(glob_objective_ptr),
        globObjectiveGradPtr(glob_objective_grad_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &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);
 
    constexpr auto symm_size = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
 
    auto t_diff_symm = diff_symmetrize(FTensor::Number<SPACE_DIM>());
 
    auto nb_gauss_pts = getGaussPts().size2();
    auto objective_ptr = boost::make_shared<VectorDouble>(nb_gauss_pts);
    auto objective_dstress =
        boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
    auto objective_dstrain =
        boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
    auto obj_grad = boost::make_shared<MatrixDouble>(SPACE_DIM, nb_gauss_pts);
 
    auto evaluate_python = [&]() {
      MoFEMFunctionBegin;
      auto &coords = OP::getCoordsAtGaussPts();
      CHKERR pythonPtr->evalObjectiveFunction(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_ptr);
      CHKERR pythonPtr->evalObjectiveGradientStress(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_dstress);
      CHKERR pythonPtr->evalObjectiveGradientStrain(
          coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
          objective_dstrain);
 
      auto t_grad_u =
          getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(commPtr->matGradPtr));
      auto t_D =
          getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
      auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(jacPtr));
      auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(diffJacPtr));
      auto t_cof = getFTensor0FromVec(*(cofVals));
      auto t_d_grad = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(dGradPtr));
 
      auto t_obj = getFTensor0FromVec(*objective_ptr);
      auto t_obj_dstress =
          getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstress);
      auto t_obj_dstrain =
          getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstrain);
 
      auto vol = OP::getMeasure();
      auto t_w = getFTensor0IntegrationWeight();
      for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
 
        auto t_det = determinantTensor(t_jac);
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_inv_jac;
        CHKERR invertTensor(t_jac, t_det, t_inv_jac);
 
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_inv_jac;
        t_diff_inv_jac(i, j) =
            -(t_inv_jac(i, I) * t_diff_jac(I, J)) * t_inv_jac(J, j);
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_grad;
        t_diff_grad(i, j) = t_grad_u(i, k) * t_diff_inv_jac(k, j);
 
        FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_d_strain;
        t_d_strain(i, j) = t_diff_symm(i, j, k, l) * (
 
                                                         t_d_grad(k, l)
 
                                                         +
 
                                                         t_diff_grad(k, l)
 
                                                     );
 
        auto alpha = t_w * vol;
 
        (*globObjectivePtr) += alpha * t_obj;
        (*globObjectiveGradPtr) +=
            alpha *
            (
 
                t_obj_dstress(i, j) * (t_D(i, j, k, l) * t_d_strain(k, l))
 
                +
 
                t_obj_dstrain(i, j) * t_d_strain(i, j)
 
                +
 
                t_obj * t_cof
 
            );
 
        ++t_w;
        ++t_jac;
        ++t_diff_jac;
        ++t_cof;
 
        ++t_obj;
        ++t_obj_dstress;
        ++t_obj_dstrain;
 
        ++t_grad_u;
        ++t_d_grad;
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR evaluate_python();
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
  boost::shared_ptr<HookeOps::CommonData> commPtr;
  boost::shared_ptr<MatrixDouble> jacPtr;
  boost::shared_ptr<MatrixDouble> diffJacPtr;
  boost::shared_ptr<VectorDouble> cofVals;
  boost::shared_ptr<MatrixDouble> dGradPtr;
  boost::shared_ptr<MatrixDouble> uPtr;
 
  boost::shared_ptr<double> globObjectivePtr;
  boost::shared_ptr<double> globObjectiveGradPtr;
};
 
MoFEMErrorCode Example::calculateGradient(PetscReal *objective_function_value,
                                          Vec objective_function_gradient,
                                          Vec adjoint_vector) {
  MOFEM_LOG_CHANNEL("WORLD");
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  auto ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
 
  auto get_essential_fe = [this]() {
    auto post_proc_rhs = boost::make_shared<FEMethod>();
    auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
                                                           post_proc_rhs, 0)();
      MoFEMFunctionReturnHot(0);
    };
    post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
    return post_proc_rhs;
  };
 
  auto get_fd_direvative_fe = [&]() {
    auto fe = boost::make_shared<DomainEle>(mField);
    fe->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe->getOpPtrVector();
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
    // Add RHS operators for internal forces
    constexpr bool debug = false;
    if constexpr (debug) {
      auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
          mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
      pip.push_back(
          new OpAdJointTestOp<SPACE_DIM, I, DomainBaseOp>("U", common_ptr));
    } else {
      HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
          mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
    }
    return fe;
  };
 
  auto calulate_fd_residual = [&](auto eps, auto diff_vec, auto fd_vec) {
    MoFEMFunctionBegin;
 
    constexpr bool debug = false;
 
    auto geom_norm = [](MoFEM::Interface &mField) {
      MoFEMFunctionBegin;
      auto field_blas = mField.getInterface<FieldBlas>();
      double nrm2 = 0.0;
      auto norm2_field = [&](const double val) {
        nrm2 += val * val;
        return val;
      };
      CHKERR field_blas->fieldLambdaOnValues(norm2_field, "GEOMETRY");
      MPI_Allreduce(MPI_IN_PLACE, &nrm2, 1, MPI_DOUBLE, MPI_SUM,
                    mField.get_comm());
      MOFEM_LOG("WORLD", Sev::inform) << "Geometry norm: " << sqrt(nrm2);
      MoFEMFunctionReturn(0);
    };
 
    if constexpr (debug)
      CHKERR geom_norm(mField);
 
    auto initial_current_geometry = createDMVector(adjointDM);
    CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
        "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
        initial_current_geometry, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecAssemblyBegin(initial_current_geometry);
    CHKERR VecAssemblyEnd(initial_current_geometry);
 
    if constexpr (debug)
      CHKERR geom_norm(mField);
 
    auto perturb_geometry = [&](auto eps, auto diff_vec) {
      MoFEMFunctionBegin;
      auto current_geometry = vectorDuplicate(initial_current_geometry);
      CHKERR VecCopy(initial_current_geometry, current_geometry);
      CHKERR VecAXPY(current_geometry, eps, diff_vec);
      CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
          "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
          current_geometry, INSERT_VALUES, SCATTER_REVERSE);
      MoFEMFunctionReturn(0);
    };
 
    auto fe = get_fd_direvative_fe();
    auto fp = vectorDuplicate(diff_vec);
    auto fm = vectorDuplicate(diff_vec);
    auto calc_impl = [&](auto f, auto eps) {
      MoFEMFunctionBegin;
      CHKERR VecZeroEntries(f);
      fe->f = f;
      CHKERR perturb_geometry(eps, diff_vec);
      CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
                                      simple->getDomainFEName(), fe);
      CHKERR VecAssemblyBegin(f);
      CHKERR VecAssemblyEnd(f);
      CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
      auto post_proc_rhs = get_essential_fe();
      post_proc_rhs->f = f;
      CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
                                             post_proc_rhs.get());
      MoFEMFunctionReturn(0);
    };
    CHKERR calc_impl(fp, eps);
    CHKERR calc_impl(fm, -eps);
    CHKERR VecWAXPY(fd_vec, -1.0, fm, fp);
    CHKERR VecScale(fd_vec, 1.0 / (2.0 * eps));
 
    CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
        "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
        initial_current_geometry, INSERT_VALUES, SCATTER_REVERSE);
 
    if constexpr (debug)
      CHKERR geom_norm(mField);
 
    MoFEMFunctionReturn(0);
  };
 
  auto get_direvative_fe = [&](auto diff_vec) {
    auto fe_adjoint = boost::make_shared<DomainEle>(mField);
    fe_adjoint->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe_adjoint->getOpPtrVector();
 
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
    auto cof_ptr = boost::make_shared<VectorDouble>();
 
    using OpCoFactor =
        AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
 
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
 
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "GEOMETRY", jac_ptr));
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "U", diff_jac_ptr,
        diff_vec)); // Note: that vector is stored on displacemnt vector, that
                    // why is used here
 
    pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
 
    // Add RHS operators for internal forces
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
        mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
    pip.push_back(new OpAdJointGradTimesSymTensor<SPACE_DIM, I, DomainBaseOp>(
        "U", common_ptr, jac_ptr, diff_jac_ptr, cof_ptr));
 
    return fe_adjoint;
  };
 
  auto get_objective_fe = [&](auto diff_vec, auto grad_vec,
                              auto glob_objective_ptr,
                              auto glob_objective_grad_ptr) {
    auto fe_adjoint = boost::make_shared<DomainEle>(mField);
    fe_adjoint->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data +  p_data - 1;
    };
    auto &pip = fe_adjoint->getOpPtrVector();
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
 
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
    auto cof_ptr = boost::make_shared<VectorDouble>();
    auto d_grad_ptr = boost::make_shared<MatrixDouble>();
    auto u_ptr = boost::make_shared<MatrixDouble>();
 
    using OpCoFactor =
        AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "GEOMETRY", jac_ptr));
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "U", diff_jac_ptr,
        diff_vec)); // Note: that vector is stored on displacemnt vector, that
                    // why is used here
    pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "U", d_grad_ptr, grad_vec));
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
 
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
        mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
    pip.push_back(new OpAdJointObjective(
        pythonPtr, common_ptr, jac_ptr, diff_jac_ptr, cof_ptr, d_grad_ptr,
        u_ptr, glob_objective_ptr, glob_objective_grad_ptr));
 
    return fe_adjoint;
  };
 
  auto dm = simple->getDM();
  auto f = createDMVector(dm);
  auto d = vectorDuplicate(f);
  auto dm_diff_vec = vectorDuplicate(d);
 
  auto adjoint_fe = get_direvative_fe(dm_diff_vec);
  auto glob_objective_ptr = boost::make_shared<double>(0.0);
  auto glob_objective_grad_ptr = boost::make_shared<double>(0.0);
  auto objective_fe = get_objective_fe(dm_diff_vec, d, glob_objective_ptr,
                                       glob_objective_grad_ptr);
 
  auto set_variance_of_geometry = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    CHKERR DMoFEMMeshToLocalVector(adjointDM, mod_vec, INSERT_VALUES,
                                   SCATTER_REVERSE);
    CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
        simple->getProblemName(), "U", "ADJOINT_FIELD", RowColData::ROW,
        dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateBegin(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_variance_internal_forces = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    CHKERR VecZeroEntries(f);
    CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
    adjoint_fe->f = f;
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), adjoint_fe);
    CHKERR VecAssemblyBegin(f);
    CHKERR VecAssemblyEnd(f);
    CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
    auto post_proc_rhs = get_essential_fe();
    post_proc_rhs->f = f;
    CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
                                           post_proc_rhs.get());
    CHKERR VecScale(f, -1.0);
 
#ifndef NDEBUG
    constexpr bool debug = false;
    if constexpr (debug) {
      double norm0;
      CHKERR VecNorm(f, NORM_2, &norm0);
      auto fd_check = vectorDuplicate(f);
      double eps = 1e-5;
      CHKERR calulate_fd_residual(eps, dm_diff_vec, fd_check);
      double nrm;
      CHKERR VecAXPY(fd_check, -1.0, f);
      CHKERR VecNorm(fd_check, NORM_2, &nrm);
      MOFEM_LOG("WORLD", Sev::inform)
          << " FD check for internal forces [ " << mode << " ]: " << nrm
          << " / " << norm0 << " ( " << (nrm / norm0) << " )";
    }
#endif
 
    MoFEMFunctionReturn(0);
  };
 
  auto calulate_variance_of_displacemnte = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    CHKERR KSPSolve(kspElastic, f, d);
    CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_variance_of_objective_function = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    *glob_objective_ptr = 0.0;
    *glob_objective_grad_ptr = 0.0;
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                    objective_fe);
    CHKERR VecSetValue(objective_function_gradient, mode,
                       *glob_objective_grad_ptr, ADD_VALUES);
    std::array<double, 2> array = {*glob_objective_ptr,
                                   *glob_objective_grad_ptr};
    MPI_Allreduce(MPI_IN_PLACE, array.data(), 2, MPI_DOUBLE, MPI_SUM,
                  mField.get_comm());
    *glob_objective_ptr = array[0];
    *glob_objective_grad_ptr = array[1];
    (*objective_function_value) += *glob_objective_ptr;
    MOFEM_LOG_CHANNEL("WORLD");
    MOFEM_LOG("WORLD", Sev::verbose) << " Objective gradient [ " << mode
                                     << " ]: " << *glob_objective_grad_ptr;
    MoFEMFunctionReturn(0);
  };
 
  CHKERR VecZeroEntries(objective_function_gradient);
  CHKERR VecZeroEntries(adjoint_vector);
  *objective_function_value = 0.0;
  int mode = 0;
  for (auto mod_vec : modeVecs) {
 
    CHKERR set_variance_of_geometry(mode, mod_vec);
    CHKERR calculate_variance_internal_forces(mode, mod_vec);
    CHKERR calulate_variance_of_displacemnte(mode, mod_vec);
    CHKERR calculate_variance_of_objective_function(mode, mod_vec);
 
    CHKERR VecAXPY(adjoint_vector, *glob_objective_grad_ptr, dm_diff_vec);
    ++mode;
  }
 
  (*objective_function_value) /= modeVecs.size();
  MOFEM_LOG("WORLD", Sev::verbose)
      << "Objective function: " << *glob_objective_ptr;
 
  CHKERR VecAssemblyBegin(objective_function_gradient);
  CHKERR VecAssemblyEnd(objective_function_gradient);
 
  CHKERR VecAssemblyBegin(adjoint_vector);
  CHKERR VecAssemblyEnd(adjoint_vector);
  CHKERR VecGhostUpdateBegin(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
 
  MoFEMFunctionReturn(0);
}
//! [calculateGradient]
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
  Py_Initialize();
  np::initialize();
 
  // Initialisation of 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);
  }
}
 
struct ObjectiveFunctionDataImpl : public ObjectiveFunctionData {
  ObjectiveFunctionDataImpl() = default;
  virtual ~ObjectiveFunctionDataImpl() = default;
 
  MoFEMErrorCode initPython(const std::string py_file);
 
  MoFEMErrorCode
  evalObjectiveFunction(MatrixDouble &coords,
                        boost::shared_ptr<MatrixDouble> u_ptr,
                        boost::shared_ptr<MatrixDouble> stress_ptr,
                        boost::shared_ptr<MatrixDouble> strain_ptr,
                        boost::shared_ptr<VectorDouble> o_ptr);
 
  MoFEMErrorCode
  evalObjectiveGradientStress(MatrixDouble &coords,
                              boost::shared_ptr<MatrixDouble> u_ptr,
                              boost::shared_ptr<MatrixDouble> stress_ptr,
                              boost::shared_ptr<MatrixDouble> strain_ptr,
                              boost::shared_ptr<MatrixDouble> o_ptr);
 
  MoFEMErrorCode
  evalObjectiveGradientStrain(MatrixDouble &coords,
                              boost::shared_ptr<MatrixDouble> u_ptr,
                              boost::shared_ptr<MatrixDouble> stress_ptr,
                              boost::shared_ptr<MatrixDouble> strain_ptr,
                              boost::shared_ptr<MatrixDouble> o_ptr);
 
  MoFEMErrorCode
  evalObjectiveGradientU(MatrixDouble &coords,
                              boost::shared_ptr<MatrixDouble> u_ptr,
                              boost::shared_ptr<MatrixDouble> stress_ptr,
                              boost::shared_ptr<MatrixDouble> strain_ptr,
                              boost::shared_ptr<MatrixDouble> o_ptr);
 
  MoFEMErrorCode numberOfModes(int block_id, int &modes);
 
  MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords,
                            std::array<double, 3> &centroid,
                            std::array<double, 6> &bbodx, MatrixDouble &o_ptr);
 
private:
  bp::object mainNamespace;
  bp::object objectiveFunction;
  bp::object objectiveGradientStress;
  bp::object objectiveGradientStrain;
  bp::object objectiveGradientU;
  bp::object objectNumberOfModes;
  bp::object objectBlockModes;
 
  MoFEMErrorCode objectiveFunctionImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  MoFEMErrorCode objectiveGradientStressImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  MoFEMErrorCode objectiveGradientStrainImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  MoFEMErrorCode objectiveGradientUImpl(
 
      np::ndarray coords, np::ndarray u,
 
      np::ndarray stress, np::ndarray strain, np::ndarray &o
 
  );
 
  MoFEMErrorCode blockModesImpl(
 
      int block_id, np::ndarray coords, np::ndarray centroid, np::ndarray bbodx,
      np::ndarray &o_ptr
 
  );
 
  np::ndarray convertToNumPy(std::vector<double> &data, int rows,
                             int nb_gauss_pts);
 
  np::ndarray convertToNumPy(double *ptr, int s);
 
  MatrixDouble copyToFull(MatrixDouble &s);
  void copyToSymmetric(double *ptr, MatrixDouble &s);
};
 
boost::shared_ptr<ObjectiveFunctionData>
create_python_objective_function(std::string py_file) {
  auto ptr = boost::make_shared<ObjectiveFunctionDataImpl>();
  CHK_THROW_MESSAGE(ptr->initPython(py_file), "init python");
  return ptr;
}
 
MoFEMErrorCode
ObjectiveFunctionDataImpl::initPython(const std::string py_file) {
  MoFEMFunctionBegin;
  try {
 
    // create main module
    auto main_module = bp::import("__main__");
    mainNamespace = main_module.attr("__dict__");
    bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
    // create a reference to python function
    objectiveFunction = mainNamespace["f"];
    objectiveGradientStress = mainNamespace["f_stress"];
    objectiveGradientStrain = mainNamespace["f_strain"];
    objectiveGradientU = mainNamespace["f_u"];
    objectNumberOfModes = mainNamespace["number_of_modes"];
    objectBlockModes = mainNamespace["block_modes"];
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MatrixDouble ObjectiveFunctionDataImpl::copyToFull(MatrixDouble &s) {
  MatrixDouble f(9, s.size2());
  f.clear();
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  auto t_f = getFTensor2FromMat<3, 3>(f);
  auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
  for (int ii = 0; ii != s.size2(); ++ii) {
    t_f(i, j) = t_s(i, j);
    ++t_f;
    ++t_s;
  }
  return f;
};
 
void ObjectiveFunctionDataImpl::copyToSymmetric(double *ptr, MatrixDouble &s) {
  FTENSOR_INDEX(SPACE_DIM, i);
  FTENSOR_INDEX(SPACE_DIM, j);
  auto t_f = FTensor::Tensor2<FTensor::PackPtr<double *, 1>, 3, 3>{
 
      ptr + 0 * s.size2(), ptr + 1 * s.size2(),
      ptr + 2 * s.size2(), ptr + 3 * s.size2(),
      ptr + 4 * s.size2(), ptr + 5 * s.size2(),
      ptr + 6 * s.size2(), ptr + 7 * s.size2(),
      ptr + 8 * s.size2()
 
  };
  auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
  for (int ii = 0; ii != s.size2(); ++ii) {
    t_s(i, j) = (t_f(i, j) || t_f(j, i)) / 2.0;
    ++t_f;
    ++t_s;
  }
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::evalObjectiveFunction(
    MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
    boost::shared_ptr<MatrixDouble> stress_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    boost::shared_ptr<VectorDouble> o_ptr) {
  MoFEMFunctionBegin;
  try {
 
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
    auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
 
    auto full_stress = copyToFull(*(stress_ptr));
    auto full_strain = copyToFull(*(strain_ptr));
 
    auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
                                    full_stress.size2());
    auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
                                    full_strain.size2());
    np::ndarray np_output = np::empty(bp::make_tuple(strain_ptr->size2()),
                                      np::dtype::get_builtin<double>());
    CHKERR objectiveFunctionImpl(np_coords, np_u, np_stress, np_strain,
                                 np_output);
    o_ptr->resize(stress_ptr->size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    std::copy(val_ptr, val_ptr + strain_ptr->size2(), o_ptr->data().begin());
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::evalObjectiveGradientStress(
    MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
    boost::shared_ptr<MatrixDouble> stress_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    boost::shared_ptr<MatrixDouble> o_ptr) {
  MoFEMFunctionBegin;
  try {
 
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
    auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
 
    auto full_stress = copyToFull(*(stress_ptr));
    auto full_strain = copyToFull(*(strain_ptr));
 
    auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
                                    full_stress.size2());
    auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
                                    full_strain.size2());
    np::ndarray np_output =
        np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
                  np::dtype::get_builtin<double>());
    CHKERR objectiveGradientStressImpl(np_coords, np_u, np_stress, np_strain,
                                       np_output);
    o_ptr->resize(stress_ptr->size1(), stress_ptr->size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    copyToSymmetric(val_ptr, *(o_ptr));
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::evalObjectiveGradientStrain(
    MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
    boost::shared_ptr<MatrixDouble> stress_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    boost::shared_ptr<MatrixDouble> o_ptr) {
  MoFEMFunctionBegin;
  try {
 
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
    auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
 
    auto full_stress = copyToFull(*(stress_ptr));
    auto full_strain = copyToFull(*(strain_ptr));
    auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
                                    full_stress.size2());
    auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
                                    full_strain.size2());
 
    np::ndarray np_output =
        np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
                  np::dtype::get_builtin<double>());
    CHKERR objectiveGradientStrainImpl(np_coords, np_u, np_stress, np_strain,
                                       np_output);
    o_ptr->resize(strain_ptr->size1(), strain_ptr->size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    copyToSymmetric(val_ptr, *(o_ptr));
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::evalObjectiveGradientU(
    MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
    boost::shared_ptr<MatrixDouble> stress_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    boost::shared_ptr<MatrixDouble> o_ptr) {
  MoFEMFunctionBegin;
  try {
 
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
    auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
 
    auto full_stress = copyToFull(*(stress_ptr));
    auto full_strain = copyToFull(*(strain_ptr));
    auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
                                    full_stress.size2());
    auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
                                    full_strain.size2());
 
    np::ndarray np_output =
        np::empty(bp::make_tuple(u_ptr->size1(), u_ptr->size2()),
                  np::dtype::get_builtin<double>());
    CHKERR objectiveGradientStrainImpl(np_coords, np_u, np_stress, np_strain,
                                       np_output);
    o_ptr->resize(u_ptr->size1(), u_ptr->size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    std::copy(val_ptr, val_ptr + u_ptr->size1() * u_ptr->size2(),
              o_ptr->data().begin());
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::blockModes(
    int block_id, MatrixDouble &coords, std::array<double, 3> &centroid,
    std::array<double, 6> &bbodx, MatrixDouble &o_ptr) {
  MoFEMFunctionBegin;
  try {
 
    int nb_modes = bp::extract<int>(objectNumberOfModes(block_id));
 
    auto np_coords =
        convertToNumPy(coords.data(), coords.size1(), coords.size2());
 
    auto np_centroid = convertToNumPy(centroid.data(), 3);
    auto np_bbodx = convertToNumPy(bbodx.data(), 6);
 
    np::ndarray np_output =
        np::empty(bp::make_tuple(coords.size1(), nb_modes, 3),
                  np::dtype::get_builtin<double>());
    CHKERR blockModesImpl(block_id, np_coords, np_centroid, np_bbodx,
                          np_output);
 
    o_ptr.resize(nb_modes, coords.size1() * coords.size2(), false);
    double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
    std::copy(val_ptr, val_ptr + coords.size1() * coords.size2() * nb_modes,
              o_ptr.data().begin());
 
  } catch (bp::error_already_set const &) {
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveFunctionImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  MoFEMFunctionBegin;
  try {
 
    // call python function
    o = bp::extract<np::ndarray>(objectiveFunction(coords, u, stress, strain));
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveGradientStressImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  MoFEMFunctionBegin;
  try {
 
    // call python function
    o = bp::extract<np::ndarray>(
        objectiveGradientStress(coords, u, stress, strain));
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveGradientStrainImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  MoFEMFunctionBegin;
  try {
 
    // call python function
    o = bp::extract<np::ndarray>(
        objectiveGradientStrain(coords, u, stress, strain));
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::objectiveGradientUImpl(
 
    np::ndarray coords, np::ndarray u,
 
    np::ndarray stress, np::ndarray strain, np::ndarray &o
 
) {
  MoFEMFunctionBegin;
  try {
 
    // call python function
    o = bp::extract<np::ndarray>(objectiveGradientU(coords, u, stress, strain));
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::numberOfModes(int block_id,
                                                        int &modes) {
  MoFEMFunctionBegin;
  try {
 
    modes = bp::extract<int>(objectNumberOfModes(block_id));
 
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ObjectiveFunctionDataImpl::blockModesImpl(int block_id,
                                                         np::ndarray coords,
                                                         np::ndarray centroid,
                                                         np::ndarray bbodx,
                                                         np::ndarray &o) {
  MoFEMFunctionBegin;
  try {
    // call python function
    o = bp::extract<np::ndarray>(
        objectBlockModes(block_id, coords, centroid, bbodx));
  } catch (bp::error_already_set const &) {
    // print all other errors to stderr
    PyErr_Print();
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
  }
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Converts a std::vector<double> to a NumPy ndarray.
 *
 * This function wraps the given vector data into a NumPy array with the
 * specified number of rows and Gauss points. The resulting ndarray shares
 * memory with the input vector, so changes to one will affect the other.
 *
 * @param data Reference to the vector containing double values to be converted.
 * @param rows Number of rows in the resulting NumPy array.
 * @param nb_gauss_pts Number of Gauss points (columns) in the resulting NumPy
 * array.
 * @return np::ndarray NumPy array view of the input data.
 *
 * @note
 * - `size` specifies the shape of the resulting ndarray as a tuple (rows,
 * nb_gauss_pts).
 * - `stride` specifies the step size in bytes to move to the next element in
 * memory. Here, it is set to sizeof(double), indicating contiguous storage for
 * each element.
 */
inline np::ndarray
ObjectiveFunctionDataImpl::convertToNumPy(std::vector<double> &data, int rows,
                                          int nb_gauss_pts) {
  auto dtype = np::dtype::get_builtin<double>();
  auto size = bp::make_tuple(rows, nb_gauss_pts);
  auto stride = bp::make_tuple(nb_gauss_pts * sizeof(double), sizeof(double));
  return (np::from_data(data.data(), dtype, size, stride, bp::object()));
}
 
inline np::ndarray ObjectiveFunctionDataImpl::convertToNumPy(double *ptr,
                                                             int s) {
  auto dtype = np::dtype::get_builtin<double>();
  auto size = bp::make_tuple(s);
  auto stride = bp::make_tuple(sizeof(double));
  return (np::from_data(ptr, dtype, size, stride, bp::object()));
}
 
Range get_range_from_block(MoFEM::Interface &m_field,
                           const std::string block_name, int dim) {
  Range r;
 
  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
  auto bcs = mesh_mng->getCubitMeshsetPtr(
 
      std::regex((boost::format("%s(.*)") % block_name).str())
 
  );
 
  for (auto bc : bcs) {
    Range faces;
    CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
                                                       faces, true),
                   "get meshset ents");
    r.merge(faces);
  }
 
  for (auto dd = dim - 1; dd >= 0; --dd) {
    if (dd >= 0) {
      Range ents;
      CHK_MOAB_THROW(m_field.get_moab().get_adjacencies(r, dd, false, ents,
                                                        moab::Interface::UNION),
                     "get adjs");
      r.merge(ents);
    } else {
      Range verts;
      CHK_MOAB_THROW(m_field.get_moab().get_connectivity(r, verts),
                     "get verts");
      r.merge(verts);
    }
    CHK_MOAB_THROW(
        m_field.getInterface<CommInterface>()->synchroniseEntities(r), "comm");
  }
 
  return r;
};
 
MoFEMErrorCode save_range(moab::Interface &moab, const std::string name,
                          const Range r) {
  MoFEMFunctionBegin;
  auto out_meshset = get_temp_meshset_ptr(moab);
  CHKERR moab.add_entities(*out_meshset, r);
  CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
  MoFEMFunctionReturn(0);
};