#ifndef SIMULATIONSTRUCTS_HPP
#define SIMULATIONSTRUCTS_HPP
#include <fstream>

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();
}
//const static IOFormat CSVFormat(StreamPrecision, DontAlignCols, ", ", "\n");
//template<class Matrix>
//void writeToCSV(const std::string &filename, Matrix &matrix)
//{
//    ofstream file(filename.c_str());
//    file << matrix.format(CSVFormat);
//}
} // 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:
    inline static std::string jsonDefaultFileName = "SimulationJob.json";
    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> &sourceToDestinationViMap,
             SimulationJob &destination_simulationJob) const
  {
      std::unordered_map<VertexIndex, std::unordered_set<int>> destination_fixedVertices;
      for (const auto &source_fixedVertex : this->constrainedVertices) {
          destination_fixedVertices[sourceToDestinationViMap.at(source_fixedVertex.first)]
              = source_fixedVertex.second;
      }

      std::unordered_map<VertexIndex, Vector6d> destination_nodalForces;
      for (const auto &source_nodalForces : this->nodalExternalForces) {
          destination_nodalForces[sourceToDestinationViMap.at(source_nodalForces.first)]
              = source_nodalForces.second;
      }

      std::unordered_map<VertexIndex, Eigen::Vector3d> destination_forcedDisplacements;
      for (const auto &source_forcedDisplacements : this->nodalForcedDisplacements) {
          destination_forcedDisplacements[sourceToDestinationViMap.at(
              source_forcedDisplacements.first)]
              = source_forcedDisplacements.second;
      }

      destination_simulationJob.constrainedVertices = destination_fixedVertices;
      destination_simulationJob.nodalExternalForces = destination_nodalForces;
      destination_simulationJob.label = this->getLabel();
      destination_simulationJob.nodalForcedDisplacements = destination_forcedDisplacements;
  }
  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)
  {
      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 (const 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 (const 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::CurveNetworkNodeVectorQuantity *externalForcesVectors
              = polyscope::getCurveNetwork(meshLabel)->addNodeVectorQuantity("External force_"
                                                                                 + label,
                                                                             externalForces);

          const std::array<double, 3> color_loads{1.0, 0, 0};
          externalForcesVectors->setVectorColor(
              glm::vec3(color_loads[0], color_loads[1], color_loads[2]));
          externalForcesVectors->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);
      }

      // per node external moments
      bool hasExternalMoments = false;
      for (const auto &forcePair : nodalExternalForces) {
          const Vector6d &load = forcePair.second;
          if (load.getRotation().norm() != 0) {
              hasExternalMoments = true;
              break;
          }
      }
      if (hasExternalMoments) {
          polyscope::getCurveNetwork(meshLabel)->removeQuantity("External moment_" + 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
*/
    inline const static std::string defaultJsonFilename{"SimulationResults.json"};
    bool converged{false};
    std::shared_ptr<SimulationJob> pJob;
    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() { pJob = 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 + pJob->pMesh->getLabel() + deliminator + pJob->getLabel();
    }

    bool saveDeformedModel(const std::string &outputFolder = std::string())
    {
        VCGEdgeMesh m;
        vcg::tri::Append<VCGEdgeMesh, SimulationMesh>::MeshCopy(m, *pJob->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::cout << "Saving results to:" << outputFolderPath << std::endl;
        std::filesystem::path simulationJobOutputFolderPath
            = std::filesystem::path(outputFolderPath).append("SimulationJob");
        std::filesystem::create_directories(simulationJobOutputFolderPath);
        pJob->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);

        nlohmann::json json;

        json[GET_VARIABLE_NAME(internalPotentialEnergy)] = internalPotentialEnergy;

        std::filesystem::path jsonFilePath(
            std::filesystem::path(resultsFolderPath).append(defaultJsonFilename));
        std::ofstream jsonFile(jsonFilePath.string());
        jsonFile << json;
        jsonFile.close();

        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)
    {
        pJob->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json").string());
        load(loadFromPath);
    }
    void load(const std::filesystem::path &loadFromPath, const std::shared_ptr<SimulationJob> &pJob)
    {
        this->pJob = pJob;
        load(loadFromPath);
    }
    static double computeDistance(
        const SimulationResults &resultsA,
        const SimulationResults &resultsB,
        const std::unordered_map<VertexIndex, VertexIndex> &resultsAToResultsBViMap)
    {
        double distance = 0;
        for (std::pair<int, int> resultsAToResultsBViPair : resultsAToResultsBViMap) {
            const double vertexToVertexDistance
                = (resultsA.displacements[resultsAToResultsBViPair.first].getTranslation()
                   - resultsB.displacements[resultsAToResultsBViPair.second].getTranslation())
                      .norm();
            distance += vertexToVertexDistance;
        }
        return distance;
    }

    double computeDistance(const SimulationResults &other,
                           const std::unordered_map<VertexIndex, VertexIndex> &thisToOtherViMap) const
    {
        return computeDistance(*this, other, thisToOtherViMap);
    }

#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;
        }
        pJob->unregister(getLabel());
        polyscope::removeCurveNetwork(getLabel());
    }
    polyscope::CurveNetwork *registerForDrawing(
        const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
        const bool &shouldEnable = true,
        const double &desiredRadius = 0.001,
        //      const double &desiredRadius = 0.0001,
        const bool &shouldShowFrames = false) const
    {
        PolyscopeInterface::init();
        const std::shared_ptr<SimulationMesh> &mesh = pJob->pMesh;
        polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh;
        if (!polyscope::hasStructure(getLabel())) {
            const auto verts = mesh->getEigenVertices();
            const auto edges = mesh->getEigenEdges();
            polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(),
                                                                               verts,
                                                                               edges);
        } else {
            polyscopeHandle_deformedEdmeMesh = polyscope::getCurveNetwork(getLabel());
        }
        polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable);
        polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, 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);

        Eigen::MatrixX3d framesX_initial(mesh->VN(), 3);
        Eigen::MatrixX3d framesY_initial(mesh->VN(), 3);
        Eigen::MatrixX3d framesZ_initial(mesh->VN(), 3);

        //      std::unordered_set<int> interfaceNodes{1, 3, 5, 7, 9, 11};
        //                  std::unordered_set<int> interfaceNodes{3, 7, 11, 15, 19, 23};
        //      std::unordered_set<int> interfaceNodes{};
        for (VertexIndex vi = 0; vi < mesh->VN(); vi++) {
            const Vector6d &nodalDisplacement = displacements[vi];
            nodalDisplacements.row(vi) = Eigen::Vector3d(nodalDisplacement[0],
                                                         nodalDisplacement[1],
                                                         nodalDisplacement[2]);
            //           Eigen::Quaternion<double> Rx(Eigen::AngleAxis(nodalDisplacement[2],Eigen::Vector3d(1, 0, 0)));
            //           Eigen::Quaternion<double> Ry(Eigen::AngleAxis(nodalDisplacement[4],Eigen::Vector3d(0, 1, 0)));
            //           Eigen::Quaternion<double> Rz(Eigen::AngleAxis(nodalDisplacement[5],Eigen::Vector3d(0, 0, 1)));
            //           Eigen::Quaternion<double> R=Rx*Ry*Rz;
            //          if (interfaceNodes.contains(vi)) {
            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]);
            framesX_initial.row(vi) = Eigen::Vector3d(1, 0, 0);
            framesY_initial.row(vi) = Eigen::Vector3d(0, 1, 0);
            framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 1);
            //          } else {
            //              framesX.row(vi) = Eigen::Vector3d(0, 0, 0);
            //              framesY.row(vi) = Eigen::Vector3d(0, 0, 0);
            //              framesZ.row(vi) = Eigen::Vector3d(0, 0, 0);
            //              framesX_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
            //              framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
            //              framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
            //          }
        }
        polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices()
                                                              + nodalDisplacements);

        const double frameRadius_default = 0.035;
        const double frameLength_default = 0.035;
        const bool shouldEnable_default = true;
        //if(showFramesOn.contains(vi)){
        auto polyscopeHandle_frameX = polyscopeHandle_deformedEdmeMesh
                                          ->addNodeVectorQuantity("FrameX", framesX);
        polyscopeHandle_frameX->setVectorLengthScale(frameLength_default);
        polyscopeHandle_frameX->setVectorRadius(frameRadius_default);
        polyscopeHandle_frameX->setVectorColor(
            /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
        auto polyscopeHandle_frameY = polyscopeHandle_deformedEdmeMesh
                                          ->addNodeVectorQuantity("FrameY", framesY);
        polyscopeHandle_frameY->setVectorLengthScale(frameLength_default);
        polyscopeHandle_frameY->setVectorRadius(frameRadius_default);
        polyscopeHandle_frameY->setVectorColor(
            /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
        auto polyscopeHandle_frameZ = polyscopeHandle_deformedEdmeMesh
                                          ->addNodeVectorQuantity("FrameZ", framesZ);
        polyscopeHandle_frameZ->setVectorLengthScale(frameLength_default);
        polyscopeHandle_frameZ->setVectorRadius(frameRadius_default);
        polyscopeHandle_frameZ->setVectorColor(
            /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));

        auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel());
        if (!polyscopeHandle_initialMesh) {
            polyscopeHandle_initialMesh = mesh->registerForDrawing();
        }

        //      auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh
        //                                                ->addNodeVectorQuantity("FrameX", framesX_initial);
        //      polyscopeHandle_frameX_initial->setVectorLengthScale(frameLength_default);
        //      polyscopeHandle_frameX_initial->setVectorRadius(frameRadius_default);
        //      polyscopeHandle_frameX_initial->setVectorColor(
        //          /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
        //      auto polyscopeHandle_frameY_initial = polyscopeHandle_initialMesh
        //                                                ->addNodeVectorQuantity("FrameY", framesY_initial);
        //      polyscopeHandle_frameY_initial->setVectorLengthScale(frameLength_default);
        //      polyscopeHandle_frameY_initial->setVectorRadius(frameRadius_default);
        //      polyscopeHandle_frameY_initial->setVectorColor(
        //          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
        //      auto polyscopeHandle_frameZ_initial = polyscopeHandle_initialMesh
        //                                                ->addNodeVectorQuantity("FrameZ", framesZ_initial);
        //      polyscopeHandle_frameZ_initial->setVectorLengthScale(frameLength_default);
        //      polyscopeHandle_frameZ_initial->setVectorRadius(frameRadius_default);
        //      polyscopeHandle_frameZ_initial->setVectorColor(
        //          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
        //      //}
        pJob->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;
        //        }
        return polyscopeHandle_deformedEdmeMesh;
    }
#endif
private:
    bool load(const std::filesystem::path &loadFromPath)
    {
        converged = true; //assuming it has converged
        assert(pJob != nullptr);
        //load job
        //Use the first .eigenBin file for loading the displacements
        for (auto const &entry : std::filesystem::recursive_directory_iterator(loadFromPath)) {
            if (std::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(pJob->pMesh->VN());
        for (int vi = 0; vi < pJob->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());
        }

        const std::filesystem::path jsonFilepath = std::filesystem::path(loadFromPath)
                                                       .append(defaultJsonFilename);
        if (!std::filesystem::exists(jsonFilepath)) {
            std::cerr << "Simulation results could not be loaded because filepath does "
                         "not exist:"
                      << jsonFilepath << std::endl;
            return false;
        }
        std::ifstream ifs(jsonFilepath);
        nlohmann::json json;
        ifs >> json;
        if (json.contains(GET_VARIABLE_NAME(internalPotentialEnergy))) {
            internalPotentialEnergy = json.at(GET_VARIABLE_NAME(internalPotentialEnergy));
        }
        return true;
    }
};

#endif // SIMULATIONHISTORY_HPP