diff --git a/beam.hpp b/beam.hpp
index 4cffe70..dc485a4 100644
--- a/beam.hpp
+++ b/beam.hpp
@@ -3,6 +3,7 @@
 #include <assert.h>
 #include <cmath>
 #include <iostream>
+#include <string>
 
 struct RectangularBeamDimensions {
   inline static std::string name{"Rectangular"};
diff --git a/beamformfinder.cpp b/beamformfinder.cpp
index b7bb397..87e2a2a 100644
--- a/beamformfinder.cpp
+++ b/beamformfinder.cpp
@@ -19,7 +19,7 @@ void FormFinder::runUnitTests() {
   //  const size_t spanGridSize = 11;
   //  mesh.createSpanGrid(spanGridSize);
   beam.loadPly(
-      std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply"));
+      std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply").string());
   std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
   fixedVertices[0] = std::unordered_set<DoFType>{0, 1, 2, 3};
   fixedVertices[beam.VN() - 1] = std::unordered_set<DoFType>{1, 2};
@@ -50,7 +50,7 @@ void FormFinder::runUnitTests() {
   simpleBeam_simulationResults.save();
   const std::string simpleBeamGroundOfTruthBinaryFilename =
       std::filesystem::path(groundOfTruthFolder)
-          .append("SimpleBeam_displacements.eigenBin");
+          .append("SimpleBeam_displacements.eigenBin").string();
   assert(std::filesystem::exists(
       std::filesystem::path(simpleBeamGroundOfTruthBinaryFilename)));
   Eigen::MatrixXd simpleBeam_groundOfTruthDisplacements;
@@ -211,7 +211,7 @@ void FormFinder::reset() {
   history.clear();
   constrainedVertices.clear();
   rigidSupports.clear();
-  pMesh.release();
+  pMesh.reset();
   plotYValues.clear();
   plotHandle.reset();
   checkedForMaximumMoment = false;
@@ -219,6 +219,7 @@ void FormFinder::reset() {
   externalMomentsNorm = 0;
   mSettings.drawingStep = 1;
   Dt = mSettings.Dtini;
+  numOfDampings=0;
 }
 
 VectorType FormFinder::computeDisplacementDifferenceDerivative(
@@ -793,8 +794,8 @@ void FormFinder::updateResidualForcesOnTheFly(
           std::vector<std::pair<int, Vector6d>>(4, {-1, Vector6d()}));
   //  omp_lock_t writelock;
   //  omp_init_lock(&writelock);
-#pragma omp parallel for schedule(static) num_threads(8)
-  for (size_t ei = 0; ei < pMesh->EN(); ei++) {
+//#pragma omp parallel for //schedule(static) num_threads(8)
+  for (int ei = 0; ei < pMesh->EN(); ei++) {
     const EdgeType &e = pMesh->edge[ei];
     const SimulationMesh::VertexType &ev_j = *e.cV(0);
     const SimulationMesh::VertexType &ev_jplus1 = *e.cV(1);
@@ -1790,14 +1791,15 @@ FormFinder::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
   //    }
   //  }
 
+  pMesh.reset();
   pMesh = std::make_unique<SimulationMesh>(*pJob->pMesh);
-  if (mSettings.beVerbose) {
+  if (mSettings.beVerbose  ) {
     std::cout << "Executing simulation for mesh with:" << pMesh->VN()
               << " nodes and " << pMesh->EN() << " elements." << std::endl;
   }
   computeRigidSupports();
 
-  if (mSettings.shouldDraw) {
+  if (mSettings.shouldDraw ) {
     initPolyscope();
     polyscope::registerCurveNetwork(
         meshPolyscopeLabel, pMesh->getEigenVertices(), pMesh->getEigenEdges());
@@ -1969,7 +1971,7 @@ mesh->currentTotalPotentialEnergykN*/
   } else if (std::isnan(pMesh->currentTotalKineticEnergy)) {
     results.converged = false;
 
-  } else if (mSettings.beVerbose) {
+  } else if (mSettings.beVerbose ) {
     auto t2 = std::chrono::high_resolution_clock::now();
     double simulationDuration =
         std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
diff --git a/csvfile.hpp b/csvfile.hpp
index a53f575..9bd0511 100644
--- a/csvfile.hpp
+++ b/csvfile.hpp
@@ -23,8 +23,16 @@ public:
       : fs_(), is_first_(true), separator_(separator), escape_seq_("\""),
         special_chars_("\"") {
     fs_.exceptions(std::ios::failbit | std::ios::badbit);
-    overwrite ? fs_.open(filename, std::ios::trunc)
-              : fs_.open(filename, std::ios::app);
+    if (!std::filesystem::exists(std::filesystem::path("../OptimizationResults")
+        .append("statistics.csv")
+    )) {
+        std::ofstream outfile(std::filesystem::path("../OptimizationResults")
+        .append("statistics.csv")
+    .string());
+        outfile.close();
+}
+    overwrite  ? fs_.open(filename, std::ios::trunc)
+                   : fs_.open(filename, std::ios::app);
   }
 
   ~csvfile() {
diff --git a/edgemesh.cpp b/edgemesh.cpp
index c925f70..a32257a 100644
--- a/edgemesh.cpp
+++ b/edgemesh.cpp
@@ -183,21 +183,11 @@ bool VCGEdgeMesh::loadPly(const std::string plyFilename) {
   assert(std::filesystem::exists(usedPath));
   this->Clear();
   const bool useDefaultImporter = false;
-  if (useDefaultImporter) {
-    if (!loadUsingDefaultLoader(usedPath)) {
-      return false;
-    }
-
-    eigenEdgeNormals.resize(EN(), 3);
-    for (int i = 0; i < EN(); i++) {
-      eigenEdgeNormals.row(i) = Eigen::Vector3d(0, 1, 0);
-    }
-  } else {
     if (!loadUsingNanoply(usedPath)) {
       std::cerr << "Error: Unable to open " + usedPath << std::endl;
       return false;
     }
-  }
+
   getEdges(eigenEdges);
   getVertices(eigenVertices);
   vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
diff --git a/edgemesh.hpp b/edgemesh.hpp
index 613deed..e420a4c 100644
--- a/edgemesh.hpp
+++ b/edgemesh.hpp
@@ -53,7 +53,6 @@ public:
 
   Eigen::MatrixX3d getNormals() const;
 
-  bool loadUsingDefaultLoader(const std::string &plyFilename);
   bool hasProperty(const std::vector<nanoply::PlyProperty> &v,
                    const std::string &propertyName);
 
diff --git a/elementalmesh.cpp b/elementalmesh.cpp
index edce87b..1647c19 100644
--- a/elementalmesh.cpp
+++ b/elementalmesh.cpp
@@ -38,6 +38,12 @@ SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) {
   eigenVertices = mesh.getEigenVertices();
 }
 
+SimulationMesh::~SimulationMesh()
+{
+    vcg::tri::Allocator<VCGEdgeMesh>::DeletePerEdgeAttribute<Element>(*this, elements);
+    vcg::tri::Allocator<VCGEdgeMesh>::DeletePerVertexAttribute<Node>(*this,nodes);
+}
+
 SimulationMesh::SimulationMesh(FlatPattern &pattern) {
   vcg::tri::MeshAssert<FlatPattern>::VertexNormalNormalized(pattern);
 
diff --git a/elementalmesh.hpp b/elementalmesh.hpp
index cd88e4e..e04c495 100644
--- a/elementalmesh.hpp
+++ b/elementalmesh.hpp
@@ -23,6 +23,7 @@ private:
 public:
   PerEdgeAttributeHandle<Element> elements;
   PerVertexAttributeHandle<Node> nodes;
+  ~SimulationMesh();
   SimulationMesh(FlatPattern &pattern);
   SimulationMesh(ConstVCGEdgeMesh &edgeMesh);
   SimulationMesh(SimulationMesh &elementalMesh);
diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp
index fc24434..02ea8e6 100644
--- a/simulationhistoryplotter.hpp
+++ b/simulationhistoryplotter.hpp
@@ -75,7 +75,7 @@ struct SimulationResultsReporter {
 
       createPlots(simulationResult.history, simulationResultPath.string(),
                   graphSuffix);
-      writeStatistics(simulationResult, simulationResultPath);
+      writeStatistics(simulationResult, simulationResultPath.string());
     }
   }
   static void createPlot(const std::string &xLabel, const std::string &yLabel,
diff --git a/simulationresult.hpp b/simulationresult.hpp
index ae43a04..1e60842 100644
--- a/simulationresult.hpp
+++ b/simulationresult.hpp
@@ -84,14 +84,14 @@ public:
   std::string toString() const {
     nlohmann::json json;
     if (!constrainedVertices.empty()) {
-      json[jsonLabels.constrainedVertices] = constrainedVertices;
+      json[jsonLabel_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;
+      json[jsonLabel_nodalForces] = arrForces;
     }
 
     return json.dump();
@@ -143,7 +143,6 @@ public:
                     }
                   });
 
-    auto cn = polyscope::getCurveNetwork(meshLabel);
     if (!nodeColors.empty()) {
       polyscope::getCurveNetwork(meshLabel)
           ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors)
@@ -194,8 +193,7 @@ public:
     if (json.contains(jsonLabels.constrainedVertices)) {
       constrainedVertices =
           //      auto conV =
-          std::unordered_map<VertexIndex, std::unordered_set<int>>(
-              json[jsonLabels.constrainedVertices]);
+              json[jsonLabels.constrainedVertices].get<std::unordered_map<VertexIndex, std::unordered_set<int>>>();
       std::cout << "Loaded constrained vertices. Number of constrained "
                    "vertices found:"
                 << constrainedVertices.size() << std::endl;
@@ -236,7 +234,7 @@ public:
         std::filesystem::absolute(
             std::filesystem::canonical(
                 std::filesystem::path(pathFolderDirectory)))
-            .append(pMesh->getLabel() + ".ply");
+            .append(pMesh->getLabel() + ".ply").string();
     returnValue = pMesh->savePly(meshFilename);
     nlohmann::json json;
     json[jsonLabels.meshFilename] = meshFilename;
@@ -260,7 +258,7 @@ public:
 
     std::string jsonFilename(
         std::filesystem::path(pathFolderDirectory)
-            .append(label + "_" + pMesh->getLabel() + "_simulationJob.json"));
+            .append(label + "_" + pMesh->getLabel() + "_simulationJob.json").string());
     std::ofstream jsonFile(jsonFilename);
     jsonFile << json;
     std::cout << "Saved simulation job as:" << jsonFilename << std::endl;
diff --git a/simulationresult.hpp.orig b/simulationresult.hpp.orig
new file mode 100644
index 0000000..4a65a58
--- /dev/null
+++ b/simulationresult.hpp.orig
@@ -0,0 +1,397 @@
+#ifndef SIMULATIONHISTORY_HPP
+#define SIMULATIONHISTORY_HPP
+
+#include "elementalmesh.hpp"
+#include "nlohmann/json.hpp"
+
+struct SimulationHistory {
+  SimulationHistory() {}
+
+  size_t numberOfSteps{0};
+  std::string label;
+  std::vector<double> residualForces;
+  std::vector<double> kineticEnergy;
+  std::vector<double> potentialEnergies;
+  std::vector<size_t> redMarks;
+  std::vector<double> greenMarks;
+
+  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(log(mesh.currentTotalKineticEnergy));
+    //    potentialEnergy.push_back(mesh.totalPotentialEnergykN);
+    residualForces.push_back(mesh.totalResidualForcesNorm);
+  }
+
+  void clear() {
+    residualForces.clear();
+    kineticEnergy.clear();
+    potentialEnergies.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 constrainedVertices{"fixed vertices"};
+    inline static std::string nodalForces{"forces"};
+    inline static std::string label{"label"};
+    inline static std::string meshLabel{
+        "label"}; // 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); }
+  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 registerForDrawing(const std::string &meshLabel) const {
+    initPolyscope();
+    if (meshLabel.empty()) {
+      assert(false);
+      std::cerr << "Expects a mesh label on which to draw the simulation job."
+                << std::endl;
+      return;
+    }
+
+    auto structs = polyscope::state::structures;
+    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) {
+      nodeColors[fixedVertex.first] = {0, 0, 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;
+                    }
+                  });
+
+    auto cn = polyscope::getCurveNetwork(meshLabel);
+    if (!nodeColors.empty()) {
+      polyscope::getCurveNetwork(meshLabel)
+          ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors)
+          ->setEnabled(false);
+    }
+
+    // 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(false);
+    }
+  }
+
+  bool load(const std::string &jsonFilename) {
+    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;
+    }
+
+    std::cout << "Loading json file:" << jsonFilename << std::endl;
+    nlohmann::json json;
+    std::ifstream ifs(jsonFilename);
+    json << ifs;
+
+    pMesh = std::make_shared<SimulationMesh>();
+    if (json.contains(jsonLabels.meshFilename)) {
+      pMesh->loadPly(json[jsonLabels.meshFilename]);
+    }
+
+    if (json.contains(jsonLabels.constrainedVertices)) {
+      constrainedVertices =
+          //      auto conV =
+<<<<<<< HEAD
+          std::unordered_map<VertexIndex, std::unordered_set<int>>(
+              json[jsonLabels.constrainedVertices]);
+=======
+              json[jsonLabel_constrainedVertices].get<std::unordered_map<VertexIndex, std::unordered_set<int>>>();
+>>>>>>> 2599d472614f403f8249060272ce288033a919d4
+      std::cout << "Loaded constrained vertices. Number of constrained "
+                   "vertices found:"
+                << constrainedVertices.size() << std::endl;
+    }
+
+<<<<<<< HEAD
+    if (json.contains(jsonLabels.nodalForces)) {
+      auto f = std::unordered_map<VertexIndex, std::array<double, 6>>(
+          json[jsonLabels.nodalForces]);
+=======
+    if (json.contains(jsonLabel_nodalForces)) {
+      auto f (
+          json[jsonLabel_nodalForces].get<std::unordered_map<VertexIndex, std::array<double, 6>>>());
+>>>>>>> 2599d472614f403f8249060272ce288033a919d4
+      for (const auto &forces : f) {
+        nodalExternalForces[forces.first] = Vector6d(forces.second);
+      }
+      std::cout << "Loaded forces. Number of forces found:"
+                << nodalExternalForces.size() << std::endl;
+    }
+
+    if (json.contains(jsonLabels.label)) {
+      label = json[jsonLabels.label];
+    }
+
+    if (json.contains(jsonLabels.meshLabel)) {
+      pMesh->setLabel(json[jsonLabels.meshLabel]);
+    }
+
+    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;
+    const std::string meshFilename =
+        std::filesystem::absolute(
+            std::filesystem::canonical(
+                std::filesystem::path(pathFolderDirectory)))
+            .append(pMesh->getLabel() + ".ply").string();
+    returnValue = pMesh->savePly(meshFilename);
+    nlohmann::json json;
+    json[jsonLabels.meshFilename] = meshFilename;
+    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 (!label.empty()) {
+      json[jsonLabels.label] = label;
+    }
+    if (!pMesh->getLabel().empty()) {
+      json[jsonLabels.meshLabel] = pMesh->getLabel();
+    }
+
+    std::string jsonFilename(
+        std::filesystem::path(pathFolderDirectory)
+<<<<<<< HEAD
+            .append(label + "_" + pMesh->getLabel() + "_simulationJob.json"));
+=======
+            .append(pMesh->getLabel() + "_simScenario.json").string());
+>>>>>>> 2599d472614f403f8249060272ce288033a919d4
+    std::ofstream jsonFile(jsonFilename);
+    jsonFile << json;
+    std::cout << "Saved simulation job as:" << jsonFilename << std::endl;
+
+    return returnValue;
+  }
+};
+
+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
+
+struct SimulationResults {
+  bool converged{true};
+  std::shared_ptr<SimulationJob> job;
+  SimulationHistory history;
+  std::vector<Vector6d> displacements;
+  double executionTime{0};
+  std::string labelPrefix{"deformed"};
+  inline static char deliminator{' '};
+
+  void setLabelPrefix(const std::string &lp) {
+    labelPrefix += deliminator + lp;
+  }
+  std::string getLabel() const {
+    return labelPrefix + deliminator + job->pMesh->getLabel() + deliminator +
+           job->getLabel();
+  }
+  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;
+    }
+    polyscope::removeCurveNetwork(getLabel());
+  }
+  void registerForDrawing() const {
+    polyscope::options::groundPlaneEnabled = false;
+    polyscope::view::upDir = polyscope::view::UpDir::ZUp;
+    const std::string branchName = "Branch:Polyscope";
+    polyscope::options::programName = branchName;
+    if (!polyscope::state::initialized) {
+      polyscope::init();
+    } /* else {
+       polyscope::removeAllStructures();
+     }*/
+    //    const std::string undeformedMeshName = "Undeformed_" + label;
+    //    polyscope::registerCurveNetwork(
+    //        undeformedMeshName, mesh->getEigenVertices(),
+    //        mesh->getEigenEdges());
+
+    const std::shared_ptr<SimulationMesh> &mesh = job->pMesh;
+    polyscope::registerCurveNetwork(getLabel(), mesh->getEigenVertices(),
+                                    mesh->getEigenEdges())
+        ->setRadius(0.0007, false);
+    Eigen::MatrixX3d nodalDisplacements(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]);
+    }
+    polyscope::getCurveNetwork(getLabel())
+        ->updateNodePositions(mesh->getEigenVertices() + nodalDisplacements);
+
+    job->registerForDrawing(getLabel());
+  }
+
+  void saveDeformedModel() {
+    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]);
+    }
+    m.savePly(getLabel() + ".ply");
+  }
+
+  void save() {
+    const std::string filename(getLabel() + "_displacements.eigenBin");
+
+    Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements);
+    Eigen::write_binary(filename, m);
+  }
+
+  // The comparison of the results happens comparing the 6-dof nodal
+  // displacements
+  bool isEqual(const Eigen::MatrixXd &nodalDisplacements) {
+    assert(nodalDisplacements.cols() == 6);
+    Eigen::MatrixXd eigenDisplacements =
+        Utilities::toEigenMatrix(this->displacements);
+    const double errorNorm = (eigenDisplacements - nodalDisplacements).norm();
+    return errorNorm < 1e-10;
+    //    return eigenDisplacements.isApprox(nodalDisplacements);
+  }
+};
+
+#endif // SIMULATIONHISTORY_HPP
diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp
index 04a85ff..01a2d29 100644
--- a/topologyenumerator.cpp
+++ b/topologyenumerator.cpp
@@ -130,7 +130,7 @@ void TopologyEnumerator::computeValidPatterns(
   if (debugIsOn) {
     // Export all valid edges in a ply
     patternAllValidEdges.savePly(
-        std::filesystem::path(resultsPath).append("allValidEdges.ply"));
+        std::filesystem::path(resultsPath).append("allValidEdges.ply").string());
   }
   //  statistics.numberOfValidEdges = validEdges.size();
 
@@ -357,7 +357,7 @@ std::vector<vcg::Point2i> TopologyEnumerator::getValidEdges(
           patternDuplicateEdges, p0, p1);
     }
     patternDuplicateEdges.savePly(
-        std::filesystem::path(duplicateEdgesPath).append("duplicateEdges.ply"));
+        std::filesystem::path(duplicateEdgesPath).append("duplicateEdges.ply").string());
   }
   statistics.numberOfDuplicateEdges = duplicateEdges.size();
 
diff --git a/topologyenumerator.hpp b/topologyenumerator.hpp
index 47d8cff..1b3202b 100644
--- a/topologyenumerator.hpp
+++ b/topologyenumerator.hpp
@@ -156,7 +156,7 @@ private:
                       const size_t &numberOfDesiredEdges) const;
   std::vector<vcg::Point2i>
   getValidEdges(const std::vector<size_t> &numberOfNodesPerSlot,
-                const std::filesystem::__cxx11::path &resultsPath,
+                const std::filesystem::path &resultsPath,
                 const FlatPatternGeometry &patternGeometryAllEdges,
                 const std::vector<vcg::Point2i> &allPossibleEdges);
   std::unordered_set<size_t> computeDuplicateEdges();
diff --git a/vcgtrimesh.cpp b/vcgtrimesh.cpp
index e68b312..6acbc35 100644
--- a/vcgtrimesh.cpp
+++ b/vcgtrimesh.cpp
@@ -11,7 +11,7 @@ void VCGTriMesh::loadFromPlyFile(const std::string &filename) {
   mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_EDGEINDEX;
   mask |= nanoply::NanoPlyWrapper<VCGTriMesh>::IO_FACEINDEX;
   if (nanoply::NanoPlyWrapper<VCGTriMesh>::LoadModel(
-          std::filesystem::absolute(filename).c_str(), *this, mask) != 0) {
+          std::filesystem::absolute(filename).string().c_str(), *this, mask) != 0) {
     std::cout << "Could not load tri mesh" << std::endl;
   }
   vcg::tri::UpdateTopology<VCGTriMesh>::FaceFace(*this);
@@ -57,7 +57,7 @@ bool VCGTriMesh::savePly(const std::string plyFilename) {
 VCGTriMesh::VCGTriMesh() {}
 
 VCGTriMesh::VCGTriMesh(const std::string &filename) {
-  const std::string extension = std::filesystem::path(filename).extension();
+  const std::string extension = std::filesystem::path(filename).extension().string();
   if (extension == ".ply") {
     loadFromPlyFile(filename);
   } else if (extension == ".obj") {