Refactoring. Removed the settings of CMAKE_SYSTEM to windows

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.8)
 project(ReducedModelOptimization)
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
-        message(STATUS "The compiler ${CMAKE_CXX_COMPILER}")
+message(STATUS "The compiler ${CMAKE_CXX_COMPILER}")
 
 #Add the project cmake scripts to the module path
 list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
@@ -26,23 +26,12 @@ endif()
 ##Create directory for the external libraries
 file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR})
 
-##Polyscope
 if(${CMAKE_BUILD_TYPE} STREQUAL "Release")
         set(USE_POLYSCOPE FALSE)
 else()
         set(USE_POLYSCOPE TRUE)
 endif()
-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}
-)
-add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR})
-add_compile_definitions(POLYSCOPE_DEFINED)
-endif()
-
+set(MYSOURCES_STATIC_LINK FALSE)
 
 set(MYSOURCES_SOURCE_DIR "/home/iason/Coding/Libraries/MySources")
 if (EXISTS ${MYSOURCES_SOURCE_DIR})
@@ -57,6 +46,20 @@ download_project(PROJ               MYSOURCES
 endif()
 add_subdirectory(${MYSOURCES_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}/MySourcesBinDir)
 
+##Polyscope
+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})
+endif()
+add_compile_definitions(POLYSCOPE_DEFINED)
+endif()
+
 #dlib
 set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin)
 file(MAKE_DIRECTORY ${DLIB_BIN_DIR})
@@ -101,8 +104,13 @@ target_include_directories(${PROJECT_NAME}
 
 
 target_link_directories(${PROJECT_NAME} PRIVATE ${MYSOURCES_SOURCE_DIR}/boost_graph/libs/)
-if(${USE_POLYSCOPE})
-        target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen  dlib::dlib MySources)
-else()
+if(${MYSOURCES_STATIC_LINK})
+        message("Linking statically")
         target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen dlib::dlib MySources)
+else()
+        if(${USE_POLYSCOPE})
+                message("Using polyscope")
+                target_link_libraries(${PROJECT_NAME} polyscope)
+        endif()
+        target_link_libraries(${PROJECT_NAME} Eigen3::Eigen dlib::dlib MySources)
 endif()

src/main.cpp

@@ -40,33 +40,45 @@ int main(int argc, char *argv[]) {
   // Set the optization settings
   ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000};
   ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000};
+  ReducedPatternOptimization::xRange beamI{"I", 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{"HexSize", 0.05, 0.95};
   ReducedPatternOptimization::xRange innerHexagonAngle{"HexAngle", -30.0, 30.0};
   ReducedPatternOptimization::Settings settings_optimization;
-  settings_optimization.xRanges = {beamE,beamA,beamJ,beamI2,beamI3,
-                                   innerHexagonSize, innerHexagonAngle};
+  settings_optimization.xRanges
+      = {/*beamE,*/ beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle};
   const bool input_numberOfFunctionCallsDefined = argc >= 4;
-  settings_optimization.numberOfFunctionCalls =
-      input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100;
+  settings_optimization.numberOfFunctionCalls = input_numberOfFunctionCallsDefined
+                                                    ? std::atoi(argv[3])
+                                                    : 100;
   settings_optimization.normalizationStrategy
       = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
-
-  settings_optimization.normalizationParameter = 0.0003;
-  settings_optimization.solverAccuracy = 0.001;
+  settings_optimization.splitGeometryMaterialOptimization = false;
+  settings_optimization.translationNormalizationParameter = 1e-3;
+  settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0);
+  //  settings_optimization.translationNormalizationParameter = 1e-15;
+  //  settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(1e-15);
+  //  settings_optimization.solverAccuracy = 1e-3;
+  settings_optimization.solverAccuracy = 1e-1;
   settings_optimization.objectiveWeights.translational = std::atof(argv[4]);
   settings_optimization.objectiveWeights.rotational = 2 - std::atof(argv[4]);
 
-  // Optimize pair
+  std::string xConcatNames;
+  for (const auto &x : settings_optimization.xRanges) {
+      xConcatNames.append(x.label + "_");
+  }
+  xConcatNames.pop_back();
+
+  // Optimize pairthere
   const std::string pairName = fullPattern.getLabel() + "@" + reducedPattern.getLabel();
   const std::string optimizationName = pairName + "("
                                        + std::to_string(settings_optimization.numberOfFunctionCalls)
                                        + "_"
                                        + to_string_with_precision(
                                            settings_optimization.objectiveWeights.translational)
-                                       + ")";
+                                       + ")" + "_" + xConcatNames;
   const bool input_resultDirectoryDefined = argc >= 6;
   const std::string optimizationResultsDirectory = input_resultDirectoryDefined
                                                        ? argv[5]
@@ -91,22 +103,21 @@ int main(int argc, char *argv[]) {
   }
 
   ReducedPatternOptimization::Results optimizationResults;
-  bool optimizationAlreadyComputed = optimizationResultFolderExists;
-  //  bool optimizationAlreadyComputed = false;
-  //  if (optimizationResultFolderExists) {
-  //      const bool resultsWereSuccessfullyLoaded = optimizationResults.load(resultsOutputDir);
-  //      if (resultsWereSuccessfullyLoaded && optimizationResults.settings == settings_optimization) {
-  //      }
-  //  }
+  constexpr bool shouldReoptimize = true;
+  bool optimizationAlreadyComputed = false;
+  if (!shouldReoptimize && optimizationResultFolderExists) {
+      const bool resultsWereSuccessfullyLoaded = optimizationResults.load(resultsOutputDir);
+      if (resultsWereSuccessfullyLoaded && optimizationResults.settings == settings_optimization) {
+          optimizationAlreadyComputed = true;
+      }
+  }
 
   if (!optimizationAlreadyComputed) {
       auto start = std::chrono::system_clock::now();
       const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
       assert(interfaceNodeIndex == numberOfNodesPerSlot[0] + numberOfNodesPerSlot[3]);
       ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
-      optimizer.initializePatterns(fullPattern,
-                                   reducedPattern,
-                                   settings_optimization.xRanges.size());
+      optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges);
       optimizer.optimize(settings_optimization, optimizationResults);
       optimizationResults.label = optimizationName;
       optimizationResults.baseTriangleFullPattern.copy(fullPattern);
@@ -118,17 +129,17 @@ int main(int argc, char *argv[]) {
       if (optimizationResults.wasSuccessful) {
           resultsOutputDir = convergedJobsDirPath.string();
           csvFile csv_results({}, false);
-          //          csvFile csv_results(std::filesystem::path(resultsOutputDir).append("optimizationDistances.csv"), false);
           csv_results << "Name";
           optimizationResults.writeHeaderTo(csv_results);
           settings_optimization.writeHeaderTo(csv_results);
           csv_results << endrow;
           csv_results << std::to_string(fullPattern.EN()) + "#" + pairName;
-          optimizationResults.writeResultsTo(settings_optimization, csv_results);
+          optimizationResults.writeResultsTo(csv_results);
           settings_optimization.writeSettingsTo(csv_results);
           csv_results << endrow;
       } else {
           resultsOutputDir = crashedJobsDirPath.string();
+          return 1;
       }
       optimizationResults.save(resultsOutputDir, true);
   }

src/reducedmodeloptimizer.cpp

@@ -20,7 +20,9 @@ struct GlobalOptimizationVariables {
     std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
         fullPatternInterfaceViPairs;
     matplot::line_handle gPlotHandle;
+    std::vector<size_t> objectiveValueHistory_iteration;
     std::vector<double> objectiveValueHistory;
+    std::vector<double> plotColors;
     Eigen::VectorXd initialParameters;
     std::vector<int> simulationScenarioIndices;
     double minY{DBL_MAX};
@@ -40,6 +42,10 @@ struct GlobalOptimizationVariables {
     double desiredMaxDisplacementValue;
     double desiredMaxRotationAngle;
     std::string currentScenarioName;
+    std::vector<std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                                   const double &newValue)>>
+        updateReducedPatternFunctions_material;
+
 } global;
 
 double ReducedModelOptimizer::computeDisplacementError(
@@ -52,6 +58,8 @@ double ReducedModelOptimizer::computeDisplacementError(
     const double rawError = computeRawTranslationalError(fullPatternDisplacements,
                                                          reducedPatternDisplacements,
                                                          reducedToFullInterfaceViMap);
+    //    std::cout << "raw trans error:" << rawError << std::endl;
+    //    std::cout << "raw trans error:" << normalizationFactor << std::endl;
     return rawError / normalizationFactor;
 }
 
@@ -120,6 +128,9 @@ double ReducedModelOptimizer::computeError(
                                    simulationResults_reducedPattern.displacements,
                                    reducedToFullInterfaceViMap,
                                    normalizationFactor_translationalDisplacement);
+
+    //    std::cout << "normalization factor:" << normalizationFactor_rotationalDisplacement << std::endl;
+    //    std::cout << "trans error:" << translationalError << std::endl;
     const double rotationalError
         = computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion,
                                  simulationResults_reducedPattern.rotationalDisplacementQuaternion,
@@ -129,18 +140,19 @@ double ReducedModelOptimizer::computeError(
            + global.optimizationSettings.objectiveWeights.rotational * rotationalError;
 }
 
-double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3,
-                                        double innerHexagonSize,
-                                        double innerHexagonRotationAngle) {
-    std::vector<double> x{E,A,J,I2,I3, innerHexagonSize, innerHexagonRotationAngle};
-    return ReducedModelOptimizer::objective(x.size(), x.data());
-}
+//double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3,
+//                                        double innerHexagonSize,
+//                                        double innerHexagonRotationAngle) {
+//    std::vector<double> x{E,A,J,I2,I3, innerHexagonSize, innerHexagonRotationAngle};
+//    return ReducedModelOptimizer::objective(x.size(), x.data());
+//}
 
-double ReducedModelOptimizer::objective(long n, const double *x) {
-    //  std::cout.precision(17);
-    //        for (int i = 0; i < n; i++) {
-    //    std::cout << x[i] << " ";
-    //        }
+double ReducedModelOptimizer::objective(const dlib::matrix<double, 0, 1> &x)
+{
+    //    std::cout.precision(17);
+    //    for (int i = 0; i < x.size(); i++) {
+    //        std::cout << x(i) << " ";
+    //    }
     //    std::cout << std::endl;
 
     //    std::cout << x[n - 2] << " " << x[n - 1] << std::endl;
@@ -150,13 +162,16 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
     //            << e.firstBendingConstFactor << " " <<
     //            e.secondBendingConstFactor
     //            << std::endl;
-    updateMesh(n, x);
+    const int n = x.size();
+    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh
+        = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh;
+    function_updateReducedPattern(x, pReducedPatternSimulationMesh);
     //    global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing();
     //    global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing();
     //    polyscope::show();
     //    global.reducedPatternSimulationJobs[0]->pMesh->unregister();
 
-    //  std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
+    //  std::cout << e.axialConstFact|A,I2,I3,J,r,thetaor << " " << e.torsionConstFactor << " "
     //            << e.firstBendingConstFactor << " " <<
     //            e.secondBendingConstFactor
     //            << std::endl;
@@ -168,6 +183,9 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
     for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
         const std::shared_ptr<SimulationJob> &reducedJob
             = global.reducedPatternSimulationJobs[simulationScenarioIndex];
+        //#ifdef POLYSCOPE_DEFINED
+        //        std::cout << reducedJob->getLabel() << ":" << std::endl;
+        //#endif
         SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob);
         //    std::string filename;
         if (!reducedModelResults.converged) {
@@ -177,7 +195,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
             std::cout << "Failed simulation" << std::endl;
 #endif
         } else {
-            const bool usePotentialEnergy = false;
+            constexpr bool usePotentialEnergy = false;
             double simulationScenarioError;
             if (usePotentialEnergy) {
                 simulationScenarioError = std::abs(
@@ -191,17 +209,17 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
                     global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
                     global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
             }
-#ifdef POLYSCOPE_DEFINED
-//            std::cout << reducedJob->getLabel() << " sim error:" << simulationScenarioError
-//                      << std::endl;
-//                        reducedJob->pMesh->registerForDrawing(Colors::reducedInitial);
-//                        reducedModelResults.registerForDrawing(Colors::reducedDeformed);
-//                        global.fullPatternResults[simulationScenarioIndex].registerForDrawing(
-//                            Colors::fullDeformed);
-//                        polyscope::show();
-//                        reducedModelResults.unregister();
-//                        global.fullPatternResults[simulationScenarioIndex].unregister();
-#endif
+            //#ifdef POLYSCOPE_DEFINED
+            //            reducedJob->pMesh->registerForDrawing(Colors::reducedInitial);
+            //            reducedModelResults.registerForDrawing(Colors::reducedDeformed);
+            //            global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed);
+            //            global.fullPatternResults[simulationScenarioIndex].registerForDrawing(
+            //                Colors::fullDeformed);
+            //            polyscope::show();
+            //            reducedModelResults.unregister();
+            //            global.pFullPatternSimulationMesh->unregister();
+            //            global.fullPatternResults[simulationScenarioIndex].unregister();
+            //#endif
             //        if (global.optimizationSettings.normalizationStrategy !=
             //        NormalizationStrategy::Epsilon &&
             //            simulationScenarioError > 1) {
@@ -231,35 +249,47 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
             //                        global.reducedToFullInterfaceViMap, true);
             //                    polyscope::removeAllStructures();
             //            #endif // POLYSCOPE_DEFINED
-            totalError += simulationScenarioError;
+            totalError += simulationScenarioError * simulationScenarioError;
         }
     }
     //  std::cout << error << std::endl;
+//    global.objectiveValueHistory.push_back(totalError);
+//    global.plotColors.push_back(10);
+    ++global.numberOfFunctionCalls;
     if (totalError < global.minY) {
         global.minY = totalError;
-        global.minX.assign(x, x + n);
+        global.objectiveValueHistory.push_back(totalError);
+        global.objectiveValueHistory_iteration.push_back(global.numberOfFunctionCalls);
+        //        std::cout << "New best:" << totalError << std::endl;
+        //        //        global.minX.assign(x.begin(), x.begin() + n);
+        //        std::cout.precision(17);
+        //        for (int i = 0; i < x.size(); i++) {
+        //            std::cout << x(i) << " ";
+        //        }
+        //        std::cout << std::endl;
+        //        global.objectiveValueHistoryY.push_back(std::log(totalError));
+        //        global.objectiveValueHistoryX.push_back(global.numberOfFunctionCalls);
+        //        global.plotColors.push_back(0.1);
+        //        auto xPlot = matplot::linspace(0,
+        //                                       global.objectiveValueHistoryY.size(),
+        //                                       global.objectiveValueHistoryY.size());
+        //        global.gPlotHandle = matplot::scatter(global.objectiveValueHistoryX,
+        //                                              global.objectiveValueHistoryY,
+        //                                              4,
+        //                                              global.plotColors);
+        //        matplot::show();
+        //                SimulationResultsReporter::createPlot("Number of Steps",
+        //                                              "Objective value",
+        //                                              global.objectiveValueHistoryY);
     }
 #ifdef POLYSCOPE_DEFINED
-    ++global.numberOfFunctionCalls;
     if (global.optimizationSettings.numberOfFunctionCalls >= 100
         && global.numberOfFunctionCalls % (global.optimizationSettings.numberOfFunctionCalls / 100)
                == 0) {
         std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl;
     }
 #endif
-//     compute error and return it
-//      global.objectiveValueHistory.push_back(totalError);
-//      auto xPlot = matplot::linspace(0, global.objectiveValueHistory.size(),
-//                                     global.objectiveValueHistory.size());
-//      std::vector<double> colors(global.gObjectiveValueHistory.size(), 2);
-//      if (global.g_firstRoundIterationIndex != 0) {
-//        for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(),
-//                 [](double &c) { c = 0.7; });
-//      }
-//      global.gPlotHandle = matplot::scatter(xPlot, global.objectiveValueHistory);
-//      SimulationResultsReporter::createPlot("Number of Steps", "Objective value",
-//                                            global.objectiveValueHistory);
-
+    //     compute error and return it
     return totalError;
 }
 
@@ -271,18 +301,19 @@ void ReducedModelOptimizer::createSimulationMeshes(
         std::cerr << "Error: A rectangular cross section is expected." << std::endl;
         terminate();
     }
+    const double ym = 1 * 1e9;
     // Full pattern
     pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
     pFullPatternSimulationMesh->setBeamCrossSection(
         CrossSectionType{0.002, 0.002});
-    pFullPatternSimulationMesh->setBeamMaterial(0.3, 2.3465 * 1e9);
+    pFullPatternSimulationMesh->setBeamMaterial(0.3, ym);
 
     // Reduced pattern
     pReducedPatternSimulationMesh =
         std::make_shared<SimulationMesh>(reducedModel);
     pReducedPatternSimulationMesh->setBeamCrossSection(
         CrossSectionType{0.002, 0.002});
-    pReducedPatternSimulationMesh->setBeamMaterial(0.3, 2.3465 * 1e9);
+    pReducedPatternSimulationMesh->setBeamMaterial(0.3, ym);
 }
 
 void ReducedModelOptimizer::createSimulationMeshes(
@@ -461,7 +492,7 @@ ReducedModelOptimizer::ReducedModelOptimizer(const std::vector<size_t> &numberOf
 
 void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
                                                PatternGeometry &reducedPattern,
-                                               const int &optimizationParameters)
+                                               const std::vector<xRange> &optimizationParameters)
 {
     assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN());
     fullPatternNumberOfEdges = fullPattern.EN();
@@ -477,82 +508,157 @@ void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
 
     computeMaps(copyFullPattern, copyReducedPattern);
     createSimulationMeshes(copyFullPattern, copyReducedPattern);
+    initializeUpdateReducedPatternFunctions();
     initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters);
 }
 
-void updateMesh(long n, const double *x) {
-    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
-        global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
-            ->pMesh;
+void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions()
+{
+    function_updateReducedPattern_geometry = [&](const dlib::matrix<double, 0, 1> &x,
+                                                 std::shared_ptr<SimulationMesh>
+                                                     &pReducedPatternSimulationMesh) {
+        const int n = x.size();
+        assert(n >= 2);
 
-    const double E=global.initialParameters(0)*x[0];
-    const double A=global.initialParameters(1) * x[1];
-    const double beamWidth=std::sqrt(A);
-        const double beamHeight=beamWidth;
+        //    CoordType center_barycentric(1, 0, 0);
+        //    CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5);
+        //    CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric)
+        //                                        * x[n - 2]);
 
-const double J=global.initialParameters(2) * x[2];
-const double I2 = global.initialParameters(3) * x[3];
-const double I3=global.initialParameters(4) * x[4];
-    for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
-        Element &e = pReducedPatternSimulationMesh->elements[ei];
-        e.setDimensions(
-            RectangularBeamDimensions(beamWidth,
-                                      beamHeight));
-        e.setMaterial(ElementMaterial(e.material.poissonsRatio,
-                                      E));
-            e.J = J;
-            e.I2 = I2;
-            e.I3 = I3;
-    }
-    assert(pReducedPatternSimulationMesh->EN() == 12);
-    assert(n >= 2);
-
-    //    CoordType center_barycentric(1, 0, 0);
-    //    CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5);
-    //    CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric)
-    //                                        * x[n - 2]);
-
-    CoordType movableVertex_barycentric(1 - x[n - 2], x[n - 2] / 2, x[n - 2] / 2);
-    CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0)
-                                                      * movableVertex_barycentric[0]
-                                                  + global.baseTriangle.cP(1)
-                                                        * movableVertex_barycentric[1]
-                                                  + global.baseTriangle.cP(2)
-                                                        * movableVertex_barycentric[2];
-    baseTriangleMovableVertexPosition
-        = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, vcg::math::ToRad(x[n - 1]))
-          * baseTriangleMovableVertexPosition;
-
-    for (int rotationCounter = 0;
-         rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {
-        pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
+        CoordType movableVertex_barycentric(1 - x(n - 2), x(n - 2) / 2, x(n - 2) / 2);
+        CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0)
+                                                          * movableVertex_barycentric[0]
+                                                      + global.baseTriangle.cP(1)
+                                                            * movableVertex_barycentric[1]
+                                                      + global.baseTriangle.cP(2)
+                                                            * movableVertex_barycentric[2];
+        baseTriangleMovableVertexPosition
             = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
-                                  vcg::math::ToRad(60.0 * rotationCounter))
+                                  vcg::math::ToRad(x(n - 1)))
               * baseTriangleMovableVertexPosition;
-    }
 
-    pReducedPatternSimulationMesh->reset();
-    //#ifdef POLYSCOPE_DEFINED
-    //    pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
-    //    pReducedPatternSimulationMesh->registerForDrawing();
-    //    std::cout << "Angle:" + std::to_string(x[n - 1]) + " size:" + std::to_string(x[n - 2])
-    //              << std::endl;
-    //    std::cout << "Verts:" << pReducedPatternSimulationMesh->VN() << std::endl;
-    //    polyscope::show();
-    //#endif
+        for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize;
+             rotationCounter++) {
+            pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
+                = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
+                                      vcg::math::ToRad(60.0 * rotationCounter))
+                  * baseTriangleMovableVertexPosition;
+        }
+    };
+
+    function_updateReducedPattern_material_E =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newE) {
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.setMaterial(ElementMaterial(e.material.poissonsRatio, newE));
+            }
+        };
+    function_updateReducedPattern_material_A =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newA) {
+            const double beamWidth = std::sqrt(newA);
+            const double beamHeight = beamWidth;
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight));
+            }
+        };
+    function_updateReducedPattern_material_I =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newI) {
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.inertia.I2 = newI;
+                e.inertia.I3 = newI;
+            }
+        };
+    function_updateReducedPattern_material_I2 =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newI2) {
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.inertia.I2 = newI2;
+            }
+        };
+    function_updateReducedPattern_material_I3 =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newI3) {
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.inertia.I3 = newI3;
+            }
+        };
+    function_updateReducedPattern_material_J =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newJ) {
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.inertia.J = newJ;
+            }
+        };
 }
 
-
 void ReducedModelOptimizer::initializeOptimizationParameters(
-    const std::shared_ptr<SimulationMesh> &mesh,const int& optimizationParamters) {
-    global.numberOfOptimizationParameters = optimizationParamters;
+    const std::shared_ptr<SimulationMesh> &mesh, const std::vector<xRange> &optimizationParameters)
+{
+    global.numberOfOptimizationParameters = optimizationParameters.size();
     global.initialParameters.resize(global.numberOfOptimizationParameters);
 
-      global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
-      global.initialParameters(1) = mesh->elements[0].A;
-      global.initialParameters(2) = mesh->elements[0].J;
-      global.initialParameters(3) = mesh->elements[0].I2;
-      global.initialParameters(4) = mesh->elements[0].I3;
+    for (int optimizationParameterIndex = 0;
+         optimizationParameterIndex < optimizationParameters.size();
+         optimizationParameterIndex++) {
+        const xRange &xOptimizationParameter = optimizationParameters[optimizationParameterIndex];
+        if (xOptimizationParameter.label == "E") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0]
+                                                                       .material.youngsModulus;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_E);
+        } else if (xOptimizationParameter.label == "A") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0].A;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_A);
+        } else if (xOptimizationParameter.label == "I") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I2;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_I);
+        } else if (xOptimizationParameter.label == "I2") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I2;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_I2);
+        } else if (xOptimizationParameter.label == "I3") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I3;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_I3);
+        } else if (xOptimizationParameter.label == "J") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.J;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_J);
+        }
+        //NOTE:Assuming that I2,I3 and J are be passed after E and A because the update of E and A changes I2,I3,J.
+        //        case "HexSize":
+        //            updateReducedPatternFunctions_material.push_back(
+        //                function_updateReducedPattern_material_E);
+        //            break;
+        //        case "HexAngle":
+        //            updateReducedPatternFunctions_material.push_back(
+        //                function_updateReducedPattern_material_E);
+        //            break;
+    }
+
+    if (global.updateReducedPatternFunctions_material.size() != 0) {
+        function_updateReducedPattern_material =
+            [&](const dlib::matrix<double, 0, 1> &x,
+                std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+                for (int optimizationParameterIndex = 0;
+                     optimizationParameterIndex
+                     < global.updateReducedPatternFunctions_material.size();
+                     optimizationParameterIndex++) {
+                    global.updateReducedPatternFunctions_material[optimizationParameterIndex](
+                        pReducedPatternSimulationMesh,
+                        global.initialParameters(optimizationParameterIndex)
+                            * x(optimizationParameterIndex));
+                }
+            };
+    } else {
+        function_updateReducedPattern_material =
+            [](const dlib::matrix<double, 0, 1> &x,
+               std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {};
+    }
 }
 
 void ReducedModelOptimizer::computeReducedModelSimulationJob(
@@ -715,19 +821,18 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
     results.optimalXNameValuePairs.resize(settings.xRanges.size());
         std::vector<double> optimalX(settings.xRanges.size());
     for (int xVariableIndex = 0; xVariableIndex < settings.xRanges.size(); xVariableIndex++) {
-        if (xVariableIndex < 5) {
+        if (xVariableIndex < settings.xRanges.size() - 2) {
             results.optimalXNameValuePairs[xVariableIndex]
                 = std::make_pair(settings.xRanges[xVariableIndex].label,
-                                 global.minX[xVariableIndex]
+                                 optimizationResult_dlib.x(xVariableIndex)
                                      * global.initialParameters(xVariableIndex));
         } else {
             //Hex size and angle are pure values (not multipliers of the initial values)
             results.optimalXNameValuePairs[xVariableIndex]
                 = std::make_pair(settings.xRanges[xVariableIndex].label,
-                                 global.minX[xVariableIndex]);
+                                 optimizationResult_dlib.x(xVariableIndex));
         }
 
-        assert(global.minX[xVariableIndex] == optimizationResult_dlib.x(xVariableIndex));
         optimalX[xVariableIndex] = optimizationResult_dlib.x(xVariableIndex);
 #ifdef POLYSCOPE_DEFINED
         std::cout << results.optimalXNameValuePairs[xVariableIndex].first << ":"
@@ -739,7 +844,10 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
 #endif
 
     // Compute obj value per simulation scenario and the raw objective value
-    updateMesh(optimalX.size(), optimalX.data());
+    //    updateMeshFunction(optimalX);
+    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh
+        = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh;
+    function_updateReducedPattern(optimizationResult_dlib.x, pReducedPatternSimulationMesh);
     //    results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios);
     //TODO:use push_back it will make the code more readable
     LinearSimulationModel simulator;
@@ -754,6 +862,10 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
         global.simulationScenarioIndices.size());
     results.objectiveValue.perSimulationScenario_total.resize(
         global.simulationScenarioIndices.size());
+    //#ifdef POLYSCOPE_DEFINED
+    //    global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed);
+    //#endif
+
     results.perScenario_fullPatternPotentialEnergy.resize(global.simulationScenarioIndices.size());
     for (int i = 0; i < global.simulationScenarioIndices.size(); i++) {
         const int simulationScenarioIndex = global.simulationScenarioIndices[i];
@@ -828,17 +940,14 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
         std::cout << "Total Error value:" << results.objectiveValue.perSimulationScenario_total[i]
                   << std::endl;
         std::cout << std::endl;
-#endif
-    }
 
-    const bool printDebugInfo = false;
-    if (printDebugInfo) {
-        std::cout << "Finished optimizing." << endl;
-    std::cout << "Total optimal objective value:" << results.objectiveValue.total << std::endl;
-    assert(global.minY == optimizationResult_dlib.y);
-    if (global.minY != optimizationResult_dlib.y) {
-        std::cerr << "ERROR in objective value" << std::endl;
-    }
+        //        reducedModelResults.registerForDrawing(Colors::reducedDeformed);
+        //        global.fullPatternResults[simulationScenarioIndex].registerForDrawing(Colors::fullDeformed);
+        //        polyscope::show();
+        //        reducedModelResults.unregister();
+        //        global.fullPatternResults[simulationScenarioIndex].unregister();
+
+#endif
     }
 
     for (int simulationScenarioIndex : global.simulationScenarioIndices) {
@@ -852,6 +961,9 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
         //        global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing();
         //        global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel(temp);
     }
+    results.objectiveValueHistory = global.objectiveValueHistory;
+    results.objectiveValueHistory_iteration = global.objectiveValueHistory_iteration;
+    //    results.draw();
 }
 
 std::vector<std::pair<BaseSimulationScenario, double>>
@@ -890,10 +1002,15 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces(
         fullPatternSimulationScenarioMaxMagnitudes
             = static_cast<std::vector<std::pair<BaseSimulationScenario, double>>>(
                 json.at("maxMagn"));
-    } else {
+        const bool shouldRecompute = fullPatternSimulationScenarioMaxMagnitudes.size()
+                                     != desiredBaseSimulationScenarioIndices.size();
+        if (!shouldRecompute) {
+            return fullPatternSimulationScenarioMaxMagnitudes;
+        }
+    }
 #endif
-        fullPatternSimulationScenarioMaxMagnitudes = computeFullPatternMaxSimulationForces(
-            desiredBaseSimulationScenarioIndices);
+    fullPatternSimulationScenarioMaxMagnitudes = computeFullPatternMaxSimulationForces(
+        desiredBaseSimulationScenarioIndices);
 
 #ifdef POLYSCOPE_DEFINED
         nlohmann::json json;
@@ -902,34 +1019,193 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces(
         std::filesystem::create_directories(forceMagnitudesDirectoryPath);
         std::ofstream jsonFile(patternMaxForceMagnitudesFilePath.string());
         jsonFile << json;
-    }
-#endif
-    assert(fullPatternSimulationScenarioMaxMagnitudes.size()
-           == desiredBaseSimulationScenarioIndices.size());
 
-    return fullPatternSimulationScenarioMaxMagnitudes;
+#endif
+        assert(fullPatternSimulationScenarioMaxMagnitudes.size()
+               == desiredBaseSimulationScenarioIndices.size());
+
+        return fullPatternSimulationScenarioMaxMagnitudes;
 }
 
 void ReducedModelOptimizer::runOptimization(const Settings &settings,
                                             ReducedPatternOptimization::Results &results)
 {
     global.objectiveValueHistory.clear();
-    dlib::matrix<double, 0, 1> xMin(global.numberOfOptimizationParameters);
-    dlib::matrix<double, 0, 1> xMax(global.numberOfOptimizationParameters);
-    for (int i = 0; i < global.numberOfOptimizationParameters; i++) {
-        xMin(i) = settings.xRanges[i].min;
-        xMax(i) = settings.xRanges[i].max;
-    }
+    global.objectiveValueHistory_iteration.clear();
+    global.objectiveValueHistory.reserve(settings.numberOfFunctionCalls / 100);
+    global.objectiveValueHistory_iteration.reserve(settings.numberOfFunctionCalls / 100);
 
-    dlib::function_evaluation result_dlib;
-    double (*objF)(double, double, double, double, double,double,double) = &objective;
-    result_dlib = dlib::find_min_global(objF,
-                                        xMin,
-                                        xMax,
-                                        dlib::max_function_calls(settings.numberOfFunctionCalls),
-                                        std::chrono::hours(24 * 365 * 290),
-                                        settings.solverAccuracy);
-    getResults(result_dlib, settings, results);
+#if POLYSCOPE_DEFINED
+    global.plotColors.reserve(settings.numberOfFunctionCalls);
+#endif
+    double (*objF)(const dlib::matrix<double, 0, 1> &) = &objective;
+
+    dlib::function_evaluation optimalResult;
+    const auto hexAngleParameterIt = std::find_if(settings.xRanges.begin(),
+                                                  settings.xRanges.end(),
+                                                  [](const xRange &x) {
+                                                      return x.label == "HexAngle";
+                                                  });
+    const bool hasHexAngleParameter = hexAngleParameterIt != settings.xRanges.end();
+    const auto hexSizeParameterIt = std::find_if(settings.xRanges.begin(),
+                                                 settings.xRanges.end(),
+                                                 [](const xRange &x) {
+                                                     return x.label == "HexSize";
+                                                 });
+    const bool hasHexSizeParameter = hexSizeParameterIt != settings.xRanges.end();
+    const bool hasBothGeometricalParameters = hasHexAngleParameter && hasHexSizeParameter;
+    assert(hasBothGeometricalParameters || (!hasHexAngleParameter && !hasHexSizeParameter));
+    const int numberOfGeometryOptimizationParameters = hasBothGeometricalParameters ? 2 : 0;
+    const int numberOfMaterialOptimizationParameters = global.numberOfOptimizationParameters
+                                                       - numberOfGeometryOptimizationParameters;
+    const bool hasGeometryAndMaterialParameters = hasBothGeometricalParameters
+                                                  && numberOfMaterialOptimizationParameters != 0;
+    if (settings.splitGeometryMaterialOptimization && hasGeometryAndMaterialParameters) {
+        //Geometry optimization of the reduced pattern
+        dlib::matrix<double, 0, 1> xGeometryMin(numberOfGeometryOptimizationParameters);
+        dlib::matrix<double, 0, 1> xGeometryMax(numberOfGeometryOptimizationParameters);
+        //        const int hexAngleParameterIndex = std::distance(settings.xRanges.begin(),
+        //                                                         hexAngleParameterIt);
+        //        const int hexSizeParameterIndex = std::distance(settings.xRanges.begin(),
+        //                                                        hexSizeParameterIt);
+        xGeometryMin(0) = hexSizeParameterIt->min;
+        xGeometryMax(0) = hexSizeParameterIt->max;
+        xGeometryMin(1) = hexAngleParameterIt->min;
+        xGeometryMax(1) = hexAngleParameterIt->max;
+
+        //Set reduced pattern update functions to be used during optimization
+        function_updateReducedPattern =
+            [&](const dlib::matrix<double, 0, 1> &x,
+                std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+                function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
+                pReducedPatternSimulationMesh->reset();
+            };
+        const int optimizationFunctionCalls_geometry
+            = global.numberOfOptimizationParameters == numberOfGeometryOptimizationParameters
+                  ? settings.numberOfFunctionCalls
+                  : static_cast<int>(settings.numberOfFunctionCalls / 3.0);
+        dlib::function_evaluation result_dlib_geometry
+            = dlib::find_min_global(objF,
+                                    xGeometryMin,
+                                    xGeometryMax,
+                                    dlib::max_function_calls(optimizationFunctionCalls_geometry),
+                                    std::chrono::hours(24 * 365 * 290),
+                                    settings.solverAccuracy);
+        //Material optimization of the reduced pattern
+        //        std::cout << "opt size:" << result_dlib_geometry.x(0) << std::endl;
+        //        std::cout << "opt size:" << global.minX[0] << std::endl;
+        dlib::function_evaluation result_dlib_material;
+        if (numberOfMaterialOptimizationParameters != 0) {
+            dlib::matrix<double, 0, 1> xMaterialMin(numberOfMaterialOptimizationParameters);
+            dlib::matrix<double, 0, 1> xMaterialMax(numberOfMaterialOptimizationParameters);
+            for (int i = 0; i < numberOfMaterialOptimizationParameters; i++) {
+                xMaterialMin(i) = settings.xRanges[i].min;
+                xMaterialMax(i) = settings.xRanges[i].max;
+            }
+            function_updateReducedPattern =
+                [&](const dlib::matrix<double, 0, 1> &x,
+                    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+                    function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
+                    function_updateReducedPattern_geometry(
+                        result_dlib_geometry.x,
+                        pReducedPatternSimulationMesh); //FIXME: I shouldn't need to call this
+                    pReducedPatternSimulationMesh->reset();
+                };
+            result_dlib_material = dlib::find_min_global(objF,
+                                                         xMaterialMin,
+                                                         xMaterialMax,
+                                                         dlib::max_function_calls(static_cast<int>(
+                                                             2.0 * settings.numberOfFunctionCalls
+                                                             / 3.0)),
+                                                         std::chrono::hours(24 * 365 * 290),
+                                                         settings.solverAccuracy);
+        }
+        //    constexpr bool useBOBYQA = false;
+        //    if (useBOBYQA) {
+        //        const size_t npt = 2 * global.numberOfOptimizationParameters;
+        //        //      ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
+        //        const double rhobeg = 0.1;
+        //        //  const double rhobeg = 10;
+        //        const double rhoend = rhobeg * 1e-6;
+        //        //  const size_t maxFun = 10 * (x.size() ^ 2);
+        //        const size_t bobyqa_maxFunctionCalls = 200;
+        //        dlib::find_min_bobyqa(objF,
+        //                              result_dlib.x,
+        //                              npt,
+        //                              xMaterialMin,
+        //                              xMaterialMax,
+        //                              rhobeg,
+        //                              rhoend,
+        //                              bobyqa_maxFunctionCalls);
+        //    }
+
+        optimalResult.x.set_size(global.numberOfOptimizationParameters);
+        std::copy(result_dlib_material.x.begin(),
+                  result_dlib_material.x.end(),
+                  optimalResult.x.begin());
+        //                  debug_x.begin());
+        std::copy(result_dlib_geometry.x.begin(),
+                  result_dlib_geometry.x.end(),
+                  optimalResult.x.begin() + numberOfMaterialOptimizationParameters);
+        //        std::cout << "opt x:";
+        //        for (const auto optx : optimalResult.x) {
+        //            std::cout << optx << " ";
+        //        }
+        //        std::cout << std::endl;
+        function_updateReducedPattern =
+            [&](const dlib::matrix<double, 0, 1> &x,
+                std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+                function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
+                function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
+                pReducedPatternSimulationMesh->reset();
+            };
+        //        const auto meshLabel = global.reducedPatternSimulationJobs[0]->pMesh->getLabel();
+        //        global.reducedPatternSimulationJobs[0]->pMesh->setLabel("pre");
+        //        global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::reducedInitial);
+        //        polyscope::show();
+        //        global.reducedPatternSimulationJobs[0]->pMesh->unregister();
+        optimalResult.y = objective(optimalResult.x);
+    } else {
+        dlib::matrix<double, 0, 1> xMin(global.numberOfOptimizationParameters);
+        dlib::matrix<double, 0, 1> xMax(global.numberOfOptimizationParameters);
+        for (int i = 0; i < global.numberOfOptimizationParameters; i++) {
+            xMin(i) = settings.xRanges[i].min;
+            xMax(i) = settings.xRanges[i].max;
+        }
+        if (numberOfGeometryOptimizationParameters != 0
+            && numberOfMaterialOptimizationParameters != 0) {
+            function_updateReducedPattern =
+                [&](const dlib::matrix<double, 0, 1> &x,
+                    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+                    function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
+                    function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
+                    pReducedPatternSimulationMesh->reset();
+                };
+        } else if (numberOfGeometryOptimizationParameters == 0) {
+            function_updateReducedPattern =
+                [&](const dlib::matrix<double, 0, 1> &x,
+                    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+                    function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
+                    pReducedPatternSimulationMesh->reset();
+                };
+
+        } else {
+            function_updateReducedPattern =
+                [&](const dlib::matrix<double, 0, 1> &x,
+                    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+                    function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
+                    pReducedPatternSimulationMesh->reset();
+                };
+        }
+        optimalResult = dlib::find_min_global(objF,
+                                              xMin,
+                                              xMax,
+                                              dlib::max_function_calls(
+                                                  settings.numberOfFunctionCalls),
+                                              std::chrono::hours(24 * 365 * 290),
+                                              settings.solverAccuracy);
+    }
+    getResults(optimalResult, settings, results);
 }
 
 void ReducedModelOptimizer::constructAxialSimulationScenario(
@@ -1117,12 +1393,12 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni
             - desiredRotationAngle);
         saveJobToPath = "../convergingJobs";
     }
-    std::filesystem::path outputPath(std::filesystem::path(saveJobToPath)
-                                         .append(job.pMesh->getLabel())
-                                         .append("mag_" + global.currentScenarioName));
-    std::filesystem::create_directories(outputPath);
-    job.save(outputPath);
-    settings.save(outputPath);
+    //    std::filesystem::path outputPath(std::filesystem::path(saveJobToPath)
+    //                                         .append(job.pMesh->getLabel())
+    //                                         .append("mag_" + global.currentScenarioName));
+    //    std::filesystem::create_directories(outputPath);
+    //    job.save(outputPath);
+    //    settings.save(outputPath);
 
     std::cout << "Force:" << forceMagnitude << " Error is:" << vcg::math::ToDeg(error) << std::endl;
 #endif
@@ -1186,7 +1462,6 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
 {
     double forceMagnitude = 1;
     double minimumError;
-    double translationalOptimizationEpsilon;
     bool wasSuccessful = false;
     global.constructScenarioFunction = constructBaseScenarioFunctions[scenario];
     const double baseTriangleHeight = vcg::Distance(global.baseTriangle.cP(0),
@@ -1194,6 +1469,7 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
                                                      + global.baseTriangle.cP(2))
                                                         / 2);
     std::function<double(const double &)> objectiveFunction;
+    double translationalOptimizationEpsilon{baseTriangleHeight * 0.001};
     double objectiveEpsilon = translationalOptimizationEpsilon;
     objectiveFunction = &fullPatternMaxSimulationForceTranslationalObjective;
     global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
@@ -1359,12 +1635,13 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
         for (int simulationScenarioIndex : global.simulationScenarioIndices) {
             if (global.optimizationSettings.normalizationStrategy ==
                 Settings::NormalizationStrategy::Epsilon) {
-                const double epsilon_translationalDisplacement = global.optimizationSettings
-                                                                     .normalizationParameter;
+                const double epsilon_translationalDisplacement
+                    = global.optimizationSettings.translationNormalizationParameter;
                 global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
                     = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex],
                                epsilon_translationalDisplacement);
-                const double epsilon_rotationalDisplacement = vcg::math::ToRad(3.0);
+                const double epsilon_rotationalDisplacement = global.optimizationSettings
+                                                                  .rotationNormalizationParameter;
                 global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
                     = std::max(fullPatternAngularDistance[simulationScenarioIndex],
                                epsilon_rotationalDisplacement);
@@ -1415,10 +1692,11 @@ void ReducedModelOptimizer::optimize(
     results.baseTriangle = global.baseTriangle;
 
     DRMSimulationModel::Settings simulationSettings;
-    simulationSettings.maxDRMIterations = 200000;
-    simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
+//    simulationSettings.maxDRMIterations = 200000;
+//    simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
     simulationSettings.viscousDampingFactor = 5e-3;
     simulationSettings.useKineticDamping = true;
+
 //    simulationSettings.averageResidualForcesCriterionThreshold = 1e-5;
 //    simulationSettings.viscousDampingFactor = 1e-3;
 //    simulationSettings.beVerbose = true;
@@ -1426,7 +1704,7 @@ void ReducedModelOptimizer::optimize(
 //    simulationSettings.isDebugMode = true;
 //    simulationSettings.debugModeStep = 100000;
 #ifdef POLYSCOPE_DEFINED
-    const bool drawFullPatternSimulationResults = false;
+    constexpr bool drawFullPatternSimulationResults = false;
     if (drawFullPatternSimulationResults) {
         global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
             ReducedPatternOptimization::Colors::fullInitial);

src/reducedmodeloptimizer.hpp

@@ -67,9 +67,10 @@ public:
 
     SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
 
-    void initializePatterns(PatternGeometry &fullPattern,
-                            PatternGeometry &reducedPatterm,
-                            const int &optimizationParameters);
+    void initializePatterns(
+        PatternGeometry &fullPattern,
+        PatternGeometry &reducedPatterm,
+        const std::vector<ReducedPatternOptimization::xRange> &optimizationParameters);
 
     static void runSimulation(const std::string &filename, std::vector<double> &x);
 
@@ -177,6 +178,33 @@ public:
         const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
             &oppositeInterfaceViPairs,
         SimulationJob &job);
+    static std::function<void(const dlib::matrix<double, 0, 1> &x,
+                              std::shared_ptr<SimulationMesh> &m)>
+        function_updateReducedPattern;
+    static std::function<void(const dlib::matrix<double, 0, 1> &x,
+                              std::shared_ptr<SimulationMesh> &m)>
+        function_updateReducedPattern_material;
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                              const double &newE)>
+        function_updateReducedPattern_material_E;
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                              const double &newA)>
+        function_updateReducedPattern_material_A;
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                              const double &newI)>
+        function_updateReducedPattern_material_I;
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                              const double &newI2)>
+        function_updateReducedPattern_material_I2;
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                              const double &newI3)>
+        function_updateReducedPattern_material_I3;
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                              const double &newJ)>
+        function_updateReducedPattern_material_J;
+    static std::function<void(const dlib::matrix<double, 0, 1> &x,
+                              std::shared_ptr<SimulationMesh> &m)>
+        function_updateReducedPattern_geometry;
 
 private:
     static void computeDesiredReducedModelDisplacements(
@@ -191,8 +219,9 @@ private:
             &maxForceMagnitudes);
     void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern);
     void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
-    static void initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh,
-                                                 const int &optimizationParamters);
+    static void initializeOptimizationParameters(
+        const std::shared_ptr<SimulationMesh> &mesh,
+        const std::vector<ReducedPatternOptimization::xRange> &optimizationParamters);
 
     static double objective(long n, const double *x);
     DRMSimulationModel simulator;
@@ -211,6 +240,32 @@ private:
     getFullPatternMaxSimulationForces(
         const std::vector<ReducedPatternOptimization::BaseSimulationScenario>
             &desiredBaseSimulationScenarioIndices);
+    static double objective(const dlib::matrix<double, 0, 1> &x);
+    static void initializeUpdateReducedPatternFunctions();
 };
-void updateMesh(long n, const double *x);
+inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
+    ReducedModelOptimizer::function_updateReducedPattern;
+inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
+    ReducedModelOptimizer::function_updateReducedPattern_material;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newE)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_E;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newA)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_A;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newI)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_I;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newI2)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_I2;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newI3)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_I3;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newJ)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_J;
+inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
+    ReducedModelOptimizer::function_updateReducedPattern_geometry;
+
 #endif // REDUCEDMODELOPTIMIZER_HPP