Added unit tests for form finder. Refactoring.

beamformfinder.cpp

@@ -8,6 +8,197 @@
 #include <execution>
 #include <omp.h>
 
+void FormFinder::runUnitTests() {
+  const std::filesystem::path groundOfTruthFolder{
+      "/home/iason/Coding/Libraries/MySources/formFinder_unitTestFiles"};
+
+  FormFinder formFinder;
+  formFinder.setTotalResidualForcesNormThreshold(1);
+
+  // First example of the paper
+  VCGEdgeMesh beam;
+  //  const size_t spanGridSize = 11;
+  //  mesh.createSpanGrid(spanGridSize);
+  beam.loadFromPly(
+      std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply"));
+  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};
+  std::unordered_map<VertexIndex, Vector6d> nodalForces{
+      {2, Vector6d({0, 1000, 1000, 0, 0, 0})}};
+  // Forced displacements
+  std::unordered_map<size_t, Eigen::Vector3d> nodalForcedDisplacements;
+  nodalForcedDisplacements.insert({beam.VN() - 1, Eigen::Vector3d(-0.2, 0, 0)});
+
+  SimulationJob beamSimulationJob{
+      std::make_shared<SimulationMesh>(beam),
+      //      SimulationJob::constructFixedVerticesSpanGrid(spanGridSize,
+      //      mesh.VN()),
+      fixedVertices, nodalForces, nodalForcedDisplacements};
+  beamSimulationJob.mesh->setBeamMaterial(0.3, 200);
+  assert(CrossSectionType::name == CylindricalBeamDimensions::name);
+
+  beamSimulationJob.mesh->setBeamCrossSection(CrossSectionType{0.03, 0.026});
+  SimulationResults simpleBeam_simulationResults =
+      formFinder.executeSimulation(beamSimulationJob, false, true);
+  simpleBeam_simulationResults.label = "SimpleBeam";
+  simpleBeam_simulationResults.save();
+  const std::string simpleBeamGroundOfTruthBinaryFilename =
+      std::filesystem::path(groundOfTruthFolder)
+          .append("SimpleBeam_displacements.eigenBin");
+  assert(std::filesystem::exists(
+      std::filesystem::path(simpleBeamGroundOfTruthBinaryFilename)));
+  Eigen::MatrixXd simpleBeam_groundOfTruthDisplacements;
+  Eigen::read_binary(simpleBeamGroundOfTruthBinaryFilename,
+                     simpleBeam_groundOfTruthDisplacements);
+  if (!simpleBeam_simulationResults.isEqual(
+          simpleBeam_groundOfTruthDisplacements)) {
+    std::cerr << "Error:Beam simulation produces wrong results." << std::endl;
+    return;
+  }
+
+  // Second example of the paper
+  VCGEdgeMesh shortSpanGrid;
+  //  const size_t spanGridSize = 11;
+  //  mesh.createSpanGrid(spanGridSize);
+  shortSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder)
+                                .append("shortSpanGridshell.ply")
+                                .string());
+
+  fixedVertices.clear();
+  //// Corner nodes
+  fixedVertices[0] = std::unordered_set<DoFType>{2};
+  fixedVertices[4] = std::unordered_set<DoFType>{2};
+  fixedVertices[16] = std::unordered_set<DoFType>{2};
+  fixedVertices[20] = std::unordered_set<DoFType>{2};
+  //// Center node
+  fixedVertices[10] = std::unordered_set<DoFType>{0, 1, 3, 4, 5};
+
+  nodalForcedDisplacements.clear();
+  nodalForcedDisplacements.insert({{0, Eigen::Vector3d(0.1, 0.1, 0)},
+                                   {4, Eigen::Vector3d(-0.1, 0.1, 0)},
+                                   {16, Eigen::Vector3d(0.1, -0.1, 0)},
+                                   {20, Eigen::Vector3d(-0.1, -0.1, 0)}});
+
+  SimulationJob shortSpanGridshellSimulationJob{
+      std::make_shared<SimulationMesh>(shortSpanGrid),
+      fixedVertices,
+      {},
+      nodalForcedDisplacements};
+  shortSpanGridshellSimulationJob.mesh->setBeamMaterial(0.3, 200);
+  assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
+  shortSpanGridshellSimulationJob.mesh->setBeamCrossSection(
+      CrossSectionType{0.03, 0.026});
+  SimulationResults shortSpanGridshellSimulationResults =
+      formFinder.executeSimulation(shortSpanGridshellSimulationJob, false,
+                                   false);
+  shortSpanGridshellSimulationResults.label = "ShortSpanGridshell";
+  shortSpanGridshellSimulationResults.save();
+
+  const std::string groundOfTruthBinaryFilename =
+      std::filesystem::path(groundOfTruthFolder)
+          .append("ShortSpanGridshell_displacements.eigenBin")
+          .string();
+  assert(std::filesystem::exists(
+      std::filesystem::path(groundOfTruthBinaryFilename)));
+  Eigen::MatrixXd shortSpanGridshell_groundOfTruthDisplacements;
+  Eigen::read_binary(groundOfTruthBinaryFilename,
+                     shortSpanGridshell_groundOfTruthDisplacements);
+  //  shortSpanGridshellSimulationResults.registerForDrawing(
+  //      shortSpanGridshellSimulationJob);
+  //  polyscope::show();
+  assert(shortSpanGridshellSimulationResults.isEqual(
+      shortSpanGridshell_groundOfTruthDisplacements));
+  if (!shortSpanGridshellSimulationResults.isEqual(
+          shortSpanGridshell_groundOfTruthDisplacements)) {
+    std::cerr << "Error:Short span simulation produces wrong results."
+              << std::endl;
+    return;
+  }
+
+  // Third example of the paper
+  VCGEdgeMesh longSpanGrid;
+  longSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder)
+                               .append("longSpanGridshell.ply")
+                               .string());
+  const size_t spanGridSize = 11;
+
+  fixedVertices.clear();
+  for (size_t vi = 0; vi < spanGridSize - 1; vi++) {
+    fixedVertices[vi] = {0, 2};
+  }
+  for (size_t vi = longSpanGrid.VN() - 1 - (spanGridSize - 2);
+       vi < longSpanGrid.VN(); vi++) {
+    fixedVertices[vi] = {0, 2};
+  }
+  for (size_t vi = spanGridSize - 1;
+       vi < longSpanGrid.VN() - 1 - (spanGridSize - 2) - spanGridSize + 1;
+       vi += spanGridSize + 1) {
+    fixedVertices[vi] = {1, 2};
+    fixedVertices[vi + spanGridSize] = {1, 2};
+  }
+
+  nodalForcedDisplacements.clear();
+  const size_t horizontalOffset = std::floor((spanGridSize - 2) / 2);
+  nodalForcedDisplacements.insert(
+      {{horizontalOffset, Eigen::Vector3d(0, 0.3, 0)},
+       {horizontalOffset + 1, Eigen::Vector3d(0, 0.3, 0)},
+       {spanGridSize * (spanGridSize + 1) - 2 + horizontalOffset,
+        Eigen::Vector3d(0, -0.3, 0)},
+       {spanGridSize * (spanGridSize + 1) - 2 + horizontalOffset + 1,
+        Eigen::Vector3d(0, -0.3, 0)},
+       {std::floor(spanGridSize / 2) * (spanGridSize + 1) - 2,
+        Eigen::Vector3d(0.3, 0, 0)},
+       {(std::floor(spanGridSize / 2) + 1) * (spanGridSize + 1) - 2,
+        Eigen::Vector3d(0.3, 0, 0)},
+       {std::floor(spanGridSize / 2) * (spanGridSize + 1) - 2 + spanGridSize,
+        Eigen::Vector3d(-0.3, 0, 0)},
+       {(std::floor(spanGridSize / 2) + 1) * (spanGridSize + 1) - 2 +
+            spanGridSize,
+        Eigen::Vector3d(-0.3, 0, 0)}});
+
+  SimulationJob longSpanGridshell_simulationJob{
+      std::make_shared<SimulationMesh>(longSpanGrid),
+      fixedVertices,
+      {},
+      nodalForcedDisplacements};
+  longSpanGridshell_simulationJob.mesh->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(
+      CrossSectionType{0.03, 0.026});
+  SimulationResults longSpanGridshell_simulationResults =
+      formFinder.executeSimulation(longSpanGridshell_simulationJob, false,
+                                   false);
+  longSpanGridshell_simulationResults.label = "LongSpanGridshell";
+  longSpanGridshell_simulationResults.save();
+
+  const std::string longSpan_groundOfTruthBinaryFilename =
+      std::filesystem::path(groundOfTruthFolder)
+          .append("LongSpanGridshell_displacements.eigenBin")
+          .string();
+  assert(std::filesystem::exists(
+      std::filesystem::path(longSpan_groundOfTruthBinaryFilename)));
+  Eigen::MatrixXd longSpanGridshell_groundOfTruthDisplacements;
+  Eigen::read_binary(longSpan_groundOfTruthBinaryFilename,
+                     longSpanGridshell_groundOfTruthDisplacements);
+  assert(longSpanGridshell_simulationResults.isEqual(
+      longSpanGridshell_groundOfTruthDisplacements));
+  if (!longSpanGridshell_simulationResults.isEqual(
+          longSpanGridshell_groundOfTruthDisplacements)) {
+    std::cerr << "Error:Long span simulation produces wrong results."
+              << std::endl;
+    return;
+  }
+
+  std::cout << "Form finder unit tests succesufully passed." << std::endl;
+
+  //  polyscope::show();
+}
+
 void FormFinder::setTotalResidualForcesNormThreshold(double value) {
   totalResidualForcesNormThreshold = value;
 }
@@ -21,6 +212,7 @@ void FormFinder::reset() {
   mesh.release();
   plotYValues.clear();
   plotHandle.reset();
+  checkedForMaximumMoment = false;
   shouldUseKineticEnergyThreshold = false;
   externalMomentsNorm = 0;
 }
@@ -791,6 +983,7 @@ void FormFinder::updateResidualForcesOnTheFly(
     const Vector6d &nodeResidualForce = force.residual;
     //    sumOfResidualForces = sumOfResidualForces + nodeResidualForce;
     const double residualForceNorm = nodeResidualForce.norm();
+    const bool shouldTrigger = std::isnan(residualForceNorm);
     totalResidualForcesNorm += residualForceNorm;
   }
   mesh->previousTotalResidualForcesNorm = mesh->totalResidualForcesNorm;
@@ -1094,16 +1287,15 @@ void FormFinder::updateNodeNormal(
      * algorithm if the kinetic energy of the system drops below the set
      * threshold.
      * */
-    static bool wasExecuted{false};
     const bool viHasMoments =
         node.force.external[3] != 0 || node.force.external[4] != 0;
-    if (!wasExecuted && viHasMoments) {
+    if (!checkedForMaximumMoment && viHasMoments) {
       shouldUseKineticEnergyThreshold = true;
       std::cout
           << "Maximum moment reached.The Kinetic energy of the system will "
              "be used as a convergence criterion"
           << std::endl;
-      wasExecuted = true;
+      checkedForMaximumMoment = true;
     }
 
   } else {
@@ -1359,92 +1551,17 @@ FormFinder::computeResults(const SimulationMesh &initialMesh) {
 void FormFinder::draw(const std::string &screenshotsFolder = {}) {
   // update positions
   //  polyscope::getCurveNetwork("Undeformed edge mesh")->setEnabled(false);
-  polyscope::getCurveNetwork("Deformed edge mesh")
+  polyscope::getCurveNetwork(meshPolyscopeLabel)
       ->updateNodePositions(mesh->getEigenVertices());
-  // Per node kinetic energies
+
+  // Vertex quantities
   std::vector<double> kineticEnergies(mesh->VN());
-  for (const VertexType &v : mesh->vert) {
-    kineticEnergies[mesh->getIndex(v)] = mesh->nodes[v].kineticEnergy;
-  }
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeScalarQuantity("Kinetic Energy", kineticEnergies);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("Kinetic Energy")
-      ->setEnabled(false);
-
-  // Per node normals
   std::vector<std::array<double, 3>> nodeNormals(mesh->VN());
-  for (const VertexType &v : mesh->vert) {
-    const VectorType n = v.cN();
-    nodeNormals[mesh->getIndex(v)] = {n[0], n[1], n[2]};
-  }
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeVectorQuantity("Node normals", nodeNormals)
-      ->setEnabled(true);
-
-  // per node external forces
-  //  std::vector<std::array<double, 3>> externalForces(mesh->VN());
-  //  for (const VertexType &v : mesh->vert) {
-  //    const Vector6d nodeForce =
-  //        mesh->nodes[v].force.external.value_or(Vector6d(0));
-  //    externalForces[mesh->getIndex(v)] = {nodeForce[0], nodeForce[1],
-  //                                         nodeForce[2]};
-  //  }
-  //  polyscope::getCurveNetwork("Deformed edge mesh")
-  //      ->addNodeVectorQuantity("External force", externalForces);
-  //  polyscope::getCurveNetwork("Deformed edge mesh")
-  //      ->getQuantity("External force")
-  //      ->setEnabled(true);
-
   std::vector<std::array<double, 3>> internalForces(mesh->VN());
   std::vector<std::array<double, 3>> externalForces(mesh->VN());
   std::vector<std::array<double, 3>> externalMoments(mesh->VN());
   std::vector<double> internalForcesNorm(mesh->VN());
-  for (const VertexType &v : mesh->vert) {
-    // per node internal forces
-    const Vector6d nodeforce = mesh->nodes[v].force.internal * (-1);
-    internalForces[mesh->getIndex(v)] = {nodeforce[0], nodeforce[1],
-                                         nodeforce[2]};
-    internalForcesNorm[mesh->getIndex(v)] = nodeforce.norm();
-    // External force
-    const Vector6d nodeExternalForce = mesh->nodes[v].force.external;
-    externalForces[mesh->getIndex(v)] = {
-        nodeExternalForce[0], nodeExternalForce[1], nodeExternalForce[2]};
-    externalMoments[mesh->getIndex(v)] = {nodeExternalForce[3],
-                                          nodeExternalForce[4], 0};
-  }
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeVectorQuantity("Internal force", internalForces);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("Internal force")
-      ->setEnabled(false);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeVectorQuantity("External force", externalForces);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("External force")
-      ->setEnabled(true);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeVectorQuantity("External Moments", externalMoments)
-      ->setEnabled(true);
-
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeScalarQuantity("Internal force norm", internalForcesNorm);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("Internal force norm")
-      ->setEnabled(true);
-  // per node internal forces
   std::vector<std::array<double, 3>> internalAxialForces(mesh->VN());
-  for (const VertexType &v : mesh->vert) {
-    const Vector6d nodeforce = mesh->nodes[v].force.internalAxial * (-1);
-    internalAxialForces[mesh->getIndex(v)] = {nodeforce[0], nodeforce[1],
-                                              nodeforce[2]};
-  }
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeVectorQuantity("Internal Axial force", internalAxialForces);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("Internal Axial force")
-      ->setEnabled(false);
-  // per node internal first bending force
   std::vector<std::array<double, 3>> internalFirstBendingTranslationForces(
       mesh->VN());
   std::vector<std::array<double, 3>> internalFirstBendingRotationForces(
@@ -1456,7 +1573,26 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
   std::vector<double> nRs(mesh->VN());
   std::vector<double> theta2(mesh->VN());
   std::vector<double> theta3(mesh->VN());
+  std::vector<std::array<double, 3>> residualForces(mesh->VN());
+  std::vector<double> residualForcesNorm(mesh->VN());
+  std::vector<double> accelerationX(mesh->VN());
   for (const VertexType &v : mesh->vert) {
+    kineticEnergies[mesh->getIndex(v)] = mesh->nodes[v].kineticEnergy;
+    const VectorType n = v.cN();
+    nodeNormals[mesh->getIndex(v)] = {n[0], n[1], n[2]};
+    // per node internal forces
+    const Vector6d nodeforce = mesh->nodes[v].force.internal * (-1);
+    internalForces[mesh->getIndex(v)] = {nodeforce[0], nodeforce[1],
+                                         nodeforce[2]};
+    internalForcesNorm[mesh->getIndex(v)] = nodeforce.norm();
+    // External force
+    const Vector6d nodeExternalForce = mesh->nodes[v].force.external;
+    externalForces[mesh->getIndex(v)] = {
+        nodeExternalForce[0], nodeExternalForce[1], nodeExternalForce[2]};
+    externalMoments[mesh->getIndex(v)] = {nodeExternalForce[3],
+                                          nodeExternalForce[4], 0};
+    internalAxialForces[mesh->getIndex(v)] = {nodeforce[0], nodeforce[1],
+                                              nodeforce[2]};
     const Node &node = mesh->nodes[v];
     const Vector6d nodeForceFirst = node.force.internalFirstBending * (-1);
     internalFirstBendingTranslationForces[mesh->getIndex(v)] = {
@@ -1470,116 +1606,89 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
     internalSecondBendingRotationForces[mesh->getIndex(v)] = {
         nodeForceSecond[3], nodeForceSecond[4], 0};
     nRs[mesh->getIndex(v)] = node.nR;
-    const std::vector<EdgePointer> incidentEdges = node.incidentElements;
-    const EdgeIndex ei = mesh->getIndex(incidentEdges.back());
-    //    theta2[mesh->getIndex(v)] =
-    //    node.rotationalDisplacements.at(ei).theta2;
-    //    theta3[mesh->getIndex(v)] =
-    //    node.rotationalDisplacements.at(ei).theta3;
-  }
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeVectorQuantity("First bending force-Translation",
-                              internalFirstBendingTranslationForces);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("First bending force-Translation")
-      ->setEnabled(false);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeVectorQuantity("First bending force-Rotation",
-                              internalFirstBendingRotationForces);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("First bending force-Rotation")
-      ->setEnabled(false);
-
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeVectorQuantity("Second bending force-Translation",
-                              internalSecondBendingTranslationForces);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("Second bending force-Translation")
-      ->setEnabled(false);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeVectorQuantity("Second bending force-Rotation",
-                              internalSecondBendingRotationForces);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("Second bending force-Rotation")
-      ->setEnabled(false);
-
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeScalarQuantity("nR", nRs);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("nR")
-      ->setEnabled(false);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeScalarQuantity("theta3", theta3);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("theta3")
-      ->setEnabled(false);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeScalarQuantity("theta2", theta2);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("theta2")
-      ->setEnabled(false);
-
-  // per node residual forces
-  std::vector<std::array<double, 3>> residualForces(mesh->VN());
-  std::vector<double> residualForcesNorm(mesh->VN());
-  for (const VertexType &v : mesh->vert) {
     const Vector6d nodeResidualForce = mesh->nodes[v].force.residual;
     residualForces[mesh->getIndex(v)] = {
         nodeResidualForce[0], nodeResidualForce[1], nodeResidualForce[2]};
     residualForcesNorm[mesh->getIndex(v)] = nodeResidualForce.norm();
-  }
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeVectorQuantity("Residual force", residualForces);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("Residual force")
-      ->setEnabled(false);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addNodeScalarQuantity("Residual force norm", residualForcesNorm);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("Residual force norm")
-      ->setEnabled(false);
-
-  // per node x acceleration
-  std::vector<double> accelerationX(mesh->VN());
-  for (const VertexType &v : mesh->vert) {
     accelerationX[mesh->getIndex(v)] = mesh->nodes[v].acceleration[0];
   }
-  polyscope::getCurveNetwork("Deformed edge mesh")
+  polyscope::getCurveNetwork(meshPolyscopeLabel)
-      ->addNodeScalarQuantity("Node acceleration x", accelerationX);
+      ->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);
 
-  // per element t1
+  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  std::vector<std::array<double, 3>> t1s(mesh->EN());
+      ->addNodeVectorQuantity("Second bending force-Translation",
-  for (const EdgeType &e : mesh->edge) {
+                              internalSecondBendingTranslationForces)
-    const VectorType &t1 = mesh->elements[e].frame.t1;
-    t1s[mesh->getIndex(e)] = {t1[0], t1[1], t1[2]};
-  }
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addEdgeVectorQuantity("t1Check", t1s);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("t1Check")
       ->setEnabled(false);
-  // per element t2
+  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  std::vector<std::array<double, 3>> t2s(mesh->EN());
+      ->addNodeVectorQuantity("Second bending force-Rotation",
-  for (const EdgeType &e : mesh->edge) {
+                              internalSecondBendingRotationForces);
-    const VectorType &t2 = mesh->elements[e].frame.t2;
+  polyscope::getCurveNetwork(meshPolyscopeLabel)
-    t2s[mesh->getIndex(e)] = {t2[0], t2[1], t2[2]};
+      ->getQuantity("Second bending force-Rotation")
-  }
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addEdgeVectorQuantity("t2", t2s);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("t2")
       ->setEnabled(false);
-  // per element t3
+  polyscope::getCurveNetwork(meshPolyscopeLabel)
-  std::vector<std::array<double, 3>> t3s(mesh->EN());
+      ->addNodeScalarQuantity("nR", nRs)
-  for (const EdgeType &e : mesh->edge) {
-    const VectorType &t3 = mesh->elements[e].frame.t3;
-    t3s[mesh->getIndex(e)] = {t3[0], t3[1], t3[2]};
-  }
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->addEdgeVectorQuantity("t3", t3s);
-  polyscope::getCurveNetwork("Deformed edge mesh")
-      ->getQuantity("t3")
       ->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()->addNodeScalarQuantity("Node acceleration x",
+                                                      accelerationX);
+
+  // Edge quantities
+  std::vector<double> A(mesh->EN());
+  std::vector<double> J(mesh->EN());
+  std::vector<double> I2(mesh->EN());
+  std::vector<double> I3(mesh->EN());
+  for (const EdgeType &e : mesh->edge) {
+    const size_t ei = mesh->getIndex(e);
+    A[ei] = mesh->elements[e].properties.A;
+    J[ei] = mesh->elements[e].properties.J;
+    I2[ei] = mesh->elements[e].properties.I2;
+    I3[ei] = mesh->elements[e].properties.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 = [&]() {
@@ -1618,15 +1727,23 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
   }
 }
 
-SimulationResults
+SimulationResults FormFinder::executeSimulation(
-FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose,
+    const SimulationJob &job, const bool &beVerbose, const bool &shouldDraw,
-                              const bool &shouldDraw,
+    const bool &shouldCreatePlots, const size_t &maxDRMIterations) {
-                              const size_t &maxDRMIterations) {
   assert(job.mesh.operator bool());
   auto t1 = std::chrono::high_resolution_clock::now();
   reset();
   mShouldDraw = shouldDraw;
 
+  //  if(!job.nodalExternalForces.empty()){
+  //      double externalForcesNorm=0;
+  //      for(const auto& externalForce:job.nodalExternalForces)
+  //      {
+  //          externalForcesNorm+=externalForce.norm();
+  //      }
+  //      setTotalResidualForcesNormThreshold(externalForcesNorm);
+  //  }
+
   constrainedVertices = job.fixedVertices;
   if (!job.nodalForcedDisplacements.empty()) {
     for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair :
@@ -1649,13 +1766,12 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose,
   }
   computeRigidSupports();
 
-  const std::string deformedMeshName = "Deformed edge mesh";
+  if (mShouldDraw || true) {
-  if (mShouldDraw) {
     if (!polyscope::state::initialized) {
       initPolyscope();
     }
-    polyscope::registerCurveNetwork(deformedMeshName, mesh->getEigenVertices(),
+    polyscope::registerCurveNetwork(
-                                    mesh->getEigenEdges());
+        meshPolyscopeLabel, mesh->getEigenVertices(), mesh->getEigenEdges());
     //    registerWorldAxes();
   }
   for (auto fixedVertex : job.fixedVertices) {
@@ -1670,7 +1786,7 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose,
   }
   updateNodalExternalForces(job.nodalExternalForces, constrainedVertices);
 
-  while (mCurrentSimulationStep < maxDRMIterations) {
+  while (maxDRMIterations == 0 || mCurrentSimulationStep < maxDRMIterations) {
     //  while (true) {
     updateNormalDerivatives();
     updateT1Derivatives();
@@ -1701,24 +1817,40 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose,
     if (mCurrentSimulationStep != 0) {
       history.stepPulse(*mesh);
     }
-    if (mShouldDraw &&
+    if (mCurrentSimulationStep > 100000) {
-        mCurrentSimulationStep % mDrawingStep == 0 /* &&
+      shouldUseKineticEnergyThreshold = true;
+    }
+    if (mCurrentSimulationStep>500000 ||(mShouldDraw  &&
+                                            mCurrentSimulationStep % mDrawingStep == 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;
-      const auto yValues = history.residualForces;
+      auto yValues = history.residualForces;
       auto xPlot = matplot::linspace(0, yValues.size(), yValues.size());
       plotHandle = matplot::scatter(xPlot, yValues);
-      const std::string label = "Residual forces";
+      std::string label = "Residual forces";
       plotHandle->legend_string(label);
-      draw();
+
+      //      yValues = history.kineticEnergy;
+      //      plotHandle = matplot::scatter(xPlot, yValues);
+      //      label = "Log of Kinetic energy";
+      //      plotHandle->legend_string(label);
+
+      //      shouldUseKineticEnergyThreshold = true;
     }
     //    t = t + Dt;
     mCurrentSimulationStep++;
@@ -1728,8 +1860,7 @@ currentSimulationStep > maxDRMIterations*/) {
     //              << std::endl;
     if (mesh->previousTotalKineticEnergy > mesh->currentTotalKineticEnergy) {
       if (/*mesh.totalPotentialEnergykN < 10 ||*/
-          //          mesh->totalResidualForcesNorm <
+          mesh->totalResidualForcesNorm < totalResidualForcesNormThreshold ||
-          //          totalResidualForcesNormThreshold ||
           (shouldUseKineticEnergyThreshold &&
            mesh->currentTotalKineticEnergy < totalKineticEnergyThreshold)) {
         if (beVerbose) {
@@ -1744,8 +1875,7 @@ currentSimulationStep > maxDRMIterations*/) {
                     << " kNm" << std::endl;
         }
         if (shouldUseKineticEnergyThreshold) {
-          std::cout << "Warning: Since maximum applied external moment reached "
+          std::cout << "Warning: The kinetic energy of the system was "
-                       "the kinetic energy of the system was "
                        " used as a convergence criterion"
                     << std::endl;
         }
@@ -1775,197 +1905,15 @@ currentSimulationStep > maxDRMIterations*/) {
   //    mesh.printVertexCoordinates(mesh.VN() / 2);
   if (mShouldDraw) {
     draw();
-    //  if (!polyscope::hasCurveNetwork(deformedMeshName)) {
+  }
-    polyscope::removeCurveNetwork(deformedMeshName);
+  if (polyscope::hasCurveNetwork(meshPolyscopeLabel)) {
-    //  }
+    polyscope::removeCurveNetwork(meshPolyscopeLabel);
   }
   SimulationResults results = computeResults(*job.mesh);
+  results.label = job.mesh->getLabel();
+  if (shouldCreatePlots) {
+    SimulationResultsReporter reporter;
+    reporter.reportResults({results}, "Results", job.mesh->getLabel());
+  }
   return results;
 }
-
-void FormFinder::runUnitTests() {
-  const std::filesystem::path groundOfTruthFolder{
-      "/home/iason/Coding/Libraries/MySources/formFinder_unitTestFiles"};
-
-  FormFinder formFinder;
-  formFinder.setTotalResidualForcesNormThreshold(1);
-
-  // First example of the paper
-  VCGEdgeMesh beam;
-  //  const size_t spanGridSize = 11;
-  //  mesh.createSpanGrid(spanGridSize);
-  beam.loadFromPly(
-      std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply"));
-  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};
-  std::unordered_map<VertexIndex, Vector6d> nodalForces{
-      {2, Vector6d({0, 1000, 1000, 0, 0, 0})}};
-  // Forced displacements
-  std::unordered_map<size_t, Eigen::Vector3d> nodalForcedDisplacements;
-  nodalForcedDisplacements.insert({beam.VN() - 1, Eigen::Vector3d(-0.2, 0, 0)});
-
-  SimulationJob beamSimulationJob{
-      std::make_shared<SimulationMesh>(beam),
-      //      SimulationJob::constructFixedVerticesSpanGrid(spanGridSize,
-      //      mesh.VN()),
-      fixedVertices, nodalForces, nodalForcedDisplacements};
-  beamSimulationJob.mesh->setBeamMaterial(0.3, 200);
-  assert(CrossSectionType::name == CylindricalBeamDimensions::name);
-
-  beamSimulationJob.mesh->setBeamCrossSection(CrossSectionType{0.03, 0.026});
-  SimulationResults simpleBeam_simulationResults =
-      formFinder.executeSimulation(beamSimulationJob, false, true);
-  simpleBeam_simulationResults.label = "SimpleBeam";
-  simpleBeam_simulationResults.save();
-  const std::string simpleBeamGroundOfTruthBinaryFilename =
-      std::filesystem::path(groundOfTruthFolder)
-          .append("SimpleBeam_displacements.eigenBin");
-  assert(std::filesystem::exists(
-      std::filesystem::path(simpleBeamGroundOfTruthBinaryFilename)));
-  Eigen::MatrixXd simpleBeam_groundOfTruthDisplacements;
-  Eigen::read_binary(simpleBeamGroundOfTruthBinaryFilename,
-                     simpleBeam_groundOfTruthDisplacements);
-  if (!simpleBeam_simulationResults.isEqual(
-          simpleBeam_groundOfTruthDisplacements)) {
-    std::cerr << "Error:Beam simulation produces wrong results." << std::endl;
-    return;
-  }
-
-  // Second example of the paper
-  VCGEdgeMesh shortSpanGrid;
-  //  const size_t spanGridSize = 11;
-  //  mesh.createSpanGrid(spanGridSize);
-  shortSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder)
-                                .append("shortSpanGridshell.ply")
-                                .string());
-
-  fixedVertices.clear();
-  //// Corner nodes
-  fixedVertices[0] = std::unordered_set<DoFType>{2};
-  fixedVertices[4] = std::unordered_set<DoFType>{2};
-  fixedVertices[16] = std::unordered_set<DoFType>{2};
-  fixedVertices[20] = std::unordered_set<DoFType>{2};
-  //// Center node
-  fixedVertices[10] = std::unordered_set<DoFType>{0, 1, 3, 4, 5};
-
-  nodalForcedDisplacements.clear();
-  nodalForcedDisplacements.insert({{0, Eigen::Vector3d(0.1, 0.1, 0)},
-                                   {4, Eigen::Vector3d(-0.1, 0.1, 0)},
-                                   {16, Eigen::Vector3d(0.1, -0.1, 0)},
-                                   {20, Eigen::Vector3d(-0.1, -0.1, 0)}});
-
-  SimulationJob shortSpanGridshellSimulationJob{
-      std::make_shared<SimulationMesh>(shortSpanGrid),
-      fixedVertices,
-      {},
-      nodalForcedDisplacements};
-  shortSpanGridshellSimulationJob.mesh->setBeamMaterial(0.3, 200);
-  assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
-  shortSpanGridshellSimulationJob.mesh->setBeamCrossSection(
-      CrossSectionType{0.03, 0.026});
-  SimulationResults shortSpanGridshellSimulationResults =
-      formFinder.executeSimulation(shortSpanGridshellSimulationJob, false,
-                                   false);
-  shortSpanGridshellSimulationResults.label = "ShortSpanGridshell";
-  shortSpanGridshellSimulationResults.save();
-
-  const std::string groundOfTruthBinaryFilename =
-      std::filesystem::path(groundOfTruthFolder)
-          .append("ShortSpanGridshell_displacements.eigenBin")
-          .string();
-  assert(std::filesystem::exists(
-      std::filesystem::path(groundOfTruthBinaryFilename)));
-  Eigen::MatrixXd shortSpanGridshell_groundOfTruthDisplacements;
-  Eigen::read_binary(groundOfTruthBinaryFilename,
-                     shortSpanGridshell_groundOfTruthDisplacements);
-  //  shortSpanGridshellSimulationResults.registerForDrawing(
-  //      shortSpanGridshellSimulationJob);
-  //  polyscope::show();
-  assert(shortSpanGridshellSimulationResults.isEqual(
-      shortSpanGridshell_groundOfTruthDisplacements));
-  if (!shortSpanGridshellSimulationResults.isEqual(
-          shortSpanGridshell_groundOfTruthDisplacements)) {
-    std::cerr << "Error:Short span simulation produces wrong results."
-              << std::endl;
-    return;
-  }
-
-  // Third example of the paper
-  VCGEdgeMesh longSpanGrid;
-  longSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder)
-                               .append("longSpanGridshell.ply")
-                               .string());
-  const size_t spanGridSize = 11;
-
-  fixedVertices.clear();
-  for (size_t vi = 0; vi < spanGridSize - 1; vi++) {
-    fixedVertices[vi] = {0, 2};
-  }
-  for (size_t vi = longSpanGrid.VN() - 1 - (spanGridSize - 2);
-       vi < longSpanGrid.VN(); vi++) {
-    fixedVertices[vi] = {0, 2};
-  }
-  for (size_t vi = spanGridSize - 1;
-       vi < longSpanGrid.VN() - 1 - (spanGridSize - 2) - spanGridSize + 1;
-       vi += spanGridSize + 1) {
-    fixedVertices[vi] = {1, 2};
-    fixedVertices[vi + spanGridSize] = {1, 2};
-  }
-
-  nodalForcedDisplacements.clear();
-  const size_t horizontalOffset = std::floor((spanGridSize - 2) / 2);
-  nodalForcedDisplacements.insert(
-      {{horizontalOffset, Eigen::Vector3d(0, 0.3, 0)},
-       {horizontalOffset + 1, Eigen::Vector3d(0, 0.3, 0)},
-       {spanGridSize * (spanGridSize + 1) - 2 + horizontalOffset,
-        Eigen::Vector3d(0, -0.3, 0)},
-       {spanGridSize * (spanGridSize + 1) - 2 + horizontalOffset + 1,
-        Eigen::Vector3d(0, -0.3, 0)},
-       {std::floor(spanGridSize / 2) * (spanGridSize + 1) - 2,
-        Eigen::Vector3d(0.3, 0, 0)},
-       {(std::floor(spanGridSize / 2) + 1) * (spanGridSize + 1) - 2,
-        Eigen::Vector3d(0.3, 0, 0)},
-       {std::floor(spanGridSize / 2) * (spanGridSize + 1) - 2 + spanGridSize,
-        Eigen::Vector3d(-0.3, 0, 0)},
-       {(std::floor(spanGridSize / 2) + 1) * (spanGridSize + 1) - 2 +
-            spanGridSize,
-        Eigen::Vector3d(-0.3, 0, 0)}});
-
-  SimulationJob longSpanGridshell_simulationJob{
-      std::make_shared<SimulationMesh>(longSpanGrid),
-      fixedVertices,
-      {},
-      nodalForcedDisplacements};
-  longSpanGridshell_simulationJob.mesh->setBeamMaterial(0.3, 200);
-  assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
-  longSpanGridshell_simulationJob.mesh->setBeamCrossSection(
-      CrossSectionType{0.03, 0.026});
-  SimulationResults longSpanGridshell_simulationResults =
-      formFinder.executeSimulation(longSpanGridshell_simulationJob, false,
-                                   false);
-  longSpanGridshell_simulationResults.label = "LongSpanGridshell";
-  longSpanGridshell_simulationResults.save();
-
-  const std::string longSpan_groundOfTruthBinaryFilename =
-      std::filesystem::path(groundOfTruthFolder)
-          .append("LongSpanGridshell_displacements.eigenBin")
-          .string();
-  assert(std::filesystem::exists(
-      std::filesystem::path(longSpan_groundOfTruthBinaryFilename)));
-  Eigen::MatrixXd longSpanGridshell_groundOfTruthDisplacements;
-  Eigen::read_binary(longSpan_groundOfTruthBinaryFilename,
-                     longSpanGridshell_groundOfTruthDisplacements);
-  assert(longSpanGridshell_simulationResults.isEqual(
-      longSpanGridshell_groundOfTruthDisplacements));
-  if (!longSpanGridshell_simulationResults.isEqual(
-          longSpanGridshell_groundOfTruthDisplacements)) {
-    std::cerr << "Error:Long span simulation produces wrong results."
-              << std::endl;
-    return;
-  }
-
-  std::cout << "Form finder unit tests succesufully passed." << std::endl;
-
-  //  polyscope::show();
-}

beamformfinder.hpp

@@ -27,16 +27,18 @@ private:
   const double Dtini{0.1};
   double Dt{Dtini};
   const double xi{0.9969};
-  double totalResidualForcesNormThreshold{0.001};
+  double totalResidualForcesNormThreshold{1e-5};
-  bool shouldUseKineticEnergyThreshold{true};
+  bool shouldUseKineticEnergyThreshold{false};
-  const double totalKineticEnergyThreshold{1e-11};
+  bool checkedForMaximumMoment{false};
+  const double totalKineticEnergyThreshold{1e-9};
   double externalMomentsNorm{0};
   size_t mCurrentSimulationStep{0};
-  int mDrawingStep{5};
+  int mDrawingStep{1};
   bool mShouldDraw{false};
   matplot::line_handle plotHandle;
   std::vector<double> plotYValues;
 
+  const std::string meshPolyscopeLabel{"Simulation mesh"};
   std::unique_ptr<SimulationMesh> mesh;
   std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
       constrainedVertices;
@@ -186,9 +188,9 @@ public:
   FormFinder();
   SimulationResults executeSimulation(
       const SimulationJob &job, const bool &beVerbose = false,
-      const bool &shouldDraw = false,
+      const bool &shouldDraw = false, const bool &createPlots = false,
       const size_t &maxDRMIterations = FormFinder::maxDRMIterations);
-  inline static const size_t maxDRMIterations{500000};
+  inline static const size_t maxDRMIterations{0};
 
   static void runUnitTests();
   void setTotalResidualForcesNormThreshold(double value);

edgemesh.hpp

@@ -21,7 +21,7 @@ struct VCGEdgeMeshUsedTypes
 class VCGEdgeMeshVertexType
     : public vcg::Vertex<VCGEdgeMeshUsedTypes, vcg::vertex::Coord3d,
                          vcg::vertex::Normal3d, vcg::vertex::BitFlags,
-                         vcg::vertex::VEAdj> {};
+                         vcg::vertex::Color4b, vcg::vertex::VEAdj> {};
 class VCGEdgeMeshEdgeType
     : public vcg::Edge<VCGEdgeMeshUsedTypes, vcg::edge::VertexRef,
                        vcg::edge::BitFlags, vcg::edge::EEAdj,

flatpattern.cpp

@@ -277,6 +277,14 @@ void FlatPattern::createFan(const size_t &fanSize) {
   }
 
   removeDuplicateVertices();
+  //  const double precision = 1e-4;
+  //    for (size_t vi = 0; vi < VN(); vi++) {
+  //      vert[vi].P()[0] = std::round(vert[vi].P()[0] * (1 / precision)) *
+  //      precision; vert[vi].P()[1] = std::round(vert[vi].P()[1] * (1 /
+  //      precision)) * precision; vert[vi].P()[2] = std::round(vert[vi].P()[2]
+  //      * (1 / precision)) * precision;
+  //    }
+
   updateEigenEdgeAndVertices();
 }
 

trianglepatterngeometry.cpp

@@ -80,7 +80,13 @@ double FlatPatternGeometry::getTriangleEdgeSize() const {
 
 FlatPatternGeometry::FlatPatternGeometry() {}
 
-std::vector<vcg::Point3d> FlatPatternGeometry::getVertices() const {}
+std::vector<vcg::Point3d> FlatPatternGeometry::getVertices() const {
+  std::vector<VCGEdgeMesh::CoordType> verts(VN());
+  for (size_t vi = 0; vi < VN(); vi++) {
+    verts[vi] = vert[vi].cP();
+  }
+  return verts;
+}
 
 FlatPatternGeometry
 FlatPatternGeometry::createTile(FlatPatternGeometry &pattern) {
@@ -235,7 +241,9 @@ std::vector<vcg::Point3d> FlatPatternGeometry::constructVertexVector(
     }
   }
 
-  if (numberOfNodesPerSlot[4] != 0) {
+  if (numberOfNodesPerSlot[4] == 1) {
+    vertices.push_back(vcg::Point3d(0, -triangleHeight, 0));
+  } else if (numberOfNodesPerSlot[4] != 0) {
     const double offset1 = 1.0 / (numberOfNodesPerSlot[4] + 1);
     const vcg::Point3d edgeVector1(triangleV2 - triangleV1);
     for (size_t vertexIndex = 0; vertexIndex < numberOfNodesPerSlot[4];