MySources/simulation_structs.hpp

#ifndef SIMULATIONSTRUCTS_HPP
#define SIMULATIONSTRUCTS_HPP

namespace Eigen {
template <class Matrix>
void write_binary(const std::string &filename, const Matrix &matrix) {
  std::ofstream out(filename,
                    std::ios::out | std::ios::binary | std::ios::trunc);
  typename Matrix::Index rows = matrix.rows(), cols = matrix.cols();
  out.write((char *)(&rows), sizeof(typename Matrix::Index));
  out.write((char *)(&cols), sizeof(typename Matrix::Index));
  out.write((char *)matrix.data(),
            rows * cols * sizeof(typename Matrix::Scalar));
  out.close();
}
template <class Matrix>
void read_binary(const std::string &filename, Matrix &matrix) {
  std::ifstream in(filename, std::ios::in | std::ios::binary);
  typename Matrix::Index rows = 0, cols = 0;
  in.read((char *)(&rows), sizeof(typename Matrix::Index));
  in.read((char *)(&cols), sizeof(typename Matrix::Index));
  matrix.resize(rows, cols);
  in.read((char *)matrix.data(), rows * cols * sizeof(typename Matrix::Scalar));
  in.close();
}
} // namespace Eigen

#include "simulationmesh.hpp"
#include "nlohmann/json.hpp"
#include <string>
#include <vector>

struct SimulationHistory {
  SimulationHistory() {}

  size_t numberOfSteps{0};
  std::string label;
  std::vector<double> logResidualForces;
  std::vector<double> kineticEnergy;
  std::vector<double> potentialEnergies;
  std::vector<size_t> redMarks;
  std::vector<double> greenMarks;
  std::vector<double> residualForcesMovingAverage;
  std::vector<double> sumOfNormalizedDisplacementNorms;
  //  std::vector<double> residualForcesMovingAverageDerivativesLog;

  void markRed(const size_t &stepNumber) { redMarks.push_back(stepNumber); }

  void markGreen(const size_t &stepNumber) { greenMarks.push_back(stepNumber); }

  void stepPulse(const SimulationMesh &mesh)
  {
      kineticEnergy.push_back(std::log10(mesh.currentTotalKineticEnergy));
      //    potentialEnergy.push_back(mesh.totalPotentialEnergykN);
      logResidualForces.push_back(std::log10(mesh.totalResidualForcesNorm));
      residualForcesMovingAverage.push_back(std::log10(mesh.residualForcesMovingAverage));
      sumOfNormalizedDisplacementNorms.push_back(std::log10(mesh.sumOfNormalizedDisplacementNorms));
      //    residualForcesMovingAverageDerivativesLog.push_back(
      //        std::log(mesh.residualForcesMovingAverageDerivativeNorm));
  }

  void clear()
  {
      logResidualForces.clear();
      kineticEnergy.clear();
      potentialEnergies.clear();
      residualForcesMovingAverage.clear();
      sumOfNormalizedDisplacementNorms.clear();
      //      residualForcesMovingAverageDerivativesLog.clear();
  }
};


namespace nlohmann {

template <> struct adl_serializer<std::unordered_map<VertexIndex, Vector6d>> {
  static void to_json(json &j,
                      const std::unordered_map<VertexIndex, Vector6d> &value) {
    // calls the "to_json" method in T's namespace
  }

  static void from_json(const nlohmann::json &j,
                        std::unordered_map<VertexIndex, Vector6d> &m) {
    std::cout << "Entered." << std::endl;
    for (const auto &p : j) {
      m.emplace(p.at(0).template get<VertexIndex>(),
                p.at(1).template get<std::array<double, 6>>());
    }
  }
};
} // namespace nlohmann

class SimulationJob {
  //  const std::unordered_map<VertexIndex, VectorType> nodalForcedNormals;
  // json labels
  struct JSONLabels {
    inline static std::string meshFilename{"mesh filename"};
    inline static std::string forcedDisplacements{"forced displacements"};
    inline static std::string constrainedVertices{"fixed vertices"};
    inline static std::string nodalForces{"forces"};
    inline static std::string label{"label"};
    inline static std::string meshLabel{
        "meshLabel"}; // TODO: should be in the savePly function of the
                      // simulation mesh class
  } jsonLabels;

public:
  std::shared_ptr<SimulationMesh> pMesh;
  std::string label{"empty_job"};
  std::unordered_map<VertexIndex, std::unordered_set<int>> constrainedVertices;
  std::unordered_map<VertexIndex, Vector6d> nodalExternalForces;
  std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
  SimulationJob(
      const std::shared_ptr<SimulationMesh> &m, const std::string &label,
      const std::unordered_map<VertexIndex, std::unordered_set<int>> &cv,
      const std::unordered_map<VertexIndex, Vector6d> &ef = {},
      const std::unordered_map<VertexIndex, Eigen::Vector3d> &fd = {})
      : pMesh(m), label(label), constrainedVertices(cv),
        nodalExternalForces(ef), nodalForcedDisplacements(fd) {}
  SimulationJob() {}
  SimulationJob(const std::string &jsonFilename) { load(jsonFilename); }

  bool isEmpty()
  {
      return constrainedVertices.empty() && nodalExternalForces.empty()
             && nodalForcedDisplacements.empty() && pMesh == nullptr;
  }
  void remap(const std::unordered_map<size_t, size_t> &thisToOtherViMap,
             SimulationJob &simulationJobMapped) const
  {
      assert(simulationJobMapped.pMesh->VN() != 0);
      std::unordered_map<VertexIndex, std::unordered_set<int>> reducedModelFixedVertices;
      for (auto fullModelFixedVertex : this->constrainedVertices) {
          reducedModelFixedVertices[thisToOtherViMap.at(fullModelFixedVertex.first)]
              = fullModelFixedVertex.second;
      }

      std::unordered_map<VertexIndex, Vector6d> reducedModelNodalForces;
      for (auto fullModelNodalForce : this->nodalExternalForces) {
          reducedModelNodalForces[thisToOtherViMap.at(fullModelNodalForce.first)]
              = fullModelNodalForce.second;
      }

      std::unordered_map<VertexIndex, Eigen::Vector3d> reducedNodalForcedDisplacements;
      for (auto fullForcedDisplacement : this->nodalForcedDisplacements) {
          reducedNodalForcedDisplacements[thisToOtherViMap.at(fullForcedDisplacement.first)]
              = fullForcedDisplacement.second;
      }

      simulationJobMapped.constrainedVertices = reducedModelFixedVertices;
      simulationJobMapped.nodalExternalForces = reducedModelNodalForces;
      simulationJobMapped.label = this->getLabel();
      simulationJobMapped.nodalForcedDisplacements = reducedNodalForcedDisplacements;
  }
  SimulationJob getCopy() const {
    SimulationJob jobCopy;
    jobCopy.pMesh = std::make_shared<SimulationMesh>();
    jobCopy.pMesh->copy(*pMesh);
    jobCopy.label = label;
    jobCopy.constrainedVertices = constrainedVertices;
    jobCopy.nodalExternalForces = nodalExternalForces;
    jobCopy.nodalForcedDisplacements = nodalForcedDisplacements;

    return jobCopy;
  }
  std::string getLabel() const { return label; }

  std::string toString() const {
    nlohmann::json json;
    if (!constrainedVertices.empty()) {
      json[jsonLabels.constrainedVertices] = constrainedVertices;
    }
    if (!nodalExternalForces.empty()) {
      std::unordered_map<VertexIndex, std::array<double, 6>> arrForces;
      for (const auto &f : nodalExternalForces) {
        arrForces[f.first] = f.second;
      }
      json[jsonLabels.nodalForces] = arrForces;
    }

    return json.dump();
  }

  void clear()
  {
      label = "empty_job";
      constrainedVertices.clear();
      nodalExternalForces.clear();
      nodalForcedDisplacements.clear();
      if (pMesh.use_count() == 1) {
          std::cout << "Job mesh is deleted" << std::endl;
      }
      pMesh.reset();
  }

  bool load(const std::string &jsonFilename, const bool &shouldLoadMesh = true)
  {
      std::cout << jsonFilename << std::endl;
      label = "empty_job";
      constrainedVertices.clear();
      nodalExternalForces.clear();
      nodalForcedDisplacements.clear();
      const bool beVerbose = false;
      if (std::filesystem::path(jsonFilename).extension() != ".json") {
          std::cerr << "A json file is expected as input. The given file has the "
                       "following extension:"
                    << std::filesystem::path(jsonFilename).extension() << std::endl;
          assert(false);
          return false;
      }

      if (!std::filesystem::exists(std::filesystem::path(jsonFilename))) {
          std::cerr << "The json file does not exist. Json file provided:" << jsonFilename
                    << std::endl;
          assert(false);
          return false;
      }

      if (beVerbose) {
          std::cout << "Loading json file:" << jsonFilename << std::endl;
      }
      nlohmann::json json;
      std::ifstream ifs(jsonFilename);
      ifs >> json;

      if (shouldLoadMesh) {
          pMesh.reset();
          pMesh = std::make_shared<SimulationMesh>();
          if (json.contains(jsonLabels.meshFilename)) {
              const std::string relativeFilepath = json[jsonLabels.meshFilename];
              const auto meshFilepath = std::filesystem::path(
                                            std::filesystem::path(jsonFilename).parent_path())
                                            .append(relativeFilepath);
              pMesh->load(meshFilepath.string());
              pMesh->setLabel(json[jsonLabels.meshLabel]); // FIXME: This should be exported using
                                                           // nanoply but nanoply might not be able
                                                           // to write a string(??)
          }

          if (json.contains(jsonLabels.meshLabel)) {
              pMesh->setLabel(json[jsonLabels.meshLabel]);
          }
      }

      if (json.contains(jsonLabels.constrainedVertices)) {
          constrainedVertices =
              //      auto conV =
              json[jsonLabels.constrainedVertices]
                  .get<std::unordered_map<VertexIndex, std::unordered_set<int>>>();
          if (beVerbose) {
              std::cout << "Loaded constrained vertices. Number of constrained "
                           "vertices found:"
                        << constrainedVertices.size() << std::endl;
          }
      }

      if (json.contains(jsonLabels.nodalForces)) {
          auto f = 
              json[jsonLabels.nodalForces].get<std::unordered_map<VertexIndex, std::array<double, 6>>>();
          for (const auto &force : f) {
              nodalExternalForces[force.first] = Vector6d(force.second);
          }
          if (beVerbose) {
              std::cout << "Loaded forces. Number of forces found:" << nodalExternalForces.size()
                        << std::endl;
          }
      }

      if (json.contains(jsonLabels.forcedDisplacements)) {
          //      auto conV =
          std::unordered_map<VertexIndex, std::array<double, 3>> forcedDisp
              = json[jsonLabels.forcedDisplacements]
                    .get<std::unordered_map<VertexIndex, std::array<double, 3>>>();

          for (const auto &fd : forcedDisp) {
              nodalForcedDisplacements[fd.first] = Eigen::Vector3d(fd.second[0],
                                                                   fd.second[1],
                                                                   fd.second[2]);
          }
          if (beVerbose) {
              std::cout << "Loaded forced displacements. Number of forced displaced"
                           "vertices found:"
                        << nodalForcedDisplacements.size() << std::endl;
          }
      }

      if (json.contains(jsonLabels.label)) {
          label = json[jsonLabels.label];
      }

      return true;
  }

  bool save(const std::string &folderDirectory) const
  {
      const std::filesystem::path pathFolderDirectory(folderDirectory);
      if (!std::filesystem::is_directory(pathFolderDirectory)) {
          std::cerr << "A folder directory is expected for saving the simulation "
                       "job. Exiting.."
                    << std::endl;
          return false;
      }

      bool returnValue = true;
      std::string jsonFilename(
          std::filesystem::path(pathFolderDirectory)
              //            .append(label + "_" + pMesh->getLabel() + "_simulationJob.json")
              .append("SimulationJob.json")
              .string());

      const std::string meshFilename = std::filesystem::absolute(
                                           std::filesystem::canonical(
                                               std::filesystem::path(pathFolderDirectory)))
                                           .append(pMesh->getLabel() + ".ply")
                                           .string();
      returnValue = pMesh->save(meshFilename);
      nlohmann::json json;
json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(meshFilename),std::filesystem::path(jsonFilename).parent_path()).string();
      json[jsonLabels.meshLabel]
          = pMesh->getLabel(); // FIXME: This should be exported using nanoply but
                               // nanoply might not be able to write a string(??)
      if (!constrainedVertices.empty()) {
          json[jsonLabels.constrainedVertices] = constrainedVertices;
      }
      if (!nodalExternalForces.empty()) {
          std::unordered_map<VertexIndex, std::array<double, 6>> arrForces;
          for (const auto &f : nodalExternalForces) {
              arrForces[f.first] = f.second;
          }
          json[jsonLabels.nodalForces] = arrForces;
      }
      if (!nodalForcedDisplacements.empty()) {
          std::unordered_map<VertexIndex, std::array<double, 3>> forcedDisp;
          for (const auto &fd : nodalForcedDisplacements) {
              forcedDisp[fd.first] = {fd.second[0], fd.second[1], fd.second[2]};
          }

          json[jsonLabels.forcedDisplacements] = forcedDisp;
      }

      if (!label.empty()) {
          json[jsonLabels.label] = label;
      }
      if (!pMesh->getLabel().empty()) {
          json[jsonLabels.meshLabel] = pMesh->getLabel();
      }
      std::ofstream jsonFile(jsonFilename);
      jsonFile << json;
      jsonFile.close();
      //    std::cout << "Saved simulation job as:" << jsonFilename << std::endl;

      return returnValue;
  }
#ifdef POLYSCOPE_DEFINED
  void registerForDrawing(const std::string &meshLabel, const bool &shouldEnable = false) const
  {
      PolyscopeInterface::init();
      if (meshLabel.empty()) {
          assert(false);
          std::cerr << "Expects a mesh label on which to draw the simulation job." << std::endl;
          return;
      }

      if (!polyscope::hasCurveNetwork(meshLabel)) {
          assert(false);
          std::cerr << "Expects mesh already being registered to draw the "
                       "simulation job. No struct named "
                           + meshLabel
                    << std::endl;
          return;
      }
      std::vector<std::array<double, 3>> nodeColors(pMesh->VN());
      for (auto fixedVertex : constrainedVertices) {
          const bool hasRotationalDoFConstrained = fixedVertex.second.contains(3)
                                                   || fixedVertex.second.contains(4)
                                                   || fixedVertex.second.contains(5);
          const bool hasTranslationalDoFConstrained = fixedVertex.second.contains(0)
                                                      || fixedVertex.second.contains(1)
                                                      || fixedVertex.second.contains(2);
          if (hasTranslationalDoFConstrained && !hasRotationalDoFConstrained) {
              nodeColors[fixedVertex.first] = {0, 0, 1};
          } else if (!hasTranslationalDoFConstrained && hasRotationalDoFConstrained) {
              nodeColors[fixedVertex.first] = {0, 1, 0};
          } else {
              nodeColors[fixedVertex.first] = {0, 1, 1};
          }
      }
      if (!nodalForcedDisplacements.empty()) {
          for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair : nodalForcedDisplacements) {
              const VertexIndex vi = viDisplPair.first;
              nodeColors[vi][0] += 1;
              nodeColors[vi][0] /= 2;
              nodeColors[vi][1] += 0;
              nodeColors[vi][1] /= 2;
              nodeColors[vi][2] += 0;
              nodeColors[vi][2] /= 2;
          }
      }
      std::for_each(nodeColors.begin(), nodeColors.end(), [](std::array<double, 3> &color) {
          const double norm = sqrt(std::pow(color[0], 2) + std::pow(color[1], 2)
                                   + std::pow(color[2], 2));
          if (norm > std::pow(10, -7)) {
              color[0] /= norm;
              color[1] /= norm;
              color[2] /= norm;
          }
      });

      if (!nodeColors.empty()) {
          polyscope::getCurveNetwork(meshLabel)
              ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors)
              ->setEnabled(shouldEnable);
      }

      // per node external forces
      std::vector<std::array<double, 3>> externalForces(pMesh->VN());
      for (const auto &forcePair : nodalExternalForces) {
          auto index = forcePair.first;
          auto force = forcePair.second;
          externalForces[index] = {force[0], force[1], force[2]};
      }

      if (!externalForces.empty()) {
          polyscope::getCurveNetwork(meshLabel)
              ->addNodeVectorQuantity("External force_" + label, externalForces)
              ->setEnabled(shouldEnable);
      }
      // per node external moments
      bool hasExternalMoments = false;
      std::vector<std::array<double, 3>> externalMoments(pMesh->VN());
      for (const auto &forcePair : nodalExternalForces) {
          auto index = forcePair.first;
          const Vector6d &load = forcePair.second;
          if (load.getRotation().norm() != 0) {
              hasExternalMoments = true;
          }
          externalMoments[index] = {load[3], load[4], load[5]};
      }

      if (hasExternalMoments) {
          polyscope::getCurveNetwork(meshLabel)
              ->addNodeVectorQuantity("External moment_" + label, externalMoments)
              ->setEnabled(shouldEnable);
      }
  }
  void unregister(const std::string &meshLabel)
  {
      if (polyscope::getCurveNetwork(meshLabel) == nullptr) {
          return;
      }
      if (!nodalExternalForces.empty()) {
          polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force_" + label);
      }
      if (!constrainedVertices.empty()) {
          polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label);
      }
  }
#endif // POLYSCOPE_DEFINED
};
struct SimulationResults
{
    /*TODO: remove rotationalDisplacementQuaternion since the last three components of the displacments
 * vector contains the same info using euler angles
*/
    bool converged{false};
    std::shared_ptr<SimulationJob> job;
    SimulationHistory history;
    std::vector<Vector6d> debug_drmDisplacements;
    std::vector<Eigen::Quaternion<double>> debug_q_f1;     //per vertex
    std::vector<Eigen::Quaternion<double>> debug_q_normal; //per vertex
    std::vector<Eigen::Quaternion<double>> debug_q_nr;     //per vertex
    std::vector<Vector6d> displacements;
    std::vector<Eigen::Quaternion<double>> rotationalDisplacementQuaternion; //per vertex
    double internalPotentialEnergy{0};
    double executionTime{0};
    std::string labelPrefix{"deformed"};
    inline static char deliminator{' '};
    SimulationResults() { job = std::make_shared<SimulationJob>(); }

    std::vector<VectorType> getTranslationalDisplacements() const
    {
        std::vector<VectorType> translationalDisplacements(displacements.size());
        std::transform(displacements.begin(),
                       displacements.end(),
                       translationalDisplacements.begin(),
                       [&](const Vector6d &d) { return VectorType(d[0], d[1], d[2]); });

        return translationalDisplacements;
    }

    void setLabelPrefix(const std::string &lp) { labelPrefix += deliminator + lp; }
    std::string getLabel() const
    {
        return labelPrefix + deliminator + job->pMesh->getLabel() + deliminator + job->getLabel();
    }

    bool saveDeformedModel(const std::string &outputFolder = std::string())
    {
        VCGEdgeMesh m;
        vcg::tri::Append<VCGEdgeMesh, SimulationMesh>::MeshCopy(m, *job->pMesh);

        for (int vi = 0; vi < m.VN(); vi++) {
            m.vert[vi].P() = m.vert[vi].P()
                             + CoordType(displacements[vi][0],
                                         displacements[vi][1],
                                         displacements[vi][2]);
        }
        return m.save(std::filesystem::path(outputFolder).append(getLabel() + ".ply").string());
    }

    void save(const std::string &outputFolder = std::string())
    {
        const std::filesystem::path outputFolderPath = outputFolder.empty()
                                                           ? std::filesystem::current_path()
                                                           : std::filesystem::path(outputFolder);
        std::filesystem::path simulationJobOutputFolderPath
            = std::filesystem::path(outputFolderPath).append("SimulationJob");
        std::filesystem::create_directories(simulationJobOutputFolderPath);
        job->save(simulationJobOutputFolderPath.string());
        const std::string filename(getLabel() + "_displacements.eigenBin");

        Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements);
        const std::filesystem::path resultsFolderPath(
            std::filesystem::path(outputFolder).append("Results"));
        std::filesystem::create_directories(resultsFolderPath);
        Eigen::write_binary(std::filesystem::path(resultsFolderPath).append(filename).string(), m);

        saveDeformedModel(resultsFolderPath.string());
    }

    // The comparison of the results happens comparing the 6-dof nodal
    // displacements
    bool isEqual(const Eigen::MatrixXd &nodalDisplacements, double &error)
    {
        assert(nodalDisplacements.cols() == 6);
        Eigen::MatrixXd eigenDisplacements = Utilities::toEigenMatrix(this->displacements);
        const double errorNorm = (eigenDisplacements - nodalDisplacements).norm();
        error = errorNorm;
        return errorNorm < 1e-10;
        //    return eigenDisplacements.isApprox(nodalDisplacements);
    }

    void load(const std::filesystem::path &loadFromPath, const std::filesystem::path &loadJobFrom)
    {
        //load job
        job->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json").string());
        //Use the first .eigenBin file for loading the displacements
        for (auto const &entry : std::filesystem::recursive_directory_iterator(loadFromPath)) {
            if (filesystem::is_regular_file(entry) && entry.path().extension() == ".eigenBin") {
                Eigen::MatrixXd displacements_eigen;
                Eigen::read_binary(entry.path().string(), displacements_eigen);
                displacements = Utilities::fromEigenMatrix(displacements_eigen);
                break;
            }
        }

        rotationalDisplacementQuaternion.resize(job->pMesh->VN());
        for (int vi = 0; vi < job->pMesh->VN(); vi++) {
            rotationalDisplacementQuaternion[vi] = Eigen::AngleAxisd(displacements[vi][3],
                                                                     Eigen::Vector3d::UnitX())
                                                   * Eigen::AngleAxisd(displacements[vi][4],
                                                                       Eigen::Vector3d::UnitY())
                                                   * Eigen::AngleAxisd(displacements[vi][5],
                                                                       Eigen::Vector3d::UnitZ());
        }
    }

#ifdef POLYSCOPE_DEFINED
  void unregister() const {
    if (!polyscope::hasCurveNetwork(getLabel())) {
      std::cerr << "No curve network registered with a name: " << getLabel()
                << std::endl;
      std::cerr << "Nothing to remove." << std::endl;
      return;
    }
    job->unregister(getLabel());
    polyscope::removeCurveNetwork(getLabel());
  }
  void registerForDrawing(const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
                          const bool &shouldEnable = true,
                          const bool &shouldShowFrames = false) const
  {
      PolyscopeInterface::init();
      const std::shared_ptr<SimulationMesh> &mesh = job->pMesh;
      polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh;
      if (!polyscope::hasStructure(getLabel())) {
          polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(),
                                                                             mesh->getEigenVertices(),
                                                                             mesh->getEigenEdges());
      }
      polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable);
      polyscopeHandle_deformedEdmeMesh->setRadius(0.002, true);
      if (desiredColor.has_value()) {
          const glm::vec3 desiredColor_glm(desiredColor.value()[0],
                                           desiredColor.value()[1],
                                           desiredColor.value()[2]);
          polyscopeHandle_deformedEdmeMesh->setColor(desiredColor_glm);
      }
      Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3);
      Eigen::MatrixX3d framesX(mesh->VN(), 3);
      Eigen::MatrixX3d framesY(mesh->VN(), 3);
      Eigen::MatrixX3d framesZ(mesh->VN(), 3);
      for (VertexIndex vi = 0; vi < mesh->VN(); vi++) {
          const Vector6d &nodalDisplacement = displacements[vi];
          nodalDisplacements.row(vi) = Eigen::Vector3d(nodalDisplacement[0],
                                                       nodalDisplacement[1],
                                                       nodalDisplacement[2]);
          auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1);
          auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0);
          auto deformedFrameX = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(1, 0, 0);
          framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0], deformedFrameX[1], deformedFrameX[2]);
          framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0], deformedFrameY[1], deformedFrameY[2]);
          framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0], deformedNormal[1], deformedNormal[2]);
      }
      polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices()
                                                            + nodalDisplacements);

      //if(showFramesOn.contains(vi)){
      auto polyscopeHandle_frameX = polyscopeHandle_deformedEdmeMesh
                                        ->addNodeVectorQuantity("FrameX", framesX);
      polyscopeHandle_frameX->setVectorLengthScale(0.01);
      polyscopeHandle_frameX->setVectorRadius(0.01);
      polyscopeHandle_frameX->setVectorColor(
          /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
      auto polyscopeHandle_frameY = polyscopeHandle_deformedEdmeMesh
                                        ->addNodeVectorQuantity("FrameY", framesY);
      polyscopeHandle_frameY->setVectorLengthScale(0.01);
      polyscopeHandle_frameY->setVectorRadius(0.01);
      polyscopeHandle_frameY->setVectorColor(
          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
      auto polyscopeHandle_frameZ = polyscopeHandle_deformedEdmeMesh
                                        ->addNodeVectorQuantity("FrameZ", framesZ);
      polyscopeHandle_frameZ->setVectorLengthScale(0.01);
      polyscopeHandle_frameZ->setVectorRadius(0.01);
      polyscopeHandle_frameZ->setVectorColor(
          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
      //}
      job->registerForDrawing(getLabel());
      static bool wasExecuted = false;
      if (!wasExecuted) {
          std::function<void()> callback = [&]() {
              static bool showFrames = shouldShowFrames;

              if (ImGui::Checkbox("Show Frames", &showFrames) && showFrames) {
                  polyscopeHandle_frameX->setEnabled(showFrames);
                  polyscopeHandle_frameY->setEnabled(showFrames);
                  polyscopeHandle_frameZ->setEnabled(showFrames);
              }
          };

          PolyscopeInterface::addUserCallback(callback);
          wasExecuted = true;
      }
  }
#endif
};

#endif // SIMULATIONHISTORY_HPP