Removed file

Hardcoded the 22 hexagon test-surface
Refactoring. Replaced PatternGeometry with ReducedModel everywhere
2022-05-06 16:40:31 +03:00 · 2022-05-06 16:40:05 +03:00 · 2022-05-06 16:38:50 +03:00 · 2022-05-06 16:36:47 +03:00 · 2022-05-06 16:28:28 +03:00 · 2022-05-06 16:27:37 +03:00
25 changed files with 5261 additions and 1467 deletions
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 2.8)
+cmake_minimum_required(VERSION 3.0)
 project(MySources)

 file(GLOB MySourcesFiles ${CMAKE_CURRENT_LIST_DIR}/*.hpp ${CMAKE_CURRENT_LIST_DIR}/*.cpp)
@ -18,6 +18,24 @@ if(NOT EXISTS ${EXTERNAL_DEPS_DIR})
 endif()
 ##Create directory for the external libraries
 file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR})
+
+##dlib
+#set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib)
+#download_project(PROJ               DLIB
+#                 GIT_REPOSITORY     https://github.com/davisking/dlib.git
+#                 GIT_TAG	    master
+#                 BINARY_DIR         ${DLIB_BIN_DIR}
+#                 PREFIX             ${EXTERNAL_DEPS_DIR}
+#                 ${UPDATE_DISCONNECTED_IF_AVAILABLE}
+#)
+#add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BIN_DIR})
+#if(${MYSOURCES_STATIC_LINK})
+#        target_link_libraries(${PROJECT_NAME} PUBLIC -static dlib::dlib)
+#else()
+#        target_link_libraries(${PROJECT_NAME} PUBLIC dlib::dlib)
+#endif()
+#add_compile_definitions(DLIB_DEFINED)
+
 if(${USE_POLYSCOPE})
        set(POLYSCOPE_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/polyscope)
        download_project(PROJ               POLYSCOPE
@ -30,6 +48,10 @@ if(${USE_POLYSCOPE})
        add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR})
        add_compile_definitions(POLYSCOPE_DEFINED)
        target_sources(${PROJECT_NAME} PUBLIC ${POLYSCOPE_SOURCE_DIR}/deps/imgui/imgui/misc/cpp/imgui_stdlib.cpp)
+        if(${USE_POLYSCOPE})
+                message("Using polyscope here")
+                target_link_libraries(${PROJECT_NAME} PUBLIC polyscope)
+        endif()
 endif()

 ##vcglib devel branch
@ -40,6 +62,9 @@ download_project(PROJ		    vcglib_devel
                 ${UPDATE_DISCONNECTED_IF_AVAILABLE}
 )
 add_subdirectory(${vcglib_devel_SOURCE_DIR} ${vcglib_devel_BINARY_DIR})
+
+#add_subdirectory("/home/iason/Coding/Libraries/vcglib" ${CMAKE_CURRENT_BINARY_DIR}/vcglib)
+#target_include_directories(${PROJECT_NAME} PUBLIC "/home/iason/Coding/Libraries/vcglib")
 target_sources(${PROJECT_NAME} PUBLIC ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)

 ##matplot++ lib
@ -56,7 +81,8 @@ add_subdirectory(${MATPLOTPLUSPLUS_SOURCE_DIR} ${MATPLOTPLUSPLUS_BINARY_DIR})
 ##threed-beam-fea
 set(threed-beam-fea_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/threed-beam-fea)
 download_project(PROJ                threed-beam-fea
-                 GIT_REPOSITORY      https://github.com/IasonManolas/threed-beam-fea.git
+#                 GIT_REPOSITORY      https://github.com/IasonManolas/threed-beam-fea.git
+                 GIT_REPOSITORY      https://gitea-s2i2s.isti.cnr.it/manolas/threed-beam-fea.git
                 GIT_TAG             master
                 BINARY_DIR		${threed-beam-fea_BINARY_DIR}
                 PREFIX             ${EXTERNAL_DEPS_DIR}
@ -76,6 +102,7 @@ download_project(PROJ                TBB
 option(TBB_BUILD_TESTS           "Build TBB tests and enable testing infrastructure" OFF)
 add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR})
 link_directories(${TBB_BINARY_DIR})
+
 ###Eigen 3 NOTE: Eigen is required on the system the code is ran
 find_package(Eigen3 3.3 REQUIRED)
 if(MSVC)
@ -83,43 +110,54 @@ if(MSVC)
 endif(MSVC)
 #link_directories(${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)

+#set_target_properties(${PROJECT_NAME} PROPERTIES POSITION_INDEPENDENT_CODE TRUE)
 if(${MYSOURCES_STATIC_LINK})
-        message("Linking statically")
-        target_link_libraries(${PROJECT_NAME} PUBLIC -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a #[[tbb_static]] pthread gfortran quadmath)
-else()
-#    set_target_properties(${PROJECT_NAME} PROPERTIES POSITION_INDEPENDENT_CODE TRUE)
-        target_link_libraries(${PROJECT_NAME} PUBLIC Eigen3::Eigen matplot ThreedBeamFEA tbb pthread)
-        if(${USE_POLYSCOPE})
-        message("Using polyscope")
-                target_link_libraries(${PROJECT_NAME} PUBLIC polyscope)
-        endif()
+        message("STATIC LINK MY SOURCES:" ${MYSOURCES_STATIC_LINK})
 endif()
+if(${MYSOURCES_STATIC_LINK})
+        message("Linking statically here")
+        target_link_libraries(${PROJECT_NAME} PUBLIC -static Eigen3::Eigen matplot ThreedBeamFEA tbb_static pthread gfortran quadmath)
+else()
+        target_link_libraries(${PROJECT_NAME} PUBLIC Eigen3::Eigen matplot ThreedBeamFEA tbb pthread)
+endif()
+
 target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
 target_include_directories(${PROJECT_NAME}
    PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph
    PUBLIC ${vcglib_devel_SOURCE_DIR}
-#    PUBLIC ${threed-beam-fea_SOURCE_DIR}
-        PUBLIC ThreedBeamFEA
-        PUBLIC matplot
+    PUBLIC ${threed-beam-fea_SOURCE_DIR}
+    PUBLIC ${MATPLOTPLUSPLUS_SOURCE_DIR}/source
    )

-if(USE_ENSMALLEN)
-        ##ENSMALLEN
-        set(ENSMALLEN_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/ensmallen)
-        download_project(PROJ			ENSMALLEN
-                         GIT_REPOSITORY		https://github.com/mlpack/ensmallen.git
-                         GIT_TAG		master
-                         BINARY_DIR		${ENSMALLEN_BINARY_DIR}
-                         PREFIX             	${EXTERNAL_DEPS_DIR}
-                         ${UPDATE_DISCONNECTED_IF_AVAILABLE}
-        )
-        add_subdirectory(${ENSMALLEN_SOURCE_DIR} ${ENSMALLEN_BINARY_DIR})
-        set(ARMADILLO_SOURCE_DIR "/home/iason/Coding/Libraries/armadillo")
-        add_subdirectory(${ARMADILLO_SOURCE_DIR} ${EXTERNAL_DEPS_DIR}/armadillo_build)
-        if(${MYSOURCES_STATIC_LINK})
-                target_link_libraries(${PROJECT_NAME} PUBLIC "-static" ensmallen ${EXTERNAL_DEPS_DIR}/armadillo_build/libarmadillo.a)
-        else()
-                target_link_libraries(${PROJECT_NAME} PUBLIC  armadillo ensmallen)
-        endif()
-        target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include ${ENSMALLEN_SOURCE_DIR})
-endif()
+
+#if(USE_ENSMALLEN)
+set(ARMADILLO_SOURCE_DIR "/home/iason/Coding/Libraries/armadillo")
+set(ARMADILLO_BIN_DIR "/home/iason/Coding/Libraries/armadillo/build")
+add_subdirectory(${ARMADILLO_SOURCE_DIR} ${ARMADILLO_BIN_DIR})
+target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include)
+add_compile_definitions(ARMA_DONT_USE_WRAPPER)
+target_link_libraries(${PROJECT_NAME} PUBLIC #[["/home/iason/Coding/Libraries/armadillo/build/libarmadillo.a"]] blas lapack)
+find_package(Armadillo REQUIRED)
+target_link_libraries(${PROJECT_NAME} PUBLIC ${ARMADILLO_LIBRARIES})
+
+##ENSMALLEN
+set(ENSMALLEN_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/ensmallen)
+download_project(PROJ			ENSMALLEN
+                 GIT_REPOSITORY		https://github.com/mlpack/ensmallen.git
+                 GIT_TAG		master
+                 BINARY_DIR		${ENSMALLEN_BINARY_DIR}
+                 PREFIX             	${EXTERNAL_DEPS_DIR}
+                 ${UPDATE_DISCONNECTED_IF_AVAILABLE}
+)
+#        add_subdirectory(${ENSMALLEN_SOURCE_DIR} ${ENSMALLEN_BINARY_DIR})
+#        target_link_libraries(${PROJECT_NAME} INTERFACE ensmallen)
+target_include_directories(${PROJECT_NAME}
+PUBLIC ${ENSMALLEN_SOURCE_DIR}/include)
+add_compile_definitions(USE_ENSMALLEN)
+#endif()
+
+#find_package(OpenMP REQUIRED)
+#target_link_libraries(${PROJECT_NAME} OpenMP::OpenMP_CXX)
+#target_include_directories(${PROJECT_NAME} OpenMP::OpenMP_CXX)
+
+#message(STATUS "OpenMP_C_INCLUDE_DIRS: ${OpenMP_CXX_INCLUDE_DIRS}")
--- a/CMakeLists.txt.user
+++ b/CMakeLists.txt.user
@ -1,355 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<!DOCTYPE QtCreatorProject>
-<!-- Written by QtCreator 4.14.2, 2021-07-15T10:27:01. -->
-<qtcreator>
- <data>
-  <variable>EnvironmentId</variable>
-  <value type="QByteArray">{4ab9621f-3cb4-4d5c-a5a5-8f7ea53a6f0f}</value>
- </data>
- <data>
-  <variable>ProjectExplorer.Project.ActiveTarget</variable>
-  <value type="int">0</value>
- </data>
- <data>
-  <variable>ProjectExplorer.Project.EditorSettings</variable>
-  <valuemap type="QVariantMap">
-   <value type="bool" key="EditorConfiguration.AutoIndent">false</value>
-   <value type="bool" key="EditorConfiguration.AutoSpacesForTabs">false</value>
-   <value type="bool" key="EditorConfiguration.CamelCaseNavigation">true</value>
-   <valuemap type="QVariantMap" key="EditorConfiguration.CodeStyle.0">
-    <value type="QString" key="language">Cpp</value>
-    <valuemap type="QVariantMap" key="value">
-     <value type="QByteArray" key="CurrentPreferences">CppGlobal</value>
-    </valuemap>
-   </valuemap>
-   <valuemap type="QVariantMap" key="EditorConfiguration.CodeStyle.1">
-    <value type="QString" key="language">QmlJS</value>
-    <valuemap type="QVariantMap" key="value">
-     <value type="QByteArray" key="CurrentPreferences">QmlJSGlobal</value>
-    </valuemap>
-   </valuemap>
-   <value type="int" key="EditorConfiguration.CodeStyle.Count">2</value>
-   <value type="QByteArray" key="EditorConfiguration.Codec">UTF-8</value>
-   <value type="bool" key="EditorConfiguration.ConstrainTooltips">false</value>
-   <value type="int" key="EditorConfiguration.IndentSize">4</value>
-   <value type="bool" key="EditorConfiguration.KeyboardTooltips">false</value>
-   <value type="int" key="EditorConfiguration.MarginColumn">80</value>
-   <value type="bool" key="EditorConfiguration.MouseHiding">true</value>
-   <value type="bool" key="EditorConfiguration.MouseNavigation">true</value>
-   <value type="int" key="EditorConfiguration.PaddingMode">1</value>
-   <value type="bool" key="EditorConfiguration.ScrollWheelZooming">true</value>
-   <value type="bool" key="EditorConfiguration.ShowMargin">false</value>
-   <value type="int" key="EditorConfiguration.SmartBackspaceBehavior">0</value>
-   <value type="bool" key="EditorConfiguration.SmartSelectionChanging">true</value>
-   <value type="bool" key="EditorConfiguration.SpacesForTabs">true</value>
-   <value type="int" key="EditorConfiguration.TabKeyBehavior">0</value>
-   <value type="int" key="EditorConfiguration.TabSize">8</value>
-   <value type="bool" key="EditorConfiguration.UseGlobal">true</value>
-   <value type="int" key="EditorConfiguration.Utf8BomBehavior">1</value>
-   <value type="bool" key="EditorConfiguration.addFinalNewLine">true</value>
-   <value type="bool" key="EditorConfiguration.cleanIndentation">true</value>
-   <value type="bool" key="EditorConfiguration.cleanWhitespace">true</value>
-   <value type="QString" key="EditorConfiguration.ignoreFileTypes">*.md, *.MD, Makefile</value>
-   <value type="bool" key="EditorConfiguration.inEntireDocument">false</value>
-   <value type="bool" key="EditorConfiguration.skipTrailingWhitespace">true</value>
-  </valuemap>
- </data>
- <data>
-  <variable>ProjectExplorer.Project.PluginSettings</variable>
-  <valuemap type="QVariantMap">
-   <valuemap type="QVariantMap" key="AutoTest.ActiveFrameworks">
-    <value type="bool" key="AutoTest.Framework.Boost">true</value>
-    <value type="bool" key="AutoTest.Framework.Catch">true</value>
-    <value type="bool" key="AutoTest.Framework.GTest">true</value>
-    <value type="bool" key="AutoTest.Framework.QtQuickTest">true</value>
-    <value type="bool" key="AutoTest.Framework.QtTest">true</value>
-   </valuemap>
-   <valuemap type="QVariantMap" key="AutoTest.CheckStates"/>
-   <value type="int" key="AutoTest.RunAfterBuild">0</value>
-   <value type="bool" key="AutoTest.UseGlobal">true</value>
-   <valuemap type="QVariantMap" key="ClangTools">
-    <value type="bool" key="ClangTools.AnalyzeOpenFiles">true</value>
-    <value type="bool" key="ClangTools.BuildBeforeAnalysis">true</value>
-    <value type="QString" key="ClangTools.DiagnosticConfig">Builtin.DefaultTidyAndClazy</value>
-    <value type="int" key="ClangTools.ParallelJobs">6</value>
-    <valuelist type="QVariantList" key="ClangTools.SelectedDirs"/>
-    <valuelist type="QVariantList" key="ClangTools.SelectedFiles"/>
-    <valuelist type="QVariantList" key="ClangTools.SuppressedDiagnostics"/>
-    <value type="bool" key="ClangTools.UseGlobalSettings">true</value>
-   </valuemap>
-  </valuemap>
- </data>
- <data>
-  <variable>ProjectExplorer.Project.Target.0</variable>
-  <valuemap type="QVariantMap">
-   <value type="QString" key="DeviceType">Desktop</value>
-   <value type="QString" key="ProjectExplorer.ProjectConfiguration.DefaultDisplayName">Clang</value>
-   <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Clang</value>
-   <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">{08d5b07f-8010-48d8-a202-51411d8039c8}</value>
-   <value type="int" key="ProjectExplorer.Target.ActiveBuildConfiguration">0</value>
-   <value type="int" key="ProjectExplorer.Target.ActiveDeployConfiguration">0</value>
-   <value type="int" key="ProjectExplorer.Target.ActiveRunConfiguration">0</value>
-   <valuemap type="QVariantMap" key="ProjectExplorer.Target.BuildConfiguration.0">
-    <value type="QString" key="CMake.Initial.Parameters">-GUnix Makefiles
-DCMAKE_BUILD_TYPE:String=Debug
-DQT_QMAKE_EXECUTABLE:STRING=%{Qt:qmakeExecutable}
-DCMAKE_PREFIX_PATH:STRING=%{Qt:QT_INSTALL_PREFIX}
-DCMAKE_C_COMPILER:STRING=%{Compiler:Executable:C}
-DCMAKE_CXX_COMPILER:STRING=%{Compiler:Executable:Cxx}</value>
-    <value type="QString" key="ProjectExplorer.BuildConfiguration.BuildDirectory">/home/iason/Coding/Libraries/MySources_viscousDampingBranch/../build-MySources_viscousDampingBranch-Clang-Debug</value>
-    <valuemap type="QVariantMap" key="ProjectExplorer.BuildConfiguration.BuildStepList.0">
-     <valuemap type="QVariantMap" key="ProjectExplorer.BuildStepList.Step.0">
-      <valuelist type="QVariantList" key="CMakeProjectManager.MakeStep.BuildTargets">
-       <value type="QString">all</value>
-      </valuelist>
-      <value type="bool" key="ProjectExplorer.BuildStep.Enabled">true</value>
-      <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Build</value>
-      <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.MakeStep</value>
-     </valuemap>
-     <value type="int" key="ProjectExplorer.BuildStepList.StepsCount">1</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DefaultDisplayName">Build</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Build</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">ProjectExplorer.BuildSteps.Build</value>
-    </valuemap>
-    <valuemap type="QVariantMap" key="ProjectExplorer.BuildConfiguration.BuildStepList.1">
-     <valuemap type="QVariantMap" key="ProjectExplorer.BuildStepList.Step.0">
-      <valuelist type="QVariantList" key="CMakeProjectManager.MakeStep.BuildTargets">
-       <value type="QString">clean</value>
-      </valuelist>
-      <value type="bool" key="ProjectExplorer.BuildStep.Enabled">true</value>
-      <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Build</value>
-      <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.MakeStep</value>
-     </valuemap>
-     <value type="int" key="ProjectExplorer.BuildStepList.StepsCount">1</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DefaultDisplayName">Clean</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Clean</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">ProjectExplorer.BuildSteps.Clean</value>
-    </valuemap>
-    <value type="int" key="ProjectExplorer.BuildConfiguration.BuildStepListCount">2</value>
-    <value type="bool" key="ProjectExplorer.BuildConfiguration.ClearSystemEnvironment">false</value>
-    <valuelist type="QVariantList" key="ProjectExplorer.BuildConfiguration.CustomParsers"/>
-    <valuelist type="QVariantList" key="ProjectExplorer.BuildConfiguration.UserEnvironmentChanges"/>
-    <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Debug</value>
-    <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.CMakeBuildConfiguration</value>
-   </valuemap>
-   <valuemap type="QVariantMap" key="ProjectExplorer.Target.BuildConfiguration.1">
-    <value type="QString" key="CMake.Initial.Parameters">-GUnix Makefiles
-DCMAKE_BUILD_TYPE:String=Release
-DQT_QMAKE_EXECUTABLE:STRING=%{Qt:qmakeExecutable}
-DCMAKE_PREFIX_PATH:STRING=%{Qt:QT_INSTALL_PREFIX}
-DCMAKE_C_COMPILER:STRING=%{Compiler:Executable:C}
-DCMAKE_CXX_COMPILER:STRING=%{Compiler:Executable:Cxx}</value>
-    <value type="QString" key="ProjectExplorer.BuildConfiguration.BuildDirectory">/home/iason/Coding/Libraries/MySources_viscousDampingBranch/../build-MySources_viscousDampingBranch-Clang-Release</value>
-    <valuemap type="QVariantMap" key="ProjectExplorer.BuildConfiguration.BuildStepList.0">
-     <valuemap type="QVariantMap" key="ProjectExplorer.BuildStepList.Step.0">
-      <valuelist type="QVariantList" key="CMakeProjectManager.MakeStep.BuildTargets">
-       <value type="QString">all</value>
-      </valuelist>
-      <value type="bool" key="ProjectExplorer.BuildStep.Enabled">true</value>
-      <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.MakeStep</value>
-     </valuemap>
-     <value type="int" key="ProjectExplorer.BuildStepList.StepsCount">1</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DefaultDisplayName">Build</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Build</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">ProjectExplorer.BuildSteps.Build</value>
-    </valuemap>
-    <valuemap type="QVariantMap" key="ProjectExplorer.BuildConfiguration.BuildStepList.1">
-     <valuemap type="QVariantMap" key="ProjectExplorer.BuildStepList.Step.0">
-      <valuelist type="QVariantList" key="CMakeProjectManager.MakeStep.BuildTargets">
-       <value type="QString">clean</value>
-      </valuelist>
-      <value type="bool" key="ProjectExplorer.BuildStep.Enabled">true</value>
-      <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.MakeStep</value>
-     </valuemap>
-     <value type="int" key="ProjectExplorer.BuildStepList.StepsCount">1</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DefaultDisplayName">Clean</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Clean</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">ProjectExplorer.BuildSteps.Clean</value>
-    </valuemap>
-    <value type="int" key="ProjectExplorer.BuildConfiguration.BuildStepListCount">2</value>
-    <value type="bool" key="ProjectExplorer.BuildConfiguration.ClearSystemEnvironment">false</value>
-    <valuelist type="QVariantList" key="ProjectExplorer.BuildConfiguration.CustomParsers"/>
-    <valuelist type="QVariantList" key="ProjectExplorer.BuildConfiguration.UserEnvironmentChanges"/>
-    <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Release</value>
-    <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.CMakeBuildConfiguration</value>
-   </valuemap>
-   <valuemap type="QVariantMap" key="ProjectExplorer.Target.BuildConfiguration.2">
-    <value type="QString" key="CMake.Initial.Parameters">-GUnix Makefiles
-DCMAKE_BUILD_TYPE:String=RelWithDebInfo
-DQT_QMAKE_EXECUTABLE:STRING=%{Qt:qmakeExecutable}
-DCMAKE_PREFIX_PATH:STRING=%{Qt:QT_INSTALL_PREFIX}
-DCMAKE_C_COMPILER:STRING=%{Compiler:Executable:C}
-DCMAKE_CXX_COMPILER:STRING=%{Compiler:Executable:Cxx}</value>
-    <value type="QString" key="ProjectExplorer.BuildConfiguration.BuildDirectory">/home/iason/Coding/Libraries/MySources_viscousDampingBranch/../build-MySources_viscousDampingBranch-Clang-RelWithDebInfo</value>
-    <valuemap type="QVariantMap" key="ProjectExplorer.BuildConfiguration.BuildStepList.0">
-     <valuemap type="QVariantMap" key="ProjectExplorer.BuildStepList.Step.0">
-      <valuelist type="QVariantList" key="CMakeProjectManager.MakeStep.BuildTargets">
-       <value type="QString">all</value>
-      </valuelist>
-      <value type="bool" key="ProjectExplorer.BuildStep.Enabled">true</value>
-      <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.MakeStep</value>
-     </valuemap>
-     <value type="int" key="ProjectExplorer.BuildStepList.StepsCount">1</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DefaultDisplayName">Build</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Build</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">ProjectExplorer.BuildSteps.Build</value>
-    </valuemap>
-    <valuemap type="QVariantMap" key="ProjectExplorer.BuildConfiguration.BuildStepList.1">
-     <valuemap type="QVariantMap" key="ProjectExplorer.BuildStepList.Step.0">
-      <valuelist type="QVariantList" key="CMakeProjectManager.MakeStep.BuildTargets">
-       <value type="QString">clean</value>
-      </valuelist>
-      <value type="bool" key="ProjectExplorer.BuildStep.Enabled">true</value>
-      <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.MakeStep</value>
-     </valuemap>
-     <value type="int" key="ProjectExplorer.BuildStepList.StepsCount">1</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DefaultDisplayName">Clean</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Clean</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">ProjectExplorer.BuildSteps.Clean</value>
-    </valuemap>
-    <value type="int" key="ProjectExplorer.BuildConfiguration.BuildStepListCount">2</value>
-    <value type="bool" key="ProjectExplorer.BuildConfiguration.ClearSystemEnvironment">false</value>
-    <valuelist type="QVariantList" key="ProjectExplorer.BuildConfiguration.CustomParsers"/>
-    <valuelist type="QVariantList" key="ProjectExplorer.BuildConfiguration.UserEnvironmentChanges"/>
-    <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Release with Debug Information</value>
-    <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.CMakeBuildConfiguration</value>
-   </valuemap>
-   <valuemap type="QVariantMap" key="ProjectExplorer.Target.BuildConfiguration.3">
-    <value type="QString" key="CMake.Initial.Parameters">-GUnix Makefiles
-DCMAKE_BUILD_TYPE:String=MinSizeRel
-DQT_QMAKE_EXECUTABLE:STRING=%{Qt:qmakeExecutable}
-DCMAKE_PREFIX_PATH:STRING=%{Qt:QT_INSTALL_PREFIX}
-DCMAKE_C_COMPILER:STRING=%{Compiler:Executable:C}
-DCMAKE_CXX_COMPILER:STRING=%{Compiler:Executable:Cxx}</value>
-    <value type="QString" key="ProjectExplorer.BuildConfiguration.BuildDirectory">/home/iason/Coding/Libraries/MySources_viscousDampingBranch/../build-MySources_viscousDampingBranch-Clang-MinSizeRel</value>
-    <valuemap type="QVariantMap" key="ProjectExplorer.BuildConfiguration.BuildStepList.0">
-     <valuemap type="QVariantMap" key="ProjectExplorer.BuildStepList.Step.0">
-      <valuelist type="QVariantList" key="CMakeProjectManager.MakeStep.BuildTargets">
-       <value type="QString">all</value>
-      </valuelist>
-      <value type="bool" key="ProjectExplorer.BuildStep.Enabled">true</value>
-      <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.MakeStep</value>
-     </valuemap>
-     <value type="int" key="ProjectExplorer.BuildStepList.StepsCount">1</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DefaultDisplayName">Build</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Build</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">ProjectExplorer.BuildSteps.Build</value>
-    </valuemap>
-    <valuemap type="QVariantMap" key="ProjectExplorer.BuildConfiguration.BuildStepList.1">
-     <valuemap type="QVariantMap" key="ProjectExplorer.BuildStepList.Step.0">
-      <valuelist type="QVariantList" key="CMakeProjectManager.MakeStep.BuildTargets">
-       <value type="QString">clean</value>
-      </valuelist>
-      <value type="bool" key="ProjectExplorer.BuildStep.Enabled">true</value>
-      <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.MakeStep</value>
-     </valuemap>
-     <value type="int" key="ProjectExplorer.BuildStepList.StepsCount">1</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DefaultDisplayName">Clean</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Clean</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">ProjectExplorer.BuildSteps.Clean</value>
-    </valuemap>
-    <value type="int" key="ProjectExplorer.BuildConfiguration.BuildStepListCount">2</value>
-    <value type="bool" key="ProjectExplorer.BuildConfiguration.ClearSystemEnvironment">false</value>
-    <valuelist type="QVariantList" key="ProjectExplorer.BuildConfiguration.CustomParsers"/>
-    <valuelist type="QVariantList" key="ProjectExplorer.BuildConfiguration.UserEnvironmentChanges"/>
-    <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Minimum Size Release</value>
-    <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">CMakeProjectManager.CMakeBuildConfiguration</value>
-   </valuemap>
-   <value type="int" key="ProjectExplorer.Target.BuildConfigurationCount">4</value>
-   <valuemap type="QVariantMap" key="ProjectExplorer.Target.DeployConfiguration.0">
-    <valuemap type="QVariantMap" key="ProjectExplorer.BuildConfiguration.BuildStepList.0">
-     <value type="int" key="ProjectExplorer.BuildStepList.StepsCount">0</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DefaultDisplayName">Deploy</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Deploy</value>
-     <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">ProjectExplorer.BuildSteps.Deploy</value>
-    </valuemap>
-    <value type="int" key="ProjectExplorer.BuildConfiguration.BuildStepListCount">1</value>
-    <valuemap type="QVariantMap" key="ProjectExplorer.DeployConfiguration.CustomData"/>
-    <value type="bool" key="ProjectExplorer.DeployConfiguration.CustomDataEnabled">false</value>
-    <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">ProjectExplorer.DefaultDeployConfiguration</value>
-   </valuemap>
-   <value type="int" key="ProjectExplorer.Target.DeployConfigurationCount">1</value>
-   <valuemap type="QVariantMap" key="ProjectExplorer.Target.RunConfiguration.0">
-    <value type="QString" key="Analyzer.Perf.CallgraphMode">dwarf</value>
-    <valuelist type="QVariantList" key="Analyzer.Perf.Events">
-     <value type="QString">cpu-cycles</value>
-    </valuelist>
-    <valuelist type="QVariantList" key="Analyzer.Perf.ExtraArguments"/>
-    <value type="int" key="Analyzer.Perf.Frequency">250</value>
-    <valuelist type="QVariantList" key="Analyzer.Perf.RecordArguments">
-     <value type="QString">-e</value>
-     <value type="QString">cpu-cycles</value>
-     <value type="QString">--call-graph</value>
-     <value type="QString">dwarf,4096</value>
-     <value type="QString">-F</value>
-     <value type="QString">250</value>
-    </valuelist>
-    <value type="QString" key="Analyzer.Perf.SampleMode">-F</value>
-    <value type="bool" key="Analyzer.Perf.Settings.UseGlobalSettings">true</value>
-    <value type="int" key="Analyzer.Perf.StackSize">4096</value>
-    <value type="bool" key="Analyzer.QmlProfiler.AggregateTraces">false</value>
-    <value type="bool" key="Analyzer.QmlProfiler.FlushEnabled">false</value>
-    <value type="uint" key="Analyzer.QmlProfiler.FlushInterval">1000</value>
-    <value type="QString" key="Analyzer.QmlProfiler.LastTraceFile"></value>
-    <value type="bool" key="Analyzer.QmlProfiler.Settings.UseGlobalSettings">true</value>
-    <valuelist type="QVariantList" key="Analyzer.Valgrind.AddedSuppressionFiles"/>
-    <value type="bool" key="Analyzer.Valgrind.Callgrind.CollectBusEvents">false</value>
-    <value type="bool" key="Analyzer.Valgrind.Callgrind.CollectSystime">false</value>
-    <value type="bool" key="Analyzer.Valgrind.Callgrind.EnableBranchSim">false</value>
-    <value type="bool" key="Analyzer.Valgrind.Callgrind.EnableCacheSim">false</value>
-    <value type="bool" key="Analyzer.Valgrind.Callgrind.EnableEventToolTips">true</value>
-    <value type="double" key="Analyzer.Valgrind.Callgrind.MinimumCostRatio">0.01</value>
-    <value type="double" key="Analyzer.Valgrind.Callgrind.VisualisationMinimumCostRatio">10</value>
-    <value type="bool" key="Analyzer.Valgrind.FilterExternalIssues">true</value>
-    <value type="QString" key="Analyzer.Valgrind.KCachegrindExecutable">kcachegrind</value>
-    <value type="int" key="Analyzer.Valgrind.LeakCheckOnFinish">1</value>
-    <value type="int" key="Analyzer.Valgrind.NumCallers">25</value>
-    <valuelist type="QVariantList" key="Analyzer.Valgrind.RemovedSuppressionFiles"/>
-    <value type="int" key="Analyzer.Valgrind.SelfModifyingCodeDetection">1</value>
-    <value type="bool" key="Analyzer.Valgrind.Settings.UseGlobalSettings">true</value>
-    <value type="bool" key="Analyzer.Valgrind.ShowReachable">false</value>
-    <value type="bool" key="Analyzer.Valgrind.TrackOrigins">true</value>
-    <value type="QString" key="Analyzer.Valgrind.ValgrindExecutable">valgrind</value>
-    <valuelist type="QVariantList" key="Analyzer.Valgrind.VisibleErrorKinds">
-     <value type="int">0</value>
-     <value type="int">1</value>
-     <value type="int">2</value>
-     <value type="int">3</value>
-     <value type="int">4</value>
-     <value type="int">5</value>
-     <value type="int">6</value>
-     <value type="int">7</value>
-     <value type="int">8</value>
-     <value type="int">9</value>
-     <value type="int">10</value>
-     <value type="int">11</value>
-     <value type="int">12</value>
-     <value type="int">13</value>
-     <value type="int">14</value>
-    </valuelist>
-    <valuelist type="QVariantList" key="CustomOutputParsers"/>
-    <value type="int" key="PE.EnvironmentAspect.Base">2</value>
-    <valuelist type="QVariantList" key="PE.EnvironmentAspect.Changes"/>
-    <value type="QString" key="ProjectExplorer.ProjectConfiguration.Id">ProjectExplorer.CustomExecutableRunConfiguration</value>
-    <value type="QString" key="ProjectExplorer.RunConfiguration.BuildKey"></value>
-    <value type="bool" key="RunConfiguration.UseCppDebugger">false</value>
-    <value type="bool" key="RunConfiguration.UseCppDebuggerAuto">true</value>
-    <value type="bool" key="RunConfiguration.UseQmlDebugger">false</value>
-    <value type="bool" key="RunConfiguration.UseQmlDebuggerAuto">true</value>
-   </valuemap>
-   <value type="int" key="ProjectExplorer.Target.RunConfigurationCount">1</value>
-  </valuemap>
- </data>
- <data>
-  <variable>ProjectExplorer.Project.TargetCount</variable>
-  <value type="int">1</value>
- </data>
- <data>
-  <variable>ProjectExplorer.Project.Updater.FileVersion</variable>
-  <value type="int">22</value>
- </data>
- <data>
-  <variable>Version</variable>
-  <value type="int">22</value>
- </data>
-</qtcreator>
--- a/csvfile.hpp
+++ b/csvfile.hpp
@ -65,10 +65,9 @@ class csvFile {
      return write(val);
  }

-  template<typename T>
-  static std::vector<std::vector<T>> parse(const std::filesystem::path &csvFilepath)
+  static std::vector<std::vector<std::string>> parse(const std::filesystem::path &csvFilepath)
  {
-      std::vector<std::vector<T>> resultCSV;
+      std::vector<std::vector<std::string>> resultCSV;
      if (!std::filesystem::exists(csvFilepath)) {
          std::cerr << "The file does not exist:" << csvFilepath.string() << std::endl;
          return resultCSV;
@ -79,21 +78,17 @@ class csvFile {
          std::cerr << "Can't open file:" << csvFilepath.string() << std::endl;
          return resultCSV;
      }
-      std::vector<T> row;
+      std::vector<std::string> row;
      std::string line;
      using Tokenizer = boost::tokenizer<boost::escaped_list_separator<char>>;
      while (std::getline(inputfile, line)) {
          Tokenizer tokenizer(line);
-          row.resize(std::distance(tokenizer.begin(), tokenizer.end()));
-          std::transform(tokenizer.begin(), tokenizer.end(), row.begin(), [](const std::string &el) {
-              return boost::lexical_cast<T>(el);
-          });
-          //          std::cout << std::endl;
-          //          row.assign(tokenizer.begin(), tokenizer.end());
-          //          for (const auto &el : row) {
-          //              std::cout << el << " ";
-          //          }
-          //          std::cout << std::endl;
+          const int numOfCols = std::distance(tokenizer.begin(), tokenizer.end());
+          row.resize(numOfCols);
+          std::copy(tokenizer.begin(), tokenizer.end(), row.begin());
+          //          std::transform(tokenizer.begin(), tokenizer.end(), row.begin(), [](const std::string &el) {
+          //              return boost::lexical_cast<T>(el);
+          //          });
          resultCSV.push_back(row);
      }

--- a/drmsimulationmodel.cpp
+++ b/drmsimulationmodel.cpp
--- a/drmsimulationmodel.hpp
+++ b/drmsimulationmodel.hpp
@ -16,7 +16,7 @@
 #include <ensmallen.hpp>
 #endif

-struct SimulationJob;
+class SimulationJob;

 enum DoF { Ux = 0, Uy, Uz, Nx, Ny, Nr, NumDoF };
 using DoFType = int;
@ -33,6 +33,8 @@ class DRMSimulationModel : public SimulationModel
 public:
    struct Settings
    {
+        bool useTranslationalKineticEnergyForKineticDamping{true};
+        bool useTotalRotationalKineticEnergyForKineticDamping{false};
        bool shouldDraw{false};
        bool beVerbose{false};
        bool shouldCreatePlots{false};
@ -41,11 +43,11 @@ public:
        //        double residualForcesMovingAverageNormThreshold{1e-8};
        double Dtini{0.1};
        double xi{0.9969};
-        std::optional<double> shouldUseTranslationalKineticEnergyThreshold;
+        std::optional<double> translationalKineticEnergyThreshold;
        int gradualForcedDisplacementSteps{50};
        //        int desiredGradualExternalLoadsSteps{1};
        double gamma{0.8};
-        double totalResidualForcesNormThreshold{1e-8};
+        double totalResidualForcesNormThreshold{1e-20};
        double totalExternalForcesNormPercentageTermination{1e-5};
        std::optional<int> maxDRMIterations;
        std::optional<int> debugModeStep;
@ -55,7 +57,6 @@ public:
        //        std::optional<int> intermediateResultsSaveStep;
        std::optional<bool> saveIntermediateBestStates;
        std::optional<double> viscousDampingFactor;
-        bool useKineticDamping{true};
        Settings() {}
        void save(const std::filesystem::path &jsonFilePath) const;
        bool load(const std::filesystem::path &filePath);
@ -153,7 +154,7 @@ private:
        const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);

 #ifdef POLYSCOPE_DEFINED
-    void draw(const std::string &screenshotsFolder = {});
+    void draw(const std::shared_ptr<SimulationJob> &pJob, const std::string &screenshotsFolder = {});
 #endif
  void
  updateNodalInternalForce(Vector6d &nodalInternalForce,
@ -254,17 +255,21 @@ private:

  void reset(const std::shared_ptr<SimulationJob> &pJob);

+  std::vector<std::array<Vector6d, 4>> computeInternalForces(
+      const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
+
  public:
  DRMSimulationModel();
  SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
                                      const Settings &settings,
                                      const SimulationResults &solutionGuess = SimulationResults());
-#ifdef USE_ENSMALLEN
-  double EvaluateWithGradient(const arma::mat &x, arma::mat &g);
-  void setJob(const std::shared_ptr<SimulationJob> &pJob);
-  SimulationMesh *getDeformedMesh(const arma::mat &x);
-
-#endif
+  //#ifdef USE_ENSMALLEN
+  //  std::shared_ptr<SimulationJob> pJob;
+  //  double EvaluateWithGradient(const arma::mat &x, arma::mat &g);
+  //  void setJob(const std::shared_ptr<SimulationJob> &pJob);
+  //  SimulationMesh *getDeformedMesh(const arma::mat &x, const std::shared_ptr<SimulationJob> &pJob);
+  //  double Evaluate(const arma::mat &x);
+  //#endif

  static void runUnitTests();
 };
--- a/edgemesh.cpp
+++ b/edgemesh.cpp
@ -1,14 +1,15 @@
 #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>

 Eigen::MatrixX2i VCGEdgeMesh::getEigenEdges() const { return eigenEdges; }

-std::vector<vcg::Point2i> VCGEdgeMesh::getEdges()
+std::vector<vcg::Point2i> VCGEdgeMesh::computeEdges()
 {
-    getEdges(eigenEdges);
+    computeEdges(eigenEdges);
    std::vector<vcg::Point2i> edges(eigenEdges.rows());
    for (int ei = 0; ei < eigenEdges.rows(); ei++) {
        edges[ei] = vcg::Point2i(eigenEdges(ei, 0), eigenEdges(ei, 1));
@ -16,16 +17,17 @@ std::vector<vcg::Point2i> VCGEdgeMesh::getEdges()
    return edges;
 }

-Eigen::MatrixX3d VCGEdgeMesh::getEigenVertices() {
-    getVertices(eigenVertices);
-  return eigenVertices;
+Eigen::MatrixX3d VCGEdgeMesh::getEigenVertices() const
+{
+    //    getVertices(eigenVertices);
+    return eigenVertices;
 }

 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()) {
@ -40,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;
    }
@ -182,8 +187,8 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth,
  vcg::tri::Allocator<VCGEdgeMesh>::DeleteVertex(
      *this, vert[(desiredHeight + 1) * (desiredWidth + 1) - 1]);
  vcg::tri::Allocator<VCGEdgeMesh>::CompactVertexVector(*this);
-  getEdges(eigenEdges);
-  getVertices(eigenVertices);
+  computeEdges(eigenEdges);
+  computeVertices(eigenVertices);
  //  vcg::tri::Allocator<VCGEdgeMesh>::CompactEdgeVector(*this);

  //  const size_t numberOfEdges =
@ -214,8 +219,8 @@ bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath)
        return false;
    }

-    getEdges(eigenEdges);
-    getVertices(eigenVertices);
+    computeEdges(eigenEdges);
+    computeVertices(eigenVertices);
    vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
    label = std::filesystem::path(meshFilePath).stem().string();

@ -246,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));
@ -258,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;
@ -299,25 +308,24 @@ Eigen::MatrixX3d VCGEdgeMesh::getNormals() const {

  return vertexNormals;
 }
-void VCGEdgeMesh::getEdges(Eigen::MatrixX3d &edgeStartingPoints,
-                           Eigen::MatrixX3d &edgeEndingPoints) const {
-  edgeStartingPoints.resize(EN(), 3);
-  edgeEndingPoints.resize(EN(), 3);
-  for (int edgeIndex = 0; edgeIndex < EN(); edgeIndex++) {
-    const VCGEdgeMesh::EdgeType &edge = this->edge[edgeIndex];
-    edgeStartingPoints.row(edgeIndex) =
-        edge.cP(0).ToEigenVector<Eigen::Vector3d>();
-    edgeEndingPoints.row(edgeIndex) =
-        edge.cP(1).ToEigenVector<Eigen::Vector3d>();
-  }
+void VCGEdgeMesh::computeEdges(Eigen::MatrixX3d &edgeStartingPoints,
+                               Eigen::MatrixX3d &edgeEndingPoints) const
+{
+    edgeStartingPoints.resize(EN(), 3);
+    edgeEndingPoints.resize(EN(), 3);
+    for (int edgeIndex = 0; edgeIndex < EN(); edgeIndex++) {
+        const VCGEdgeMesh::EdgeType &edge = this->edge[edgeIndex];
+        edgeStartingPoints.row(edgeIndex) = edge.cP(0).ToEigenVector<Eigen::Vector3d>();
+        edgeEndingPoints.row(edgeIndex) = edge.cP(1).ToEigenVector<Eigen::Vector3d>();
+    }
 }

 VCGEdgeMesh::VCGEdgeMesh() {}

 void VCGEdgeMesh::updateEigenEdgeAndVertices() {
 #ifdef POLYSCOPE_DEFINED
-  getEdges(eigenEdges);
-  getVertices(eigenVertices);
+    computeEdges(eigenEdges);
+    computeVertices(eigenVertices);
 #endif
 }

@ -325,13 +333,15 @@ bool VCGEdgeMesh::copy(VCGEdgeMesh &mesh) {
  vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::MeshCopy(*this, mesh);
  label = mesh.getLabel();
  eigenEdges = mesh.getEigenEdges();
-  if (eigenEdges.rows() == 0) {
-    getEdges(eigenEdges);
-  }
+  //  assert(eigenEdges.rows() != 0);
+  //  if (eigenEdges.rows() == 0) {
+  //    getEdges(eigenEdges);
+  //  }
  eigenVertices = mesh.getEigenVertices();
-  if (eigenVertices.rows() == 0) {
-    getVertices(eigenVertices);
-  }
+  //  assert(eigenVertices.rows() != 0);
+  //  if (eigenVertices.rows() == 0) {
+  //    getVertices(eigenVertices);
+  //  }
  vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);

  return true;
@ -356,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();
 }
@ -399,7 +409,8 @@ void VCGEdgeMesh::deleteDanglingVertices(vcg::tri::Allocator<VCGEdgeMesh>::Point
    updateEigenEdgeAndVertices();
 }

-void VCGEdgeMesh::getVertices(Eigen::MatrixX3d &vertices) {
+void VCGEdgeMesh::computeVertices(Eigen::MatrixX3d &vertices)
+{
    vertices = Eigen::MatrixX3d();
    vertices.resize(VN(), 3);
    for (int vi = 0; vi < VN(); vi++) {
@ -412,21 +423,21 @@ void VCGEdgeMesh::getVertices(Eigen::MatrixX3d &vertices) {
    }
 }

-void VCGEdgeMesh::getEdges(Eigen::MatrixX2i &edges) {
-  edges = Eigen::MatrixX2i();
-  edges.resize(EN(), 2);
-  for (int edgeIndex = 0; edgeIndex < EN(); edgeIndex++) {
-    const VCGEdgeMesh::EdgeType &edge = this->edge[edgeIndex];
-    assert(!edge.IsD());
-    auto vp0 = edge.cV(0);
-    auto vp1 = edge.cV(1);
-    assert(vcg::tri::IsValidPointer(*this, vp0) &&
-           vcg::tri::IsValidPointer(*this, vp1));
-    const size_t vi0 = vcg::tri::Index<VCGEdgeMesh>(*this, vp0);
-    const size_t vi1 = vcg::tri::Index<VCGEdgeMesh>(*this, vp1);
-    assert(vi0 != -1 && vi1 != -1);
-    edges.row(edgeIndex) = Eigen::Vector2i(vi0, vi1);
-  }
+void VCGEdgeMesh::computeEdges(Eigen::MatrixX2i &edges)
+{
+    edges = Eigen::MatrixX2i();
+    edges.resize(EN(), 2);
+    for (int edgeIndex = 0; edgeIndex < EN(); edgeIndex++) {
+        const VCGEdgeMesh::EdgeType &edge = this->edge[edgeIndex];
+        assert(!edge.IsD());
+        auto vp0 = edge.cV(0);
+        auto vp1 = edge.cV(1);
+        assert(vcg::tri::IsValidPointer(*this, vp0) && vcg::tri::IsValidPointer(*this, vp1));
+        const size_t vi0 = vcg::tri::Index<VCGEdgeMesh>(*this, vp0);
+        const size_t vi1 = vcg::tri::Index<VCGEdgeMesh>(*this, vp1);
+        assert(vi0 != -1 && vi1 != -1);
+        edges.row(edgeIndex) = Eigen::Vector2i(vi0, vi1);
+    }
 }

 void VCGEdgeMesh::printVertexCoordinates(const size_t &vi) const {
@ -435,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,
@ -456,7 +483,7 @@ polyscope::CurveNetwork *VCGEdgeMesh::registerForDrawing(
        const glm::vec3 desiredColor_glm(desiredColor.value()[0],
                                         desiredColor.value()[1],
                                         desiredColor.value()[2]);
-        polyscopeHandle_edgeMesh->setColor(desiredColor_glm);
+        polyscopeHandle_edgeMesh->setColor(/*glm::normalize(*/ desiredColor_glm /*)*/);
    }

    return polyscopeHandle_edgeMesh;
--- 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>
@ -40,8 +39,8 @@ protected:
    Eigen::MatrixX3d eigenVertices;
    Eigen::MatrixX3d eigenEdgeNormals;

-    void getEdges(Eigen::MatrixX2i &edges);
-    void getVertices(Eigen::MatrixX3d &vertices);
+    void computeEdges(Eigen::MatrixX2i &edges);
+    void computeVertices(Eigen::MatrixX3d &vertices);

 public:
    VCGEdgeMesh();
@ -65,7 +64,7 @@ public:
    virtual void deleteDanglingVertices(
        vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu);

-    void getEdges(Eigen::MatrixX3d &edgeStartingPoints, Eigen::MatrixX3d &edgeEndingPoints) const;
+    void computeEdges(Eigen::MatrixX3d &edgeStartingPoints, Eigen::MatrixX3d &edgeEndingPoints) const;

    Eigen::MatrixX3d getNormals() const;

@ -81,9 +80,9 @@ public:
    void createSpiral(const float &degreesOfArm, const size_t &numberOfSamples);

    Eigen::MatrixX2i getEigenEdges() const;
-    std::vector<vcg::Point2i> getEdges();
+    std::vector<vcg::Point2i> computeEdges();

-    Eigen::MatrixX3d getEigenVertices();
+    Eigen::MatrixX3d getEigenVertices() const;
    Eigen::MatrixX3d getEigenEdgeNormals() const;
    void printVertexCoordinates(const size_t &vi) const;
 #ifdef POLYSCOPE_DEFINED
@ -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;
--- a/reducedmodel.cpp
+++ b/reducedmodel.cpp
@ -0,0 +1,94 @@
+#include "reducedmodel.hpp"
+
+ReducedModel::ReducedModel() : PatternGeometry({1, 0, 0, 2, 1, 2, 1}, {{0, 3}})
+
+{
+    interfaceNodeIndex = 3;
+    label = "ReducedModel";
+    //    constexpr double initialHexSize = 0.3;
+    //    CoordType movableVertex_barycentric(1 - initialHexSize, initialHexSize / 2, initialHexSize / 2);
+    //    const vcg::Triangle3<double> &baseTriangle = getBaseTriangle();
+    //    vert[0].P() = baseTriangle.cP(0) * movableVertex_barycentric[0]
+    //                  + baseTriangle.cP(1) * movableVertex_barycentric[1]
+    //                  + baseTriangle.cP(2) * movableVertex_barycentric[2];
+}
+
+void ReducedModel::updateBaseTriangleGeometry_theta(const double &newTheta_deg)
+{
+    const CoordType baseTriangleHexagonVertexPosition = vert[0].cP();
+
+    const CoordType thetaRotatedHexagonBaseTriangleVertexPosition
+        = vcg::RotationMatrix(PatternGeometry::DefaultNormal, vcg::math::ToRad(newTheta_deg))
+          * baseTriangleHexagonVertexPosition;
+    vert[0].P() = thetaRotatedHexagonBaseTriangleVertexPosition;
+
+    //    constexpr int fanSize = 6;
+    //    for (int rotationCounter = 0; rotationCounter < /*ReducedModelOptimizer::*/ fanSize;
+    //         rotationCounter++) {
+    //        vert[2 * rotationCounter].P() = vcg::RotationMatrix(PatternGeometry::DefaultNormal,
+    //                                                            vcg::math::ToRad(60.0 * rotationCounter))
+    //                                        * thetaRotatedHexagonBaseTriangleVertexPosition;
+    //    }
+    updateEigenEdgeAndVertices();
+}
+
+void ReducedModel::updateBaseTriangleGeometry_R(const double &newR)
+{
+    const CoordType barycentricCoordinates_hexagonBaseTriangleVertex(1 - newR, newR / 2, newR / 2);
+    const vcg::Triangle3<double> &baseTriangle = getBaseTriangle();
+    const CoordType hexagonBaseTriangleVertexPosition
+        = baseTriangle.cP(0) * barycentricCoordinates_hexagonBaseTriangleVertex[0]
+          + baseTriangle.cP(1) * barycentricCoordinates_hexagonBaseTriangleVertex[1]
+          + baseTriangle.cP(2) * barycentricCoordinates_hexagonBaseTriangleVertex[2];
+    vert[0].P() = hexagonBaseTriangleVertexPosition;
+
+    //    constexpr int fanSize = 6;
+    //    for (int rotationCounter = 0; rotationCounter < /*ReducedModelOptimizer::*/ fanSize;
+    //         rotationCounter++) {
+    //        vert[2 * rotationCounter].P() = vcg::RotationMatrix(PatternGeometry::DefaultNormal,
+    //                                                            vcg::math::ToRad(60.0 * rotationCounter))
+    //                                        * hexagonBaseTriangleVertexPosition;
+    //        //        std::cout << vert[2 * rotationCounter].P()[0] << " " << vert[2 * rotationCounter].P()[1]
+    //        //                  << " " << vert[2 * rotationCounter].P()[2] << std::endl;
+    //    }
+    updateEigenEdgeAndVertices();
+}
+
+void ReducedModel::updateBaseTriangleGeometry_theta(
+    const double &newTheta_deg, std::shared_ptr<VCGEdgeMesh> &pReducedModel_edgeMesh)
+{
+    std::dynamic_pointer_cast<ReducedModel>(pReducedModel_edgeMesh)
+        ->updateBaseTriangleGeometry_theta(newTheta_deg);
+}
+
+void ReducedModel::updateBaseTriangleGeometry_R(const double &newR,
+                                                std::shared_ptr<VCGEdgeMesh> &pReducedModel_edgeMesh)
+{
+    std::dynamic_pointer_cast<ReducedModel>(pReducedModel_edgeMesh)
+        ->updateBaseTriangleGeometry_R(newR);
+}
+
+void ReducedModel::createFan(const size_t &fanSize)
+{
+    deleteDanglingVertices();
+    PatternGeometry rotatedPattern;
+    vcg::tri::Append<PatternGeometry, PatternGeometry>::MeshCopy(rotatedPattern, *this);
+    for (int rotationCounter = 1; rotationCounter < fanSize; rotationCounter++) {
+        vcg::Matrix44d R;
+        auto rotationAxis = PatternGeometry::DefaultNormal;
+        R.SetRotateDeg(360.0 / fanSize, rotationAxis);
+        vcg::tri::UpdatePosition<PatternGeometry>::Matrix(rotatedPattern, R);
+        vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(*this, rotatedPattern);
+        //For the reduced model we also need to connect to neighbouring base triangles of the fan
+        if (rotationCounter == fanSize - 1) {
+            vcg::tri::Allocator<PatternGeometry>::AddEdge(*this,
+                                                          0,
+                                                          this->VN() - rotatedPattern.VN());
+        }
+        vcg::tri::Allocator<PatternGeometry>::AddEdge(*this,
+                                                      this->VN() - 2 * rotatedPattern.VN(),
+                                                      this->VN() - rotatedPattern.VN());
+        removeDuplicateVertices();
+        updateEigenEdgeAndVertices();
+    }
+}
--- a/reducedmodel.hpp
+++ b/reducedmodel.hpp
@ -0,0 +1,20 @@
+#ifndef REDUCEDMODEL_HPP
+#define REDUCEDMODEL_HPP
+
+#include "trianglepatterngeometry.hpp"
+
+class ReducedModel : public PatternGeometry
+{
+
+public:
+    ReducedModel();
+    void updateBaseTriangleGeometry_theta(const double &newTheta_deg);
+    void updateBaseTriangleGeometry_R(const double &newR);
+    static void updateBaseTriangleGeometry_theta(
+        const double &newTheta_deg, std::shared_ptr<VCGEdgeMesh> &pReducedModel_edgeMesh);
+    static void updateBaseTriangleGeometry_R(const double &newR,
+                                             std::shared_ptr<VCGEdgeMesh> &pReducedModel_edgeMesh);
+    void createFan(const size_t &fanSize = 6) override;
+};
+
+#endif // REDUCEDMODEL_HPP
--- a/reducedmodelevaluator.cpp
+++ b/reducedmodelevaluator.cpp
@ -1,6 +1,9 @@
 #include "reducedmodelevaluator.hpp"
 #include "hexagonremesher.hpp"
+#include "reducedmodel.hpp"
+#include "reducedmodeloptimizer.hpp"
 #include "trianglepatterngeometry.hpp"
+#include <execution>
 #include <filesystem>

 using FullPatternVertexIndex = VertexIndex;
@ -8,89 +11,301 @@ 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 = [&]() {
+    pTileIntoSurface = [&]() {
+        std::istringstream inputStream_tileIntoTriSurface(tileIntoSurfaceFileContent);
        VCGTriMesh tileInto_triMesh;
-        const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filename);
+        const bool surfaceLoadSuccessfull = tileInto_triMesh.load(inputStream_tileIntoTriSurface);
+        tileInto_triMesh.setLabel("instantMeshes_plane_34");
        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);
+}
+
+//double ReducedModelEvaluator::evaluateOptimizationSettings(
+//    const ReducedModelOptimization::Settings &optimizationSettings,
+//    const std::vector<std::shared_ptr<PatternGeometry>> &pPatterns,
+//    std::vector<ReducedModelEvaluator::Results> &patternEvaluationResults)
+//{
+//    assert(!pPatterns.empty());
+//    double optimizationError = 0;
+//    auto start = std::chrono::high_resolution_clock::now();
+//    std::vector<double> averageNormalizedError(pPatterns.size(), 0);
+//    patternEvaluationResults.clear();
+//    patternEvaluationResults.resize(pPatterns.size());
+//    std::for_each(
+//        //        std::execution::par_unseq,
+//        pPatterns.begin(),
+//        pPatterns.end(),
+//        [&](const std::shared_ptr<ConstPatternGeometry> &pPattern) {
+//            //            std::cout << "Optimizing " << pPattern->getLabel() << std::endl;
+//            ReducedModelOptimization::Results optimizationResults;
+//            ReducedModelOptimizer optimizer;
+//            optimizer.optimize(*pPattern, optimizationSettings, optimizationResults);
+//            const auto evaluationResults
+//                = ReducedModelEvaluator::evaluateReducedModel(optimizationResults,
+//                                                              tileIntoSurfaceFilePath,
+//                                                              scenariosDirPath,
+//                                                              fullPatternTessellatedResultsDirPath);
+//            const double averageNormalizedErrorOfPattern
+//                = std::reduce(evaluationResults.distances_normalizedDrm2reduced.begin(),
+//                              evaluationResults.distances_normalizedDrm2reduced.end())
+//                  / evaluationResults.distances_normalizedDrm2reduced.size();
+//            const int patternIndex = &pPattern - &patterns[0];
+//            averageNormalizedError[patternIndex] = averageNormalizedErrorOfPattern;
+//            patternsEvaluationResults[patternIndex] = evaluationResults;
+//        });
+//    const double strategyAverageNormalizedError = std::reduce(std::execution::par_unseq,
+//                                                              averageNormalizedError.begin(),
+//                                                              averageNormalizedError.end())
+//                                                  / pPointers.size();
+//    const auto totalDuration_min = std::chrono::duration_cast<std::chrono::seconds>(
+//                                       std::chrono::high_resolution_clock::now() - start)
+//                                       .count()
+//                                   / 60.0;
+//    std::cout << "Optimized pattern(s) in:" << totalDuration_min << " minutes." << std::endl;
+//    std::cout << "Average time per pattern:" << totalDuration_min / patternsPointers.size()
+//              << " minutes." << std::endl;
+//    std::cout << "Objective value:" << strategyAverageNormalizedError << std::endl;
+//    return strategyAverageNormalizedError;
+//    //    std::cout << "After:" << ++numberOfOptimizationRoundsExecuted << " iterations." << std::endl;
+//}
+
+ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
+    ReducedModelOptimization::Results &optimizationResult)
+{
+    const std::filesystem::path scenariosDirectoryPath
+        = "/home/iason/Coding/Projects/Approximating shapes with flat "
+          "patterns/ReducedModelEvaluator/Scenarios";
+    const std::filesystem::path fullPatternTessellatedResultsDirectoryPath
+        = "/home/iason/Coding/Projects/Approximating shapes with flat "
+          "patterns/ReducedModelEvaluator/TessellatedResults";
+
+    return evaluateReducedModel(optimizationResult,
+                                scenariosDirectoryPath,
+                                fullPatternTessellatedResultsDirectoryPath);
+}
+//void ReducedModelEvaluator::printResults(const Results &evaluationResults,
+//                                         const std::string &resultsLabel,
+//                                         const std::filesystem::path &resultsOutputPath)
+//{
+//    const bool outputPathIsDirectory = !resultsOutputPath.empty()
+//                                       && !resultsOutputPath.has_extension();
+//    const bool outputPathIsFile = !resultsOutputPath.empty() && resultsOutputPath.has_extension();
+//    assert(outputPathIsDirectory && outputPathIsFile);
+//    if (outputPathIsDirectory) {
+//        std::filesystem::create_directories(resultsOutputPath);
+//    }
+
+//#else
+//    std::filesystem::path csvOutputFilePath;
+//    bool shouldOverwrite = false;
+//    if (outputPathIsDirectory) {
+//        csvOutputFilePath = std::filesystem::path(resultsOutputPath)
+//                                .append("distances_" + resultsLabel + ".csv")
+//                                .string();
+//        shouldOverwrite = true;
+//    } else if (outputPathIsFile) {
+//        csvOutputFilePath = resultsOutputPath;
+//    }
+
+//    csvFile csvOutput(csvOutputFilePath, shouldOverwrite);
+//    printResults(evaluationResults, resultsLabel, csvOutput);
+//}
+
+void ReducedModelEvaluator::printResults(const Results &evaluationResults,
+                                         const std::string &resultsLabel)
+{
+    csvFile csvOutputToCout({}, true);
+    Settings exportSettings;
+    exportSettings.exportingDirection = Vertical;
+    exportSettings.shouldWriteHeader = false;
+    exportSettings.resultsLabel = resultsLabel;
+    printResults(evaluationResults, exportSettings, csvOutputToCout);
+}
+
+void ReducedModelEvaluator::printResults(const Results &evaluationResults,
+                                         const Settings &settings,
+                                         csvFile &csvOutput)
+{
+    if (settings.shouldWriteHeader) {
+        csvOutput << csvExportingDataStrings[settings.exportingData];
+        printHeader(settings, csvOutput);
+        //    csvOutput << "Average error";
+        //    csvOutput<<"Cumulative error";
+        csvOutput << endrow;
+    }
+    if (!settings.resultsLabel.empty()) {
+        csvOutput << settings.resultsLabel;
+    }
+    if (settings.exportingDirection == Vertical) {
+        printResultsVertically(evaluationResults, csvOutput);
+    } else {
+        printResultsHorizontally(evaluationResults, csvOutput);
+    }
+}
+
+void ReducedModelEvaluator::printHeader(const Settings &settings, csvFile &csvOutput)
+{
+    if (settings.exportingDirection == Horizontal) {
+        //        csvOutput << "Job label";
+        for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
+             jobIndex++) {
+            const std::string &jobLabel = ReducedModelEvaluator::scenariosTestSetLabels[jobIndex];
+            csvOutput << jobLabel;
+        }
+    } else {
+        std::cout << "Vertical header not defined" << std::endl;
+        assert(false);
+        std::terminate();
+    }
+}
+
+void ReducedModelEvaluator::printResultsHorizontally(const Results &evaluationResults,
+                                                     csvFile &csvOutput)
+{
+    //print header
+    //print raw error
+    constexpr bool shouldPrintRawError = false;
+    if (shouldPrintRawError) {
+        //        csvOutput << "drm2Reduced";
+        double sumOfFull2Reduced = 0;
+        int numOfNonEvaluatedScenarios = 0;
+        for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
+             jobIndex++) {
+            const double &distance_fullDrmToReduced = evaluationResults
+                                                          .distances_drm2reduced[jobIndex];
+            if (distance_fullDrmToReduced == -1) {
+                csvOutput << "notEvaluated";
+                numOfNonEvaluatedScenarios++;
+            } else {
+                csvOutput << distance_fullDrmToReduced;
+                sumOfFull2Reduced += distance_fullDrmToReduced;
+            }
+        }
+        //    const int numOfEvaluatedScenarios = ReducedModelEvaluator::scenariosTestSetLabels.size()
+        //                                        - numOfNonEvaluatedScenarios;
+        //    const double averageDistance_full2Reduced = sumOfFull2Reduced / numOfEvaluatedScenarios;
+        //    csvOutput << averageDistance_full2Reduced;
+        //        csvOutput << endrow;
+    }
+    //print normalized error
+    //    csvOutput << "norm_drm2Reduced";
+    double sumOfNormalizedFull2Reduced = 0;
+    for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
+         jobIndex++) {
+        const double &distance_normalizedFullDrmToReduced
+            = evaluationResults.distances_normalizedDrm2reduced[jobIndex];
+        if (distance_normalizedFullDrmToReduced == -1) {
+            csvOutput << "notEvaluated";
+        } else {
+            csvOutput << distance_normalizedFullDrmToReduced;
+            sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
+        }
+    }
+    //    const double averageDistance_normalizedFull2Reduced = sumOfNormalizedFull2Reduced
+    //                                                          / numOfEvaluatedScenarios;
+    //    csvOutput << averageDistance_normalizedFull2Reduced;
+    //    csvOutput << endrow;
+}
+
+void ReducedModelEvaluator::printResultsVertically(const Results &evaluationResults,
+                                                   csvFile &csvOutput)
+{
+#ifdef POLYSCOPE_DEFINED
+    csvOutput << "drm2Reduced"
+              << "norm_drm2Reduced";
+#else
+    csvOutput << "Job Label"
+              << "drm2Reduced"
+              << "norm_drm2Reduced";
+
+#endif
+    csvOutput << endrow;
+    double sumOfFull2Reduced = 0;
+    double sumOfNormalizedFull2Reduced = 0;
+    int numOfNonEvaluatedScenarios = 0;
+    for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
+         jobIndex++) {
+        const double &distance_fullDrmToReduced = evaluationResults.distances_drm2reduced[jobIndex];
+        const double &distance_normalizedFullDrmToReduced
+            = evaluationResults.distances_normalizedDrm2reduced[jobIndex];
+#ifndef POLYSCOPE_DEFINED
+        const std::string &jobLabel = ReducedModelEvaluator::scenariosTestSetLabels[jobIndex];
+        csvOutput << jobLabel;
+#endif
+        if (distance_fullDrmToReduced == -1) {
+            csvOutput << "notEvaluated"
+                      << "notEvaluated";
+            numOfNonEvaluatedScenarios++;
+        } else {
+            csvOutput << distance_fullDrmToReduced << distance_normalizedFullDrmToReduced;
+            sumOfFull2Reduced += distance_fullDrmToReduced;
+            sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
+        }
+        csvOutput << endrow;
+        //        sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
+    }
+    const int numOfEvaluatedScenarios = ReducedModelEvaluator::scenariosTestSetLabels.size()
+                                        - numOfNonEvaluatedScenarios;
+    const double averageDistance_full2Reduced = sumOfFull2Reduced / numOfEvaluatedScenarios;
+    const double averageDistance_normalizedFull2Reduced = sumOfNormalizedFull2Reduced
+                                                          / numOfEvaluatedScenarios;
+#ifndef POLYSCOPE_DEFINED
+    csvOutput << "Average error";
+#endif
+    csvOutput << averageDistance_full2Reduced << averageDistance_normalizedFull2Reduced;
+    csvOutput << endrow;
+#ifndef POLYSCOPE_DEFINED
+    csvOutput << "Cumulative error";
+#endif
+    csvOutput << sumOfFull2Reduced << sumOfNormalizedFull2Reduced;
+    csvOutput << endrow;
+    csvOutput << endrow;
+}
+
+//void ReducedModelEvaluator::createFullAndReducedPatternTessellations(){
+
+//}
+
+//ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
+//    std::vector<ReducedModelOptimization::Results> &optimizationResults,
+//    const std::filesystem::path &tileInto_triMesh_filePath,
+//    const std::filesystem::path &scenariosDirectoryPath,
+//    //    const std::filesystem::path &reducedPatternFilePath,
+//    const std::filesystem::path &fullPatternTessellatedResultsDirectoryPath)
+//{
+//    //Load surface
+//    std::shared_ptr<VCGPolyMesh> pTileIntoSurface = [&]() {
+//        VCGTriMesh tileInto_triMesh;
+//        const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filePath);
+//        assert(surfaceLoadSuccessfull);
+//        return PolygonalRemeshing::remeshWithPolygons(tileInto_triMesh);
+//    }();
+//    const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResult.baseTriangle.cP(0),
+//                                                             (optimizationResult.baseTriangle.cP(1)
+//                                                              + optimizationResult.baseTriangle.cP(
+//                                                                  2))
+//                                                                 / 2);
+//    pTileIntoSurface->moveToCenter();
+//    const double scaleFactor = optimizedBaseTriangleHeight
+//                               / pTileIntoSurface->getAverageFaceRadius();
+//    vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
+//    //Tile full pattern into surface
+
+//}
+
+ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
+    ReducedModelOptimization::Results &optimizationResult,
+    const std::filesystem::path &scenariosDirectoryPath,
+    const std::filesystem::path &fullPatternTessellatedResultsDirectoryPath)
+{
+    //    const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResult.baseTriangle.cP(0),
+    //                                                             (optimizationResult.baseTriangle.cP(1)
+    //                                                              + optimizationResult.baseTriangle.cP(
+    //                                                                  2))
+    //                                                                 / 2);
    pTileIntoSurface->moveToCenter();
-    const double scaleFactor = optimizedBaseTriangleHeight
+    const double scaleFactor = optimizationResult.settings.targetBaseTriangleSize
                               / pTileIntoSurface->getAverageFaceRadius();
    vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
    //    tileIntoSurface.registerForDrawing();
@ -98,7 +313,7 @@ std::vector<double> ReducedModelEvaluator::evaluateReducedModel(

    //Tile full pattern into surface
    std::vector<PatternGeometry> fullPatterns(1);
-    fullPatterns[0].copy(optimizationResults.baseTriangleFullPattern);
+    fullPatterns[0].copy(optimizationResult.baseTriangleFullPattern);
    ////        Base triangle pattern might contain dangling vertices.Remove those
    fullPatterns[0].interfaceNodeIndex = 3;
    fullPatterns[0].deleteDanglingVertices();
@ -115,15 +330,18 @@ std::vector<double> ReducedModelEvaluator::evaluateReducedModel(
    pTiledFullPattern->setLabel("Tiled_full_patterns");
    //    pTiledFullPattern->registerForDrawing();
    //Tile reduced pattern into surface
-    PatternGeometry reducedPattern;
-    reducedPattern.load(reducedPatternFilePath);
-    reducedPattern.deleteDanglingVertices();
-    assert(reducedPattern.interfaceNodeIndex == 1);
+    //    PatternGeometry reducedPattern;
+    ReducedModel reducedModel;
+    //    reducedModel.registerForDrawing();
+    //    polyscope::show();
+    reducedModel.deleteDanglingVertices();
+    //    reducedPattern.interfaceNodeIndex = 1;
+    //    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]);
+    reducedPatterns[0].copy(reducedModel);
+    const auto reducedPatternBaseTriangle = reducedModel.computeBaseTriangle();
+    ReducedModelOptimization::Results::applyOptimizationResults_reducedModel_nonFanned(
+        optimizationResult, reducedPatternBaseTriangle, reducedPatterns[0]);

    std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
    std::vector<size_t> tileIntoEdgeToTiledReducedVi;
@ -135,7 +353,8 @@ std::vector<double> ReducedModelEvaluator::evaluateReducedModel(
                                       tileIntoEdgeToTiledReducedVi,
                                       perPatternIndexTiledReducedPatternEdgeIndices);
    pTiledReducedPattern->setLabel("Tiled_reduced_patterns");
-    //        pTiledReducedPattern->registerForDrawing();
+    //    pTiledReducedPattern->registerForDrawing();
+    //    polyscope::show();

    std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
        fullToReducedViMap; //of only the common vertices
@ -160,12 +379,12 @@ std::vector<double> ReducedModelEvaluator::evaluateReducedModel(
    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) {
+    if (optimizationResult.fullPatternYoungsModulus == 0) {
        std::cerr << "Full pattern's young modulus not found." << std::endl;
        std::terminate();
    }
    pTiledFullPattern_simulationMesh->setBeamMaterial(0.3,
-                                                      optimizationResults.fullPatternYoungsModulus);
+                                                      optimizationResult.fullPatternYoungsModulus);
    pTiledFullPattern_simulationMesh->reset();

    ////Tessellated reduced pattern simulation mesh
@ -174,144 +393,179 @@ std::vector<double> ReducedModelEvaluator::evaluateReducedModel(
    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();
+        optimizationResult, pTiledReducedPattern_simulationMesh);
+    //    pTiledReducedPattern_simulationMesh->reset();
+    Results evaluationResults;
+    evaluationResults.distances_drm2reduced.fill(-1);
+    evaluationResults.distances_normalizedDrm2reduced.fill(-1);
+    DRMSimulationModel::Settings drmSimulationSettings;
+    drmSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
+    //    drmSimulationSettings.load(drmSettingsFilePath);
+    drmSimulationSettings.beVerbose = true;
+    drmSimulationSettings.maxDRMIterations = 5e6;
+    drmSimulationSettings.debugModeStep = 100000;
+    drmSimulationSettings.translationalKineticEnergyThreshold = 1e-15;
+    drmSimulationSettings.linearGuessForceScaleFactor = 0.8;
+    drmSimulationSettings.viscousDampingFactor = 7e-3;
+    drmSimulationSettings.xi = 0.9999;
+    //    drmSimulationSettings.Dtini = 5.86;
+    drmSimulationSettings.gamma = 0.25;
+#ifdef POLYSCOPE_DEFINED
+    //    drmSimulationSettings.shouldDraw = true;
+    drmSimulationSettings.shouldCreatePlots = true;
+#endif
+    constexpr bool shouldRerunFullPatternSimulation = false;
+    //    for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
+    std::for_each(
+        //#ifndef POLYSCOPE_DEFINED
+        //        std::execution::par_unseq,
+        //#endif
+        scenariosTestSetLabels.begin(),
+        scenariosTestSetLabels.end(),
+        [&](const std::string &jobLabel) {
+            //check if reduced model scenario exists
+            //        const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
+            const std::filesystem::path tiledReducedPatternJobFilePath
+                = std::filesystem::path(scenariosDirectoryPath)
+                      .append(pTileIntoSurface->getLabel())
+                      .append(jobLabel)
+                      .append("ReducedJob")
+                      .append(SimulationJob::jsonDefaultFileName);
+            if (!std::filesystem::exists(tiledReducedPatternJobFilePath)) {
+                std::cerr << "Scenario " << jobLabel
+                          << " not found in:" << tiledReducedPatternJobFilePath << std::endl;
+                //            continue; //if not move on to the next scenario
+                return;
+            }
+            //Map the reduced job to the job on the pattern tessellation
+            //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();
+            const std::filesystem::path surfaceFolderPath
+                = std::filesystem::path(fullPatternTessellatedResultsDirectoryPath)
+                      .append(pTileIntoSurface->getLabel());
+            const std::string scenarioLabel = pJob_tiledFullPattern->getLabel();
+            const std::filesystem::path scenarioDirectoryPath
+                = std::filesystem::path(surfaceFolderPath).append(scenarioLabel);
+            //Save reduced job
+            constexpr bool exportReducedJob = false;
+            if (exportReducedJob) {
+                const std::filesystem::path reducedJobDirectoryPath
+                    = std::filesystem::path(scenarioDirectoryPath).append("ReducedJob");
+                std::filesystem::create_directories(reducedJobDirectoryPath);
+                pJob_tiledReducedPattern->save(reducedJobDirectoryPath);
+            }
+            //Check if the drm simulation of the full pattern has already been computed
+            ////Full
+            const std::string &patternLabel = [&]() {
+                const std::string patternLabel = optimizationResult.baseTriangleFullPattern
+                                                     .getLabel();
+                if (patternLabel.find("_") == std::string::npos) {
+                    return std::to_string(optimizationResult.baseTriangleFullPattern.EN()) + "_"
+                           + patternLabel;
+                } else {
+                    return patternLabel;
+                }
+            }();
+            const auto fullResultsFolderPath = std::filesystem::path(scenarioDirectoryPath)
+                                                   .append(patternLabel)
+                                                   .append("Results");
+            if (shouldRerunFullPatternSimulation && std::filesystem::exists(fullResultsFolderPath)) {
+                std::filesystem::remove_all(fullResultsFolderPath);
            }
-        }();
-        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;
-        }
+            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);
+                //#ifdef POLYSCOPE_DEFINED
+                //            std::array<double, 3> resultsColor({28.0, 99.0, 227.0});
+                //            simulationResults_tiledFullPattern_drm.registerForDrawing(resultsColor);
+                //            std::ifstream ifs("CameraSettings.json");
+                //            nlohmann::json json;
+                //            ifs >> json;
+                //            polyscope::view::setCameraFromJson(json.dump(), false);
+                //            //            polyscope::show();
+                //            const std::string cameraJson = polyscope::view::getCameraJson();
+                //            std::filesystem::path jsonFilePath("CameraSettings.json");
+                //            std::ofstream jsonFile_cameraSettings(jsonFilePath.string());
+                //            jsonFile_cameraSettings << cameraJson;
+                //            jsonFile_cameraSettings.close();
+                //            std::filesystem::create_directories("screenshots");
+                //            const std::string screenshotOutputFilePath
+                //                = (std::filesystem::current_path()
+                //                       .append("screenshots")
+                //                       .append(optimizationResult.label + "_" + pJob_tiledFullPattern->getLabel()))
+                //                      .string()
+                //                  + ".png";
+                //            //            std::cout << "Saving image to:" << screenshotOutputFilePath << std::endl;
+                //            polyscope::screenshot(screenshotOutputFilePath, false);
+                //            simulationResults_tiledFullPattern_drm.unregister();
+                //#endif

-        LinearSimulationModel linearSimulationModel;
-        SimulationResults simulationResults_tiledReducedPattern
-            = linearSimulationModel.executeSimulation(pJob_tiledReducedPattern);
+                simulationResults_tiledFullPattern_drm.converged = true;
+            } else {
+                std::cout << "Drm results not found in:" << fullResultsFolderPath << std::endl;
+                //Full
+                std::cout << "Executing:" << jobLabel << std::endl;
+                DRMSimulationModel drmSimulationModel;
+                simulationResults_tiledFullPattern_drm
+                    = drmSimulationModel.executeSimulation(pJob_tiledFullPattern,
+                                                           drmSimulationSettings);
+                simulationResults_tiledFullPattern_drm.setLabelPrefix("DRM");
+            }
+            if (!simulationResults_tiledFullPattern_drm.converged) {
+                std::cerr << "Full pattern simulation failed." << std::endl;
+                std::cerr << "Not saving results" << std::endl;
+                //            continue;
+                return;
+            }
+            std::filesystem::create_directories(fullResultsFolderPath);
+            const std::filesystem::path drmResultsOutputPath
+                = std::filesystem::path(scenarioDirectoryPath).append(patternLabel);
+            simulationResults_tiledFullPattern_drm.save(drmResultsOutputPath);

-        //        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
+            LinearSimulationModel linearSimulationModel;
+            SimulationResults simulationResults_tiledReducedPattern
+                = linearSimulationModel.executeSimulation(pJob_tiledReducedPattern);

-    //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;
-    }
+            //            simulationResults_tiledReducedPattern.registerForDrawing();
+            //            simulationResults_tiledFullPattern_drm.registerForDrawing();
+            //            polyscope::show();

-    std::cout << "Average normalized error per scenario:"
-              << sumOfNormalizedFull2Reduced / scenariosTestSetLabels.size() << std::endl;
-    return distances_normalizedDrm2reduced;
-    //#endif
+            //compute the full2reduced 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;
+            const int jobIndex = &jobLabel - &scenariosTestSetLabels[0];
+            evaluationResults.distances_drm2reduced[jobIndex] = distance_fullDrmToReduced;
+            evaluationResults.distances_normalizedDrm2reduced[jobIndex]
+                = distance_normalizedFullDrmToReduced;
+        });
+
+    return evaluationResults;
 }
--- a/reducedmodelevaluator.hpp
+++ b/reducedmodelevaluator.hpp
@ -6,9 +6,212 @@
 class ReducedModelEvaluator
 {
 public:
+    enum CSVExportingDirection { Vertical = 0, Horizontal };
+    enum CSVExportingData {
+        raw_drm2Reduced = 0,
+        norm_drm2Reduced,
+        raw_and_norm_drm2Reduced,
+        NumberOfDataTypes
+    };
+    inline static std::array<std::string, NumberOfDataTypes>
+        csvExportingDataStrings{"raw_drm2Reduced", "norm_drm2Reduced", "raw_and_norm_drm2Reduced"};
+    struct Settings
+    {
+        CSVExportingDirection exportingDirection{Horizontal};
+        CSVExportingData exportingData{norm_drm2Reduced};
+        bool shouldWriteHeader{true};
+        std::string resultsLabel;
+    };
+
+    inline static constexpr int NumberOfEvaluationScenarios{22};
+    struct Results
+    {
+        std::array<double, NumberOfEvaluationScenarios> distances_drm2reduced;
+        std::array<double, NumberOfEvaluationScenarios> distances_normalizedDrm2reduced;
+        std::array<std::string, NumberOfEvaluationScenarios> evaluationScenarioLabels;
+    };
    ReducedModelEvaluator();
-    static std::vector<double> evaluateReducedModel(
-        ReducedModelOptimization::Results &optimizationResults);
+    Results evaluateReducedModel(
+        ReducedModelOptimization::Results &optimizationResult,
+        const std::filesystem::path &scenariosDirectoryPath,
+        //        const std::filesystem::path &reducedPatternFilePath,
+        const std::filesystem::path &fullPatternTessellatedResultsDirectoryPath);
+    Results evaluateReducedModel(ReducedModelOptimization::Results &optimizationResult);
+    static void printResultsVertically(const ReducedModelEvaluator::Results &evaluationResults,
+                                       csvFile &csvOutput);
+    static void printResults(const ReducedModelEvaluator::Results &evaluationResults,
+                             const std::string &resultsLabel);
+
+    inline static std::array<std::string, NumberOfEvaluationScenarios> 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",
+        //#ifdef POLYSCOPE_DEFINED
+        //                               "22Hex_shear_2N",
+        //                               "22Hex_shear_5N",
+        //                               "22Hex_axial_10N",
+        //                               "22Hex_axial_20N",
+        //#else
+        //                                       "notUsed_22Hex_shear_2N",
+        //                                       "notUsed_22Hex_shear_5N",
+        //                                       "notUsed_22Hex_axial_10N",
+        //                                       "notUsed_22Hex_axial_20N",
+        //#endif
+        "22Hex_cylinder_0.05N",
+        "22Hex_cylinder_0.1N",
+        "22Hex_s_0.05N",
+        //                               "22Hex_s_0.1N"
+    };
+    static void printResultsHorizontally(const Results &evaluationResults, csvFile &csvOutput);
+    static void printResults(const Results &evaluationResults,
+                             const Settings &settings,
+                             csvFile &csvOutput);
+    static void printHeader(const Settings &settings, csvFile &csvOutput);
+    //    static double evaluateOptimizationSettings(
+    //        const ReducedModelOptimization::Settings &optimizationSettings,
+    //        const std::vector<std::shared_ptr<PatternGeometry>> &pPatterns,
+    //        std::vector<ReducedModelEvaluator::Results> &patternEvaluationResults);
+    std::shared_ptr<VCGPolyMesh> pTileIntoSurface;
+    inline static constexpr char *tileIntoSurfaceFileContent = R"~(OFF
+46 66 0
+-0.0745923 0.03573945 0 
+-0.07464622 0.02191801 0 
+-0.06264956 0.02878203 0 
+-0.08658896 0.02887542 0 
+-0.06270347 0.01496059 0 
+-0.06259564 0.04260346 0 
+-0.08664289 0.01505398 0 
+-0.08653505 0.04269686 0 
+-0.0507068 0.02182462 0 
+-0.07453838 0.04956088 0 
+-0.05065288 0.03564605 0 
+-0.09858564 0.02201139 0 
+-0.09853172 0.03583283 0 
+-0.06254172 0.0564249 0 
+-0.05059896 0.04946748 0 
+-0.08648112 0.0565183 0 
+-0.09847781 0.04965428 0 
+-0.03871013 0.02868864 0 
+-0.07448447 0.06338232 0 
+-0.03865621 0.04251007 0 
+-0.1105284 0.02896881 0 
+-0.1104745 0.04279025 0 
+-0.03876406 0.0148672 0 
+-0.05054504 0.06328891 0 
+-0.0624878 0.07024634 0 
+-0.03860229 0.0563315 0 
+-0.1105823 0.01514738 0 
+-0.08642721 0.07033974 0 
+-0.09842389 0.06347572 0 
+-0.1104206 0.05661169 0 
+-0.07443054 0.07720376 0 
+-0.05049112 0.07711035 0 
+-0.03854837 0.07015293 0 
+-0.06243387 0.08406778 0 
+-0.08637329 0.08416118 0 
+-0.09836997 0.07729716 0 
+-0.1103666 0.07043312 0 
+-0.07437663 0.09102518 0 
+-0.03849445 0.08397438 0 
+-0.0504372 0.0909318 0 
+-0.06237995 0.09788923 0 
+-0.08631937 0.09798261 0 
+-0.09831604 0.09111859 0 
+-0.1103127 0.08425456 0 
+-0.03844052 0.09779582 0 
+-0.1102588 0.09807601 0 
+3 0 1 2 
+3 3 1 0 
+3 4 2 1 
+3 5 0 2 
+3 3 6 1 
+3 3 0 7 
+3 8 2 4 
+3 5 9 0 
+3 10 5 2 
+3 3 11 6 
+3 7 0 9 
+3 12 3 7 
+3 10 2 8 
+3 5 13 9 
+3 10 14 5 
+3 12 11 3 
+3 7 9 15 
+3 12 7 16 
+3 10 8 17 
+3 14 13 5 
+3 18 9 13 
+3 19 14 10 
+3 12 20 11 
+3 18 15 9 
+3 16 7 15 
+3 21 12 16 
+3 22 17 8 
+3 19 10 17 
+3 23 13 14 
+3 18 13 24 
+3 19 25 14 
+3 21 20 12 
+3 26 11 20 
+3 18 27 15 
+3 16 15 28 
+3 21 16 29 
+3 23 24 13 
+3 23 14 25 
+3 30 18 24 
+3 30 27 18 
+3 15 27 28 
+3 29 16 28 
+3 23 31 24 
+3 23 25 32 
+3 30 24 33 
+3 30 34 27 
+3 35 28 27 
+3 29 28 36 
+3 23 32 31 
+3 24 31 33 
+3 30 33 37 
+3 30 37 34 
+3 35 27 34 
+3 35 36 28 
+3 32 38 31 
+3 33 31 39 
+3 40 37 33 
+3 34 37 41 
+3 35 34 42 
+3 35 43 36 
+3 38 39 31 
+3 40 33 39 
+3 34 41 42 
+3 35 42 43 
+3 44 39 38 
+3 45 43 42 
+)~";
 };

 #endif // REDUCEDMODELEVALUATOR_HPP
--- a/reducedmodeloptimizer.cpp
+++ b/reducedmodeloptimizer.cpp
--- a/reducedmodeloptimizer.hpp
+++ b/reducedmodeloptimizer.hpp
@ -0,0 +1,362 @@
+#ifndef REDUCEDMODELOPTIMIZER_HPP
+#define REDUCEDMODELOPTIMIZER_HPP
+
+#include "csvfile.hpp"
+#include "drmsimulationmodel.hpp"
+#include "edgemesh.hpp"
+#include "linearsimulationmodel.hpp"
+#ifdef POLYSCOPE_DEFINED
+#include "matplot/matplot.h"
+#endif
+#include "reducedmodel.hpp"
+#include "reducedmodeloptimizer_structs.hpp"
+#include "simulationmesh.hpp"
+#include <Eigen/Dense>
+#ifdef DLIB_DEFINED
+#include <dlib/global_optimization.h>
+#include <dlib/optimization.h>
+#endif
+
+#ifdef POLYSCOPE_DEFINED
+#include "polyscope/color_management.h"
+#endif // POLYSCOPE_DEFINED
+using FullPatternVertexIndex = VertexIndex;
+using ReducedPatternVertexIndex = VertexIndex;
+
+class ReducedModelOptimizer
+{
+public:
+    struct OptimizationState
+    {
+        std::vector<SimulationResults> fullPatternResults;
+        std::vector<double> translationalDisplacementNormalizationValues;
+        std::vector<double> rotationalDisplacementNormalizationValues;
+        std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
+        std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
+        std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+            reducedToFullInterfaceViMap;
+        std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            fullPatternOppositeInterfaceViPairs;
+        matplot::line_handle gPlotHandle;
+        std::vector<size_t> objectiveValueHistory_iteration;
+        std::vector<double> objectiveValueHistory;
+        std::vector<double> plotColors;
+        std::array<double,
+                   ReducedModelOptimization::OptimizationParameterIndex::NumberOfOptimizationVariables>
+            parametersInitialValue;
+        std::array<double,
+                   ReducedModelOptimization::OptimizationParameterIndex::NumberOfOptimizationVariables>
+            optimizationInitialValue;
+        std::vector<int> simulationScenarioIndices;
+        double minY{DBL_MAX};
+        std::vector<double> minX;
+        int numOfSimulationCrashes{false};
+        int numberOfFunctionCalls{0};
+        //Variables for finding the full pattern simulation forces
+        std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh;
+        std::array<std::function<void(const double &newValue,
+                                      std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>,
+                   7>
+            functions_updateReducedPatternParameter;
+        std::vector<double> xMin;
+        std::vector<double> xMax;
+        std::vector<double> scenarioWeights;
+        std::vector<ReducedModelOptimization::Settings::ObjectiveWeights> objectiveWeights;
+    };
+
+private:
+    OptimizationState optimizationState;
+    vcg::Triangle3<double> baseTriangle;
+    std::function<void(const double &,
+                       const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>> &,
+                       SimulationJob &)>
+        constructScenarioFunction;
+    std::shared_ptr<SimulationMesh> m_pReducedPatternSimulationMesh;
+    std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh;
+    std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
+        m_fullToReducedInterfaceViMap;
+    std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+        m_fullPatternOppositeInterfaceViPairs;
+    std::unordered_map<size_t, size_t> nodeToSlot;
+    std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
+    std::string optimizationNotes;
+    std::array<std::function<void(
+                   const double &,
+                   const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>> &,
+                   SimulationJob &)>,
+               ReducedModelOptimization::NumberOfBaseSimulationScenarios>
+        constructBaseScenarioFunctions;
+    std::vector<bool> scenarioIsSymmetrical;
+    int fullPatternNumberOfEdges;
+    constexpr static double youngsModulus{1 * 1e9};
+    std::string fullPatternLabel;
+    //    ReducedModelOptimization::Settings optimizationSettings;
+
+public:
+    struct FunctionEvaluation
+    {
+        FunctionEvaluation() = default;
+        FunctionEvaluation(const std::vector<double> &x, double y) : x(x), y(y) {}
+
+        std::vector<double> x;
+        double y = std::numeric_limits<double>::quiet_NaN();
+    };
+
+    //    struct ParameterLabels
+    //    {
+    //        inline const static std::string E = {"E"};
+    //        inline const static std::string A = {"A"};
+    //        inline const static std::string I2 ={"I2"};
+    //        inline const static std::string I3 ={"I3"};
+    //        inline const static std::string J = {"J"};
+    //        inline const static std::string th= {"Theta"};
+    //        inline const static std::string R = {"R"};
+    //    };
+    //    inline constexpr static ParameterLabels parameterLabels();
+
+    inline static std::array<std::string, ReducedModelOptimization::NumberOfOptimizationVariables>
+        parameterLabels = {"E", "A", "I2", "I3", "J", "Theta", "R"};
+    constexpr static std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
+        simulationScenariosResolution = {12, 12, 22, 22, 22, 22};
+    constexpr static std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
+        baseScenarioWeights = {1, 1, 2, 3, 2};
+    inline static int totalNumberOfSimulationScenarios
+        = std::accumulate(simulationScenariosResolution.begin(),
+                          simulationScenariosResolution.end(),
+                          0);
+    inline static int fanSize{6};
+    inline static double initialHexagonSize{0.3};
+    int interfaceNodeIndex;
+    double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;
+
+    ReducedModelOptimizer();
+    static void computeReducedModelSimulationJob(
+        const SimulationJob &simulationJobOfFullModel,
+        const std::unordered_map<size_t, size_t> &fullToReducedMap,
+        SimulationJob &simulationJobOfReducedModel);
+
+    SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
+
+    static void runSimulation(const std::string &filename, std::vector<double> &x);
+
+    static std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationJobs(
+        const std::shared_ptr<SimulationMesh> &pMesh,
+        const std::unordered_map<size_t, size_t> &fullPatternOppositeInterfaceViMap);
+
+    static void createSimulationMeshes(
+        PatternGeometry &fullModel,
+        PatternGeometry &reducedModel,
+        std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
+        std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh);
+    void computeMaps(const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
+                     PatternGeometry &fullPattern,
+                     ReducedModel &reducedPattern,
+                     std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+                         &reducedToFullInterfaceViMap,
+                     std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
+                         &fullToReducedInterfaceViMap,
+                     std::vector<std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex>>
+                         &fullPatternOppositeInterfaceViMap);
+    static void visualizeResults(
+        const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs,
+        const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs,
+        const std::vector<ReducedModelOptimization::BaseSimulationScenario> &simulationScenarios,
+        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+            &reducedToFullInterfaceViMap);
+    static void registerResultsForDrawing(
+        const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
+        const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
+        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+            &reducedToFullInterfaceViMap);
+
+    static double computeRawTranslationalError(
+        const std::vector<Vector6d> &fullPatternDisplacements,
+        const std::vector<Vector6d> &reducedPatternDisplacements,
+        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+            &reducedToFullInterfaceViMap);
+
+    static double computeDisplacementError(
+        const std::vector<Vector6d> &fullPatternDisplacements,
+        const std::vector<Vector6d> &reducedPatternDisplacements,
+        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+            &reducedToFullInterfaceViMap,
+        const double &normalizationFactor);
+
+    static double computeRawRotationalError(
+        const std::vector<Eigen::Quaterniond> &rotatedQuaternion_fullPattern,
+        const std::vector<Eigen::Quaterniond> &rotatedQuaternion_reducedPattern,
+        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+            &reducedToFullInterfaceViMap);
+
+    static double computeRotationalError(const std::vector<Eigen::Quaterniond> &rotatedQuaternion_fullPattern,
+        const std::vector<Eigen::Quaterniond> &rotatedQuaternion_reducedPattern,
+        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+            &reducedToFullInterfaceViMap,
+        const double &normalizationFactor);
+    static double computeError(
+        const SimulationResults &simulationResults_fullPattern,
+        const SimulationResults &simulationResults_reducedPattern,
+        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+            &reducedToFullInterfaceViMap,
+        const double &normalizationFactor_translationalDisplacement,
+        const double &normalizationFactor_rotationalDisplacement,
+        const double &scenarioWeight,
+        const ReducedModelOptimization::Settings::ObjectiveWeights &objectiveWeights);
+    static void constructAxialSimulationScenario(
+        const double &forceMagnitude,
+        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            &oppositeInterfaceViPairs,
+        SimulationJob &job);
+
+    static void constructShearSimulationScenario(
+        const double &forceMagnitude,
+        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            &oppositeInterfaceViPairs,
+        SimulationJob &job);
+
+    static void constructBendingSimulationScenario(
+        const double &forceMagnitude,
+        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            &oppositeInterfaceViPairs,
+        SimulationJob &job);
+
+    static void constructDomeSimulationScenario(
+        const double &forceMagnitude,
+        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            &oppositeInterfaceViPairs,
+        SimulationJob &job);
+
+    static void constructSaddleSimulationScenario(
+        const double &forceMagnitude,
+        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            &oppositeInterfaceViPairs,
+        SimulationJob &job);
+    static void constructSSimulationScenario(
+        const double &forceMagnitude,
+        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            &oppositeInterfaceViPairs,
+        SimulationJob &job);
+    static std::function<void(const std::vector<double> &x,
+                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+        function_updateReducedPattern;
+    static std::function<void(const double &newE,
+                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+        function_updateReducedPattern_material_E;
+    static std::function<void(const double &newA,
+                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+        function_updateReducedPattern_material_A;
+    static std::function<void(const double &newI,
+                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+        function_updateReducedPattern_material_I;
+    static std::function<void(const double &newI2,
+                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+        function_updateReducedPattern_material_I2;
+    static std::function<void(const double &newI3,
+                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+        function_updateReducedPattern_material_I3;
+    static std::function<void(const double &newJ,
+                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+        function_updateReducedPattern_material_J;
+    static double objective(const std::vector<double> &x);
+    void initializeUpdateReducedPatternFunctions();
+    //    static double objective(const double &xValue);
+
+    ReducedModelOptimization::Results optimize(
+        ConstPatternGeometry &fullPattern,
+        const ReducedModelOptimization::Settings &optimizationSettings);
+
+    void optimize(ConstPatternGeometry &fullPattern,
+                  ReducedModelOptimization::Settings &optimizationSettings,
+                  ReducedModelOptimization::Results &optimizationResults);
+    static double objective(const std::vector<double> &x,
+                            ReducedModelOptimizer::OptimizationState &optimizationState);
+
+private:
+    void optimize(
+        ReducedModelOptimization::Settings &optimizationSettings,
+        ReducedModelOptimization::Results &results,
+        const std::vector<ReducedModelOptimization::BaseSimulationScenario> &simulationScenarios
+        = std::vector<ReducedModelOptimization::BaseSimulationScenario>(
+            {ReducedModelOptimization::Axial,
+             ReducedModelOptimization::Shear,
+             ReducedModelOptimization::Bending,
+             ReducedModelOptimization::Dome,
+             ReducedModelOptimization::Saddle,
+             ReducedModelOptimization::S}));
+
+    void initializePatterns(PatternGeometry &fullPattern,
+                            ReducedModel &reducedPattern,
+                            const std::array<ReducedModelOptimization::xRange,
+                                             ReducedModelOptimization::NumberOfOptimizationVariables>
+                                &optimizationParameters);
+    static void computeDesiredReducedModelDisplacements(
+        const SimulationResults &fullModelResults,
+        const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
+        Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
+    void runOptimization(const ReducedModelOptimization::Settings &settings,
+                         ReducedModelOptimization::Results &results);
+    void computeMaps(PatternGeometry &fullModel, ReducedModel &reducedModel);
+    void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
+    void initializeOptimizationParameters(
+        const std::shared_ptr<SimulationMesh> &mesh,
+        const std::array<ReducedModelOptimization::xRange,
+                         ReducedModelOptimization::NumberOfOptimizationVariables>
+            &optimizationParamters);
+
+    DRMSimulationModel simulator;
+    void computeObjectiveValueNormalizationFactors(
+        const ReducedModelOptimization::Settings &optimizationSettings);
+
+    void getResults(const FunctionEvaluation &optimalObjective,
+                    const ReducedModelOptimization::Settings &settings,
+                    ReducedModelOptimization::Results &results);
+    double computeFullPatternMaxSimulationForce(
+        const ReducedModelOptimization::BaseSimulationScenario &scenario) const;
+
+#ifdef DLIB_DEFINED
+    static double objective(const dlib::matrix<double, 0, 1> &x);
+#endif
+    std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
+    computeFullPatternMaxSimulationForces(
+        const std::vector<ReducedModelOptimization::BaseSimulationScenario>
+            &desiredBaseSimulationScenario) const;
+    std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationJobs(
+        const std::shared_ptr<SimulationMesh> &pMesh,
+        const std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
+            &baseScenarioMaxForceMagnitudes) const;
+    std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
+    getFullPatternMaxSimulationForces(
+        const std::vector<ReducedModelOptimization::BaseSimulationScenario>
+            &desiredBaseSimulationScenarioIndices,
+        const std::filesystem::path &intermediateResultsDirectoryPath,
+        const bool &recomputeForceMagnitudes);
+    std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
+    getFullPatternMaxSimulationForces();
+    void computeScenarioWeights(
+        const std::vector<ReducedModelOptimization::BaseSimulationScenario> &baseSimulationScenarios,
+        const ReducedModelOptimization::Settings &optimizationSettings);
+};
+inline std::function<void(const double &newE,
+                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_E;
+inline std::function<void(const double &newA,
+                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_A;
+inline std::function<void(const double &newI,
+                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_I;
+inline std::function<void(const double &newI2,
+                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_I2;
+inline std::function<void(const double &newI3,
+                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_I3;
+inline std::function<void(const double &newJ,
+                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_J;
+inline std::function<void(const std::vector<double> &x, std::shared_ptr<SimulationMesh> &m)>
+    ReducedModelOptimizer::function_updateReducedPattern;
+
+extern ReducedModelOptimizer::OptimizationState global;
+
+#endif // REDUCEDMODELOPTIMIZER_HPP
--- a/reducedmodeloptimizer_structs.hpp
+++ b/reducedmodeloptimizer_structs.hpp
@ -15,13 +15,15 @@ enum BaseSimulationScenario {
    Bending,
    Dome,
    Saddle,
+    S,
    NumberOfBaseSimulationScenarios
 };
 inline static std::vector<std::string> baseSimulationScenarioNames{"Axial",
                                                                   "Shear",
                                                                   "Bending",
                                                                   "Dome",
-                                                                   "Saddle"};
+                                                                   "Saddle",
+                                                                   "S"};

 struct Colors
 {
@ -107,36 +109,136 @@ inline int getParameterIndex(const std::string &s)
 struct Settings
 {
    inline static std::string defaultFilename{"OptimizationSettings.json"};
+    //    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,
+    //                                  1.8073431893126763,
+    //                                  0.2864731720436702,
+    //                                  0.14982243562639147,
+    //                                  0.18514829631059054};//median
+    //    std::array<double, NumberOfBaseSimulationScenarios>
+    //        baseScenarioMaxMagnitudes{1.1725844893199244,
+    //                                  0.3464275389927846,
+    //                                  0.09527915004635197,
+    //                                  0.06100757786262501,
+    //                                  0.10631914784812076}; //5_1565 0.03axial
+    //    std::array<double, NumberOfBaseSimulationScenarios>
+    //        baseScenarioMaxMagnitudes{/*15*/ 0 /*1.711973658196369*/,
+    //                                  1.878077115238504,
+    //                                  0.8,
+    //                                  0.15851675178327318,
+    //                                  0.8,
+    //                                  /*1.711973658196369*/ 0}; //custom
+    //    std::array<double, NumberOfBaseSimulationScenarios>
+    //        baseScenarioMaxMagnitudes{0,
+    //                                  14.531854387244818,
+    //                                  0.38321932238436796,
+    //                                  0.21381267870193282,
+    //                                  0.28901381608791094,
+    //                                  1.711973658196369}; //9_14423
+    //    std::array<double, NumberOfBaseSimulationScenarios>
+    //        baseScenarioMaxMagnitudes{1.1725844893199244,
+    //                                  0.3464275389927846,
+    //                                  0.09527915004635197,
+    //                                  0.06100757786262501,
+    //                                  0.10631914784812076}; //5_1565 0.03axial
+    //    std::array<double, NumberOfBaseSimulationScenarios>
+    //        baseScenarioMagnitudes{9.82679, 0.138652, 0.247242, 0.739443, 0.00675865}; //Hyperparam opt
+    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{
+    //        30, 8, 0.4421382884449713, 0.22758433903942452, 0.3247935583883217};
+
+    //    std::array<double, NumberOfBaseSimulationScenarios>
+    //        baseScenarioMagnitudes{10 * 6.310485381644259,
+    //                               10 * 1.7100142258819078,
+    //                               10 * 0.18857048204421728,
+    //                               10 * 0.10813697502645818,
+    //                               10 * 0.11982893539207524}; //6_338
+    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMagnitudes{7.72224,
+    //                                                                               7.72224,
+    //                                                                               0.89468,
+    //                                                                               0.445912,
+    //                                                                               0.625905};
+    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMagnitudes{0.407714,
+    //                                                                               22.3524,
+    //                                                                               0.703164,
+    //                                                                               0.0226138,
+    //                                                                               0.161316};
+    //    std::array<double, NumberOfBaseSimulationScenarios>
+    //        baseScenarioMagnitudes{2, 1, 0.4, 0.2, 0.2}; //8_15444 magnitudes from randomBending0
+
+    //        std::array<double, NumberOfBaseSimulationScenarios>
+    //            baseScenarioMagnitudes{1.0600375226325425,
+    //                                   0.6381040280710403,
+    //                                   0.17201755995098306,
+    //                                   0.0706601822149856,
+    //                                   0.13578373479448072}; //8_15444 magnitudes from displacements
+    //    std::array<double, NumberOfBaseSimulationScenarios>
+    //        baseScenarioMagnitudes{10 * 1.0600375226325425,
+    //                               10 * 0.6381040280710403,
+    //                               10 * 0.17201755995098306,
+    //                               10 * 0.0706601822149856,
+    //                               10 * 0.13578373479448072}; //8_15444 magnitudes from displacements
+
+    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes;
    std::vector<std::vector<OptimizationParameterIndex>> optimizationStrategy = {
-        {E, A, I2, I3, J, R, Theta}};
-    //        {A, I2, I3, J, R, Theta}};
+        //        {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{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},
+    std::array<xRange, NumberOfOptimizationVariables> variablesRanges{xRange{"E", 1e-3, 1e3},
+                                                                      xRange{"A", 1e-3, 1e3},
+                                                                      xRange{"I2", 1e-3, 1e3},
+                                                                      xRange{"I3", 1e-3, 1e3},
+                                                                      xRange{"J", 1e-3, 1e3},
                                                                      xRange{"R", 0.05, 0.95},
                                                                      xRange{"Theta", -30, 30}};
-    int numberOfFunctionCalls{100000};
-    double solverAccuracy{1e-3};
-    double translationEpsilon{3e-3};
+    struct SettingsPSO
+    {
+        int numberOfParticles{200};
+#ifdef USE_PSO
+        inline static std::string optimizerName{"pso"};
+#else
+        inline static std::string optimizerName{"sa"};
+#endif
+    } pso;
+
+    struct SettingsDlibGlobal
+    {
+        int numberOfFunctionCalls{100000};
+    } dlib;
+
+    double solverAccuracy{1e-2};
+    double translationEpsilon{4e-3};
+    //    double translationEpsilon{0};
+    //    double angularDistanceEpsilon{vcg::math::ToRad(2.0)};
    double angularDistanceEpsilon{vcg::math::ToRad(0.0)};
+    double targetBaseTriangleSize{0.03};
+    std::filesystem::path intermediateResultsDirectoryPath;
    struct ObjectiveWeights
    {
        double translational{1.2};
        double rotational{0.8};
-        bool operator==(const ObjectiveWeights &other) const
-        {
-            return this->translational == other.translational
-                   && this->rotational == other.rotational;
-        }
+        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
    {
@ -161,10 +263,8 @@ struct Settings
        inline static std::string ObjectiveWeights{"ObjectiveWeight"};
    };

-    void save(const std::filesystem::path &saveToPath)
+    nlohmann::json toJson() const
    {
-        assert(std::filesystem::is_directory(saveToPath));
-
        nlohmann::json json;
        json[GET_VARIABLE_NAME(optimizationStrategy)] = optimizationStrategy;
        if (optimizationVariablesGroupsWeights.has_value()) {
@ -191,7 +291,6 @@ struct Settings
        //                = xRange.toString();
        //        }

-        json[GET_VARIABLE_NAME(numberOfFunctionCalls)] = numberOfFunctionCalls;
        json[GET_VARIABLE_NAME(solverAccuracy)] = solverAccuracy;
        //Objective weights
        std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios> objectiveWeightsPairs;
@ -204,9 +303,24 @@ struct Settings
                       });
        json[JsonKeys::ObjectiveWeights] = objectiveWeightsPairs;
        json[GET_VARIABLE_NAME(translationEpsilon)] = translationEpsilon;
-        json[GET_VARIABLE_NAME(angularDistanceEpsilon)] = vcg::math::ToDeg(
-            angularDistanceEpsilon);
+        json[GET_VARIABLE_NAME(angularDistanceEpsilon)] = vcg::math::ToDeg(angularDistanceEpsilon);
+        json[GET_VARIABLE_NAME(targetBaseTriangleSize)] = targetBaseTriangleSize;
+        json[GET_VARIABLE_NAME(baseScenarioMaxMagnitudes)] = baseScenarioMaxMagnitudes;

+#ifdef USE_ENSMALLEN
+#ifdef USE_PSO
+        json[GET_VARIABLE_NAME(pso.numberOfParticles)] = pso.numberOfParticles;
+#endif
+        json[GET_VARIABLE_NAME(pso.optimizerName)] = pso.optimizerName;
+#else
+        json[GET_VARIABLE_NAME(dlib.numberOfFunctionCalls)] = dlib.numberOfFunctionCalls;
+#endif
+        return json;
+    }
+    void save(const std::filesystem::path &saveToPath)
+    {
+        assert(std::filesystem::is_directory(saveToPath));
+        nlohmann::json json = toJson();
        std::filesystem::path jsonFilePath(
            std::filesystem::path(saveToPath).append(defaultFilename));
        std::ofstream jsonFile(jsonFilePath.string());
@ -214,20 +328,17 @@ struct Settings
        jsonFile.close();
    }

-    bool load(const std::filesystem::path &jsonFilepath)
+    bool load(const std::filesystem::path &jsonFilePath)
    {
-        assert(std::filesystem::is_directory(loadFromPath));
-        //Load optimal X
-        nlohmann::json json;
-        std::filesystem::path jsonFilepath(
-            std::filesystem::path(loadFromPath).append(JsonKeys::filename));
-        if (!std::filesystem::exists(jsonFilepath)) {
+        if (!std::filesystem::exists(jsonFilePath)) {
            std::cerr << "Optimization settings could not be loaded because input filepath does "
                         "not exist:"
-                      << jsonFilepath << std::endl;
+                      << jsonFilePath << std::endl;
+            assert(false);
            return false;
        }
-        std::ifstream ifs(jsonFilepath);
+        std::ifstream ifs(jsonFilePath);
+        nlohmann::json json;
        ifs >> json;

        if (json.contains(GET_VARIABLE_NAME(optimizationStrategy))) {
@ -260,8 +371,8 @@ struct Settings
            }
        }

-        if (json.contains(GET_VARIABLE_NAME(numberOfFunctionCalls))) {
-            numberOfFunctionCalls = json.at(GET_VARIABLE_NAME(numberOfFunctionCalls));
+        if (json.contains(GET_VARIABLE_NAME(dlib.numberOfFunctionCalls))) {
+            dlib.numberOfFunctionCalls = json.at(GET_VARIABLE_NAME(dlib.numberOfFunctionCalls));
        }
        if (json.contains(GET_VARIABLE_NAME(solverAccuracy))) {
            solverAccuracy = json.at(GET_VARIABLE_NAME(solverAccuracy));
@ -288,33 +399,22 @@ struct Settings
                static_cast<double>(json[GET_VARIABLE_NAME(angularDistanceEpsilon)]));
        }

+        if (json.contains(GET_VARIABLE_NAME(targetBaseTriangleSize))) {
+            targetBaseTriangleSize = static_cast<double>(
+                json[GET_VARIABLE_NAME(targetBaseTriangleSize)]);
+        }
+
+        if (json.contains(GET_VARIABLE_NAME(pso.numberOfParticles))) {
+            pso.numberOfParticles = static_cast<int>(json[GET_VARIABLE_NAME(pso.numberOfParticles)]);
+        }
+
        //        perBaseScenarioObjectiveWeights = json.at(JsonKeys::ObjectiveWeights);
        //        objectiveWeights.translational = json.at(JsonKeys::ObjectiveWeights);
        //        objectiveWeights.rotational = 2 - objectiveWeights.translational;
        return true;
    }

-    std::string toString() const
-    {
-        std::string settingsString;
-        if (!variablesRanges.empty()) {
-            std::string xRangesString;
-            for (const xRange &range : variablesRanges) {
-                xRangesString += range.toString() + " ";
-            }
-            settingsString += xRangesString;
-        }
-        settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls)
-                          + " Accuracy=" + std::to_string(solverAccuracy);
-        for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
-             baseScenario++) {
-            +" W_t=" + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].translational)
-                + " W_r="
-                + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].rotational);
-        }
-
-        return settingsString;
-    }
+    std::string toString() const { return toJson().dump(); }

    void writeHeaderTo(csvFile &os) const
    {
@ -341,7 +441,7 @@ struct Settings
                os << range.min;
            }
        }
-        os << numberOfFunctionCalls;
+        os << dlib.numberOfFunctionCalls;
        os << solverAccuracy;
        os << std::to_string(translationEpsilon);
        os << std::to_string(vcg::math::ToDeg(angularDistanceEpsilon));
@ -378,12 +478,11 @@ struct Settings
                          settings1.perBaseScenarioObjectiveWeights.end(),
                          settings2.perBaseScenarioObjectiveWeights.begin())
                .first
-            != settings1.perBaseScenarioObjectiveWeights.end();
-      return settings1.numberOfFunctionCalls == settings2.numberOfFunctionCalls
+            == settings1.perBaseScenarioObjectiveWeights.end();
+      return settings1.dlib.numberOfFunctionCalls == settings2.dlib.numberOfFunctionCalls
             && settings1.variablesRanges == settings2.variablesRanges
             && settings1.solverAccuracy == settings2.solverAccuracy && haveTheSameObjectiveWeights
-             && settings1.translationEpsilon
-                    == settings2.translationEpsilon;
+             && settings1.translationEpsilon == settings2.translationEpsilon;
  }

  inline void updateMeshWithOptimalVariables(const std::vector<double> &x, SimulationMesh &m)
@ -454,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
@ -467,23 +566,26 @@ 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;

+      inline static std::string DefaultFileName{"OptimizationResults.json"};
+
      struct CSVExportingSettings
      {
          bool exportPE{false};
          bool exportIterationOfMinima{false};
-          bool exportRawObjectiveValue{true};
+          bool exportRawObjectiveValue{false};
          CSVExportingSettings() {}
      };

      const CSVExportingSettings exportSettings;
      struct JsonKeys
      {
-          inline static std::string filename{"OptimizationResults.json"};
          inline static std::string optimizationVariables{"OptimizationVariables"};
          inline static std::string baseTriangle{"BaseTriangle"};
          inline static std::string Label{"Label"};
@ -493,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
@ -536,18 +685,18 @@ 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(
-              std::filesystem::path(saveToPath).append(JsonKeys::filename));
+              std::filesystem::path(saveToPath).append(DefaultFileName));
          std::ofstream jsonFile_optimizationResults(jsonFilePath.string());
          jsonFile_optimizationResults << json_optimizationResults;

@ -582,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);

@ -597,19 +761,22 @@ struct Settings
 #endif
      }

-      bool load(const std::string &loadFromPath, const bool &shouldLoadDebugFiles = false)
+      bool load(const std::filesystem::path &loadFromPath, const bool &shouldLoadDebugFiles = false)
      {
          assert(std::filesystem::is_directory(loadFromPath));
          std::filesystem::path jsonFilepath(
-              std::filesystem::path(loadFromPath).append(JsonKeys::filename));
+              std::filesystem::path(loadFromPath).append(DefaultFileName));
          if (!std::filesystem::exists(jsonFilepath)) {
+              std::cerr << "Input path does not exist:" << loadFromPath << std::endl;
              return false;
          }
          //Load optimal X
          nlohmann::json json_optimizationResults;
-          std::ifstream ifs(std::filesystem::path(loadFromPath).append(JsonKeys::filename));
+          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);
@ -635,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);
@ -659,40 +830,30 @@ struct Settings
              const std::filesystem::path simulationJobsFolderPath(
                  std::filesystem::path(loadFromPath).append("SimulationJobs"));

-              std::vector<std::filesystem::path> sortedByName;
-              for (auto &entry : std::filesystem::directory_iterator(simulationJobsFolderPath))
-                  sortedByName.push_back(entry.path());
+              const std::vector<std::filesystem::path> scenariosSortedByName = [&]() {
+                  std::vector<std::filesystem::path> sortedByName;
+                  for (auto &entry : std::filesystem::directory_iterator(simulationJobsFolderPath))
+                      sortedByName.push_back(entry.path());

-              std::sort(sortedByName.begin(), sortedByName.end(), &Utilities::compareNat);
-
-              for (const auto &simulationScenarioPath : sortedByName) {
+                  std::sort(sortedByName.begin(), sortedByName.end(), &Utilities::compareNat);
+                  return sortedByName;
+              }();
+              for (const auto &simulationScenarioPath : scenariosSortedByName) {
                  if (!std::filesystem::is_directory(simulationScenarioPath)) {
                      continue;
                  }
-                  // Load full pattern files
-                  for (const auto &fileEntry : filesystem::directory_iterator(
-                           std::filesystem::path(simulationScenarioPath).append("Full"))) {
-                      const auto filepath = fileEntry.path();
-                      if (filepath.extension() == ".json") {
-                          SimulationJob job;
-                          job.load(filepath.string());
-                          job.pMesh->setBeamMaterial(0.3, fullPatternYoungsModulus);
-                          fullPatternSimulationJobs.push_back(std::make_shared<SimulationJob>(job));
-                      }
-                  }
-
+                  // Load full job
                  const auto fullJobFilepath = Utilities::getFilepathWithExtension(
                      std::filesystem::path(simulationScenarioPath).append("Full"), ".json");
                  SimulationJob fullJob;
                  fullJob.load(fullJobFilepath.string());
                  fullJob.pMesh->setBeamMaterial(0.3, fullPatternYoungsModulus);
                  fullPatternSimulationJobs.push_back(std::make_shared<SimulationJob>(fullJob));
-
+                  // Load reduced job
                  const auto reducedJobFilepath = Utilities::getFilepathWithExtension(
                      std::filesystem::path(simulationScenarioPath).append("Reduced"), ".json");
                  SimulationJob reducedJob;
                  reducedJob.load(reducedJobFilepath.string());
-                  //                          applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh);
                  applyOptimizationResults_elements(*this, reducedJob.pMesh);
                  reducedPatternSimulationJobs.push_back(
                      std::make_shared<SimulationJob>(reducedJob));
@ -704,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)
@ -716,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,
@ -781,8 +942,8 @@ struct Settings
              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];
+              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+                  Element &e = pReducedPattern_simulationMesh->elements[ei];
                  e.setDimensions(elementDimensions);
              }
          }
@ -792,8 +953,8 @@ struct Settings
          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];
+              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+                  Element &e = pReducedPattern_simulationMesh->elements[ei];
                  e.setMaterial(elementMaterial);
              }
          }
@ -801,8 +962,8 @@ struct Settings
          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];
+              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+                  Element &e = pReducedPattern_simulationMesh->elements[ei];
                  e.dimensions.inertia.J = J;
              }
          }
@ -810,8 +971,8 @@ struct Settings
          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];
+              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+                  Element &e = pReducedPattern_simulationMesh->elements[ei];
                  e.dimensions.inertia.I2 = I2;
              }
          }
@ -819,76 +980,86 @@ struct Settings
          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];
+              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+                  Element &e = pReducedPattern_simulationMesh->elements[ei];
                  e.dimensions.inertia.I3 = I3;
              }
          }
-          pTiledReducedPattern_simulationMesh->reset();
+          pReducedPattern_simulationMesh->reset();
      }

 #if POLYSCOPE_DEFINED
-    void draw() const {
-        PolyscopeInterface::init();
-        DRMSimulationModel drmSimulator;
-        LinearSimulationModel linearSimulator;
-        assert(fullPatternSimulationJobs.size() == reducedPatternSimulationJobs.size());
-        fullPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::fullInitial);
-        reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::reducedInitial, false);
+      void draw(const std::vector<int> &desiredSimulationScenariosIndices = std::vector<int>()) const
+      {
+          PolyscopeInterface::init();
+          DRMSimulationModel drmSimulator;
+          LinearSimulationModel linearSimulator;
+          assert(fullPatternSimulationJobs.size() == reducedPatternSimulationJobs.size());
+          fullPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::fullInitial);
+          reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::reducedInitial, false);

-        const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
-        for (int simulationJobIndex = 0; simulationJobIndex < numberOfSimulationJobs;
-             simulationJobIndex++) {
-            // Drawing of full pattern results
-            const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
-                = fullPatternSimulationJobs[simulationJobIndex];
-            pFullPatternSimulationJob->registerForDrawing(
-                fullPatternSimulationJobs[0]->pMesh->getLabel());
-            DRMSimulationModel::Settings drmSettings;
-            SimulationResults fullModelResults
-                = drmSimulator.executeSimulation(pFullPatternSimulationJob, drmSettings);
-            fullModelResults.registerForDrawing(Colors::fullDeformed, true);
-            //        SimulationResults fullModelLinearResults =
-            //            linearSimulator.executeSimulation(pFullPatternSimulationJob);
-            //        fullModelLinearResults.setLabelPrefix("linear");
-            //        fullModelLinearResults.registerForDrawing(Colors::fullDeformed,false);
-            // Drawing of reduced pattern results
-            const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob
-                = reducedPatternSimulationJobs[simulationJobIndex];
-            //            pReducedPatternSimulationJob->pMesh->registerForDrawing();
-            //            polyscope::show();
-            //                    SimulationResults reducedModelResults = drmSimulator.executeSimulation(
-            //                        pReducedPatternSimulationJob);
-            //        reducedModelResults.registerForDrawing(Colors::reducedDeformed, false);
-            //            SimulationResults reducedModelResults
-            //                = drmSimulator.executeSimulation(pReducedPatternSimulationJob,
-            //                                                 DRMSimulationModel::Settings());
-            SimulationResults reducedModelResults = linearSimulator.executeSimulation(
-                pReducedPatternSimulationJob);
-            reducedModelResults.setLabelPrefix("linear");
-            reducedModelResults.registerForDrawing(Colors::reducedDeformed, true);
-            polyscope::options::programName = fullPatternSimulationJobs[0]->pMesh->getLabel();
-            //            polyscope::view::resetCameraToHomeView();
-            polyscope::show();
-            // Save a screensh
-            const std::string screenshotFilename
-                = "/home/iason/Coding/Projects/Approximating shapes with flat "
-                  "patterns/RodModelOptimizationForPatterns/Results/Images/"
-                  + fullPatternSimulationJobs[0]->pMesh->getLabel() + "_"
-                  + pFullPatternSimulationJob->getLabel();
-            polyscope::screenshot(screenshotFilename, false);
-            pFullPatternSimulationJob->unregister(fullPatternSimulationJobs[0]->pMesh->getLabel());
-            fullModelResults.unregister();
-            //        reducedModelResults.unregister();
-            reducedModelResults.unregister();
-            //        fullModelLinearResults.unregister();
-            //            double error = computeError(reducedModelLinearResults.displacements,
-            //                                        fullModelResults.displacements);
-            //            std::cout << "Error of simulation scenario "
-            //                      << simulationScenarioStrings[simulationScenarioIndex] << " is " << error
-            //                      << std::endl;
+          const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
+          const std::vector<int> scenariosToDraw = [&]() {
+              if (desiredSimulationScenariosIndices.empty()) {
+                  std::vector<int> v(numberOfSimulationJobs);
+                  std::iota(v.begin(), v.end(), 0); //draw all
+                  return v;
+              } else {
+                  return desiredSimulationScenariosIndices;
+              }
+          }();
+
+          for (const int &simulationJobIndex : scenariosToDraw) {
+              // Drawing of full pattern results
+              const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
+                  = fullPatternSimulationJobs[simulationJobIndex];
+              pFullPatternSimulationJob->registerForDrawing(
+                  fullPatternSimulationJobs[0]->pMesh->getLabel());
+              DRMSimulationModel::Settings drmSettings;
+              SimulationResults fullModelResults
+                  = drmSimulator.executeSimulation(pFullPatternSimulationJob, drmSettings);
+              fullModelResults.registerForDrawing(Colors::fullDeformed, true);
+              //        SimulationResults fullModelLinearResults =
+              //            linearSimulator.executeSimulation(pFullPatternSimulationJob);
+              //        fullModelLinearResults.setLabelPrefix("linear");
+              //        fullModelLinearResults.registerForDrawing(Colors::fullDeformed,false);
+              // Drawing of reduced pattern results
+              const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob
+                  = reducedPatternSimulationJobs[simulationJobIndex];
+              //            pReducedPatternSimulationJob->pMesh->registerForDrawing();
+              //            polyscope::show();
+              //                    SimulationResults reducedModelResults = drmSimulator.executeSimulation(
+              //                        pReducedPatternSimulationJob);
+              //        reducedModelResults.registerForDrawing(Colors::reducedDeformed, false);
+              //            SimulationResults reducedModelResults
+              //                = drmSimulator.executeSimulation(pReducedPatternSimulationJob,
+              //                                                 DRMSimulationModel::Settings());
+              SimulationResults reducedModelResults = linearSimulator.executeSimulation(
+                  pReducedPatternSimulationJob);
+              reducedModelResults.setLabelPrefix("linear");
+              reducedModelResults.registerForDrawing(Colors::reducedDeformed, true);
+              polyscope::options::programName = fullPatternSimulationJobs[0]->pMesh->getLabel();
+              //            polyscope::view::resetCameraToHomeView();
+              polyscope::show();
+              // Save a screensh
+              const std::string screenshotFilename
+                  = "/home/iason/Coding/Projects/Approximating shapes with flat "
+                    "patterns/RodModelOptimizationForPatterns/Results/Images/"
+                    + fullPatternSimulationJobs[0]->pMesh->getLabel() + "_"
+                    + pFullPatternSimulationJob->getLabel();
+              polyscope::screenshot(screenshotFilename, false);
+              pFullPatternSimulationJob->unregister(fullPatternSimulationJobs[0]->pMesh->getLabel());
+              fullModelResults.unregister();
+              //        reducedModelResults.unregister();
+              reducedModelResults.unregister();
+              //        fullModelLinearResults.unregister();
+              //            double error = computeError(reducedModelLinearResults.displacements,
+              //                                        fullModelResults.displacements);
+              //            std::cout << "Error of simulation scenario "
+              //                      << simulationScenarioStrings[simulationScenarioIndex] << " is " << error
+              //                      << std::endl;
+          }
      }
-    }
 #endif // POLYSCOPE_DEFINED
    void saveMeshFiles() const {
      const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
@ -1021,5 +1192,9 @@ struct Settings
  };
  enum SimulationModelTag { DRM, Linear };

-  }    // namespace ReducedPatternOptimization
+  inline bool Settings::ObjectiveWeights::operator==(const ObjectiveWeights &other) const
+  {
+      return this->translational == other.translational && this->rotational == other.rotational;
+  }
+  }    // namespace ReducedModelOptimization
 #endif // REDUCEDMODELOPTIMIZER_STRUCTS_HPP
--- a/simulation_structs.hpp
+++ b/simulation_structs.hpp
@ -1,6 +1,12 @@
 #ifndef SIMULATIONSTRUCTS_HPP
 #define SIMULATIONSTRUCTS_HPP
+#include "csvfile.hpp"
+#include "nlohmann/json.hpp"
+#include "simulationmesh.hpp"
+#include "utilities.hpp"
 #include <fstream>
+#include <string>
+#include <vector>

 namespace Eigen {
 template <class Matrix>
@ -33,10 +39,6 @@ void read_binary(const std::string &filename, Matrix &matrix) {
 //}
 } // namespace Eigen

-#include "simulationmesh.hpp"
-#include "nlohmann/json.hpp"
-#include <string>
-#include <vector>

 struct SimulationHistory {
  SimulationHistory() {}
@ -49,8 +51,13 @@ struct SimulationHistory {
  std::vector<size_t> redMarks;
  std::vector<double> greenMarks;
  std::vector<double> residualForcesMovingAverage;
-  std::vector<double> sumOfNormalizedDisplacementNorms;
-  //  std::vector<double> residualForcesMovingAverageDerivativesLog;
+  std::vector<double> perVertexAverageNormalizedDisplacementNorm;
+  std::vector<double> residualForcesMovingAverageDerivativesLog;
+  //internal forces
+  std::vector<double> logOfSumOfAxialForcesNorm;
+  std::vector<double> logOfSumOfTorsionForcesNorm;
+  std::vector<double> logOfSumOfFirstBendingForcesNorm;
+  std::vector<double> logOfSumOfSecondBendingForcesNorm;

  void markRed(const size_t &stepNumber) { redMarks.push_back(stepNumber); }

@ -62,9 +69,47 @@ struct SimulationHistory {
      //    potentialEnergy.push_back(mesh.totalPotentialEnergykN);
      logResidualForces.push_back(std::log10(mesh.totalResidualForcesNorm));
      residualForcesMovingAverage.push_back(std::log10(mesh.residualForcesMovingAverage));
-      sumOfNormalizedDisplacementNorms.push_back(std::log10(mesh.sumOfNormalizedDisplacementNorms));
-      //    residualForcesMovingAverageDerivativesLog.push_back(
-      //        std::log(mesh.residualForcesMovingAverageDerivativeNorm));
+      perVertexAverageNormalizedDisplacementNorm.push_back(
+          mesh.perVertexAverageNormalizedDisplacementNorm);
+      residualForcesMovingAverageDerivativesLog.push_back(
+          std::log(mesh.residualForcesMovingAverageDerivativeNorm));
+
+      //Internal forces
+      const double axialSumOfNorms = std::accumulate(mesh.nodes._handle->data.begin(),
+                                                     mesh.nodes._handle->data.end(),
+                                                     0.0,
+                                                     [](const double &sum, const Node &node) {
+                                                         return sum
+                                                                + node.force.internalAxial.norm();
+                                                     });
+      logOfSumOfAxialForcesNorm.push_back(std::log(axialSumOfNorms));
+
+      const double torsionSumOfNorms
+          = std::accumulate(mesh.nodes._handle->data.begin(),
+                            mesh.nodes._handle->data.end(),
+                            0.0,
+                            [](const double &sum, const Node &node) {
+                                return sum + node.force.internalTorsion.norm();
+                            });
+      logOfSumOfTorsionForcesNorm.push_back(std::log(torsionSumOfNorms));
+
+      const double firstBendingSumOfNorms
+          = std::accumulate(mesh.nodes._handle->data.begin(),
+                            mesh.nodes._handle->data.end(),
+                            0.0,
+                            [](const double &sum, const Node &node) {
+                                return sum + node.force.internalFirstBending.norm();
+                            });
+      logOfSumOfFirstBendingForcesNorm.push_back(std::log(firstBendingSumOfNorms));
+
+      const double secondBendingSumOfNorms
+          = std::accumulate(mesh.nodes._handle->data.begin(),
+                            mesh.nodes._handle->data.end(),
+                            0.0,
+                            [](const double &sum, const Node &node) {
+                                return sum + node.force.internalSecondBending.norm();
+                            });
+      logOfSumOfSecondBendingForcesNorm.push_back(std::log(secondBendingSumOfNorms));
  }

  void clear()
@ -73,7 +118,7 @@ struct SimulationHistory {
      kineticEnergy.clear();
      potentialEnergies.clear();
      residualForcesMovingAverage.clear();
-      sumOfNormalizedDisplacementNorms.clear();
+      perVertexAverageNormalizedDisplacementNorm.clear();
      //      residualForcesMovingAverageDerivativesLog.clear();
  }
 };
@ -190,6 +235,10 @@ public:

    return json.dump();
  }
+  bool operator==(const SimulationJob &otherSimulationJob)
+  {
+      return this->toString() == otherSimulationJob.toString();
+  }

  void clear()
  {
@ -416,11 +465,11 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
          }
      });

-      //      if (!nodeColors.empty()) {
-      //          polyscope::getCurveNetwork(meshLabel)
-      //              ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors)
-      //              ->setEnabled(shouldEnable);
-      //      }
+      if (!nodeColors.empty()) {
+          polyscope::getCurveNetwork(meshLabel)
+              ->addNodeColorQuantity("Boundary conditions", nodeColors)
+              ->setEnabled(shouldEnable);
+      }

      // per node external forces
      std::vector<std::array<double, 3>> externalForces(pMesh->VN());
@ -431,10 +480,11 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
      }

      if (!externalForces.empty()) {
+          const std::string polyscopeLabel_externalForces = "External force";
+          polyscope::getCurveNetwork(meshLabel)->removeQuantity(polyscopeLabel_externalForces);
          polyscope::CurveNetworkNodeVectorQuantity *externalForcesVectors
-              = polyscope::getCurveNetwork(meshLabel)->addNodeVectorQuantity("External force_"
-                                                                                 + label,
-                                                                             externalForces);
+              = polyscope::getCurveNetwork(meshLabel)
+                    ->addNodeVectorQuantity(polyscopeLabel_externalForces, externalForces);

          const std::array<double, 3> color_loads{1.0, 0, 0};
          externalForcesVectors->setVectorColor(
@ -455,7 +505,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m

      if (hasExternalMoments) {
          polyscope::getCurveNetwork(meshLabel)
-              ->addNodeVectorQuantity("External moment_" + label, externalMoments)
+              ->addNodeVectorQuantity("External moment", externalMoments)
              ->setEnabled(shouldEnable);
      }
  }
@ -465,10 +515,10 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
          return;
      }
      if (!nodalExternalForces.empty()) {
-          polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force_" + label);
+          polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force");
      }
      if (!constrainedVertices.empty()) {
-          polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label);
+          polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions");
      }

      // per node external moments
@ -481,7 +531,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
          }
      }
      if (hasExternalMoments) {
-          polyscope::getCurveNetwork(meshLabel)->removeQuantity("External moment_" + label);
+          polyscope::getCurveNetwork(meshLabel)->removeQuantity("External moment");
      }
  }
 #endif // POLYSCOPE_DEFINED
@ -503,6 +553,7 @@ struct SimulationResults
    std::vector<Eigen::Quaternion<double>> rotationalDisplacementQuaternion; //per vertex
    double internalPotentialEnergy{0};
    double executionTime{0};
+    std::vector<std::array<Vector6d, 4>> perVertexInternalForces; //axial,torsion,bending1,bending2
    std::string labelPrefix{"deformed"};
    inline static char deliminator{' '};
    SimulationResults() { pJob = std::make_shared<SimulationJob>(); }
@ -538,6 +589,53 @@ struct SimulationResults
        return m.save(std::filesystem::path(outputFolder).append(getLabel() + ".ply").string());
    }

+    void saveInternalForces(const std::filesystem::path &outputDirPath)
+    {
+        std::cout << "out to:" << outputDirPath << std::endl;
+        const std::filesystem::path internalForcesDirPath = std::filesystem::path(outputDirPath);
+        std::filesystem::create_directories(internalForcesDirPath);
+        csvFile csv_axial6d(std::filesystem::path(internalForcesDirPath)
+                                .append("forces_axial_6d.csv"),
+                            true);
+        csvFile csv_axialMagn(std::filesystem::path(internalForcesDirPath)
+                                  .append("forces_axial_magn.csv"),
+                              true);
+        csvFile csv_torsion6d(std::filesystem::path(internalForcesDirPath)
+                                  .append("forces_torsion_6d.csv"),
+                              true);
+        csvFile csv_torsionMagn(std::filesystem::path(internalForcesDirPath)
+                                    .append("forces_torsion_magn.csv"),
+                                true);
+        csvFile csv_firstBending6d(std::filesystem::path(internalForcesDirPath)
+                                       .append("forces_firstBending_6d.csv"),
+                                   true);
+        csvFile csv_firstBendingMagn(std::filesystem::path(internalForcesDirPath)
+                                         .append("forces_firstBending_magn.csv"),
+                                     true);
+        csvFile csv_secondBending6d(std::filesystem::path(internalForcesDirPath)
+                                        .append("forces_secondBending_6d.csv"),
+                                    true);
+        csvFile csv_secondBendingMagn(std::filesystem::path(internalForcesDirPath)
+                                          .append("forces_secondBending_magn.csv"),
+                                      true);
+        for (const std::array<Vector6d, 4> &internalForce : perVertexInternalForces) {
+            for (int dofi = 0; dofi < 6; dofi++) {
+                csv_axial6d << internalForce[0][dofi];
+                csv_torsion6d << internalForce[1][dofi];
+                csv_firstBending6d << internalForce[2][dofi];
+                csv_secondBending6d << internalForce[3][dofi];
+            }
+            csv_axial6d << endrow;
+            csv_torsion6d << endrow;
+            csv_firstBending6d << endrow;
+            csv_secondBending6d << endrow;
+            csv_axialMagn << internalForce[0].norm() << endrow;
+            csv_torsionMagn << internalForce[1].norm() << endrow;
+            csv_firstBendingMagn << internalForce[2].norm() << endrow;
+            csv_secondBendingMagn << internalForce[3].norm() << endrow;
+        }
+    }
+
    void save(const std::string &outputFolder = std::string())
    {
        const std::filesystem::path outputFolderPath = outputFolder.empty()
@ -559,7 +657,27 @@ struct SimulationResults
        nlohmann::json json;

        json[GET_VARIABLE_NAME(internalPotentialEnergy)] = internalPotentialEnergy;
-
+        //Write internal forces
+        if (!perVertexInternalForces.empty()) {
+            std::vector<Vector6d> internalForces_axial(perVertexInternalForces.size());
+            std::vector<Vector6d> internalForces_torsion(perVertexInternalForces.size());
+            std::vector<Vector6d> internalForces_firstBending(perVertexInternalForces.size());
+            std::vector<Vector6d> internalForces_secondBending(perVertexInternalForces.size());
+            for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
+                internalForces_axial[vi] = perVertexInternalForces[vi][0];
+                internalForces_torsion[vi] = perVertexInternalForces[vi][1];
+                internalForces_firstBending[vi] = perVertexInternalForces[vi][2];
+                internalForces_secondBending[vi] = perVertexInternalForces[vi][3];
+            }
+            json[std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_axial"]
+                = internalForces_axial;
+            json[std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_torsion"]
+                = internalForces_torsion;
+            json[std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_firstBending"]
+                = internalForces_firstBending;
+            json[std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_secondBending"]
+                = internalForces_secondBending;
+        }
        std::filesystem::path jsonFilePath(
            std::filesystem::path(resultsFolderPath).append(defaultJsonFilename));
        std::ofstream jsonFile(jsonFilePath.string());
@ -649,7 +767,7 @@ struct SimulationResults
            const glm::vec3 desiredColor_glm(desiredColor.value()[0],
                                             desiredColor.value()[1],
                                             desiredColor.value()[2]);
-            polyscopeHandle_deformedEdmeMesh->setColor(desiredColor_glm);
+            polyscopeHandle_deformedEdmeMesh->setColor(/*glm::normalize(*/ desiredColor_glm /*)*/);
        }
        Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3);
        Eigen::MatrixX3d framesX(mesh->VN(), 3);
@ -723,10 +841,10 @@ struct SimulationResults
        polyscopeHandle_frameZ->setVectorColor(
            /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));

-        auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel());
-        if (!polyscopeHandle_initialMesh) {
-            polyscopeHandle_initialMesh = mesh->registerForDrawing();
-        }
+        //        if (!polyscope::hasCurveNetwork(mesh->getLabel())) {
+        //            const std::array<double, 3> initialColor({0, 0, 0});
+        //            /*auto polyscopeHandle_initialMesh =*/mesh->registerForDrawing(initialColor);
+        //        }

        //      auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh
        //                                                ->addNodeVectorQuantity("FrameX", framesX_initial);
@ -747,7 +865,7 @@ struct SimulationResults
        //      polyscopeHandle_frameZ_initial->setVectorColor(
        //          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
        //      //}
-        pJob->registerForDrawing(getLabel());
+        pJob->registerForDrawing(getLabel(), false);
        //        static bool wasExecuted =false;
        //        if (!wasExecuted) {
        //            std::function<void()> callback = [&]() {
@ -794,17 +912,52 @@ private:

        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;
-        if (json.contains(GET_VARIABLE_NAME(internalPotentialEnergy))) {
-            internalPotentialEnergy = json.at(GET_VARIABLE_NAME(internalPotentialEnergy));
+        if (std::filesystem::exists(jsonFilepath)) {
+            std::ifstream ifs(jsonFilepath);
+            nlohmann::json json;
+            ifs >> json;
+            //        if (json.contains(GET_VARIABLE_NAME(internalPotentialEnergy))) {
+            //            internalPotentialEnergy = json.at(GET_VARIABLE_NAME(internalPotentialEnergy));
+            //        }
+            if (json.contains(std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_axial")) {
+                perVertexInternalForces.resize(pJob->pMesh->VN());
+                std::vector<Vector6d> perVertexInternalForces_axial
+                    = static_cast<std::vector<Vector6d>>(json.at(
+                        std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_axial"));
+                for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
+                    perVertexInternalForces[vi][0] = perVertexInternalForces_axial[vi];
+                }
+            }
+            if (json.contains(std::string(GET_VARIABLE_NAME(perVertexInternalForces))
+                              + "_torsion")) {
+                perVertexInternalForces.resize(pJob->pMesh->VN());
+                std::vector<Vector6d> perVertexInternalForces_axial
+                    = static_cast<std::vector<Vector6d>>(json.at(
+                        std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_torsion"));
+                for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
+                    perVertexInternalForces[vi][0] = perVertexInternalForces_axial[vi];
+                }
+            }
+            if (json.contains(std::string(GET_VARIABLE_NAME(perVertexInternalForces))
+                              + "_firstBending")) {
+                perVertexInternalForces.resize(pJob->pMesh->VN());
+                std::vector<Vector6d> perVertexInternalForces_axial
+                    = static_cast<std::vector<Vector6d>>(json.at(
+                        std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_firstBending"));
+                for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
+                    perVertexInternalForces[vi][0] = perVertexInternalForces_axial[vi];
+                }
+            }
+            if (json.contains(std::string(GET_VARIABLE_NAME(perVertexInternalForces))
+                              + "_secondBending")) {
+                perVertexInternalForces.resize(pJob->pMesh->VN());
+                std::vector<Vector6d> perVertexInternalForces_axial
+                    = static_cast<std::vector<Vector6d>>(json.at(
+                        std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_secondBending"));
+                for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
+                    perVertexInternalForces[vi][0] = perVertexInternalForces_axial[vi];
+                }
+            }
        }
        return true;
    }
--- 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(),
--- a/simulationmesh.cpp
+++ b/simulationmesh.cpp
@ -64,11 +64,11 @@ SimulationMesh::SimulationMesh(VCGTriMesh &triMesh)
    label = triMesh.getLabel();
    //    eigenEdges = triMesh.getEigenEdges();
    //    if (eigenEdges.rows() == 0) {
-    getEdges(eigenEdges);
+    computeEdges(eigenEdges);
    //    }
    //    eigenVertices = triMesh.getEigenVertices();
    //    if (eigenVertices.rows() == 0) {
-    getVertices(eigenVertices);
+    computeVertices(eigenVertices);
    //    }
    vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);

--- 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(); }
--- a/topologyenumerator.cpp
+++ b/topologyenumerator.cpp
@ -115,7 +115,7 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
                                                         patternGeometryAllEdges,
                                                         allPossibleEdges);
    PatternGeometry patternAllValidEdges;
-    patternAllValidEdges.add(patternGeometryAllEdges.getVertices(), validEdges);
+    patternAllValidEdges.add(patternGeometryAllEdges.computeVertices(), validEdges);
    if (debugIsOn) {
        // Export all valid edges in a ply
        patternAllValidEdges.save(
@ -138,12 +138,12 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
                const size_t ei1 = *setIt;
                vcg::tri::Allocator<PatternGeometry>::AddEdge(
                    intersectingEdgePair,
-                    patternGeometryAllEdges.getVertices()[validEdges[ei0][0]],
-                    patternGeometryAllEdges.getVertices()[validEdges[ei0][1]]);
+                    patternGeometryAllEdges.computeVertices()[validEdges[ei0][0]],
+                    patternGeometryAllEdges.computeVertices()[validEdges[ei0][1]]);
                vcg::tri::Allocator<PatternGeometry>::AddEdge(
                    intersectingEdgePair,
-                    patternGeometryAllEdges.getVertices()[validEdges[ei1][0]],
-                    patternGeometryAllEdges.getVertices()[validEdges[ei1][1]]);
+                    patternGeometryAllEdges.computeVertices()[validEdges[ei1][0]],
+                    patternGeometryAllEdges.computeVertices()[validEdges[ei1][1]]);
                intersectingEdgePair.save(std::filesystem::path(intersectingEdgesPath)
                                              .append(std::to_string(mapIt->first) + "_"
                                                      + std::to_string(*setIt) + ".ply")
@ -192,7 +192,7 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
            computeValidPatterns(numberOfNodesPerSlot,
                                 numberOfEdges,
                                 perEdgeResultPath,
-                                 patternGeometryAllEdges.getVertices(),
+                                 patternGeometryAllEdges.computeVertices(),
                                 intersectingEdges,
                                 validEdges,
                                 interfaceNodes);
@ -213,7 +213,7 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
        computeValidPatterns(numberOfNodesPerSlot,
                             numberOfDesiredEdges,
                             perEdgeResultPath,
-                             patternGeometryAllEdges.getVertices(),
+                             patternGeometryAllEdges.computeVertices(),
                             intersectingEdges,
                             validEdges,
                             interfaceNodes);
@ -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
--- a/trianglepatterngeometry.cpp
+++ b/trianglepatterngeometry.cpp
@ -78,12 +78,13 @@ double PatternGeometry::getTriangleEdgeSize() const { return triangleEdgeSize; }

 PatternGeometry::PatternGeometry() {}

-std::vector<vcg::Point3d> PatternGeometry::getVertices() const {
-  std::vector<VCGEdgeMesh::CoordType> verts(VN());
-  for (size_t vi = 0; vi < VN(); vi++) {
-    verts[vi] = vert[vi].cP();
-  }
-  return verts;
+std::vector<vcg::Point3d> PatternGeometry::computeVertices() const
+{
+    std::vector<VCGEdgeMesh::CoordType> verts(VN());
+    for (size_t vi = 0; vi < VN(); vi++) {
+        verts[vi] = vert[vi].cP();
+    }
+    return verts;
 }

 PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern)
@ -140,13 +141,19 @@ PatternGeometry PatternGeometry::createFan(PatternGeometry &pattern) {
  return fan;
 }

-PatternGeometry::PatternGeometry(PatternGeometry &other) {
-  vcg::tri::Append<PatternGeometry, PatternGeometry>::MeshCopy(*this, other);
-  this->vertices = other.getVertices();
-  baseTriangle = other.getBaseTriangle();
-  baseTriangleHeight = computeBaseTriangleHeight();
-  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
-  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
+void PatternGeometry::updateBaseTriangle()
+{
+    baseTriangle = computeBaseTriangle();
+}
+
+PatternGeometry::PatternGeometry(PatternGeometry &other)
+{
+    vcg::tri::Append<PatternGeometry, PatternGeometry>::MeshCopy(*this, other);
+    this->vertices = other.computeVertices();
+    baseTriangle = other.getBaseTriangle();
+    baseTriangleHeight = computeBaseTriangleHeight();
+    vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
+    vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
 }

 bool PatternGeometry::load(const std::filesystem::__cxx11::path &meshFilePath)
@ -154,6 +161,7 @@ bool PatternGeometry::load(const std::filesystem::__cxx11::path &meshFilePath)
    if (!VCGEdgeMesh::load(meshFilePath)) {
        return false;
    }
+    addNormals();
    vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
    baseTriangleHeight = computeBaseTriangleHeight();
    baseTriangle = computeBaseTriangle();
@ -165,7 +173,7 @@ bool PatternGeometry::load(const std::filesystem::__cxx11::path &meshFilePath)
 void PatternGeometry::add(const std::vector<vcg::Point3d> &vertices) {
  this->vertices = vertices;
  std::for_each(vertices.begin(), vertices.end(), [&](const vcg::Point3d &p) {
-    vcg::tri::Allocator<PatternGeometry>::AddVertex(*this, p);
+      vcg::tri::Allocator<PatternGeometry>::AddVertex(*this, p, DefaultNormal);
  });
  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
@ -419,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;
        }
@ -432,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;
@ -753,27 +776,28 @@ PatternGeometry::getIntersectingEdges(
  return intersectingEdges;
 }

-PatternGeometry::PatternGeometry(const std::string &filename,
-                                 bool addNormalsIfAbsent) {
-  if (!std::filesystem::exists(std::filesystem::path(filename))) {
-    assert(false);
-    std::cerr << "No flat pattern with name " << filename << std::endl;
-    return;
-  }
-  if (!load(filename)) {
-    assert(false);
-    std::cerr << "File could not be loaded  " << filename << std::endl;
-    return;
-  }
-  if (addNormalsIfAbsent) {
-    addNormals();
-  }
+PatternGeometry::PatternGeometry(const std::filesystem::path &patternFilePath,
+                                 bool addNormalsIfAbsent)
+{
+    if (!std::filesystem::exists(std::filesystem::path(patternFilePath))) {
+        assert(false);
+        std::cerr << "No flat pattern with name " << patternFilePath << std::endl;
+        return;
+    }
+    if (!load(patternFilePath)) {
+        assert(false);
+        std::cerr << "File could not be loaded  " << patternFilePath << std::endl;
+        return;
+    }
+    if (addNormalsIfAbsent) {
+        addNormals();
+    }

-  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
-  baseTriangleHeight = computeBaseTriangleHeight();
-  baseTriangle = computeBaseTriangle();
+    vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
+    baseTriangleHeight = computeBaseTriangleHeight();
+    baseTriangle = computeBaseTriangle();

-  updateEigenEdgeAndVertices();
+    updateEigenEdgeAndVertices();
 }

 double PatternGeometry::computeBaseTriangleHeight() const
@ -781,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];
    }
 }

@ -831,6 +859,8 @@ PatternGeometry::PatternGeometry(
  addNormals();
  baseTriangleHeight = computeBaseTriangleHeight();
  baseTriangle = computeBaseTriangle();
+  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
+  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
  updateEigenEdgeAndVertices();
 }

@ -866,8 +896,8 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
    assert(vcg::tri::HasFEAdjacency(tileInto));
    assert(vcg::tri::HasFVAdjacency(tileInto));
    for (const VCGPolyMesh::FaceType &f : tileInto.face) {
-        const int facePatternIndex = perSurfaceFacePatternIndices[tileInto.getIndex(f)];
-        if (facePatternIndex == -1) {
+        const int patternIndex = perSurfaceFacePatternIndices[tileInto.getIndex(f)];
+        if (patternIndex == -1) {
            continue;
        }
        CoordType centerOfFace(0, 0, 0);
@ -883,7 +913,7 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
        for (int &vi : firstInFanConnectToNeighbor_vi) {
            vi += pTiledPattern->VN();
        }
-        ConstPatternGeometry &pattern = patterns[facePatternIndex];
+        ConstPatternGeometry &pattern = patterns[patternIndex];

        for (size_t vi = 0; vi < f.VN(); vi++) {
            auto ep = f.FEp(vi);
@ -905,12 +935,15 @@ 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(
-                    meshTrianglePoints[0] * barycentricCoordinates[facePatternIndex][vi][0]
-                    + meshTrianglePoints[1] * barycentricCoordinates[facePatternIndex][vi][1]
-                    + meshTrianglePoints[2] * barycentricCoordinates[facePatternIndex][vi][2]);
+                    meshTrianglePoints[0] * barycentricCoordinates[patternIndex][vi][0]
+                    + meshTrianglePoints[1] * barycentricCoordinates[patternIndex][vi][1]
+                    + meshTrianglePoints[2] * barycentricCoordinates[patternIndex][vi][2]);
            }

            for (VertexType &v : transformedPattern.vert) {
@ -922,8 +955,12 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
                                                                     transformedPattern,
                                                                     remap);
            for (size_t ei = 0; ei < pattern.EN(); ei++) {
-                perPatternIndexToTiledPatternEdgeIndex[facePatternIndex].push_back(remap.edge[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
@ -937,7 +974,7 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
                            firstInFanConnectToNeighbor_vi[connectToNeighborIndex],
                            pTiledPattern->VN() - pattern.VN()
                                + connectToNeighborsVi[connectToNeighborIndex]);
-                        perPatternIndexToTiledPatternEdgeIndex[facePatternIndex].push_back(
+                        perPatternIndexToTiledPatternEdgeIndex[patternIndex].push_back(
                            pTiledPattern->getIndex(*eIt));
                    }
                }
@ -951,15 +988,11 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
                                + connectToNeighborsVi[connectToNeighborIndex],
                            pTiledPattern->VN() - pattern.VN()
                                + connectToNeighborsVi[connectToNeighborIndex]);
-                        perPatternIndexToTiledPatternEdgeIndex[facePatternIndex].push_back(
+                        perPatternIndexToTiledPatternEdgeIndex[patternIndex].push_back(
                            pTiledPattern->getIndex(*eIt));
                    }
                }
            }
-            //              tiledPattern.updateEigenEdgeAndVertices();
-            //              tiledPattern.registerForDrawing();
-
-            //              polyscope::show();
        }
    }
    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> pu_vertices;
@ -978,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,
@ -986,6 +1020,7 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
                                                        return sum + v.size();
                                                    });

+    const int en = pTiledPattern->EN();
    assert(pTiledPattern->EN() == sumOfEdgeIndices);

    tileIntoEdgesToTiledVi.clear();
@ -1007,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;
 }

@ -1068,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);
@ -1089,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();
    }
--- 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;
@ -40,6 +39,7 @@ private:
                                                           const int &vi);

  public:
+  inline static VectorType DefaultNormal{0.0, 0.0, 1.0};
  PatternGeometry();
  /*The following function should be a copy constructor with
   * a const argument but this is impossible due to the
@ -50,25 +50,26 @@ private:
  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);
-  void add(const std::vector<size_t> &numberOfNodesPerSlot,
-           const std::vector<vcg::Point2i> &edges);
-  static std::vector<vcg::Point3d>
-  constructVertexVector(const std::vector<size_t> &numberOfNodesPerSlot,
-                        const size_t &fanSize, const double &triangleEdgeSize);
+  void add(const std::vector<size_t> &numberOfNodesPerSlot, const std::vector<vcg::Point2i> &edges);
+  static std::vector<vcg::Point3d> constructVertexVector(
+      const std::vector<size_t> &numberOfNodesPerSlot,
+      const size_t &fanSize,
+      const double &triangleEdgeSize);
  bool hasDanglingEdges(const std::vector<size_t> &numberOfNodesPerSlot);
  bool hasValenceGreaterThan(const std::vector<size_t> &numberOfNodesPerSlot,
                             const size_t &valenceThreshold);
-  std::vector<vcg::Point3d> getVertices() const;
+  std::vector<vcg::Point3d> computeVertices() const;
  static PatternGeometry createFan(PatternGeometry &pattern);
  static PatternGeometry createTile(PatternGeometry &pattern);
  double getTriangleEdgeSize() const;
  bool hasUntiledDanglingEdges();
-  std::unordered_map<size_t, std::unordered_set<size_t>>
-  getIntersectingEdges(size_t &numberOfIntersectingEdgePairs) const;
+  std::unordered_map<size_t, std::unordered_set<size_t>> getIntersectingEdges(
+      size_t &numberOfIntersectingEdgePairs) const;

-  static size_t binomialCoefficient(size_t n, size_t m) {
-    assert(n >= m);
-    return tgamma(n + 1) / (tgamma(m + 1) * tgamma(n - m + 1));
+  static size_t binomialCoefficient(size_t n, size_t m)
+  {
+      assert(n >= m);
+      return tgamma(n + 1) / (tgamma(m + 1) * tgamma(n - m + 1));
  }
  bool isFullyConnectedWhenFanned();

@ -83,7 +84,7 @@ private:

  PatternGeometry(const std::vector<size_t> &numberOfNodesPerSlot,
                  const std::vector<vcg::Point2i> &edges);
-  PatternGeometry(const std::string &filename, bool addNormalsIfAbsent = true);
+  PatternGeometry(const std::filesystem::path &patternFilePath, bool addNormalsIfAbsent = true);

  bool createHoneycombAtom();
  void copy(PatternGeometry &copyFrom);
@ -93,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);
--- a/utilities.cpp
+++ b/utilities.cpp
@ -0,0 +1,108 @@
+#include "utilities.hpp"
+#include "matplot/matplot.h"
+
+void Utilities::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);
+    }
+}
+#ifdef POLYSCOPE_DEFINED
+#include "polyscope/curve_network.h"
+#include "polyscope/pick.h"
+#include "polyscope/polyscope.h"
+#include <functional>
+
+void PolyscopeInterface::mainCallback()
+{
+    ImGui::PushItemWidth(100); // Make ui elements 100 pixels wide,
+        // instead of full width. Must have
+        // matching PopItemWidth() below.
+
+    for (std::function<void()> &userCallback : globalPolyscopeData.userCallbacks) {
+        userCallback();
+    }
+
+    ImGui::PopItemWidth();
+}
+
+void PolyscopeInterface::addUserCallback(const std::function<void()> &userCallback)
+{
+    globalPolyscopeData.userCallbacks.push_back(userCallback);
+}
+
+void PolyscopeInterface::deinitPolyscope()
+{
+    if (!polyscope::state::initialized) {
+        return;
+    }
+
+    polyscope::render::engine->shutdownImGui();
+}
+
+void PolyscopeInterface::init()
+{
+    if (polyscope::state::initialized) {
+        return;
+    }
+    polyscope::init();
+    polyscope::options::groundPlaneEnabled = false;
+    polyscope::view::upDir = polyscope::view::UpDir::ZUp;
+
+    polyscope::state::userCallback = &mainCallback;
+    polyscope::options::autocenterStructures = false;
+    polyscope::options::autoscaleStructures = false;
+}
+
+std::pair<PolyscopeInterface::PolyscopeLabel, size_t> PolyscopeInterface::getSelection()
+{
+    std::pair<polyscope::Structure *, size_t> selection = polyscope::pick::getSelection();
+    if (selection.first == nullptr) {
+        return std::make_pair(std::string(), 0);
+    }
+    return std::make_pair(selection.first->name, selection.second);
+}
+
+void PolyscopeInterface::registerWorldAxes()
+{
+    PolyscopeInterface::init();
+
+    Eigen::MatrixX3d axesPositions(4, 3);
+    axesPositions.row(0) = Eigen::Vector3d(0, 0, 0);
+    //    axesPositions.row(1) = Eigen::Vector3d(polyscope::state::lengthScale, 0, 0);
+    //    axesPositions.row(2) = Eigen::Vector3d(0, polyscope::state::lengthScale, 0);
+    //    axesPositions.row(3) = Eigen::Vector3d(0, 0, polyscope::state::lengthScale);
+    axesPositions.row(1) = Eigen::Vector3d(1, 0, 0);
+    axesPositions.row(2) = Eigen::Vector3d(0, 1, 0);
+    axesPositions.row(3) = Eigen::Vector3d(0, 0, 1);
+
+    Eigen::MatrixX2i axesEdges(3, 2);
+    axesEdges.row(0) = Eigen::Vector2i(0, 1);
+    axesEdges.row(1) = Eigen::Vector2i(0, 2);
+    axesEdges.row(2) = Eigen::Vector2i(0, 3);
+    Eigen::MatrixX3d axesColors(3, 3);
+    axesColors.row(0) = Eigen::Vector3d(1, 0, 0);
+    axesColors.row(1) = Eigen::Vector3d(0, 1, 0);
+    axesColors.row(2) = Eigen::Vector3d(0, 0, 1);
+
+    const std::string worldAxesName = "World Axes";
+    polyscope::registerCurveNetwork(worldAxesName, axesPositions, axesEdges);
+    polyscope::getCurveNetwork(worldAxesName)->setRadius(0.0001, false);
+    const std::string worldAxesColorName = worldAxesName + " Color";
+    polyscope::getCurveNetwork(worldAxesName)
+        ->addEdgeColorQuantity(worldAxesColorName, axesColors)
+        ->setEnabled(true);
+}
+
+#endif
--- a/utilities.hpp
+++ b/utilities.hpp
@ -10,6 +10,7 @@
 #include <iterator>
 #include <numeric>
 #include <regex>
+#include <string_view>

 #define GET_VARIABLE_NAME(Variable) (#Variable)

@ -142,6 +143,15 @@ struct Vector6d : public std::array<double, 6> {
 };

 namespace Utilities {
+template<typename T>
+std::string to_string_with_precision(const T a_value, const int n = 2)
+{
+    std::ostringstream out;
+    out.precision(n);
+    out << std::fixed << a_value;
+    return out.str();
+}
+
 inline bool compareNat(const std::string &a, const std::string &b)
 {
    if (a.empty())
@ -314,6 +324,7 @@ inline std::vector<Vector6d> fromEigenMatrix(const Eigen::MatrixXd &m)
 inline std::filesystem::path getFilepathWithExtension(const std::filesystem::path &folderPath,
                                                      const std::string &extension)
 {
+    assert(std::filesystem::exists(folderPath));
    for (const std::filesystem::directory_entry &dirEntry :
         std::filesystem::directory_iterator(folderPath)) {
        if (dirEntry.is_regular_file() && std::filesystem::path(dirEntry).extension() == extension) {
@ -324,13 +335,17 @@ inline std::filesystem::path getFilepathWithExtension(const std::filesystem::pat
    return "";
 }

+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 = {});
+
 } // namespace Utilities

 #ifdef POLYSCOPE_DEFINED
-#include "polyscope/curve_network.h"
-#include "polyscope/pick.h"
-#include "polyscope/polyscope.h"
-#include <functional>

 namespace PolyscopeInterface {
 inline struct GlobalPolyscopeData
@ -338,86 +353,17 @@ inline struct GlobalPolyscopeData
    std::vector<std::function<void()>> userCallbacks;
 } globalPolyscopeData;

-inline void mainCallback()
-{
-    ImGui::PushItemWidth(100); // Make ui elements 100 pixels wide,
-                               // instead of full width. Must have
-                               // matching PopItemWidth() below.
+void mainCallback();

-    for (std::function<void()> &userCallback : globalPolyscopeData.userCallbacks) {
-        userCallback();
-    }
+void addUserCallback(const std::function<void()> &userCallback);

-    ImGui::PopItemWidth();
-}
+void deinitPolyscope();

-inline void addUserCallback(const std::function<void()> &userCallback)
-{
-    globalPolyscopeData.userCallbacks.push_back(userCallback);
-}
-
-inline void deinitPolyscope()
-{
-    if (!polyscope::state::initialized) {
-        return;
-    }
-
-    polyscope::render::engine->shutdownImGui();
-}
-
-inline void init()
-{
-    if (polyscope::state::initialized) {
-        return;
-    }
-    polyscope::init();
-    polyscope::options::groundPlaneEnabled = false;
-    polyscope::view::upDir = polyscope::view::UpDir::ZUp;
-
-    polyscope::state::userCallback = &mainCallback;
-    polyscope::options::autocenterStructures = false;
-    polyscope::options::autoscaleStructures = false;
-}
+void init();
 using PolyscopeLabel = std::string;
-inline std::pair<PolyscopeLabel, size_t> getSelection()
-{
-    std::pair<polyscope::Structure *, size_t> selection = polyscope::pick::getSelection();
-    if (selection.first == nullptr) {
-        return std::make_pair(std::string(), 0);
-    }
-    return std::make_pair(selection.first->name, selection.second);
-}
+std::pair<PolyscopeLabel, size_t> getSelection();

-inline void registerWorldAxes()
-{
-    PolyscopeInterface::init();
-
-    Eigen::MatrixX3d axesPositions(4, 3);
-    axesPositions.row(0) = Eigen::Vector3d(0, 0, 0);
-    //    axesPositions.row(1) = Eigen::Vector3d(polyscope::state::lengthScale, 0, 0);
-    //    axesPositions.row(2) = Eigen::Vector3d(0, polyscope::state::lengthScale, 0);
-    //    axesPositions.row(3) = Eigen::Vector3d(0, 0, polyscope::state::lengthScale);
-    axesPositions.row(1) = Eigen::Vector3d(1, 0, 0);
-    axesPositions.row(2) = Eigen::Vector3d(0, 1, 0);
-    axesPositions.row(3) = Eigen::Vector3d(0, 0, 1);
-
-    Eigen::MatrixX2i axesEdges(3, 2);
-    axesEdges.row(0) = Eigen::Vector2i(0, 1);
-    axesEdges.row(1) = Eigen::Vector2i(0, 2);
-    axesEdges.row(2) = Eigen::Vector2i(0, 3);
-    Eigen::MatrixX3d axesColors(3, 3);
-    axesColors.row(0) = Eigen::Vector3d(1, 0, 0);
-    axesColors.row(1) = Eigen::Vector3d(0, 1, 0);
-    axesColors.row(2) = Eigen::Vector3d(0, 0, 1);
-
-    const std::string worldAxesName = "World Axes";
-    polyscope::registerCurveNetwork(worldAxesName, axesPositions, axesEdges);
-    polyscope::getCurveNetwork(worldAxesName)->setRadius(0.0001, false);
-    const std::string worldAxesColorName = worldAxesName + " Color";
-    polyscope::getCurveNetwork(worldAxesName)
-        ->addEdgeColorQuantity(worldAxesColorName, axesColors)
-        ->setEnabled(true);
-}
+void registerWorldAxes();
 } // namespace PolyscopeInterface

 #endif
@ -440,15 +386,6 @@ template <typename T> std::string toString(const T &v) {
         std::to_string(v[2]) + ")";
 }

-template<typename T>
-std::string to_string_with_precision(const T a_value, const int n = 2)
-{
-    std::ostringstream out;
-    out.precision(n);
-    out << std::fixed << a_value;
-    return out.str();
-}
-
 template<typename T>
 size_t computeHashUnordered(const std::vector<T> &v)
 {
--- a/vcgtrimesh.cpp
+++ b/vcgtrimesh.cpp
@ -6,6 +6,9 @@
 #include <wrap/io_trimesh/export.h>
 #include <wrap/io_trimesh/import.h>
 //#include <wrap/nanoply/include/nanoplyWrapper.hpp>
+#ifdef POLYSCOPE_DEFINED
+#include <polyscope/curve_network.h>
+#endif

 bool VCGTriMesh::load(const std::filesystem::__cxx11::path &meshFilePath) {
    assert(std::filesystem::exists(meshFilePath));
@ -36,6 +39,32 @@ bool VCGTriMesh::load(const std::filesystem::__cxx11::path &meshFilePath) {
    return true;
 }

+bool VCGTriMesh::load(std::istringstream &offInputStream)
+{
+    Clear();
+    //  assert(plyFileHasAllRequiredFields(plyFilename));
+    // Load the ply file
+    int mask = 0;
+    mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
+    mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
+    mask |= vcg::tri::io::Mask::IOM_VERTCOLOR;
+    mask |= vcg::tri::io::Mask::IOM_EDGEINDEX;
+    const bool openingFromStreamErrorCode
+        = vcg::tri::io::ImporterOFF<VCGTriMesh>::OpenStream(*this, offInputStream, mask);
+    if (openingFromStreamErrorCode != 0) {
+        std::cerr << "Error reading from stream:"
+                  << vcg::tri::io::ImporterOFF<VCGTriMesh>::ErrorMsg(openingFromStreamErrorCode)
+                  << std::endl;
+        return false;
+    }
+    vcg::tri::UpdateTopology<VCGTriMesh>::AllocateEdge(*this);
+    vcg::tri::UpdateTopology<VCGTriMesh>::FaceFace(*this);
+    vcg::tri::UpdateTopology<VCGTriMesh>::VertexFace(*this);
+    vcg::tri::UpdateTopology<VCGTriMesh>::VertexEdge(*this);
+    vcg::tri::UpdateNormal<VCGTriMesh>::PerVertexNormalized(*this);
+    return true;
+}
+
 Eigen::MatrixX3d VCGTriMesh::getVertices() const
 {
    //    vcg::tri::Allocator<VCGTriMesh>::CompactVertexVector(m);
--- a/vcgtrimesh.hpp
+++ b/vcgtrimesh.hpp
@ -40,6 +40,7 @@ public:
  VCGTriMesh();
  VCGTriMesh(const std::string &filename);
  bool load(const std::filesystem::path &meshFilePath) override;
+  bool load(std::istringstream &offInputStream);
  Eigen::MatrixX3d getVertices() const;
  Eigen::MatrixX3i getFaces() const;
  bool save(const std::string plyFilename);
Author	SHA1	Message	Date
iasonmanolas	3e5dc29bc1	Removed file	2022-05-06 16:40:31 +03:00
iasonmanolas	b2bd12d8df	Hardcoded the 22 hexagon test-surface	2022-05-06 16:40:05 +03:00
iasonmanolas	90dc2ac081	Refactoring. Replaced PatternGeometry with ReducedModel everywhere	2022-05-06 16:38:50 +03:00
iasonmanolas	7f11ea8a47	refactoring. Added functionality to export internal forces	2022-05-06 16:36:47 +03:00
iasonmanolas	f531b16b19	Added function for loading tri mes using istringstream	2022-05-06 16:28:28 +03:00
iasonmanolas	398df24056	Moved implementations to cpp file	2022-05-06 16:27:37 +03:00
iasonmanolas	0cb175e72e	Added functions for changing the geometry and material props of the reduced model	2022-05-06 16:27:00 +03:00
iasonmanolas	a7fdf431dd	Refactoring	2022-05-06 16:26:05 +03:00
iasonmanolas	e9707e2cfb	Refactoring	2022-05-06 16:16:43 +03:00
iasonmanolas	436ece0d88	Refactoring	2022-05-06 16:12:42 +03:00
iasonmanolas	3d0de46922	Added the S scenario. Refactoring	2022-02-18 17:50:16 +02:00
iasonmanolas	96807c3b85	Moved global variable to struct when ensmallen is used. And left the global variable when dlib is used for optimization. Added the S scenario. Removed static keyword from several functions	2022-02-18 17:49:17 +02:00
iasonmanolas	5c4f8c0bd5	Added Settings struct. Expanded printing of evaluation results.Refactoring	2022-02-18 17:46:28 +02:00
iasonmanolas	45eed0e3da	Moved string with precision function to Utilities namespace	2022-02-18 17:45:14 +02:00
iasonmanolas	bef2ae8860	Renamed getEdges and getVertices to computeEdges and computeVertices	2022-02-18 17:44:33 +02:00
iasonmanolas	ec7160e637	Renamed getEdges and getVertices to computeEdges and computeVertices	2022-02-18 17:43:07 +02:00
iasonmanolas	6c16d8d48e	Renamed getEdges and getVertices to computeEdges and computeVertices	2022-02-18 17:42:47 +02:00
iasonmanolas	5f863af7ce	Refactoring	2022-02-18 17:41:37 +02:00
iasonmanolas	e1f66515d5	Renamed getEdges and getVertices to computeEdges and computeVertices	2022-02-18 17:40:56 +02:00
iasonmanolas	9e121beade	Made parse function of csvFile struct return a vector<vector<string>> instead of being a template function	2022-02-18 17:39:02 +02:00
iasonmanolas	33facf05ab	Fixed bug	2022-02-01 13:08:07 +02:00
iasonmanolas	e54dc0a27c	Changed hardcoded directory of tessellated results	2022-02-01 13:06:53 +02:00
iasonmanolas	7f543ef21a	Added computation of internal forces in the converged state	2022-02-01 13:06:16 +02:00
iasonmanolas	366727ced6	Refactoring	2022-01-28 20:06:56 +02:00
iasonmanolas	068626f299	Refactoring	2022-01-28 20:01:44 +02:00
iasonmanolas	c513abc45b	Added dlib to dependencies since I move the reducedmodeloptimizer to MySources. Refacoting	2022-01-28 19:59:42 +02:00
iasonmanolas	8e33ba35aa	Moved the reduced model optimizer class to MySources	2022-01-27 14:45:54 +02:00
iasonmanolas	38c2535a34	Added reduced model class that handles the reduced model in terms of construction	2022-01-27 14:45:16 +02:00
iasonmanolas	643e8b35be	Added reduced model class that handles the reduced model in terms of construction	2022-01-27 14:44:30 +02:00