diff --git a/beam.hpp b/beam.hpp
index f629662..1ec22a9 100644
--- a/beam.hpp
+++ b/beam.hpp
@@ -6,9 +6,9 @@
 
 struct RectangularBeamDimensions {
   inline static std::string name{"Rectangular"};
-  float b;
-  float h;
-  RectangularBeamDimensions(const float &width, const float &height)
+  double b;
+  double h;
+  RectangularBeamDimensions(const double &width, const double &height)
       : b(width), h(height) {
     assert(width > 0 && height > 0);
   }
@@ -17,12 +17,12 @@ struct RectangularBeamDimensions {
 
 struct CylindricalBeamDimensions {
   inline static std::string name{"Cylindrical"};
-  float od; // Cylinder outside diameter
-  float
+  double od; // Cylinder outside diameter
+  double
       id; // Cylinder inside diameter
           // https://www.engineeringtoolbox.com/area-moment-inertia-d_1328.html
-  CylindricalBeamDimensions(const float &outsideDiameter,
-                            const float &insideDiameter)
+  CylindricalBeamDimensions(const double &outsideDiameter,
+                            const double &insideDiameter)
       : od(outsideDiameter), id(insideDiameter) {
     assert(outsideDiameter > 0 && insideDiameter > 0 &&
            outsideDiameter > insideDiameter);
@@ -31,14 +31,20 @@ struct CylindricalBeamDimensions {
 };
 
 struct ElementMaterial {
-  float poissonsRatio; // NOTE: if I make this double the unit
-                       // tests produced different results and thus fail
+  double poissonsRatio; // NOTE: if I make this double the unit
+                        // tests produced different results and thus fail
   double youngsModulus;
-  ElementMaterial(const float &poissonsRatio, const double &youngsModulus)
+  ElementMaterial(const double &poissonsRatio, const double &youngsModulus)
       : poissonsRatio(poissonsRatio), youngsModulus(youngsModulus) {
     assert(poissonsRatio <= 0.5 && poissonsRatio >= -1);
   }
   ElementMaterial() : poissonsRatio(0.3), youngsModulus(200) {}
+  std::string toString() const {
+    return std::string("Material:") + std::string("\nPoisson's ratio=") +
+           std::to_string(poissonsRatio) +
+           std::string("\nYoung's Modulus(GPa)=") +
+           std::to_string(youngsModulus / 1e9);
+  }
 };
 
 #endif // BEAM_HPP
diff --git a/beamformfinder.cpp b/beamformfinder.cpp
index f0e22ea..975257d 100644
--- a/beamformfinder.cpp
+++ b/beamformfinder.cpp
@@ -1865,9 +1865,12 @@ FormFinder::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
       std::cout << "Non terminating simulation found. Saved simulation job to:"
                 << dir << std::endl;
       std::cout << "Exiting.." << std::endl;
-      //      FormFinder debug;
-      //      debug.executeSimulation(pJob, true, true, true);
-      //      std::terminate();
+      FormFinder debug;
+      FormFinder::Settings settings;
+      settings.shouldDraw = true;
+      settings.beVerbose = true;
+      debug.executeSimulation(pJob, settings);
+      std::terminate();
       break;
     }
 
diff --git a/elementalmesh.cpp b/elementalmesh.cpp
index e8da8e1..5f02a8e 100644
--- a/elementalmesh.cpp
+++ b/elementalmesh.cpp
@@ -333,11 +333,11 @@ bool SimulationMesh::savePly(const std::string &plyFilename) {
   nanoply::NanoPlyWrapper<VCGEdgeMesh>::CustomAttributeDescriptor customAttrib;
   customAttrib.GetMeshAttrib(plyFilename);
 
-  dimensions = getBeamDimensions();
+  std::vector<CrossSectionType> dimensions = getBeamDimensions();
   customAttrib.AddEdgeAttribDescriptor<CrossSectionType, float, 2>(
       plyPropertyBeamDimensionsID, nanoply::NNP_LIST_INT8_FLOAT32,
       dimensions.data());
-  material = getBeamMaterial();
+  std::vector<ElementMaterial> material = getBeamMaterial();
   customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>(
       plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT32,
       material.data());
diff --git a/elementalmesh.hpp b/elementalmesh.hpp
index 25cc456..31909cc 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:
@@ -20,8 +20,8 @@ private:
   const std::string plyPropertyBeamDimensionsID{"beam_dimensions"};
   const std::string plyPropertyBeamMaterialID{"beam_material"};
   const std::string plyPropertyBeamProperties{"beam_properties"};
-  std::vector<ElementMaterial> material;
-  std::vector<CrossSectionType> dimensions;
+  //  std::vector<ElementMaterial> material;
+  //  std::vector<CrossSectionType> dimensions;
 
 public:
   PerEdgeAttributeHandle<Element> elements;
diff --git a/formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin b/formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin
index 6649189..cd8534c 100644
Binary files a/formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin and b/formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin differ
diff --git a/formFinder_unitTestFiles/ShortSpanGridshell_displacements.eigenBin b/formFinder_unitTestFiles/ShortSpanGridshell_displacements.eigenBin
index 19a07ed..ed58e09 100644
Binary files a/formFinder_unitTestFiles/ShortSpanGridshell_displacements.eigenBin and b/formFinder_unitTestFiles/ShortSpanGridshell_displacements.eigenBin differ
diff --git a/formFinder_unitTestFiles/SimpleBeam_displacements.eigenBin b/formFinder_unitTestFiles/SimpleBeam_displacements.eigenBin
index 61a9a23..335be93 100644
Binary files a/formFinder_unitTestFiles/SimpleBeam_displacements.eigenBin and b/formFinder_unitTestFiles/SimpleBeam_displacements.eigenBin differ