Refactoring

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)
@@ -26,16 +29,18 @@ if(${USE_POLYSCOPE})
         )
         add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR})
         add_compile_definitions(POLYSCOPE_DEFINED)
+        target_sources(${PROJECT_NAME} PUBLIC ${POLYSCOPE_SOURCE_DIR}/deps/imgui/imgui/misc/cpp/imgui_stdlib.cpp)
 endif()
 
 ##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)
@@ -77,23 +82,26 @@ 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)
-target_include_directories(${PROJECT_NAME}
-    PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph
-    PUBLIC ${vcglib_devel_SOURCE_DIR}
-    PUBLIC ${threed-beam-fea_SOURCE_DIR}
-    )
+
 if(${MYSOURCES_STATIC_LINK})
         message("Linking statically")
-        target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread gfortran quadmath)
+        target_link_libraries(${PROJECT_NAME} PUBLIC -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a #[[tbb_static]] pthread gfortran quadmath)
 else()
-        target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot ThreedBeamFEA tbb pthread)
+#    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})
-                target_link_libraries(${PROJECT_NAME} 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 ThreedBeamFEA
+        PUBLIC matplot
+    )
 
 if(USE_ENSMALLEN)
         ##ENSMALLEN
@@ -113,5 +121,5 @@ if(USE_ENSMALLEN)
         else()
                 target_link_libraries(${PROJECT_NAME} PUBLIC  armadillo ensmallen)
         endif()
-        target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include)
+        target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include ${ENSMALLEN_SOURCE_DIR})
 endif()

drmsimulationmodel.cpp

@@ -1889,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();
 }
 

drmsimulationmodel.hpp

@@ -152,7 +152,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,

edgemesh.cpp

@@ -25,7 +25,7 @@ Eigen::MatrixX3d VCGEdgeMesh::getEigenEdgeNormals() const {
   return eigenEdgeNormals;
 }
 
-bool VCGEdgeMesh::save(const string &plyFilename)
+bool VCGEdgeMesh::save(const std::string &plyFilename)
 {
     std::string filename = plyFilename;
     if (filename.empty()) {
@@ -197,36 +197,36 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth,
   return true;
 }
 
-bool VCGEdgeMesh::load(const string &plyFilename) {
+bool VCGEdgeMesh::load(const std::string &plyFilename)
+{
+    std::string usedPath = plyFilename;
+    if (std::filesystem::path(plyFilename).is_relative()) {
+        usedPath = std::filesystem::absolute(plyFilename).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));
-  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(plyFilename).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 << plyFilename << " 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(plyFilename).stem().string();
+    return true;
 }
 
 //bool VCGEdgeMesh::loadUsingNanoply(const std::string &plyFilename) {

linearsimulationmodel.cpp

@@ -108,8 +108,6 @@ SimulationResults LinearSimulationModel::getResults(
         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");
 
@@ -214,7 +212,11 @@ SimulationResults LinearSimulationModel::executeSimulation(
     const std::shared_ptr<SimulationJob> &pSimulationJob)
 {
     assert(pSimulationJob->pMesh->VN() != 0);
-    if (!simulator.wasInitialized()) {
+    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);

reducedmodeloptimizer_structs.hpp

@@ -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,120 +60,321 @@ 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;
+    }
 };
 
-enum OptimizationParameterIndex { E, A, I2, I3, J, R, Theta, NumberOfOptimizationParameters };
+enum OptimizationParameterIndex { E, A, I2, I3, J, R, Theta, NumberOfOptimizationVariables };
+
+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;
+    }
+}
 
 struct Settings
 {
-    std::filesystem::path intermediateResultsDirectoryPath;
-    std::vector<std::vector<OptimizationParameterIndex>> optimizationVariables;
+    inline static std::string defaultFilename{"OptimizationSettings.json"};
+    std::vector<std::vector<OptimizationParameterIndex>> optimizationStrategy;
+    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 };
-    std::vector<xRange> parameterRanges;
-    inline static vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
-    int numberOfFunctionCalls{100000};
-    double solverAccuracy{1e-2};
+    inline static std::vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
     NormalizationStrategy normalizationStrategy{Epsilon};
-    double translationNormalizationParameter{0.003};
-    double rotationNormalizationParameter{3};
+    std::array<xRange, NumberOfOptimizationVariables> variablesRanges;
+    int numberOfFunctionCalls{100000};
+    double solverAccuracy{1e-3};
+    double translationNormalizationEpsilon{1e-4};
+    double rotationNormalizationEpsilon{vcg::math::ToRad(3.0)};
     struct ObjectiveWeights
     {
         double translational{1};
         double rotational{1};
-    } objectiveWeights;
+        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;
+    }
 
     struct JsonKeys
     {
-        inline static std::string filename{"OptimizationSettings.json"};
+        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 TranslationalObjectiveWeight{"TransObjWeight"};
-        inline static std::string OptimizationParameters{"OptimizationParameters"};
+        inline static std::string ObjectiveWeights{"ObjectiveWeight"};
     };
 
+    void setDefault()
+    {
+        ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000};
+        ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000};
+        ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000};
+        ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000};
+        ReducedPatternOptimization::xRange beamJ{"J", 0.001, 1000};
+        ReducedPatternOptimization::xRange innerHexagonSize{"R", 0.05, 0.95};
+        ReducedPatternOptimization::xRange innerHexagonAngle{"Theta", -30.0, 30.0};
+        variablesRanges = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle};
+        numberOfFunctionCalls = 100000;
+
+        enum OptimizationParameterComparisonScenarioIndex {
+            AllVar,
+            GeoYM,
+            MatGeo,
+            YMMat_Geo,
+            YM_MatGeo,
+            MatGeo_YM,
+            Geo_YM_Mat,
+            YM_Geo_Mat,
+            Geo_Mat,
+            YMGeo_Mat,
+            NumberOfScenarios
+        };
+        const std::vector<
+            std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>>
+            optimizationParameters = [&]() {
+                std::vector<
+                    std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>>
+                    optimizationParameters(NumberOfScenarios);
+                using namespace ReducedPatternOptimization;
+                optimizationParameters[AllVar] = {{E, A, I2, I3, J, R, Theta}};
+                optimizationParameters[GeoYM] = {{R, Theta, E}};
+                optimizationParameters[MatGeo] = {{A, I2, I3, J, R, Theta}};
+                optimizationParameters[YMMat_Geo] = {{E, A, I2, I3, J}, {R, Theta}};
+                optimizationParameters[YM_MatGeo] = {{E}, {A, I2, I3, J, R, Theta}};
+                optimizationParameters[MatGeo_YM] = {{A, I2, I3, J, R, Theta}, {E}};
+                optimizationParameters[Geo_YM_Mat] = {{R, Theta}, {E}, {A, I2, I3, J}};
+                optimizationParameters[YM_Geo_Mat] = {{E}, {R, Theta}, {A, I2, I3, J}};
+                optimizationParameters[Geo_Mat] = {{R, Theta}, {A, I2, I3, J}};
+                optimizationParameters[YMGeo_Mat] = {{E, R, Theta}, {A, I2, I3, J}};
+                return optimizationParameters;
+            }();
+
+        constexpr OptimizationParameterComparisonScenarioIndex scenario = MatGeo;
+        optimizationStrategy = optimizationParameters[scenario];
+        if (scenario == YMGeo_Mat) {
+            optimizationVariablesGroupsWeights = {0.15, 0.85};
+        }
+
+        normalizationStrategy = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
+        //#ifdef POLYSCOPE_DEFINED
+        //#else
+        translationNormalizationEpsilon = 4e-3;
+        rotationNormalizationEpsilon = vcg::math::ToRad(3.0);
+
+        //     translationNormalizationParameter = 0;
+        //     rotationNormalizationParameter = 0;
+        //     translationNormalizationParameter = 1e-3;
+        //     rotationNormalizationParameter = vcg::math::ToRad(3.0);
+        //     translationNormalizationParameter = 5e-2;
+        //     rotationNormalizationParameter = vcg::math::ToRad(3.0);
+        //#endif
+        solverAccuracy = 1e-2;
+        //TODO refactor that
+        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].translational = 1.5;
+        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].rotational
+            = 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].translational;
+        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].translational = 1.5;
+        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].rotational
+            = 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].translational;
+        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].translational = 1.2;
+        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].rotational
+            = 2
+              - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].translational;
+        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].translational = 1.95;
+        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].rotational
+            = 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].translational;
+        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].translational = 1.2;
+        perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].rotational
+            = 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].translational;
+    }
+
     void save(const std::filesystem::path &saveToPath)
     {
         assert(std::filesystem::is_directory(saveToPath));
 
         nlohmann::json json;
-        //write x ranges
-        for (size_t xRangeIndex = 0; xRangeIndex < parameterRanges.size(); xRangeIndex++) {
-            const auto &xRange = parameterRanges[xRangeIndex];
-            json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)]
-                = xRange.toString();
+        json[JsonKeys::OptimizationStrategy] = optimizationStrategy;
+        if (optimizationVariablesGroupsWeights.has_value()) {
+            json[JsonKeys::OptimizationStrategyGroupWeights] = 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();
+        //        }
 
         json[JsonKeys::NumberOfFunctionCalls] = numberOfFunctionCalls;
         json[JsonKeys::SolverAccuracy] = solverAccuracy;
-        json[JsonKeys::TranslationalObjectiveWeight] = objectiveWeights.translational;
-        json[JsonKeys::OptimizationParameters] = optimizationVariables;
+        //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;
 
         std::filesystem::path jsonFilePath(
-            std::filesystem::path(saveToPath).append(JsonKeys::filename));
+            std::filesystem::path(saveToPath).append(defaultFilename));
         std::ofstream jsonFile(jsonFilePath.string());
         jsonFile << json;
     }
 
-    bool load(const std::filesystem::path &loadFromPath)
+    bool load(const std::filesystem::path &jsonFilepath)
     {
-        assert(std::filesystem::is_directory(loadFromPath));
-        //Load optimal X
         nlohmann::json json;
-        std::filesystem::path jsonFilepath(
-            std::filesystem::path(loadFromPath).append(JsonKeys::filename));
         if (!std::filesystem::exists(jsonFilepath)) {
+            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;
 
+        if (json.contains(JsonKeys::OptimizationStrategy)) {
+            optimizationStrategy
+                = std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>(
+                    (json.at(JsonKeys::OptimizationStrategy)));
+        }
+        if (json.contains(JsonKeys::OptimizationStrategyGroupWeights)) {
+            optimizationVariablesGroupsWeights = std::vector<double>(
+                json[JsonKeys::OptimizationStrategyGroupWeights]);
+        }
         //read x ranges
-        size_t xRangeIndex = 0;
-        while (true) {
-            const std::string jsonXRangeKey = JsonKeys::OptimizationVariables + "_"
-                                              + std::to_string(xRangeIndex);
-            if (!json.contains(jsonXRangeKey)) {
-                break;
+        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;
             }
-            xRange x;
-            x.fromString(json.at(jsonXRangeKey));
-            parameterRanges.push_back(x);
-            xRangeIndex++;
         }
 
-        optimizationVariables
-            = std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>(
-                (json.at(JsonKeys::OptimizationParameters)));
         numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls);
         solverAccuracy = json.at(JsonKeys::SolverAccuracy);
-        objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight);
-        objectiveWeights.rotational = 2 - objectiveWeights.translational;
+        //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};
+                           });
+        }
+
+        //        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 (!parameterRanges.empty()) {
+        if (!variablesRanges.empty()) {
             std::string xRangesString;
-            for (const xRange &range : parameterRanges) {
+            for (const xRange &range : variablesRanges) {
                 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);
+                          + " 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 (!parameterRanges.empty()) {
-            for (const xRange &range : parameterRanges) {
+        if (!variablesRanges.empty()) {
+            for (const xRange &range : variablesRanges) {
                 os << range.label + " max";
                 os << range.label + " min";
             }
@@ -181,16 +382,15 @@ struct Settings
         os << "Function Calls";
         os << "Solution Accuracy";
         os << "Normalization strategy";
-        os << "Trans weight";
-        os << "Rot weight";
+        os << JsonKeys::ObjectiveWeights;
         os << "Optimization parameters";
         //        os << std::endl;
     }
 
     void writeSettingsTo(csvFile &os) const
     {
-        if (!parameterRanges.empty()) {
-            for (const xRange &range : parameterRanges) {
+        if (!variablesRanges.empty()) {
+            for (const xRange &range : variablesRanges) {
                 os << range.max;
                 os << range.min;
             }
@@ -198,13 +398,22 @@ struct Settings
         os << numberOfFunctionCalls;
         os << solverAccuracy;
         os << normalizationStrategyStrings[normalizationStrategy] + "_"
-                  + std::to_string(translationNormalizationParameter);
-        os << objectiveWeights.translational;
-        os << objectiveWeights.rotational;
+                  + std::to_string(translationNormalizationEpsilon);
+        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 : optimizationVariables) {
+        for (const std::vector<OptimizationParameterIndex> &v : optimizationStrategy) {
             std::vector<int> vi;
             vi.reserve(v.size());
             for (const OptimizationParameterIndex &parameter : v) {
@@ -218,12 +427,17 @@ struct Settings
 
   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.parameterRanges == settings2.parameterRanges
-             && settings1.solverAccuracy == settings2.solverAccuracy
-             && settings1.objectiveWeights.translational == settings2.objectiveWeights.translational
-             && settings1.translationNormalizationParameter
-                    == settings2.translationNormalizationParameter;
+             && settings1.variablesRanges == settings2.variablesRanges
+             && settings1.solverAccuracy == settings2.solverAccuracy && haveTheSameObjectiveWeights
+             && settings1.translationNormalizationEpsilon
+                    == settings2.translationNormalizationEpsilon;
   }
 
   inline void updateMeshWithOptimalVariables(const std::vector<double> &x, SimulationMesh &m)
@@ -306,6 +520,7 @@ struct Settings
           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;
@@ -332,7 +547,7 @@ struct Settings
           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)) {
@@ -405,7 +620,8 @@ struct Settings
                       = 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");
@@ -435,7 +651,7 @@ struct Settings
 #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(
@@ -496,40 +712,36 @@ struct Settings
           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(
-                           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;
       }
@@ -543,9 +755,16 @@ struct Settings
           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"));
+          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);
       }
@@ -589,7 +808,7 @@ struct Settings
 
       static void applyOptimizationResults_elements(
           const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
-          const shared_ptr<SimulationMesh> &pTiledReducedPattern_simulationMesh)
+          const std::shared_ptr<SimulationMesh> &pTiledReducedPattern_simulationMesh)
       {
           assert(CrossSectionType::name == RectangularBeamDimensions::name);
           // Set reduced parameters fron the optimization results

simulation_structs.hpp

@@ -133,32 +133,32 @@ public:
       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;
@@ -377,7 +377,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);
@@ -393,7 +393,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;
@@ -712,21 +713,21 @@ struct SimulationResults
         //          /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
         //      //}
         pJob->registerForDrawing(getLabel());
-        static bool wasExecuted = false;
-        if (!wasExecuted) {
-            std::function<void()> callback = [&]() {
-                static bool showFrames = shouldShowFrames;
+        //        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);
-                }
-            };
+        //                if (ImGui::Checkbox("Show Frames", &showFrames) && showFrames) {
+        //                    polyscopeHandle_frameX->setEnabled(showFrames);
+        //                    polyscopeHandle_frameY->setEnabled(showFrames);
+        //                    polyscopeHandle_frameZ->setEnabled(showFrames);
+        //                }
+        //            };
 
-            PolyscopeInterface::addUserCallback(callback);
-            wasExecuted = true;
-        }
+        //            PolyscopeInterface::addUserCallback(callback);
+        //            wasExecuted = true;
+        //        }
         return polyscopeHandle_deformedEdmeMesh;
     }
 #endif
@@ -738,7 +739,7 @@ private:
         //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);

simulationhistoryplotter.hpp

@@ -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()) {

simulationmesh.cpp

@@ -284,85 +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 (!VCGEdgeMesh::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();
-  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;
+    //  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)

simulationmesh.hpp

@@ -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();
 

trianglepatterngeometry.cpp

@@ -145,7 +145,7 @@ PatternGeometry::PatternGeometry(PatternGeometry &other) {
   vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
 }
 
-bool PatternGeometry::load(const string &plyFilename)
+bool PatternGeometry::load(const std::string &plyFilename)
 {
     if (!VCGEdgeMesh::load(plyFilename)) {
         return false;

utilities.hpp

@@ -4,9 +4,11 @@
 #include <Eigen/Dense>
 #include <algorithm>
 #include <array>
+#include <chrono>
 #include <filesystem>
 #include <fstream>
 #include <iterator>
+#include <numeric>
 #include <regex>
 
 #define GET_VARIABLE_NAME(Variable) (#Variable)
@@ -140,16 +142,37 @@ struct Vector6d : public std::array<double, 6> {
 };
 
 namespace Utilities {
-//inline std::vector<std::string> split(std::string text, char delim)
-//{
-//    std::string line;
-//    std::vector<std::string> vec;
-//    std::stringstream ss(text);
-//    while (std::getline(ss, line, delim)) {
-//        vec.push_back(line);
-//    }
-//    return vec;
-//}
+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]));
+    }
+
+    // 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;
+
+    // Numbers are the same --> remove numbers and recurse
+    std::string anew, bnew;
+    std::getline(issa, anew);
+    std::getline(issb, bnew);
+    return (compareNat(anew, bnew));
+}
 
 inline std::string_view leftTrimSpaces(const std::string_view& str)
 {
@@ -175,6 +198,21 @@ std::string_view trimmedString=str;
     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;
@@ -273,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
@@ -293,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();
     }
 
@@ -418,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