Added loading and saving of simulation jobs. Saving of the const factors of each element when exporting a ply file. Refactoring

beam.hpp

@@ -1,6 +1,5 @@
 #ifndef BEAM_HPP
 #define BEAM_HPP
-#include <Eigen/Dense>
 #include <assert.h>
 #include <cmath>
 #include <iostream>

beamformfinder.cpp

@@ -31,17 +31,16 @@ void FormFinder::runUnitTests() {
   nodalForcedDisplacements.insert({beam.VN() - 1, Eigen::Vector3d(-0.2, 0, 0)});
 
   SimulationJob beamSimulationJob{
-      std::make_shared<SimulationMesh>(beam),
+      std::make_shared<SimulationMesh>(beam), "First paper example",
       //      SimulationJob::constructFixedVerticesSpanGrid(spanGridSize,
       //      mesh.VN()),
       fixedVertices, nodalForces, nodalForcedDisplacements};
-  beamSimulationJob.mesh->setBeamMaterial(0.3, 200);
+  beamSimulationJob.pMesh->setBeamMaterial(0.3, 200);
   assert(CrossSectionType::name == CylindricalBeamDimensions::name);
 
-  beamSimulationJob.mesh->setBeamCrossSection(CrossSectionType{0.03, 0.026});
+  beamSimulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026});
-  SimulationResults simpleBeam_simulationResults =
+  SimulationResults simpleBeam_simulationResults = formFinder.executeSimulation(
-      formFinder.executeSimulation(beamSimulationJob, false, true);
+      std::make_shared<SimulationJob>(beamSimulationJob), false, true);
-  simpleBeam_simulationResults.label = "SimpleBeam";
   simpleBeam_simulationResults.save();
   const std::string simpleBeamGroundOfTruthBinaryFilename =
       std::filesystem::path(groundOfTruthFolder)
@@ -82,17 +81,18 @@ void FormFinder::runUnitTests() {
 
   SimulationJob shortSpanGridshellSimulationJob{
       std::make_shared<SimulationMesh>(shortSpanGrid),
+      "Short span gridshell",
       fixedVertices,
       {},
       nodalForcedDisplacements};
-  shortSpanGridshellSimulationJob.mesh->setBeamMaterial(0.3, 200);
+  shortSpanGridshellSimulationJob.pMesh->setBeamMaterial(0.3, 200);
   assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
-  shortSpanGridshellSimulationJob.mesh->setBeamCrossSection(
+  shortSpanGridshellSimulationJob.pMesh->setBeamCrossSection(
       CrossSectionType{0.03, 0.026});
   SimulationResults shortSpanGridshellSimulationResults =
-      formFinder.executeSimulation(shortSpanGridshellSimulationJob, false,
+      formFinder.executeSimulation(
-                                   false);
+          std::make_shared<SimulationJob>(shortSpanGridshellSimulationJob),
-  shortSpanGridshellSimulationResults.label = "ShortSpanGridshell";
+          false, false);
   shortSpanGridshellSimulationResults.save();
 
   const std::string groundOfTruthBinaryFilename =
@@ -159,21 +159,22 @@ void FormFinder::runUnitTests() {
 
   SimulationJob longSpanGridshell_simulationJob{
       std::make_shared<SimulationMesh>(longSpanGrid),
+      "long span gridshell",
       fixedVertices,
       {},
       nodalForcedDisplacements};
-  longSpanGridshell_simulationJob.mesh->setBeamMaterial(0.3, 200);
+  longSpanGridshell_simulationJob.pMesh->setBeamMaterial(0.3, 200);
   if (typeid(CrossSectionType) != typeid(CylindricalBeamDimensions)) {
     std::cerr << "A cylindrical cross section is expected for running the "
                  "paper examples."
               << std::endl;
   }
-  longSpanGridshell_simulationJob.mesh->setBeamCrossSection(
+  longSpanGridshell_simulationJob.pMesh->setBeamCrossSection(
       CrossSectionType{0.03, 0.026});
   SimulationResults longSpanGridshell_simulationResults =
-      formFinder.executeSimulation(longSpanGridshell_simulationJob, false,
+      formFinder.executeSimulation(
-                                   false);
+          std::make_shared<SimulationJob>(longSpanGridshell_simulationJob),
-  longSpanGridshell_simulationResults.label = "LongSpanGridshell";
+          false, false);
   longSpanGridshell_simulationResults.save();
 
   const std::string longSpan_groundOfTruthBinaryFilename =
@@ -200,11 +201,11 @@ void FormFinder::runUnitTests() {
 }
 
 void FormFinder::setTotalResidualForcesNormThreshold(double value) {
-  totalResidualForcesNormThreshold = value;
+  settings.totalResidualForcesNormThreshold = value;
 }
 
 void FormFinder::reset() {
-  Dt = Dtini;
+  settings.Dt = settings.Dtini;
   mCurrentSimulationStep = 0;
   history.clear();
   constrainedVertices.clear();
@@ -215,6 +216,7 @@ void FormFinder::reset() {
   checkedForMaximumMoment = false;
   shouldUseKineticEnergyThreshold = false;
   externalMomentsNorm = 0;
+  settings.drawingStep = 1;
 }
 
 VectorType FormFinder::computeDisplacementDifferenceDerivative(
@@ -1160,24 +1162,24 @@ void FormFinder::updateNodalMasses() {
       assert(rotationalSumSk_J != 0);
     }
     mesh->nodes[v].translationalMass =
-        gamma * pow(Dtini, 2) * 2 * translationalSumSk;
+        gamma * pow(settings.Dtini, 2) * 2 * translationalSumSk;
     mesh->nodes[v].rotationalMass_I2 =
-        gamma * pow(Dtini, 2) * 8 * rotationalSumSk_I2;
+        gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I2;
     mesh->nodes[v].rotationalMass_I3 =
-        gamma * pow(Dtini, 2) * 8 * rotationalSumSk_I3;
+        gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I3;
     mesh->nodes[v].rotationalMass_J =
-        gamma * pow(Dtini, 2) * 8 * rotationalSumSk_J;
+        gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_J;
 
-    assert(std::pow(Dtini, 2.0) * translationalSumSk /
+    assert(std::pow(settings.Dtini, 2.0) * translationalSumSk /
                mesh->nodes[v].translationalMass <
            2);
-    assert(std::pow(Dtini, 2.0) * rotationalSumSk_I2 /
+    assert(std::pow(settings.Dtini, 2.0) * rotationalSumSk_I2 /
                mesh->nodes[v].rotationalMass_I2 <
            2);
-    assert(std::pow(Dtini, 2.0) * rotationalSumSk_I3 /
+    assert(std::pow(settings.Dtini, 2.0) * rotationalSumSk_I3 /
                mesh->nodes[v].rotationalMass_I3 <
            2);
-    assert(std::pow(Dtini, 2.0) * rotationalSumSk_J /
+    assert(std::pow(settings.Dtini, 2.0) * rotationalSumSk_J /
                mesh->nodes[v].rotationalMass_J <
            2);
   }
@@ -1209,7 +1211,7 @@ void FormFinder::updateNodalVelocities() {
   for (VertexType &v : mesh->vert) {
     const VertexIndex vi = mesh->getIndex(v);
     Node &node = mesh->nodes[v];
-    node.velocity = node.velocity + node.acceleration * Dt;
+    node.velocity = node.velocity + node.acceleration * settings.Dt;
   }
   updateKineticEnergy();
 }
@@ -1217,7 +1219,14 @@ void FormFinder::updateNodalVelocities() {
 void FormFinder::updateNodalDisplacements() {
   for (VertexType &v : mesh->vert) {
     Node &node = mesh->nodes[v];
-    node.displacements = node.displacements + node.velocity * Dt;
+    node.displacements = node.displacements + node.velocity * settings.Dt;
+    if (settings.beVerbose) {
+      std::cout << "Node " << node.vi << ":" << endl;
+      std::cout << node.displacements[0] << " " << node.displacements[0] << " "
+                << node.displacements[1] << " " << node.displacements[2] << " "
+                << node.displacements[3] << " " << node.displacements[4] << " "
+                << node.displacements[5] << std::endl;
+    }
   }
 }
 
@@ -1331,7 +1340,7 @@ void FormFinder::applyDisplacements(
     updateNodeNormal(v, fixedVertices);
   }
   updateElementalFrames();
-  if (mShouldDraw) {
+  if (settings.shouldDraw || true) {
     mesh->updateEigenEdgeAndVertices();
   }
 }
@@ -1451,13 +1460,13 @@ void FormFinder::updatePositionsOnTheFly(
       //      assert(rotationalSumSk_J != 0);
     }
     mesh->nodes[v].translationalMass =
-        gamma * pow(Dtini, 2) * 2 * translationalSumSk;
+        gamma * pow(settings.Dtini, 2) * 2 * translationalSumSk;
     mesh->nodes[v].rotationalMass_I2 =
-        gamma * pow(Dtini, 2) * 8 * rotationalSumSk_I2;
+        gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I2;
     mesh->nodes[v].rotationalMass_I3 =
-        gamma * pow(Dtini, 2) * 8 * rotationalSumSk_I3;
+        gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I3;
     mesh->nodes[v].rotationalMass_J =
-        gamma * pow(Dtini, 2) * 8 * rotationalSumSk_J;
+        gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_J;
 
     //    assert(std::pow(Dtini, 2.0) * translationalSumSk /
     //               mesh->nodes[v].translationalMass <
@@ -1494,13 +1503,13 @@ void FormFinder::updatePositionsOnTheFly(
 
   for (VertexType &v : mesh->vert) {
     Node &node = mesh->nodes[v];
-    node.velocity = node.velocity + node.acceleration * Dt;
+    node.velocity = node.velocity + node.acceleration * settings.Dt;
   }
   updateKineticEnergy();
 
   for (VertexType &v : mesh->vert) {
     Node &node = mesh->nodes[v];
-    node.displacements = node.displacements + node.velocity * Dt;
+    node.displacements = node.displacements + node.velocity * settings.Dt;
   }
 
   for (VertexType &v : mesh->vert) {
@@ -1535,17 +1544,14 @@ void FormFinder::updatePositionsOnTheFly(
 }
 
 SimulationResults
-FormFinder::computeResults(const SimulationMesh &initialMesh) {
+FormFinder::computeResults(const std::shared_ptr<SimulationJob> &pJob) {
-  const size_t numberOfVertices = initialMesh.VN();
+  std::vector<Vector6d> displacements(mesh->VN());
-  std::vector<Vector6d> displacements(numberOfVertices);
+  for (size_t vi = 0; vi < mesh->VN(); vi++) {
-
-  for (size_t vi = 0; vi < numberOfVertices; vi++) {
     displacements[vi] = mesh->nodes[vi].displacements;
   }
-
   history.numberOfSteps = mCurrentSimulationStep;
 
-  return SimulationResults{history, displacements};
+  return SimulationResults{pJob, history, displacements};
 }
 
 void FormFinder::draw(const std::string &screenshotsFolder = {}) {
@@ -1667,8 +1673,8 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
   polyscope::getCurveNetwork(meshPolyscopeLabel)
       ->addNodeScalarQuantity("Residual force norm", residualForcesNorm)
       ->setEnabled(false);
-  polyscope::getCurveNetwork()->addNodeScalarQuantity("Node acceleration x",
+  polyscope::getCurveNetwork(meshPolyscopeLabel)
-                                                      accelerationX);
+      ->addNodeScalarQuantity("Node acceleration x", accelerationX);
 
   // Edge quantities
   std::vector<double> A(mesh->EN());
@@ -1700,9 +1706,9 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
                                // matching PopItemWidth() below.
 
     ImGui::InputInt("Simulation drawing step",
-                    &mDrawingStep); // set a int variable
+                    &settings.drawingStep); // set a int variable
     ImGui::Checkbox("Enable drawing",
-                    &mShouldDraw); // set a int variable
+                    &settings.shouldDraw); // set a int variable
     ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep);
 
     ImGui::PopItemWidth();
@@ -1727,13 +1733,16 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
   }
 }
 
-SimulationResults FormFinder::executeSimulation(
+SimulationResults
-    const SimulationJob &job, const bool &beVerbose, const bool &shouldDraw,
+FormFinder::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
-    const bool &shouldCreatePlots, const size_t &maxDRMIterations) {
+                              const bool &shouldDraw, const bool &beVerbose,
-  assert(job.mesh.operator bool());
+                              const bool &shouldCreatePlots,
+                              const size_t &maxDRMIterations) {
+  assert(pJob->pMesh.operator bool());
   auto t1 = std::chrono::high_resolution_clock::now();
   reset();
-  mShouldDraw = shouldDraw;
+  settings.shouldDraw = shouldDraw;
+  settings.beVerbose = beVerbose;
 
   //  if(!job.nodalExternalForces.empty()){
   //      double externalForcesNorm=0;
@@ -1744,29 +1753,29 @@ SimulationResults FormFinder::executeSimulation(
   //      setTotalResidualForcesNormThreshold(externalForcesNorm);
   //  }
 
-  constrainedVertices = job.fixedVertices;
+  constrainedVertices = pJob->constrainedVertices;
-  if (!job.nodalForcedDisplacements.empty()) {
+  if (!pJob->nodalForcedDisplacements.empty()) {
     for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair :
-         job.nodalForcedDisplacements) {
+         pJob->nodalForcedDisplacements) {
       const VertexIndex vi = viDisplPair.first;
       constrainedVertices[vi].insert({0, 1, 2});
     }
   }
-  if (!job.nodalForcedNormals.empty()) {
+  //  if (!pJob->nodalForcedNormals.empty()) {
-    for (std::pair<VertexIndex, Eigen::Vector3d> viNormalPair :
+  //    for (std::pair<VertexIndex, Eigen::Vector3d> viNormalPair :
-         job.nodalForcedDisplacements) {
+  //         pJob->nodalForcedDisplacements) {
-      assert(viNormalPair.second[2] >= 0);
+  //      assert(viNormalPair.second[2] >= 0);
-    }
+  //    }
-  }
+  //  }
 
-  mesh = std::make_unique<SimulationMesh>(*job.mesh);
+  mesh = std::make_unique<SimulationMesh>(*pJob->pMesh);
   if (beVerbose) {
     std::cout << "Executing simulation for mesh with:" << mesh->VN()
               << " nodes and " << mesh->EN() << " elements." << std::endl;
   }
   computeRigidSupports();
 
-  if (mShouldDraw || true) {
+  if (settings.shouldDraw || true) {
     if (!polyscope::state::initialized) {
       initPolyscope();
     }
@@ -1774,18 +1783,17 @@ SimulationResults FormFinder::executeSimulation(
         meshPolyscopeLabel, mesh->getEigenVertices(), mesh->getEigenEdges());
     //    registerWorldAxes();
   }
-  for (auto fixedVertex : job.fixedVertices) {
+  for (auto fixedVertex : pJob->constrainedVertices) {
     assert(fixedVertex.first < mesh->VN());
   }
 
   updateElementalFrames();
   updateNodalMasses();
-  if (!job.nodalForcedDisplacements.empty() &&
+  if (!pJob->nodalForcedDisplacements.empty() &&
-      job.nodalExternalForces.empty()) {
+      pJob->nodalExternalForces.empty()) {
     shouldApplyInitialDistortion = true;
   }
-  updateNodalExternalForces(job.nodalExternalForces, constrainedVertices);
+  updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
-
   while (maxDRMIterations == 0 || mCurrentSimulationStep < maxDRMIterations) {
     //  while (true) {
     updateNormalDerivatives();
@@ -1804,14 +1812,14 @@ SimulationResults FormFinder::executeSimulation(
     updateNodalVelocities();
     updateNodalDisplacements();
     applyDisplacements(constrainedVertices);
-    if (!job.nodalForcedDisplacements.empty()) {
+    if (!pJob->nodalForcedDisplacements.empty()) {
       applyForcedDisplacements(
 
-          job.nodalForcedDisplacements);
+          pJob->nodalForcedDisplacements);
-    }
-    if (!job.nodalForcedNormals.empty()) {
-      applyForcedNormals(job.nodalForcedNormals);
     }
+    //    if (!pJob->nodalForcedNormals.empty()) {
+    //      applyForcedNormals(pJob->nodalForcedNormals);
+    //    }
     updateElementalLengths();
 
     if (mCurrentSimulationStep != 0) {
@@ -1820,30 +1828,36 @@ SimulationResults FormFinder::executeSimulation(
     if (mCurrentSimulationStep > 100000) {
       shouldUseKineticEnergyThreshold = true;
     }
-    if (mCurrentSimulationStep>500000 ||(mShouldDraw  &&
+
-                                            mCurrentSimulationStep % mDrawingStep == 0) /* &&
+    if (mCurrentSimulationStep > 500000) {
-currentSimulationStep > maxDRMIterations*/) {
+      std::time_t now = time(0);
+      std::tm *ltm = std::localtime(&now);
+      std::string dir =
+          "../ProblematicSimulationJobs/" + std::to_string(ltm->tm_hour) + ":" +
+          std::to_string(ltm->tm_min) + "_" + std::to_string(ltm->tm_mday) +
+          "_" + std::to_string(1 + ltm->tm_mon) + "_" +
+          std::to_string(1900 + ltm->tm_year);
+      std::filesystem::create_directories(dir);
+      pJob->save(dir);
+
+      std::cout << "Non terminating simulation found. Saved simulation job to:"
+                << dir << std::endl;
+      std::cout << "Exiting.." << std::endl;
+      FormFinder debug;
+      debug.executeSimulation(pJob, true, true, true);
+      std::terminate();
+    }
+
+    if (settings.shouldDraw &&
+        mCurrentSimulationStep % settings.drawingStep == 0) /* &&
+currentSimulationStep > maxDRMIterations*/
+    {
       //      std::string saveTo = std::filesystem::current_path()
       //                               .append("Debugging_files")
       //                               .append("Screenshots")
       //                               .string();
       //      draw(saveTo);
       draw();
-      static bool wasExecuted = false;
-      if (mCurrentSimulationStep > 500000) {
-        FormFinder debug;
-        debug.executeSimulation(job, true, true, true);
-        wasExecuted = true;
-      }
-      std::cout << "Residual forces norm: " << mesh->totalResidualForcesNorm
-                << std::endl;
-      std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy
-                << std::endl;
-      auto yValues = history.residualForces;
-      auto xPlot = matplot::linspace(0, yValues.size(), yValues.size());
-      plotHandle = matplot::scatter(xPlot, yValues);
-      std::string label = "Residual forces";
-      plotHandle->legend_string(label);
 
       //      yValues = history.kineticEnergy;
       //      plotHandle = matplot::scatter(xPlot, yValues);
@@ -1852,6 +1866,22 @@ currentSimulationStep > maxDRMIterations*/) {
 
       //      shouldUseKineticEnergyThreshold = true;
     }
+    //    if (mCurrentSimulationStep % 100 == 0 && mCurrentSimulationStep >
+    //    100000) {
+    //      std::cout << "Current simulation step:" << mCurrentSimulationStep
+    //                << std::endl;
+    //      std::cout << "Residual forces norm: " <<
+    //      mesh->totalResidualForcesNorm
+    //                << std::endl;
+    //      std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy
+    //                << std::endl;
+    //      //      auto yValues = history.residualForces;
+    //      //      auto xPlot = matplot::linspace(0, yValues.size(),
+    //      yValues.size());
+    //      //      plotHandle = matplot::scatter(xPlot, yValues);
+    //      //      std::string label = "Residual forces";
+    //      //      plotHandle->legend_string(label);
+    //    }
     //    t = t + Dt;
     mCurrentSimulationStep++;
     //    std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy
@@ -1860,9 +1890,11 @@ currentSimulationStep > maxDRMIterations*/) {
     //              << std::endl;
     if (mesh->previousTotalKineticEnergy > mesh->currentTotalKineticEnergy) {
       if (/*mesh.totalPotentialEnergykN < 10 ||*/
-          mesh->totalResidualForcesNorm < totalResidualForcesNormThreshold ||
+          mesh->totalResidualForcesNorm <
+              settings.totalResidualForcesNormThreshold ||
           (shouldUseKineticEnergyThreshold &&
-           mesh->currentTotalKineticEnergy < totalKineticEnergyThreshold)) {
+           mesh->currentTotalKineticEnergy <
+               settings.totalKineticEnergyThreshold)) {
         if (beVerbose) {
           std::cout << "Convergence after " << mCurrentSimulationStep
                     << " steps" << std::endl;
@@ -1883,7 +1915,7 @@ currentSimulationStep > maxDRMIterations*/) {
         //        }
       }
       resetVelocities();
-      Dt = Dt * xi;
+      settings.Dt = settings.Dt * settings.xi;
     }
   }
   if (mCurrentSimulationStep == maxDRMIterations) {
@@ -1903,17 +1935,16 @@ currentSimulationStep > maxDRMIterations*/) {
               << "s" << std::endl;
   }
   //    mesh.printVertexCoordinates(mesh.VN() / 2);
-  if (mShouldDraw) {
+  if (settings.shouldDraw) {
     draw();
   }
   if (polyscope::hasCurveNetwork(meshPolyscopeLabel)) {
     polyscope::removeCurveNetwork(meshPolyscopeLabel);
   }
-  SimulationResults results = computeResults(*job.mesh);
+  SimulationResults results = computeResults(pJob);
-  results.label = job.mesh->getLabel();
   if (shouldCreatePlots) {
     SimulationResultsReporter reporter;
-    reporter.reportResults({results}, "Results", job.mesh->getLabel());
+    reporter.reportResults({results}, "Results", pJob->pMesh->getLabel());
   }
   return results;
 }

beamformfinder.hpp

@@ -22,19 +22,24 @@ struct DifferentiateWithRespectTo {
   const DoFType &dofi;
 };
 class FormFinder {
+  struct Settings {
+    bool shouldDraw{false};
+    int drawingStep{1};
+    const double totalKineticEnergyThreshold{1e-9};
+    double totalResidualForcesNormThreshold{1e-3};
+    const double Dtini{0.000001};
+    double Dt{Dtini};
+    const double xi{0.9969};
+    bool beVerbose;
+  };
+
+  Settings settings;
 
 private:
-  const double Dtini{0.1};
-  double Dt{Dtini};
-  const double xi{0.9969};
-  double totalResidualForcesNormThreshold{1e-5};
   bool shouldUseKineticEnergyThreshold{false};
   bool checkedForMaximumMoment{false};
-  const double totalKineticEnergyThreshold{1e-9};
   double externalMomentsNorm{0};
   size_t mCurrentSimulationStep{0};
-  int mDrawingStep{1};
-  bool mShouldDraw{false};
   matplot::line_handle plotHandle;
   std::vector<double> plotYValues;
 
@@ -88,7 +93,7 @@ private:
 
   void resetVelocities();
 
-  SimulationResults computeResults(const SimulationMesh &initialiMesh);
+  SimulationResults computeResults(const std::shared_ptr<SimulationJob> &pJob);
 
   void updateNodePosition(
       VertexType &v,
@@ -187,8 +192,9 @@ private:
 public:
   FormFinder();
   SimulationResults executeSimulation(
-      const SimulationJob &job, const bool &beVerbose = false,
+      const std::shared_ptr<SimulationJob> &pJob,
-      const bool &shouldDraw = false, const bool &createPlots = false,
+      const bool &shouldDraw = false, const bool &beVerbose = false,
+      const bool &createPlots = false,
       const size_t &maxDRMIterations = FormFinder::maxDRMIterations);
   inline static const size_t maxDRMIterations{0};
 

edgemesh.cpp

@@ -175,6 +175,7 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth,
 }
 
 bool VCGEdgeMesh::loadFromPly(const std::string plyFilename) {
+
   std::string usedPath = plyFilename;
   if (std::filesystem::path(plyFilename).is_relative()) {
     usedPath = std::filesystem::absolute(plyFilename).string();
@@ -200,6 +201,7 @@ bool VCGEdgeMesh::loadFromPly(const std::string plyFilename) {
   getEdges(eigenEdges);
   getVertices(eigenVertices);
   vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
+  label = std::filesystem::path(plyFilename).stem().string();
 
   std::cout << plyFilename << " was loaded successfuly." << std::endl;
   std::cout << "Mesh has " << EN() << " edges." << std::endl;
@@ -347,5 +349,8 @@ void VCGEdgeMesh::printVertexCoordinates(const size_t &vi) const {
 
 void VCGEdgeMesh::registerForDrawing() const {
   initPolyscope();
-  polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges());
+  const double drawingRadius = 0.0007;
+  polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges())
+      ->setRadius(drawingRadius, false);
+  std::cout << "Registered:" << label << std::endl;
 }

edgemesh.hpp

@@ -46,7 +46,7 @@ public:
     return vcg::tri::Index<VCGEdgeMesh>(*this, meshElement);
   }
   void updateEigenEdgeAndVertices();
-  void copy(VCGEdgeMesh &mesh);
+  virtual void copy(VCGEdgeMesh &mesh);
 
   void getEdges(Eigen::MatrixX3d &edgeStartingPoints,
                 Eigen::MatrixX3d &edgeEndingPoints) const;

elementalmesh.cpp

@@ -1,5 +1,12 @@
 #include "elementalmesh.hpp"
 
+SimulationMesh::SimulationMesh() {
+  elements = vcg::tri::Allocator<VCGEdgeMesh>::GetPerEdgeAttribute<Element>(
+      *this, std::string("Elements"));
+  nodes = vcg::tri::Allocator<VCGEdgeMesh>::GetPerVertexAttribute<Node>(
+      *this, std::string("Nodes"));
+}
+
 SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) {
   vcg::tri::MeshAssert<VCGEdgeMesh>::VertexNormalNormalized(mesh);
   //  bool containsNormals = true;
@@ -85,7 +92,6 @@ void SimulationMesh::initializeNodes() {
     Node &node = nodes[v];
     node.vi = vi;
     node.initialLocation = v.cP();
-    node.previousLocation = v.cP();
     node.initialNormal = v.cN();
     node.derivativeOfNormal.resize(6, VectorType(0, 0, 0));
 
@@ -126,6 +132,52 @@ void SimulationMesh::initializeNodes() {
   }
 }
 
+void SimulationMesh::reset() {
+  for (const EdgeType &e : edge) {
+    Element &element = elements[e];
+    element.ei = getIndex(e);
+    const VCGEdgeMesh::CoordType p0 = e.cP(0);
+    const VCGEdgeMesh::CoordType p1 = e.cP(1);
+    const vcg::Segment3<double> s(p0, p1);
+    element.initialLength = s.Length();
+    element.length = element.initialLength;
+    element.updateConstFactors();
+  }
+
+  for (const VertexType &v : vert) {
+    Node &node = nodes[v];
+    node.vi = getIndex(v);
+    node.initialLocation = v.cP();
+    node.initialNormal = v.cN();
+    node.derivativeOfNormal.resize(6, VectorType(0, 0, 0));
+
+    node.displacements[3] =
+        v.cN()[0]; // initialize nx diplacement with vertex normal x
+                   //        component
+    node.displacements[4] =
+        v.cN()[1]; // initialize ny(in the paper) diplacement with vertex
+    //        normal
+    // y component.
+
+    const EdgeType &referenceElement = *getReferenceElement(v);
+    const VectorType t01 =
+        computeT1Vector(referenceElement.cP(0), referenceElement.cP(1));
+    const VectorType f01 = (t01 - (v.cN() * (t01.dot(v.cN())))).Normalize();
+
+    for (const VCGEdgeMesh::EdgePointer &ep : nodes[v].incidentElements) {
+      assert(referenceElement.cV(0) == ep->cV(0) ||
+             referenceElement.cV(0) == ep->cV(1) ||
+             referenceElement.cV(1) == ep->cV(0) ||
+             referenceElement.cV(1) == ep->cV(1));
+      const VectorType t1 = computeT1Vector(*ep);
+      const VectorType f1 = t1 - (v.cN() * (t1.dot(v.cN()))).Normalize();
+      const EdgeIndex ei = getIndex(ep);
+      const double alphaAngle = computeAngle(f01, f1, v.cN());
+      node.alphaAngles[ei] = alphaAngle;
+    }
+  }
+}
+
 void SimulationMesh::initializeElements() {
   computeElementalProperties();
   for (const EdgeType &e : edge) {
@@ -208,45 +260,93 @@ std::vector<ElementMaterial> SimulationMesh::getBeamMaterial() {
   return beamMaterial;
 }
 
-bool SimulationMesh::loadFromPly(const string &plyFilename) {
+bool SimulationMesh::loadPly(const string &plyFilename) {
   this->Clear();
   //  assert(plyFileHasAllRequiredFields(plyFilename));
   // Load the ply file
-  VCGEdgeMesh::PerEdgeAttributeHandle<CrossSectionType> handleBeamDimensions;
+  VCGEdgeMesh::PerEdgeAttributeHandle<CrossSectionType> handleBeamDimensions =
-  VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial> handleBeamMaterial;
+      vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<
-  handleBeamDimensions =
+          CrossSectionType>(*this, plyPropertyBeamDimensionsID);
-      vcg::tri::Allocator<VCGEdgeMesh>::AddPerEdgeAttribute<CrossSectionType>(
+  VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial> handleBeamMaterial =
-          *this, plyPropertyBeamDimensionsID);
+      vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<ElementMaterial>(
-  handleBeamMaterial =
-      vcg::tri::Allocator<VCGEdgeMesh>::AddPerEdgeAttribute<ElementMaterial>(
           *this, plyPropertyBeamMaterialID);
-  nanoply::NanoPlyWrapper<VCGEdgeMesh>::CustomAttributeDescriptor customAttrib;
+  nanoply::NanoPlyWrapper<SimulationMesh>::CustomAttributeDescriptor
+      customAttrib;
   customAttrib.GetMeshAttrib(plyFilename);
   customAttrib.AddEdgeAttribDescriptor<CrossSectionType, float, 2>(
-      plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr);
+      plyPropertyBeamDimensionsID, nanoply::NNP_LIST_INT8_FLOAT32, nullptr);
+  /*FIXME: Since I allow CrossSectionType to take two types I should export the
+   * type as well such that that when loaded the correct type of cross section
+   * is used.
+   */
   customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>(
-      plyPropertyBeamMaterialID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr);
+      plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT32, nullptr);
+
+  VCGEdgeMesh::PerEdgeAttributeHandle<std::array<double, 6>>
+      handleBeamProperties =
+          vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<
+              std::array<double, 6>>(*this, plyPropertyBeamProperties);
+  customAttrib.AddEdgeAttribDescriptor<std::array<double, 6>, double, 6>(
+      plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
+
+  //  VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial>
+  //  handleBeamRigidityContants;
+  //  customAttrib.AddEdgeAttribDescriptor<vcg::Point4f, float, 4>(
+  //      plyPropertyBeamRigidityConstantsID, nanoply::NNP_LIST_INT8_FLOAT32,
+  //      nullptr);
   unsigned int mask = 0;
   mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
   mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
   mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
   mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
-  if (nanoply::NanoPlyWrapper<SimulationMesh>::LoadModel(plyFilename.c_str(),
+  if (nanoply::NanoPlyWrapper<SimulationMesh>::LoadModel(
-                                                         *this, mask) != 0) {
+          plyFilename.c_str(), *this, mask, customAttrib) != 0) {
     return false;
   }
+
+  elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(
+      *this, std::string("Elements"));
+  nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(
+      *this, std::string("Nodes"));
+  vcg::tri::UpdateTopology<SimulationMesh>::VertexEdge(*this);
+  initializeNodes();
+  initializeElements();
+  updateEigenEdgeAndVertices();
+
+  if (!handleBeamProperties._handle->data.empty()) {
+    for (size_t ei = 0; ei < EN(); ei++) {
+      elements[ei].properties = Element::Properties(handleBeamProperties[ei]);
+      elements[ei].updateConstFactors();
+    }
+  }
+
   return true;
 }
 
 bool SimulationMesh::savePly(const std::string &plyFilename) {
   nanoply::NanoPlyWrapper<VCGEdgeMesh>::CustomAttributeDescriptor customAttrib;
   customAttrib.GetMeshAttrib(plyFilename);
+
+  dimensions = getBeamDimensions();
   customAttrib.AddEdgeAttribDescriptor<CrossSectionType, float, 2>(
-      plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32,
+      plyPropertyBeamDimensionsID, nanoply::NNP_LIST_INT8_FLOAT32,
-      getBeamDimensions().data());
+      dimensions.data());
+  material = getBeamMaterial();
   customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>(
-      plyPropertyBeamMaterialID, nanoply::NNP_LIST_UINT8_FLOAT32,
+      plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT32,
-      getBeamMaterial().data());
+      material.data());
+  std::vector<std::array<double, 6>> beamProperties(EN());
+  for (size_t ei = 0; ei < EN(); ei++) {
+    auto props = elements[ei].properties.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;
@@ -254,7 +354,7 @@ bool SimulationMesh::savePly(const std::string &plyFilename) {
   mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
   mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
   if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::SaveModel(
-          plyFilename.c_str(), *this, mask, customAttrib, false) != 0) {
+          plyFilename.c_str(), *this, mask, customAttrib, false) != 1) {
     return false;
   }
   return true;
@@ -354,12 +454,24 @@ Element::Properties::Properties(const CrossSectionType &dimensions,
   computeDimensionsProperties(dimensions);
   computeMaterialProperties(material);
 }
-void Element::updateConstFactors() {
 
+Element::Properties::Properties(const std::array<double, 6> &arr) {
+  assert(arr.size() == 6);
+  E = arr[0];
+  G = arr[1];
+  A = arr[2];
+  I2 = arr[3];
+  I3 = arr[4];
+  J = arr[5];
+}
+
+std::array<double, 6> Element::Properties::toArray() const {
+  return std::array<double, 6>({E, G, A, I2, I3, J});
+}
+
+void Element::updateConstFactors() {
   axialConstFactor = properties.E * properties.A / initialLength;
   torsionConstFactor = properties.G * properties.J / initialLength;
   firstBendingConstFactor = 2 * properties.E * properties.I2 / initialLength;
   secondBendingConstFactor = 2 * properties.E * properties.I3 / initialLength;
-  int i = 0;
-  i++;
 }

elementalmesh.hpp

@@ -7,18 +7,21 @@
 
 struct Element;
 struct Node;
-// using CrossSectionType = RectangularBeamDimensions;
+using CrossSectionType = RectangularBeamDimensions;
-using CrossSectionType = CylindricalBeamDimensions;
+// using CrossSectionType = CylindricalBeamDimensions;
 
 class SimulationMesh : public VCGEdgeMesh {
 private:
   void computeElementalProperties();
-  void initializeNodes();
   void initializeElements();
+  void initializeNodes();
   EdgePointer getReferenceElement(const VertexType &v);
 
   const std::string plyPropertyBeamDimensionsID{"beam_dimensions"};
   const std::string plyPropertyBeamMaterialID{"beam_material"};
+  const std::string plyPropertyBeamProperties{"beam_properties"};
+  std::vector<ElementMaterial> material;
+  std::vector<CrossSectionType> dimensions;
 
 public:
   PerEdgeAttributeHandle<Element> elements;
@@ -30,7 +33,7 @@ public:
 
   std::vector<VCGEdgeMesh::EdgePointer>
   getIncidentElements(const VCGEdgeMesh::VertexType &v);
-  bool loadFromPly(const string &plyFilename);
+  bool loadPly(const string &plyFilename);
   std::vector<CrossSectionType> getBeamDimensions();
   std::vector<ElementMaterial> getBeamMaterial();
 
@@ -42,6 +45,8 @@ public:
   bool savePly(const std::string &plyFilename);
   void setBeamCrossSection(const CrossSectionType &beamDimensions);
   void setBeamMaterial(const double &pr, const double &ym);
+  void reset();
+  SimulationMesh();
 };
 
 struct Element {
@@ -63,7 +68,11 @@ struct Element {
     void setMaterial(const ElementMaterial &material);
     Properties(const CrossSectionType &dimensions,
                const ElementMaterial &material);
+    Properties(const std::array<double, 6> &arr);
     Properties() {}
+
+  public:
+    std::array<double, 6> toArray() const;
   };
 
   struct LocalFrame {

flatpattern.cpp

@@ -4,7 +4,7 @@
 
 FlatPattern::FlatPattern() {}
 
-FlatPattern::FlatPattern(const string &filename, bool addNormalsIfAbsent) {
+FlatPattern::FlatPattern(const std::string &filename, bool addNormalsIfAbsent) {
   assert(std::filesystem::exists(std::filesystem::path(filename)));
   loadFromPly(filename);
   if (addNormalsIfAbsent) {
@@ -17,7 +17,8 @@ FlatPattern::FlatPattern(const string &filename, bool addNormalsIfAbsent) {
   }
 
   vcg::tri::UpdateTopology<FlatPattern>::VertexEdge(*this);
-  scale();
+  const double baseTriangleCentralEdgeSize =
+      (vert[0].cP() - vert[3].cP()).Norm();
 
   updateEigenEdgeAndVertices();
 }
@@ -32,7 +33,7 @@ FlatPattern::FlatPattern(const std::vector<size_t> &numberOfNodesPerSlot,
       v.N() = CoordType(0, 0, 1);
     }
   }
-  scale();
+  baseTriangleHeight = (vert[0].cP() - vert[3].cP()).Norm();
   updateEigenEdgeAndVertices();
 }
 
@@ -84,6 +85,11 @@ bool FlatPattern::createHoneycombAtom() {
   return true;
 }
 
+void FlatPattern::copy(FlatPattern &pattern) {
+  VCGEdgeMesh::copy(pattern);
+  baseTriangleHeight = pattern.getBaseTriangleHeight();
+}
+
 void FlatPattern::deleteDanglingEdges() {
   for (VertexType &v : vert) {
     std::vector<VCGEdgeMesh::EdgePointer> incidentElements;
@@ -101,7 +107,8 @@ void FlatPattern::deleteDanglingEdges() {
   vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
 }
 
-void FlatPattern::scale() {
+void FlatPattern::scale(const double &desiredBaseTriangleCentralEdgeSize) {
+  this->baseTriangleHeight = desiredBaseTriangleCentralEdgeSize;
   const double baseTriangleCentralEdgeSize =
       (vert[0].cP() - vert[3].cP()).Norm();
   const double scaleRatio =
@@ -296,6 +303,8 @@ void FlatPattern::removeDuplicateVertices() {
   vcg::tri::UpdateTopology<FlatPattern>::VertexEdge(*this);
 }
 
+double FlatPattern::getBaseTriangleHeight() const { return baseTriangleHeight; }
+
 void FlatPattern::tilePattern(VCGEdgeMesh &pattern, VCGTriMesh &tileInto) {
 
   const size_t middleIndex = 3;

flatpattern.hpp

@@ -11,6 +11,7 @@ public:
               const std::vector<vcg::Point2i> &edges);
 
   bool createHoneycombAtom();
+  void copy(FlatPattern &pattern);
 
   void tilePattern(VCGEdgeMesh &pattern, VCGTriMesh &tileInto);
   void tilePattern(VCGEdgeMesh &pattern, VCGPolyMesh &tileInto,
@@ -21,13 +22,14 @@ public:
   void deleteDanglingVertices(
       vcg::tri::Allocator<FlatPattern>::PointerUpdater<VertexPointer> &pu);
   void deleteDanglingVertices();
+  void scale(const double &desiredBaseTriangleCentralEdgeSize);
+
+  double getBaseTriangleHeight() const;
 
 private:
   void deleteDanglingEdges();
   void removeDuplicateVertices();
-  void scale();
+  double baseTriangleHeight;
-
-  const double desiredBaseTriangleCentralEdgeSize{0.03};
 };
 
 #endif // FLATPATTERN_HPP

simulationhistoryplotter.hpp

@@ -59,11 +59,9 @@ struct SimulationResultsReporter {
     file.close();
   }
 
-  void reportResults(
+  void reportResults(const std::vector<SimulationResults> &results,
-      const std::vector<SimulationResults> &results,
+                     const std::string &reportFolderPath,
-      const std::string &reportFolderPath,
+                     const std::string &graphSuffix = std::string()) {
-      const std::string &graphSuffix = std::string(),
-      const SimulationResults &groundOfTruthResults = SimulationResults()) {
     if (results.empty()) {
       return;
     }
@@ -73,7 +71,7 @@ struct SimulationResultsReporter {
     for (const SimulationResults &simulationResult : results) {
       const auto simulationResultPath =
           std::filesystem::path(reportFolderPath)
-              .append(simulationResult.label);
+              .append(simulationResult.getLabel());
       std::filesystem::create_directory(simulationResultPath.string());
 
       createPlots(simulationResult.history, simulationResultPath.string(),

simulationresult.hpp

@@ -2,6 +2,7 @@
 #define SIMULATIONHISTORY_HPP
 
 #include "elementalmesh.hpp"
+#include "nlohmann/json.hpp"
 
 struct SimulationHistory {
   SimulationHistory() {}
@@ -31,23 +32,85 @@ struct SimulationHistory {
   }
 };
 
-struct SimulationJob {
+namespace nlohmann {
-  std::shared_ptr<SimulationMesh> mesh;
+
-  std::unordered_map<VertexIndex, std::unordered_set<int>> fixedVertices;
+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
+  inline static std::string jsonLabel_meshFilename{"mesh filename"};
+  inline static std::string jsonLabel_constrainedVertices{"fixed vertices"};
+  inline static std::string jsonLabel_nodalForces{"forces"};
+
+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;
-  std::unordered_map<VertexIndex, VectorType> nodalForcedNormals;
+  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; }
 
-  void registerForDrawing() const {
+  std::string toString() const {
+    nlohmann::json json;
+    if (!constrainedVertices.empty()) {
+      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[jsonLabel_nodalForces] = arrForces;
+    }
+
+    return json.dump();
+  }
+
+  void registerForDrawing(const std::string &meshLabel) const {
     initPolyscope();
-    const std::string polyscopeName = mesh->getLabel() + "_Simulation Job";
+    if (meshLabel.empty()) {
-    polyscope::registerCurveNetwork(polyscopeName, mesh->getEigenVertices(),
+      assert(false);
-                                    mesh->getEigenEdges());
+      std::cerr << "Expects a mesh label on which to draw the simulation job."
-    //  polyscope::registerPointCloud("Undeformed edge mesh PC",
+                << std::endl;
-    //                                job.edgeMesh.getEigenVertices());
+      return;
+    }
 
-    std::vector<std::array<double, 3>> nodeColors(mesh->VN());
+    auto structs = polyscope::state::structures;
-    for (auto fixedVertex : fixedVertices) {
+    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()) {
@@ -75,15 +138,13 @@ struct SimulationJob {
                   });
 
     if (!nodeColors.empty()) {
-      polyscope::getCurveNetwork(polyscopeName)
+      polyscope::getCurveNetwork(meshLabel)
-          ->addNodeColorQuantity("Boundary conditions", nodeColors);
+          ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors)
-      polyscope::getCurveNetwork(polyscopeName)
+          ->setEnabled(false);
-          ->getQuantity("Boundary conditions")
-          ->setEnabled(true);
     }
 
     // per node external forces
-    std::vector<std::array<double, 3>> externalForces(mesh->VN());
+    std::vector<std::array<double, 3>> externalForces(pMesh->VN());
     for (const auto &forcePair : nodalExternalForces) {
       auto index = forcePair.first;
       auto force = forcePair.second;
@@ -91,13 +152,99 @@ struct SimulationJob {
     }
 
     if (!externalForces.empty()) {
-      polyscope::getCurveNetwork(polyscopeName)
+      polyscope::getCurveNetwork(meshLabel)
-          ->addNodeVectorQuantity("External force", externalForces);
+          ->addNodeVectorQuantity("External force_" + label, externalForces)
-      polyscope::getCurveNetwork(polyscopeName)
+          ->setEnabled(false);
-          ->getQuantity("External force")
-          ->setEnabled(true);
     }
   }
+
+  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(jsonLabel_meshFilename)) {
+      pMesh->loadPly(json[jsonLabel_meshFilename]);
+    }
+
+    if (json.contains(jsonLabel_constrainedVertices)) {
+      constrainedVertices =
+          //      auto conV =
+          std::unordered_map<VertexIndex, std::unordered_set<int>>(
+              json[jsonLabel_constrainedVertices]);
+      std::cout << "Loaded constrained vertices. Number of constrained "
+                   "vertices found:"
+                << constrainedVertices.size() << std::endl;
+    }
+
+    if (json.contains(jsonLabel_nodalForces)) {
+      auto f = std::unordered_map<VertexIndex, std::array<double, 6>>(
+          json[jsonLabel_nodalForces]);
+      for (const auto &forces : f) {
+        nodalExternalForces[forces.first] = Vector6d(forces.second);
+      }
+      std::cout << "Loaded forces. Number of forces found:"
+                << nodalExternalForces.size() << std::endl;
+    }
+
+    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");
+    returnValue = pMesh->savePly(meshFilename);
+    nlohmann::json json;
+    json[jsonLabel_meshFilename] = meshFilename;
+    if (!constrainedVertices.empty()) {
+      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[jsonLabel_nodalForces] = arrForces;
+    }
+
+    std::string jsonFilename(
+        std::filesystem::path(pathFolderDirectory)
+            .append(pMesh->getLabel() + "_simScenario.json"));
+    std::ofstream jsonFile(jsonFilename);
+    jsonFile << json;
+    std::cout << "Saved simulation job as:" << jsonFilename << std::endl;
+
+    return returnValue;
+  }
 };
 
 namespace Eigen {
@@ -125,12 +272,30 @@ void read_binary(const std::string &filename, Matrix &matrix) {
 } // namespace Eigen
 
 struct SimulationResults {
+  std::shared_ptr<SimulationJob> job;
   SimulationHistory history;
   std::vector<Vector6d> displacements;
   double executionTime{0};
-  std::string label{"NoLabel"};
+  std::string labelPrefix{"deformed"};
+  inline static char deliminator{' '};
 
-  void registerForDrawing(const SimulationJob &job) const {
+  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";
@@ -140,87 +305,41 @@ struct SimulationResults {
     } /* else {
        polyscope::removeAllStructures();
      }*/
-    const std::string undeformedMeshName = "Undeformed_" + label;
+    //    const std::string undeformedMeshName = "Undeformed_" + label;
-    polyscope::registerCurveNetwork(undeformedMeshName,
+    //    polyscope::registerCurveNetwork(
-                                    job.mesh->getEigenVertices(),
+    //        undeformedMeshName, mesh->getEigenVertices(),
-                                    job.mesh->getEigenEdges());
+    //        mesh->getEigenEdges());
 
-    const std::string deformedMeshName = "Deformed_" + label;
+    const std::shared_ptr<SimulationMesh> &mesh = job->pMesh;
-    polyscope::registerCurveNetwork(deformedMeshName,
+    polyscope::registerCurveNetwork(getLabel(), mesh->getEigenVertices(),
-                                    job.mesh->getEigenVertices(),
+                                    mesh->getEigenEdges())
-                                    job.mesh->getEigenEdges());
+        ->setRadius(0.0007, false);
-    Eigen::MatrixX3d nodalDisplacements(job.mesh->VN(), 3);
+    Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3);
-    for (VertexIndex vi = 0; vi < job.mesh->VN(); vi++) {
+    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(deformedMeshName)
+    polyscope::getCurveNetwork(getLabel())
-        ->updateNodePositions(job.mesh->getEigenVertices() +
+        ->updateNodePositions(mesh->getEigenVertices() + nodalDisplacements);
-                              nodalDisplacements);
 
-    std::vector<std::array<double, 3>> nodeColors(job.mesh->VN());
+    job->registerForDrawing(getLabel());
-    for (auto fixedVertex : job.fixedVertices) {
+  }
-      nodeColors[fixedVertex.first] = {0, 0, 1};
+
+  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]);
     }
-    if (!job.nodalForcedDisplacements.empty()) {
+    m.savePly(getLabel() + ".ply");
-      for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair :
-           job.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;
-                    }
-                  });
-    // per node external forces
-    std::vector<std::array<double, 3>> externalForces(job.mesh->VN());
-    std::vector<std::array<double, 3>> externalMoments(job.mesh->VN());
-    for (const auto &forcePair : job.nodalExternalForces) {
-      auto index = forcePair.first;
-      auto force = forcePair.second;
-      externalForces[index] = {force[0], force[1], force[2]};
-      externalMoments[index] = {force[3], force[4], 0};
-    }
-
-    polyscope::getCurveNetwork(undeformedMeshName)
-        ->addNodeColorQuantity("Boundary conditions", nodeColors)
-        ->setEnabled(true);
-
-    polyscope::getCurveNetwork(undeformedMeshName)
-        ->addNodeVectorQuantity("External force", externalForces)
-        ->setEnabled(true);
-
-    polyscope::getCurveNetwork(undeformedMeshName)
-        ->addNodeVectorQuantity("External moments", externalMoments)
-        ->setEnabled(true);
-
-    polyscope::getCurveNetwork(deformedMeshName)
-        ->addNodeColorQuantity("Boundary conditions", nodeColors)
-        ->setEnabled(false);
-
-    polyscope::getCurveNetwork(deformedMeshName)
-        ->addNodeVectorQuantity("External force", externalForces)
-        ->setEnabled(true);
-
-    //    polyscope::show();
   }
 
   void save() {
-    const std::string filename(label + "_displacements.eigenBin");
+    const std::string filename(getLabel() + "_displacements.eigenBin");
 
     Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements);
     Eigen::write_binary(filename, m);

trianglepatterngeometry.cpp

@@ -169,6 +169,7 @@ void FlatPatternGeometry::add(const std::vector<vcg::Point3d> &vertices) {
   });
   vcg::tri::UpdateTopology<FlatPatternGeometry>::VertexEdge(*this);
   vcg::tri::UpdateTopology<FlatPatternGeometry>::EdgeEdge(*this);
+  updateEigenEdgeAndVertices();
 }
 
 void FlatPatternGeometry::add(const std::vector<vcg::Point2i> &edges) {
@@ -177,12 +178,14 @@ void FlatPatternGeometry::add(const std::vector<vcg::Point2i> &edges) {
   });
   vcg::tri::UpdateTopology<FlatPatternGeometry>::VertexEdge(*this);
   vcg::tri::UpdateTopology<FlatPatternGeometry>::EdgeEdge(*this);
+  updateEigenEdgeAndVertices();
 }
 
 void FlatPatternGeometry::add(const std::vector<vcg::Point3d> &vertices,
                               const std::vector<vcg::Point2i> &edges) {
   add(vertices);
   add(edges);
+  updateEigenEdgeAndVertices();
 }
 
 void FlatPatternGeometry::add(const std::vector<size_t> &numberOfNodesPerSlot,

utilities.hpp

@@ -19,6 +19,8 @@ struct Vector6d : public std::array<double, 6> {
     }
   }
 
+  Vector6d(const std::array<double, 6> &arr) : std::array<double, 6>(arr) {}
+
   Vector6d(const std::initializer_list<double> &initList) {
     std::copy(initList.begin(), initList.end(), this->begin());
   }
@@ -155,78 +157,6 @@ inline Eigen::MatrixXd toEigenMatrix(const std::vector<Vector6d> &v) {
 
 // namespace ConfigurationFile {
 
-// inline void getPlyFilename(const std::string jsonFilepath,
-//                           std::string &plyFilename) {
-//  std::ifstream inFile(jsonFilepath);
-//  std::string jsonContents((std::istreambuf_iterator<char>(inFile)),
-//                           std::istreambuf_iterator<char>());
-//  nlohmann::json jsonFile(nlohmann::json::parse(jsonContents));
-
-//  if (jsonFile.contains("plyFilename")) {
-//    plyFilename = jsonFile["plyFilename"];
-//  }
-//}
-
-// inline void
-// getNodalForces(const std::string jsonFilepath,
-//               std::unordered_map<size_t, Eigen::Vector3d> &nodalForces) {
-//  std::ifstream inFile(jsonFilepath);
-//  std::string jsonContents((std::istreambuf_iterator<char>(inFile)),
-//                           std::istreambuf_iterator<char>());
-//  nlohmann::json jsonFile(nlohmann::json::parse(jsonContents));
-
-//  nodalForces.clear();
-//  if (jsonFile.contains("forces")) {
-//    std::vector<std::vector<double>> forces = jsonFile["forces"];
-//    for (size_t forceIndex = 0; forceIndex < forces.size(); forceIndex++) {
-//      const BeamFormFinder::NodalForce nf{
-//          static_cast<size_t>(forces[forceIndex][0]),
-//          static_cast<size_t>(forces[forceIndex][1]), forces[forceIndex][2]};
-//      const size_t vertexIndex = forces[forceIndex][0];
-//      const Eigen::Vector3d forceVector(
-//          forces[forceIndex][1], forces[forceIndex][2],
-//          forces[forceIndex][3]);
-//      assert(forceIndex >= 0 && forceVector.norm() >= 0);
-//      nodalForces[vertexIndex] = forceVector;
-//    }
-//  }
-//}
-
-// inline void getFixedVertices(const std::string jsonFilepath,
-//                             std::vector<size_t> &fixedVertices) {
-//  std::ifstream inFile(jsonFilepath);
-//  std::string jsonContents((std::istreambuf_iterator<char>(inFile)),
-//                           std::istreambuf_iterator<char>());
-//  nlohmann::json jsonFile(nlohmann::json::parse(jsonContents));
-
-//  fixedVertices.clear();
-//  if (jsonFile.contains("fixedVertices")) {
-//    fixedVertices =
-//    static_cast<std::vector<size_t>>(jsonFile["fixedVertices"]);
-//  }
-//}
-
-// struct SimulationScenario {
-//  std::string edgeMeshFilename;
-//  std::vector<size_t> fixedVertices;
-//  std::unordered_map<size_t, Eigen::Vector3d> nodalForces;
-//};
-
-// void to_json(nlohmann::json &json, const SimulationScenario &scenario) {
-//  json["plyFilename"] = scenario.edgeMeshFilename;
-//  if (!scenario.fixedVertices.empty()) {
-//    json["fixedVertices"] = scenario.fixedVertices;
-//  }
-//  if (!scenario.nodalForces.empty()) {
-//    std::vector<std::tuple<size_t, double, double, double>> forces;
-//    std::transform(scenario.nodalForces.begin(), scenario.nodalForces.end(),
-//                   std::back_inserter(forces),
-//                   [](const std::pair<size_t, Eigen::Vector3d> &f) {
-//                     return std::tuple<size_t, double, double, double>{
-//                         f.first, f.second[0], f.second[1], f.second[2]};
-//                   });
-//    json["forces"] = forces;
-//  }
 //}
 //} // namespace ConfigurationFile
 #include "polyscope/curve_network.h"
@@ -263,7 +193,7 @@ inline void registerWorldAxes() {
 
   const std::string worldAxesName = "World Axes";
   polyscope::registerCurveNetwork(worldAxesName, axesPositions, axesEdges);
-  polyscope::getCurveNetwork(worldAxesName)->setRadius(0.001);
+  polyscope::getCurveNetwork(worldAxesName)->setRadius(0.0001, false);
   const std::string worldAxesColorName = worldAxesName + " Color";
   polyscope::getCurveNetwork(worldAxesName)
       ->addEdgeColorQuantity(worldAxesColorName, axesColors)