Decreased the range of the E,A,I2,I3,J of optimization variables to [0.01,100]

CMakeLists.txt

@@ -1,7 +1,10 @@
 cmake_minimum_required(VERSION 2.8)
 project(MySources)
-set(CMAKE_CXX_STANDARD 20)
-set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+file(GLOB MySourcesFiles ${CMAKE_CURRENT_LIST_DIR}/*.hpp ${CMAKE_CURRENT_LIST_DIR}/*.cpp)
+add_library(${PROJECT_NAME} ${MySourcesFiles} )
+target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20)
+
 #Download external dependencies NOTE: If the user has one of these libs it shouldn't be downloaded again.
 include(${CMAKE_CURRENT_LIST_DIR}/cmake/DownloadProject.cmake)
 if (CMAKE_VERSION VERSION_LESS 3.2)
@@ -9,101 +12,114 @@ if (CMAKE_VERSION VERSION_LESS 3.2)
 else()
     set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1")
 endif()
-
 set(EXTERNAL_DEPS_DIR "/home/iason/Coding/build/external dependencies")
 if(NOT EXISTS ${EXTERNAL_DEPS_DIR})
     set(EXTERNAL_DEPS_DIR ${CMAKE_CURRENT_BINARY_DIR})
 endif()
 ##Create directory for the external libraries
 file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR})
-#message(${POLYSCOPE_ALREADY_COMPILED})
-
 if(${USE_POLYSCOPE})
-download_project(PROJ               POLYSCOPE
-                 GIT_REPOSITORY     https://github.com/nmwsharp/polyscope.git
-                 GIT_TAG	    master
-                 PREFIX             ${EXTERNAL_DEPS_DIR}
-                 ${UPDATE_DISCONNECTED_IF_AVAILABLE}
-)
-if(NOT EXISTS ${POLYSCOPE_BINARY_DIR})
-add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR})
-add_compile_definitions(POLYSCOPE_DEFINED)
+        set(POLYSCOPE_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/polyscope)
+        download_project(PROJ               POLYSCOPE
+                         GIT_REPOSITORY     https://github.com/nmwsharp/polyscope.git
+                         GIT_TAG	    master
+                         PREFIX             ${EXTERNAL_DEPS_DIR}
+                         BINARY_DIR		${POLYSCOPE_BINARY_DIR}
+                         ${UPDATE_DISCONNECTED_IF_AVAILABLE}
+        )
+        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)
 endif()
-endif()
-
-#dlib
-set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin)
-file(MAKE_DIRECTORY ${DLIB_BIN_DIR})
-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_BINARY_DIR})
 
 ##vcglib devel branch
 download_project(PROJ		    vcglib_devel
-                 GIT_REPOSITORY     https://github.com/IasonManolas/vcglib.git
+                 GIT_REPOSITORY     https://gitea-s2i2s.isti.cnr.it/manolas/vcglib.git
                  GIT_TAG	    devel
                  PREFIX             ${EXTERNAL_DEPS_DIR}
                  ${UPDATE_DISCONNECTED_IF_AVAILABLE}
 )
+add_subdirectory(${vcglib_devel_SOURCE_DIR} ${vcglib_devel_BINARY_DIR})
+target_sources(${PROJECT_NAME} PUBLIC ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)
 
 ##matplot++ lib
+set(MATPLOTPLUSPLUS_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/matplot)
 download_project(PROJ               MATPLOTPLUSPLUS
                  GIT_REPOSITORY     https://github.com/alandefreitas/matplotplusplus
                  GIT_TAG	    master
+                 BINARY_DIR		${MATPLOTPLUSPLUS_BINARY_DIR}
                  PREFIX             ${EXTERNAL_DEPS_DIR}
                  ${UPDATE_DISCONNECTED_IF_AVAILABLE}
                  )
 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}
                  ${UPDATE_DISCONNECTED_IF_AVAILABLE}
 )
 add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR})
 
 ##TBB
+set(TBB_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/tbb)
 download_project(PROJ                TBB
                  GIT_REPOSITORY      https://github.com/wjakob/tbb.git
                  GIT_TAG             master
                  PREFIX             ${EXTERNAL_DEPS_DIR}
+                 BINARY_DIR		${TBB_BINARY_DIR}
                  ${UPDATE_DISCONNECTED_IF_AVAILABLE}
 )
 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)
         add_compile_definitions(_HAS_STD_BYTE=0)
 endif(MSVC)
-
 #link_directories(${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
-file(GLOB MySourcesFiles ${CMAKE_CURRENT_LIST_DIR}/*.hpp ${CMAKE_CURRENT_LIST_DIR}/*.cpp)
-add_library(${PROJECT_NAME} ${MySourcesFiles} ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)
-
 
+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()
+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 ${MySourcesFiles}
+        PUBLIC matplot
     )
 
-if(${MYSOURCES_STATIC_LINK})
-        target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread)
-elseif(${USE_POLYSCOPE})
-        target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib glad ThreedBeamFEA tbb pthread polyscope)
-else()
-        target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib glad ThreedBeamFEA tbb pthread)
+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()
-target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)

beam.hpp

@@ -9,15 +9,36 @@ struct RectangularBeamDimensions {
   inline static std::string name{"Rectangular"};
   double b;
   double h;
-  RectangularBeamDimensions(const double &width, const double &height)
-      : b(width), h(height) {
-    assert(width > 0 && height > 0);
+
+  double A{0}; // cross sectional area
+
+  struct MomentsOfInertia
+  {
+      double I2{0}; // second moment of inertia
+      double I3{0}; // third moment of inertia
+      double J{0};  // torsional constant (polar moment of inertia)
+  } inertia;
+
+  RectangularBeamDimensions(const double &width, const double &height) : b(width), h(height)
+  {
+      assert(width > 0 && height > 0);
+      updateProperties();
   }
-  RectangularBeamDimensions() : b(0.002), h(0.002) {}
+  RectangularBeamDimensions() : b(0.002), h(0.002) { updateProperties(); }
   std::string toString() const {
     return std::string("b=") + std::to_string(b) + std::string(" h=") +
            std::to_string(h);
   }
+
+  void updateProperties()
+  {
+      A = b * h;
+      inertia.I2 = b * std::pow(h, 3) / 12;
+      inertia.I3 = h * std::pow(b, 3) / 12;
+      inertia.J = inertia.I2 + inertia.I3;
+  }
+  static void computeMomentsOfInertia(const RectangularBeamDimensions &dimensions,
+                                      MomentsOfInertia &inertia);
 };
 
 struct CylindricalBeamDimensions {
@@ -39,18 +60,21 @@ struct ElementMaterial
 {
     double poissonsRatio;
     double youngsModulus;
+    double G;
     ElementMaterial(const double &poissonsRatio, const double &youngsModulus)
         : poissonsRatio(poissonsRatio), youngsModulus(youngsModulus)
     {
         assert(poissonsRatio <= 0.5 && poissonsRatio >= -1);
+        updateProperties();
     }
-    ElementMaterial() : poissonsRatio(0.3), youngsModulus(200 * 1e9) {}
+    ElementMaterial() : poissonsRatio(0.3), youngsModulus(200 * 1e9) { updateProperties(); }
     std::string toString() const
     {
         return std::string("Material:") + std::string("\nPoisson's ratio=")
                + std::to_string(poissonsRatio) + std::string("\nYoung's Modulus(GPa)=")
                + std::to_string(youngsModulus / 1e9);
     }
+    void updateProperties() { G = youngsModulus / (2 * (1 + poissonsRatio)); }
 };
 
 #endif // BEAM_HPP

csvfile.hpp

@@ -21,23 +21,24 @@ class csvFile {
   const std::string escape_seq_;
   const std::string special_chars_;
 
-public:
-  csvFile(const std::string filename, const bool &overwrite,
+  public:
+  csvFile(const std::filesystem::path &filename,
+          const bool &overwrite,
           const std::string separator = ",")
-      : fs_(), is_first_(true), separator_(separator), escape_seq_("\""),
-        special_chars_("\"") {
-    fs_.exceptions(std::ios::failbit | std::ios::badbit);
-    if (filename.empty()) {
-      fs_.copyfmt(std::cout);
-      fs_.clear(std::cout.rdstate());
-      fs_.basic_ios<char>::rdbuf(std::cout.rdbuf());
-    } else {
-        if (!std::filesystem::exists(std::filesystem::path(filename))) {
-            std::ofstream outfile(filename);
-            outfile.close();
-        }
-        overwrite ? fs_.open(filename, std::ios::trunc) : fs_.open(filename, std::ios::app);
-    }
+      : fs_(), is_first_(true), separator_(separator), escape_seq_("\""), special_chars_("\"")
+  {
+      fs_.exceptions(std::ios::failbit | std::ios::badbit);
+      if (filename.empty()) {
+          fs_.copyfmt(std::cout);
+          fs_.clear(std::cout.rdstate());
+          fs_.basic_ios<char>::rdbuf(std::cout.rdbuf());
+      } else {
+          if (!std::filesystem::exists(filename)) {
+              std::ofstream outfile(filename);
+              outfile.close();
+          }
+          overwrite ? fs_.open(filename, std::ios::trunc) : fs_.open(filename, std::ios::app);
+      }
   }
 
   ~csvFile() {

drmsimulationmodel.cpp

@@ -217,28 +217,16 @@ void DRMSimulationModel::runUnitTests()
 
     //  polyscope::show();
 }
-
-void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings)
+void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob)
 {
-    mSettings = settings;
-    mCurrentSimulationStep = 0;
-    history.clear();
-    history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
-    constrainedVertices.clear();
+    //#ifdef USE_ENSMALLEN
+    //    this->pJob = pJob;
+    //#endif
     pMesh.reset();
     pMesh = std::make_unique<SimulationMesh>(*pJob->pMesh);
     vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
-    plotYValues.clear();
-    plotHandle.reset();
-    checkedForMaximumMoment = false;
-    externalMomentsNorm = 0;
-    Dt = settings.Dtini;
-    numOfDampings = 0;
-    shouldTemporarilyDampForces = false;
-    externalLoadStep = 1;
-    isVertexConstrained.clear();
-    minTotalResidualForcesNorm = std::numeric_limits<double>::max();
 
+    constrainedVertices.clear();
     constrainedVertices = pJob->constrainedVertices;
     if (!pJob->nodalForcedDisplacements.empty()) {
         for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair : pJob->nodalForcedDisplacements) {
@@ -251,6 +239,26 @@ void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const
     for (auto fixedVertex : constrainedVertices) {
         isVertexConstrained[fixedVertex.first] = true;
     }
+}
+
+void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings)
+{
+    mSettings = settings;
+    mCurrentSimulationStep = 0;
+    history.clear();
+    history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel();
+    plotYValues.clear();
+    plotHandle.reset();
+    checkedForMaximumMoment = false;
+    externalMomentsNorm = 0;
+    Dt = settings.Dtini;
+    numOfDampings = 0;
+    shouldTemporarilyDampForces = false;
+    externalLoadStep = 1;
+    isVertexConstrained.clear();
+    minTotalResidualForcesNorm = std::numeric_limits<double>::max();
+
+    reset(pJob);
 
 #ifdef POLYSCOPE_DEFINED
     if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) {
@@ -788,21 +796,21 @@ double DRMSimulationModel::computeTotalInternalPotentialEnergy()
 
         const EdgeIndex ei = pMesh->getIndex(e);
         const double e_k = element.length - element.initialLength;
-        const double axialPE = pow(e_k, 2) * element.material.youngsModulus * element.A
+        const double axialPE = pow(e_k, 2) * element.material.youngsModulus * element.dimensions.A
                                / (2 * element.initialLength);
         const double theta1Diff = theta1_jplus1 - theta1_j;
-        const double torsionalPE = element.G * element.inertia.J * pow(theta1Diff, 2)
-                                   / (2 * element.initialLength);
+        const double torsionalPE = element.material.G * element.dimensions.inertia.J
+                                   * pow(theta1Diff, 2) / (2 * element.initialLength);
         const double firstBendingPE_inBracketsTerm = 4 * pow(theta2_j, 2)
                                                      + 4 * theta2_j * theta2_jplus1
                                                      + 4 * pow(theta2_jplus1, 2);
         const double firstBendingPE = firstBendingPE_inBracketsTerm * element.material.youngsModulus
-                                      * element.inertia.I2 / (2 * element.initialLength);
+                                      * element.dimensions.inertia.I2 / (2 * element.initialLength);
         const double secondBendingPE_inBracketsTerm = 4 * pow(theta3_j, 2)
                                                       + 4 * theta3_j * theta3_jplus1
                                                       + 4 * pow(theta3_jplus1, 2);
         const double secondBendingPE = secondBendingPE_inBracketsTerm
-                                       * element.material.youngsModulus * element.inertia.I3
+                                       * element.material.youngsModulus * element.dimensions.inertia.I3
                                        / (2 * element.initialLength);
 
         totalInternalPotentialEnergy += axialPE + torsionalPE + firstBendingPE + secondBendingPE;
@@ -881,25 +889,26 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
                 const Node &edgeNode = pMesh->nodes[ev];
                 internalForcesContributionFromThisEdge[evi].first = edgeNode.vi;
 
-            const VertexPointer &rev_j = edgeNode.referenceElement->cV(0);
-            const VertexPointer &rev_jplus1 = edgeNode.referenceElement->cV(1);
-            // refElemOtherVertex can be precomputed
-            const VertexPointer &refElemOtherVertex = rev_j == &ev ? rev_jplus1 : rev_j;
-            const Node &refElemOtherVertexNode = pMesh->nodes[refElemOtherVertex];
-            if (edgeNode.referenceElement != &e) {
-                internalForcesContributionFromThisEdge[evi + 2].first = refElemOtherVertexNode.vi;
-            }
-            const size_t vi = edgeNode.vi;
-            //      #pragma omp parallel for schedule(static) num_threads(6)
-            for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
-                const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi]
-                                                    && fixedVertices.at(edgeNode.vi).contains(dofi);
-                if (!isDofConstrainedFor_ev) {
-                    DifferentiateWithRespectTo dui{ev, dofi};
-                    // Axial force computation
-                    const double e_k = element.length - element.initialLength;
-                    const double e_kDeriv = computeDerivativeElementLength(e, dui);
-                    const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial;
+                const VertexPointer &rev_j = edgeNode.referenceElement->cV(0);
+                const VertexPointer &rev_jplus1 = edgeNode.referenceElement->cV(1);
+                // refElemOtherVertex can be precomputed
+                const VertexPointer &refElemOtherVertex = rev_j == &ev ? rev_jplus1 : rev_j;
+                const Node &refElemOtherVertexNode = pMesh->nodes[refElemOtherVertex];
+                if (edgeNode.referenceElement != &e) {
+                    internalForcesContributionFromThisEdge[evi + 2].first = refElemOtherVertexNode.vi;
+                }
+                const size_t vi = edgeNode.vi;
+                //      #pragma omp parallel for schedule(static) num_threads(6)
+                for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
+                    const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi]
+                                                        && fixedVertices.at(edgeNode.vi)
+                                                               .contains(dofi);
+                    if (!isDofConstrainedFor_ev) {
+                        DifferentiateWithRespectTo dui{ev, dofi};
+                        // Axial force computation
+                        const double e_k = element.length - element.initialLength;
+                        const double e_kDeriv = computeDerivativeElementLength(e, dui);
+                        const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial;
 
 #ifdef POLYSCOPE_DEFINED
                         if (std::isnan(axialForce_dofi)) {
@@ -976,7 +985,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
                         internalForcesContributionFromThisEdge[evi].second[dofi]
                             = axialForce_dofi + firstBendingForce_dofi + secondBendingForce_dofi
                               + torsionalForce_dofi;
-                }
+                    }
                     if (edgeNode.referenceElement != &e) {
                         const bool isDofConstrainedFor_refElemOtherVertex
                             = isVertexConstrained[refElemOtherVertexNode.vi]
@@ -1064,7 +1073,6 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
         }
 #endif
     }
-    pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm;
     pMesh->totalResidualForcesNorm = totalResidualForcesNorm;
     pMesh->averageResidualForcesNorm = totalResidualForcesNorm / pMesh->VN();
     //    pMesh->averageResidualForcesNorm = sumOfResidualForces.norm() / pMesh->VN();
@@ -1253,14 +1261,15 @@ void DRMSimulationModel::updateNodalMasses()
         for (const EdgePointer &ep : pMesh->nodes[v].incidentElements) {
             const size_t ei = pMesh->getIndex(ep);
             const Element &elem = pMesh->elements[ep];
-            const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length;
+            const double SkTranslational = elem.material.youngsModulus * elem.dimensions.A
+                                           / elem.length;
             translationalSumSk += SkTranslational;
             const double lengthToThe3 = std::pow(elem.length, 3);
-            const long double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2
+            const long double SkRotational_I2 = elem.material.youngsModulus * elem.dimensions.inertia.I2
                                                 / lengthToThe3;
-            const long double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3
+            const long double SkRotational_I3 = elem.material.youngsModulus * elem.dimensions.inertia.I3
                                                 / lengthToThe3;
-            const long double SkRotational_J = elem.material.youngsModulus * elem.inertia.J
+            const long double SkRotational_J = elem.material.youngsModulus * elem.dimensions.inertia.J
                                                / lengthToThe3;
             rotationalSumSk_I2 += SkRotational_I2;
             rotationalSumSk_I3 += SkRotational_I3;
@@ -1619,14 +1628,15 @@ void DRMSimulationModel::updatePositionsOnTheFly(
         double rotationalSumSk_J = 0;
         for (const EdgePointer &ep : pMesh->nodes[v].incidentElements) {
             const Element &elem = pMesh->elements[ep];
-            const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length;
+            const double SkTranslational = elem.material.youngsModulus * elem.dimensions.A
+                                           / elem.length;
             translationalSumSk += SkTranslational;
             const double lengthToThe3 = std::pow(elem.length, 3);
-            const double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2
+            const double SkRotational_I2 = elem.material.youngsModulus * elem.dimensions.inertia.I2
                                            / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
-            const double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3
+            const double SkRotational_I3 = elem.material.youngsModulus * elem.dimensions.inertia.I3
                                            / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
-            const double SkRotational_J = elem.material.youngsModulus * elem.inertia.J
+            const double SkRotational_J = elem.material.youngsModulus * elem.dimensions.inertia.J
                                           / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
             rotationalSumSk_I2 += SkRotational_I2;
             rotationalSumSk_I3 += SkRotational_I3;
@@ -1823,10 +1833,13 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptr<Simul
             results.debug_q_normal[vi] = q_normal;
             results.debug_q_nr[vi] = q_nr_nInit;
             results.rotationalDisplacementQuaternion[vi]
+                //Eigen::Quaterniond R
                 = (q_normal
                    * (q_f1_nInit * q_nr_nInit)); //q_f1_nDeformed * q_nr_nDeformed * q_normal;
             //Update the displacement vector to contain the euler angles
-            const Eigen::Vector3d eulerAngles = results.rotationalDisplacementQuaternion[vi]
+            const Eigen::Vector3d eulerAngles = results
+                                                    .rotationalDisplacementQuaternion[vi]
+                                                    //            R
                                                     .toRotationMatrix()
                                                     .eulerAngles(0, 1, 2);
             results.displacements[vi][3] = eulerAngles[0];
@@ -1876,7 +1889,7 @@ void DRMSimulationModel::printCurrentState() const
 void DRMSimulationModel::printDebugInfo() const
 {
     std::cout << pMesh->elements[0].rigidity.toString() << std::endl;
-    std::cout << "Number of dampings:" << numOfDampings << endl;
+    std::cout << "Number of dampings:" << numOfDampings << std::endl;
     printCurrentState();
 }
 
@@ -1999,10 +2012,10 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
     std::vector<double> I3(pMesh->EN());
     for (const EdgeType &e : pMesh->edge) {
         const size_t ei = pMesh->getIndex(e);
-        A[ei] = pMesh->elements[e].A;
-        J[ei] = pMesh->elements[e].inertia.J;
-        I2[ei] = pMesh->elements[e].inertia.I2;
-        I3[ei] = pMesh->elements[e].inertia.I3;
+        A[ei] = pMesh->elements[e].dimensions.A;
+        J[ei] = pMesh->elements[e].dimensions.inertia.J;
+        I2[ei] = pMesh->elements[e].dimensions.inertia.I2;
+        I3[ei] = pMesh->elements[e].dimensions.inertia.I3;
     }
 
     polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A);
@@ -2656,54 +2669,119 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
     return results;
 }
 
-void DRMSimulationModel::Settings::save(const std::filesystem::path &folderPath) const
-{
-    std::filesystem::create_directories(folderPath);
-    nlohmann::json json;
+//#ifdef USE_ENSMALLEN
 
-    json[jsonLabels.gamma] = gamma;
-    json[jsonLabels.shouldDraw] = shouldDraw ? "true" : "false";
-    json[jsonLabels.shouldCreatePlots] = shouldCreatePlots ? "true" : "false";
-    json[jsonLabels.beVerbose] = beVerbose ? "true" : "false";
-    json[jsonLabels.Dtini] = Dtini;
-    json[jsonLabels.xi] = xi;
+//void DRMSimulationModel::setJob(const std::shared_ptr<SimulationJob> &pJob)
+//{
+//    reset(pJob);
+
+//    updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
+//    //    arma::mat externalForces(pMesh->VN() * NumDoF);
+//    //    for (int vi = 0; vi < pMesh->VN(); vi++) {
+//    //        const Vector6d &nodeExternalForce = pMesh.vert[vi].node.externalForce;
+//    //        for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
+//    //            externalForces(vi * NumDoF + dofi) = nodeExternalForce[dofi];
+//    //        }
+//    //    }
+//}
+
+//SimulationMesh *DRMSimulationModel::getDeformedMesh(const arma::mat &x)
+//{
+//    for (int vi = 0; vi < pMesh->VN(); vi++) {
+//        Node &node = pMesh->nodes[vi];
+//        for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
+//            node.displacements[dofi] = x(vi * DoF::NumDoF + dofi);
+//        }
+//    }
+
+//    applyDisplacements(constrainedVertices);
+//    if (!pJob->nodalForcedDisplacements.empty()) {
+//        applyForcedDisplacements(pJob->nodalForcedDisplacements);
+//    }
+
+//    return pMesh.get();
+//}
+
+//double DRMSimulationModel::EvaluateWithGradient(const arma::mat &x, arma::mat &g)
+//{
+//    for (int vi = 0; vi < pMesh->VN(); vi++) {
+//        Node &node = pMesh->nodes[vi];
+//        for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
+//            node.displacements[dofi] = x(vi * DoF::NumDoF + dofi);
+//        }
+//    }
+
+//    applyDisplacements(constrainedVertices);
+//    if (!pJob->nodalForcedDisplacements.empty()) {
+//        applyForcedDisplacements(pJob->nodalForcedDisplacements);
+//    }
+//    updateElementalLengths();
+//    updateNormalDerivatives();
+//    updateT1Derivatives();
+//    updateRDerivatives();
+//    updateT2Derivatives();
+//    updateT3Derivatives();
+//    updateResidualForcesOnTheFly(constrainedVertices);
+//    for (int vi = 0; vi < pMesh->VN(); vi++) {
+//        Node &node = pMesh->nodes[vi];
+//        for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
+//            g(vi * NumDoF + dofi) = node.force.residual[dofi];
+//        }
+//    }
+//    mCurrentSimulationStep++;
+//    const double PE = computeTotalPotentialEnergy();
+//    std::cout << "PE:" << PE << std::endl;
+//    return PE;
+//}
+
+//#endif
+
+void DRMSimulationModel::Settings::save(const std::filesystem::path &jsonFilePath) const
+{
+    nlohmann::json json;
+    json[GET_VARIABLE_NAME(shouldDraw)] = shouldDraw;
+    json[GET_VARIABLE_NAME(beVerbose)] = beVerbose;
+    json[GET_VARIABLE_NAME(shouldCreatePlots)] = shouldCreatePlots;
+    json[GET_VARIABLE_NAME(Dtini)] = Dtini;
+    json[GET_VARIABLE_NAME(xi)] = xi;
+    json[GET_VARIABLE_NAME(gamma)] = gamma;
+    json[GET_VARIABLE_NAME(useKineticDamping)] = totalResidualForcesNormThreshold;
     if (maxDRMIterations.has_value()) {
-        json[jsonLabels.maxDRMIterations] = maxDRMIterations.value();
+        json[GET_VARIABLE_NAME(jsonLabels.maxDRMIterations)] = maxDRMIterations.value();
     }
     if (debugModeStep.has_value()) {
-        json[jsonLabels.debugModeStep] = debugModeStep.value();
+        json[GET_VARIABLE_NAME(debugModeStep)] = debugModeStep.value();
     }
     if (totalTranslationalKineticEnergyThreshold.has_value()) {
-        json[jsonLabels.totalTranslationalKineticEnergyThreshold]
+        json[GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold)]
             = totalTranslationalKineticEnergyThreshold.value();
     }
-
     if (averageResidualForcesCriterionThreshold.has_value()) {
-        json[jsonLabels.averageResidualForcesCriterionThreshold]
-            = averageResidualForcesCriterionThreshold.value();
-    }
-
-    if (averageResidualForcesCriterionThreshold.has_value()) {
-        json[jsonLabels.averageResidualForcesCriterionThreshold]
+        json[GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold)]
             = averageResidualForcesCriterionThreshold.value();
     }
     if (linearGuessForceScaleFactor.has_value()) {
-        json[jsonLabels.linearGuessForceScaleFactor] = linearGuessForceScaleFactor.value();
+        json[GET_VARIABLE_NAME(linearGuessForceScaleFactor)] = linearGuessForceScaleFactor.value();
     }
-
+    if (viscousDampingFactor.has_value()) {
+        json[GET_VARIABLE_NAME(viscousDampingFactor)] = viscousDampingFactor.value();
+    }
+    json[GET_VARIABLE_NAME(useKineticDamping)] = useKineticDamping;
     //    if (intermediateResultsSaveStep.has_value()) {
     //        json[GET_VARIABLE_NAME(intermediateResultsSaveStep)] = intermediateResultsSaveStep.value();
     //    }
 
-    if (saveIntermediateBestStates.has_value()) {
-        json[GET_VARIABLE_NAME(saveIntermediateBestStates)] = saveIntermediateBestStates.value()
-                                                                  ? "true"
-                                                                  : "false";
-    }
+    //    if (saveIntermediateBestStates.has_value()) {
+    //        json[GET_VARIABLE_NAME(saveIntermediateBestStates)] = saveIntermediateBestStates.value()
+    //                                                                  ? "true"
+    //                                                                  : "false";
+    //    }
 
-    const std::string jsonFilename = "drmSettings.json";
-    std::ofstream jsonFile(std::filesystem::path(folderPath).append(jsonFilename));
+    std::filesystem::create_directories(jsonFilePath.parent_path());
+    std::ofstream jsonFile(jsonFilePath);
+    std::cout << "Saving DRM settings to:" << jsonFilePath << std::endl;
     jsonFile << json;
+    jsonFile.close();
 }
 
 bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePath)
@@ -2727,41 +2805,40 @@ bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePat
     std::ifstream ifs(jsonFilePath.string());
     ifs >> json;
 
-    if (json.contains(jsonLabels.shouldDraw)) {
-        shouldDraw = json.at(jsonLabels.shouldDraw) == "true" ? true : false;
+    if (json.contains(GET_VARIABLE_NAME(shouldDraw))) {
+        shouldDraw = json.at(GET_VARIABLE_NAME(shouldDraw)) == "true" ? true : false;
     }
-    if (json.contains(jsonLabels.beVerbose)) {
-        beVerbose = json.at(jsonLabels.beVerbose) == "true" ? true : false;
+    if (json.contains(GET_VARIABLE_NAME(beVerbose))) {
+        beVerbose = json.at(GET_VARIABLE_NAME(beVerbose)) == "true" ? true : false;
     }
-    if (json.contains(jsonLabels.shouldCreatePlots)) {
-        shouldCreatePlots = json.at(jsonLabels.shouldCreatePlots) == "true" ? true : false;
+    if (json.contains(GET_VARIABLE_NAME(shouldCreatePlots))) {
+        shouldCreatePlots = json.at(GET_VARIABLE_NAME(shouldCreatePlots)) == "true" ? true : false;
     }
-    if (json.contains(jsonLabels.Dtini)) {
-        Dtini = json.at(jsonLabels.Dtini);
+    if (json.contains(GET_VARIABLE_NAME(Dtini))) {
+        Dtini = json.at(GET_VARIABLE_NAME(Dtini));
     }
-    if (json.contains(jsonLabels.xi)) {
-        xi = json.at(jsonLabels.xi);
+    if (json.contains(GET_VARIABLE_NAME(xi))) {
+        xi = json.at(GET_VARIABLE_NAME(xi));
     }
-    if (json.contains(jsonLabels.maxDRMIterations)) {
-        maxDRMIterations = json.at(jsonLabels.maxDRMIterations);
+    if (json.contains(GET_VARIABLE_NAME(maxDRMIterations))) {
+        maxDRMIterations = json.at(GET_VARIABLE_NAME(maxDRMIterations));
     }
-    if (json.contains(jsonLabels.debugModeStep)) {
-        debugModeStep = json.at(jsonLabels.debugModeStep);
+    if (json.contains(GET_VARIABLE_NAME(debugModeStep))) {
+        debugModeStep = json.at(GET_VARIABLE_NAME(debugModeStep));
     }
-    if (json.contains(jsonLabels.totalTranslationalKineticEnergyThreshold)) {
+    if (json.contains(GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold))) {
         totalTranslationalKineticEnergyThreshold = json.at(
-            jsonLabels.totalTranslationalKineticEnergyThreshold);
+            GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold));
     }
-    if (json.contains(jsonLabels.gamma)) {
-        gamma = json.at(jsonLabels.gamma);
+    if (json.contains(GET_VARIABLE_NAME(gamma))) {
+        gamma = json.at(GET_VARIABLE_NAME(gamma));
     }
-    if (json.contains(jsonLabels.averageResidualForcesCriterionThreshold)) {
+    if (json.contains(GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold))) {
         averageResidualForcesCriterionThreshold = json.at(
-            jsonLabels.averageResidualForcesCriterionThreshold);
+            GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold));
     }
-
-    if (json.contains(jsonLabels.linearGuessForceScaleFactor)) {
-        linearGuessForceScaleFactor = json.at(jsonLabels.linearGuessForceScaleFactor);
+    if (json.contains(GET_VARIABLE_NAME(linearGuessForceScaleFactor))) {
+        linearGuessForceScaleFactor = json.at(GET_VARIABLE_NAME(linearGuessForceScaleFactor));
     }
 
     //    if (json.contains(GET_VARIABLE_NAME(intermediateResultsSaveStep))) {
@@ -2774,6 +2851,7 @@ bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePat
                                          ? true
                                          : false;
     }
+    //    std::cout << json.dump() << std::endl;
 
     return true;
 }

drmsimulationmodel.hpp

@@ -1,13 +1,21 @@
 #ifndef BEAMFORMFINDER_HPP
 #define BEAMFORMFINDER_HPP
 
-#include "matplot/matplot.h"
+//#ifdef USE_MATPLOT
+//#include "matplot.h"
+#include <matplot/matplot.h>
+//#endif
 #include "simulationmesh.hpp"
 #include "simulationmodel.hpp"
 #include <Eigen/Dense>
 #include <filesystem>
 #include <unordered_set>
 
+#ifdef USE_ENSMALLEN
+#include <armadillo>
+#include <ensmallen.hpp>
+#endif
+
 struct SimulationJob;
 
 enum DoF { Ux = 0, Uy, Uz, Nx, Ny, Nr, NumDoF };
@@ -25,7 +33,6 @@ class DRMSimulationModel : public SimulationModel
 public:
     struct Settings
     {
-        //        bool isDebugMode{false};
         bool shouldDraw{false};
         bool beVerbose{false};
         bool shouldCreatePlots{false};
@@ -38,8 +45,8 @@ public:
         int gradualForcedDisplacementSteps{50};
         //        int desiredGradualExternalLoadsSteps{1};
         double gamma{0.8};
-        double totalResidualForcesNormThreshold{1e-3};
-        double totalExternalForcesNormPercentageTermination{1e-3};
+        double totalResidualForcesNormThreshold{1e-8};
+        double totalExternalForcesNormPercentageTermination{1e-8};
         std::optional<int> maxDRMIterations;
         std::optional<int> debugModeStep;
         std::optional<double> totalTranslationalKineticEnergyThreshold;
@@ -47,11 +54,11 @@ public:
         std::optional<double> linearGuessForceScaleFactor;
         //        std::optional<int> intermediateResultsSaveStep;
         std::optional<bool> saveIntermediateBestStates;
-        Settings() {}
-        void save(const std::filesystem::path &folderPath) const;
-        bool load(const std::filesystem::path &filePath);
         std::optional<double> viscousDampingFactor;
         bool useKineticDamping{true};
+        Settings() {}
+        void save(const std::filesystem::path &jsonFilePath) const;
+        bool load(const std::filesystem::path &filePath);
         struct JsonLabels
         {
             const std::string shouldDraw{"shouldDraw"};
@@ -59,14 +66,16 @@ public:
             const std::string shouldCreatePlots{"shouldCreatePlots"};
             const std::string Dtini{"DtIni"};
             const std::string xi{"xi"};
+            const std::string gamma{"gamma"};
+            const std::string totalResidualForcesNormThreshold;
             const std::string maxDRMIterations{"maxIterations"};
             const std::string debugModeStep{"debugModeStep"};
             const std::string totalTranslationalKineticEnergyThreshold{
                 "totalTranslationaKineticEnergyThreshold"};
-            const std::string gamma{"gamma"};
             const std::string averageResidualForcesCriterionThreshold{
                 "averageResidualForcesThreshold"};
             const std::string linearGuessForceScaleFactor{"linearGuessForceScaleFactor"};
+            const std::string viscousDampingFactor{"viscousDampingFactor"};
         };
         static JsonLabels jsonLabels;
     };
@@ -94,7 +103,9 @@ private:
     //  Eigen::Tensor<double, 4> theta3Derivatives;
     //  std::unordered_map<MyKeyType, double, key_hash> theta3Derivatives;
     bool shouldApplyInitialDistortion = false;
-
+    //#ifdef USE_ENSMALLEN
+    //    std::shared_ptr<SimulationJob> pJob;
+    //#endif
     void reset(const std::shared_ptr<SimulationJob> &pJob, const Settings &settings);
     void updateNodalInternalForces(
         const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
@@ -142,7 +153,7 @@ private:
         const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
 
 #ifdef POLYSCOPE_DEFINED
-  void draw(const string &screenshotsFolder= {});
+    void draw(const std::string &screenshotsFolder = {});
 #endif
   void
   updateNodalInternalForce(Vector6d &nodalInternalForce,
@@ -241,11 +252,19 @@ private:
 
   virtual SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob);
 
+  void reset(const std::shared_ptr<SimulationJob> &pJob);
+
   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
 
   static void runUnitTests();
 };

edgemesh.cpp

@@ -25,13 +25,13 @@ Eigen::MatrixX3d VCGEdgeMesh::getEigenEdgeNormals() const {
   return eigenEdgeNormals;
 }
 
-bool VCGEdgeMesh::save(const string &plyFilename)
+bool VCGEdgeMesh::save(const std::filesystem::__cxx11::path &meshFilePath)
 {
-    std::string filename = plyFilename;
+    std::string filename = meshFilePath;
     if (filename.empty()) {
         filename = std::filesystem::current_path().append(getLabel() + ".ply").string();
-    } else if (std::filesystem::is_directory(std::filesystem::path(plyFilename))) {
-        filename = std::filesystem::path(plyFilename).append(getLabel() + ".ply").string();
+    } else if (std::filesystem::is_directory(std::filesystem::path(meshFilePath))) {
+        filename = std::filesystem::path(meshFilePath).append(getLabel() + ".ply").string();
     }
     assert(std::filesystem::path(filename).extension().string() == ".ply");
     unsigned int mask = 0;
@@ -197,36 +197,36 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth,
   return true;
 }
 
-bool VCGEdgeMesh::load(const string &plyFilename) {
+bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath)
+{
+    std::string usedPath = meshFilePath;
+    if (std::filesystem::path(meshFilePath).is_relative()) {
+        usedPath = std::filesystem::absolute(meshFilePath).string();
+    }
+    assert(std::filesystem::exists(usedPath));
+    Clear();
+    //  if (!loadUsingNanoply(usedPath)) {
+    //    std::cerr << "Error: Unable to open " + usedPath << std::endl;
+    //    return false;
+    //  }
+    if (!loadUsingDefaultLoader(usedPath)) {
+        std::cerr << "Error: Unable to open " + usedPath << std::endl;
+        return false;
+    }
 
-  std::string usedPath = plyFilename;
-  if (std::filesystem::path(plyFilename).is_relative()) {
-    usedPath = std::filesystem::absolute(plyFilename).string();
-  }
-  assert(std::filesystem::exists(usedPath));
-  this->Clear();
-  //  if (!loadUsingNanoply(usedPath)) {
-  //    std::cerr << "Error: Unable to open " + usedPath << std::endl;
-  //    return false;
-  //  }
-  if (!loadUsingDefaultLoader(usedPath)) {
-      std::cerr << "Error: Unable to open " + usedPath << std::endl;
-      return false;
-  }
+    getEdges(eigenEdges);
+    getVertices(eigenVertices);
+    vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
+    label = std::filesystem::path(meshFilePath).stem().string();
 
-  getEdges(eigenEdges);
-  getVertices(eigenVertices);
-  vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
-  label = std::filesystem::path(plyFilename).stem().string();
+    const bool printDebugInfo = false;
+    if (printDebugInfo) {
+        std::cout << meshFilePath << " was loaded successfuly." << std::endl;
+        std::cout << "Mesh has " << EN() << " edges." << std::endl;
+    }
 
-  const bool printDebugInfo = false;
-  if (printDebugInfo) {
-    std::cout << plyFilename << " was loaded successfuly." << std::endl;
-    std::cout << "Mesh has " << EN() << " edges." << std::endl;
-  }
-
-  label=std::filesystem::path(plyFilename).stem().string();
-  return true;
+    label = std::filesystem::path(meshFilePath).stem().string();
+    return true;
 }
 
 //bool VCGEdgeMesh::loadUsingNanoply(const std::string &plyFilename) {
@@ -248,7 +248,7 @@ bool VCGEdgeMesh::load(const string &plyFilename) {
 
 bool VCGEdgeMesh::loadUsingDefaultLoader(const std::string &plyFilename)
 {
-    this->Clear();
+    Clear();
     //  assert(plyFileHasAllRequiredFields(plyFilename));
     // Load the ply file
     int mask = 0;
@@ -337,10 +337,35 @@ bool VCGEdgeMesh::copy(VCGEdgeMesh &mesh) {
   return true;
 }
 
+void VCGEdgeMesh::set(const std::vector<double> &vertexPositions, const std::vector<int> &edges)
+{
+    Clear();
+    for (int ei = 0; ei < edges.size(); ei += 2) {
+        const int vi0 = edges[ei];
+        const int vi1 = edges[ei + 1];
+        const int vi0_offset = 3 * vi0;
+        const int vi1_offset = 3 * vi1;
+        const CoordType p0(vertexPositions[vi0_offset],
+                           vertexPositions[vi0_offset + 1],
+                           vertexPositions[vi0_offset + 2]);
+        const CoordType p1(vertexPositions[vi1_offset],
+                           vertexPositions[vi1_offset + 1],
+                           vertexPositions[vi1_offset + 2]);
+        auto eIt = vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(*this, p0, p1);
+        CoordType n(0, 0, 1);
+        eIt->cV(0)->N() = n;
+        eIt->cV(1)->N() = n;
+    }
+    removeDuplicateVertices();
+
+    updateEigenEdgeAndVertices();
+}
+
 void VCGEdgeMesh::removeDuplicateVertices()
 {
     vcg::tri::Clean<VCGEdgeMesh>::MergeCloseVertex(*this, 0.000061524);
     vcg::tri::Allocator<VCGEdgeMesh>::CompactVertexVector(*this);
+    vcg::tri::Allocator<VCGEdgeMesh>::CompactEdgeVector(*this);
     vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
 }
 
@@ -418,8 +443,12 @@ polyscope::CurveNetwork *VCGEdgeMesh::registerForDrawing(
 {
     PolyscopeInterface::init();
     const double drawingRadius = desiredRadius;
+    updateEigenEdgeAndVertices();
     polyscope::CurveNetwork *polyscopeHandle_edgeMesh
         = polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges());
+    //    std::cout << "EDGES:" << polyscopeHandle_edgeMesh->nEdges() << std::endl;
+    assert(polyscopeHandle_edgeMesh->nEdges() == getEigenEdges().rows()
+           && polyscopeHandle_edgeMesh->nNodes() == getEigenVertices().rows());
 
     polyscopeHandle_edgeMesh->setEnabled(shouldEnable);
     polyscopeHandle_edgeMesh->setRadius(drawingRadius, true);

edgemesh.hpp

@@ -58,9 +58,12 @@ public:
              * */
     bool copy(VCGEdgeMesh &mesh);
 
+    void set(const std::vector<double> &vertexPositions, const std::vector<int> &edges);
+
     void removeDuplicateVertices();
-    void deleteDanglingVertices();
-    void deleteDanglingVertices(vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu);
+    virtual void deleteDanglingVertices();
+    virtual void deleteDanglingVertices(
+        vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu);
 
     void getEdges(Eigen::MatrixX3d &edgeStartingPoints, Eigen::MatrixX3d &edgeEndingPoints) const;
 
@@ -70,8 +73,8 @@ public:
 
     //    bool loadUsingNanoply(const std::string &plyFilename);
 
-    bool load(const std::string &plyFilename) override;
-    bool save(const std::string &plyFilename = std::string()) override;
+    bool load(const std::filesystem::path &meshFilePath) override;
+    bool save(const std::filesystem::path &meshFilePath = std::filesystem::path()) override;
 
     bool createSpanGrid(const size_t squareGridDimension);
     bool createSpanGrid(const size_t desiredWidth, const size_t desiredHeight);

hexagonremesher.hpp

@@ -11,7 +11,7 @@
 
 namespace PolygonalRemeshing {
 //std::shared_ptr<VCGPolyMesh>
-std::shared_ptr<VCGPolyMesh> remeshWithPolygons(VCGTriMesh &surface)
+inline std::shared_ptr<VCGPolyMesh> remeshWithPolygons(VCGTriMesh &surface)
 {
     //    const std::string tileIntoFilePath(
     //        "/home/iason/Coding/build/PatternTillingReducedModel/RelWithDebInfo/plane.ply");

linearsimulationmodel.cpp

@@ -0,0 +1,274 @@
+#include "linearsimulationmodel.hpp"
+//#include "gsl/gsl"
+
+std::vector<fea::Elem> LinearSimulationModel::getFeaElements(
+    const std::shared_ptr<SimulationMesh> &pMesh)
+{
+    const int numberOfEdges = pMesh->EN();
+    std::vector<fea::Elem> elements(numberOfEdges);
+    for (int edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
+        const SimulationMesh::CoordType &evn0 = pMesh->edge[edgeIndex].cV(0)->cN();
+        const SimulationMesh::CoordType &evn1 = pMesh->edge[edgeIndex].cV(1)->cN();
+        const std::vector<double> nAverageVector{(evn0[0] + evn1[0]) / 2,
+                                                 (evn0[1] + evn1[1]) / 2,
+                                                 (evn0[2] + evn1[2]) / 2};
+        const Element &element = pMesh->elements[edgeIndex];
+        const double E = element.material.youngsModulus;
+        fea::Props feaProperties(E * element.dimensions.A,
+                                 E * element.dimensions.inertia.I3,
+                                 E * element.dimensions.inertia.I2,
+                                 element.material.G * element.dimensions.inertia.J,
+                                 nAverageVector);
+        const int vi0 = pMesh->getIndex(pMesh->edge[edgeIndex].cV(0));
+        const int vi1 = pMesh->getIndex(pMesh->edge[edgeIndex].cV(1));
+        elements[edgeIndex] = fea::Elem(vi0, vi1, feaProperties);
+    }
+
+    return elements;
+}
+
+std::vector<fea::Node> LinearSimulationModel::getFeaNodes(
+    const std::shared_ptr<SimulationMesh> &pMesh)
+{
+    const int numberOfNodes = pMesh->VN();
+    std::vector<fea::Node> feaNodes(numberOfNodes);
+    for (int vi = 0; vi < numberOfNodes; vi++) {
+        const CoordType &p = pMesh->vert[vi].cP();
+        feaNodes[vi] = fea::Node(p[0], p[1], p[2]);
+    }
+
+    return feaNodes;
+}
+
+std::vector<fea::BC> LinearSimulationModel::getFeaFixedNodes(
+    const std::shared_ptr<SimulationJob> &job)
+{
+    std::vector<fea::BC> boundaryConditions;
+    //        boundaryConditions.reserve(job->constrainedVertices.size() * 6
+    //                                   + job->nodalForcedDisplacements.size() * 3);
+    for (auto fixedVertex : job->constrainedVertices) {
+        const int vertexIndex = fixedVertex.first;
+        for (int dofIndex : fixedVertex.second) {
+            boundaryConditions.emplace_back(
+                fea::BC(vertexIndex, static_cast<fea::DOF>(dofIndex), 0));
+        }
+    }
+
+    for (auto forcedDisplacement : job->nodalForcedDisplacements) {
+        const int vi = forcedDisplacement.first;
+        for (int dofIndex = 0; dofIndex < 3; dofIndex++) {
+            boundaryConditions.emplace_back(
+                fea::BC(vi, static_cast<fea::DOF>(dofIndex), forcedDisplacement.second[dofIndex]));
+        }
+    }
+
+    return boundaryConditions;
+}
+
+std::vector<fea::Force> LinearSimulationModel::getFeaNodalForces(
+    const std::shared_ptr<SimulationJob> &job)
+{
+    std::vector<fea::Force> nodalForces;
+    nodalForces.reserve(job->nodalExternalForces.size() * 6);
+
+    for (auto nodalForce : job->nodalExternalForces) {
+        for (int dofIndex = 0; dofIndex < 6; dofIndex++) {
+            if (nodalForce.second[dofIndex] == 0) {
+                continue;
+            }
+            fea::Force f(nodalForce.first, dofIndex, nodalForce.second[dofIndex]);
+            nodalForces.emplace_back(f);
+        }
+    }
+
+    return nodalForces;
+}
+
+SimulationResults LinearSimulationModel::getResults(
+    const fea::Summary &feaSummary, const std::shared_ptr<SimulationJob> &simulationJob)
+{
+    SimulationResults results;
+    results.converged = feaSummary.converged;
+    if (!results.converged) {
+        return results;
+    }
+
+    results.executionTime = feaSummary.total_time_in_ms * 1000;
+    // displacements
+    results.displacements.resize(feaSummary.num_nodes);
+    results.rotationalDisplacementQuaternion.resize(feaSummary.num_nodes);
+    for (int vi = 0; vi < feaSummary.num_nodes; vi++) {
+        results.displacements[vi] = Vector6d(feaSummary.nodal_displacements[vi]);
+
+        const Vector6d &nodalDisplacement = results.displacements[vi];
+        Eigen::Quaternion<double> q_nx;
+        q_nx = Eigen::AngleAxis<double>(nodalDisplacement[3], Eigen::Vector3d(1, 0, 0));
+        Eigen::Quaternion<double> q_ny;
+        q_ny = Eigen::AngleAxis<double>(nodalDisplacement[4], Eigen::Vector3d(0, 1, 0));
+        Eigen::Quaternion<double> q_nz;
+        q_nz = Eigen::AngleAxis<double>(nodalDisplacement[5], Eigen::Vector3d(0, 0, 1));
+        results.rotationalDisplacementQuaternion[vi] = q_nx * q_ny * q_nz;
+    }
+    results.setLabelPrefix("Linear");
+
+    //  // Convert forces
+    //  // Convert to vector of eigen matrices of the form force component-> per
+    //  // Edge
+    //  const int numDof = 6;
+    //  const size_t numberOfEdges = feaSummary.element_forces.size();
+    //  edgeForces =
+    //      std::vector<Eigen::VectorXd>(numDof, Eigen::VectorXd(2 *
+    //      numberOfEdges));
+    //  for (gsl::index edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
+    //    for (gsl::index forceComponentIndex = 0; forceComponentIndex < numDof;
+    //         forceComponentIndex++) {
+    //      (edgeForces[forceComponentIndex])(2 * edgeIndex) =
+    //          feaSummary.element_forces[edgeIndex][forceComponentIndex];
+    //      (edgeForces[forceComponentIndex])(2 * edgeIndex + 1) =
+    //          feaSummary.element_forces[edgeIndex][numDof +
+    //          forceComponentIndex];
+    //    }
+    //  }
+
+    for (int ei = 0; ei < feaSummary.num_elems; ei++) {
+        const std::vector<double> &elementForces = feaSummary.element_forces[ei];
+        const Element &element = simulationJob->pMesh->elements[ei];
+        //Axial
+        const double elementPotentialEnergy_axial_v0 = std::pow(elementForces[0], 2)
+                                                       * element.initialLength
+                                                       / (4 * element.material.youngsModulus
+                                                          * element.dimensions.A);
+        const double elementPotentialEnergy_axial_v1 = std::pow(elementForces[6], 2)
+                                                       * element.initialLength
+                                                       / (4 * element.material.youngsModulus
+                                                          * element.dimensions.A);
+        const double elementPotentialEnergy_axial = elementPotentialEnergy_axial_v0
+                                                    + elementPotentialEnergy_axial_v1;
+        //Shear
+        const double elementPotentialEnergy_shearY_v0 = std::pow(elementForces[1], 2)
+                                                        * element.initialLength
+                                                        / (4 * element.dimensions.A
+                                                           * element.material.G);
+        const double elementPotentialEnergy_shearZ_v0 = std::pow(elementForces[2], 2)
+                                                        * element.initialLength
+                                                        / (4 * element.dimensions.A
+                                                           * element.material.G);
+        const double elementPotentialEnergy_shearY_v1 = std::pow(elementForces[7], 2)
+                                                        * element.initialLength
+                                                        / (4 * element.dimensions.A
+                                                           * element.material.G);
+        const double elementPotentialEnergy_shearZ_v1 = std::pow(elementForces[8], 2)
+                                                        * element.initialLength
+                                                        / (4 * element.dimensions.A
+                                                           * element.material.G);
+        const double elementPotentialEnergy_shear = elementPotentialEnergy_shearY_v0
+                                                    + elementPotentialEnergy_shearZ_v0
+                                                    + elementPotentialEnergy_shearY_v1
+                                                    + elementPotentialEnergy_shearZ_v1;
+        //Bending
+        const double elementPotentialEnergy_bendingY_v0 = std::pow(elementForces[4], 2)
+                                                          * element.initialLength
+                                                          / (4 * element.material.youngsModulus
+                                                             * element.dimensions.inertia.I2);
+        const double elementPotentialEnergy_bendingZ_v0 = std::pow(elementForces[5], 2)
+                                                          * element.initialLength
+                                                          / (4 * element.material.youngsModulus
+                                                             * element.dimensions.inertia.I3);
+        const double elementPotentialEnergy_bendingY_v1 = std::pow(elementForces[10], 2)
+                                                          * element.initialLength
+                                                          / (4 * element.material.youngsModulus
+                                                             * element.dimensions.inertia.I2);
+        const double elementPotentialEnergy_bendingZ_v1 = std::pow(elementForces[11], 2)
+                                                          * element.initialLength
+                                                          / (4 * element.material.youngsModulus
+                                                             * element.dimensions.inertia.I3);
+        const double elementPotentialEnergy_bending = elementPotentialEnergy_bendingY_v0
+                                                      + elementPotentialEnergy_bendingZ_v0
+                                                      + elementPotentialEnergy_bendingY_v1
+                                                      + elementPotentialEnergy_bendingZ_v1;
+        //Torsion
+        const double elementPotentialEnergy_torsion_v0 = std::pow(elementForces[3], 2)
+                                                         * element.initialLength
+                                                         / (4 * element.dimensions.inertia.J
+                                                            * element.material.G);
+        const double elementPotentialEnergy_torsion_v1 = std::pow(elementForces[9], 2)
+                                                         * element.initialLength
+                                                         / (4 * element.dimensions.inertia.J
+                                                            * element.material.G);
+        const double elementPotentialEnergy_torsion = elementPotentialEnergy_torsion_v0
+                                                      + elementPotentialEnergy_torsion_v1;
+        const double elementInternalPotentialEnergy = elementPotentialEnergy_axial
+                                                      + elementPotentialEnergy_shear
+                                                      + elementPotentialEnergy_bending
+                                                      + elementPotentialEnergy_torsion;
+        results.internalPotentialEnergy += elementInternalPotentialEnergy;
+    }
+    results.pJob = simulationJob;
+
+    return results;
+}
+
+SimulationResults LinearSimulationModel::executeSimulation(
+    const std::shared_ptr<SimulationJob> &pSimulationJob)
+{
+    assert(pSimulationJob->pMesh->VN() != 0);
+    const bool wasInitializedWithDifferentStructure = pSimulationJob->pMesh->EN()
+                                                          != simulator.structure.elems.size()
+                                                      || pSimulationJob->pMesh->VN()
+                                                             != simulator.structure.nodes.size();
+    if (!simulator.wasInitialized() || wasInitializedWithDifferentStructure) {
+        setStructure(pSimulationJob->pMesh);
+    }
+    //  printInfo(job);
+
+    // create the default options
+    fea::Options opts;
+    opts.save_elemental_forces = false;
+    opts.save_nodal_displacements = false;
+    opts.save_nodal_forces = false;
+    opts.save_report = false;
+    opts.save_tie_forces = false;
+    //  if (!elementalForcesOutputFilepath.empty()) {
+    //    opts.save_elemental_forces = true;
+    //    opts.elemental_forces_filename = elementalForcesOutputFilepath;
+    //  }
+    //  if (!nodalDisplacementsOutputFilepath.empty()) {
+    //    opts.save_nodal_displacements = true;
+    //    opts.nodal_displacements_filename = nodalDisplacementsOutputFilepath;
+    //  }
+
+    // have the program output status updates
+    opts.verbose = false;
+
+    // form an empty vector of ties
+    std::vector<fea::Tie> ties;
+
+    // also create an empty list of equations
+    std::vector<fea::Equation> equations;
+    auto fixedVertices = getFeaFixedNodes(pSimulationJob);
+    auto nodalForces = getFeaNodalForces(pSimulationJob);
+    fea::Summary feaResults = simulator.solve(fixedVertices, nodalForces, ties, equations, opts);
+
+    SimulationResults results = getResults(feaResults, pSimulationJob);
+    results.pJob = pSimulationJob;
+    return results;
+}
+
+void LinearSimulationModel::setStructure(const std::shared_ptr<SimulationMesh> &pMesh)
+{
+    fea::BeamStructure structure(getFeaNodes(pMesh), getFeaElements(pMesh));
+    simulator.initialize(structure);
+}
+
+void LinearSimulationModel::printInfo(const fea::BeamStructure &job)
+{
+    std::cout << "Details regarding the fea::Job:" << std::endl;
+    std::cout << "Nodes:" << std::endl;
+    for (fea::Node n : job.nodes) {
+        std::cout << n << std::endl;
+    }
+    std::cout << "Elements:" << std::endl;
+    for (Eigen::Vector2i e : job.elems) {
+        std::cout << e << std::endl;
+    }
+}

linearsimulationmodel.hpp

@@ -12,265 +12,21 @@ public:
     LinearSimulationModel(){
 
     }
-    static std::vector<fea::Elem>
-    getFeaElements(const std::shared_ptr<SimulationJob> &job) {
-        const int numberOfEdges = job->pMesh->EN();
-        std::vector<fea::Elem> elements(numberOfEdges);
-        for (int edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
-                const SimulationMesh::CoordType &evn0 =
-                        job->pMesh->edge[edgeIndex].cV(0)->cN();
-                const SimulationMesh::CoordType &evn1 = job->pMesh->edge[edgeIndex].cV(1)->cN();
-                const std::vector<double> nAverageVector{(evn0[0] + evn1[0]) / 2,
-                                                         (evn0[1] + evn1[1]) / 2,
-                                                         (evn0[2] + evn1[2]) / 2};
-                const Element &element = job->pMesh->elements[edgeIndex];
-                const double E = element.material.youngsModulus;
-                fea::Props feaProperties(E * element.A,
-                                         E * element.inertia.I3,
-                                         E * element.inertia.I2,
-                                         element.G * element.inertia.J,
-                                         nAverageVector);
-                const int vi0 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(0));
-                const int vi1 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(1));
-                elements[edgeIndex] = fea::Elem(vi0, vi1, feaProperties);
-            }
+    static std::vector<fea::Elem> getFeaElements(const std::shared_ptr<SimulationMesh> &pMesh);
 
-        return elements;
-    }
+    static std::vector<fea::Node> getFeaNodes(const std::shared_ptr<SimulationMesh> &pMesh);
+    static std::vector<fea::BC> getFeaFixedNodes(const std::shared_ptr<SimulationJob> &job);
 
-    static std::vector<fea::Node>
-    getFeaNodes(const std::shared_ptr<SimulationJob> &job) {
-        const int numberOfNodes = job->pMesh->VN();
-        std::vector<fea::Node> feaNodes(numberOfNodes);
-        for (int vi = 0; vi < numberOfNodes; vi++) {
-                const CoordType &p = job->pMesh->vert[vi].cP();
-                feaNodes[vi] = fea::Node(p[0], p[1], p[2]);
-            }
-
-        return feaNodes;
-    }
-    static std::vector<fea::BC>
-    getFeaFixedNodes(const std::shared_ptr<SimulationJob> &job) {
-        std::vector<fea::BC> boundaryConditions;
-        //        boundaryConditions.reserve(job->constrainedVertices.size() * 6
-        //                                   + job->nodalForcedDisplacements.size() * 3);
-        for (auto fixedVertex : job->constrainedVertices) {
-                const int vertexIndex = fixedVertex.first;
-                for (int dofIndex : fixedVertex.second) {
-                    boundaryConditions.emplace_back(
-                        fea::BC(vertexIndex, static_cast<fea::DOF>(dofIndex), 0));
-                }
-        }
-
-        for (auto forcedDisplacement : job->nodalForcedDisplacements) {
-            const int vi = forcedDisplacement.first;
-            for (int dofIndex = 0; dofIndex < 3; dofIndex++) {
-                boundaryConditions.emplace_back(fea::BC(vi,
-                                                        static_cast<fea::DOF>(dofIndex),
-                                                        forcedDisplacement.second[dofIndex]));
-            }
-        }
-
-        return boundaryConditions;
-    }
-
-    static std::vector<fea::Force>
-    getFeaNodalForces(const std::shared_ptr<SimulationJob> &job) {
-        std::vector<fea::Force> nodalForces;
-        nodalForces.reserve(job->nodalExternalForces.size() * 6);
-
-        for (auto nodalForce : job->nodalExternalForces) {
-                for (int dofIndex = 0; dofIndex < 6; dofIndex++) {
-                        if (nodalForce.second[dofIndex] == 0) {
-                                continue;
-                            }
-                        fea::Force f(nodalForce.first, dofIndex, nodalForce.second[dofIndex]);
-                        nodalForces.emplace_back(f);
-                    }
-            }
-
-        return nodalForces;
-    }
+    static std::vector<fea::Force> getFeaNodalForces(const std::shared_ptr<SimulationJob> &job);
     static SimulationResults getResults(const fea::Summary &feaSummary,
-                                        const std::shared_ptr<SimulationJob> &simulationJob)
-    {
-        SimulationResults results;
-        results.converged = feaSummary.converged;
-        if (!results.converged) {
-            return results;
-        }
+                                        const std::shared_ptr<SimulationJob> &simulationJob);
 
-        results.executionTime = feaSummary.total_time_in_ms * 1000;
-        // displacements
-        results.displacements.resize(feaSummary.num_nodes);
-        results.rotationalDisplacementQuaternion.resize(feaSummary.num_nodes);
-        for (int vi = 0; vi < feaSummary.num_nodes; vi++) {
-            results.displacements[vi] = Vector6d(feaSummary.nodal_displacements[vi]);
-
-            const Vector6d &nodalDisplacement = results.displacements[vi];
-            Eigen::Quaternion<double> q_nx;
-            q_nx = Eigen::AngleAxis<double>(nodalDisplacement[3], Eigen::Vector3d(1, 0, 0));
-            Eigen::Quaternion<double> q_ny;
-            q_ny = Eigen::AngleAxis<double>(nodalDisplacement[4], Eigen::Vector3d(0, 1, 0));
-            Eigen::Quaternion<double> q_nz;
-            q_nz = Eigen::AngleAxis<double>(nodalDisplacement[5], Eigen::Vector3d(0, 0, 1));
-            results.rotationalDisplacementQuaternion[vi] = q_nx * q_ny * q_nz;
-            //            results.rotationalDisplacementQuaternion[vi] = q_nz * q_ny * q_nx;
-            //            results.rotationalDisplacementQuaternion[vi] = q_nz * q_nx * q_ny;
-        }
-        results.setLabelPrefix("Linear");
-
-        //  // Convert forces
-        //  // Convert to vector of eigen matrices of the form force component-> per
-        //  // Edge
-        //  const int numDof = 6;
-        //  const size_t numberOfEdges = feaSummary.element_forces.size();
-        //  edgeForces =
-        //      std::vector<Eigen::VectorXd>(numDof, Eigen::VectorXd(2 *
-        //      numberOfEdges));
-        //  for (gsl::index edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
-        //    for (gsl::index forceComponentIndex = 0; forceComponentIndex < numDof;
-        //         forceComponentIndex++) {
-        //      (edgeForces[forceComponentIndex])(2 * edgeIndex) =
-        //          feaSummary.element_forces[edgeIndex][forceComponentIndex];
-        //      (edgeForces[forceComponentIndex])(2 * edgeIndex + 1) =
-        //          feaSummary.element_forces[edgeIndex][numDof +
-        //          forceComponentIndex];
-        //    }
-        //  }
-
-        for (int ei = 0; ei < feaSummary.num_elems; ei++) {
-            const std::vector<double> &elementForces = feaSummary.element_forces[ei];
-            const Element &element = simulationJob->pMesh->elements[ei];
-            //Axial
-            const double elementPotentialEnergy_axial_v0 = std::pow(elementForces[0], 2)
-                                                           * element.initialLength
-                                                           / (4 * element.material.youngsModulus
-                                                              * element.A);
-            const double elementPotentialEnergy_axial_v1 = std::pow(elementForces[6], 2)
-                                                           * element.initialLength
-                                                           / (4 * element.material.youngsModulus
-                                                              * element.A);
-            const double elementPotentialEnergy_axial = elementPotentialEnergy_axial_v0
-                                                        + elementPotentialEnergy_axial_v1;
-            //Shear
-            const double elementPotentialEnergy_shearY_v0 = std::pow(elementForces[1], 2)
-                                                            * element.initialLength
-                                                            / (4 * element.A * element.G);
-            const double elementPotentialEnergy_shearZ_v0 = std::pow(elementForces[2], 2)
-                                                            * element.initialLength
-                                                            / (4 * element.A * element.G);
-            const double elementPotentialEnergy_shearY_v1 = std::pow(elementForces[7], 2)
-                                                            * element.initialLength
-                                                            / (4 * element.A * element.G);
-            const double elementPotentialEnergy_shearZ_v1 = std::pow(elementForces[8], 2)
-                                                            * element.initialLength
-                                                            / (4 * element.A * element.G);
-            const double elementPotentialEnergy_shear = elementPotentialEnergy_shearY_v0
-                                                        + elementPotentialEnergy_shearZ_v0
-                                                        + elementPotentialEnergy_shearY_v1
-                                                        + elementPotentialEnergy_shearZ_v1;
-            //Bending
-            const double elementPotentialEnergy_bendingY_v0 = std::pow(elementForces[4], 2)
-                                                              * element.initialLength
-                                                              / (4 * element.material.youngsModulus
-                                                                 * element.inertia.I2);
-            const double elementPotentialEnergy_bendingZ_v0 = std::pow(elementForces[5], 2)
-                                                              * element.initialLength
-                                                              / (4 * element.material.youngsModulus
-                                                                 * element.inertia.I3);
-            const double elementPotentialEnergy_bendingY_v1 = std::pow(elementForces[10], 2)
-                                                              * element.initialLength
-                                                              / (4 * element.material.youngsModulus
-                                                                 * element.inertia.I2);
-            const double elementPotentialEnergy_bendingZ_v1 = std::pow(elementForces[11], 2)
-                                                              * element.initialLength
-                                                              / (4 * element.material.youngsModulus
-                                                                 * element.inertia.I3);
-            const double elementPotentialEnergy_bending = elementPotentialEnergy_bendingY_v0
-                                                          + elementPotentialEnergy_bendingZ_v0
-                                                          + elementPotentialEnergy_bendingY_v1
-                                                          + elementPotentialEnergy_bendingZ_v1;
-            //Torsion
-            const double elementPotentialEnergy_torsion_v0 = std::pow(elementForces[3], 2)
-                                                             * element.initialLength
-                                                             / (4 * element.inertia.J * element.G);
-            const double elementPotentialEnergy_torsion_v1 = std::pow(elementForces[9], 2)
-                                                             * element.initialLength
-                                                             / (4 * element.inertia.J * element.G);
-            const double elementPotentialEnergy_torsion = elementPotentialEnergy_torsion_v0
-                                                          + elementPotentialEnergy_torsion_v1;
-            const double elementInternalPotentialEnergy = elementPotentialEnergy_axial
-                                                          + elementPotentialEnergy_shear
-                                                          + elementPotentialEnergy_bending
-                                                          + elementPotentialEnergy_torsion;
-            results.internalPotentialEnergy += elementInternalPotentialEnergy;
-        }
-        results.pJob = simulationJob;
-
-        return results;
-    }
-
-  SimulationResults
-  executeSimulation(const std::shared_ptr<SimulationJob> &simulationJob) {
-      assert(simulationJob->pMesh->VN() != 0);
-      fea::Job job(getFeaNodes(simulationJob), getFeaElements(simulationJob));
-      //  printInfo(job);
-
-      // create the default options
-      fea::Options opts;
-      opts.save_elemental_forces = false;
-      opts.save_nodal_displacements = false;
-      opts.save_nodal_forces = false;
-      opts.save_report = false;
-      opts.save_tie_forces = false;
-      //  if (!elementalForcesOutputFilepath.empty()) {
-      //    opts.save_elemental_forces = true;
-      //    opts.elemental_forces_filename = elementalForcesOutputFilepath;
-      //  }
-      //  if (!nodalDisplacementsOutputFilepath.empty()) {
-      //    opts.save_nodal_displacements = true;
-      //    opts.nodal_displacements_filename = nodalDisplacementsOutputFilepath;
-      //  }
-
-      // have the program output status updates
-      opts.verbose = false;
-
-      // form an empty vector of ties
-      std::vector<fea::Tie> ties;
-
-      // also create an empty list of equations
-      std::vector<fea::Equation> equations;
-      auto fixedVertices = getFeaFixedNodes(simulationJob);
-      auto nodalForces = getFeaNodalForces(simulationJob);
-      fea::Summary feaResults =
-              fea::solve(job, fixedVertices, nodalForces, ties, equations, opts);
-
-      SimulationResults results = getResults(feaResults, simulationJob);
-      results.pJob = simulationJob;
-      return results;
-  }
-  //  SimulationResults getResults() const;
-  //  void setResultsNodalDisplacementCSVFilepath(const std::string
-  //  &outputPath); void setResultsElementalForcesCSVFilepath(const std::string
-  //  &outputPath);
+    SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &simulationJob);
+    void setStructure(const std::shared_ptr<SimulationMesh> &pMesh);
 
 private:
-  //  std::string nodalDisplacementsOutputFilepath{"nodal_displacement.csv"};
-  //  std::string elementalForcesOutputFilepath{"elemental_forces.csv"};
-  //  SimulationResults results;
-
-  static void printInfo(const fea::Job &job) {
-      std::cout << "Details regarding the fea::Job:" << std::endl;
-      std::cout << "Nodes:" << std::endl;
-      for (fea::Node n : job.nodes) {
-              std::cout << n << std::endl;
-          }
-      std::cout << "Elements:" << std::endl;
-      for (Eigen::Vector2i e : job.elems) {
-              std::cout << e << std::endl;
-          }
-  }
+    fea::ThreedBeamFEA simulator;
+    static void printInfo(const fea::BeamStructure &job);
 };
 
 #endif // LINEARSIMULATIONMODEL_HPP

mesh.hpp

@@ -1,7 +1,7 @@
 #ifndef MESH_HPP
 #define MESH_HPP
 
-#include <Eigen/Core>
+#include <filesystem>
 #include <string>
 
 class Mesh
@@ -11,8 +11,8 @@ protected:
 
 public:
   virtual ~Mesh() = default;
-  virtual bool load(const std::string &filePath) { return false; }
-  virtual bool save(const std::string &filePath) { return false; }
+  virtual bool load(const std::filesystem::path &meshFilePath) { return false; }
+  virtual bool save(const std::filesystem::path &meshFilePath) { return false; }
 
   std::string getLabel() const;
   void setLabel(const std::string &newLabel);

polymesh.hpp

@@ -5,6 +5,7 @@
 #include <filesystem>
 #include <wrap/io_trimesh/export_obj.h>
 #include <wrap/io_trimesh/export_ply.h>
+#include <wrap/io_trimesh/import.h>
 
 #ifdef POLYSCOPE_DEFINED
 #include <polyscope/surface_mesh.h>
@@ -64,14 +65,25 @@ class VCGPolyMesh
       public Mesh
 {
 public:
-    //  virtual bool load(const std::string &filename) override {
-    //    int mask;
-    //    vcg::tri::io::Importer<VCGPolyMesh>::LoadMask(filename.c_str(), mask);
-    //    int error = vcg::tri::io::Importer<VCGPolyMesh>::Open(
-    //        *this, filename.c_str(), mask);
-    //    if (error != 0) {
-    //      return false;
-    //    }
+    virtual bool load(const std::filesystem::__cxx11::path &meshFilePath) override
+    {
+        int mask;
+        vcg::tri::io::Importer<VCGPolyMesh>::LoadMask(meshFilePath.c_str(), mask);
+        int error = vcg::tri::io::Importer<VCGPolyMesh>::Open(*this, meshFilePath.c_str(), mask);
+        if (error != 0) {
+            std::cerr << "Could not load polygonal mesh:" << meshFilePath << std::endl;
+            return false;
+        }
+
+        vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
+        vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
+        vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
+        vcg::tri::UpdateNormal<VCGPolyMesh>::PerVertexNormalized(*this);
+        vcg::tri::Clean<VCGPolyMesh>::RemoveUnreferencedVertex(*this);
+        //finally remove valence 1 vertices on the border
+        //        vcg::PolygonalAlgorithm<PolyMeshType>::RemoveValence2Vertices(dual);
+        return true;
+    }
     //    //    vcg::tri::io::ImporterOBJ<VCGPolyMesh>::Open();
     //    //    unsigned int mask = 0;
     //    //    mask |= nanoply::NanoPlyWrapper<VCGPolyMesh>::IO_VERTCOORD;
@@ -84,11 +96,6 @@ public:
     //    //            0) {
     //    //      std::cout << "Could not load tri mesh" << std::endl;
     //    //    }
-    //    vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
-    //    vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
-    //    //    vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
-    //    vcg::tri::UpdateNormal<VCGPolyMesh>::PerVertexNormalized(*this);
-    //    return true;
     //  }
     Eigen::MatrixX3d getVertices() const
     {
@@ -113,6 +120,34 @@ public:
       return faces;
   }
 
+  //  bool load(const std::filesystem::__cxx11::path &meshFilePath)
+  //  {
+  //      const std::string extension = ".ply";
+  //      std::filesystem::path filePath = meshFilePath;
+  //      assert(std::filesystem::path(filePath).extension().string() == extension);
+  //      unsigned int mask = 0;
+  //      mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
+  //      mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
+  //      mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
+  //      mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
+  //      if (vcg::tri::io::Importer<VCGPolyMesh>::Open(*this, filePath.c_str()) != 0) {
+  //          return false;
+  //      }
+  //      label = meshFilePath.filename();
+  //      return true;
+  //  }
+
+  bool save(const std::filesystem::__cxx11::path &meshFilePath = std::filesystem::path())
+  {
+      if (meshFilePath.extension() == ".obj") {
+          return saveOBJ(meshFilePath);
+      } else if (meshFilePath.extension() == ".ply") {
+          return savePLY(meshFilePath);
+      }
+
+      return false;
+  }
+
   bool saveOBJ(const std::filesystem::path &objFilePath = std::filesystem::path())
   {
       const std::string extension = ".obj";
@@ -128,7 +163,7 @@ public:
       mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
       mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
       mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
-      if (vcg::tri::io::ExporterOBJ<VCGPolyMesh>::Save(*this, filePath.c_str(), mask) != 0) {
+      if (vcg::tri::io::ExporterOBJ<VCGPolyMesh>::Save(*this, filePath.string().c_str(), mask) != 0) {
           return false;
       }
       return true;
@@ -149,7 +184,7 @@ public:
       mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
       mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
       mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
-      if (vcg::tri::io::ExporterPLY<VCGPolyMesh>::Save(*this, filePath.c_str(), mask, false) != 0) {
+      if (vcg::tri::io::ExporterPLY<VCGPolyMesh>::Save(*this, filePath.string().c_str(), mask, false) != 0) {
           return false;
       }
       return true;

reducedmodelevaluator.cpp

@@ -0,0 +1,317 @@
+#include "reducedmodelevaluator.hpp"
+#include "hexagonremesher.hpp"
+#include "trianglepatterngeometry.hpp"
+#include <filesystem>
+
+using FullPatternVertexIndex = VertexIndex;
+using ReducedPatternVertexIndex = VertexIndex;
+
+ReducedModelEvaluator::ReducedModelEvaluator()
+{
+
+}
+
+std::vector<double> ReducedModelEvaluator::evaluateReducedModel(
+		ReducedModelOptimization::Results &optimizationResults)
+{
+	//    std::shared_ptr<VCGPolyMesh> pTileIntoSurface = std::make_shared<VCGPolyMesh>();
+	//    std::filesystem::path tileIntoSurfaceFilePath{
+	//        "/home/iason/Coding/Projects/Approximating shapes with flat "
+	//        "patterns/ReducedModelOptimization/TestSet/TileIntoSurface/plane_34Polygons.ply"};
+	//    assert(std::filesystem::exists(tileIntoSurfaceFilePath));
+	//    pTileIntoSurface->load(tileIntoSurfaceFilePath);
+
+	//Set required file paths
+	const std::filesystem::path tileInto_triMesh_filename
+		= "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/"
+		"instantMeshes_plane_34.ply";
+	const std::filesystem::path reducedPatternFilePath
+		= "/home/iason/Coding/Projects/Approximating shapes with flat "
+		"patterns/ReducedModelOptimization/TestSet/ReducedPatterns/single_reduced.ply";
+	const std::filesystem::path intermediateResultsDirectoryPath
+		= "/home/iason/Coding/build/ReducedModelOptimization/IntermediateResults";
+	//    const std::filesystem::path drmSettingsFilePath
+	//        = "/home/iason/Coding/Projects/Approximating shapes with flat "
+	//          "patterns/ReducedModelOptimization/DefaultSettings/DRMSettings/"
+	//          "defaultDRMSimulationSettings.json";
+	DRMSimulationModel::Settings drmSimulationSettings;
+	//    drmSimulationSettings.load(drmSettingsFilePath);
+	drmSimulationSettings.linearGuessForceScaleFactor = 0.7;
+	drmSimulationSettings.debugModeStep = 10000;
+	drmSimulationSettings.beVerbose = true;
+	constexpr bool shouldRerunFullPatternSimulation = false;
+	const std::vector<std::string> scenariosTestSetLabels{"22Hex_randomBending0",
+		"22Hex_randomBending1",
+		"22Hex_randomBending2",
+		"22Hex_randomBending3",
+		"22Hex_randomBending4",
+		"22Hex_randomBending5",
+		"22Hex_randomBending6",
+		"22Hex_randomBending7",
+		"22Hex_randomBending8",
+		"22Hex_randomBending9",
+		"22Hex_randomBending10",
+		"22Hex_randomBending11",
+		"22Hex_randomBending12",
+		"22Hex_randomBending13",
+		"22Hex_randomBending14",
+		"22Hex_randomBending15",
+		"22Hex_randomBending16",
+		"22Hex_randomBending17",
+		"22Hex_randomBending18",
+		"22Hex_randomBending19",
+		"22Hex_randomBending20",
+		"22Hex_bending_0.005N",
+		"22Hex_bending_0.01N",
+		"22Hex_bending_0.03N",
+		"22Hex_bending_0.05N",
+		"22Hex_pullOppositeVerts_0.05N",
+		"22Hex_pullOppositeVerts_0.1N",
+		"22Hex_pullOppositeVerts_0.3N",
+		"22Hex_shear_2N",
+		"22Hex_shear_5N",
+		"22Hex_axial_10N",
+		"22Hex_axial_20N",
+		"22Hex_cylinder_0.05N",
+		"22Hex_cylinder_0.1N",
+		"22Hex_s_0.05N",
+		"22Hex_s_0.1N"};
+
+	//Load surface
+	std::shared_ptr<VCGPolyMesh> pTileIntoSurface = [&]() {
+		VCGTriMesh tileInto_triMesh;
+		const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filename);
+		assert(surfaceLoadSuccessfull);
+		return PolygonalRemeshing::remeshWithPolygons(tileInto_triMesh);
+	}();
+	const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResults.baseTriangle.cP(0),
+			(optimizationResults.baseTriangle.cP(1)
+			 + optimizationResults.baseTriangle.cP(
+				 2))
+			/ 2);
+	pTileIntoSurface->moveToCenter();
+	const double scaleFactor = optimizedBaseTriangleHeight
+		/ pTileIntoSurface->getAverageFaceRadius();
+	vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
+	//    tileIntoSurface.registerForDrawing();
+	//    polyscope::show();
+
+	//Tile full pattern into surface
+	std::vector<PatternGeometry> fullPatterns(1);
+	fullPatterns[0].copy(optimizationResults.baseTriangleFullPattern);
+	////        Base triangle pattern might contain dangling vertices.Remove those
+	fullPatterns[0].interfaceNodeIndex = 3;
+	fullPatterns[0].deleteDanglingVertices();
+	std::vector<int> perSurfaceFacePatternIndices(pTileIntoSurface->FN(), 0);
+	std::vector<std::vector<size_t>> perPatternIndexTiledFullPatternEdgeIndices;
+	std::vector<size_t> tileIntoEdgeToTiledFullVi;
+	std::shared_ptr<PatternGeometry> pTiledFullPattern
+		= PatternGeometry::tilePattern(fullPatterns,
+				{},
+				*pTileIntoSurface,
+				perSurfaceFacePatternIndices,
+				tileIntoEdgeToTiledFullVi,
+				perPatternIndexTiledFullPatternEdgeIndices);
+	pTiledFullPattern->setLabel("Tiled_full_patterns");
+	//    pTiledFullPattern->registerForDrawing();
+	//Tile reduced pattern into surface
+	PatternGeometry reducedPattern;
+	reducedPattern.load(reducedPatternFilePath);
+	reducedPattern.deleteDanglingVertices();
+	assert(reducedPattern.interfaceNodeIndex == 1);
+	std::vector<PatternGeometry> reducedPatterns(1);
+	reducedPatterns[0].copy(reducedPattern);
+	const auto reducedPatternBaseTriangle = reducedPattern.computeBaseTriangle();
+	ReducedModelOptimization::Results::applyOptimizationResults_innerHexagon(
+			optimizationResults, reducedPatternBaseTriangle, reducedPatterns[0]);
+
+	std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
+	std::vector<size_t> tileIntoEdgeToTiledReducedVi;
+	std::shared_ptr<PatternGeometry> pTiledReducedPattern
+		= PatternGeometry::tilePattern(reducedPatterns,
+				{0},
+				*pTileIntoSurface,
+				perSurfaceFacePatternIndices,
+				tileIntoEdgeToTiledReducedVi,
+				perPatternIndexTiledReducedPatternEdgeIndices);
+	pTiledReducedPattern->setLabel("Tiled_reduced_patterns");
+	//        pTiledReducedPattern->registerForDrawing();
+
+	std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
+		fullToReducedViMap; //of only the common vertices
+	std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+		reducedToFullViMap; //of only the common vertices
+	for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) {
+		fullToReducedViMap[tileIntoEdgeToTiledFullVi[ei]] = tileIntoEdgeToTiledReducedVi[ei];
+	}
+	constructInverseMap(fullToReducedViMap, reducedToFullViMap);
+
+	std::vector<size_t> tilledFullPatternInterfaceVi;
+	tilledFullPatternInterfaceVi.clear();
+	tilledFullPatternInterfaceVi.resize(fullToReducedViMap.size());
+	size_t viIndex = 0;
+	for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
+		tilledFullPatternInterfaceVi[viIndex++] = fullToReducedPair.first;
+	}
+
+	//Create simulation meshes
+	////Tessellated full pattern simulation mesh
+	std::shared_ptr<SimulationMesh> pTiledFullPattern_simulationMesh;
+	pTiledFullPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledFullPattern);
+	//NOTE: Those should be derived from the optimization results instead of hardcoding them
+	pTiledFullPattern_simulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002});
+	if (optimizationResults.fullPatternYoungsModulus == 0) {
+		std::cerr << "Full pattern's young modulus not found." << std::endl;
+		std::terminate();
+	}
+	pTiledFullPattern_simulationMesh->setBeamMaterial(0.3,
+			optimizationResults.fullPatternYoungsModulus);
+	pTiledFullPattern_simulationMesh->reset();
+
+	////Tessellated reduced pattern simulation mesh
+	std::shared_ptr<SimulationMesh> pTiledReducedPattern_simulationMesh;
+	pTiledReducedPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledReducedPattern);
+	const std::vector<size_t> &tiledPatternElementIndicesForReducedPattern
+		= perPatternIndexTiledReducedPatternEdgeIndices[0];
+	ReducedModelOptimization::Results::applyOptimizationResults_elements(
+			optimizationResults, pTiledReducedPattern_simulationMesh);
+	pTiledReducedPattern_simulationMesh->reset();
+	std::vector<double> distances_drm2reduced(scenariosTestSetLabels.size(), 0);
+	std::vector<double> distances_normalizedDrm2reduced(scenariosTestSetLabels.size(), 0);
+	for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
+		const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
+		const std::filesystem::path scenariosDirectoryPath
+			= "/home/iason/Coding/build/PatternTillingReducedModel/Scenarios/";
+		const std::filesystem::path tiledReducedPatternJobFilePath
+			= std::filesystem::path(scenariosDirectoryPath)
+			.append(jobLabel)
+			.append("SimulationJobs")
+			.append("Reduced")
+			.append(SimulationJob::jsonDefaultFileName);
+		//set jobs
+		std::shared_ptr<SimulationJob> pJob_tiledReducedPattern;
+		pJob_tiledReducedPattern = std::make_shared<SimulationJob>(SimulationJob());
+		pJob_tiledReducedPattern->load(tiledReducedPatternJobFilePath, false);
+		pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
+		std::shared_ptr<SimulationJob> pJob_tiledFullPattern;
+		pJob_tiledFullPattern = std::make_shared<SimulationJob>(SimulationJob());
+		pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
+		pJob_tiledReducedPattern->remap(reducedToFullViMap, *pJob_tiledFullPattern);
+		//        pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern->getLabel());
+		//        pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel());
+		//    polyscope::show();
+		//Save reduced job
+		const std::filesystem::path tesellatedResultsFolderPath
+			= std::filesystem::path(intermediateResultsDirectoryPath).append("TessellatedResults");
+		const std::filesystem::path surfaceFolderPath = std::filesystem::path(
+				tesellatedResultsFolderPath)
+			.append(pTileIntoSurface->getLabel());
+		const std::string scenarioLabel = pJob_tiledFullPattern->getLabel();
+		const std::filesystem::path scenarioDirectoryPath = std::filesystem::path(surfaceFolderPath)
+			.append(scenarioLabel);
+		const std::filesystem::path reducedJobDirectoryPath
+			= std::filesystem::path(scenarioDirectoryPath).append("ReducedJob");
+		std::filesystem::create_directories(reducedJobDirectoryPath);
+		pJob_tiledReducedPattern->save(reducedJobDirectoryPath);
+		//Run scenario
+
+		////Full
+
+		const std::string patternLabel = [&]() {
+			const std::string patternLabel = optimizationResults.baseTriangleFullPattern.getLabel();
+			const int numberOfOccurences = std::count_if(patternLabel.begin(),
+					patternLabel.end(),
+					[](char c) { return c == '#'; });
+			if (numberOfOccurences == 0) {
+				return std::to_string(optimizationResults.baseTriangleFullPattern.EN()) + "#"
+					+ optimizationResults.baseTriangleFullPattern.getLabel();
+			} else if (numberOfOccurences == 1) {
+				return optimizationResults.baseTriangleFullPattern.getLabel();
+			}
+		}();
+		const auto fullResultsFolderPath
+			= std::filesystem::path(scenarioDirectoryPath).append(patternLabel).append("Results");
+		if (shouldRerunFullPatternSimulation && std::filesystem::exists(fullResultsFolderPath)) {
+			std::filesystem::remove_all(fullResultsFolderPath);
+		}
+
+		const std::filesystem::path fullPatternJobFolderPath = std::filesystem::path(
+				scenarioDirectoryPath)
+			.append(patternLabel)
+			.append("SimulationJob");
+		SimulationResults simulationResults_tiledFullPattern_drm;
+		if (std::filesystem::exists(fullResultsFolderPath)) {
+			//Load full pattern results
+			assert(std::filesystem::exists(fullPatternJobFolderPath));
+			simulationResults_tiledFullPattern_drm.load(fullResultsFolderPath,
+					fullPatternJobFolderPath);
+			simulationResults_tiledFullPattern_drm.converged = true;
+		} else {
+			//Full
+			std::cout << "Executing:" << jobLabel << std::endl;
+			DRMSimulationModel drmSimulationModel;
+			simulationResults_tiledFullPattern_drm
+				= drmSimulationModel.executeSimulation(pJob_tiledFullPattern, drmSimulationSettings);
+			simulationResults_tiledFullPattern_drm.setLabelPrefix("DRM");
+		}
+		std::filesystem::create_directories(fullResultsFolderPath);
+		simulationResults_tiledFullPattern_drm.save(
+				std::filesystem::path(scenarioDirectoryPath).append(patternLabel));
+		if (!simulationResults_tiledFullPattern_drm.converged) {
+			std::cerr << "Full pattern simulation failed." << std::endl;
+		}
+
+		LinearSimulationModel linearSimulationModel;
+		SimulationResults simulationResults_tiledReducedPattern
+			= linearSimulationModel.executeSimulation(pJob_tiledReducedPattern);
+
+		//        simulationResults_tiledReducedPattern.registerForDrawing();
+		//        simulationResults_tiledFullPattern_drm.registerForDrawing();
+		//        polyscope::show();
+		//measure distance
+		const double distance_fullDrmToReduced
+			= simulationResults_tiledFullPattern_drm
+			.computeDistance(simulationResults_tiledReducedPattern, fullToReducedViMap);
+		double distance_fullSumOfAllVerts = 0;
+		for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
+			distance_fullSumOfAllVerts += simulationResults_tiledFullPattern_drm
+				.displacements[fullToReducedPair.first]
+				.getTranslation()
+				.norm();
+		}
+		const double distance_normalizedFullDrmToReduced = distance_fullDrmToReduced
+			/ distance_fullSumOfAllVerts;
+		distances_drm2reduced[jobIndex] = distance_fullDrmToReduced;
+		distances_normalizedDrm2reduced[jobIndex] = distance_normalizedFullDrmToReduced;
+	}
+	//#ifndef POLYSCOPE_DEFINED
+	//    return distances_drm2reduced;
+	//#else
+
+	//report distance
+	//        csvFile csv_results("", flse);
+	csvFile csv_results({}, false);
+	//    csvFile csv_results(std::filesystem::path(dirPath_thisOptimization)
+	//                            .append("results.csv")
+	//                            .string(),
+	//                        false);
+	csv_results /*<< "Job Label"*/
+		<< "drm2Reduced"
+		<< "norm_drm2Reduced";
+	csv_results << endrow;
+	double sumOfNormalizedFull2Reduced = 0;
+	for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
+		const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
+		const double &distance_fullDrmToReduced = distances_drm2reduced[jobIndex];
+		const double &distance_normalizedFullDrmToReduced = distances_normalizedDrm2reduced[jobIndex];
+		csv_results /*<< jobLabel*/ << distance_fullDrmToReduced
+			<< distance_normalizedFullDrmToReduced;
+		csv_results << endrow;
+		sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
+	}
+
+	std::cout << "Average normalized error per scenario:"
+		<< sumOfNormalizedFull2Reduced / scenariosTestSetLabels.size() << std::endl;
+	return distances_normalizedDrm2reduced;
+	//#endif
+}

reducedmodelevaluator.hpp

@@ -0,0 +1,14 @@
+#ifndef REDUCEDMODELEVALUATOR_HPP
+#define REDUCEDMODELEVALUATOR_HPP
+
+#include "reducedmodeloptimizer_structs.hpp"
+
+class ReducedModelEvaluator
+{
+public:
+    ReducedModelEvaluator();
+    static std::vector<double> evaluateReducedModel(
+        ReducedModelOptimization::Results &optimizationResults);
+};
+
+#endif // REDUCEDMODELEVALUATOR_HPP

reducedmodeloptimizer_structs.hpp

@@ -7,7 +7,7 @@
 #include "unordered_map"
 #include <string>
 
-namespace ReducedPatternOptimization {
+namespace ReducedModelOptimization {
 
 enum BaseSimulationScenario {
     Axial,
@@ -34,9 +34,9 @@ struct Colors
 
 struct xRange
 {
-    std::string label;
-    double min;
-    double max;
+    std::string label{};
+    double min{0};
+    double max{0};
 
     inline bool operator<(const xRange &other) { return label < other.label; }
     std::string toString() const
@@ -60,151 +60,331 @@ struct xRange
     {
         return label == xrange.label && min == xrange.min && max == xrange.max;
     }
+
+    std::tuple<std::string, double, double> toTuple() const
+    {
+        return std::make_tuple(label, min, max);
+    }
+
+    void set(const std::tuple<std::string, double, double> &inputTuple)
+    {
+        if (std::get<1>(inputTuple) > std::get<2>(inputTuple)) {
+            std::cerr
+                << "Invalid xRange tuple. Second argument(min) cant be smaller than the third(max)"
+                << std::endl;
+            std::terminate();
+            //            return;
+        }
+        std::tie(label, min, max) = inputTuple;
+    }
 };
 
-  struct Settings
-  {
-      enum NormalizationStrategy { NonNormalized, Epsilon };
-      std::vector<xRange> xRanges;
-      inline static vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
-      int numberOfFunctionCalls{100000};
-      double solverAccuracy{1e-2};
-      NormalizationStrategy normalizationStrategy{Epsilon};
-      double translationNormalizationParameter{0.003};
-      double rotationNormalizationParameter{3};
-      bool splitGeometryMaterialOptimization{false};
-      struct ObjectiveWeights
-      {
-          double translational{1};
-          double rotational{1};
-      } objectiveWeights;
+enum OptimizationParameterIndex { E, A, I2, I3, J, R, Theta, NumberOfOptimizationVariables };
 
-      struct JsonKeys
-      {
-          inline static std::string filename{"OptimizationSettings.json"};
-          inline static std::string OptimizationVariables{"OptimizationVariables"};
-          inline static std::string NumberOfFunctionCalls{"NumberOfFunctionCalls"};
-          inline static std::string SolverAccuracy{"SolverAccuracy"};
-          inline static std::string TranslationalObjectiveWeight{"TransObjWeight"};
-      };
+inline int getParameterIndex(const std::string &s)
+{
+    if ("E" == s) {
+        return E;
+    } else if ("A" == s) {
+        return A;
+    } else if ("I2" == s) {
+        return I2;
+    } else if ("I3" == s) {
+        return I3;
+    } else if ("J" == s) {
+        return J;
+    } else if ("R" == s || "HexSize" == s) {
+        return R;
+    } else if ("Theta" == s || "HexAngle" == s) {
+        return Theta;
+    } else {
+        std::cerr << "Input is not recognized as a valid optimization variable index:" << s
+                  << std::endl;
+        return -1;
+    }
+}
 
-      void save(const std::filesystem::path &saveToPath)
-      {
-          assert(std::filesystem::is_directory(saveToPath));
+struct Settings
+{
+    inline static std::string defaultFilename{"OptimizationSettings.json"};
+    std::vector<std::vector<OptimizationParameterIndex>> optimizationStrategy = {
+        {E, A, I2, I3, J, R, Theta}};
+    //        {A, I2, I3, J, R, Theta}};
+    std::optional<std::vector<double>>
+        optimizationVariablesGroupsWeights; //TODO:should be removed in the future if not a splitting strategy is used for the optimization
+    enum NormalizationStrategy { NonNormalized, Epsilon };
+    inline static std::vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
+    NormalizationStrategy normalizationStrategy{Epsilon};
+    std::array<xRange, NumberOfOptimizationVariables> variablesRanges{xRange{"E", 0.01, 100},
+                                                                      xRange{"A", 0.01, 100},
+                                                                      xRange{"I2", 0.01, 100},
+                                                                      xRange{"I3", 0.01, 100},
+                                                                      xRange{"J", 0.01, 100},
+                                                                      xRange{"R", 0.05, 0.95},
+                                                                      xRange{"Theta", -30, 30}};
+    int numberOfFunctionCalls{100000};
+    double solverAccuracy{1e-3};
+    double translationEpsilon{3e-4};
+    //    double translationEpsilon{0};
+    double angularDistanceEpsilon{vcg::math::ToRad(0.0)};
+    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;
+        }
+    };
+    std::array<ObjectiveWeights, NumberOfBaseSimulationScenarios> perBaseScenarioObjectiveWeights;
+    std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios>
+    convertObjectiveWeightsToPairs() const
+    {
+        std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios> objectiveWeightsPairs;
+        for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
+             baseScenario++) {
+            objectiveWeightsPairs[baseScenario]
+                = std::make_pair(perBaseScenarioObjectiveWeights[baseScenario].translational,
+                                 perBaseScenarioObjectiveWeights[baseScenario].rotational);
+        }
+        return objectiveWeightsPairs;
+    }
 
-          nlohmann::json json;
-          //write x ranges
-          for (size_t xRangeIndex = 0; xRangeIndex < xRanges.size(); xRangeIndex++) {
-              const auto &xRange = xRanges[xRangeIndex];
-              json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)]
-                  = xRange.toString();
-          }
+    struct JsonKeys
+    {
+        inline static std::string OptimizationStrategy{"OptimizationStrategy"};
+        inline static std::string OptimizationStrategyGroupWeights{
+            "OptimizationStrategyGroupWeights"};
+        inline static std::string OptimizationVariables{"OptimizationVariables"};
+        inline static std::string NumberOfFunctionCalls{"NumberOfFunctionCalls"};
+        inline static std::string SolverAccuracy{"SolverAccuracy"};
+        inline static std::string ObjectiveWeights{"ObjectiveWeight"};
+    };
 
-          json[JsonKeys::NumberOfFunctionCalls] = numberOfFunctionCalls;
-          json[JsonKeys::SolverAccuracy] = solverAccuracy;
-          json[JsonKeys::TranslationalObjectiveWeight] = objectiveWeights.translational;
+    void save(const std::filesystem::path &saveToPath)
+    {
+        assert(std::filesystem::is_directory(saveToPath));
 
-          std::filesystem::path jsonFilePath(
-              std::filesystem::path(saveToPath).append(JsonKeys::filename));
-          std::ofstream jsonFile(jsonFilePath.string());
-          jsonFile << json;
-      }
+        nlohmann::json json;
+        json[GET_VARIABLE_NAME(optimizationStrategy)] = optimizationStrategy;
+        if (optimizationVariablesGroupsWeights.has_value()) {
+            json[GET_VARIABLE_NAME(ptimizationStrategyGroupWeights)]
+                = optimizationVariablesGroupsWeights.value();
+        }
+        //write x ranges
+        const std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
+            xRangesAsTuples = [=]() {
+                std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
+                    xRangesAsTuples;
+                for (int optimizationParameterIndex = E;
+                     optimizationParameterIndex != NumberOfOptimizationVariables;
+                     optimizationParameterIndex++) {
+                    xRangesAsTuples[optimizationParameterIndex]
+                        = variablesRanges[optimizationParameterIndex].toTuple();
+                }
+                return xRangesAsTuples;
+            }();
+        json[JsonKeys::OptimizationVariables] = xRangesAsTuples;
+        //        for (size_t xRangeIndex = 0; xRangeIndex < variablesRanges.size(); xRangeIndex++) {
+        //            const auto &xRange = variablesRanges[xRangeIndex];
+        //            json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)]
+        //                = xRange.toString();
+        //        }
 
-      bool load(const std::filesystem::path &loadFromPath)
-      {
-          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)) {
-              return false;
-          }
-          std::ifstream ifs(jsonFilepath);
-          ifs >> json;
+        json[GET_VARIABLE_NAME(numberOfFunctionCalls)] = numberOfFunctionCalls;
+        json[GET_VARIABLE_NAME(solverAccuracy)] = solverAccuracy;
+        //Objective weights
+        std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios> objectiveWeightsPairs;
+        std::transform(perBaseScenarioObjectiveWeights.begin(),
+                       perBaseScenarioObjectiveWeights.end(),
+                       objectiveWeightsPairs.begin(),
+                       [](const ObjectiveWeights &objectiveWeights) {
+                           return std::make_pair(objectiveWeights.translational,
+                                                 objectiveWeights.rotational);
+                       });
+        json[JsonKeys::ObjectiveWeights] = objectiveWeightsPairs;
+        json[GET_VARIABLE_NAME(translationEpsilon)] = translationEpsilon;
+        json[GET_VARIABLE_NAME(angularDistanceEpsilon)] = vcg::math::ToDeg(
+            angularDistanceEpsilon);
 
-          //read x ranges
-          size_t xRangeIndex = 0;
-          while (true) {
-              const std::string jsonXRangeKey = JsonKeys::OptimizationVariables + "_"
-                                                + std::to_string(xRangeIndex);
-              if (!json.contains(jsonXRangeKey)) {
-                  break;
-              }
-              xRange x;
-              x.fromString(json.at(jsonXRangeKey));
-              xRanges.push_back(x);
-              xRangeIndex++;
-          }
+        std::filesystem::path jsonFilePath(
+            std::filesystem::path(saveToPath).append(defaultFilename));
+        std::ofstream jsonFile(jsonFilePath.string());
+        jsonFile << json;
+        jsonFile.close();
+    }
 
-          numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls);
-          solverAccuracy = json.at(JsonKeys::SolverAccuracy);
-          objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight);
-          objectiveWeights.rotational = 2 - objectiveWeights.translational;
-          return true;
-      }
+    bool load(const std::filesystem::path &loadFromPath)
+    {
+        assert(std::filesystem::is_directory(loadFromPath));
+        //Load optimal X
+        nlohmann::json json;
+        std::filesystem::path jsonFilepath(
+            std::filesystem::path(loadFromPath).append(defaultFilename));
+        if (!std::filesystem::exists(jsonFilepath)) {
+            std::cerr << "Optimization settings could not be loaded because input filepath does "
+                         "not exist:"
+                      << jsonFilepath << std::endl;
+            return false;
+        }
+        std::ifstream ifs(jsonFilepath);
+        ifs >> json;
 
-      std::string toString() const
-      {
-          std::string settingsString;
-          if (!xRanges.empty()) {
-              std::string xRangesString;
-              for (const xRange &range : xRanges) {
-                  xRangesString += range.toString() + " ";
-              }
-              settingsString += xRangesString;
-          }
-          settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls)
-                            + " Accuracy=" + std::to_string(solverAccuracy)
-                            + " TransWeight=" + std::to_string(objectiveWeights.translational)
-                            + " RotWeight=" + std::to_string(objectiveWeights.rotational);
+        if (json.contains(GET_VARIABLE_NAME(optimizationStrategy))) {
+            optimizationStrategy
+                = std::vector<std::vector<ReducedModelOptimization::OptimizationParameterIndex>>(
+                    (json.at(GET_VARIABLE_NAME(optimizationStrategy))));
+        }
+        if (json.contains(GET_VARIABLE_NAME(optimizationStrategyGroupWeights))) {
+            optimizationVariablesGroupsWeights = std::vector<double>(
+                json[GET_VARIABLE_NAME(optimizationStrategyGroupWeights)]);
+        }
+        //read x ranges
+        if (json.contains(JsonKeys::OptimizationVariables)) {
+            const std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
+                xRangesAsTuples = json.at(JsonKeys::OptimizationVariables);
+            for (const auto &rangeTuple : xRangesAsTuples) {
+                variablesRanges[getParameterIndex(std::get<0>(rangeTuple))].set(rangeTuple);
+            }
+        } else { //NOTE:legacy compatibility
+            size_t xRangeIndex = 0;
+            while (true) {
+                const std::string jsonXRangeKey = JsonKeys::OptimizationVariables + "_"
+                                                  + std::to_string(xRangeIndex++);
+                if (!json.contains(jsonXRangeKey)) {
+                    break;
+                }
+                xRange x;
+                x.fromString(json.at(jsonXRangeKey));
+                variablesRanges[getParameterIndex(x.label)] = x;
+            }
+        }
 
-          return settingsString;
-      }
+        if (json.contains(GET_VARIABLE_NAME(numberOfFunctionCalls))) {
+            numberOfFunctionCalls = json.at(GET_VARIABLE_NAME(numberOfFunctionCalls));
+        }
+        if (json.contains(GET_VARIABLE_NAME(solverAccuracy))) {
+            solverAccuracy = json.at(GET_VARIABLE_NAME(solverAccuracy));
+        }
+        //Objective weights
+        if (json.contains(JsonKeys::ObjectiveWeights)) {
+            std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios>
+                objectiveWeightsPairs = json.at(JsonKeys::ObjectiveWeights);
+            std::transform(objectiveWeightsPairs.begin(),
+                           objectiveWeightsPairs.end(),
+                           perBaseScenarioObjectiveWeights.begin(),
+                           [](const std::pair<double, double> &objectiveWeightsPair) {
+                               return ObjectiveWeights{objectiveWeightsPair.first,
+                                                       objectiveWeightsPair.second};
+                           });
+        }
+        if (json.contains(GET_VARIABLE_NAME(translationalNormalizationEpsilon))) {
+            translationEpsilon
+                = json[GET_VARIABLE_NAME(translationalNormalizationEpsilon)];
+        }
 
-      void writeHeaderTo(csvFile &os) const
-      {
-          if (!xRanges.empty()) {
-              for (const xRange &range : xRanges) {
-                  os << range.label + " max";
-                  os << range.label + " min";
-              }
-          }
-          os << "Function Calls";
-          os << "Solution Accuracy";
-          os << "Normalization strategy";
-          os << "Trans weight";
-          os << "Rot weight";
-          os << "Splitted geo from mat opt";
-          //        os << std::endl;
-      }
+        if (json.contains(GET_VARIABLE_NAME(angularDistanceEpsilon))) {
+            angularDistanceEpsilon = vcg::math::ToRad(
+                static_cast<double>(json[GET_VARIABLE_NAME(angularDistanceEpsilon)]));
+        }
 
-      void writeSettingsTo(csvFile &os) const
-      {
-          if (!xRanges.empty()) {
-              for (const xRange &range : xRanges) {
-                  os << range.max;
-                  os << range.min;
-              }
-          }
-          os << numberOfFunctionCalls;
-          os << solverAccuracy;
-          os << normalizationStrategyStrings[normalizationStrategy] + "_"
-                    + std::to_string(translationNormalizationParameter);
-          os << objectiveWeights.translational;
-          os << objectiveWeights.rotational;
-          os << (splitGeometryMaterialOptimization == true ? "true" : "false");
-      }
-  };
+        //        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;
+    }
+
+    void writeHeaderTo(csvFile &os) const
+    {
+        if (!variablesRanges.empty()) {
+            for (const xRange &range : variablesRanges) {
+                os << range.label + " max";
+                os << range.label + " min";
+            }
+        }
+        os << "Function Calls";
+        os << "Solution Accuracy";
+        os << "Normalization trans epsilon";
+        os << "Normalization rot epsilon(deg)";
+        os << JsonKeys::ObjectiveWeights;
+        os << "Optimization parameters";
+        //        os << std::endl;
+    }
+
+    void writeSettingsTo(csvFile &os) const
+    {
+        if (!variablesRanges.empty()) {
+            for (const xRange &range : variablesRanges) {
+                os << range.max;
+                os << range.min;
+            }
+        }
+        os << numberOfFunctionCalls;
+        os << solverAccuracy;
+        os << std::to_string(translationEpsilon);
+        os << std::to_string(vcg::math::ToDeg(angularDistanceEpsilon));
+        std::string objectiveWeightsString;
+        objectiveWeightsString += "{";
+        for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
+             baseScenario++) {
+            objectiveWeightsString
+                += "{" + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].translational)
+                   + "," + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].rotational)
+                   + "}";
+        }
+        objectiveWeightsString += "}";
+        os << objectiveWeightsString;
+
+        //export optimization parameters
+        std::vector<std::vector<int>> vv;
+        for (const std::vector<OptimizationParameterIndex> &v : optimizationStrategy) {
+            std::vector<int> vi;
+            vi.reserve(v.size());
+            for (const OptimizationParameterIndex &parameter : v) {
+                vi.emplace_back(parameter);
+            }
+            vv.push_back(vi);
+        }
+        os << Utilities::toString(vv);
+    }
+};
 
   inline bool operator==(const Settings &settings1, const Settings &settings2) noexcept
   {
+      const bool haveTheSameObjectiveWeights
+          = std::mismatch(settings1.perBaseScenarioObjectiveWeights.begin(),
+                          settings1.perBaseScenarioObjectiveWeights.end(),
+                          settings2.perBaseScenarioObjectiveWeights.begin())
+                .first
+            != settings1.perBaseScenarioObjectiveWeights.end();
       return settings1.numberOfFunctionCalls == settings2.numberOfFunctionCalls
-             && settings1.xRanges == settings2.xRanges
-             && settings1.solverAccuracy == settings2.solverAccuracy
-             && settings1.objectiveWeights.translational == settings2.objectiveWeights.translational
-             && settings1.translationNormalizationParameter
-                    == settings2.translationNormalizationParameter;
+             && settings1.variablesRanges == settings2.variablesRanges
+             && settings1.solverAccuracy == settings2.solverAccuracy && haveTheSameObjectiveWeights
+             && settings1.translationEpsilon
+                    == settings2.translationEpsilon;
   }
 
   inline void updateMeshWithOptimalVariables(const std::vector<double> &x, SimulationMesh &m)
@@ -222,9 +402,9 @@ struct xRange
           Element &e = m.elements[ei];
           e.setDimensions(CrossSectionType(beamWidth, beamHeight));
           e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
-          e.inertia.J = J;
-          e.inertia.I2 = I2;
-          e.inertia.I3 = I3;
+          e.dimensions.inertia.J = J;
+          e.dimensions.inertia.I2 = I2;
+          e.dimensions.inertia.I3 = I3;
       }
 
       CoordType center_barycentric(1, 0, 0);
@@ -271,6 +451,7 @@ struct xRange
       std::vector<std::pair<std::string, double>> optimalXNameValuePairs;
       std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
       std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
+      double fullPatternYoungsModulus{0};
 
       PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern
       vcg::Triangle3<double> baseTriangle;
@@ -286,6 +467,7 @@ struct xRange
           std::vector<double> perSimulationScenario_translational;
           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;
@@ -294,8 +476,8 @@ struct xRange
       struct CSVExportingSettings
       {
           bool exportPE{false};
-          bool exportIterationOfMinima{true};
-          bool exportRawObjectiveValue{false};
+          bool exportIterationOfMinima{false};
+          bool exportRawObjectiveValue{true};
           CSVExportingSettings() {}
       };
 
@@ -309,9 +491,10 @@ struct xRange
           inline static std::string FullPatternLabel{"FullPatternLabel"};
           inline static std::string Settings{"OptimizationSettings"};
           inline static std::string FullPatternPotentialEnergies{"PE"};
+          inline static std::string fullPatternYoungsModulus{"youngsModulus"};
       };
 
-      void save(const string &saveToPath, const bool shouldExportDebugFiles = false)
+      void save(const std::string &saveToPath, const bool shouldExportDebugFiles = false)
       {
           //clear directory
           if (std::filesystem::exists(saveToPath)) {
@@ -362,6 +545,7 @@ struct xRange
                   = 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));
@@ -383,7 +567,8 @@ struct xRange
                       = fullPatternSimulationJobs[simulationScenarioIndex];
                   std::filesystem::path simulationJobFolderPath(
                       std::filesystem::path(simulationJobsPath)
-                          .append(pFullPatternSimulationJob->label));
+                          .append(std::to_string(simulationScenarioIndex) + "_"
+                                  + pFullPatternSimulationJob->label));
                   std::filesystem::create_directories(simulationJobFolderPath);
                   const auto fullPatternDirectoryPath
                       = std::filesystem::path(simulationJobFolderPath).append("Full");
@@ -413,12 +598,13 @@ struct xRange
 #endif
       }
 
-      bool load(const string &loadFromPath, const bool &shouldLoadDebugFiles = false)
+      bool load(const std::string &loadFromPath, const bool &shouldLoadDebugFiles = false)
       {
           assert(std::filesystem::is_directory(loadFromPath));
           std::filesystem::path jsonFilepath(
               std::filesystem::path(loadFromPath).append(JsonKeys::filename));
           if (!std::filesystem::exists(jsonFilepath)) {
+              std::cerr << "Input path does not exist:" << loadFromPath << std::endl;
               return false;
           }
           //Load optimal X
@@ -465,59 +651,81 @@ struct xRange
           baseTriangle.P2(0) = CoordType(baseTriangleVertices[6],
                                          baseTriangleVertices[7],
                                          baseTriangleVertices[8]);
+          if (json_optimizationResults.contains(JsonKeys::fullPatternYoungsModulus)) {
+              fullPatternYoungsModulus = json_optimizationResults.at(
+                  JsonKeys::fullPatternYoungsModulus);
+          } else {
+              fullPatternYoungsModulus = 1 * 1e9;
+          }
           if (shouldLoadDebugFiles) {
               const std::filesystem::path simulationJobsFolderPath(
                   std::filesystem::path(loadFromPath).append("SimulationJobs"));
-              for (const auto &directoryEntry :
-                   filesystem::directory_iterator(simulationJobsFolderPath)) {
-                  const auto simulationScenarioPath = directoryEntry.path();
+
+              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) {
                   if (!std::filesystem::is_directory(simulationScenarioPath)) {
                       continue;
                   }
                   // Load full pattern files
-                  for (const auto &fileEntry : filesystem::directory_iterator(
+                  for (const auto &fileEntry : std::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 reduced pattern files and apply optimal parameters
-                  for (const auto &fileEntry : filesystem::directory_iterator(
-                           std::filesystem::path(simulationScenarioPath).append("Reduced"))) {
-                      const auto filepath = fileEntry.path();
-                      if (filepath.extension() == ".json") {
-                          SimulationJob job;
-                          job.load(filepath.string());
-                          //                          applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh);
-                          applyOptimizationResults_elements(*this, job.pMesh);
-                          reducedPatternSimulationJobs.push_back(
-                              std::make_shared<SimulationJob>(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));
+
+                  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));
               }
           }
-          settings.load(loadFromPath);
+          settings.load(std::filesystem::path(loadFromPath).append(Settings::defaultFilename));
 
           return true;
       }
 
       template<typename MeshType>
       static void applyOptimizationResults_innerHexagon(
-          const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
+          const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
           const vcg::Triangle3<double> &patternBaseTriangle,
           MeshType &reducedPattern)
       {
           std::unordered_map<std::string, double>
               optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(),
                                 reducedPattern_optimizationResults.optimalXNameValuePairs.end());
-          assert(optimalXVariables.contains("HexSize") && optimalXVariables.contains("HexAngle"));
-          applyOptimizationResults_innerHexagon(optimalXVariables.at("HexSize"),
-                                                optimalXVariables.at("HexAngle"),
+          assert(
+              (optimalXVariables.contains("R") && optimalXVariables.contains("Theta"))
+              || (optimalXVariables.contains("HexSize") && optimalXVariables.contains("HexAngle")));
+          if (optimalXVariables.contains("HexSize")) {
+              applyOptimizationResults_innerHexagon(optimalXVariables.at("HexSize"),
+                                                    optimalXVariables.at("HexAngle"),
+                                                    patternBaseTriangle,
+                                                    reducedPattern);
+              return;
+          }
+          applyOptimizationResults_innerHexagon(optimalXVariables.at("R"),
+                                                optimalXVariables.at("Theta"),
                                                 patternBaseTriangle,
                                                 reducedPattern);
       }
@@ -560,8 +768,8 @@ struct xRange
       }
 
       static void applyOptimizationResults_elements(
-          const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
-          const shared_ptr<SimulationMesh> &pTiledReducedPattern_simulationMesh)
+          const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
+          const std::shared_ptr<SimulationMesh> &pReducedPattern_simulationMesh)
       {
           assert(CrossSectionType::name == RectangularBeamDimensions::name);
           // Set reduced parameters fron the optimization results
@@ -575,8 +783,8 @@ struct xRange
               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);
               }
           }
@@ -586,8 +794,8 @@ struct xRange
           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);
               }
           }
@@ -595,30 +803,30 @@ struct xRange
           const std::string JLabel = "J";
           if (optimalXVariables.contains(JLabel)) {
               const double J = optimalXVariables.at(JLabel);
-              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
-                  e.inertia.J = J;
+              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+                  Element &e = pReducedPattern_simulationMesh->elements[ei];
+                  e.dimensions.inertia.J = J;
               }
           }
 
           const std::string I2Label = "I2";
           if (optimalXVariables.contains(I2Label)) {
               const double I2 = optimalXVariables.at(I2Label);
-              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
-                  e.inertia.I2 = I2;
+              for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
+                  Element &e = pReducedPattern_simulationMesh->elements[ei];
+                  e.dimensions.inertia.I2 = I2;
               }
           }
 
           const std::string I3Label = "I3";
           if (optimalXVariables.contains(I3Label)) {
               const double I3 = optimalXVariables.at(I3Label);
-              for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
-                  Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
-                  e.inertia.I3 = I3;
+              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
@@ -715,7 +923,7 @@ struct xRange
       }
     }
 
-    void writeHeaderTo(csvFile &os)
+    void writeHeaderTo(csvFile &os) const
     {
         if (exportSettings.exportRawObjectiveValue) {
             os << "Total raw Obj value";
@@ -751,13 +959,20 @@ struct xRange
         //        os << "notes";
     }
 
+    void writeHeaderTo(std::vector<csvFile *> &vectorOfPointersToOutputStreams) const
+    {
+        for (csvFile *outputStream : vectorOfPointersToOutputStreams) {
+            writeHeaderTo(*outputStream);
+        }
+    }
+
     void writeResultsTo(csvFile &os) const
     {
         if (exportSettings.exportRawObjectiveValue) {
             os << objectiveValue.totalRaw;
         }
         os << objectiveValue.total;
-        if (exportSettings.exportIterationOfMinima) {
+        if (exportSettings.exportIterationOfMinima && !objectiveValueHistory_iteration.empty()) {
             os << objectiveValueHistory_iteration.back();
         }
         for (int simulationScenarioIndex = 0;
@@ -787,6 +1002,13 @@ struct xRange
         //        os << notes;
     }
 
+    void writeResultsTo(std::vector<csvFile *> &vectorOfPointersToOutputStreams) const
+    {
+        for (csvFile *&outputStream : vectorOfPointersToOutputStreams) {
+            writeResultsTo(*outputStream);
+        }
+    }
+
     void writeMinimaInfoTo(csvFile &outputCsv)
     {
         outputCsv << "Iteration";

simulation_structs.hpp

@@ -1,5 +1,6 @@
 #ifndef SIMULATIONSTRUCTS_HPP
 #define SIMULATIONSTRUCTS_HPP
+#include <fstream>
 
 namespace Eigen {
 template <class Matrix>
@@ -23,13 +24,13 @@ void read_binary(const std::string &filename, Matrix &matrix) {
   in.read((char *)matrix.data(), rows * cols * sizeof(typename Matrix::Scalar));
   in.close();
 }
-const static IOFormat CSVFormat(StreamPrecision, DontAlignCols, ", ", "\n");
-template<class Matrix>
-void writeToCSV(const std::string &filename, Matrix &matrix)
-{
-    ofstream file(filename.c_str());
-    file << matrix.format(CSVFormat);
-}
+//const static IOFormat CSVFormat(StreamPrecision, DontAlignCols, ", ", "\n");
+//template<class Matrix>
+//void writeToCSV(const std::string &filename, Matrix &matrix)
+//{
+//    ofstream file(filename.c_str());
+//    file << matrix.format(CSVFormat);
+//}
 } // namespace Eigen
 
 #include "simulationmesh.hpp"
@@ -112,52 +113,54 @@ class SimulationJob {
   } jsonLabels;
 
 public:
-  std::shared_ptr<SimulationMesh> pMesh;
-  std::string label{"empty_job"};
-  std::unordered_map<VertexIndex, std::unordered_set<int>> constrainedVertices;
-  std::unordered_map<VertexIndex, Vector6d> nodalExternalForces;
-  std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
-  SimulationJob(
-      const std::shared_ptr<SimulationMesh> &m, const std::string &label,
-      const std::unordered_map<VertexIndex, std::unordered_set<int>> &cv,
-      const std::unordered_map<VertexIndex, Vector6d> &ef = {},
-      const std::unordered_map<VertexIndex, Eigen::Vector3d> &fd = {})
-      : pMesh(m), label(label), constrainedVertices(cv),
-        nodalExternalForces(ef), nodalForcedDisplacements(fd) {}
-  SimulationJob() {}
-  SimulationJob(const std::string &jsonFilename) { load(jsonFilename); }
+    inline static std::string jsonDefaultFileName = "SimulationJob.json";
+    std::shared_ptr<SimulationMesh> pMesh;
+    std::string label{"empty_job"};
+    std::unordered_map<VertexIndex, std::unordered_set<int>> constrainedVertices;
+    std::unordered_map<VertexIndex, Vector6d> nodalExternalForces;
+    std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
+    SimulationJob(const std::shared_ptr<SimulationMesh> &m,
+                  const std::string &label,
+                  const std::unordered_map<VertexIndex, std::unordered_set<int>> &cv,
+                  const std::unordered_map<VertexIndex, Vector6d> &ef = {},
+                  const std::unordered_map<VertexIndex, Eigen::Vector3d> &fd = {})
+        : pMesh(m), label(label), constrainedVertices(cv), nodalExternalForces(ef),
+          nodalForcedDisplacements(fd)
+    {}
+    SimulationJob() {}
+    SimulationJob(const std::string &jsonFilename) { load(jsonFilename); }
 
-  bool isEmpty()
-  {
-      return constrainedVertices.empty() && nodalExternalForces.empty()
-             && nodalForcedDisplacements.empty() && pMesh == nullptr;
+    bool isEmpty()
+    {
+        return constrainedVertices.empty() && nodalExternalForces.empty()
+               && nodalForcedDisplacements.empty() && pMesh == nullptr;
   }
-  void remap(const std::unordered_map<size_t, size_t> &thisToOtherViMap,
-             SimulationJob &simulationJobMapped) const
+  void remap(const std::unordered_map<size_t, size_t> &sourceToDestinationViMap,
+             SimulationJob &destination_simulationJob) const
   {
-      assert(simulationJobMapped.pMesh->VN() != 0);
-      std::unordered_map<VertexIndex, std::unordered_set<int>> reducedModelFixedVertices;
-      for (auto fullModelFixedVertex : this->constrainedVertices) {
-          reducedModelFixedVertices[thisToOtherViMap.at(fullModelFixedVertex.first)]
-              = fullModelFixedVertex.second;
+      std::unordered_map<VertexIndex, std::unordered_set<int>> destination_fixedVertices;
+      for (const auto &source_fixedVertex : this->constrainedVertices) {
+          destination_fixedVertices[sourceToDestinationViMap.at(source_fixedVertex.first)]
+              = source_fixedVertex.second;
       }
 
-      std::unordered_map<VertexIndex, Vector6d> reducedModelNodalForces;
-      for (auto fullModelNodalForce : this->nodalExternalForces) {
-          reducedModelNodalForces[thisToOtherViMap.at(fullModelNodalForce.first)]
-              = fullModelNodalForce.second;
+      std::unordered_map<VertexIndex, Vector6d> destination_nodalForces;
+      for (const auto &source_nodalForces : this->nodalExternalForces) {
+          destination_nodalForces[sourceToDestinationViMap.at(source_nodalForces.first)]
+              = source_nodalForces.second;
       }
 
-      std::unordered_map<VertexIndex, Eigen::Vector3d> reducedNodalForcedDisplacements;
-      for (auto fullForcedDisplacement : this->nodalForcedDisplacements) {
-          reducedNodalForcedDisplacements[thisToOtherViMap.at(fullForcedDisplacement.first)]
-              = fullForcedDisplacement.second;
+      std::unordered_map<VertexIndex, Eigen::Vector3d> destination_forcedDisplacements;
+      for (const auto &source_forcedDisplacements : this->nodalForcedDisplacements) {
+          destination_forcedDisplacements[sourceToDestinationViMap.at(
+              source_forcedDisplacements.first)]
+              = source_forcedDisplacements.second;
       }
 
-      simulationJobMapped.constrainedVertices = reducedModelFixedVertices;
-      simulationJobMapped.nodalExternalForces = reducedModelNodalForces;
-      simulationJobMapped.label = this->getLabel();
-      simulationJobMapped.nodalForcedDisplacements = reducedNodalForcedDisplacements;
+      destination_simulationJob.constrainedVertices = destination_fixedVertices;
+      destination_simulationJob.nodalExternalForces = destination_nodalForces;
+      destination_simulationJob.label = this->getLabel();
+      destination_simulationJob.nodalForcedDisplacements = destination_forcedDisplacements;
   }
   SimulationJob getCopy() const {
     SimulationJob jobCopy;
@@ -376,7 +379,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
           return;
       }
       std::vector<std::array<double, 3>> nodeColors(pMesh->VN());
-      for (auto fixedVertex : constrainedVertices) {
+      for (const auto &fixedVertex : constrainedVertices) {
           const bool hasRotationalDoFConstrained = fixedVertex.second.contains(3)
                                                    || fixedVertex.second.contains(4)
                                                    || fixedVertex.second.contains(5);
@@ -392,7 +395,8 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
           }
       }
       if (!nodalForcedDisplacements.empty()) {
-          for (std::pair<VertexIndex, Eigen::Vector3d> viDisplPair : nodalForcedDisplacements) {
+          for (const std::pair<VertexIndex, Eigen::Vector3d> &viDisplPair :
+               nodalForcedDisplacements) {
               const VertexIndex vi = viDisplPair.first;
               nodeColors[vi][0] += 1;
               nodeColors[vi][0] /= 2;
@@ -487,6 +491,7 @@ struct SimulationResults
     /*TODO: remove rotationalDisplacementQuaternion since the last three components of the displacments
  * vector contains the same info using euler angles
 */
+    inline const static std::string defaultJsonFilename{"SimulationResults.json"};
     bool converged{false};
     std::shared_ptr<SimulationJob> pJob;
     SimulationHistory history;
@@ -538,7 +543,7 @@ struct SimulationResults
         const std::filesystem::path outputFolderPath = outputFolder.empty()
                                                            ? std::filesystem::current_path()
                                                            : std::filesystem::path(outputFolder);
-        std::cout << "Saving results to:" << outputFolderPath << std::endl;
+        //        std::cout << "Saving results to:" << outputFolderPath << std::endl;
         std::filesystem::path simulationJobOutputFolderPath
             = std::filesystem::path(outputFolderPath).append("SimulationJob");
         std::filesystem::create_directories(simulationJobOutputFolderPath);
@@ -551,6 +556,16 @@ struct SimulationResults
         std::filesystem::create_directories(resultsFolderPath);
         Eigen::write_binary(std::filesystem::path(resultsFolderPath).append(filename).string(), m);
 
+        nlohmann::json json;
+
+        json[GET_VARIABLE_NAME(internalPotentialEnergy)] = internalPotentialEnergy;
+
+        std::filesystem::path jsonFilePath(
+            std::filesystem::path(resultsFolderPath).append(defaultJsonFilename));
+        std::ofstream jsonFile(jsonFilePath.string());
+        jsonFile << json;
+        jsonFile.close();
+
         saveDeformedModel(resultsFolderPath.string());
     }
 
@@ -568,11 +583,198 @@ struct SimulationResults
 
     void load(const std::filesystem::path &loadFromPath, const std::filesystem::path &loadJobFrom)
     {
-        //load job
         pJob->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json").string());
+        load(loadFromPath);
+    }
+    void load(const std::filesystem::path &loadFromPath, const std::shared_ptr<SimulationJob> &pJob)
+    {
+        this->pJob = pJob;
+        load(loadFromPath);
+    }
+    static double computeDistance(
+        const SimulationResults &resultsA,
+        const SimulationResults &resultsB,
+        const std::unordered_map<VertexIndex, VertexIndex> &resultsAToResultsBViMap)
+    {
+        double distance = 0;
+        for (std::pair<int, int> resultsAToResultsBViPair : resultsAToResultsBViMap) {
+            const double vertexToVertexDistance
+                = (resultsA.displacements[resultsAToResultsBViPair.first].getTranslation()
+                   - resultsB.displacements[resultsAToResultsBViPair.second].getTranslation())
+                      .norm();
+            distance += vertexToVertexDistance;
+        }
+        return distance;
+    }
+
+    double computeDistance(const SimulationResults &other,
+                           const std::unordered_map<VertexIndex, VertexIndex> &thisToOtherViMap) const
+    {
+        return computeDistance(*this, other, thisToOtherViMap);
+    }
+
+#ifdef POLYSCOPE_DEFINED
+    void unregister() const
+    {
+        if (!polyscope::hasCurveNetwork(getLabel())) {
+            std::cerr << "No curve network registered with a name: " << getLabel() << std::endl;
+            std::cerr << "Nothing to remove." << std::endl;
+            return;
+        }
+        pJob->unregister(getLabel());
+        polyscope::removeCurveNetwork(getLabel());
+    }
+    polyscope::CurveNetwork *registerForDrawing(
+        const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
+        const bool &shouldEnable = true,
+        const double &desiredRadius = 0.001,
+        //      const double &desiredRadius = 0.0001,
+        const bool &shouldShowFrames = false) const
+    {
+        PolyscopeInterface::init();
+        const std::shared_ptr<SimulationMesh> &mesh = pJob->pMesh;
+        polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh;
+        if (!polyscope::hasStructure(getLabel())) {
+            const auto verts = mesh->getEigenVertices();
+            const auto edges = mesh->getEigenEdges();
+            polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(),
+                                                                               verts,
+                                                                               edges);
+        } else {
+            polyscopeHandle_deformedEdmeMesh = polyscope::getCurveNetwork(getLabel());
+        }
+        polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable);
+        polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, true);
+        if (desiredColor.has_value()) {
+            const glm::vec3 desiredColor_glm(desiredColor.value()[0],
+                                             desiredColor.value()[1],
+                                             desiredColor.value()[2]);
+            polyscopeHandle_deformedEdmeMesh->setColor(desiredColor_glm);
+        }
+        Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3);
+        Eigen::MatrixX3d framesX(mesh->VN(), 3);
+        Eigen::MatrixX3d framesY(mesh->VN(), 3);
+        Eigen::MatrixX3d framesZ(mesh->VN(), 3);
+
+        Eigen::MatrixX3d framesX_initial(mesh->VN(), 3);
+        Eigen::MatrixX3d framesY_initial(mesh->VN(), 3);
+        Eigen::MatrixX3d framesZ_initial(mesh->VN(), 3);
+
+        //      std::unordered_set<int> interfaceNodes{1, 3, 5, 7, 9, 11};
+        //                  std::unordered_set<int> interfaceNodes{3, 7, 11, 15, 19, 23};
+        //      std::unordered_set<int> interfaceNodes{};
+        for (VertexIndex vi = 0; vi < mesh->VN(); vi++) {
+            const Vector6d &nodalDisplacement = displacements[vi];
+            nodalDisplacements.row(vi) = Eigen::Vector3d(nodalDisplacement[0],
+                                                         nodalDisplacement[1],
+                                                         nodalDisplacement[2]);
+            //           Eigen::Quaternion<double> Rx(Eigen::AngleAxis(nodalDisplacement[2],Eigen::Vector3d(1, 0, 0)));
+            //           Eigen::Quaternion<double> Ry(Eigen::AngleAxis(nodalDisplacement[4],Eigen::Vector3d(0, 1, 0)));
+            //           Eigen::Quaternion<double> Rz(Eigen::AngleAxis(nodalDisplacement[5],Eigen::Vector3d(0, 0, 1)));
+            //           Eigen::Quaternion<double> R=Rx*Ry*Rz;
+            //          if (interfaceNodes.contains(vi)) {
+            auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1);
+            auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0);
+            auto deformedFrameX = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(1, 0, 0);
+            framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0],
+                                              deformedFrameX[1],
+                                              deformedFrameX[2]);
+            framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0],
+                                              deformedFrameY[1],
+                                              deformedFrameY[2]);
+            framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0],
+                                              deformedNormal[1],
+                                              deformedNormal[2]);
+            framesX_initial.row(vi) = Eigen::Vector3d(1, 0, 0);
+            framesY_initial.row(vi) = Eigen::Vector3d(0, 1, 0);
+            framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 1);
+            //          } else {
+            //              framesX.row(vi) = Eigen::Vector3d(0, 0, 0);
+            //              framesY.row(vi) = Eigen::Vector3d(0, 0, 0);
+            //              framesZ.row(vi) = Eigen::Vector3d(0, 0, 0);
+            //              framesX_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
+            //              framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
+            //              framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
+            //          }
+        }
+        polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices()
+                                                              + nodalDisplacements);
+
+        const double frameRadius_default = 0.035;
+        const double frameLength_default = 0.035;
+        const bool shouldEnable_default = true;
+        //if(showFramesOn.contains(vi)){
+        auto polyscopeHandle_frameX = polyscopeHandle_deformedEdmeMesh
+                                          ->addNodeVectorQuantity("FrameX", framesX);
+        polyscopeHandle_frameX->setVectorLengthScale(frameLength_default);
+        polyscopeHandle_frameX->setVectorRadius(frameRadius_default);
+        polyscopeHandle_frameX->setVectorColor(
+            /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
+        auto polyscopeHandle_frameY = polyscopeHandle_deformedEdmeMesh
+                                          ->addNodeVectorQuantity("FrameY", framesY);
+        polyscopeHandle_frameY->setVectorLengthScale(frameLength_default);
+        polyscopeHandle_frameY->setVectorRadius(frameRadius_default);
+        polyscopeHandle_frameY->setVectorColor(
+            /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
+        auto polyscopeHandle_frameZ = polyscopeHandle_deformedEdmeMesh
+                                          ->addNodeVectorQuantity("FrameZ", framesZ);
+        polyscopeHandle_frameZ->setVectorLengthScale(frameLength_default);
+        polyscopeHandle_frameZ->setVectorRadius(frameRadius_default);
+        polyscopeHandle_frameZ->setVectorColor(
+            /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
+
+        auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel());
+        if (!polyscopeHandle_initialMesh) {
+            polyscopeHandle_initialMesh = mesh->registerForDrawing();
+        }
+
+        //      auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh
+        //                                                ->addNodeVectorQuantity("FrameX", framesX_initial);
+        //      polyscopeHandle_frameX_initial->setVectorLengthScale(frameLength_default);
+        //      polyscopeHandle_frameX_initial->setVectorRadius(frameRadius_default);
+        //      polyscopeHandle_frameX_initial->setVectorColor(
+        //          /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
+        //      auto polyscopeHandle_frameY_initial = polyscopeHandle_initialMesh
+        //                                                ->addNodeVectorQuantity("FrameY", framesY_initial);
+        //      polyscopeHandle_frameY_initial->setVectorLengthScale(frameLength_default);
+        //      polyscopeHandle_frameY_initial->setVectorRadius(frameRadius_default);
+        //      polyscopeHandle_frameY_initial->setVectorColor(
+        //          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
+        //      auto polyscopeHandle_frameZ_initial = polyscopeHandle_initialMesh
+        //                                                ->addNodeVectorQuantity("FrameZ", framesZ_initial);
+        //      polyscopeHandle_frameZ_initial->setVectorLengthScale(frameLength_default);
+        //      polyscopeHandle_frameZ_initial->setVectorRadius(frameRadius_default);
+        //      polyscopeHandle_frameZ_initial->setVectorColor(
+        //          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
+        //      //}
+        pJob->registerForDrawing(getLabel());
+        //        static bool wasExecuted =false;
+        //        if (!wasExecuted) {
+        //            std::function<void()> callback = [&]() {
+        //                static bool showFrames = shouldShowFrames;
+
+        //                if (ImGui::Checkbox("Show Frames", &showFrames) && showFrames) {
+        //                    polyscopeHandle_frameX->setEnabled(showFrames);
+        //                    polyscopeHandle_frameY->setEnabled(showFrames);
+        //                    polyscopeHandle_frameZ->setEnabled(showFrames);
+        //                }
+        //            };
+
+        //            PolyscopeInterface::addUserCallback(callback);
+        //            wasExecuted = true;
+        //        }
+        return polyscopeHandle_deformedEdmeMesh;
+    }
+#endif
+private:
+    bool load(const std::filesystem::path &loadFromPath)
+    {
+        converged = true; //assuming it has converged
+        assert(pJob != nullptr);
+        //load job
         //Use the first .eigenBin file for loading the displacements
         for (auto const &entry : std::filesystem::recursive_directory_iterator(loadFromPath)) {
-            if (filesystem::is_regular_file(entry) && entry.path().extension() == ".eigenBin") {
+            if (std::filesystem::is_regular_file(entry) && entry.path().extension() == ".eigenBin") {
                 Eigen::MatrixXd displacements_eigen;
                 Eigen::read_binary(entry.path().string(), displacements_eigen);
                 displacements = Utilities::fromEigenMatrix(displacements_eigen);
@@ -589,152 +791,23 @@ struct SimulationResults
                                                    * Eigen::AngleAxisd(displacements[vi][5],
                                                                        Eigen::Vector3d::UnitZ());
         }
+
+        const std::filesystem::path jsonFilepath = std::filesystem::path(loadFromPath)
+                                                       .append(defaultJsonFilename);
+        if (!std::filesystem::exists(jsonFilepath)) {
+            std::cerr << "Simulation results could not be loaded because filepath does "
+                         "not exist:"
+                      << jsonFilepath << std::endl;
+            return false;
+        }
+        std::ifstream ifs(jsonFilepath);
+        nlohmann::json json;
+        ifs >> json;
+        if (json.contains(GET_VARIABLE_NAME(internalPotentialEnergy))) {
+            internalPotentialEnergy = json.at(GET_VARIABLE_NAME(internalPotentialEnergy));
+        }
+        return true;
     }
-
-#ifdef POLYSCOPE_DEFINED
-  void unregister() const {
-    if (!polyscope::hasCurveNetwork(getLabel())) {
-      std::cerr << "No curve network registered with a name: " << getLabel()
-                << std::endl;
-      std::cerr << "Nothing to remove." << std::endl;
-      return;
-    }
-    pJob->unregister(getLabel());
-    polyscope::removeCurveNetwork(getLabel());
-  }
-  polyscope::CurveNetwork *registerForDrawing(
-      const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
-      const bool &shouldEnable = true,
-      const double &desiredRadius = 0.001,
-      //      const double &desiredRadius = 0.0001,
-      const bool &shouldShowFrames = false) const
-  {
-      PolyscopeInterface::init();
-      const std::shared_ptr<SimulationMesh> &mesh = pJob->pMesh;
-      polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh;
-      if (!polyscope::hasStructure(getLabel())) {
-          polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(),
-                                                                             mesh->getEigenVertices(),
-                                                                             mesh->getEigenEdges());
-      }
-      polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable);
-      polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, true);
-      if (desiredColor.has_value()) {
-          const glm::vec3 desiredColor_glm(desiredColor.value()[0],
-                                           desiredColor.value()[1],
-                                           desiredColor.value()[2]);
-          polyscopeHandle_deformedEdmeMesh->setColor(desiredColor_glm);
-      }
-      Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3);
-      Eigen::MatrixX3d framesX(mesh->VN(), 3);
-      Eigen::MatrixX3d framesY(mesh->VN(), 3);
-      Eigen::MatrixX3d framesZ(mesh->VN(), 3);
-
-      Eigen::MatrixX3d framesX_initial(mesh->VN(), 3);
-      Eigen::MatrixX3d framesY_initial(mesh->VN(), 3);
-      Eigen::MatrixX3d framesZ_initial(mesh->VN(), 3);
-
-      //      std::unordered_set<int> interfaceNodes{1, 3, 5, 7, 9, 11};
-      std::unordered_set<int> interfaceNodes{3, 7, 11, 15, 19, 23};
-      for (VertexIndex vi = 0; vi < mesh->VN(); vi++) {
-          const Vector6d &nodalDisplacement = displacements[vi];
-          nodalDisplacements.row(vi) = Eigen::Vector3d(nodalDisplacement[0],
-                                                       nodalDisplacement[1],
-                                                       nodalDisplacement[2]);
-          if (interfaceNodes.contains(vi)) {
-              auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1);
-              auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0);
-              auto deformedFrameX = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(1, 0, 0);
-              framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0],
-                                                deformedFrameX[1],
-                                                deformedFrameX[2]);
-              framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0],
-                                                deformedFrameY[1],
-                                                deformedFrameY[2]);
-              framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0],
-                                                deformedNormal[1],
-                                                deformedNormal[2]);
-              framesX_initial.row(vi) = Eigen::Vector3d(1, 0, 0);
-              framesY_initial.row(vi) = Eigen::Vector3d(0, 1, 0);
-              framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 1);
-          } else {
-              framesX.row(vi) = Eigen::Vector3d(0, 0, 0);
-              framesY.row(vi) = Eigen::Vector3d(0, 0, 0);
-              framesZ.row(vi) = Eigen::Vector3d(0, 0, 0);
-              framesX_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
-              framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
-              framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0);
-          }
-      }
-      polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices()
-                                                            + nodalDisplacements);
-
-      const double frameRadius_default = 0.035;
-      const double frameLength_default = 0.035;
-      const bool shouldEnable_default = true;
-      //if(showFramesOn.contains(vi)){
-      auto polyscopeHandle_frameX = polyscopeHandle_deformedEdmeMesh
-                                        ->addNodeVectorQuantity("FrameX", framesX);
-      polyscopeHandle_frameX->setVectorLengthScale(frameLength_default);
-      polyscopeHandle_frameX->setVectorRadius(frameRadius_default);
-      polyscopeHandle_frameX->setVectorColor(
-          /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
-      auto polyscopeHandle_frameY = polyscopeHandle_deformedEdmeMesh
-                                        ->addNodeVectorQuantity("FrameY", framesY);
-      polyscopeHandle_frameY->setVectorLengthScale(frameLength_default);
-      polyscopeHandle_frameY->setVectorRadius(frameRadius_default);
-      polyscopeHandle_frameY->setVectorColor(
-          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
-      auto polyscopeHandle_frameZ = polyscopeHandle_deformedEdmeMesh
-                                        ->addNodeVectorQuantity("FrameZ", framesZ);
-      polyscopeHandle_frameZ->setVectorLengthScale(frameLength_default);
-      polyscopeHandle_frameZ->setVectorRadius(frameRadius_default);
-      polyscopeHandle_frameZ->setVectorColor(
-          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
-
-      //      auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel());
-      //      if (!polyscopeHandle_initialMesh) {
-      //          polyscopeHandle_initialMesh = mesh->registerForDrawing();
-      //      }
-
-      //      auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh
-      //                                                ->addNodeVectorQuantity("FrameX", framesX_initial);
-      //      polyscopeHandle_frameX_initial->setVectorLengthScale(frameLength_default);
-      //      polyscopeHandle_frameX_initial->setVectorRadius(frameRadius_default);
-      //      polyscopeHandle_frameX_initial->setVectorColor(
-      //          /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0));
-      //      auto polyscopeHandle_frameY_initial = polyscopeHandle_initialMesh
-      //                                                ->addNodeVectorQuantity("FrameY", framesY_initial);
-      //      polyscopeHandle_frameY_initial->setVectorLengthScale(frameLength_default);
-      //      polyscopeHandle_frameY_initial->setVectorRadius(frameRadius_default);
-      //      polyscopeHandle_frameY_initial->setVectorColor(
-      //          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0));
-      //      auto polyscopeHandle_frameZ_initial = polyscopeHandle_initialMesh
-      //                                                ->addNodeVectorQuantity("FrameZ", framesZ_initial);
-      //      polyscopeHandle_frameZ_initial->setVectorLengthScale(frameLength_default);
-      //      polyscopeHandle_frameZ_initial->setVectorRadius(frameRadius_default);
-      //      polyscopeHandle_frameZ_initial->setVectorColor(
-      //          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
-      //      //}
-      pJob->registerForDrawing(getLabel());
-      static bool wasExecuted = false;
-      if (!wasExecuted) {
-          std::function<void()> callback = [&]() {
-              static bool showFrames = shouldShowFrames;
-
-              if (ImGui::Checkbox("Show Frames", &showFrames) && showFrames) {
-                  polyscopeHandle_frameX->setEnabled(showFrames);
-                  polyscopeHandle_frameY->setEnabled(showFrames);
-                  polyscopeHandle_frameZ->setEnabled(showFrames);
-              }
-          };
-
-          PolyscopeInterface::addUserCallback(callback);
-          wasExecuted = true;
-      }
-      return polyscopeHandle_deformedEdmeMesh;
-  }
-#endif
 };
 
 #endif // SIMULATIONHISTORY_HPP

simulationhistoryplotter.hpp

@@ -1,8 +1,8 @@
 #ifndef SIMULATIONHISTORYPLOTTER_HPP
 #define SIMULATIONHISTORYPLOTTER_HPP
 
-#include "simulationmesh.hpp"
 #include "simulation_structs.hpp"
+#include "simulationmesh.hpp"
 #include "utilities.hpp"
 #include <algorithm>
 #include <matplot/matplot.h>
@@ -15,27 +15,25 @@ struct SimulationResultsReporter {
 
   void writeStatistics(const SimulationResults &results,
                        const std::string &reportFolderPath) {
+      std::ofstream file;
+      file.open(std::filesystem::path(reportFolderPath).append("results.txt").string());
 
-    ofstream file;
-    file.open(
-        std::filesystem::path(reportFolderPath).append("results.txt").string());
+      const size_t numberOfSteps = results.history.numberOfSteps;
+      file << "Number of steps " << numberOfSteps << "\n";
 
-    const size_t numberOfSteps = results.history.numberOfSteps;
-    file << "Number of steps " << numberOfSteps << "\n";
+      //    file << "Force threshold used " << 1000 << "\n";
 
-    //    file << "Force threshold used " << 1000 << "\n";
-
-    //    assert(numberOfSteps == results.history.potentialEnergy.size() &&
-    //           numberOfSteps == results.history.residualForces.size());
-    // Write kinetic energies
-    const SimulationHistory &history = results.history;
-    if (!history.kineticEnergy.empty()) {
-      file << "Kinetic energies"
-           << "\n";
-      for (size_t step = 0; step < numberOfSteps; step++) {
-        file << history.kineticEnergy[step] << "\n";
-      }
-      file << "\n";
+      //    assert(numberOfSteps == results.history.potentialEnergy.size() &&
+      //           numberOfSteps == results.history.residualForces.size());
+      // Write kinetic energies
+      const SimulationHistory &history = results.history;
+      if (!history.kineticEnergy.empty()) {
+          file << "Kinetic energies"
+               << "\n";
+          for (size_t step = 0; step < numberOfSteps; step++) {
+              file << history.kineticEnergy[step] << "\n";
+          }
+          file << "\n";
     }
 
     if (!history.logResidualForces.empty()) {
@@ -114,93 +112,81 @@ struct SimulationResultsReporter {
       if (YvaluesToPlot.size() < 2) {
           return;
       }
-      auto x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size());
-      std::vector<double> colors(x.size(), 0.5);
-      std::vector<double> markerSizes(x.size(), 5);
+      std::vector<double> colors(YvaluesToPlot.size(), 0.5);
+      std::vector<double> markerSizes(YvaluesToPlot.size(), 5);
       if (!markPoints.empty()) {
           for (const auto pointIndex : markPoints) {
               colors[pointIndex] = 0.9;
               markerSizes[pointIndex] = 14;
           }
       }
+      std::vector<double> x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size());
       createPlot(xLabel, yLabel, x, YvaluesToPlot, markerSizes, colors, saveTo);
-      //      matplot::figure(true);
-      //      matplot::hold(matplot::on);
-
-      //        ->marker_indices(history.redMarks)
-      //          ->marker_indices(truncatedRedMarks)
-      //        .marker_color("r")
-      //                ->marker_size(1)
-      //      ;
-      //      auto greenMarksY = matplot::transform(history.greenMarks,
-      //                                            [&](auto x) { return YvaluesToPlot[x]; });
-      //      matplot::scatter(history.greenMarks, greenMarksY)->color("green").marker_size(10);
-      //      matplot::hold(matplot::off);
   }
 
   void createPlots(const SimulationHistory &history,
                    const std::string &reportFolderPath,
-                   const std::string &graphSuffix) {
-    const auto graphsFolder =
-        std::filesystem::path(reportFolderPath).append("Graphs");
-    std::filesystem::remove_all(graphsFolder);
-    std::filesystem::create_directory(graphsFolder.string());
+                   const std::string &graphSuffix)
+  {
+      const auto graphsFolder = std::filesystem::path(reportFolderPath).append("Graphs");
+      std::filesystem::remove_all(graphsFolder);
+      std::filesystem::create_directory(graphsFolder.string());
 
-    if (!history.kineticEnergy.empty()) {
-        createPlot("Number of Iterations",
-                   "Log of Kinetic Energy log",
-                   history.kineticEnergy,
-                   std::filesystem::path(graphsFolder)
-                       .append("KineticEnergyLog_" + graphSuffix + ".png")
-                       .string(),
-                   history.redMarks);
-    }
+      if (!history.kineticEnergy.empty()) {
+          createPlot("Number of Iterations",
+                     "Log of Kinetic Energy log",
+                     history.kineticEnergy,
+                     std::filesystem::path(graphsFolder)
+                         .append("KineticEnergyLog_" + graphSuffix + ".png")
+                         .string(),
+                     history.redMarks);
+      }
 
-    if (!history.logResidualForces.empty()) {
-        createPlot("Number of Iterations",
-                   "Residual Forces norm log",
-                   history.logResidualForces,
-                   std::filesystem::path(graphsFolder)
-                       .append("ResidualForcesLog_" + graphSuffix + ".png")
-                       .string(),
-                   history.redMarks);
-    }
+      if (!history.logResidualForces.empty()) {
+          createPlot("Number of Iterations",
+                     "Residual Forces norm log",
+                     history.logResidualForces,
+                     std::filesystem::path(graphsFolder)
+                         .append("ResidualForcesLog_" + graphSuffix + ".png")
+                         .string(),
+                     history.redMarks);
+      }
 
-    if (!history.potentialEnergies.empty()) {
-        createPlot("Number of Iterations",
-                   "Potential energy",
-                   history.potentialEnergies,
-                   std::filesystem::path(graphsFolder)
-                       .append("PotentialEnergy_" + graphSuffix + ".png")
-                       .string(),
-                   history.redMarks);
-    }
-    if (!history.residualForcesMovingAverage.empty()) {
-        createPlot("Number of Iterations",
-                   "Residual forces moving average",
-                   history.residualForcesMovingAverage,
-                   std::filesystem::path(graphsFolder)
-                       .append("ResidualForcesMovingAverage_" + graphSuffix + ".png")
-                       .string(),
-                   history.redMarks);
-    }
-    //    if (!history.residualForcesMovingAverageDerivativesLog.empty()) {
-    //        createPlot("Number of Iterations",
-    //                   "Residual forces moving average derivative log",
-    //                   history.residualForcesMovingAverageDerivativesLog,
-    //                   std::filesystem::path(graphsFolder)
-    //                       .append("ResidualForcesMovingAverageDerivativesLog_" + graphSuffix + ".png")
-    //                       .string());
-    //    }
-    if (!history.sumOfNormalizedDisplacementNorms.empty()) {
-        createPlot("Number of Iterations",
-                   "Sum of normalized displacement norms",
-                   history.sumOfNormalizedDisplacementNorms,
-                   std::filesystem::path(graphsFolder)
-                       .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png")
-                       .string(),
-                   history.redMarks);
-    }
+      if (!history.potentialEnergies.empty()) {
+          createPlot("Number of Iterations",
+                     "Potential energy",
+                     history.potentialEnergies,
+                     std::filesystem::path(graphsFolder)
+                         .append("PotentialEnergy_" + graphSuffix + ".png")
+                         .string(),
+                     history.redMarks);
+      }
+      if (!history.residualForcesMovingAverage.empty()) {
+          createPlot("Number of Iterations",
+                     "Residual forces moving average",
+                     history.residualForcesMovingAverage,
+                     std::filesystem::path(graphsFolder)
+                         .append("ResidualForcesMovingAverage_" + graphSuffix + ".png")
+                         .string(),
+                     history.redMarks);
+      }
+      //    if (!history.residualForcesMovingAverageDerivativesLog.empty()) {
+      //        createPlot("Number of Iterations",
+      //                   "Residual forces moving average derivative log",
+      //                   history.residualForcesMovingAverageDerivativesLog,
+      //                   std::filesystem::path(graphsFolder)
+      //                       .append("ResidualForcesMovingAverageDerivativesLog_" + graphSuffix + ".png")
+      //                       .string());
+      //    }
+      if (!history.sumOfNormalizedDisplacementNorms.empty()) {
+          createPlot("Number of Iterations",
+                     "Sum of normalized displacement norms",
+                     history.sumOfNormalizedDisplacementNorms,
+                     std::filesystem::path(graphsFolder)
+                         .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png")
+                         .string(),
+                     history.redMarks);
+      }
   }
 };
 

simulationmesh.cpp

@@ -15,8 +15,8 @@ SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) {
 
   elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(
       *this, std::string("Elements"));
-  nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(
-      *this, std::string("Nodes"));
+  nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(*this,
+                                                                           std::string("Nodes"));
   initializeNodes();
   initializeElements();
 }
@@ -198,7 +198,6 @@ void SimulationMesh::reset() {
 }
 
 void SimulationMesh::initializeElements() {
-  computeElementalProperties();
   for (const EdgeType &e : edge) {
     Element &element = elements[e];
     element.ei = getIndex(e);
@@ -243,8 +242,6 @@ void SimulationMesh::updateElementalLengths() {
     const vcg::Segment3<double> s(p0, p1);
     const double elementLength = s.Length();
     elements[e].length = elementLength;
-    int i = 0;
-    i++;
   }
 }
 
@@ -260,7 +257,6 @@ void SimulationMesh::setBeamCrossSection(
     const CrossSectionType &beamDimensions) {
   for (size_t ei = 0; ei < EN(); ei++) {
     elements[ei].dimensions = beamDimensions;
-    elements[ei].computeDimensionsProperties(beamDimensions);
     elements[ei].updateRigidity();
   }
 }
@@ -288,84 +284,85 @@ std::vector<ElementMaterial> SimulationMesh::getBeamMaterial() {
   return beamMaterial;
 }
 
-bool SimulationMesh::load(const string &plyFilename) {
-  this->Clear();
-  //  assert(plyFileHasAllRequiredFields(plyFilename));
-  // Load the ply file
-  //  VCGEdgeMesh::PerEdgeAttributeHandle<CrossSectionType> handleBeamDimensions =
-  //      vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<
-  //          CrossSectionType>(*this, plyPropertyBeamDimensionsID);
-  //  VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial> handleBeamMaterial =
-  //      vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<ElementMaterial>(
-  //          *this, plyPropertyBeamMaterialID);
-  //  nanoply::NanoPlyWrapper<SimulationMesh>::CustomAttributeDescriptor
-  //      customAttrib;
-  //  customAttrib.GetMeshAttrib(plyFilename);
-  //  customAttrib.AddEdgeAttribDescriptor<CrossSectionType, double, 2>(
-  //      plyPropertyBeamDimensionsID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
-  //  /*FIXME: Since I allow CrossSectionType to take two types I should export the
-  //   * type as well such that that when loaded the correct type of cross section
-  //   * is used.
-  //   */
-  //  customAttrib.AddEdgeAttribDescriptor<vcg::Point2d, double, 2>(
-  //      plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
+bool SimulationMesh::load(const std::string &plyFilename)
+{
+    this->Clear();
+    //  assert(plyFileHasAllRequiredFields(plyFilename));
+    // Load the ply file
+    //  VCGEdgeMesh::PerEdgeAttributeHandle<CrossSectionType> handleBeamDimensions =
+    //      vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<
+    //          CrossSectionType>(*this, plyPropertyBeamDimensionsID);
+    //  VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial> handleBeamMaterial =
+    //      vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<ElementMaterial>(
+    //          *this, plyPropertyBeamMaterialID);
+    //  nanoply::NanoPlyWrapper<SimulationMesh>::CustomAttributeDescriptor
+    //      customAttrib;
+    //  customAttrib.GetMeshAttrib(plyFilename);
+    //  customAttrib.AddEdgeAttribDescriptor<CrossSectionType, double, 2>(
+    //      plyPropertyBeamDimensionsID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
+    //  /*FIXME: Since I allow CrossSectionType to take two types I should export the
+    //   * type as well such that that when loaded the correct type of cross section
+    //   * is used.
+    //   */
+    //  customAttrib.AddEdgeAttribDescriptor<vcg::Point2d, double, 2>(
+    //      plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
 
-  //  VCGEdgeMesh::PerEdgeAttributeHandle<std::array<double, 6>>
-  //      handleBeamProperties =
-  //          vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<
-  //              std::array<double, 6>>(*this, plyPropertyBeamProperties);
-  //  customAttrib.AddEdgeAttribDescriptor<std::array<double, 6>, double, 6>(
-  //      plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
+    //  VCGEdgeMesh::PerEdgeAttributeHandle<std::array<double, 6>>
+    //      handleBeamProperties =
+    //          vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<
+    //              std::array<double, 6>>(*this, plyPropertyBeamProperties);
+    //  customAttrib.AddEdgeAttribDescriptor<std::array<double, 6>, double, 6>(
+    //      plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
 
-  //  VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial>
-  //  handleBeamRigidityContants;
-  //  customAttrib.AddEdgeAttribDescriptor<vcg::Point4f, float, 4>(
-  //      plyPropertyBeamRigidityConstantsID, nanoply::NNP_LIST_INT8_FLOAT32,
-  //      nullptr);
-  //  unsigned int mask = 0;
-  //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
-  //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
-  //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
-  //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
-  //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_MESHATTRIB;
-  if (!load(plyFilename)) {
-      return false;
-  }
-
-  //  elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(
-  //      *this, std::string("Elements"));
-  //  nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(
-  //      *this, std::string("Nodes"));
-  vcg::tri::UpdateTopology<SimulationMesh>::VertexEdge(*this);
-  initializeNodes();
-  initializeElements();
-  updateEigenEdgeAndVertices();
-
-  //  if (!handleBeamProperties._handle->data.empty()) {
-  //    for (size_t ei = 0; ei < EN(); ei++) {
-  //      elements[ei] =
-  //      Element::Properties(handleBeamProperties[ei]);
-  //      elements[ei].updateRigidity();
-  //    }
-  //  }
-  //  for (size_t ei = 0; ei < EN(); ei++) {
-  //    elements[ei].setDimensions(handleBeamDimensions[ei]);
-  //    elements[ei].setMaterial(handleBeamMaterial[ei]);
-  //    elements[ei].updateRigidity();
-  //  }
-
-  bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001;
-  if (normalsAreAbsent) {
-    CoordType normalVector(0, 0, 1);
-    std::cout << "Warning: Normals are missing from " << plyFilename
-              << ". Added normal vector:" << toString(normalVector)
-              << std::endl;
-    for (auto &v : vert) {
-      v.N() = normalVector;
+    //  VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial>
+    //  handleBeamRigidityContants;
+    //  customAttrib.AddEdgeAttribDescriptor<vcg::Point4f, float, 4>(
+    //      plyPropertyBeamRigidityConstantsID, nanoply::NNP_LIST_INT8_FLOAT32,
+    //      nullptr);
+    //  unsigned int mask = 0;
+    //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
+    //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
+    //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
+    //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
+    //  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_MESHATTRIB;
+    if (!VCGEdgeMesh::load(plyFilename)) {
+        return false;
     }
-  }
 
-  return true;
+    //  elements = vcg::tri::Allocator<SimulationMesh>::GetPerEdgeAttribute<Element>(
+    //      *this, std::string("Elements"));
+    //  nodes = vcg::tri::Allocator<SimulationMesh>::GetPerVertexAttribute<Node>(
+    //      *this, std::string("Nodes"));
+    vcg::tri::UpdateTopology<SimulationMesh>::VertexEdge(*this);
+    initializeNodes();
+    initializeElements();
+    setBeamMaterial(0.3, 1 * 1e9);
+    updateEigenEdgeAndVertices();
+
+    //  if (!handleBeamProperties._handle->data.empty()) {
+    //    for (size_t ei = 0; ei < EN(); ei++) {
+    //      elements[ei] =
+    //      Element::Properties(handleBeamProperties[ei]);
+    //      elements[ei].updateRigidity();
+    //    }
+    //  }
+    //  for (size_t ei = 0; ei < EN(); ei++) {
+    //    elements[ei].setDimensions(handleBeamDimensions[ei]);
+    //    elements[ei].setMaterial(handleBeamMaterial[ei]);
+    //    elements[ei].updateRigidity();
+    //  }
+
+    bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001;
+    if (normalsAreAbsent) {
+        CoordType normalVector(0, 0, 1);
+        std::cout << "Warning: Normals are missing from " << plyFilename
+                  << ". Added normal vector:" << toString(normalVector) << std::endl;
+        for (auto &v : vert) {
+            v.N() = normalVector;
+        }
+    }
+
+    return true;
 }
 
 bool SimulationMesh::save(const std::string &plyFilename)
@@ -398,7 +395,7 @@ bool SimulationMesh::save(const std::string &plyFilename)
     //        != 1) {
     //        return false;
     //    }
-    if (!save(plyFilename)) {
+    if (!VCGEdgeMesh::save(plyFilename)) {
         return false;
     }
     return true;
@@ -456,60 +453,59 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1,
   return angle;
 }
 
-void Element::computeMaterialProperties(const ElementMaterial &material) {
-    G = material.youngsModulus / (2 * (1 + material.poissonsRatio));
-}
+//void Element::computeMaterialProperties(const ElementMaterial &material) {
+//    G = material.youngsModulus / (2 * (1 + material.poissonsRatio));
+//}
 
-void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A)
-{
-    A = dimensions.b * dimensions.h;
-}
+//void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A)
+//{
+//    A = dimensions.b * dimensions.h;
+//}
 
-void Element::computeDimensionsProperties(
-    const RectangularBeamDimensions &dimensions) {
-  assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions));
-  computeCrossSectionArea(dimensions, A);
-  computeMomentsOfInertia(dimensions, inertia);
-}
+//void Element::computeDimensionsProperties(
+//    const RectangularBeamDimensions &dimensions) {
+//  assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions));
+//  computeCrossSectionArea(dimensions, A);
+//  computeMomentsOfInertia(dimensions, dimensions.inertia);
+//}
 
-void Element::computeDimensionsProperties(
-    const CylindricalBeamDimensions &dimensions) {
-  assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
-  A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2));
-  inertia.I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64;
-  inertia.I3 = inertia.I2;
-  inertia.J = inertia.I2 + inertia.I3;
-}
+//void Element::computeDimensionsProperties(
+//    const CylindricalBeamDimensions &dimensions) {
+//  assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
+//  A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2));
+//  dimensions.inertia.I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64;
+//  dimensions.inertia.I3 = dimensions.inertia.I2;
+//  dimensions.inertia.J = dimensions.inertia.I2 + dimensions.inertia.I3;
+//}
 
 void Element::setDimensions(const CrossSectionType &dimensions) {
   this->dimensions = dimensions;
-  computeDimensionsProperties(dimensions);
+  assert(this->dimensions.A == dimensions.A);
+  //  computeDimensionsProperties(dimensions);
   updateRigidity();
 }
 
-void Element::setMaterial(const ElementMaterial &material) {
-  this->material = material;
-  computeMaterialProperties(material);
-  updateRigidity();
+void Element::setMaterial(const ElementMaterial &material)
+{
+    this->material = material;
+    //  computeMaterialProperties(material);
+    updateRigidity();
 }
 
 double Element::getMass(const double &materialDensity)
 {
-    const double beamVolume = A * length;
+    const double beamVolume = dimensions.A * length;
     return beamVolume * materialDensity;
 }
 
-void Element::computeMomentsOfInertia(const RectangularBeamDimensions &dimensions,
-                                      Element::MomentsOfInertia &inertia)
-{
-    inertia.I2 = dimensions.b * std::pow(dimensions.h, 3) / 12;
-    inertia.I3 = dimensions.h * std::pow(dimensions.b, 3) / 12;
-    inertia.J = inertia.I2 + inertia.I3;
-}
-
 void Element::updateRigidity() {
-  rigidity.axial = material.youngsModulus * A / initialLength;
-  rigidity.torsional = G * inertia.J / initialLength;
-  rigidity.firstBending = 2 * material.youngsModulus * inertia.I2 / initialLength;
-  rigidity.secondBending = 2 * material.youngsModulus * inertia.I3 / initialLength;
+    //    assert(initialLength != 0);
+    rigidity.axial = material.youngsModulus * dimensions.A / initialLength;
+    //    assert(rigidity.axial != 0);
+    rigidity.torsional = material.G * dimensions.inertia.J / initialLength;
+    //    assert(rigidity.torsional != 0);
+    rigidity.firstBending = 2 * material.youngsModulus * dimensions.inertia.I2 / initialLength;
+    //    assert(rigidity.firstBending != 0);
+    rigidity.secondBending = 2 * material.youngsModulus * dimensions.inertia.I3 / initialLength;
+    //    assert(rigidity.secondBending != 0);
 }

simulationmesh.hpp

@@ -1,9 +1,9 @@
 #ifndef SIMULATIONMESH_HPP
 #define SIMULATIONMESH_HPP
 
-#include "Eigen/Dense"
 #include "edgemesh.hpp"
 #include "trianglepatterngeometry.hpp"
+#include <Eigen/Dense>
 
 //extern bool drawGlobal;
 struct Element;
@@ -34,7 +34,7 @@ public:
 
   std::vector<VCGEdgeMesh::EdgePointer>
   getIncidentElements(const VCGEdgeMesh::VertexType &v);
-  virtual bool load(const string &plyFilename);
+  virtual bool load(const std::string &plyFilename);
   std::vector<CrossSectionType> getBeamDimensions();
   std::vector<ElementMaterial> getBeamMaterial();
 
@@ -62,26 +62,14 @@ struct Element {
 
   CrossSectionType dimensions;
   ElementMaterial material;
-  double G{0};
-  double A{0};  // cross sectional area
 
   void computeMaterialProperties(const ElementMaterial &material);
-  void computeDimensionsProperties(const RectangularBeamDimensions &dimensions);
-  void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions);
+  //  void computeDimensionsProperties(const RectangularBeamDimensions &dimensions);
+  //  void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions);
   void setDimensions(const CrossSectionType &dimensions);
   void setMaterial(const ElementMaterial &material);
   double getMass(const double &matrialDensity);
 
-  struct MomentsOfInertia
-  {
-      double I2{0}; // second moment of inertia
-      double I3{0}; // third moment of inertia
-      double J{0};  // torsional constant (polar moment of inertia)
-  } inertia;
-
-  static void computeMomentsOfInertia(const RectangularBeamDimensions &dimensions,
-                                      MomentsOfInertia &inertia);
-
   struct LocalFrame {
     VectorType t1;
     VectorType t2;

trianglepatterngeometry.cpp

@@ -86,39 +86,43 @@ std::vector<vcg::Point3d> PatternGeometry::getVertices() const {
   return verts;
 }
 
-PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern) {
+PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern)
+{
+    const size_t fanSize = PatternGeometry().getFanSize();
+    PatternGeometry fan(createFan(pattern));
+    PatternGeometry tile(fan);
 
-  const size_t fanSize = PatternGeometry().getFanSize();
-  PatternGeometry fan(createFan(pattern));
-  PatternGeometry tile(fan);
+    if (fanSize % 2 == 1) {
+        vcg::Matrix44d R;
+        auto rotationAxis = vcg::Point3d(0, 0, 1);
+        R.SetRotateDeg(180, rotationAxis);
+        vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, R);
+    }
+    vcg::Matrix44d T;
+    //    const double centerAngle = 2 * M_PI / fanSize;
+    //    const double triangleHeight = std::sin((M_PI - centerAngle) / 2)
+    //                                  * PatternGeometry().triangleEdgeSize;
+    vcg::tri::UpdateBounding<PatternGeometry>::Box(fan);
+    const double triangleHeight = fan.bbox.DimY() / 2;
+    T.SetTranslate(0, -2 * triangleHeight, 0);
+    vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, T);
+    //    fan.registerForDrawing();
+    //    polyscope::show();
 
-  if (fanSize % 2 == 1) {
-    vcg::Matrix44d R;
-    auto rotationAxis = vcg::Point3d(0, 0, 1);
-    R.SetRotateDeg(180, rotationAxis);
-    vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, R);
-  }
-  vcg::Matrix44d T;
-  const double centerAngle = 2 * M_PI / fanSize;
-  const double triangleHeight =
-      std::sin((M_PI - centerAngle) / 2) * PatternGeometry().triangleEdgeSize;
-  T.SetTranslate(0, -2 * triangleHeight, 0);
-  vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, T);
+    PatternGeometry fanOfFan = createFan(fan);
+    vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(tile, fanOfFan);
+    vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(tile, 0.0000005);
+    vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(tile);
+    vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(tile);
+    vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(tile);
 
-  PatternGeometry fanOfFan = createFan(fan);
-  vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(tile, fanOfFan);
-  vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(tile, 0.0000005);
-  vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(tile);
-  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(tile);
-  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(tile);
+    //  for (size_t vi = 0; vi < pattern.vn; vi++) {
+    //    tile.vert[vi].C() = vcg::Color4b::Blue;
+    //  }
 
-  //  for (size_t vi = 0; vi < pattern.vn; vi++) {
-  //    tile.vert[vi].C() = vcg::Color4b::Blue;
-  //  }
-
-  tile.setLabel("tilled_" + pattern.getLabel());
-  tile.updateEigenEdgeAndVertices();
-  return tile;
+    tile.setLabel("tilled_" + pattern.getLabel());
+    tile.updateEigenEdgeAndVertices();
+    return tile;
 }
 
 PatternGeometry PatternGeometry::createFan(PatternGeometry &pattern) {
@@ -145,9 +149,9 @@ PatternGeometry::PatternGeometry(PatternGeometry &other) {
   vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
 }
 
-bool PatternGeometry::load(const string &plyFilename)
+bool PatternGeometry::load(const std::filesystem::__cxx11::path &meshFilePath)
 {
-    if (!VCGEdgeMesh::load(plyFilename)) {
+    if (!VCGEdgeMesh::load(meshFilePath)) {
         return false;
     }
     vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
@@ -777,6 +781,22 @@ double PatternGeometry::computeBaseTriangleHeight() const
     return vcg::Distance(vert[0].cP(), vert[interfaceNodeIndex].cP());
 }
 
+void PatternGeometry::deleteDanglingVertices()
+{
+    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> pu;
+    VCGEdgeMesh::deleteDanglingVertices(pu);
+    if (!pu.remap.empty()) {
+        interfaceNodeIndex
+            = pu.remap[interfaceNodeIndex]; //TODO:Could this be automatically be determined?
+    }
+}
+
+void PatternGeometry::deleteDanglingVertices(
+    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu)
+{
+    VCGEdgeMesh::deleteDanglingVertices(pu);
+}
+
 vcg::Triangle3<double> PatternGeometry::getBaseTriangle() const
 {
     return baseTriangle;
@@ -814,6 +834,7 @@ PatternGeometry::PatternGeometry(
   updateEigenEdgeAndVertices();
 }
 
+//TODO: refactor such that it takes as an input a single pattern and tiles into the desired faces without requiring the each face will contain a single pattern
 std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
     std::vector<ConstPatternGeometry> &patterns,
     const std::vector<int> &connectToNeighborsVi,

trianglepatterngeometry.hpp

@@ -46,7 +46,7 @@ private:
    * vcglib interface.
    * */
   PatternGeometry(PatternGeometry &other);
-  bool load(const std::string &plyFilename) override;
+  bool load(const std::filesystem::path &meshFilePath) override;
   void add(const std::vector<vcg::Point3d> &vertices);
   void add(const std::vector<vcg::Point2i> &edges);
   void add(const std::vector<vcg::Point3d> &vertices, const std::vector<vcg::Point2i> &edges);
@@ -124,6 +124,9 @@ private:
       std::vector<std::vector<size_t>> &perPatternIndexTiledPatternEdgeIndex);
   std::unordered_set<VertexIndex> getInterfaceNodes(const std::vector<size_t> &numberOfNodesPerSlot);
   bool isInterfaceConnected(const std::unordered_set<VertexIndex> &interfaceNodes);
+  void deleteDanglingVertices() override;
+  void deleteDanglingVertices(
+      vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu) override;
 };
 
 #endif // FLATPATTERNGEOMETRY_HPP

utilities.hpp

@@ -2,12 +2,14 @@
 #define UTILITIES_H
 
 #include <Eigen/Dense>
-#include <filesystem>
-#include <fstream>
-#include <regex>
-#include <iterator>
 #include <algorithm>
 #include <array>
+#include <chrono>
+#include <filesystem>
+#include <fstream>
+#include <iterator>
+#include <numeric>
+#include <regex>
 
 #define GET_VARIABLE_NAME(Variable) (#Variable)
 
@@ -140,17 +142,121 @@ struct Vector6d : public std::array<double, 6> {
 };
 
 namespace Utilities {
-inline void parseIntegers(const std::string &str, std::vector<size_t> &result) {
-  typedef std::regex_iterator<std::string::const_iterator> re_iterator;
-  typedef re_iterator::value_type re_iterated;
+inline bool compareNat(const std::string &a, const std::string &b)
+{
+    if (a.empty())
+        return true;
+    if (b.empty())
+        return false;
+    if (std::isdigit(a[0]) && !std::isdigit(b[0]))
+        return true;
+    if (!std::isdigit(a[0]) && std::isdigit(b[0]))
+        return false;
+    if (!std::isdigit(a[0]) && !std::isdigit(b[0])) {
+        if (std::toupper(a[0]) == std::toupper(b[0]))
+            return compareNat(a.substr(1), b.substr(1));
+        return (std::toupper(a[0]) < std::toupper(b[0]));
+    }
 
-  std::regex re("(\\d+)");
+    // Both strings begin with digit --> parse both numbers
+    std::istringstream issa(a);
+    std::istringstream issb(b);
+    int ia, ib;
+    issa >> ia;
+    issb >> ib;
+    if (ia != ib)
+        return ia < ib;
 
-  re_iterator rit(str.begin(), str.end(), re);
-  re_iterator rend;
+    // Numbers are the same --> remove numbers and recurse
+    std::string anew, bnew;
+    std::getline(issa, anew);
+    std::getline(issb, bnew);
+    return (compareNat(anew, bnew));
+}
 
-  std::transform(rit, rend, std::back_inserter(result),
-                 [](const re_iterated &it) { return std::stoi(it[1]); });
+inline std::string_view leftTrimSpaces(const std::string_view& str)
+{
+std::string_view trimmedString=str;
+    const auto pos(str.find_first_not_of(" \t\n\r\f\v"));
+    trimmedString.remove_prefix(std::min(pos, trimmedString.length()));
+    return trimmedString;
+}
+
+inline std::string_view rightTrimSpaces(const std::string_view& str)
+{
+std::string_view trimmedString=str;
+    const auto pos(trimmedString.find_last_not_of(" \t\n\r\f\v"));
+    trimmedString.remove_suffix(std::min(trimmedString.length() - pos - 1,trimmedString.length()));
+    return trimmedString;
+}
+
+inline std::string_view trimLeftAndRightSpaces(std::string_view str)
+{
+std::string_view trimmedString=str;
+    trimmedString = leftTrimSpaces(trimmedString);
+    trimmedString = rightTrimSpaces(trimmedString);
+    return trimmedString;
+}
+
+template<typename InputIt>
+inline void normalize(InputIt itBegin, InputIt itEnd)
+{
+    const auto squaredSumOfElements = std::accumulate(itBegin,
+                                                      itEnd,
+                                                      0.0,
+                                                      [](const auto &sum, const auto &el) {
+                                                          return sum + el * el;
+                                                      });
+    assert(squaredSumOfElements != 0);
+    std::transform(itBegin, itEnd, itBegin, [&](auto &element) {
+        return element / std::sqrt(squaredSumOfElements);
+    });
+}
+
+inline std::vector<std::string> split(const std::string& text, std::string delim)
+{
+    std::vector<std::string> vec;
+    size_t pos = 0, prevPos = 0;
+    while (1) {
+        pos = text.find(delim, prevPos);
+        if (pos == std::string::npos) {
+            vec.push_back(text.substr(prevPos));
+            return vec;
+        }
+
+        vec.push_back(text.substr(prevPos, pos - prevPos));
+        prevPos = pos + delim.length();
+    }
+}
+
+inline std::string toString(const std::vector<std::vector<int>> &vv)
+{
+    std::string s;
+    s.append("{");
+    for (const std::vector<int> &v : vv) {
+        s.append("{");
+        for (const int &i : v) {
+            s.append(std::to_string(i) + ",");
+        }
+        s.pop_back();
+        s.append("}");
+    }
+    s.append("}");
+    return s;
+}
+inline void parseIntegers(const std::string &str, std::vector<size_t> &result)
+{
+    typedef std::regex_iterator<std::string::const_iterator> re_iterator;
+    typedef re_iterator::value_type re_iterated;
+
+    std::regex re("(\\d+)");
+
+    re_iterator rit(str.begin(), str.end(), re);
+    re_iterator rend;
+
+    std::transform(rit, rend, std::back_inserter(result), [](const re_iterated &it) {
+        return std::stoi(it[1]);
+    });
 }
 
 inline Eigen::MatrixXd toEigenMatrix(const std::vector<Vector6d> &v) {
@@ -179,7 +285,6 @@ inline std::vector<Vector6d> fromEigenMatrix(const Eigen::MatrixXd &m)
 
     return v;
 }
-
 // std::string convertToLowercase(const std::string &s) {
 //  std::string lowercase;
 //  std::transform(s.begin(), s.end(), lowercase.begin(),
@@ -206,6 +311,19 @@ inline std::vector<Vector6d> fromEigenMatrix(const Eigen::MatrixXd &m)
 //  return true;
 //}
 
+inline std::filesystem::path getFilepathWithExtension(const std::filesystem::path &folderPath,
+                                                      const std::string &extension)
+{
+    for (const std::filesystem::directory_entry &dirEntry :
+         std::filesystem::directory_iterator(folderPath)) {
+        if (dirEntry.is_regular_file() && std::filesystem::path(dirEntry).extension() == extension) {
+            return std::filesystem::path(dirEntry);
+        }
+    }
+
+    return "";
+}
+
 } // namespace Utilities
 
 #ifdef POLYSCOPE_DEFINED
@@ -226,7 +344,7 @@ inline void mainCallback()
                                // instead of full width. Must have
                                // matching PopItemWidth() below.
 
-    for (std::function<void()> userCallback : globalPolyscopeData.userCallbacks) {
+    for (std::function<void()> &userCallback : globalPolyscopeData.userCallbacks) {
         userCallback();
     }
 
@@ -257,7 +375,8 @@ inline void init()
     polyscope::view::upDir = polyscope::view::UpDir::ZUp;
 
     polyscope::state::userCallback = &mainCallback;
-    //    polyscope::options::autocenterStructures = true;
+    polyscope::options::autocenterStructures = false;
+    polyscope::options::autoscaleStructures = false;
 }
 using PolyscopeLabel = std::string;
 inline std::pair<PolyscopeLabel, size_t> getSelection()
@@ -275,9 +394,12 @@ inline void registerWorldAxes()
 
     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(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);
@@ -347,17 +469,4 @@ inline size_t computeHashOrdered(const std::vector<int> &v)
     return std::hash<std::string>{}(elementsString);
 }
 
-//inline std::filesystem::path getFilepathWithExtension(const std::filesystem::path &folderPath,
-//                                                      const std::string &extension)
-//{
-//    for (const std::filesystem::directory_entry &dirEntry :
-//         std::filesystem::directory_iterator(folderPath)) {
-//        if (dirEntry.is_regular_file() && std::filesystem::path(dirEntry).extension() == extension) {
-//            return std::filesystem::path(dirEntry);
-//        }
-//    }
-
-//    return "";
-//}
-
 #endif // UTILITIES_H

vcgtrimesh.cpp

@@ -7,8 +7,8 @@
 #include <wrap/io_trimesh/import.h>
 //#include <wrap/nanoply/include/nanoplyWrapper.hpp>
 
-bool VCGTriMesh::load(const std::string &filename) {
-    assert(std::filesystem::exists(filename));
+bool VCGTriMesh::load(const std::filesystem::__cxx11::path &meshFilePath) {
+    assert(std::filesystem::exists(meshFilePath));
     unsigned int mask = 0;
     mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
     mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
@@ -20,7 +20,7 @@ bool VCGTriMesh::load(const std::string &filename) {
     //            std::filesystem::absolute(filename).string().c_str(), *this, mask)
     //        != 0) {
     if (vcg::tri::io::Importer<VCGTriMesh>::Open(*this,
-                                                 std::filesystem::absolute(filename).string().c_str())
+                                                 std::filesystem::absolute(meshFilePath).string().c_str())
         != 0) {
         std::cout << "Could not load tri mesh" << std::endl;
         return false;
@@ -32,7 +32,7 @@ bool VCGTriMesh::load(const std::string &filename) {
     vcg::tri::UpdateTopology<VCGTriMesh>::VertexEdge(*this);
     vcg::tri::UpdateNormal<VCGTriMesh>::PerVertexNormalized(*this);
 
-    label = std::filesystem::path(filename).stem().string();
+    label = std::filesystem::path(meshFilePath).stem().string();
     return true;
 }
 

vcgtrimesh.hpp

@@ -39,7 +39,7 @@ class VCGTriMesh : public vcg::tri::TriMesh<std::vector<VCGTriMeshVertex>,
 public:
   VCGTriMesh();
   VCGTriMesh(const std::string &filename);
-  bool load(const std::string &filename) override;
+  bool load(const std::filesystem::path &meshFilePath) override;
   Eigen::MatrixX3d getVertices() const;
   Eigen::MatrixX3i getFaces() const;
   bool save(const std::string plyFilename);