diff --git a/beam.hpp b/beam.hpp
index 4851276..e965c9c 100644
--- a/beam.hpp
+++ b/beam.hpp
@@ -6,6 +6,7 @@
 #include <iostream>
 
 struct RectangularBeamDimensions {
+  inline static std::string name{"Rectangular"};
   float b;
   float h;
   RectangularBeamDimensions(const float &width, const float &height)
@@ -16,6 +17,7 @@ struct RectangularBeamDimensions {
 };
 
 struct CylindricalBeamDimensions {
+  inline static std::string name{"Cylindrical"};
   float od; // Cylinder outside diameter
   float
       id; // Cylinder inside diameter
diff --git a/beamformfinder.cpp b/beamformfinder.cpp
index 1d39e20..853855e 100644
--- a/beamformfinder.cpp
+++ b/beamformfinder.cpp
@@ -8,6 +8,10 @@
 #include <execution>
 #include <omp.h>
 
+void FormFinder::setTotalResidualForcesNormThreshold(double value) {
+  totalResidualForcesNormThreshold = value;
+}
+
 void FormFinder::reset() {
   Dt = Dtini;
   mCurrentSimulationStep = 0;
@@ -17,6 +21,8 @@ void FormFinder::reset() {
   mesh.release();
   plotYValues.clear();
   plotHandle.reset();
+  shouldUseKineticEnergyThreshold = false;
+  externalMomentsNorm = 0;
 }
 
 VectorType FormFinder::computeDisplacementDifferenceDerivative(
@@ -779,13 +785,17 @@ void FormFinder::updateResidualForcesOnTheFly(
       force.residual = force.residual + (internalForcePair.second * -1);
     }
   }
+  Vector6d sumOfResidualForces(0);
   for (size_t vi = 0; vi < mesh->VN(); vi++) {
-    const Vector6d &nodeResidualForce = mesh->nodes[vi].force.residual;
+    Node::Forces &force = mesh->nodes[vi].force;
+    const Vector6d &nodeResidualForce = force.residual;
+    //    sumOfResidualForces = sumOfResidualForces + nodeResidualForce;
     const double residualForceNorm = nodeResidualForce.norm();
     totalResidualForcesNorm += residualForceNorm;
   }
   mesh->previousTotalResidualForcesNorm = mesh->totalResidualForcesNorm;
   mesh->totalResidualForcesNorm = totalResidualForcesNorm;
+  //  mesh->totalResidualForcesNorm = sumOfResidualForces.norm();
 
   //  plotYValues.push_back(totalResidualForcesNorm);
   //  auto xPlot = matplot::linspace(0, plotYValues.size(), plotYValues.size());
@@ -797,6 +807,7 @@ void FormFinder::updateNodalExternalForces(
     const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
         &fixedVertices) {
 
+  externalMomentsNorm = 0;
   for (const std::pair<VertexIndex, Vector6d> &nodalForce : nodalForces) {
     const VertexIndex nodeIndex = nodalForce.first;
     const bool isNodeConstrained = fixedVertices.contains(nodeIndex);
@@ -810,6 +821,12 @@ void FormFinder::updateNodalExternalForces(
       }
       nodalExternalForce[dofi] = nodalForce.second[dofi];
     }
+    externalMomentsNorm += std::sqrt(pow(nodalExternalForce[3], 2) +
+                                     pow(nodalExternalForce[4], 2));
+    //    CoordType v = (mesh->vert[nodeIndex].cP() -
+    //    mesh->vert[0].cP()).Normalize(); const double forceMagnitude = 0.1;
+    //    nodalExternalForce[3] = v[0] * forceMagnitude;
+    //    nodalExternalForce[4] = v[1] * forceMagnitude;
     node.force.external = nodalExternalForce;
   }
 }
@@ -1032,7 +1049,7 @@ void FormFinder::updateNodePosition(
   node.previousLocation = previousLocation;
   v.P() = node.initialLocation + displacementVector;
   if (shouldApplyInitialDistortion && mCurrentSimulationStep < 40) {
-    VectorType desiredInitialDisplacement(0, 0, 0.001);
+    VectorType desiredInitialDisplacement(0, 0, 0.01);
     v.P() = v.P() + desiredInitialDisplacement;
   }
 }
@@ -1061,11 +1078,34 @@ void FormFinder::updateNodeNormal(
   const double &nx = v.N()[0], ny = v.N()[1], nz = v.N()[2];
   const double nxnyMagnitude = std::pow(nx, 2) + std::pow(ny, 2);
   VectorType n = v.N();
-  const bool shouldBreak = mCurrentSimulationStep == 118 && vi == 1;
+  const bool shouldBreak = mCurrentSimulationStep == 118 && vi == 3;
   if (nxnyMagnitude > 1) {
     VectorType newNormal(nx / std::sqrt(nxnyMagnitude),
                          ny / std::sqrt(nxnyMagnitude), 0);
     v.N() = newNormal;
+
+    /*If an external moment caused the normal to lay on the xy plane this means
+     * that in order to disable its effect a greater internal force is needed
+     * than what is possible (the constraint on the z of the normals imposes a
+     * constraint on the maximum internal force). Because of that the
+     * totalResidualForcesNorm can't drop below the magnitude of external moment
+     * applied on vertex vi. In order to allow termination of the simulation
+     * when the described phenomenon happens we allow the termination of the
+     * 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) {
+      shouldUseKineticEnergyThreshold = true;
+      std::cout
+          << "Maximum moment reached.The Kinetic energy of the system will "
+             "be used as a convergence criterion"
+          << std::endl;
+      wasExecuted = true;
+    }
+
   } else {
     const double nzSquared = 1.0 - nxnyMagnitude;
     const double nz = std::sqrt(nzSquared);
@@ -1358,6 +1398,7 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
 
   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
@@ -1369,6 +1410,8 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
     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);
@@ -1380,6 +1423,10 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
   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")
@@ -1651,6 +1698,9 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose,
     }
     updateElementalLengths();
 
+    if (mCurrentSimulationStep != 0) {
+      history.stepPulse(*mesh);
+    }
     if (mShouldDraw &&
         mCurrentSimulationStep % mDrawingStep == 0 /* &&
 currentSimulationStep > maxDRMIterations*/) {
@@ -1661,11 +1711,15 @@ currentSimulationStep > maxDRMIterations*/) {
       //      draw(saveTo);
       std::cout << "Residual forces norm: " << mesh->totalResidualForcesNorm
                 << std::endl;
+      std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy
+                << std::endl;
+      const auto yValues = history.residualForces;
+      auto xPlot = matplot::linspace(0, yValues.size(), yValues.size());
+      plotHandle = matplot::scatter(xPlot, yValues);
+      const std::string label = "Residual forces";
+      plotHandle->legend_string(label);
       draw();
     }
-    if (mCurrentSimulationStep != 0) {
-      history.stepPulse(*mesh);
-    }
     //    t = t + Dt;
     mCurrentSimulationStep++;
     //    std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy
@@ -1674,7 +1728,10 @@ currentSimulationStep > maxDRMIterations*/) {
     //              << std::endl;
     if (mesh->previousTotalKineticEnergy > mesh->currentTotalKineticEnergy) {
       if (/*mesh.totalPotentialEnergykN < 10 ||*/
-          mesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) {
+          //          mesh->totalResidualForcesNorm <
+          //          totalResidualForcesNormThreshold ||
+          (shouldUseKineticEnergyThreshold &&
+           mesh->currentTotalKineticEnergy < totalKineticEnergyThreshold)) {
         if (beVerbose) {
           std::cout << "Convergence after " << mCurrentSimulationStep
                     << " steps" << std::endl;
@@ -1686,6 +1743,12 @@ currentSimulationStep > maxDRMIterations*/) {
           std::cout << "Total potential:" << mesh->totalPotentialEnergykN
                     << " kNm" << std::endl;
         }
+        if (shouldUseKineticEnergyThreshold) {
+          std::cout << "Warning: Since maximum applied external moment reached "
+                       "the kinetic energy of the system was "
+                       " used as a convergence criterion"
+                    << std::endl;
+        }
         break;
         //        }
       }
@@ -1721,37 +1784,188 @@ currentSimulationStep > maxDRMIterations*/) {
 }
 
 void FormFinder::runUnitTests() {
+  const std::filesystem::path groundOfTruthFolder{
+      "/home/iason/Coding/Libraries/MySources/formFinder_unitTestFiles"};
+
   FormFinder formFinder;
-  VCGEdgeMesh mesh;
+  formFinder.setTotalResidualForcesNormThreshold(1);
+
+  // First example of the paper
+  VCGEdgeMesh beam;
   //  const size_t spanGridSize = 11;
   //  mesh.createSpanGrid(spanGridSize);
-  mesh.loadFromPly("/home/iason/Models/simple_beam_paper_example.ply");
+  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[mesh.VN() - 1] = std::unordered_set<DoFType>{1, 2};
+  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({mesh.VN() - 1, Eigen::Vector3d(-0.2, 0, 0)});
+  nodalForcedDisplacements.insert({beam.VN() - 1, Eigen::Vector3d(-0.2, 0, 0)});
 
   SimulationJob beamSimulationJob{
-      std::make_shared<SimulationMesh>(mesh),
+      std::make_shared<SimulationMesh>(beam),
       //      SimulationJob::constructFixedVerticesSpanGrid(spanGridSize,
       //      mesh.VN()),
       fixedVertices, nodalForces, nodalForcedDisplacements};
   beamSimulationJob.mesh->setBeamMaterial(0.3, 200);
-  assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
-  beamSimulationJob.mesh->setBeamCrossSection(CrossSectionType{0.03, 0.026});
-  SimulationResults beamSimulationResults =
-      formFinder.executeSimulation(beamSimulationJob, true, false);
+  assert(CrossSectionType::name == CylindricalBeamDimensions::name);
 
-  const bool testSuccessful =
-      beamSimulationResults.displacements[2][0] - (-0.09556) < 1e-5 &&
-      beamSimulationResults.displacements[2][1] - 0.40666 < 1e-5 &&
-      beamSimulationResults.displacements[2][2] - 0.40666 < 1e-5;
-  assert(testSuccessful);
-  beamSimulationResults.simulationLabel = "Beam";
-  beamSimulationResults.registerForDrawing(beamSimulationJob);
-  polyscope::show();
+  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();
 }
diff --git a/beamformfinder.hpp b/beamformfinder.hpp
index 2e6c08f..e7d31c5 100644
--- a/beamformfinder.hpp
+++ b/beamformfinder.hpp
@@ -27,7 +27,10 @@ private:
   const double Dtini{0.1};
   double Dt{Dtini};
   const double xi{0.9969};
-  const double totalResidualForcesNormThreshold{0.01};
+  double totalResidualForcesNormThreshold{0.001};
+  bool shouldUseKineticEnergyThreshold{true};
+  const double totalKineticEnergyThreshold{1e-11};
+  double externalMomentsNorm{0};
   size_t mCurrentSimulationStep{0};
   int mDrawingStep{5};
   bool mShouldDraw{false};
@@ -188,6 +191,7 @@ public:
   inline static const size_t maxDRMIterations{500000};
 
   static void runUnitTests();
+  void setTotalResidualForcesNormThreshold(double value);
 };
 
 template <typename PointType> PointType Cross(PointType p1, PointType p2) {
@@ -211,28 +215,4 @@ inline bool TriggerBreakpoint(const VertexIndex &vi, const EdgeIndex &ei) {
   return (vi == 1 && ei == 0) || (vi == 9 && ei == 9);
 }
 
-namespace Eigen {
-template <class Matrix>
-void write_binary(const std::string &filename, const Matrix &matrix) {
-  std::ofstream out(filename,
-                    std::ios::out | std::ios::binary | std::ios::trunc);
-  typename Matrix::Index rows = matrix.rows(), cols = matrix.cols();
-  out.write((char *)(&rows), sizeof(typename Matrix::Index));
-  out.write((char *)(&cols), sizeof(typename Matrix::Index));
-  out.write((char *)matrix.data(),
-            rows * cols * sizeof(typename Matrix::Scalar));
-  out.close();
-}
-template <class Matrix>
-void read_binary(const std::string &filename, Matrix &matrix) {
-  std::ifstream in(filename, std::ios::in | std::ios::binary);
-  typename Matrix::Index rows = 0, cols = 0;
-  in.read((char *)(&rows), sizeof(typename Matrix::Index));
-  in.read((char *)(&cols), sizeof(typename Matrix::Index));
-  matrix.resize(rows, cols);
-  in.read((char *)matrix.data(), rows * cols * sizeof(typename Matrix::Scalar));
-  in.close();
-}
-} // namespace Eigen
-
 #endif // BEAMFORMFINDER_HPP
diff --git a/edgemesh.cpp b/edgemesh.cpp
index 4ca26fd..fa1fa44 100644
--- a/edgemesh.cpp
+++ b/edgemesh.cpp
@@ -12,7 +12,17 @@ Eigen::MatrixX3d VCGEdgeMesh::getEigenEdgeNormals() const {
   return eigenEdgeNormals;
 }
 
-bool VCGEdgeMesh::savePly(const std::string plyFilename) {}
+bool VCGEdgeMesh::savePly(const std::string plyFilename) {
+  unsigned int mask = 0;
+  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
+  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
+  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
+  if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::SaveModel(
+          plyFilename.c_str(), *this, mask, false) != 0) {
+    return false;
+  }
+  return true;
+}
 
 void VCGEdgeMesh::GeneratedRegularSquaredPattern(
     const double angleDeg, std::vector<std::vector<vcg::Point2d>> &pattern,
@@ -90,7 +100,7 @@ bool VCGEdgeMesh::createSpanGrid(const size_t squareGridDimension) {
 }
 
 bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth,
-                             const size_t desiredHeight) {
+                                 const size_t desiredHeight) {
   std::cout << "Grid of dimensions:" << desiredWidth << "," << desiredHeight
             << std::endl;
   const VCGEdgeMesh::CoordType n(0, 0, 1);
diff --git a/elementalmesh.hpp b/elementalmesh.hpp
index 463aa98..e60ecdf 100644
--- a/elementalmesh.hpp
+++ b/elementalmesh.hpp
@@ -7,8 +7,8 @@
 
 struct Element;
 struct Node;
-using CrossSectionType = RectangularBeamDimensions;
-// using CrossSectionType = CylindricalBeamDimensions;
+// using CrossSectionType = RectangularBeamDimensions;
+using CrossSectionType = CylindricalBeamDimensions;
 
 class SimulationMesh : public VCGEdgeMesh {
 private:
diff --git a/formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin b/formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin
new file mode 100644
index 0000000..6649189
Binary files /dev/null and b/formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin differ
diff --git a/formFinder_unitTestFiles/ShortSpanGridshell_displacements.eigenBin b/formFinder_unitTestFiles/ShortSpanGridshell_displacements.eigenBin
new file mode 100644
index 0000000..19a07ed
Binary files /dev/null and b/formFinder_unitTestFiles/ShortSpanGridshell_displacements.eigenBin differ
diff --git a/formFinder_unitTestFiles/SimpleBeam_displacements.eigenBin b/formFinder_unitTestFiles/SimpleBeam_displacements.eigenBin
new file mode 100644
index 0000000..61a9a23
Binary files /dev/null and b/formFinder_unitTestFiles/SimpleBeam_displacements.eigenBin differ
diff --git a/formFinder_unitTestFiles/longSpanGridshell.ply b/formFinder_unitTestFiles/longSpanGridshell.ply
new file mode 100644
index 0000000..d56cbb7
--- /dev/null
+++ b/formFinder_unitTestFiles/longSpanGridshell.ply
@@ -0,0 +1,374 @@
+ply
+format ascii 1.0
+comment nanoply generated
+element vertex 140
+property double x
+property double y
+property double z
+property double nx
+property double ny
+property double nz
+element edge 220
+property int vertex1
+property int vertex2
+end_header
+1.000000 0.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 0.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 0.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 0.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 0.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 0.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 0.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 0.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 0.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 0.000000 0.000000 0.000000 0.000000 1.000000
+0.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+1.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+11.000000 1.000000 0.000000 0.000000 0.000000 1.000000
+0.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+1.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+11.000000 2.000000 0.000000 0.000000 0.000000 1.000000
+0.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+1.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+11.000000 3.000000 0.000000 0.000000 0.000000 1.000000
+0.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+1.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+11.000000 4.000000 0.000000 0.000000 0.000000 1.000000
+0.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+1.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+11.000000 5.000000 0.000000 0.000000 0.000000 1.000000
+0.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+1.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+11.000000 6.000000 0.000000 0.000000 0.000000 1.000000
+0.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+1.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+11.000000 7.000000 0.000000 0.000000 0.000000 1.000000
+0.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+1.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+11.000000 8.000000 0.000000 0.000000 0.000000 1.000000
+0.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+1.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+11.000000 9.000000 0.000000 0.000000 0.000000 1.000000
+0.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+1.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+11.000000 10.000000 0.000000 0.000000 0.000000 1.000000
+1.000000 11.000000 0.000000 0.000000 0.000000 1.000000
+2.000000 11.000000 0.000000 0.000000 0.000000 1.000000
+3.000000 11.000000 0.000000 0.000000 0.000000 1.000000
+4.000000 11.000000 0.000000 0.000000 0.000000 1.000000
+5.000000 11.000000 0.000000 0.000000 0.000000 1.000000
+6.000000 11.000000 0.000000 0.000000 0.000000 1.000000
+7.000000 11.000000 0.000000 0.000000 0.000000 1.000000
+8.000000 11.000000 0.000000 0.000000 0.000000 1.000000
+9.000000 11.000000 0.000000 0.000000 0.000000 1.000000
+10.000000 11.000000 0.000000 0.000000 0.000000 1.000000
+0 11
+1 12
+2 13
+3 14
+4 15
+5 16
+6 17
+7 18
+8 19
+9 20
+10 11
+11 23
+11 12
+12 24
+12 13
+13 25
+13 14
+14 26
+14 15
+15 27
+15 16
+16 28
+16 17
+17 29
+17 18
+18 30
+18 19
+19 31
+19 20
+20 32
+20 21
+22 23
+23 35
+23 24
+24 36
+24 25
+25 37
+25 26
+26 38
+26 27
+27 39
+27 28
+28 40
+28 29
+29 41
+29 30
+30 42
+30 31
+31 43
+31 32
+32 44
+32 33
+34 35
+35 47
+35 36
+36 48
+36 37
+37 49
+37 38
+38 50
+38 39
+39 51
+39 40
+40 52
+40 41
+41 53
+41 42
+42 54
+42 43
+43 55
+43 44
+44 56
+44 45
+46 47
+47 59
+47 48
+48 60
+48 49
+49 61
+49 50
+50 62
+50 51
+51 63
+51 52
+52 64
+52 53
+53 65
+53 54
+54 66
+54 55
+55 67
+55 56
+56 68
+56 57
+58 59
+59 71
+59 60
+60 72
+60 61
+61 73
+61 62
+62 74
+62 63
+63 75
+63 64
+64 76
+64 65
+65 77
+65 66
+66 78
+66 67
+67 79
+67 68
+68 80
+68 69
+70 71
+71 83
+71 72
+72 84
+72 73
+73 85
+73 74
+74 86
+74 75
+75 87
+75 76
+76 88
+76 77
+77 89
+77 78
+78 90
+78 79
+79 91
+79 80
+80 92
+80 81
+82 83
+83 95
+83 84
+84 96
+84 85
+85 97
+85 86
+86 98
+86 87
+87 99
+87 88
+88 100
+88 89
+89 101
+89 90
+90 102
+90 91
+91 103
+91 92
+92 104
+92 93
+94 95
+95 107
+95 96
+96 108
+96 97
+97 109
+97 98
+98 110
+98 99
+99 111
+99 100
+100 112
+100 101
+101 113
+101 102
+102 114
+102 103
+103 115
+103 104
+104 116
+104 105
+106 107
+107 119
+107 108
+108 120
+108 109
+109 121
+109 110
+110 122
+110 111
+111 123
+111 112
+112 124
+112 113
+113 125
+113 114
+114 126
+114 115
+115 127
+115 116
+116 128
+116 117
+118 119
+119 130
+119 120
+120 131
+120 121
+121 132
+121 122
+122 133
+122 123
+123 134
+123 124
+124 135
+124 125
+125 136
+125 126
+126 137
+126 127
+127 138
+127 128
+128 139
+128 129
diff --git a/formFinder_unitTestFiles/shortSpanGridshell.ply b/formFinder_unitTestFiles/shortSpanGridshell.ply
new file mode 100644
index 0000000..e25e0cd
--- /dev/null
+++ b/formFinder_unitTestFiles/shortSpanGridshell.ply
@@ -0,0 +1,69 @@
+ply
+format ascii 1.0
+element vertex 21           { define "vertex" element, 8 of them in file }
+property float x           { vertex contains float "x" coordinate }
+property float y           { y coordinate is also a vertex property }
+property float z           { z coordinate, too }
+property float nx
+property float ny
+property float nz
+element edge 24            { there are 6 "face" elements in the file }
+property int vertex1
+property int vertex2
+property list uchar float beam_dimensions
+property list uchar float beam_material {poissons ratio , young's modulus in GPascal}
+end_header
+
+0 0 0 0 0 1
+1 0 0 0 0 1
+2 0 0 0 0 1
+3 0 0 0 0 1
+4 0 0 0 0 1
+0 1 0 0 0 1
+2 1 0 0 0 1
+4 1 0 0 0 1
+0 2 0 0 0 1
+1 2 0 0 0 1
+2 2 0 0 0 1
+3 2 0 0 0 1
+4 2 0 0 0 1
+0 3 0 0 0 1
+2 3 0 0 0 1
+4 3 0 0 0 1
+0 4 0 0 0 1
+1 4 0 0 0 1
+2 4 0 0 0 1
+3 4 0 0 0 1
+4 4 0 0 0 1
+
+0 1 2 0.03 0.026 2 0.3 200
+1 2 2 0.03 0.026 2 0.3 200
+2 3 2 0.03 0.026 2 0.3 200
+3 4 2 0.03 0.026 2 0.3 200
+
+0 5 2 0.03 0.026 2 0.3 200
+2 6 2 0.03 0.026 2 0.3 200
+4 7 2 0.03 0.026 2 0.3 200
+
+
+8 9 2 0.03 0.026 2 0.3 200
+9 10 2 0.03 0.026 2 0.3 200
+10 11 2 0.03 0.026 2 0.3 200
+11 12 2 0.03 0.026 2 0.3 200
+
+8 5 2 0.03 0.026 2 0.3 200
+10 6 2 0.03 0.026 2 0.3 200
+12 7 2 0.03 0.026 2 0.3 200
+
+8 13 2 0.03 0.026 2 0.3 200
+10 14 2 0.03 0.026 2 0.3 200
+12 15 2 0.03 0.026 2 0.3 200
+
+16 13 2 0.03 0.026 2 0.3 200
+18 14 2 0.03 0.026 2 0.3 200
+20 15 2 0.03 0.026 2 0.3 200
+
+16 17 2 0.03 0.026 2 0.3 200
+17 18 2 0.03 0.026 2 0.3 200
+18 19 2 0.03 0.026 2 0.3 200
+19 20 2 0.03 0.026 2 0.3 200
\ No newline at end of file
diff --git a/formFinder_unitTestFiles/simpleBeam.ply b/formFinder_unitTestFiles/simpleBeam.ply
new file mode 100644
index 0000000..1d3d728
--- /dev/null
+++ b/formFinder_unitTestFiles/simpleBeam.ply
@@ -0,0 +1,27 @@
+ply
+format ascii 1.0
+element vertex 5           { define "vertex" element, 8 of them in file }
+property float x           { vertex contains float "x" coordinate }
+property float y           { y coordinate is also a vertex property }
+property float z           { z coordinate, too }
+property float nx
+property float ny
+property float nz
+
+element edge 4            { there are 6 "face" elements in the file }
+property int vertex1
+property int vertex2
+property list uchar float beam_dimensions
+property list uchar float beam_material {poissons ratio , young's modulus in GPascal}
+end_header
+
+0 0 0 0 0 1
+1 0 0 0 0 1
+2 0 0 0 0 1
+3 0 0 0 0 1
+4 0 0 0 0 1
+
+0 1 2 0.03 0.026 2 0.3 200
+1 2 2 0.03 0.026 2 0.3 200
+2 3 2 0.03 0.026 2 0.3 200
+3 4 2 0.03 0.026 2 0.3 200
\ No newline at end of file
diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp
index 6d5ffc5..8335572 100644
--- a/simulationhistoryplotter.hpp
+++ b/simulationhistoryplotter.hpp
@@ -73,7 +73,7 @@ struct SimulationResultsReporter {
     for (const SimulationResults &simulationResult : results) {
       const auto simulationResultPath =
           std::filesystem::path(reportFolderPath)
-              .append(simulationResult.simulationLabel);
+              .append(simulationResult.label);
       std::filesystem::create_directory(simulationResultPath.string());
 
       createPlots(simulationResult.history, simulationResultPath.string(),
diff --git a/simulationresult.hpp b/simulationresult.hpp
index da4492a..c992519 100644
--- a/simulationresult.hpp
+++ b/simulationresult.hpp
@@ -98,33 +98,37 @@ struct SimulationJob {
           ->setEnabled(true);
     }
   }
-  static std::unordered_map<size_t, std::unordered_set<int>>
-  constructFixedVerticesSpanGrid(const size_t &spanGridSize,
-                                 const size_t &gridVertices) {
-    std::unordered_map<size_t, std::unordered_set<int>> fixedVertices;
-    for (size_t vi = 0; vi < spanGridSize - 1; vi++) {
-      fixedVertices[vi] = {0, 1, 2};
-    }
-    for (size_t vi = gridVertices - 1 - (spanGridSize - 2); vi < gridVertices;
-         vi++) {
-      fixedVertices[vi] = {0, 1, 2};
-    }
-    for (size_t vi = spanGridSize - 1;
-         vi < gridVertices - 1 - (spanGridSize - 2) - spanGridSize + 1;
-         vi += spanGridSize + 1) {
-      fixedVertices[vi] = {0, 1, 2};
-      fixedVertices[vi + spanGridSize] = {0, 1, 2};
-    }
-
-    return fixedVertices;
-  }
 };
 
+namespace Eigen {
+template <class Matrix>
+void write_binary(const std::string &filename, const Matrix &matrix) {
+  std::ofstream out(filename,
+                    std::ios::out | std::ios::binary | std::ios::trunc);
+  typename Matrix::Index rows = matrix.rows(), cols = matrix.cols();
+  out.write((char *)(&rows), sizeof(typename Matrix::Index));
+  out.write((char *)(&cols), sizeof(typename Matrix::Index));
+  out.write((char *)matrix.data(),
+            rows * cols * sizeof(typename Matrix::Scalar));
+  out.close();
+}
+template <class Matrix>
+void read_binary(const std::string &filename, Matrix &matrix) {
+  std::ifstream in(filename, std::ios::in | std::ios::binary);
+  typename Matrix::Index rows = 0, cols = 0;
+  in.read((char *)(&rows), sizeof(typename Matrix::Index));
+  in.read((char *)(&cols), sizeof(typename Matrix::Index));
+  matrix.resize(rows, cols);
+  in.read((char *)matrix.data(), rows * cols * sizeof(typename Matrix::Scalar));
+  in.close();
+}
+} // namespace Eigen
+
 struct SimulationResults {
   SimulationHistory history;
   std::vector<Vector6d> displacements;
   double executionTime{0};
-  std::string simulationLabel{"NoLabel"};
+  std::string label{"NoLabel"};
 
   void registerForDrawing(const SimulationJob &job) const {
     polyscope::options::groundPlaneEnabled = false;
@@ -136,12 +140,12 @@ struct SimulationResults {
     } /* else {
        polyscope::removeAllStructures();
      }*/
-    const std::string undeformedMeshName = "Undeformed_" + simulationLabel;
+    const std::string undeformedMeshName = "Undeformed_" + label;
     polyscope::registerCurveNetwork(undeformedMeshName,
                                     job.mesh->getEigenVertices(),
                                     job.mesh->getEigenEdges());
 
-    const std::string deformedMeshName = "Deformed_" + simulationLabel;
+    const std::string deformedMeshName = "Deformed_" + label;
     polyscope::registerCurveNetwork(deformedMeshName,
                                     job.mesh->getEigenVertices(),
                                     job.mesh->getEigenEdges());
@@ -214,6 +218,24 @@ struct SimulationResults {
 
     //    polyscope::show();
   }
+
+  void save() {
+    const std::string filename(label + "_displacements.eigenBin");
+
+    Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements);
+    Eigen::write_binary(filename, m);
+  }
+
+  // The comparison of the results happens comparing the 6-dof nodal
+  // displacements
+  bool isEqual(const Eigen::MatrixXd &nodalDisplacements) {
+    assert(nodalDisplacements.cols() == 6);
+    Eigen::MatrixXd eigenDisplacements =
+        Utilities::toEigenMatrix(this->displacements);
+    const double errorNorm = (eigenDisplacements - nodalDisplacements).norm();
+    return errorNorm < 1e-10;
+    //    return eigenDisplacements.isApprox(nodalDisplacements);
+  }
 };
 
 #endif // SIMULATIONHISTORY_HPP
diff --git a/utilities.hpp b/utilities.hpp
index 5571974..b1d640d 100644
--- a/utilities.hpp
+++ b/utilities.hpp
@@ -6,53 +6,6 @@
 #include <fstream>
 #include <regex>
 
-namespace Utilities {
-inline void parseIntegers(const std::string &str, std::vector<size_t> &result) {
-  typedef std::regex_iterator<std::string::const_iterator> re_iterator;
-  typedef re_iterator::value_type re_iterated;
-
-  std::regex re("(\\d+)");
-
-  re_iterator rit(str.begin(), str.end(), re);
-  re_iterator rend;
-
-  std::transform(rit, rend, std::back_inserter(result),
-                 [](const re_iterated &it) { return std::stoi(it[1]); });
-}
-
-// std::string convertToLowercase(const std::string &s) {
-//  std::string lowercase;
-//  std::transform(s.begin(), s.end(), lowercase.begin(),
-//                 [](unsigned char c) { return std::tolower(c); });
-//  return lowercase;
-//}
-// bool hasExtension(const std::string &filename, const std::string &extension)
-// {
-//  const std::filesystem::path path(filename);
-//  if (!path.has_extension()) {
-//    std::cerr << "Error: No file extension found in " << filename <<
-//    std::endl; return false;
-//  }
-
-//  const std::string detectedExtension = path.extension().string();
-
-//  if (convertToLowercase(detectedExtension) != convertToLowercase(extension))
-//  {
-//    std::cerr << "Error: detected extension is " + detectedExtension +
-//                     " and not " + extension
-//              << std::endl;
-//    return false;
-//  }
-//  return true;
-//}
-
-} // namespace Utilities
-
-struct NodalForce {
-  size_t index;
-  size_t dof;
-  double magnitude;
-};
 struct Vector6d : public std::array<double, 6> {
   Vector6d() {
     for (size_t i = 0; i < 6; i++) {
@@ -145,6 +98,61 @@ struct Vector6d : public std::array<double, 6> {
   }
 };
 
+namespace Utilities {
+inline void parseIntegers(const std::string &str, std::vector<size_t> &result) {
+  typedef std::regex_iterator<std::string::const_iterator> re_iterator;
+  typedef re_iterator::value_type re_iterated;
+
+  std::regex re("(\\d+)");
+
+  re_iterator rit(str.begin(), str.end(), re);
+  re_iterator rend;
+
+  std::transform(rit, rend, std::back_inserter(result),
+                 [](const re_iterated &it) { return std::stoi(it[1]); });
+}
+
+inline Eigen::MatrixXd toEigenMatrix(const std::vector<Vector6d> &v) {
+  Eigen::MatrixXd m(v.size(), 6);
+
+  for (size_t vi = 0; vi < v.size(); vi++) {
+    const Vector6d &vec = v[vi];
+    for (size_t i = 0; i < 6; i++) {
+      m(vi, i) = vec[i];
+    }
+  }
+
+  return m;
+}
+
+// std::string convertToLowercase(const std::string &s) {
+//  std::string lowercase;
+//  std::transform(s.begin(), s.end(), lowercase.begin(),
+//                 [](unsigned char c) { return std::tolower(c); });
+//  return lowercase;
+//}
+// bool hasExtension(const std::string &filename, const std::string &extension)
+// {
+//  const std::filesystem::path path(filename);
+//  if (!path.has_extension()) {
+//    std::cerr << "Error: No file extension found in " << filename <<
+//    std::endl; return false;
+//  }
+
+//  const std::string detectedExtension = path.extension().string();
+
+//  if (convertToLowercase(detectedExtension) != convertToLowercase(extension))
+//  {
+//    std::cerr << "Error: detected extension is " + detectedExtension +
+//                     " and not " + extension
+//              << std::endl;
+//    return false;
+//  }
+//  return true;
+//}
+
+} // namespace Utilities
+
 // namespace ConfigurationFile {
 
 // inline void getPlyFilename(const std::string jsonFilepath,