diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp
index 3b68b57..64712da 100755
--- a/drmsimulationmodel.cpp
+++ b/drmsimulationmodel.cpp
@@ -889,25 +889,26 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
                 const Node &edgeNode = pMesh->nodes[ev];
                 internalForcesContributionFromThisEdge[evi].first = edgeNode.vi;
 
-            const VertexPointer &rev_j = edgeNode.referenceElement->cV(0);
-            const VertexPointer &rev_jplus1 = edgeNode.referenceElement->cV(1);
-            // refElemOtherVertex can be precomputed
-            const VertexPointer &refElemOtherVertex = rev_j == &ev ? rev_jplus1 : rev_j;
-            const Node &refElemOtherVertexNode = pMesh->nodes[refElemOtherVertex];
-            if (edgeNode.referenceElement != &e) {
-                internalForcesContributionFromThisEdge[evi + 2].first = refElemOtherVertexNode.vi;
-            }
-            const size_t vi = edgeNode.vi;
-            //      #pragma omp parallel for schedule(static) num_threads(6)
-            for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
-                const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi]
-                                                    && fixedVertices.at(edgeNode.vi).contains(dofi);
-                if (!isDofConstrainedFor_ev) {
-                    DifferentiateWithRespectTo dui{ev, dofi};
-                    // Axial force computation
-                    const double e_k = element.length - element.initialLength;
-                    const double e_kDeriv = computeDerivativeElementLength(e, dui);
-                    const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial;
+                const VertexPointer &rev_j = edgeNode.referenceElement->cV(0);
+                const VertexPointer &rev_jplus1 = edgeNode.referenceElement->cV(1);
+                // refElemOtherVertex can be precomputed
+                const VertexPointer &refElemOtherVertex = rev_j == &ev ? rev_jplus1 : rev_j;
+                const Node &refElemOtherVertexNode = pMesh->nodes[refElemOtherVertex];
+                if (edgeNode.referenceElement != &e) {
+                    internalForcesContributionFromThisEdge[evi + 2].first = refElemOtherVertexNode.vi;
+                }
+                const size_t vi = edgeNode.vi;
+                //      #pragma omp parallel for schedule(static) num_threads(6)
+                for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
+                    const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi]
+                                                        && fixedVertices.at(edgeNode.vi)
+                                                               .contains(dofi);
+                    if (!isDofConstrainedFor_ev) {
+                        DifferentiateWithRespectTo dui{ev, dofi};
+                        // Axial force computation
+                        const double e_k = element.length - element.initialLength;
+                        const double e_kDeriv = computeDerivativeElementLength(e, dui);
+                        const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial;
 
 #ifdef POLYSCOPE_DEFINED
                         if (std::isnan(axialForce_dofi)) {
@@ -984,7 +985,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
                         internalForcesContributionFromThisEdge[evi].second[dofi]
                             = axialForce_dofi + firstBendingForce_dofi + secondBendingForce_dofi
                               + torsionalForce_dofi;
-                }
+                    }
                     if (edgeNode.referenceElement != &e) {
                         const bool isDofConstrainedFor_refElemOtherVertex
                             = isVertexConstrained[refElemOtherVertexNode.vi]
@@ -1072,7 +1073,6 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
         }
 #endif
     }
-    pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm;
     pMesh->totalResidualForcesNorm = totalResidualForcesNorm;
     pMesh->averageResidualForcesNorm = totalResidualForcesNorm / pMesh->VN();
     //    pMesh->averageResidualForcesNorm = sumOfResidualForces.norm() / pMesh->VN();
@@ -2736,54 +2736,52 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
 
 //#endif
 
-void DRMSimulationModel::Settings::save(const std::filesystem::path &folderPath) const
+void DRMSimulationModel::Settings::save(const std::filesystem::path &jsonFilePath) const
 {
-    std::filesystem::create_directories(folderPath);
     nlohmann::json json;
-
-    json[jsonLabels.gamma] = gamma;
-    json[jsonLabels.shouldDraw] = shouldDraw ? "true" : "false";
-    json[jsonLabels.shouldCreatePlots] = shouldCreatePlots ? "true" : "false";
-    json[jsonLabels.beVerbose] = beVerbose ? "true" : "false";
-    json[jsonLabels.Dtini] = Dtini;
-    json[jsonLabels.xi] = xi;
+    json[GET_VARIABLE_NAME(shouldDraw)] = shouldDraw;
+    json[GET_VARIABLE_NAME(beVerbose)] = beVerbose;
+    json[GET_VARIABLE_NAME(shouldCreatePlots)] = shouldCreatePlots;
+    json[GET_VARIABLE_NAME(Dtini)] = Dtini;
+    json[GET_VARIABLE_NAME(xi)] = xi;
+    json[GET_VARIABLE_NAME(gamma)] = gamma;
+    json[GET_VARIABLE_NAME(useKineticDamping)] = totalResidualForcesNormThreshold;
     if (maxDRMIterations.has_value()) {
-        json[jsonLabels.maxDRMIterations] = maxDRMIterations.value();
+        json[GET_VARIABLE_NAME(jsonLabels.maxDRMIterations)] = maxDRMIterations.value();
     }
     if (debugModeStep.has_value()) {
-        json[jsonLabels.debugModeStep] = debugModeStep.value();
+        json[GET_VARIABLE_NAME(debugModeStep)] = debugModeStep.value();
     }
     if (totalTranslationalKineticEnergyThreshold.has_value()) {
-        json[jsonLabels.totalTranslationalKineticEnergyThreshold]
+        json[GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold)]
             = totalTranslationalKineticEnergyThreshold.value();
     }
-
     if (averageResidualForcesCriterionThreshold.has_value()) {
-        json[jsonLabels.averageResidualForcesCriterionThreshold]
-            = averageResidualForcesCriterionThreshold.value();
-    }
-
-    if (averageResidualForcesCriterionThreshold.has_value()) {
-        json[jsonLabels.averageResidualForcesCriterionThreshold]
+        json[GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold)]
             = averageResidualForcesCriterionThreshold.value();
     }
     if (linearGuessForceScaleFactor.has_value()) {
-        json[jsonLabels.linearGuessForceScaleFactor] = linearGuessForceScaleFactor.value();
+        json[GET_VARIABLE_NAME(linearGuessForceScaleFactor)] = linearGuessForceScaleFactor.value();
     }
-
+    if (viscousDampingFactor.has_value()) {
+        json[GET_VARIABLE_NAME(viscousDampingFactor)] = viscousDampingFactor.value();
+    }
+    json[GET_VARIABLE_NAME(useKineticDamping)] = useKineticDamping;
     //    if (intermediateResultsSaveStep.has_value()) {
     //        json[GET_VARIABLE_NAME(intermediateResultsSaveStep)] = intermediateResultsSaveStep.value();
     //    }
 
-    if (saveIntermediateBestStates.has_value()) {
-        json[GET_VARIABLE_NAME(saveIntermediateBestStates)] = saveIntermediateBestStates.value()
-                                                                  ? "true"
-                                                                  : "false";
-    }
+    //    if (saveIntermediateBestStates.has_value()) {
+    //        json[GET_VARIABLE_NAME(saveIntermediateBestStates)] = saveIntermediateBestStates.value()
+    //                                                                  ? "true"
+    //                                                                  : "false";
+    //    }
 
-    const std::string jsonFilename = "drmSettings.json";
-    std::ofstream jsonFile(std::filesystem::path(folderPath).append(jsonFilename));
+    std::filesystem::create_directories(jsonFilePath.parent_path());
+    std::ofstream jsonFile(jsonFilePath);
+    std::cout << "Saving DRM settings to:" << jsonFilePath << std::endl;
     jsonFile << json;
+    jsonFile.close();
 }
 
 bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePath)
@@ -2807,41 +2805,40 @@ bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePat
     std::ifstream ifs(jsonFilePath.string());
     ifs >> json;
 
-    if (json.contains(jsonLabels.shouldDraw)) {
-        shouldDraw = json.at(jsonLabels.shouldDraw) == "true" ? true : false;
+    if (json.contains(GET_VARIABLE_NAME(shouldDraw))) {
+        shouldDraw = json.at(GET_VARIABLE_NAME(shouldDraw)) == "true" ? true : false;
     }
-    if (json.contains(jsonLabels.beVerbose)) {
-        beVerbose = json.at(jsonLabels.beVerbose) == "true" ? true : false;
+    if (json.contains(GET_VARIABLE_NAME(beVerbose))) {
+        beVerbose = json.at(GET_VARIABLE_NAME(beVerbose)) == "true" ? true : false;
     }
-    if (json.contains(jsonLabels.shouldCreatePlots)) {
-        shouldCreatePlots = json.at(jsonLabels.shouldCreatePlots) == "true" ? true : false;
+    if (json.contains(GET_VARIABLE_NAME(shouldCreatePlots))) {
+        shouldCreatePlots = json.at(GET_VARIABLE_NAME(shouldCreatePlots)) == "true" ? true : false;
     }
-    if (json.contains(jsonLabels.Dtini)) {
-        Dtini = json.at(jsonLabels.Dtini);
+    if (json.contains(GET_VARIABLE_NAME(Dtini))) {
+        Dtini = json.at(GET_VARIABLE_NAME(Dtini));
     }
-    if (json.contains(jsonLabels.xi)) {
-        xi = json.at(jsonLabels.xi);
+    if (json.contains(GET_VARIABLE_NAME(xi))) {
+        xi = json.at(GET_VARIABLE_NAME(xi));
     }
-    if (json.contains(jsonLabels.maxDRMIterations)) {
-        maxDRMIterations = json.at(jsonLabels.maxDRMIterations);
+    if (json.contains(GET_VARIABLE_NAME(maxDRMIterations))) {
+        maxDRMIterations = json.at(GET_VARIABLE_NAME(maxDRMIterations));
     }
-    if (json.contains(jsonLabels.debugModeStep)) {
-        debugModeStep = json.at(jsonLabels.debugModeStep);
+    if (json.contains(GET_VARIABLE_NAME(debugModeStep))) {
+        debugModeStep = json.at(GET_VARIABLE_NAME(debugModeStep));
     }
-    if (json.contains(jsonLabels.totalTranslationalKineticEnergyThreshold)) {
+    if (json.contains(GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold))) {
         totalTranslationalKineticEnergyThreshold = json.at(
-            jsonLabels.totalTranslationalKineticEnergyThreshold);
+            GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold));
     }
-    if (json.contains(jsonLabels.gamma)) {
-        gamma = json.at(jsonLabels.gamma);
+    if (json.contains(GET_VARIABLE_NAME(gamma))) {
+        gamma = json.at(GET_VARIABLE_NAME(gamma));
     }
-    if (json.contains(jsonLabels.averageResidualForcesCriterionThreshold)) {
+    if (json.contains(GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold))) {
         averageResidualForcesCriterionThreshold = json.at(
-            jsonLabels.averageResidualForcesCriterionThreshold);
+            GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold));
     }
-
-    if (json.contains(jsonLabels.linearGuessForceScaleFactor)) {
-        linearGuessForceScaleFactor = json.at(jsonLabels.linearGuessForceScaleFactor);
+    if (json.contains(GET_VARIABLE_NAME(linearGuessForceScaleFactor))) {
+        linearGuessForceScaleFactor = json.at(GET_VARIABLE_NAME(linearGuessForceScaleFactor));
     }
 
     //    if (json.contains(GET_VARIABLE_NAME(intermediateResultsSaveStep))) {
@@ -2854,6 +2851,7 @@ bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePat
                                          ? true
                                          : false;
     }
+    //    std::cout << json.dump() << std::endl;
 
     return true;
 }
diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp
index 4d7e052..660e29b 100755
--- a/drmsimulationmodel.hpp
+++ b/drmsimulationmodel.hpp
@@ -33,7 +33,6 @@ class DRMSimulationModel : public SimulationModel
 public:
     struct Settings
     {
-        //        bool isDebugMode{false};
         bool shouldDraw{false};
         bool beVerbose{false};
         bool shouldCreatePlots{false};
@@ -46,8 +45,8 @@ public:
         int gradualForcedDisplacementSteps{50};
         //        int desiredGradualExternalLoadsSteps{1};
         double gamma{0.8};
-        double totalResidualForcesNormThreshold{1e-3};
-        double totalExternalForcesNormPercentageTermination{1e-3};
+        double totalResidualForcesNormThreshold{1e-8};
+        double totalExternalForcesNormPercentageTermination{1e-5};
         std::optional<int> maxDRMIterations;
         std::optional<int> debugModeStep;
         std::optional<double> totalTranslationalKineticEnergyThreshold;
@@ -55,11 +54,11 @@ public:
         std::optional<double> linearGuessForceScaleFactor;
         //        std::optional<int> intermediateResultsSaveStep;
         std::optional<bool> saveIntermediateBestStates;
-        Settings() {}
-        void save(const std::filesystem::path &folderPath) const;
-        bool load(const std::filesystem::path &filePath);
         std::optional<double> viscousDampingFactor;
         bool useKineticDamping{true};
+        Settings() {}
+        void save(const std::filesystem::path &jsonFilePath) const;
+        bool load(const std::filesystem::path &filePath);
         struct JsonLabels
         {
             const std::string shouldDraw{"shouldDraw"};
@@ -67,14 +66,16 @@ public:
             const std::string shouldCreatePlots{"shouldCreatePlots"};
             const std::string Dtini{"DtIni"};
             const std::string xi{"xi"};
+            const std::string gamma{"gamma"};
+            const std::string totalResidualForcesNormThreshold;
             const std::string maxDRMIterations{"maxIterations"};
             const std::string debugModeStep{"debugModeStep"};
             const std::string totalTranslationalKineticEnergyThreshold{
                 "totalTranslationaKineticEnergyThreshold"};
-            const std::string gamma{"gamma"};
             const std::string averageResidualForcesCriterionThreshold{
                 "averageResidualForcesThreshold"};
             const std::string linearGuessForceScaleFactor{"linearGuessForceScaleFactor"};
+            const std::string viscousDampingFactor{"viscousDampingFactor"};
         };
         static JsonLabels jsonLabels;
     };
diff --git a/edgemesh.cpp b/edgemesh.cpp
index 3fa7665..f479978 100755
--- a/edgemesh.cpp
+++ b/edgemesh.cpp
@@ -25,13 +25,13 @@ Eigen::MatrixX3d VCGEdgeMesh::getEigenEdgeNormals() const {
   return eigenEdgeNormals;
 }
 
-bool VCGEdgeMesh::save(const std::string &plyFilename)
+bool VCGEdgeMesh::save(const std::filesystem::__cxx11::path &meshFilePath)
 {
-    std::string filename = plyFilename;
+    std::string filename = meshFilePath;
     if (filename.empty()) {
         filename = std::filesystem::current_path().append(getLabel() + ".ply").string();
-    } else if (std::filesystem::is_directory(std::filesystem::path(plyFilename))) {
-        filename = std::filesystem::path(plyFilename).append(getLabel() + ".ply").string();
+    } else if (std::filesystem::is_directory(std::filesystem::path(meshFilePath))) {
+        filename = std::filesystem::path(meshFilePath).append(getLabel() + ".ply").string();
     }
     assert(std::filesystem::path(filename).extension().string() == ".ply");
     unsigned int mask = 0;
@@ -197,11 +197,11 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth,
   return true;
 }
 
-bool VCGEdgeMesh::load(const std::string &plyFilename)
+bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath)
 {
-    std::string usedPath = plyFilename;
-    if (std::filesystem::path(plyFilename).is_relative()) {
-        usedPath = std::filesystem::absolute(plyFilename).string();
+    std::string usedPath = meshFilePath;
+    if (std::filesystem::path(meshFilePath).is_relative()) {
+        usedPath = std::filesystem::absolute(meshFilePath).string();
     }
     assert(std::filesystem::exists(usedPath));
     Clear();
@@ -217,15 +217,15 @@ bool VCGEdgeMesh::load(const std::string &plyFilename)
     getEdges(eigenEdges);
     getVertices(eigenVertices);
     vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
-    label = std::filesystem::path(plyFilename).stem().string();
+    label = std::filesystem::path(meshFilePath).stem().string();
 
     const bool printDebugInfo = false;
     if (printDebugInfo) {
-        std::cout << plyFilename << " was loaded successfuly." << std::endl;
+        std::cout << meshFilePath << " was loaded successfuly." << std::endl;
         std::cout << "Mesh has " << EN() << " edges." << std::endl;
     }
 
-    label = std::filesystem::path(plyFilename).stem().string();
+    label = std::filesystem::path(meshFilePath).stem().string();
     return true;
 }
 
diff --git a/edgemesh.hpp b/edgemesh.hpp
index 8c5d8d2..ebd7ed8 100755
--- a/edgemesh.hpp
+++ b/edgemesh.hpp
@@ -61,8 +61,9 @@ public:
     void set(const std::vector<double> &vertexPositions, const std::vector<int> &edges);
 
     void removeDuplicateVertices();
-    void deleteDanglingVertices();
-    void deleteDanglingVertices(vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu);
+    virtual void deleteDanglingVertices();
+    virtual void deleteDanglingVertices(
+        vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu);
 
     void getEdges(Eigen::MatrixX3d &edgeStartingPoints, Eigen::MatrixX3d &edgeEndingPoints) const;
 
@@ -72,8 +73,8 @@ public:
 
     //    bool loadUsingNanoply(const std::string &plyFilename);
 
-    bool load(const std::string &plyFilename) override;
-    bool save(const std::string &plyFilename = std::string()) override;
+    bool load(const std::filesystem::path &meshFilePath) override;
+    bool save(const std::filesystem::path &meshFilePath = std::filesystem::path()) override;
 
     bool createSpanGrid(const size_t squareGridDimension);
     bool createSpanGrid(const size_t desiredWidth, const size_t desiredHeight);
diff --git a/hexagonremesher.hpp b/hexagonremesher.hpp
index 6d76967..4c2a6c2 100644
--- a/hexagonremesher.hpp
+++ b/hexagonremesher.hpp
@@ -11,7 +11,7 @@
 
 namespace PolygonalRemeshing {
 //std::shared_ptr<VCGPolyMesh>
-std::shared_ptr<VCGPolyMesh> remeshWithPolygons(VCGTriMesh &surface)
+inline std::shared_ptr<VCGPolyMesh> remeshWithPolygons(VCGTriMesh &surface)
 {
     //    const std::string tileIntoFilePath(
     //        "/home/iason/Coding/build/PatternTillingReducedModel/RelWithDebInfo/plane.ply");
diff --git a/mesh.hpp b/mesh.hpp
index 7ff46eb..31d0ab8 100755
--- a/mesh.hpp
+++ b/mesh.hpp
@@ -1,7 +1,7 @@
 #ifndef MESH_HPP
 #define MESH_HPP
 
-//#include <Eigen/Core>
+#include <filesystem>
 #include <string>
 
 class Mesh
@@ -11,8 +11,8 @@ protected:
 
 public:
   virtual ~Mesh() = default;
-  virtual bool load(const std::string &filePath) { return false; }
-  virtual bool save(const std::string &filePath) { return false; }
+  virtual bool load(const std::filesystem::path &meshFilePath) { return false; }
+  virtual bool save(const std::filesystem::path &meshFilePath) { return false; }
 
   std::string getLabel() const;
   void setLabel(const std::string &newLabel);
diff --git a/polymesh.hpp b/polymesh.hpp
index aa50d47..60ef843 100755
--- a/polymesh.hpp
+++ b/polymesh.hpp
@@ -5,6 +5,7 @@
 #include <filesystem>
 #include <wrap/io_trimesh/export_obj.h>
 #include <wrap/io_trimesh/export_ply.h>
+#include <wrap/io_trimesh/import.h>
 
 #ifdef POLYSCOPE_DEFINED
 #include <polyscope/surface_mesh.h>
@@ -64,14 +65,25 @@ class VCGPolyMesh
       public Mesh
 {
 public:
-    //  virtual bool load(const std::string &filename) override {
-    //    int mask;
-    //    vcg::tri::io::Importer<VCGPolyMesh>::LoadMask(filename.c_str(), mask);
-    //    int error = vcg::tri::io::Importer<VCGPolyMesh>::Open(
-    //        *this, filename.c_str(), mask);
-    //    if (error != 0) {
-    //      return false;
-    //    }
+    virtual bool load(const std::filesystem::__cxx11::path &meshFilePath) override
+    {
+        int mask;
+        vcg::tri::io::Importer<VCGPolyMesh>::LoadMask(meshFilePath.c_str(), mask);
+        int error = vcg::tri::io::Importer<VCGPolyMesh>::Open(*this, meshFilePath.c_str(), mask);
+        if (error != 0) {
+            std::cerr << "Could not load polygonal mesh:" << meshFilePath << std::endl;
+            return false;
+        }
+
+        vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
+        vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
+        vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
+        vcg::tri::UpdateNormal<VCGPolyMesh>::PerVertexNormalized(*this);
+        vcg::tri::Clean<VCGPolyMesh>::RemoveUnreferencedVertex(*this);
+        //finally remove valence 1 vertices on the border
+        //        vcg::PolygonalAlgorithm<PolyMeshType>::RemoveValence2Vertices(dual);
+        return true;
+    }
     //    //    vcg::tri::io::ImporterOBJ<VCGPolyMesh>::Open();
     //    //    unsigned int mask = 0;
     //    //    mask |= nanoply::NanoPlyWrapper<VCGPolyMesh>::IO_VERTCOORD;
@@ -84,11 +96,6 @@ public:
     //    //            0) {
     //    //      std::cout << "Could not load tri mesh" << std::endl;
     //    //    }
-    //    vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
-    //    vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
-    //    //    vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
-    //    vcg::tri::UpdateNormal<VCGPolyMesh>::PerVertexNormalized(*this);
-    //    return true;
     //  }
     Eigen::MatrixX3d getVertices() const
     {
@@ -113,6 +120,34 @@ public:
       return faces;
   }
 
+  //  bool load(const std::filesystem::__cxx11::path &meshFilePath)
+  //  {
+  //      const std::string extension = ".ply";
+  //      std::filesystem::path filePath = meshFilePath;
+  //      assert(std::filesystem::path(filePath).extension().string() == extension);
+  //      unsigned int mask = 0;
+  //      mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
+  //      mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
+  //      mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
+  //      mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
+  //      if (vcg::tri::io::Importer<VCGPolyMesh>::Open(*this, filePath.c_str()) != 0) {
+  //          return false;
+  //      }
+  //      label = meshFilePath.filename();
+  //      return true;
+  //  }
+
+  bool save(const std::filesystem::__cxx11::path &meshFilePath = std::filesystem::path())
+  {
+      if (meshFilePath.extension() == ".obj") {
+          return saveOBJ(meshFilePath);
+      } else if (meshFilePath.extension() == ".ply") {
+          return savePLY(meshFilePath);
+      }
+
+      return false;
+  }
+
   bool saveOBJ(const std::filesystem::path &objFilePath = std::filesystem::path())
   {
       const std::string extension = ".obj";
diff --git a/reducedmodelevaluator.cpp b/reducedmodelevaluator.cpp
new file mode 100644
index 0000000..89ee988
--- /dev/null
+++ b/reducedmodelevaluator.cpp
@@ -0,0 +1,317 @@
+#include "reducedmodelevaluator.hpp"
+#include "hexagonremesher.hpp"
+#include "trianglepatterngeometry.hpp"
+#include <filesystem>
+
+using FullPatternVertexIndex = VertexIndex;
+using ReducedPatternVertexIndex = VertexIndex;
+
+ReducedModelEvaluator::ReducedModelEvaluator()
+{
+
+}
+
+std::vector<double> ReducedModelEvaluator::evaluateReducedModel(
+    ReducedModelOptimization::Results &optimizationResults)
+{
+    //    std::shared_ptr<VCGPolyMesh> pTileIntoSurface = std::make_shared<VCGPolyMesh>();
+    //    std::filesystem::path tileIntoSurfaceFilePath{
+    //        "/home/iason/Coding/Projects/Approximating shapes with flat "
+    //        "patterns/ReducedModelOptimization/TestSet/TileIntoSurface/plane_34Polygons.ply"};
+    //    assert(std::filesystem::exists(tileIntoSurfaceFilePath));
+    //    pTileIntoSurface->load(tileIntoSurfaceFilePath);
+
+    //Set required file paths
+    const std::filesystem::path tileInto_triMesh_filename
+        = "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/"
+          "instantMeshes_plane_34.ply";
+    const std::filesystem::path reducedPatternFilePath
+        = "/home/iason/Coding/Projects/Approximating shapes with flat "
+          "patterns/ReducedModelOptimization/TestSet/ReducedPatterns/single_reduced.ply";
+    const std::filesystem::path intermediateResultsDirectoryPath
+        = "/home/iason/Coding/build/ReducedModelOptimization/IntermediateResults";
+    //    const std::filesystem::path drmSettingsFilePath
+    //        = "/home/iason/Coding/Projects/Approximating shapes with flat "
+    //          "patterns/ReducedModelOptimization/DefaultSettings/DRMSettings/"
+    //          "defaultDRMSimulationSettings.json";
+    DRMSimulationModel::Settings drmSimulationSettings;
+    //    drmSimulationSettings.load(drmSettingsFilePath);
+    drmSimulationSettings.linearGuessForceScaleFactor = 1;
+    drmSimulationSettings.debugModeStep = 1000;
+    drmSimulationSettings.beVerbose = true;
+    constexpr bool shouldRerunFullPatternSimulation = false;
+    const std::vector<std::string> scenariosTestSetLabels{"22Hex_randomBending0",
+                                                          "22Hex_randomBending1",
+                                                          "22Hex_randomBending2",
+                                                          "22Hex_randomBending3",
+                                                          "22Hex_randomBending4",
+                                                          "22Hex_randomBending5",
+                                                          "22Hex_randomBending6",
+                                                          "22Hex_randomBending7",
+                                                          "22Hex_randomBending8",
+                                                          "22Hex_randomBending9",
+                                                          "22Hex_randomBending10",
+                                                          "22Hex_randomBending11",
+                                                          "22Hex_randomBending12",
+                                                          "22Hex_randomBending13",
+                                                          "22Hex_randomBending14",
+                                                          "22Hex_randomBending15",
+                                                          "22Hex_randomBending16",
+                                                          "22Hex_randomBending17",
+                                                          "22Hex_randomBending18",
+                                                          "22Hex_randomBending19",
+                                                          "22Hex_randomBending20",
+                                                          "22Hex_bending_0.005N",
+                                                          "22Hex_bending_0.01N",
+                                                          "22Hex_bending_0.03N",
+                                                          "22Hex_bending_0.05N",
+                                                          "22Hex_pullOppositeVerts_0.05N",
+                                                          "22Hex_pullOppositeVerts_0.1N",
+                                                          "22Hex_pullOppositeVerts_0.3N",
+                                                          "22Hex_shear_2N",
+                                                          "22Hex_shear_5N",
+                                                          "22Hex_axial_10N",
+                                                          "22Hex_axial_20N",
+                                                          "22Hex_cylinder_0.05N",
+                                                          "22Hex_cylinder_0.1N",
+                                                          "22Hex_s_0.05N",
+                                                          "22Hex_s_0.1N"};
+
+    //Load surface
+    std::shared_ptr<VCGPolyMesh> pTileIntoSurface = [&]() {
+        VCGTriMesh tileInto_triMesh;
+        const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filename);
+        assert(surfaceLoadSuccessfull);
+        return PolygonalRemeshing::remeshWithPolygons(tileInto_triMesh);
+    }();
+    const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResults.baseTriangle.cP(0),
+                                                             (optimizationResults.baseTriangle.cP(1)
+                                                              + optimizationResults.baseTriangle.cP(
+                                                                  2))
+                                                                 / 2);
+    pTileIntoSurface->moveToCenter();
+    const double scaleFactor = optimizedBaseTriangleHeight
+                               / pTileIntoSurface->getAverageFaceRadius();
+    vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
+    //    tileIntoSurface.registerForDrawing();
+    //    polyscope::show();
+
+    //Tile full pattern into surface
+    std::vector<PatternGeometry> fullPatterns(1);
+    fullPatterns[0].copy(optimizationResults.baseTriangleFullPattern);
+    ////        Base triangle pattern might contain dangling vertices.Remove those
+    fullPatterns[0].interfaceNodeIndex = 3;
+    fullPatterns[0].deleteDanglingVertices();
+    std::vector<int> perSurfaceFacePatternIndices(pTileIntoSurface->FN(), 0);
+    std::vector<std::vector<size_t>> perPatternIndexTiledFullPatternEdgeIndices;
+    std::vector<size_t> tileIntoEdgeToTiledFullVi;
+    std::shared_ptr<PatternGeometry> pTiledFullPattern
+        = PatternGeometry::tilePattern(fullPatterns,
+                                       {},
+                                       *pTileIntoSurface,
+                                       perSurfaceFacePatternIndices,
+                                       tileIntoEdgeToTiledFullVi,
+                                       perPatternIndexTiledFullPatternEdgeIndices);
+    pTiledFullPattern->setLabel("Tiled_full_patterns");
+    //    pTiledFullPattern->registerForDrawing();
+    //Tile reduced pattern into surface
+    PatternGeometry reducedPattern;
+    reducedPattern.load(reducedPatternFilePath);
+    reducedPattern.deleteDanglingVertices();
+    assert(reducedPattern.interfaceNodeIndex == 1);
+    std::vector<PatternGeometry> reducedPatterns(1);
+    reducedPatterns[0].copy(reducedPattern);
+    const auto reducedPatternBaseTriangle = reducedPattern.computeBaseTriangle();
+    ReducedModelOptimization::Results::applyOptimizationResults_innerHexagon(
+        optimizationResults, reducedPatternBaseTriangle, reducedPatterns[0]);
+
+    std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
+    std::vector<size_t> tileIntoEdgeToTiledReducedVi;
+    std::shared_ptr<PatternGeometry> pTiledReducedPattern
+        = PatternGeometry::tilePattern(reducedPatterns,
+                                       {0},
+                                       *pTileIntoSurface,
+                                       perSurfaceFacePatternIndices,
+                                       tileIntoEdgeToTiledReducedVi,
+                                       perPatternIndexTiledReducedPatternEdgeIndices);
+    pTiledReducedPattern->setLabel("Tiled_reduced_patterns");
+    //        pTiledReducedPattern->registerForDrawing();
+
+    std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
+        fullToReducedViMap; //of only the common vertices
+    std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+        reducedToFullViMap; //of only the common vertices
+    for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) {
+        fullToReducedViMap[tileIntoEdgeToTiledFullVi[ei]] = tileIntoEdgeToTiledReducedVi[ei];
+    }
+    constructInverseMap(fullToReducedViMap, reducedToFullViMap);
+
+    std::vector<size_t> tilledFullPatternInterfaceVi;
+    tilledFullPatternInterfaceVi.clear();
+    tilledFullPatternInterfaceVi.resize(fullToReducedViMap.size());
+    size_t viIndex = 0;
+    for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
+        tilledFullPatternInterfaceVi[viIndex++] = fullToReducedPair.first;
+    }
+
+    //Create simulation meshes
+    ////Tessellated full pattern simulation mesh
+    std::shared_ptr<SimulationMesh> pTiledFullPattern_simulationMesh;
+    pTiledFullPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledFullPattern);
+    //NOTE: Those should be derived from the optimization results instead of hardcoding them
+    pTiledFullPattern_simulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002});
+    if (optimizationResults.fullPatternYoungsModulus == 0) {
+        std::cerr << "Full pattern's young modulus not found." << std::endl;
+        std::terminate();
+    }
+    pTiledFullPattern_simulationMesh->setBeamMaterial(0.3,
+                                                      optimizationResults.fullPatternYoungsModulus);
+    pTiledFullPattern_simulationMesh->reset();
+
+    ////Tessellated reduced pattern simulation mesh
+    std::shared_ptr<SimulationMesh> pTiledReducedPattern_simulationMesh;
+    pTiledReducedPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledReducedPattern);
+    const std::vector<size_t> &tiledPatternElementIndicesForReducedPattern
+        = perPatternIndexTiledReducedPatternEdgeIndices[0];
+    ReducedModelOptimization::Results::applyOptimizationResults_elements(
+        optimizationResults, pTiledReducedPattern_simulationMesh);
+    pTiledReducedPattern_simulationMesh->reset();
+    std::vector<double> distances_drm2reduced(scenariosTestSetLabels.size(), 0);
+    std::vector<double> distances_normalizedDrm2reduced(scenariosTestSetLabels.size(), 0);
+    for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
+        const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
+        const std::filesystem::path scenariosDirectoryPath
+            = "/home/iason/Coding/build/PatternTillingReducedModel/Scenarios/";
+        const std::filesystem::path tiledReducedPatternJobFilePath
+            = std::filesystem::path(scenariosDirectoryPath)
+                  .append(jobLabel)
+                  .append("SimulationJobs")
+                  .append("Reduced")
+                  .append(SimulationJob::jsonDefaultFileName);
+        //set jobs
+        std::shared_ptr<SimulationJob> pJob_tiledReducedPattern;
+        pJob_tiledReducedPattern = std::make_shared<SimulationJob>(SimulationJob());
+        pJob_tiledReducedPattern->load(tiledReducedPatternJobFilePath, false);
+        pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
+        std::shared_ptr<SimulationJob> pJob_tiledFullPattern;
+        pJob_tiledFullPattern = std::make_shared<SimulationJob>(SimulationJob());
+        pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
+        pJob_tiledReducedPattern->remap(reducedToFullViMap, *pJob_tiledFullPattern);
+        //        pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern->getLabel());
+        //        pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel());
+        //    polyscope::show();
+        //Save reduced job
+        const std::filesystem::path tesellatedResultsFolderPath
+            = std::filesystem::path(intermediateResultsDirectoryPath).append("TessellatedResults");
+        const std::filesystem::path surfaceFolderPath = std::filesystem::path(
+                                                            tesellatedResultsFolderPath)
+                                                            .append(pTileIntoSurface->getLabel());
+        const std::string scenarioLabel = pJob_tiledFullPattern->getLabel();
+        const std::filesystem::path scenarioDirectoryPath = std::filesystem::path(surfaceFolderPath)
+                                                                .append(scenarioLabel);
+        const std::filesystem::path reducedJobDirectoryPath
+            = std::filesystem::path(scenarioDirectoryPath).append("ReducedJob");
+        std::filesystem::create_directories(reducedJobDirectoryPath);
+        pJob_tiledReducedPattern->save(reducedJobDirectoryPath);
+        //Run scenario
+
+        ////Full
+
+        const std::string patternLabel = [&]() {
+            const std::string patternLabel = optimizationResults.baseTriangleFullPattern.getLabel();
+            const int numberOfOccurences = std::count_if(patternLabel.begin(),
+                                                         patternLabel.end(),
+                                                         [](char c) { return c == '#'; });
+            if (numberOfOccurences == 0) {
+                return std::to_string(optimizationResults.baseTriangleFullPattern.EN()) + "#"
+                       + optimizationResults.baseTriangleFullPattern.getLabel();
+            } else if (numberOfOccurences == 1) {
+                return optimizationResults.baseTriangleFullPattern.getLabel();
+            }
+        }();
+        const auto fullResultsFolderPath
+            = std::filesystem::path(scenarioDirectoryPath).append(patternLabel).append("Results");
+        if (shouldRerunFullPatternSimulation && std::filesystem::exists(fullResultsFolderPath)) {
+            std::filesystem::remove_all(fullResultsFolderPath);
+        }
+
+        const std::filesystem::path fullPatternJobFolderPath = std::filesystem::path(
+                                                                   scenarioDirectoryPath)
+                                                                   .append(patternLabel)
+                                                                   .append("SimulationJob");
+        SimulationResults simulationResults_tiledFullPattern_drm;
+        if (std::filesystem::exists(fullResultsFolderPath)) {
+            //Load full pattern results
+            assert(std::filesystem::exists(fullPatternJobFolderPath));
+            simulationResults_tiledFullPattern_drm.load(fullResultsFolderPath,
+                                                        fullPatternJobFolderPath);
+            simulationResults_tiledFullPattern_drm.converged = true;
+        } else {
+            //Full
+            std::cout << "Executing:" << jobLabel << std::endl;
+            DRMSimulationModel drmSimulationModel;
+            simulationResults_tiledFullPattern_drm
+                = drmSimulationModel.executeSimulation(pJob_tiledFullPattern, drmSimulationSettings);
+            simulationResults_tiledFullPattern_drm.setLabelPrefix("DRM");
+        }
+        std::filesystem::create_directories(fullResultsFolderPath);
+        simulationResults_tiledFullPattern_drm.save(
+            std::filesystem::path(scenarioDirectoryPath).append(patternLabel));
+        if (!simulationResults_tiledFullPattern_drm.converged) {
+            std::cerr << "Full pattern simulation failed." << std::endl;
+        }
+
+        LinearSimulationModel linearSimulationModel;
+        SimulationResults simulationResults_tiledReducedPattern
+            = linearSimulationModel.executeSimulation(pJob_tiledReducedPattern);
+
+        //        simulationResults_tiledReducedPattern.registerForDrawing();
+        //        simulationResults_tiledFullPattern_drm.registerForDrawing();
+        //        polyscope::show();
+        //measure distance
+        const double distance_fullDrmToReduced
+            = simulationResults_tiledFullPattern_drm
+                  .computeDistance(simulationResults_tiledReducedPattern, fullToReducedViMap);
+        double distance_fullSumOfAllVerts = 0;
+        for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
+            distance_fullSumOfAllVerts += simulationResults_tiledFullPattern_drm
+                                              .displacements[fullToReducedPair.first]
+                                              .getTranslation()
+                                              .norm();
+        }
+        const double distance_normalizedFullDrmToReduced = distance_fullDrmToReduced
+                                                           / distance_fullSumOfAllVerts;
+        distances_drm2reduced[jobIndex] = distance_fullDrmToReduced;
+        distances_normalizedDrm2reduced[jobIndex] = distance_normalizedFullDrmToReduced;
+    }
+    //#ifndef POLYSCOPE_DEFINED
+    //    return distances_drm2reduced;
+    //#else
+
+    //report distance
+    //        csvFile csv_results("", flse);
+    csvFile csv_results({}, false);
+    //    csvFile csv_results(std::filesystem::path(dirPath_thisOptimization)
+    //                            .append("results.csv")
+    //                            .string(),
+    //                        false);
+    csv_results /*<< "Job Label"*/
+        << "drm2Reduced"
+        << "norm_drm2Reduced";
+    csv_results << endrow;
+    double sumOfNormalizedFull2Reduced = 0;
+    for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
+        const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
+        const double &distance_fullDrmToReduced = distances_drm2reduced[jobIndex];
+        const double &distance_normalizedFullDrmToReduced = distances_normalizedDrm2reduced[jobIndex];
+        csv_results /*<< jobLabel*/ << distance_fullDrmToReduced
+                                    << distance_normalizedFullDrmToReduced;
+        csv_results << endrow;
+        sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
+    }
+
+    std::cout << "Average normalized error per scenario:"
+              << sumOfNormalizedFull2Reduced / scenariosTestSetLabels.size() << std::endl;
+    return distances_normalizedDrm2reduced;
+    //#endif
+}
diff --git a/reducedmodelevaluator.hpp b/reducedmodelevaluator.hpp
new file mode 100644
index 0000000..dd06918
--- /dev/null
+++ b/reducedmodelevaluator.hpp
@@ -0,0 +1,14 @@
+#ifndef REDUCEDMODELEVALUATOR_HPP
+#define REDUCEDMODELEVALUATOR_HPP
+
+#include "reducedmodeloptimizer_structs.hpp"
+
+class ReducedModelEvaluator
+{
+public:
+    ReducedModelEvaluator();
+    static std::vector<double> evaluateReducedModel(
+        ReducedModelOptimization::Results &optimizationResults);
+};
+
+#endif // REDUCEDMODELEVALUATOR_HPP
diff --git a/reducedmodeloptimizer_structs.hpp b/reducedmodeloptimizer_structs.hpp
index e79f126..6b2c5ea 100644
--- a/reducedmodeloptimizer_structs.hpp
+++ b/reducedmodeloptimizer_structs.hpp
@@ -7,7 +7,7 @@
 #include "unordered_map"
 #include <string>
 
-namespace ReducedPatternOptimization {
+namespace ReducedModelOptimization {
 
 enum BaseSimulationScenario {
     Axial,
@@ -107,21 +107,29 @@ inline int getParameterIndex(const std::string &s)
 struct Settings
 {
     inline static std::string defaultFilename{"OptimizationSettings.json"};
-    std::vector<std::vector<OptimizationParameterIndex>> optimizationStrategy;
+    std::vector<std::vector<OptimizationParameterIndex>> optimizationStrategy = {
+        {E, A, I2, I3, J, R, Theta}};
+    //        {A, I2, I3, J, R, Theta}};
     std::optional<std::vector<double>>
         optimizationVariablesGroupsWeights; //TODO:should be removed in the future if not a splitting strategy is used for the optimization
     enum NormalizationStrategy { NonNormalized, Epsilon };
     inline static std::vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
     NormalizationStrategy normalizationStrategy{Epsilon};
-    std::array<xRange, NumberOfOptimizationVariables> variablesRanges;
+    std::array<xRange, NumberOfOptimizationVariables> variablesRanges{xRange{"E", 0.001, 1000},
+                                                                      xRange{"A", 0.001, 1000},
+                                                                      xRange{"I2", 0.001, 1000},
+                                                                      xRange{"I3", 0.001, 1000},
+                                                                      xRange{"J", 0.001, 1000},
+                                                                      xRange{"R", 0.05, 0.95},
+                                                                      xRange{"Theta", -30, 30}};
     int numberOfFunctionCalls{100000};
     double solverAccuracy{1e-3};
-    double translationNormalizationEpsilon{1e-4};
-    double rotationNormalizationEpsilon{vcg::math::ToRad(3.0)};
+    double translationEpsilon{4e-3};
+    double angularDistanceEpsilon{vcg::math::ToRad(0.0)};
     struct ObjectiveWeights
     {
-        double translational{1};
-        double rotational{1};
+        double translational{1.2};
+        double rotational{0.8};
         bool operator==(const ObjectiveWeights &other) const
         {
             return this->translational == other.translational
@@ -155,15 +163,15 @@ struct Settings
 
     void setDefault()
     {
-        ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000};
-        ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000};
-        ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000};
-        ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000};
-        ReducedPatternOptimization::xRange beamJ{"J", 0.001, 1000};
-        ReducedPatternOptimization::xRange innerHexagonSize{"R", 0.05, 0.95};
-        ReducedPatternOptimization::xRange innerHexagonAngle{"Theta", -30.0, 30.0};
-        variablesRanges = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle};
-        numberOfFunctionCalls = 100000;
+        //        ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000};
+        //        ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000};
+        //        ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000};
+        //        ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000};
+        //        ReducedPatternOptimization::xRange beamJ{"J", 0.001, 1000};
+        //        ReducedPatternOptimization::xRange innerHexagonSize{"R", 0.05, 0.95};
+        //        ReducedPatternOptimization::xRange innerHexagonAngle{"Theta", -30.0, 30.0};
+        //        variablesRanges = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle};
+        //        numberOfFunctionCalls = 100000;
 
         enum OptimizationParameterComparisonScenarioIndex {
             AllVar,
@@ -179,12 +187,12 @@ struct Settings
             NumberOfScenarios
         };
         const std::vector<
-            std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>>
+            std::vector<std::vector<ReducedModelOptimization::OptimizationParameterIndex>>>
             optimizationParameters = [&]() {
                 std::vector<
-                    std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>>
+                    std::vector<std::vector<ReducedModelOptimization::OptimizationParameterIndex>>>
                     optimizationParameters(NumberOfScenarios);
-                using namespace ReducedPatternOptimization;
+                using namespace ReducedModelOptimization;
                 optimizationParameters[AllVar] = {{E, A, I2, I3, J, R, Theta}};
                 optimizationParameters[GeoYM] = {{R, Theta, E}};
                 optimizationParameters[MatGeo] = {{A, I2, I3, J, R, Theta}};
@@ -198,43 +206,28 @@ struct Settings
                 return optimizationParameters;
             }();
 
-        constexpr OptimizationParameterComparisonScenarioIndex scenario = MatGeo;
+        constexpr OptimizationParameterComparisonScenarioIndex scenario = AllVar;
         optimizationStrategy = optimizationParameters[scenario];
         if (scenario == YMGeo_Mat) {
             optimizationVariablesGroupsWeights = {0.15, 0.85};
         }
 
-        normalizationStrategy = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
-        //#ifdef POLYSCOPE_DEFINED
-        //#else
-        translationNormalizationEpsilon = 4e-3;
-        rotationNormalizationEpsilon = vcg::math::ToRad(3.0);
+        //        normalizationStrategy = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
+        //        translationNormalizationEpsilon = 0;
+        //        rotationNormalizationEpsilon = vcg::math::ToRad(0.0);
+        //        solverAccuracy = 1e-2;
 
-        //     translationNormalizationParameter = 0;
-        //     rotationNormalizationParameter = 0;
-        //     translationNormalizationParameter = 1e-3;
-        //     rotationNormalizationParameter = vcg::math::ToRad(3.0);
-        //     translationNormalizationParameter = 5e-2;
-        //     rotationNormalizationParameter = vcg::math::ToRad(3.0);
-        //#endif
-        solverAccuracy = 1e-2;
-        //TODO refactor that
-        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].translational = 1.5;
-        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].rotational
-            = 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].translational;
-        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].translational = 1.5;
-        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].rotational
-            = 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].translational;
-        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].translational = 1.2;
-        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].rotational
-            = 2
-              - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].translational;
-        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].translational = 1.95;
-        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].rotational
-            = 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].translational;
-        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].translational = 1.2;
-        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].rotational
-            = 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].translational;
+        //        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].translational = 1.174011;
+        //        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].translational = 1.690460;
+        //        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].translational = 1.812089;
+        //        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].translational = 0.127392;
+        //        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].translational = 0.503043;
+
+        //        for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
+        //             baseScenario++) {
+        //            perBaseScenarioObjectiveWeights[baseScenario].rotational
+        //                = 2 - perBaseScenarioObjectiveWeights[baseScenario].translational;
+        //        }
     }
 
     void save(const std::filesystem::path &saveToPath)
@@ -242,10 +235,10 @@ struct Settings
         assert(std::filesystem::is_directory(saveToPath));
 
         nlohmann::json json;
-        json[JsonKeys::OptimizationStrategy] = optimizationStrategy;
+        json[GET_VARIABLE_NAME(optimizationStrategy)] = optimizationStrategy;
         if (optimizationVariablesGroupsWeights.has_value()) {
-            json[JsonKeys::OptimizationStrategyGroupWeights] = optimizationVariablesGroupsWeights
-                                                                   .value();
+            json[GET_VARIABLE_NAME(ptimizationStrategyGroupWeights)]
+                = optimizationVariablesGroupsWeights.value();
         }
         //write x ranges
         const std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
@@ -267,8 +260,8 @@ struct Settings
         //                = xRange.toString();
         //        }
 
-        json[JsonKeys::NumberOfFunctionCalls] = numberOfFunctionCalls;
-        json[JsonKeys::SolverAccuracy] = solverAccuracy;
+        json[GET_VARIABLE_NAME(numberOfFunctionCalls)] = numberOfFunctionCalls;
+        json[GET_VARIABLE_NAME(solverAccuracy)] = solverAccuracy;
         //Objective weights
         std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios> objectiveWeightsPairs;
         std::transform(perBaseScenarioObjectiveWeights.begin(),
@@ -279,11 +272,15 @@ struct Settings
                                                  objectiveWeights.rotational);
                        });
         json[JsonKeys::ObjectiveWeights] = objectiveWeightsPairs;
+        json[GET_VARIABLE_NAME(translationEpsilon)] = translationEpsilon;
+        json[GET_VARIABLE_NAME(angularDistanceEpsilon)] = vcg::math::ToDeg(
+            angularDistanceEpsilon);
 
         std::filesystem::path jsonFilePath(
             std::filesystem::path(saveToPath).append(defaultFilename));
         std::ofstream jsonFile(jsonFilePath.string());
         jsonFile << json;
+        jsonFile.close();
     }
 
     bool load(const std::filesystem::path &jsonFilepath)
@@ -298,14 +295,14 @@ struct Settings
         std::ifstream ifs(jsonFilepath);
         ifs >> json;
 
-        if (json.contains(JsonKeys::OptimizationStrategy)) {
+        if (json.contains(GET_VARIABLE_NAME(optimizationStrategy))) {
             optimizationStrategy
-                = std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>(
-                    (json.at(JsonKeys::OptimizationStrategy)));
+                = std::vector<std::vector<ReducedModelOptimization::OptimizationParameterIndex>>(
+                    (json.at(GET_VARIABLE_NAME(optimizationStrategy))));
         }
-        if (json.contains(JsonKeys::OptimizationStrategyGroupWeights)) {
+        if (json.contains(GET_VARIABLE_NAME(optimizationStrategyGroupWeights))) {
             optimizationVariablesGroupsWeights = std::vector<double>(
-                json[JsonKeys::OptimizationStrategyGroupWeights]);
+                json[GET_VARIABLE_NAME(optimizationStrategyGroupWeights)]);
         }
         //read x ranges
         if (json.contains(JsonKeys::OptimizationVariables)) {
@@ -328,8 +325,12 @@ struct Settings
             }
         }
 
-        numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls);
-        solverAccuracy = json.at(JsonKeys::SolverAccuracy);
+        if (json.contains(GET_VARIABLE_NAME(numberOfFunctionCalls))) {
+            numberOfFunctionCalls = json.at(GET_VARIABLE_NAME(numberOfFunctionCalls));
+        }
+        if (json.contains(GET_VARIABLE_NAME(solverAccuracy))) {
+            solverAccuracy = json.at(GET_VARIABLE_NAME(solverAccuracy));
+        }
         //Objective weights
         if (json.contains(JsonKeys::ObjectiveWeights)) {
             std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios>
@@ -342,6 +343,15 @@ struct Settings
                                                        objectiveWeightsPair.second};
                            });
         }
+        if (json.contains(GET_VARIABLE_NAME(translationalNormalizationEpsilon))) {
+            translationEpsilon
+                = json[GET_VARIABLE_NAME(translationalNormalizationEpsilon)];
+        }
+
+        if (json.contains(GET_VARIABLE_NAME(angularDistanceEpsilon))) {
+            angularDistanceEpsilon = vcg::math::ToRad(
+                static_cast<double>(json[GET_VARIABLE_NAME(angularDistanceEpsilon)]));
+        }
 
         //        perBaseScenarioObjectiveWeights = json.at(JsonKeys::ObjectiveWeights);
         //        objectiveWeights.translational = json.at(JsonKeys::ObjectiveWeights);
@@ -381,7 +391,8 @@ struct Settings
         }
         os << "Function Calls";
         os << "Solution Accuracy";
-        os << "Normalization strategy";
+        os << "Normalization trans epsilon";
+        os << "Normalization rot epsilon(deg)";
         os << JsonKeys::ObjectiveWeights;
         os << "Optimization parameters";
         //        os << std::endl;
@@ -397,8 +408,8 @@ struct Settings
         }
         os << numberOfFunctionCalls;
         os << solverAccuracy;
-        os << normalizationStrategyStrings[normalizationStrategy] + "_"
-                  + std::to_string(translationNormalizationEpsilon);
+        os << std::to_string(translationEpsilon);
+        os << std::to_string(vcg::math::ToDeg(angularDistanceEpsilon));
         std::string objectiveWeightsString;
         objectiveWeightsString += "{";
         for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
@@ -436,8 +447,8 @@ struct Settings
       return settings1.numberOfFunctionCalls == settings2.numberOfFunctionCalls
              && settings1.variablesRanges == settings2.variablesRanges
              && settings1.solverAccuracy == settings2.solverAccuracy && haveTheSameObjectiveWeights
-             && settings1.translationNormalizationEpsilon
-                    == settings2.translationNormalizationEpsilon;
+             && settings1.translationEpsilon
+                    == settings2.translationEpsilon;
   }
 
   inline void updateMeshWithOptimalVariables(const std::vector<double> &x, SimulationMesh &m)
@@ -504,7 +515,7 @@ struct Settings
       std::vector<std::pair<std::string, double>> optimalXNameValuePairs;
       std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
       std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
-      double fullPatternYoungsModulus;
+      double fullPatternYoungsModulus{0};
 
       PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern
       vcg::Triangle3<double> baseTriangle;
@@ -646,7 +657,7 @@ struct Settings
 
           settings.save(saveToPath);
 
-#ifdef POLYSCOPE_DEFINED12
+#ifdef POLYSCOPE_DEFINED
           std::cout << "Saved optimization results to:" << saveToPath << std::endl;
 #endif
       }
@@ -748,14 +759,16 @@ struct Settings
 
       template<typename MeshType>
       static void applyOptimizationResults_innerHexagon(
-          const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
+          const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
           const vcg::Triangle3<double> &patternBaseTriangle,
           MeshType &reducedPattern)
       {
           std::unordered_map<std::string, double>
               optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(),
                                 reducedPattern_optimizationResults.optimalXNameValuePairs.end());
-          assert(optimalXVariables.contains("R") && optimalXVariables.contains("Theta"));
+          assert(
+              (optimalXVariables.contains("R") && optimalXVariables.contains("Theta"))
+              || (optimalXVariables.contains("HexSize") && optimalXVariables.contains("HexAngle")));
           if (optimalXVariables.contains("HexSize")) {
               applyOptimizationResults_innerHexagon(optimalXVariables.at("HexSize"),
                                                     optimalXVariables.at("HexAngle"),
@@ -807,8 +820,8 @@ struct Settings
       }
 
       static void applyOptimizationResults_elements(
-          const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
-          const std::shared_ptr<SimulationMesh> &pTiledReducedPattern_simulationMesh)
+          const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
+          const std::shared_ptr<SimulationMesh> &pReducedPattern_simulationMesh)
       {
           assert(CrossSectionType::name == RectangularBeamDimensions::name);
           // Set reduced parameters fron the optimization results
@@ -817,55 +830,37 @@ struct Settings
                                 reducedPattern_optimizationResults.optimalXNameValuePairs.end());
 
           const std::string ALabel = "A";
-          if (optimalXVariables.contains(ALabel)) {
-              const double A = optimalXVariables.at(ALabel);
-              const double beamWidth = std::sqrt(A);
-              const double beamHeight = beamWidth;
-              CrossSectionType elementDimensions(beamWidth, beamHeight);
-              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
-                  e.setDimensions(elementDimensions);
-              }
-          }
+          assert(optimalXVariables.contains(ALabel));
+          const double A = optimalXVariables.at(ALabel);
+          const double beamWidth = std::sqrt(A);
+          const double beamHeight = beamWidth;
+          CrossSectionType elementDimensions(beamWidth, beamHeight);
 
           const double poissonsRatio = 0.3;
           const std::string ymLabel = "E";
-          if (optimalXVariables.contains(ymLabel)) {
-              const double E = optimalXVariables.at(ymLabel);
-              const ElementMaterial elementMaterial(poissonsRatio, E);
-              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
-                  e.setMaterial(elementMaterial);
-              }
-          }
+          assert(optimalXVariables.contains(ymLabel));
+          const double E = optimalXVariables.at(ymLabel);
+          const ElementMaterial elementMaterial(poissonsRatio, E);
 
           const std::string JLabel = "J";
-          if (optimalXVariables.contains(JLabel)) {
-              const double J = optimalXVariables.at(JLabel);
-              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
-                  e.dimensions.inertia.J = J;
-              }
-          }
+          assert(optimalXVariables.contains(JLabel));
+          const double J = optimalXVariables.at(JLabel);
 
           const std::string I2Label = "I2";
-          if (optimalXVariables.contains(I2Label)) {
-              const double I2 = optimalXVariables.at(I2Label);
-              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
-                  e.dimensions.inertia.I2 = I2;
-              }
-          }
+          assert(optimalXVariables.contains(I2Label));
+          const double I2 = optimalXVariables.at(I2Label);
 
           const std::string I3Label = "I3";
-          if (optimalXVariables.contains(I3Label)) {
-              const double I3 = optimalXVariables.at(I3Label);
-              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
-                  e.dimensions.inertia.I3 = I3;
-              }
+          assert(optimalXVariables.contains(I3Label));
+          const double I3 = optimalXVariables.at(I3Label);
+          for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+              Element &e = pReducedPattern_simulationMesh->elements[ei];
+              e.setDimensions(elementDimensions);
+              e.setMaterial(elementMaterial);
+              e.dimensions.inertia.J = J;
+              e.dimensions.inertia.I2 = I2;
+              e.dimensions.inertia.I3 = I3;
           }
-          pTiledReducedPattern_simulationMesh->reset();
       }
 
 #if POLYSCOPE_DEFINED
diff --git a/simulation_structs.hpp b/simulation_structs.hpp
index e30240b..7f4b908 100755
--- a/simulation_structs.hpp
+++ b/simulation_structs.hpp
@@ -113,25 +113,27 @@ class SimulationJob {
   } jsonLabels;
 
 public:
-  std::shared_ptr<SimulationMesh> pMesh;
-  std::string label{"empty_job"};
-  std::unordered_map<VertexIndex, std::unordered_set<int>> constrainedVertices;
-  std::unordered_map<VertexIndex, Vector6d> nodalExternalForces;
-  std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
-  SimulationJob(
-      const std::shared_ptr<SimulationMesh> &m, const std::string &label,
-      const std::unordered_map<VertexIndex, std::unordered_set<int>> &cv,
-      const std::unordered_map<VertexIndex, Vector6d> &ef = {},
-      const std::unordered_map<VertexIndex, Eigen::Vector3d> &fd = {})
-      : pMesh(m), label(label), constrainedVertices(cv),
-        nodalExternalForces(ef), nodalForcedDisplacements(fd) {}
-  SimulationJob() {}
-  SimulationJob(const std::string &jsonFilename) { load(jsonFilename); }
+    inline static std::string jsonDefaultFileName = "SimulationJob.json";
+    std::shared_ptr<SimulationMesh> pMesh;
+    std::string label{"empty_job"};
+    std::unordered_map<VertexIndex, std::unordered_set<int>> constrainedVertices;
+    std::unordered_map<VertexIndex, Vector6d> nodalExternalForces;
+    std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
+    SimulationJob(const std::shared_ptr<SimulationMesh> &m,
+                  const std::string &label,
+                  const std::unordered_map<VertexIndex, std::unordered_set<int>> &cv,
+                  const std::unordered_map<VertexIndex, Vector6d> &ef = {},
+                  const std::unordered_map<VertexIndex, Eigen::Vector3d> &fd = {})
+        : pMesh(m), label(label), constrainedVertices(cv), nodalExternalForces(ef),
+          nodalForcedDisplacements(fd)
+    {}
+    SimulationJob() {}
+    SimulationJob(const std::string &jsonFilename) { load(jsonFilename); }
 
-  bool isEmpty()
-  {
-      return constrainedVertices.empty() && nodalExternalForces.empty()
-             && nodalForcedDisplacements.empty() && pMesh == nullptr;
+    bool isEmpty()
+    {
+        return constrainedVertices.empty() && nodalExternalForces.empty()
+               && nodalForcedDisplacements.empty() && pMesh == nullptr;
   }
   void remap(const std::unordered_map<size_t, size_t> &sourceToDestinationViMap,
              SimulationJob &destination_simulationJob) const
@@ -489,6 +491,7 @@ struct SimulationResults
     /*TODO: remove rotationalDisplacementQuaternion since the last three components of the displacments
  * vector contains the same info using euler angles
 */
+    inline const static std::string defaultJsonFilename{"SimulationResults.json"};
     bool converged{false};
     std::shared_ptr<SimulationJob> pJob;
     SimulationHistory history;
@@ -553,6 +556,16 @@ struct SimulationResults
         std::filesystem::create_directories(resultsFolderPath);
         Eigen::write_binary(std::filesystem::path(resultsFolderPath).append(filename).string(), m);
 
+        nlohmann::json json;
+
+        json[GET_VARIABLE_NAME(internalPotentialEnergy)] = internalPotentialEnergy;
+
+        std::filesystem::path jsonFilePath(
+            std::filesystem::path(resultsFolderPath).append(defaultJsonFilename));
+        std::ofstream jsonFile(jsonFilePath.string());
+        jsonFile << json;
+        jsonFile.close();
+
         saveDeformedModel(resultsFolderPath.string());
     }
 
@@ -578,6 +591,28 @@ struct SimulationResults
         this->pJob = pJob;
         load(loadFromPath);
     }
+    static double computeDistance(
+        const SimulationResults &resultsA,
+        const SimulationResults &resultsB,
+        const std::unordered_map<VertexIndex, VertexIndex> &resultsAToResultsBViMap)
+    {
+        double distance = 0;
+        for (std::pair<int, int> resultsAToResultsBViPair : resultsAToResultsBViMap) {
+            const double vertexToVertexDistance
+                = (resultsA.displacements[resultsAToResultsBViPair.first].getTranslation()
+                   - resultsB.displacements[resultsAToResultsBViPair.second].getTranslation())
+                      .norm();
+            distance += vertexToVertexDistance;
+        }
+        return distance;
+    }
+
+    double computeDistance(const SimulationResults &other,
+                           const std::unordered_map<VertexIndex, VertexIndex> &thisToOtherViMap) const
+    {
+        return computeDistance(*this, other, thisToOtherViMap);
+    }
+
 #ifdef POLYSCOPE_DEFINED
     void unregister() const
     {
@@ -732,7 +767,7 @@ struct SimulationResults
     }
 #endif
 private:
-    void load(const std::filesystem::path &loadFromPath)
+    bool load(const std::filesystem::path &loadFromPath)
     {
         converged = true; //assuming it has converged
         assert(pJob != nullptr);
@@ -756,6 +791,19 @@ private:
                                                    * Eigen::AngleAxisd(displacements[vi][5],
                                                                        Eigen::Vector3d::UnitZ());
         }
+
+        const std::filesystem::path jsonFilepath = std::filesystem::path(loadFromPath)
+                                                       .append(defaultJsonFilename);
+        if (!std::filesystem::exists(jsonFilepath)) {
+            std::cerr << "Simulation results could not be loaded because filepath does "
+                         "not exist:"
+                      << jsonFilepath << std::endl;
+            return false;
+        }
+        std::ifstream ifs(jsonFilepath);
+        nlohmann::json json;
+        ifs >> json;
+        internalPotentialEnergy = json.at(GET_VARIABLE_NAME(internalPotentialEnergy));
     }
 };
 
diff --git a/trianglepatterngeometry.cpp b/trianglepatterngeometry.cpp
index 9d30d9a..0031a23 100755
--- a/trianglepatterngeometry.cpp
+++ b/trianglepatterngeometry.cpp
@@ -86,39 +86,39 @@ std::vector<vcg::Point3d> PatternGeometry::getVertices() const {
   return verts;
 }
 
-PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern) {
+PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern)
+{
+    const size_t fanSize = PatternGeometry().getFanSize();
+    PatternGeometry fan(createFan(pattern));
+    PatternGeometry tile(fan);
 
-  const size_t fanSize = PatternGeometry().getFanSize();
-  PatternGeometry fan(createFan(pattern));
-  PatternGeometry tile(fan);
+    if (fanSize % 2 == 1) {
+        vcg::Matrix44d R;
+        auto rotationAxis = vcg::Point3d(0, 0, 1);
+        R.SetRotateDeg(180, rotationAxis);
+        vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, R);
+    }
+    vcg::Matrix44d T;
+    const double centerAngle = 2 * M_PI / fanSize;
+    const double triangleHeight = std::sin((M_PI - centerAngle) / 2)
+                                  * PatternGeometry().triangleEdgeSize;
+    T.SetTranslate(0, -2 * triangleHeight, 0);
+    vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, T);
 
-  if (fanSize % 2 == 1) {
-    vcg::Matrix44d R;
-    auto rotationAxis = vcg::Point3d(0, 0, 1);
-    R.SetRotateDeg(180, rotationAxis);
-    vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, R);
-  }
-  vcg::Matrix44d T;
-  const double centerAngle = 2 * M_PI / fanSize;
-  const double triangleHeight =
-      std::sin((M_PI - centerAngle) / 2) * PatternGeometry().triangleEdgeSize;
-  T.SetTranslate(0, -2 * triangleHeight, 0);
-  vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, T);
+    PatternGeometry fanOfFan = createFan(fan);
+    vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(tile, fanOfFan);
+    vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(tile, 0.0000005);
+    vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(tile);
+    vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(tile);
+    vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(tile);
 
-  PatternGeometry fanOfFan = createFan(fan);
-  vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(tile, fanOfFan);
-  vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(tile, 0.0000005);
-  vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(tile);
-  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(tile);
-  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(tile);
+    //  for (size_t vi = 0; vi < pattern.vn; vi++) {
+    //    tile.vert[vi].C() = vcg::Color4b::Blue;
+    //  }
 
-  //  for (size_t vi = 0; vi < pattern.vn; vi++) {
-  //    tile.vert[vi].C() = vcg::Color4b::Blue;
-  //  }
-
-  tile.setLabel("tilled_" + pattern.getLabel());
-  tile.updateEigenEdgeAndVertices();
-  return tile;
+    tile.setLabel("tilled_" + pattern.getLabel());
+    tile.updateEigenEdgeAndVertices();
+    return tile;
 }
 
 PatternGeometry PatternGeometry::createFan(PatternGeometry &pattern) {
@@ -145,9 +145,9 @@ PatternGeometry::PatternGeometry(PatternGeometry &other) {
   vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
 }
 
-bool PatternGeometry::load(const std::string &plyFilename)
+bool PatternGeometry::load(const std::filesystem::__cxx11::path &meshFilePath)
 {
-    if (!VCGEdgeMesh::load(plyFilename)) {
+    if (!VCGEdgeMesh::load(meshFilePath)) {
         return false;
     }
     vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
@@ -777,6 +777,22 @@ double PatternGeometry::computeBaseTriangleHeight() const
     return vcg::Distance(vert[0].cP(), vert[interfaceNodeIndex].cP());
 }
 
+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?
+    }
+}
+
+void PatternGeometry::deleteDanglingVertices(
+    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu)
+{
+    VCGEdgeMesh::deleteDanglingVertices(pu);
+}
+
 vcg::Triangle3<double> PatternGeometry::getBaseTriangle() const
 {
     return baseTriangle;
@@ -814,6 +830,7 @@ PatternGeometry::PatternGeometry(
   updateEigenEdgeAndVertices();
 }
 
+//TODO: refactor such that it takes as an input a single pattern and tiles into the desired faces without requiring the each face will contain a single pattern
 std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
     std::vector<ConstPatternGeometry> &patterns,
     const std::vector<int> &connectToNeighborsVi,
diff --git a/trianglepatterngeometry.hpp b/trianglepatterngeometry.hpp
index ccce648..150912f 100755
--- a/trianglepatterngeometry.hpp
+++ b/trianglepatterngeometry.hpp
@@ -46,7 +46,7 @@ private:
    * vcglib interface.
    * */
   PatternGeometry(PatternGeometry &other);
-  bool load(const std::string &plyFilename) override;
+  bool load(const std::filesystem::path &meshFilePath) override;
   void add(const std::vector<vcg::Point3d> &vertices);
   void add(const std::vector<vcg::Point2i> &edges);
   void add(const std::vector<vcg::Point3d> &vertices, const std::vector<vcg::Point2i> &edges);
@@ -124,6 +124,9 @@ private:
       std::vector<std::vector<size_t>> &perPatternIndexTiledPatternEdgeIndex);
   std::unordered_set<VertexIndex> getInterfaceNodes(const std::vector<size_t> &numberOfNodesPerSlot);
   bool isInterfaceConnected(const std::unordered_set<VertexIndex> &interfaceNodes);
+  void deleteDanglingVertices() override;
+  void deleteDanglingVertices(
+      vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu) override;
 };
 
 #endif // FLATPATTERNGEOMETRY_HPP
diff --git a/vcgtrimesh.cpp b/vcgtrimesh.cpp
index 486c619..c40696d 100755
--- a/vcgtrimesh.cpp
+++ b/vcgtrimesh.cpp
@@ -7,8 +7,8 @@
 #include <wrap/io_trimesh/import.h>
 //#include <wrap/nanoply/include/nanoplyWrapper.hpp>
 
-bool VCGTriMesh::load(const std::string &filename) {
-    assert(std::filesystem::exists(filename));
+bool VCGTriMesh::load(const std::filesystem::__cxx11::path &meshFilePath) {
+    assert(std::filesystem::exists(meshFilePath));
     unsigned int mask = 0;
     mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
     mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
@@ -20,7 +20,7 @@ bool VCGTriMesh::load(const std::string &filename) {
     //            std::filesystem::absolute(filename).string().c_str(), *this, mask)
     //        != 0) {
     if (vcg::tri::io::Importer<VCGTriMesh>::Open(*this,
-                                                 std::filesystem::absolute(filename).string().c_str())
+                                                 std::filesystem::absolute(meshFilePath).string().c_str())
         != 0) {
         std::cout << "Could not load tri mesh" << std::endl;
         return false;
@@ -32,7 +32,7 @@ bool VCGTriMesh::load(const std::string &filename) {
     vcg::tri::UpdateTopology<VCGTriMesh>::VertexEdge(*this);
     vcg::tri::UpdateNormal<VCGTriMesh>::PerVertexNormalized(*this);
 
-    label = std::filesystem::path(filename).stem().string();
+    label = std::filesystem::path(meshFilePath).stem().string();
     return true;
 }
 
diff --git a/vcgtrimesh.hpp b/vcgtrimesh.hpp
index 6f79d02..c871075 100755
--- a/vcgtrimesh.hpp
+++ b/vcgtrimesh.hpp
@@ -39,7 +39,7 @@ class VCGTriMesh : public vcg::tri::TriMesh<std::vector<VCGTriMeshVertex>,
 public:
   VCGTriMesh();
   VCGTriMesh(const std::string &filename);
-  bool load(const std::string &filename) override;
+  bool load(const std::filesystem::path &meshFilePath) override;
   Eigen::MatrixX3d getVertices() const;
   Eigen::MatrixX3i getFaces() const;
   bool save(const std::string plyFilename);