Refactoring

beamformfinder.cpp

@@ -1814,7 +1814,8 @@ FormFinder::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
     shouldApplyInitialDistortion = true;
   }
   updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
-  while (maxDRMIterations == 0 || mCurrentSimulationStep < maxDRMIterations) {
+  while (settings.maxDRMIterations == 0 ||
+         mCurrentSimulationStep < settings.maxDRMIterations) {
     //  while (true) {
     updateNormalDerivatives();
     updateT1Derivatives();
@@ -1850,33 +1851,25 @@ FormFinder::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
     }
 
     if (std::isnan(pMesh->currentTotalKineticEnergy)) {
-      if (pMesh->elements[0].dimensions.b / pMesh->elements[0].dimensions.h <
-          10) {
       std::time_t now = time(0);
       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 =
+      const std::string subDir = std::filesystem::path(dir)
-            std::filesystem::path(dir)
+                                     .append(std::to_string(ltm->tm_hour) +
-                .append(std::to_string(ltm->tm_hour) + ":" +
+                                             ":" + std::to_string(ltm->tm_min))
-                        std::to_string(ltm->tm_min))
                                      .string();
       std::filesystem::create_directories(subDir);
       pJob->save(subDir);
-        //        std::cout
+      std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl;
-        //            << "Non terminating simulation found. Saved simulation job
+      FormFinder debug;
-        //            to : "
+      FormFinder::Settings settings;
-        //            << dir << std::endl;
+      settings.shouldDraw = true;
-        //        std::terminate();
+      settings.beVerbose = true;
-      }
+      debug.executeSimulation(pJob, settings);
-      //      std::cout << "Exiting.." << std::endl;
+      std::terminate();
-      //      FormFinder debug;
-      //      FormFinder::Settings settings;
-      //      settings.shouldDraw = true;
-      //      settings.beVerbose = true;
-      //      debug.executeSimulation(pJob, settings);
       break;
     }
 
@@ -1966,10 +1959,12 @@ mesh->currentTotalPotentialEnergykN*/
 
   SimulationResults results = computeResults(pJob);
 
-  if (mCurrentSimulationStep == maxDRMIterations) {
+  if (mCurrentSimulationStep == settings.maxDRMIterations &&
+      mCurrentSimulationStep != 0) {
     std::cout << "Did not reach equilibrium before reaching the maximum number "
                  "of DRM steps ("
-              << maxDRMIterations << "). Breaking simulation" << std::endl;
+              << settings.maxDRMIterations << "). Breaking simulation"
+              << std::endl;
     results.converged = false;
   } else if (std::isnan(pMesh->currentTotalKineticEnergy)) {
     results.converged = false;

beamformfinder.hpp

@@ -204,7 +204,6 @@ public:
   SimulationResults
   executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
                     const Settings &settings = Settings());
-  inline static const size_t maxDRMIterations{0};
 
   static void runUnitTests();
   void setTotalResidualForcesNormThreshold(double value);

elementalmesh.cpp

@@ -326,6 +326,17 @@ bool SimulationMesh::loadPly(const string &plyFilename) {
     elements[ei].updateRigidity();
   }
 
+  bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001;
+  if (normalsAreAbsent) {
+    CoordType normalVector(0, 0, 1);
+    std::cout << "Warning: Normals are missing from " << plyFilename
+              << ". Added normal vector:" << toString(normalVector)
+              << std::endl;
+    for (auto &v : vert) {
+      v.N() = normalVector;
+    }
+  }
+
   return true;
 }
 

flatpattern.cpp

@@ -93,9 +93,9 @@ bool FlatPattern::createHoneycombAtom() {
   return true;
 }
 
-void FlatPattern::copy(FlatPattern &pattern) {
+void FlatPattern::copy(FlatPattern &copyFrom) {
-  VCGEdgeMesh::copy(pattern);
+  VCGEdgeMesh::copy(copyFrom);
-  baseTriangleHeight = pattern.getBaseTriangleHeight();
+  baseTriangleHeight = copyFrom.getBaseTriangleHeight();
 }
 
 void FlatPattern::deleteDanglingEdges() {

flatpattern.hpp

@@ -11,7 +11,7 @@ public:
               const std::vector<vcg::Point2i> &edges);
 
   bool createHoneycombAtom();
-  void copy(FlatPattern &pattern);
+  void copy(FlatPattern &copyFrom);
 
   void tilePattern(VCGEdgeMesh &pattern, VCGTriMesh &tileInto);
   void tilePattern(VCGEdgeMesh &pattern, VCGPolyMesh &tileInto,

simulationresult.hpp

@@ -54,9 +54,15 @@ template <> struct adl_serializer<std::unordered_map<VertexIndex, Vector6d>> {
 class SimulationJob {
   //  const std::unordered_map<VertexIndex, VectorType> nodalForcedNormals;
   // json labels
-  inline static std::string jsonLabel_meshFilename{"mesh filename"};
+  struct JSONLabels {
-  inline static std::string jsonLabel_constrainedVertices{"fixed vertices"};
+    inline static std::string meshFilename{"mesh filename"};
-  inline static std::string jsonLabel_nodalForces{"forces"};
+    inline static std::string constrainedVertices{"fixed vertices"};
+    inline static std::string nodalForces{"forces"};
+    inline static std::string label{"label"};
+    inline static std::string meshLabel{
+        "label"}; // TODO: should be in the savePly function of the simulation
+                  // mesh class
+  } jsonLabels;
 
 public:
   std::shared_ptr<SimulationMesh> pMesh;
@@ -78,14 +84,14 @@ public:
   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();
@@ -137,6 +143,7 @@ public:
                     }
                   });
 
+    auto cn = polyscope::getCurveNetwork(meshLabel);
     if (!nodeColors.empty()) {
       polyscope::getCurveNetwork(meshLabel)
           ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors)
@@ -180,23 +187,23 @@ public:
     json << ifs;
 
     pMesh = std::make_shared<SimulationMesh>();
-    if (json.contains(jsonLabel_meshFilename)) {
+    if (json.contains(jsonLabels.meshFilename)) {
-      pMesh->loadPly(json[jsonLabel_meshFilename]);
+      pMesh->loadPly(json[jsonLabels.meshFilename]);
     }
 
-    if (json.contains(jsonLabel_constrainedVertices)) {
+    if (json.contains(jsonLabels.constrainedVertices)) {
       constrainedVertices =
           //      auto conV =
           std::unordered_map<VertexIndex, std::unordered_set<int>>(
-              json[jsonLabel_constrainedVertices]);
+              json[jsonLabels.constrainedVertices]);
       std::cout << "Loaded constrained vertices. Number of constrained "
                    "vertices found:"
                 << constrainedVertices.size() << std::endl;
     }
 
-    if (json.contains(jsonLabel_nodalForces)) {
+    if (json.contains(jsonLabels.nodalForces)) {
       auto f = std::unordered_map<VertexIndex, std::array<double, 6>>(
-          json[jsonLabel_nodalForces]);
+          json[jsonLabels.nodalForces]);
       for (const auto &forces : f) {
         nodalExternalForces[forces.first] = Vector6d(forces.second);
       }
@@ -204,6 +211,14 @@ public:
                 << 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;
   }
 
@@ -224,21 +239,28 @@ public:
             .append(pMesh->getLabel() + ".ply");
     returnValue = pMesh->savePly(meshFilename);
     nlohmann::json json;
-    json[jsonLabel_meshFilename] = meshFilename;
+    json[jsonLabels.meshFilename] = meshFilename;
     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;
+    }
+
+    if (!label.empty()) {
+      json[jsonLabels.label] = label;
+    }
+    if (!pMesh->getLabel().empty()) {
+      json[jsonLabels.meshLabel] = pMesh->getLabel();
     }
 
     std::string jsonFilename(
         std::filesystem::path(pathFolderDirectory)
-            .append(pMesh->getLabel() + "_simScenario.json"));
+            .append(label + "_" + pMesh->getLabel() + "_simulationJob.json"));
     std::ofstream jsonFile(jsonFilename);
     jsonFile << json;
     std::cout << "Saved simulation job as:" << jsonFilename << std::endl;

utilities.hpp

@@ -162,6 +162,13 @@ inline Eigen::MatrixXd toEigenMatrix(const std::vector<Vector6d> &v) {
 #include "polyscope/curve_network.h"
 #include "polyscope/polyscope.h"
 
+inline void deinitPolyscope() {
+  if (!polyscope::state::initialized) {
+    return;
+  }
+
+  polyscope::shutdown();
+}
 inline void initPolyscope() {
   if (polyscope::state::initialized) {
     return;
@@ -207,4 +214,9 @@ void constructInverseMap(const T1 &map, T2 &oppositeMap) {
   }
 }
 
+template <typename T> std::string toString(const T &v) {
+  return "(" + std::to_string(v[0]) + "," + std::to_string(v[1]) + "," +
+         std::to_string(v[2]) + ")";
+}
+
 #endif // UTILITIES_H