Refactoring

beamformfinder.cpp

@@ -1,7 +1,4 @@
 #include "beamformfinder.hpp"
-#include "polyscope/curve_network.h"
-#include "polyscope/histogram.h"
-#include "polyscope/polyscope.h"
 #include "simulationhistoryplotter.hpp"
 #include <algorithm>
 #include <chrono>
@@ -18,8 +15,9 @@ void FormFinder::runUnitTests() {
   VCGEdgeMesh beam;
   //  const size_t spanGridSize = 11;
   //  mesh.createSpanGrid(spanGridSize);
-  beam.loadPly(
+  beam.loadPly(std::filesystem::path(groundOfTruthFolder)
-      std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply").string());
+                   .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 +48,8 @@ void FormFinder::runUnitTests() {
   simpleBeam_simulationResults.save();
   const std::string simpleBeamGroundOfTruthBinaryFilename =
       std::filesystem::path(groundOfTruthFolder)
-          .append("SimpleBeam_displacements.eigenBin").string();
+          .append("SimpleBeam_displacements.eigenBin")
+          .string();
   assert(std::filesystem::exists(
       std::filesystem::path(simpleBeamGroundOfTruthBinaryFilename)));
   Eigen::MatrixXd simpleBeam_groundOfTruthDisplacements;
@@ -1568,183 +1567,6 @@ FormFinder::computeResults(const std::shared_ptr<SimulationJob> &pJob) {
   return SimulationResults{true, pJob, history, displacements};
 }
 
-void FormFinder::draw(const std::string &screenshotsFolder = {}) {
-  // update positions
-  //  polyscope::getCurveNetwork("Undeformed edge mesh")->setEnabled(false);
-  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->updateNodePositions(pMesh->getEigenVertices());
-
-  // Vertex quantities
-  std::vector<double> kineticEnergies(pMesh->VN());
-  std::vector<std::array<double, 3>> nodeNormals(pMesh->VN());
-  std::vector<std::array<double, 3>> internalForces(pMesh->VN());
-  std::vector<std::array<double, 3>> externalForces(pMesh->VN());
-  std::vector<std::array<double, 3>> externalMoments(pMesh->VN());
-  std::vector<double> internalForcesNorm(pMesh->VN());
-  std::vector<std::array<double, 3>> internalAxialForces(pMesh->VN());
-  std::vector<std::array<double, 3>> internalFirstBendingTranslationForces(
-      pMesh->VN());
-  std::vector<std::array<double, 3>> internalFirstBendingRotationForces(
-      pMesh->VN());
-  std::vector<std::array<double, 3>> internalSecondBendingTranslationForces(
-      pMesh->VN());
-  std::vector<std::array<double, 3>> internalSecondBendingRotationForces(
-      pMesh->VN());
-  std::vector<double> nRs(pMesh->VN());
-  std::vector<double> theta2(pMesh->VN());
-  std::vector<double> theta3(pMesh->VN());
-  std::vector<std::array<double, 3>> residualForces(pMesh->VN());
-  std::vector<double> residualForcesNorm(pMesh->VN());
-  std::vector<double> accelerationX(pMesh->VN());
-  for (const VertexType &v : pMesh->vert) {
-    kineticEnergies[pMesh->getIndex(v)] = pMesh->nodes[v].kineticEnergy;
-    const VectorType n = v.cN();
-    nodeNormals[pMesh->getIndex(v)] = {n[0], n[1], n[2]};
-    // per node internal forces
-    const Vector6d nodeforce = pMesh->nodes[v].force.internal * (-1);
-    internalForces[pMesh->getIndex(v)] = {nodeforce[0], nodeforce[1],
-                                          nodeforce[2]};
-    internalForcesNorm[pMesh->getIndex(v)] = nodeforce.norm();
-    // External force
-    const Vector6d nodeExternalForce = pMesh->nodes[v].force.external;
-    externalForces[pMesh->getIndex(v)] = {
-        nodeExternalForce[0], nodeExternalForce[1], nodeExternalForce[2]};
-    externalMoments[pMesh->getIndex(v)] = {nodeExternalForce[3],
-                                           nodeExternalForce[4], 0};
-    internalAxialForces[pMesh->getIndex(v)] = {nodeforce[0], nodeforce[1],
-                                               nodeforce[2]};
-    const Node &node = pMesh->nodes[v];
-    const Vector6d nodeForceFirst = node.force.internalFirstBending * (-1);
-    internalFirstBendingTranslationForces[pMesh->getIndex(v)] = {
-        nodeForceFirst[0], nodeForceFirst[1], nodeForceFirst[2]};
-    internalFirstBendingRotationForces[pMesh->getIndex(v)] = {
-        nodeForceFirst[3], nodeForceFirst[4], 0};
-
-    const Vector6d nodeForceSecond = node.force.internalSecondBending * (-1);
-    internalSecondBendingTranslationForces[pMesh->getIndex(v)] = {
-        nodeForceSecond[0], nodeForceSecond[1], nodeForceSecond[2]};
-    internalSecondBendingRotationForces[pMesh->getIndex(v)] = {
-        nodeForceSecond[3], nodeForceSecond[4], 0};
-    nRs[pMesh->getIndex(v)] = node.nR;
-    const Vector6d nodeResidualForce = pMesh->nodes[v].force.residual;
-    residualForces[pMesh->getIndex(v)] = {
-        nodeResidualForce[0], nodeResidualForce[1], nodeResidualForce[2]};
-    residualForcesNorm[pMesh->getIndex(v)] = nodeResidualForce.norm();
-    accelerationX[pMesh->getIndex(v)] = pMesh->nodes[v].acceleration[0];
-  }
-  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->addNodeScalarQuantity("Kinetic Energy", kineticEnergies)
-      ->setEnabled(false);
-  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->addNodeVectorQuantity("Node normals", nodeNormals)
-      ->setEnabled(true);
-  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->addNodeVectorQuantity("Internal force", internalForces)
-      ->setEnabled(false);
-  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->addNodeVectorQuantity("External force", externalForces)
-      ->setEnabled(true);
-  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->addNodeVectorQuantity("External Moments", externalMoments)
-      ->setEnabled(true);
-  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->addNodeScalarQuantity("Internal force norm", internalForcesNorm)
-      ->setEnabled(true);
-  //  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  //      ->addNodeVectorQuantity("Internal Axial force", internalAxialForces)
-  //      ->setEnabled(false);
-  //  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  //      ->addNodeVectorQuantity("First bending force-Translation",
-  //                              internalFirstBendingTranslationForces)
-  //      ->setEnabled(false);
-  //  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  //      ->addNodeVectorQuantity("First bending force-Rotation",
-  //                              internalFirstBendingRotationForces)
-  //      ->setEnabled(false);
-
-  //  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  //      ->addNodeVectorQuantity("Second bending force-Translation",
-  //                              internalSecondBendingTranslationForces)
-  //      ->setEnabled(false);
-  //  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  //      ->addNodeVectorQuantity("Second bending force-Rotation",
-  //                              internalSecondBendingRotationForces)
-  //      ->setEnabled(false);
-  //  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  //      ->addNodeScalarQuantity("nR", nRs)
-  //      ->setEnabled(false);
-  //  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  //      ->addNodeScalarQuantity("theta3", theta3)
-  //      ->setEnabled(false);
-  //  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  //      ->addNodeScalarQuantity("theta2", theta2)
-  //      ->setEnabled(false);
-  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->addNodeVectorQuantity("Residual force", residualForces)
-      ->setEnabled(false);
-  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->addNodeScalarQuantity("Residual force norm", residualForcesNorm)
-      ->setEnabled(false);
-  //  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  //      ->addNodeScalarQuantity("Node acceleration x", accelerationX);
-
-  // Edge quantities
-  std::vector<double> A(pMesh->EN());
-  std::vector<double> J(pMesh->EN());
-  std::vector<double> I2(pMesh->EN());
-  std::vector<double> I3(pMesh->EN());
-  for (const EdgeType &e : pMesh->edge) {
-    const size_t ei = pMesh->getIndex(e);
-    A[ei] = pMesh->elements[e].A;
-    J[ei] = pMesh->elements[e].J;
-    I2[ei] = pMesh->elements[e].I2;
-    I3[ei] = pMesh->elements[e].I3;
-  }
-
-  polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A);
-  polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("J", J);
-  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->addEdgeScalarQuantity("I2", I2);
-  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->addEdgeScalarQuantity("I3", I3);
-
-  // Specify the callback
-  polyscope::state::userCallback = [&]() {
-    // Since options::openImGuiWindowForUserCallback == true by default,
-    // we can immediately start using ImGui commands to build a UI
-
-    ImGui::PushItemWidth(100); // Make ui elements 100 pixels wide,
-                               // instead of full width. Must have
-                               // matching PopItemWidth() below.
-
-    ImGui::InputInt("Simulation drawing step",
-                    &mSettings.drawingStep); // set a int variable
-    ImGui::Checkbox("Enable drawing",
-                    &mSettings.shouldDraw); // set a int variable
-    ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep);
-
-    ImGui::PopItemWidth();
-  };
-
-  if (!screenshotsFolder.empty()) {
-    static bool firstDraw = true;
-    if (firstDraw) {
-      for (const auto &entry :
-           std::filesystem::directory_iterator(screenshotsFolder))
-        std::filesystem::remove_all(entry.path());
-      //      polyscope::view::processZoom(5);
-      firstDraw = false;
-    }
-    polyscope::screenshot(
-        std::filesystem::path(screenshotsFolder)
-            .append(std::to_string(mCurrentSimulationStep) + ".png")
-            .string(),
-        false);
-  } else {
-    polyscope::show();
-  }
-}
-
 void FormFinder::printCurrentState() const {
   std::cout << "Simulation steps executed:" << mCurrentSimulationStep
             << std::endl;
@@ -1798,13 +1620,6 @@ FormFinder::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
               << " nodes and " << pMesh->EN() << " elements." << std::endl;
   }
   computeRigidSupports();
-
-  if (mSettings.shouldDraw ) {
-    initPolyscope();
-    polyscope::registerCurveNetwork(
-        meshPolyscopeLabel, pMesh->getEigenVertices(), pMesh->getEigenEdges());
-    //    registerWorldAxes();
-  }
   for (auto fixedVertex : pJob->constrainedVertices) {
     assert(fixedVertex.first < pMesh->VN());
   }
@@ -1853,25 +1668,9 @@ FormFinder::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
     }
 
     if (std::isnan(pMesh->currentTotalKineticEnergy)) {
-      std::time_t now = time(0);
+      if (mSettings.beVerbose) {
-      std::tm *ltm = std::localtime(&now);
-      std::string dir = "../ProblematicSimulationJobs/" +
-                        std::to_string(ltm->tm_mday) + "_" +
-                        std::to_string(1 + ltm->tm_mon) + "_" +
-                        std::to_string(1900 + ltm->tm_year);
-      const std::string subDir = std::filesystem::path(dir)
-                                     .append(std::to_string(ltm->tm_hour) +
-                                             ":" + std::to_string(ltm->tm_min))
-                                     .string();
-      std::filesystem::create_directories(subDir);
-      pJob->save(subDir);
         std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl;
-      FormFinder debug;
+      }
-      FormFinder::Settings settings;
-      settings.shouldDraw = true;
-      settings.beVerbose = true;
-      debug.executeSimulation(pJob, settings);
-      std::terminate();
       break;
     }
 
@@ -1884,8 +1683,6 @@ currentSimulationStep > maxDRMIterations*/
       //                               .append("Screenshots")
       //                               .string();
       //      draw(saveTo);
-      draw();
-
       //      yValues = history.kineticEnergy;
       //      plotHandle = matplot::scatter(xPlot, yValues);
       //      label = "Log of Kinetic energy";
@@ -1990,12 +1787,6 @@ mesh->currentTotalPotentialEnergykN*/
       SimulationResultsReporter reporter;
       reporter.reportResults({results}, "Results", pJob->pMesh->getLabel());
     }
-    if (mSettings.shouldDraw) {
-      draw();
-    }
-  }
-  if (polyscope::hasCurveNetwork(meshPolyscopeLabel)) {
-    polyscope::removeCurveNetwork(meshPolyscopeLabel);
   }
   return results;
 }

beamformfinder.hpp

@@ -3,8 +3,8 @@
 
 #include "elementalmesh.hpp"
 #include "matplot/matplot.h"
-#include "polyscope/curve_network.h"
+//#include "polyscope/curve_network.h"
-#include "polyscope/polyscope.h"
+//#include "polyscope/polyscope.h"
 #include "simulationresult.hpp"
 #include <Eigen/Dense>
 #include <filesystem>

csvfile.hpp

@@ -5,12 +5,12 @@
 #include <sstream>
 #include <string>
 
-class csvfile;
+class csvFile;
 
-inline static csvfile &endrow(csvfile &file);
+inline static csvFile &endrow(csvFile &file);
-inline static csvfile &flush(csvfile &file);
+inline static csvFile &flush(csvFile &file);
 
-class csvfile {
+class csvFile {
   std::ofstream fs_;
   bool is_first_;
   const std::string separator_;
@@ -18,14 +18,19 @@ class csvfile {
   const std::string special_chars_;
 
 public:
-  csvfile(const std::string filename, const bool &overwrite,
+  csvFile(const std::string filename, const bool &overwrite,
-          const std::string separator = ";")
+          const std::string separator = ",")
       : fs_(), is_first_(true), separator_(separator), escape_seq_("\""),
         special_chars_("\"") {
     fs_.exceptions(std::ios::failbit | std::ios::badbit);
-    if (!std::filesystem::exists(std::filesystem::path("../OptimizationResults")
+    if (filename.empty()) {
-        .append("statistics.csv")
+      fs_.copyfmt(std::cout);
-    )) {
+      fs_.clear(std::cout.rdstate());
+      fs_.basic_ios<char>::rdbuf(std::cout.rdbuf());
+    } else {
+      if (!std::filesystem::exists(
+              std::filesystem::path("../OptimizationResults")
+                  .append("statistics.csv"))) {
         std::ofstream outfile(std::filesystem::path("../OptimizationResults")
                                   .append("statistics.csv")
                                   .string());
@@ -34,8 +39,9 @@ public:
       overwrite ? fs_.open(filename, std::ios::trunc)
                 : fs_.open(filename, std::ios::app);
     }
+  }
 
-  ~csvfile() {
+  ~csvFile() {
     flush();
     fs_.close();
   }
@@ -47,16 +53,16 @@ public:
     is_first_ = true;
   }
 
-  csvfile &operator<<(csvfile &(*val)(csvfile &)) { return val(*this); }
+  csvFile &operator<<(csvFile &(*val)(csvFile &)) { return val(*this); }
 
-  csvfile &operator<<(const char *val) { return write(escape(val)); }
+  csvFile &operator<<(const char *val) { return write(escape(val)); }
 
-  csvfile &operator<<(const std::string &val) { return write(escape(val)); }
+  csvFile &operator<<(const std::string &val) { return write(escape(val)); }
 
-  template <typename T> csvfile &operator<<(const T &val) { return write(val); }
+  template <typename T> csvFile &operator<<(const T &val) { return write(val); }
 
 private:
-  template <typename T> csvfile &write(const T &val) {
+  template <typename T> csvFile &write(const T &val) {
     if (!is_first_) {
       fs_ << separator_;
     } else {
@@ -80,12 +86,12 @@ private:
   }
 };
 
-inline static csvfile &endrow(csvfile &file) {
+inline static csvFile &endrow(csvFile &file) {
   file.endrow();
   return file;
 }
 
-inline static csvfile &flush(csvfile &file) {
+inline static csvFile &flush(csvFile &file) {
   file.flush();
   return file;
 }

edgemesh.cpp

@@ -193,8 +193,11 @@ bool VCGEdgeMesh::loadPly(const std::string plyFilename) {
   vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
   label = std::filesystem::path(plyFilename).stem().string();
 
+  const bool printDebugInfo = false;
+  if (printDebugInfo) {
     std::cout << plyFilename << " was loaded successfuly." << std::endl;
     std::cout << "Mesh has " << EN() << " edges." << std::endl;
+  }
   return true;
 }
 
@@ -337,9 +340,27 @@ void VCGEdgeMesh::printVertexCoordinates(const size_t &vi) const {
             << " " << vert[vi].cP()[2] << std::endl;
 }
 
-void VCGEdgeMesh::registerForDrawing() const {
+#ifdef POLYSCOPE_DEFINED
+void VCGEdgeMesh::registerForDrawing(
+    const std::optional<glm::vec3> &desiredColor) const {
   initPolyscope();
-  const double drawingRadius = 0.0007;
+  const double drawingRadius = 0.002;
-  polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges())
+  auto polyscopeHandle_edgeMesh =
-      ->setRadius(drawingRadius, false);
+      polyscope::registerCurveNetwork(label, getEigenVertices(),
+                                      getEigenEdges())
+          ->setRadius(drawingRadius, true);
+  if (desiredColor.has_value()) {
+    polyscopeHandle_edgeMesh->setColor(desiredColor.value());
   }
+}
+
+void VCGEdgeMesh::unregister() const {
+  if (!polyscope::hasCurveNetwork(label)) {
+    std::cerr << "No curve network registered with a name: " << getLabel()
+              << std::endl;
+    std::cerr << "Nothing to remove." << std::endl;
+    return;
+  }
+  polyscope::removeCurveNetwork(label);
+}
+#endif

edgemesh.hpp

@@ -1,6 +1,7 @@
 #ifndef EDGEMESH_HPP
 #define EDGEMESH_HPP
 #include "beam.hpp"
+#include "externvariables.hpp"
 #include "polymesh.hpp"
 #include "utilities.hpp"
 #include "vcgtrimesh.hpp"
@@ -77,8 +78,11 @@ public:
   Eigen::MatrixX3d getEigenVertices() const;
   Eigen::MatrixX3d getEigenEdgeNormals() const;
   void printVertexCoordinates(const size_t &vi) const;
-  void registerForDrawing() const;
+#ifdef POLYSCOPE_DEFINED
-
+  void registerForDrawing(
+      const std::optional<glm::vec3> &desiredColor = std::nullopt) const;
+  void unregister() const;
+#endif
   std::string getLabel() const;
   void setLabel(const std::string &value);
 

elementalmesh.cpp

@@ -38,10 +38,11 @@ SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) {
   eigenVertices = mesh.getEigenVertices();
 }
 
-SimulationMesh::~SimulationMesh()
+SimulationMesh::~SimulationMesh() {
-{
+  vcg::tri::Allocator<VCGEdgeMesh>::DeletePerEdgeAttribute<Element>(*this,
-    vcg::tri::Allocator<VCGEdgeMesh>::DeletePerEdgeAttribute<Element>(*this, elements);
+                                                                    elements);
-    vcg::tri::Allocator<VCGEdgeMesh>::DeletePerVertexAttribute<Node>(*this,nodes);
+  vcg::tri::Allocator<VCGEdgeMesh>::DeletePerVertexAttribute<Node>(*this,
+                                                                   nodes);
 }
 
 SimulationMesh::SimulationMesh(FlatPattern &pattern) {
@@ -305,6 +306,7 @@ bool SimulationMesh::loadPly(const string &plyFilename) {
   mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
   mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
   mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
+  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_MESHATTRIB;
   if (nanoply::NanoPlyWrapper<SimulationMesh>::LoadModel(
           plyFilename.c_str(), *this, mask, customAttrib) != 0) {
     return false;
@@ -358,19 +360,6 @@ bool SimulationMesh::savePly(const std::string &plyFilename) {
   customAttrib.AddEdgeAttribDescriptor<vcg::Point2d, double, 2>(
       plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT64,
       material.data());
-  //  std::vector<std::array<double, 6>> beamProperties(EN());
-  //  for (size_t ei = 0; ei < EN(); ei++) {
-  //    auto props = elements[ei].toArray();
-  //    for (auto p : props) {
-  //      std::cout << p << " ";
-  //    }
-  //    std::cout << std::endl;
-  //    beamProperties[ei] = props;
-  //  }
-  //  customAttrib.AddEdgeAttribDescriptor<std::array<double, 6>, double, 6>(
-  //      plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64,
-  //      beamProperties.data());
-  // Load the ply file
   unsigned int mask = 0;
   mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
   mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
@@ -386,11 +375,6 @@ bool SimulationMesh::savePly(const std::string &plyFilename) {
 SimulationMesh::EdgePointer
 SimulationMesh::getReferenceElement(const VCGEdgeMesh::VertexType &v) {
   const VertexIndex vi = getIndex(v);
-  //  return nodes[v].incidentElements[0];
-  //  if (vi == 0 || vi == 1) {
-  //    return nodes[0].incidentElements[0];
-  //  }
-
   return nodes[v].incidentElements[0];
 }
 

optimizationstructs.hpp

@@ -0,0 +1,5 @@
+#ifndef OPTIMIZATIONSTRUCTS_HPP
+#define OPTIMIZATIONSTRUCTS_HPP
+#include "beamformfinder.hpp"
+
+#endif // OPTIMIZATIONSTRUCTS_HPP

simulationresult.hpp

@@ -60,8 +60,8 @@ class SimulationJob {
     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
+        "meshLabel"}; // TODO: should be in the savePly function of the
-                  // mesh class
+                      // simulation mesh class
   } jsonLabels;
 
 public:
@@ -79,24 +79,155 @@ public:
         nodalExternalForces(ef), nodalForcedDisplacements(fd) {}
   SimulationJob() {}
   SimulationJob(const std::string &jsonFilename) { load(jsonFilename); }
+
+  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[jsonLabel_constrainedVertices] = constrainedVertices;
+      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[jsonLabel_nodalForces] = arrForces;
+      json[jsonLabels.nodalForces] = arrForces;
     }
 
     return json.dump();
   }
+  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);
+    ifs >> json;
+
+    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->loadPly(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.constrainedVertices)) {
+      constrainedVertices =
+          //      auto conV =
+          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;
+    }
+
+    if (json.contains(jsonLabels.nodalForces)) {
+      auto f = std::unordered_map<VertexIndex, std::array<double, 6>>(
+          json[jsonLabels.nodalForces]);
+      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;
+    std::string jsonFilename(
+        std::filesystem::path(pathFolderDirectory)
+            .append(label + "_" + pMesh->getLabel() + "_simulationJob.json")
+            .string());
+
+    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] = std::filesystem::relative(
+        std::filesystem::path(meshFilename),
+        std::filesystem::path(
+            std::filesystem::path(jsonFilename).parent_path()));
+    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 (!label.empty()) {
+      json[jsonLabels.label] = label;
+    }
+    if (!pMesh->getLabel().empty()) {
+      json[jsonLabels.meshLabel] = pMesh->getLabel();
+    }
+    std::ofstream jsonFile(jsonFilename);
+    jsonFile << json;
+    //    std::cout << "Saved simulation job as:" << jsonFilename << std::endl;
+
+    return returnValue;
+  }
+#ifdef POLYSCOPE_DEFINED
   void registerForDrawing(const std::string &meshLabel) const {
     initPolyscope();
     if (meshLabel.empty()) {
@@ -163,108 +294,7 @@ public:
           ->setEnabled(false);
     }
   }
-
+#endif // POLYSCOPE_DEFINED
-  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 =
-              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;
-    }
-
-    if (json.contains(jsonLabels.nodalForces)) {
-      auto f = std::unordered_map<VertexIndex, std::array<double, 6>>(
-          json[jsonLabels.nodalForces]);
-      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)
-            .append(label + "_" + pMesh->getLabel() + "_simulationJob.json").string());
-    std::ofstream jsonFile(jsonFilename);
-    jsonFile << json;
-    std::cout << "Saved simulation job as:" << jsonFilename << std::endl;
-
-    return returnValue;
-  }
 };
 
 namespace Eigen {
@@ -300,6 +330,16 @@ struct SimulationResults {
   std::string labelPrefix{"deformed"};
   inline static char deliminator{' '};
 
+  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;
   }
@@ -307,45 +347,6 @@ struct SimulationResults {
     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;
@@ -376,6 +377,53 @@ struct SimulationResults {
     return errorNorm < 1e-10;
     //    return eigenDisplacements.isApprox(nodalDisplacements);
   }
+
+#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;
+    }
+    polyscope::removeCurveNetwork(getLabel());
+  }
+  void registerForDrawing(
+      const std::optional<glm::vec3> &desiredColor = std::nullopt) 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;
+    auto polyscopeHandle_deformedEdmeMesh =
+        polyscope::registerCurveNetwork(getLabel(), mesh->getEigenVertices(),
+                                        mesh->getEigenEdges())
+            ->setRadius(0.002, true);
+    if (desiredColor.has_value()) {
+      polyscopeHandle_deformedEdmeMesh->setColor(desiredColor.value());
+    }
+    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]);
+    }
+    polyscopeHandle_deformedEdmeMesh->updateNodePositions(
+        mesh->getEigenVertices() + nodalDisplacements);
+
+    job->registerForDrawing(getLabel());
+  }
+#endif
 };
 
 #endif // SIMULATIONHISTORY_HPP

utilities.hpp

@@ -98,6 +98,11 @@ struct Vector6d : public std::array<double, 6> {
       return true;
     });
   }
+
+  Eigen::Vector3d getTranslation() const {
+    return Eigen::Vector3d(this->operator[](0), this->operator[](1),
+                           this->operator[](2));
+  }
 };
 
 namespace Utilities {
@@ -155,10 +160,7 @@ inline Eigen::MatrixXd toEigenMatrix(const std::vector<Vector6d> &v) {
 
 } // namespace Utilities
 
-// namespace ConfigurationFile {
+#ifdef POLYSCOPE_DEFINED
-
-//}
-//} // namespace ConfigurationFile
 #include "polyscope/curve_network.h"
 #include "polyscope/polyscope.h"
 
@@ -204,7 +206,12 @@ inline void registerWorldAxes() {
       ->addEdgeColorQuantity(worldAxesColorName, axesColors)
       ->setEnabled(true);
 }
+#endif
 
+// namespace ConfigurationFile {
+
+//}
+//} // namespace ConfigurationFile
 template <typename T1, typename T2>
 void constructInverseMap(const T1 &map, T2 &oppositeMap) {
   assert(!map.empty());