diff --git a/edgemesh.cpp b/edgemesh.cpp
index c2b845f..87b0df2 100755
--- a/edgemesh.cpp
+++ b/edgemesh.cpp
@@ -1,5 +1,6 @@
 #include "edgemesh.hpp"
 #include "vcg/simplex/face/topology.h"
+#include <wrap/io_trimesh/import.h>
 //#include <wrap/nanoply/include/nanoplyWrapper.hpp>
 #include <wrap/io_trimesh/export.h>
 #include <wrap/io_trimesh/import.h>
@@ -26,7 +27,7 @@ Eigen::MatrixX3d VCGEdgeMesh::getEigenEdgeNormals() const {
   return eigenEdgeNormals;
 }
 
-bool VCGEdgeMesh::save(const std::filesystem::__cxx11::path &meshFilePath)
+bool VCGEdgeMesh::save(const std::filesystem::path &meshFilePath)
 {
     std::string filename = meshFilePath;
     if (filename.empty()) {
@@ -41,6 +42,9 @@ bool VCGEdgeMesh::save(const std::filesystem::__cxx11::path &meshFilePath)
     mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
     mask |= vcg::tri::io::Mask::IOM_VERTCOLOR;
     //    if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::SaveModel(filename.c_str(), *this, mask, false) != 0) {
+    if (std::filesystem::is_directory(meshFilePath.parent_path())) {
+        std::filesystem::create_directories(meshFilePath.parent_path());
+    }
     if (vcg::tri::io::Exporter<VCGEdgeMesh>::Save(*this, filename.c_str(), mask) != 0) {
         return false;
     }
@@ -247,7 +251,7 @@ bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath)
 //  return true;
 //}
 
-bool VCGEdgeMesh::loadUsingDefaultLoader(const std::string &plyFilename)
+bool VCGEdgeMesh::loadUsingDefaultLoader(const std::string &plyFilePath)
 {
     Clear();
     //  assert(plyFileHasAllRequiredFields(plyFilename));
@@ -259,7 +263,11 @@ bool VCGEdgeMesh::loadUsingDefaultLoader(const std::string &plyFilename)
     mask |= vcg::tri::io::Mask::IOM_EDGEINDEX;
     //  if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::LoadModel(plyFilename.c_str(),
     //                                                      *this, mask) != 0) {
-    if (vcg::tri::io::Importer<VCGEdgeMesh>::Open(*this, plyFilename.c_str(), mask) != 0) {
+    const int loadErrorCode = vcg::tri::io::Importer<VCGEdgeMesh>::Open(*this,
+                                                                        plyFilePath.c_str(),
+                                                                        mask);
+    if (loadErrorCode != 0) {
+        std::cerr << vcg::tri::io::Importer<VCGEdgeMesh>::ErrorMsg(loadErrorCode) << std::endl;
         return false;
     }
     return true;
@@ -358,7 +366,7 @@ void VCGEdgeMesh::set(const std::vector<double> &vertexPositions, const std::vec
         eIt->cV(0)->N() = n;
         eIt->cV(1)->N() = n;
     }
-    removeDuplicateVertices();
+    //    removeDuplicateVertices();
 
     updateEigenEdgeAndVertices();
 }
@@ -438,6 +446,22 @@ void VCGEdgeMesh::printVertexCoordinates(const size_t &vi) const {
 }
 
 #ifdef POLYSCOPE_DEFINED
+void VCGEdgeMesh::markVertices(const std::vector<size_t> &vertsToMark)
+{
+    if (vertsToMark.empty()) {
+        return;
+    }
+
+    std::vector<std::array<double, 3>> nodeColors(VN(), {0, 0, 0});
+    for (const size_t vi : vertsToMark) {
+        nodeColors[vi] = {1, 0, 0};
+    }
+
+    polyscope::getCurveNetwork(getLabel())
+        ->addNodeColorQuantity("Marked vertices" + getLabel(), nodeColors)
+        ->setEnabled(true);
+}
+
 //TODO: make const getEigenVertices is not
 polyscope::CurveNetwork *VCGEdgeMesh::registerForDrawing(
     const std::optional<std::array<double, 3>> &desiredColor,
diff --git a/edgemesh.hpp b/edgemesh.hpp
index de8cd51..73cfd35 100755
--- a/edgemesh.hpp
+++ b/edgemesh.hpp
@@ -3,11 +3,10 @@
 #include "beam.hpp"
 #include "mesh.hpp"
 #include "utilities.hpp"
+#include <array>
+#include <optional>
 #include <vcg/complex/complex.h>
 #include <vector>
-#include <wrap/io_trimesh/import.h>
-#include <optional>
-#include <array>
 
 #ifdef POLYSCOPE_DEFINED
 #include <polyscope/curve_network.h>
@@ -111,11 +110,13 @@ public:
         }
     }
 
+    void markVertices(const std::vector<size_t> &vertsToMark);
+
 private:
     void GeneratedRegularSquaredPattern(const double angleDeg,
                                         std::vector<std::vector<vcg::Point2d>> &pattern,
                                         const size_t &desiredNumberOfSamples);
-    bool loadUsingDefaultLoader(const std::string &plyFilename);
+    bool loadUsingDefaultLoader(const std::string &plyFilePath);
 };
 
 using VectorType = VCGEdgeMesh::CoordType;
diff --git a/reducedmodeloptimizer_structs.hpp b/reducedmodeloptimizer_structs.hpp
index 64c6f03..d68d8a0 100644
--- a/reducedmodeloptimizer_structs.hpp
+++ b/reducedmodeloptimizer_structs.hpp
@@ -109,17 +109,19 @@ inline int getParameterIndex(const std::string &s)
 struct Settings
 {
     inline static std::string defaultFilename{"OptimizationSettings.json"};
-    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{0,
-                                                                                  3,
-                                                                                  0.95,
-                                                                                  0.087,
-                                                                                  1.9,
-                                                                                  /*3*/ 0}; //custom
-    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{2.74429,
-    //                                                                                  5.9173,
-    //                                                                                  0.493706,
-    //                                                                                  0.0820769,
-    //                                                                                  1.57789};
+    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{0.590241,
+    //                                                                                  0.888372,
+    //                                                                                  0.368304,
+    //                                                                                  0.0127508,
+    //                                                                                  1.18079,
+    //                                                                                  0}; //final
+    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{0.590241 / 3,
+                                                                                  0.588372 / 3,
+                                                                                  0.368304,
+                                                                                  0.1,
+                                                                                  1.18 / 2,
+                                                                                  0}; //final
+    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{0, 0, 0, 0.1, 0};
 
     //    std::array<double, NumberOfBaseSimulationScenarios>
     //        baseScenarioMaxMagnitudes{20.85302947095844,
@@ -235,6 +237,8 @@ struct Settings
         bool operator==(const ObjectiveWeights &other) const;
     };
     std::array<ObjectiveWeights, NumberOfBaseSimulationScenarios> perBaseScenarioObjectiveWeights;
+    //    std::array<ObjectiveWeights, NumberOfBaseSimulationScenarios> perBaseScenarioObjectiveWeights{
+    //        {{1.95, 0.05}, {0.87, 1.13}, {0.37, 1.63}, {0.01, 1.99}, {0.94, 1.06}, {1.2, 0.8}}};
     std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios>
     convertObjectiveWeightsToPairs() const
     {
@@ -549,7 +553,7 @@ struct Settings
 
       PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern
       vcg::Triangle3<double> baseTriangle;
-      std::string notes;
+      //      std::string notes;
 
       //Data gathered for csv exporting
       struct ObjectiveValues
@@ -562,7 +566,9 @@ struct Settings
           std::vector<double> perSimulationScenario_rotational;
           std::vector<double> perSimulationScenario_total;
           std::vector<double> perSimulationScenario_total_unweighted;
+
       } objectiveValue;
+
       std::vector<double> perScenario_fullPatternPotentialEnergy;
       std::vector<double> objectiveValueHistory;
       std::vector<size_t> objectiveValueHistory_iteration;
@@ -573,7 +579,7 @@ struct Settings
       {
           bool exportPE{false};
           bool exportIterationOfMinima{false};
-          bool exportRawObjectiveValue{true};
+          bool exportRawObjectiveValue{false};
           CSVExportingSettings() {}
       };
 
@@ -589,6 +595,53 @@ struct Settings
           inline static std::string fullPatternYoungsModulus{"youngsModulus"};
       };
 
+      void saveObjectiveValuePlot(const std::filesystem::path &outputImageDirPath) const
+      {
+          std::vector<std::string> scenarioLabels(objectiveValue.perSimulationScenario_total.size());
+          const double colorAxial = 1;
+          const double colorShear = 3;
+          const double colorBending = 5;
+          const double colorDome = 0.1;
+          const double colorSaddle = 0;
+          std::vector<double> colors(objectiveValue.perSimulationScenario_total.size());
+          for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) {
+              scenarioLabels[scenarioIndex] = reducedPatternSimulationJobs[scenarioIndex]->getLabel();
+              if (scenarioLabels[scenarioIndex].rfind("Axial", 0) == 0) {
+                  colors[scenarioIndex] = colorAxial;
+
+              } else if (scenarioLabels[scenarioIndex].rfind("Shear", 0) == 0) {
+                  colors[scenarioIndex] = colorShear;
+
+              } else if (scenarioLabels[scenarioIndex].rfind("Bending", 0) == 0) {
+                  colors[scenarioIndex] = colorBending;
+
+              } else if (scenarioLabels[scenarioIndex].rfind("Dome", 0) == 0) {
+                  colors[scenarioIndex] = colorDome;
+              } else if (scenarioLabels[scenarioIndex].rfind("Saddle", 0) == 0) {
+                  colors[scenarioIndex] = colorSaddle;
+              } else {
+                  std::cerr << "Label could not be identified" << std::endl;
+              }
+          }
+          std::vector<double> y(objectiveValue.perSimulationScenario_total.size());
+          for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) {
+              y[scenarioIndex]
+                  //          = optimizationResults.objectiveValue.perSimulationScenario_rawTranslational[scenarioIndex]
+                  //            + optimizationResults.objectiveValue.perSimulationScenario_rawRotational[scenarioIndex];
+                  = objectiveValue.perSimulationScenario_total_unweighted[scenarioIndex];
+          }
+          std::vector<double> x = matplot::linspace(0, y.size() - 1, y.size());
+          std::vector<double> markerSizes(y.size(), 5);
+          Utilities::createPlot("scenario index",
+                                "error",
+                                x,
+                                y,
+                                markerSizes,
+                                colors,
+                                std::filesystem::path(outputImageDirPath)
+                                    .append("perScenarioObjectiveValues.svg"));
+      }
+
       void save(const std::string &saveToPath, const bool shouldExportDebugFiles = false)
       {
           //clear directory
@@ -632,14 +685,14 @@ struct Settings
           json_optimizationResults[JsonKeys::FullPatternLabel] = baseTriangleFullPattern.getLabel();
 
           //potential energies
-          const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
-          std::vector<double> fullPatternPE(numberOfSimulationJobs);
-          for (int simulationScenarioIndex = 0; simulationScenarioIndex < numberOfSimulationJobs;
-               simulationScenarioIndex++) {
-              fullPatternPE[simulationScenarioIndex]
-                  = perScenario_fullPatternPotentialEnergy[simulationScenarioIndex];
-          }
-          json_optimizationResults[JsonKeys::FullPatternPotentialEnergies] = fullPatternPE;
+          //          const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
+          //          std::vector<double> fullPatternPE(numberOfSimulationJobs);
+          //          for (int simulationScenarioIndex = 0; simulationScenarioIndex < numberOfSimulationJobs;
+          //               simulationScenarioIndex++) {
+          //              fullPatternPE[simulationScenarioIndex]
+          //                  = perScenario_fullPatternPotentialEnergy[simulationScenarioIndex];
+          //          }
+          //          json_optimizationResults[JsonKeys::FullPatternPotentialEnergies] = fullPatternPE;
           json_optimizationResults[JsonKeys::fullPatternYoungsModulus] = fullPatternYoungsModulus;
           ////Save to json file
           std::filesystem::path jsonFilePath(
@@ -678,13 +731,28 @@ struct Settings
                   }
                   pReducedPatternSimulationJob->save(reducedPatternDirectoryPath.string());
               }
+
+              //              constexpr bool shouldSaveObjectiveValuePlot = shouldExportDebugFiles;
+              //              if (shouldSaveObjectiveValuePlot) {
+              saveObjectiveValuePlot(saveToPath);
+              //              }
           }
+          csvFile csv_resultsLocalFile(std::filesystem::path(saveToPath).append("results.csv"),
+                                       true);
+          csv_resultsLocalFile << "Name";
+          writeHeaderTo(csv_resultsLocalFile);
+          settings.writeHeaderTo(csv_resultsLocalFile);
+          csv_resultsLocalFile << endrow;
+          csv_resultsLocalFile << baseTriangleFullPattern.getLabel();
+          writeResultsTo(csv_resultsLocalFile);
+          settings.writeSettingsTo(csv_resultsLocalFile);
+          csv_resultsLocalFile << endrow;
 
           //save minima info
-          std::filesystem::path csvFilepathMinimaInfo = std::filesystem::path(saveToPath)
-                                                            .append("minimaInfo.csv");
-          csvFile csv_minimaInfo(csvFilepathMinimaInfo, false);
-          writeMinimaInfoTo(csv_minimaInfo);
+          //          std::filesystem::path csvFilepathMinimaInfo = std::filesystem::path(saveToPath)
+          //                                                            .append("minimaInfo.csv");
+          //          csvFile csv_minimaInfo(csvFilepathMinimaInfo, false);
+          //          writeMinimaInfoTo(csv_minimaInfo);
 
           settings.save(saveToPath);
 
@@ -707,6 +775,8 @@ struct Settings
           std::ifstream ifs(std::filesystem::path(loadFromPath).append(DefaultFileName));
           ifs >> json_optimizationResults;
 
+          //          std::cout << json_optimizationResults.dump() << std::endl;
+
           label = json_optimizationResults.at(JsonKeys::Label);
           std::string optimizationVariablesString = *json_optimizationResults.find(
               JsonKeys::optimizationVariables);
@@ -732,8 +802,12 @@ struct Settings
 
           const std::string fullPatternLabel = json_optimizationResults.at(
               JsonKeys::FullPatternLabel);
-          baseTriangleFullPattern.load(
-              std::filesystem::path(loadFromPath).append(fullPatternLabel + ".ply").string());
+          if (!baseTriangleFullPattern.load(
+                  std::filesystem::path(loadFromPath).append(fullPatternLabel + ".ply").string())) {
+              !baseTriangleFullPattern.load(std::filesystem::path(loadFromPath)
+                                                .append(loadFromPath.stem().string() + ".ply")
+                                                .string());
+          }
 
           std::vector<double> baseTriangleVertices = json_optimizationResults.at(
               JsonKeys::baseTriangle);
@@ -791,7 +865,7 @@ struct Settings
       }
 
       template<typename MeshType>
-      static void applyOptimizationResults_innerHexagon(
+      static void applyOptimizationResults_reducedModel_nonFanned(
           const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
           const vcg::Triangle3<double> &patternBaseTriangle,
           MeshType &reducedPattern)
@@ -803,20 +877,20 @@ struct Settings
               (optimalXVariables.contains("R") && optimalXVariables.contains("Theta"))
               || (optimalXVariables.contains("HexSize") && optimalXVariables.contains("HexAngle")));
           if (optimalXVariables.contains("HexSize")) {
-              applyOptimizationResults_innerHexagon(optimalXVariables.at("HexSize"),
-                                                    optimalXVariables.at("HexAngle"),
-                                                    patternBaseTriangle,
-                                                    reducedPattern);
+              applyOptimizationResults_reducedModel_nonFanned(optimalXVariables.at("HexSize"),
+                                                              optimalXVariables.at("HexAngle"),
+                                                              patternBaseTriangle,
+                                                              reducedPattern);
               return;
           }
-          applyOptimizationResults_innerHexagon(optimalXVariables.at("R"),
-                                                optimalXVariables.at("Theta"),
-                                                patternBaseTriangle,
-                                                reducedPattern);
+          applyOptimizationResults_reducedModel_nonFanned(optimalXVariables.at("R"),
+                                                          optimalXVariables.at("Theta"),
+                                                          patternBaseTriangle,
+                                                          reducedPattern);
       }
 
       template<typename MeshType>
-      static void applyOptimizationResults_innerHexagon(
+      static void applyOptimizationResults_reducedModel_nonFanned(
           const double &hexSize,
           const double &hexAngle,
           const vcg::Triangle3<double> &patternBaseTriangle,
diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp
index f49612b..f4f329f 100755
--- a/simulationhistoryplotter.hpp
+++ b/simulationhistoryplotter.hpp
@@ -85,23 +85,6 @@ struct SimulationResultsReporter {
       writeStatistics(simulationResult, simulationResultPath.string());
     }
   }
-  static void createPlot(const std::string &xLabel,
-                         const std::string &yLabel,
-                         const std::vector<double> &x,
-                         const std::vector<double> &y,
-                         const std::vector<double> &markerSizes,
-                         const std::vector<double> &c,
-                         const std::string &saveTo = {})
-  {
-      //    matplot::figure(true);
-      matplot::xlabel(xLabel);
-      matplot::ylabel(yLabel);
-      matplot::grid(matplot::on);
-      matplot::scatter(x, y, markerSizes, c)->marker_face(true);
-      if (!saveTo.empty()) {
-          matplot::save(saveTo);
-      }
-  }
 
   static void createPlot(const std::string &xLabel,
                          const std::string &yLabel,
@@ -121,7 +104,7 @@ struct SimulationResultsReporter {
           }
       }
       std::vector<double> x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size());
-      createPlot(xLabel, yLabel, x, YvaluesToPlot, markerSizes, colors, saveTo);
+      Utilities::createPlot(xLabel, yLabel, x, YvaluesToPlot, markerSizes, colors, saveTo);
   }
 
   void createPlots(const SimulationHistory &history,
@@ -178,10 +161,10 @@ struct SimulationResultsReporter {
       //                       .append("ResidualForcesMovingAverageDerivativesLog_" + graphSuffix + ".png")
       //                       .string());
       //    }
-      if (!history.sumOfNormalizedDisplacementNorms.empty()) {
+      if (!history.perVertexAverageNormalizedDisplacementNorm.empty()) {
           createPlot("Number of Iterations",
                      "Sum of normalized displacement norms",
-                     history.sumOfNormalizedDisplacementNorms,
+                     history.perVertexAverageNormalizedDisplacementNorm,
                      std::filesystem::path(graphsFolder)
                          .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png")
                          .string(),
diff --git a/simulationmesh.hpp b/simulationmesh.hpp
index ecdcb52..d6c8d57 100755
--- a/simulationmesh.hpp
+++ b/simulationmesh.hpp
@@ -39,9 +39,12 @@ public:
   std::vector<ElementMaterial> getBeamMaterial();
 
   double previousTotalKineticEnergy{0};
+  double previousTranslationalKineticEnergy{0};
+  double previousTotalRotationalKineticEnergy{0};
   double previousTotalResidualForcesNorm{0};
   double currentTotalKineticEnergy{0};
   double currentTotalTranslationalKineticEnergy{0};
+  double currentTotalRotationalKineticEnergy{0};
   double totalResidualForcesNorm{0};
   double totalExternalForcesNorm{0};
   double averageResidualForcesNorm{0};
@@ -49,7 +52,7 @@ public:
   double previousTotalPotentialEnergykN{0};
   double residualForcesMovingAverageDerivativeNorm{0};
   double residualForcesMovingAverage{0};
-  double sumOfNormalizedDisplacementNorms{0};
+  double perVertexAverageNormalizedDisplacementNorm{0};
   bool save(const std::string &plyFilename = std::string());
   void setBeamCrossSection(const CrossSectionType &beamDimensions);
   void setBeamMaterial(const double &pr, const double &ym);
@@ -133,6 +136,7 @@ struct Node {
     Vector6d internal{0};
     Vector6d residual{0};
     Vector6d internalAxial{0};
+    Vector6d internalTorsion{0};
     Vector6d internalFirstBending{0};
     Vector6d internalSecondBending{0};
     bool hasExternalForce() const { return external.isZero(); }
diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp
index 21d8342..d785670 100755
--- a/topologyenumerator.cpp
+++ b/topologyenumerator.cpp
@@ -441,11 +441,12 @@ void TopologyEnumerator::computeValidPatterns(
     //  std::string previousPatternBinaryRepresentation(validEdges.size(),'0');
     size_t patternIndex = 0;
     bool validPatternsExist = false;
-    const bool exportTilledPattern = true;
-    const bool saveCompressedFormat = false;
+    constexpr bool exportTilledPattern = false;
+    constexpr bool saveCompressedFormat = false;
     do {
         patternIndex++;
-        const std::string patternName = std::to_string(patternIndex);
+        const std::string patternName = std::to_string(numberOfDesiredEdges) + "_"
+                                        + std::to_string(patternIndex);
         //    std::cout << "Pattern name:" + patternBinaryRepresentation <<
         //    std::endl; isValidPattern(patternBinaryRepresentation, validEdges,
         //                   numberOfDesiredEdges);
@@ -465,6 +466,32 @@ void TopologyEnumerator::computeValidPatterns(
         patternGeometry.add(allVertices, patternEdges);
         patternGeometry.setLabel(patternName);
 
+#ifdef POLYSCOPE_DEFINED
+        //1st example
+        //        const bool shouldBreak = patternBinaryRepresentation == "00100000100100000"; //398
+        //        const bool shouldBreak = patternBinaryRepresentation == "10000010101110110";//13036
+        //        const bool shouldBreak = patternBinaryRepresentation == "00010111000010100"; //2481
+        //        const bool shouldBreak = patternBinaryRepresentation == "10000101100110010"; //12116
+        //        const bool shouldBreak = patternBinaryRepresentation == "10010111000000110"; //13915
+        //2nd example
+        //        const bool shouldBreak = patternBinaryRepresentation == "00001011100010011"; //7_1203
+        //        const bool shouldBreak = patternBinaryRepresentation == "00110001100100111"; //4865
+        //        const bool shouldBreak = patternBinaryRepresentation == "00010000101000110"; //1380
+        //        const bool shouldBreak = patternBinaryRepresentation == "00000010100010111"; //268
+        //3rd
+        //        const bool shouldBreak = patternBinaryRepresentation == "10011011100000010"; //14272
+        //        const bool shouldBreak = patternBinaryRepresentation == "10000111100110110"; //11877
+        //        const bool shouldBreak = patternBinaryRepresentation == "00001011100010011"; //1203
+        //        const bool shouldBreak = patternBinaryRepresentation == "00010101000110000"; //12117
+
+        //        const bool shouldBreak = patternBinaryRepresentation == "10000101100110100"; //12117
+//        if (shouldBreak) {
+//            patternGeometry.registerForDrawing();
+//            polyscope::show();
+//            patternGeometry.unregister();
+//        }
+#endif
+
         // Check if pattern contains intersecting edges
         const bool isInterfaceConnected = patternGeometry.isInterfaceConnected(interfaceNodes);
         // Export the tiled ply file if it contains intersecting edges
diff --git a/trianglepatterngeometry.cpp b/trianglepatterngeometry.cpp
index 163c2fb..6a88cbb 100755
--- a/trianglepatterngeometry.cpp
+++ b/trianglepatterngeometry.cpp
@@ -427,8 +427,9 @@ bool PatternGeometry::hasAngleSmallerThanThreshold(const std::vector<size_t> &nu
         //                                      (vert[vi].cP() + tiledIncidentVectors[tiledVectorIndex])[1],
         //                                      (vert[vi].cP() + tiledIncidentVectors[tiledVectorIndex])[2]});
         //        }
+        //        polyscope::init();
         //        polyscope::registerCurveNetworkLine("temp", edgePoints);
-        //        polyscope::show();
+        //                polyscope::removeStructure("temp");
         if (tiledIncidentVectors.size() == 1) {
             continue;
         }
@@ -440,14 +441,28 @@ bool PatternGeometry::hasAngleSmallerThanThreshold(const std::vector<size_t> &nu
                        [](const VectorType &v) { return vcg::Point2d(v[0], v[1]).Angle(); });
         //sort them using theta angles
         std::sort(thetaAnglesOfIncidentVectors.begin(), thetaAnglesOfIncidentVectors.end());
+        //        polyscope::show();
+
+        //        std::vector<double> angles_theta(thetaAnglesOfIncidentVectors);
+        //        for (double &theta_rad : angles_theta) {
+        //            theta_rad = vcg::math::ToDeg(theta_rad);
+        //        }
 
         //find nodes that contain incident edges with relative angles less than the threshold
         const double angleThreshold_rad = vcg::math::ToRad(angleThreshold_degrees);
-        for (int thetaAngleIndex = 1; thetaAngleIndex < thetaAnglesOfIncidentVectors.size();
+        for (int thetaAngleIndex = 0; thetaAngleIndex < thetaAnglesOfIncidentVectors.size();
              thetaAngleIndex++) {
-            const double absAngleDifference = std::abs(
-                thetaAnglesOfIncidentVectors[thetaAngleIndex]
-                - thetaAnglesOfIncidentVectors[thetaAngleIndex - 1]);
+            const auto &va_theta
+                = thetaAnglesOfIncidentVectors[(thetaAngleIndex + 1)
+                                               % thetaAnglesOfIncidentVectors.size()];
+            const auto &vb_theta = thetaAnglesOfIncidentVectors[thetaAngleIndex];
+            //            const auto &va
+            //                = tiledIncidentVectors[(thetaAngleIndex + 1) % thetaAnglesOfIncidentVectors.size()];
+            //            const auto &vb = tiledIncidentVectors[thetaAngleIndex];
+            const double absAngleDifference = std::abs(va_theta - vb_theta);
+            //            const double debug_difDegOtherway = vcg::math::ToDeg(
+            //                std::acos((va * vb) / (va.Norm() * vb.Norm())));
+            //            const double debug_diffDeg = vcg::math::ToDeg(absAngleDifference);
             if (absAngleDifference < angleThreshold_rad
                 /*&& absAngleDifference > vcg::math::ToRad(0.01)*/) {
                 //                std::cout << "Found angDiff:" << absAngleDifference << std::endl;
@@ -790,13 +805,17 @@ double PatternGeometry::computeBaseTriangleHeight() const
     return vcg::Distance(vert[0].cP(), vert[interfaceNodeIndex].cP());
 }
 
+void PatternGeometry::updateBaseTriangleHeight()
+{
+    baseTriangleHeight = computeBaseTriangleHeight();
+}
+
 void PatternGeometry::deleteDanglingVertices()
 {
     vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> pu;
     VCGEdgeMesh::deleteDanglingVertices(pu);
     if (!pu.remap.empty()) {
-        interfaceNodeIndex
-            = pu.remap[interfaceNodeIndex]; //TODO:Could this be automatically be determined?
+        interfaceNodeIndex = pu.remap[interfaceNodeIndex];
     }
 }
 
@@ -840,6 +859,8 @@ PatternGeometry::PatternGeometry(
   addNormals();
   baseTriangleHeight = computeBaseTriangleHeight();
   baseTriangle = computeBaseTriangle();
+  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
+  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
   updateEigenEdgeAndVertices();
 }
 
@@ -914,6 +935,9 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
             fit->N() = faceNormal;
             PatternGeometry transformedPattern;
             transformedPattern.copy(pattern);
+            //            pattern.registerForDrawing();
+            //            polyscope::show();
+            //            pattern.unregister();
             //Transform the base triangle nodes to the mesh triangle using barycentric coords
             for (int vi = 0; vi < transformedPattern.VN(); vi++) {
                 transformedPattern.vert[vi].P() = CoordType(
@@ -933,6 +957,10 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
             for (size_t ei = 0; ei < pattern.EN(); ei++) {
                 perPatternIndexToTiledPatternEdgeIndex[patternIndex].push_back(remap.edge[ei]);
             }
+            //            pTiledPattern->registerForDrawing();
+            //            pTiledPattern->markVertices({remap.vert[pattern.interfaceNodeIndex]});
+            //            polyscope::show();
+            //            pTiledPattern->unregister();
             const size_t ei = tileInto.getIndex(ep);
             tileIntoEdgeToInterfaceVi[ei].push_back(remap.vert[pattern.interfaceNodeIndex]);
             //Add edges for connecting the desired vertices
@@ -965,10 +993,6 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
                     }
                 }
             }
-            //            pTiledPattern->updateEigenEdgeAndVertices();
-            //            pTiledPattern->registerForDrawing();
-
-            //            polyscope::show();
         }
     }
     vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> pu_vertices;
@@ -987,6 +1011,7 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
 
     //        tiledPatternEdges.erase(end, tiledPatternEdges.end());
     //    }
+
     const size_t sumOfEdgeIndices = std::accumulate(perPatternIndexToTiledPatternEdgeIndex.begin(),
                                                     perPatternIndexToTiledPatternEdgeIndex.end(),
                                                     0,
@@ -1017,7 +1042,7 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
     pTiledPattern->deleteDanglingVertices();
     vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(*pTiledPattern);
     pTiledPattern->updateEigenEdgeAndVertices();
-    pTiledPattern->save();
+    //    pTiledPattern->save();
     return pTiledPattern;
 }
 
@@ -1078,17 +1103,21 @@ void PatternGeometry::copy(PatternGeometry &copyFrom)
     interfaceNodeIndex = copyFrom.interfaceNodeIndex;
 }
 
-void PatternGeometry::scale(const double &desiredBaseTriangleCentralEdgeSize,
-                            const int &interfaceNodeIndex)
+void PatternGeometry::scale(const double &desiredBaseTriangleCentralEdgeSize)
 {
-    const double baseTriangleCentralEdgeSize = computeBaseTriangleHeight();
+    const double baseTriangleCentralEdgeSize = getBaseTriangleHeight();
     const double scaleRatio = desiredBaseTriangleCentralEdgeSize / baseTriangleCentralEdgeSize;
     vcg::tri::UpdatePosition<VCGEdgeMesh>::Scale(*this, scaleRatio);
     baseTriangle = computeBaseTriangle();
     baseTriangleHeight = computeBaseTriangleHeight();
+    const double debug_baseTriHeight = vcg::Distance(baseTriangle.cP(0),
+                                                     (baseTriangle.cP(1) + baseTriangle.cP(2)) / 2);
+    assert(std::abs(desiredBaseTriangleCentralEdgeSize - baseTriangleHeight) < 1e-10);
+    int i = 0;
+    i++;
 }
 
-void PatternGeometry::createFan(const std::vector<int> &connectToNeighborsVi, const size_t &fanSize)
+void PatternGeometry::createFan(const size_t &fanSize)
 {
     PatternGeometry rotatedPattern;
     vcg::tri::Append<PatternGeometry, PatternGeometry>::MeshCopy(rotatedPattern, *this);
@@ -1099,28 +1128,6 @@ void PatternGeometry::createFan(const std::vector<int> &connectToNeighborsVi, co
         vcg::tri::UpdatePosition<PatternGeometry>::Matrix(rotatedPattern, R);
         vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(*this, rotatedPattern);
         //Add edges for connecting the desired vertices
-        if (!connectToNeighborsVi.empty()) {
-            if (rotationCounter == fanSize - 1) {
-                for (int connectToNeighborIndex = 0;
-                     connectToNeighborIndex < connectToNeighborsVi.size();
-                     connectToNeighborIndex++) {
-                    vcg::tri::Allocator<PatternGeometry>::AddEdge(
-                        *this,
-                        connectToNeighborsVi[connectToNeighborIndex],
-                        this->VN() - rotatedPattern.VN()
-                            + connectToNeighborsVi[connectToNeighborIndex]);
-                }
-            }
-            for (int connectToNeighborIndex = 0;
-                 connectToNeighborIndex < connectToNeighborsVi.size();
-                 connectToNeighborIndex++) {
-                vcg::tri::Allocator<PatternGeometry>::AddEdge(
-                    *this,
-                    this->VN() - 2 * rotatedPattern.VN()
-                        + connectToNeighborsVi[connectToNeighborIndex],
-                    this->VN() - rotatedPattern.VN() + connectToNeighborsVi[connectToNeighborIndex]);
-            }
-        }
         removeDuplicateVertices();
         updateEigenEdgeAndVertices();
     }
diff --git a/trianglepatterngeometry.hpp b/trianglepatterngeometry.hpp
index 249af04..c4241d2 100755
--- a/trianglepatterngeometry.hpp
+++ b/trianglepatterngeometry.hpp
@@ -21,7 +21,6 @@ private:
 
   void addNormals();
   double baseTriangleHeight;
-  double computeBaseTriangleHeight() const;
 
   inline static size_t fanSize{6};
   std::vector<VCGEdgeMesh::CoordType> vertices;
@@ -95,20 +94,23 @@ private:
                    const int &interfaceNodeIndex,
                    const bool &shouldDeleteDanglingEdges);
 
-  void scale(const double &desiredBaseTriangleCentralEdgeSize, const int &interfaceNodeIndex);
+  void scale(const double &desiredBaseTriangleCentralEdgeSize);
 
   double getBaseTriangleHeight() const;
   vcg::Triangle3<double> computeBaseTriangle() const;
+  void updateBaseTriangle();
+  double computeBaseTriangleHeight() const;
+  void updateBaseTriangleHeight();
 
-  PatternGeometry(const std::vector<vcg::Point2d> &vertices, const std::vector<vcg::Point2i> &edges);
+  //  PatternGeometry(const std::vector<vcg::Point2d> &vertices, const std::vector<vcg::Point2i> &edges);
   //  static std::shared_ptr<PatternGeometry> tilePattern(
   //      std::vector<PatternGeometry> &pattern,
   //      const std::vector<int> &connectToNeighborsVi,
   //      const VCGPolyMesh &tileInto,
   //      std::vector<int> &tileIntoEdgesToTiledVi,
   //      std::vector<std::vector<size_t>> &perPatternEdges);
-  void createFan(const std::vector<int> &connectToNeighborsVi = std::vector<int>(),
-                 const size_t &fanSize = 6);
+  virtual void createFan(/*const std::vector<int> &connectToNeighborsVi = std::vector<int>(),*/
+                         const size_t &fanSize = 6);
   int interfaceNodeIndex{3}; //TODO: Fix this. This should be automatically computed
   bool hasAngleSmallerThanThreshold(const std::vector<size_t> &numberOfNodesPerSlot,
                                     const double &angleThreshold_degrees);
@@ -129,7 +131,6 @@ private:
   void deleteDanglingVertices() override;
   void deleteDanglingVertices(
       vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu) override;
-  void updateBaseTriangle();
 };
 
 #endif // FLATPATTERNGEOMETRY_HPP