Made normalization strategy member of the optimization settings.

CMakeLists.txt

@@ -37,16 +37,15 @@ download_project(PROJ               MYSOURCES
 
 endif()
 
-
+##Polyscope
-#Polyscope
+#download_project(PROJ               POLYSCOPE
-download_project(PROJ               POLYSCOPE
+#                 GIT_REPOSITORY     https://github.com/nmwsharp/polyscope.git
-                 GIT_REPOSITORY     https://github.com/nmwsharp/polyscope.git
+#                 GIT_TAG	    master
-                 GIT_TAG	    master
+#                 PREFIX             ${CMAKE_CURRENT_SOURCE_DIR}/build/external/
-                 PREFIX             ${CMAKE_CURRENT_SOURCE_DIR}/build/external/
+#                 ${UPDATE_DISCONNECTED_IF_AVAILABLE}
-                 ${UPDATE_DISCONNECTED_IF_AVAILABLE}
+#)
-)
+#add_subdirectory(${POLYSCOPE_SOURCE_DIR})
-add_subdirectory(${POLYSCOPE_SOURCE_DIR})
+#add_compile_definitions(POLYSCOPE_DEFINED)
-add_compile_definitions(POLYSCOPE_DEFINED)
 
 #dlib
 download_project(PROJ               DLIB
@@ -89,4 +88,5 @@ target_include_directories(${PROJECT_NAME}
         )
 
         link_directories(${MYSOURCES_SOURCE_DIR}/boost_graph/libs)
-target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen matplot dlib::dlib)
+#target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen matplot dlib::dlib)
+target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib)

src/main.cpp

@@ -2,9 +2,6 @@
 #include "csvfile.hpp"
 #include "edgemesh.hpp"
 #include "flatpattern.hpp"
-#include "polyscope/curve_network.h"
-#include "polyscope/point_cloud.h"
-#include "polyscope/polyscope.h"
 #include "reducedmodeloptimizer.hpp"
 #include "simulationhistoryplotter.hpp"
 #include "trianglepattterntopology.hpp"
@@ -16,6 +13,11 @@
 #include <string>
 #include <vcg/complex/algorithms/update/position.h>
 
+#ifdef POLYSCOPE_DEFINED
+#include "polyscope/curve_network.h"
+#include "polyscope/point_cloud.h"
+#include "polyscope/polyscope.h"
+#endif
 int main(int argc, char *argv[]) {
   if (argc < 3) {
     std::cerr << "Specify at least the two pattern filepaths to be "
@@ -49,7 +51,10 @@ int main(int argc, char *argv[]) {
   const bool input_numberOfFunctionCallsDefined = argc >= 4;
   settings_optimization.numberOfFunctionCalls =
       input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100;
-  settings_optimization.shouldNormalizeObjectiveValue = true;
+  settings_optimization.normalizationStrategy =
+      ReducedModelOptimizer::Settings::NormalizationStrategy::Epsilon;
+  settings_optimization.normalizationParameter = 0.003;
+  settings_optimization.solutionAccuracy = 0.01;
 
   // Optimize pair
   const std::string pairName =
@@ -65,21 +70,24 @@ int main(int argc, char *argv[]) {
   const bool input_resultDirectoryDefined = argc >= 5;
   std::string optimizationResultsDirectory =
       input_resultDirectoryDefined ? argv[4] : "OptimizationResults";
+  std::string resultsOutputDir;
   if (optimizationResults.numberOfSimulationCrashes != 0) {
     const auto crashedJobsDirPath =
         std::filesystem::path(optimizationResultsDirectory)
             .append("CrashedJobs")
             .append(pairName);
     std::filesystem::create_directories(crashedJobsDirPath);
-    optimizationResults.save(crashedJobsDirPath.string());
+    resultsOutputDir = crashedJobsDirPath.string();
   } else {
     std::filesystem::path convergedJobsDirPath(
         std::filesystem::path(optimizationResultsDirectory)
             .append("ConvergedJobs")
             .append(pairName));
     std::filesystem::create_directories(convergedJobsDirPath);
-    optimizationResults.save(convergedJobsDirPath.string());
+    resultsOutputDir = convergedJobsDirPath.string();
   }
+  optimizationResults.save(resultsOutputDir);
+  // Write results in csv
   csvFile csv_results({}, false);
   //    csvFile csv_results(std::filesystem::path(dirPath_thisOptimization)
   //                            .append("results.csv")

src/reducedmodeloptimizer.cpp

@@ -9,7 +9,7 @@
 struct GlobalOptimizationVariables {
   std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
   std::vector<std::vector<Vector6d>> fullPatternDisplacements;
-  std::vector<double> reducedPatternMaximumDisplacementNormSum;
+  std::vector<double> objectiveNormalizationValues;
   std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
   std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
   std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
@@ -46,16 +46,10 @@ double ReducedModelOptimizer::computeError(
   const double rawError =
       computeRawError(reducedPatternDisplacements, fullPatternDisplacements,
                       reducedToFullInterfaceViMap);
-  if (global.optimizationSettings.shouldNormalizeObjectiveValue) {
+  if (global.optimizationSettings.normalizationStrategy !=
-    assert(rawError <= normalizationFactor);
+      Settings::NormalizationStrategy::NonNormalized) {
-    if (rawError > normalizationFactor) {
+    return rawError / normalizationFactor;
-      std::cout << "Normalized error is bigger than one." << std::endl;
-    }
-    return rawError /
-           normalizationFactor; // std::max(interfaceDisplacementsNormSum,
-                                // 0.00003);
   }
-
   return rawError;
 }
 
@@ -66,24 +60,22 @@ double ReducedModelOptimizer::computeRawError(
         &reducedToFullInterfaceViMap) {
   double error = 0;
   for (const auto reducedFullViPair : reducedToFullInterfaceViMap) {
-    //    const auto pos =
+    const VertexIndex reducedModelVi = reducedFullViPair.first;
-    //        g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP();
-    //    std::cout << "Interface vi " << reducedModelVi << " is at position "
-    //              << pos[0] << " " << pos[1] << " " << pos[2] << std::endl;
     Eigen::Vector3d reducedVertexDisplacement(
-        reducedPatternDisplacements[reducedFullViPair.first][0],
+        reducedPatternDisplacements[reducedModelVi][0],
-        reducedPatternDisplacements[reducedFullViPair.first][1],
+        reducedPatternDisplacements[reducedModelVi][1],
-        reducedPatternDisplacements[reducedFullViPair.first][2]);
+        reducedPatternDisplacements[reducedModelVi][2]);
     if (!std::isfinite(reducedVertexDisplacement[0]) ||
         !std::isfinite(reducedVertexDisplacement[1]) ||
         !std::isfinite(reducedVertexDisplacement[2])) {
       std::cout << "Displacements are not finite" << std::endl;
       std::terminate();
     }
+    const VertexIndex fullModelVi = reducedFullViPair.second;
     Eigen::Vector3d fullVertexDisplacement(
-        fullPatternDisplacements[reducedFullViPair.second][0],
+        fullPatternDisplacements[fullModelVi][0],
-        fullPatternDisplacements[reducedFullViPair.second][1],
+        fullPatternDisplacements[fullModelVi][1],
-        fullPatternDisplacements[reducedFullViPair.second][2]);
+        fullPatternDisplacements[fullModelVi][2]);
     Eigen::Vector3d errorVector =
         fullVertexDisplacement - reducedVertexDisplacement;
     //    error += errorVector.squaredNorm();
@@ -179,7 +171,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
   //            << std::endl;
 
   // run simulations
-  double error = 0;
+  double totalError = 0;
   FormFinder simulator;
   FormFinder::Settings simulationSettings;
   //    simulationSettings.shouldDraw = true;
@@ -192,75 +184,54 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
         simulationSettings);
     //    std::string filename;
     if (!reducedModelResults.converged) {
-      error += std::numeric_limits<double>::max();
+      totalError += std::numeric_limits<double>::max();
       global.numOfSimulationCrashes++;
       bool beVerbose = false; // TODO: use an extern or data member for that
       if (beVerbose) {
         std::cout << "Failed simulation" << std::endl;
       }
     } else {
-      const double thisSimulationScenarioError =
+      double simulationScenarioError = computeError(
-          computeError(reducedModelResults.displacements,
+          reducedModelResults.displacements,
-                       global.fullPatternDisplacements[simulationScenarioIndex],
+          global.fullPatternDisplacements[simulationScenarioIndex],
-                       global.reducedToFullInterfaceViMap,
+          global.reducedToFullInterfaceViMap,
-                       global.reducedPatternMaximumDisplacementNormSum
+          global.objectiveNormalizationValues[simulationScenarioIndex]);
-                           [simulationScenarioIndex]);
-      if (thisSimulationScenarioError > 1) {
-        std::cout << "Simulation scenario "
-                  << simulationScenarioStrings[simulationScenarioIndex]
-                  << " results in an error bigger than one." << std::endl;
-        for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
-          std::cout << "x[" + std::to_string(parameterIndex) + "]="
-                    << x[parameterIndex] << std::endl;
-        }
 
-        //#ifdef POLYSCOPE_DEFINED
+      //        if (global.optimizationSettings.normalizationStrategy !=
-        //        ReducedModelOptimizer::visualizeResults(
+      //        NormalizationStrategy::Epsilon &&
-        //            global.fullPatternSimulationJobs[simulationScenarioIndex],
+      //            simulationScenarioError > 1) {
-        //            global.reducedPatternSimulationJobs[simulationScenarioIndex],
+      //          std::cout << "Simulation scenario "
-        //            global.reducedToFullInterfaceViMap, false);
+      //                    <<
+      //                    simulationScenarioStrings[simulationScenarioIndex]
+      //                    << " results in an error bigger than one." <<
+      //                    std::endl;
+      //          for (size_t parameterIndex = 0; parameterIndex < n;
+      //               parameterIndex++) {
+      //            std::cout << "x[" + std::to_string(parameterIndex) + "]="
+      //                      << x[parameterIndex] << std::endl;
+      //          }
+      //        }
 
-        //        ReducedModelOptimizer::visualizeResults(
+      //#ifdef POLYSCOPE_DEFINED
-        //            global.fullPatternSimulationJobs[simulationScenarioIndex],
+      //        ReducedModelOptimizer::visualizeResults(
-        //            std::make_shared<SimulationJob>(
+      //            global.fullPatternSimulationJobs[simulationScenarioIndex],
-        //                reducedPatternMaximumDisplacementSimulationJobs
+      //            global.reducedPatternSimulationJobs[simulationScenarioIndex],
-        //                    [simulationScenarioIndex]),
+      //            global.reducedToFullInterfaceViMap, false);
-        //            global.reducedToFullInterfaceViMap, true);
+
-        //        polyscope::removeAllStructures();
+      //        ReducedModelOptimizer::visualizeResults(
-        //#endif // POLYSCOPE_DEFINED
+      //            global.fullPatternSimulationJobs[simulationScenarioIndex],
-      }
+      //            std::make_shared<SimulationJob>(
-      error += thisSimulationScenarioError;
+      //                reducedPatternMaximumDisplacementSimulationJobs
+      //                    [simulationScenarioIndex]),
+      //            global.reducedToFullInterfaceViMap, true);
+      //        polyscope::removeAllStructures();
+      //#endif // POLYSCOPE_DEFINED
+      totalError += simulationScenarioError;
     }
-    //    std::ofstream out(filename, std::ios_base::app);
-    //    auto pMesh =
-    //        global
-    //            .reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
-    //            ->pMesh;
-    //    for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++)
-    //    {
-    //      out << "x[" + std::to_string(parameterIndex) + "]=" <<
-    //      x[parameterIndex]
-    //          << std::endl;
-    //    }
-    //    out << pMesh->elements[0].dimensions.toString() + "\n" +
-    //               pMesh->elements[0].material.toString() + " \nA="
-    //        << pMesh->elements[0].A << " \nratio="
-    //        << pMesh->elements[0].dimensions.b /
-    //        pMesh->elements[0].dimensions.h
-    //        << " \naxialRig:" << pMesh->elements[0].rigidity.axial
-    //        << " \ntorsionalRig:" << pMesh->elements[0].rigidity.torsional
-    //        << " \nfirstBendingRig:" <<
-    //        pMesh->elements[0].rigidity.firstBending
-    //        << " \nsecondBendingRig:" <<
-    //        pMesh->elements[0].rigidity.secondBending
-    //        << " \nscenario:" +
-    //        simulationScenarioStrings[simulationScenarioIndex] +
-    //               "\n\n";
-    //    out.close();
   }
   //  std::cout << error << std::endl;
-  if (error < global.minY) {
+  if (totalError < global.minY) {
-    global.minY = error;
+    global.minY = totalError;
     global.minX.assign(x, x + n);
   }
   // if (++global.numberOfFunctionCalls %100== 0) {
@@ -269,7 +240,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
   //}
 
   // compute error and return it
-  global.gObjectiveValueHistory.push_back(error);
+  //  global.gObjectiveValueHistory.push_back(error);
   //  auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
   //                                 gObjectiveValueHistory.size());
   //  std::vector<double> colors(gObjectiveValueHistory.size(), 2);
@@ -282,7 +253,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
   //  value",
   //                                        gObjectiveValueHistory);
 
-  return error;
+  return totalError;
 }
 
 void ReducedModelOptimizer::createSimulationMeshes(
@@ -578,11 +549,73 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob(
 
 #if POLYSCOPE_DEFINED
 void ReducedModelOptimizer::visualizeResults(
+    const std::vector<std::shared_ptr<SimulationJob>>
+        &fullPatternSimulationJobs,
+    const std::vector<std::shared_ptr<SimulationJob>>
+        &reducedPatternSimulationJobs,
+    const std::vector<SimulationScenario> &simulationScenarios,
+    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+        &reducedToFullInterfaceViMap) {
+  FormFinder simulator;
+  std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
+      fullPatternSimulationJobs[0]->pMesh;
+  pFullPatternSimulationMesh->registerForDrawing();
+  pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() +
+                                      "_undeformed.ply");
+  reducedPatternSimulationJobs[0]->pMesh->savePly(
+      reducedPatternSimulationJobs[0]->pMesh->getLabel() + "_undeformed.ply");
+  double totalError = 0;
+  for (const int simulationScenarioIndex : simulationScenarios) {
+    const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
+        fullPatternSimulationJobs[simulationScenarioIndex];
+    pFullPatternSimulationJob->registerForDrawing(
+        pFullPatternSimulationMesh->getLabel());
+    SimulationResults fullModelResults =
+        simulator.executeSimulation(pFullPatternSimulationJob);
+    fullModelResults.registerForDrawing();
+    // fullModelResults.saveDeformedModel();
+    const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
+        reducedPatternSimulationJobs[simulationScenarioIndex];
+    SimulationResults reducedModelResults =
+        simulator.executeSimulation(pReducedPatternSimulationJob);
+    double normalizationFactor = 1;
+    if (global.optimizationSettings.normalizationStrategy !=
+        Settings::NormalizationStrategy::NonNormalized) {
+      normalizationFactor =
+          global.objectiveNormalizationValues[simulationScenarioIndex];
+    }
+    reducedModelResults.saveDeformedModel();
+    fullModelResults.saveDeformedModel();
+    double error = computeError(
+        reducedModelResults.displacements, fullModelResults.displacements,
+        reducedToFullInterfaceViMap, normalizationFactor);
+    std::cout << "Error of simulation scenario "
+              << simulationScenarioStrings[simulationScenarioIndex] << " is "
+              << error << std::endl;
+    totalError += error;
+    reducedModelResults.registerForDrawing();
+    //    firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
+    //    registerWorldAxes();
+    const std::string screenshotFilename =
+        "/home/iason/Coding/Projects/Approximating shapes with flat "
+        "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
+        "Images/" +
+        pFullPatternSimulationMesh->getLabel() + "_" +
+        simulationScenarioStrings[simulationScenarioIndex];
+    polyscope::show();
+    polyscope::screenshot(screenshotFilename, false);
+    fullModelResults.unregister();
+    reducedModelResults.unregister();
+    //    firstOptimizationRoundResults[simulationScenarioIndex].unregister();
+  }
+  std::cout << "Total error:" << totalError << std::endl;
+}
+
+void ReducedModelOptimizer::registerResultsForDrawing(
     const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
     const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
     const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
-        &reducedToFullInterfaceViMap,
+        &reducedToFullInterfaceViMap) {
-    const bool &draw = true) {
   FormFinder simulator;
   std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
       pFullPatternSimulationJob->pMesh;
@@ -609,19 +642,6 @@ void ReducedModelOptimizer::visualizeResults(
       global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/);
   std::cout << "Error is " << error << std::endl;
   reducedModelResults.registerForDrawing();
-  if (draw) {
-    polyscope::show();
-    const std::string screenshotFilename =
-        "/home/iason/Coding/Projects/Approximating shapes with flat "
-        "patterns/RodModelOptimizationForPatterns/Results/"
-        "Images/" +
-        pFullPatternSimulationMesh->getLabel() + "_" + "noScenarioName";
-    //        simulationScenarioStrings[simulationScenarioIndex];
-    polyscope::screenshot(screenshotFilename, false);
-    fullModelResults.unregister();
-    reducedModelResults.unregister();
-    pFullPatternSimulationMesh->unregister();
-  }
 }
 #endif // POLYSCOPE_DEFINED
 
@@ -651,46 +671,6 @@ ReducedModelOptimizer::Results
 ReducedModelOptimizer::runOptimization(const Settings &settings) {
 
   global.gObjectiveValueHistory.clear();
-
-  // g_optimizeInnerHexagonSize ? 5: 4;
-  //  const size_t npt = (n + 1) * (n + 2) / 2;
-  //  //      ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
-  //  assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2);
-  //  assert(npt <= 2 * n + 1 && "The choice of the number of interpolation "
-  //                             "conditions is not recommended.");
-  // Set initial guess of solution
-
-  // const size_t initialGuess = 1;
-  // std::vector<double> x(n, initialGuess);
-  // if (global.optimizeInnerHexagonSize) {
-  //  x[n - 1] = global.g_innerHexagonInitialPos;
-  //}
-  /*if (!initialGuess.empty()) {
-    x = g_optimizationInitialGuess;
-  }*/ //      {0.10000000000000	001,
-  //  2, 1.9999999971613847, 6.9560343643347764};
-  //      {1, 5.9277};
-  //      {0.0001, 2, 2.000000005047502, 1.3055270196964464};
-  //      {initialGuess(0), initialGuess(1), initialGuess(2),
-  //      initialGuess(3)};
-  //  assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) {
-  //           return d >= xMax || d <= xMin;
-  //         }));
-  //  std::vector<double> xLow(x.size(), xMin);
-  //  std::vector<double> xUpper(x.size(), xMax);
-  //  if (g_optimizeInnerHexagonSize) {
-  //    xLow[n - 1] = 0.1;
-  //    xUpper[n - 1] = 0.9;
-  //  }
-  //  const double maxX = *std::max_element(
-  //      x.begin(), x.end(),
-  //      [](const double &a, const double &b) { return abs(a) < abs(b); });
-  //  const double rhobeg = std::min(0.95, 0.2 * maxX);
-  //  double rhobeg = 1;
-  //  double rhoend = rhobeg * 1e-8;
-  //  const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2;
-  //  std::vector<double> w(wSize);
-  //  const size_t maxFun = std::min(100.0 * (x.size() + 1), 1000.0);
   dlib::matrix<double, 0, 1> xMin(global.numberOfOptimizationParameters);
   dlib::matrix<double, 0, 1> xMax(global.numberOfOptimizationParameters);
   for (int i = 0; i < global.numberOfOptimizationParameters; i++) {
@@ -720,17 +700,59 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
   results.numberOfSimulationCrashes = global.numOfSimulationCrashes;
   results.x = global.minX;
   results.objectiveValue = global.minY;
-  if (settings.shouldNormalizeObjectiveValue) {
+  if (global.minY != result.y) {
-    global.optimizationSettings.shouldNormalizeObjectiveValue = false;
+    std::cerr << "Global min objective is not equal to result objective"
+              << std::endl;
+  }
+
+  // Compute raw objective value
+  if (global.optimizationSettings.normalizationStrategy !=
+      Settings::NormalizationStrategy::NonNormalized) {
+    auto temp = global.optimizationSettings.normalizationStrategy;
+    global.optimizationSettings.normalizationStrategy =
+        Settings::NormalizationStrategy::NonNormalized;
     if (global.optimizeInnerHexagonSize) {
-      results.objectiveValue = objective(4, results.x.data());
+      results.rawObjectiveValue = objective(4, results.x.data());
 
     } else {
-      results.objectiveValue =
+      results.rawObjectiveValue =
           objective(results.x[0], results.x[1], results.x[2]);
     }
-    global.optimizationSettings.shouldNormalizeObjectiveValue = true;
+    global.optimizationSettings.normalizationStrategy = temp;
+
+  } else {
+    results.rawObjectiveValue = results.objectiveValue;
   }
+
+  // Compute obj value per simulation scenario
+  updateMesh(results.x.size(), results.x.data());
+  results.objectiveValuePerSimulationScenario.resize(
+      NumberOfSimulationScenarios);
+  FormFinder::Settings simulationSettings;
+  FormFinder simulator;
+  for (int simulationScenarioIndex = 0;
+       simulationScenarioIndex < NumberOfSimulationScenarios;
+       simulationScenarioIndex++) {
+    SimulationResults reducedModelResults = simulator.executeSimulation(
+        global.reducedPatternSimulationJobs[simulationScenarioIndex],
+        simulationSettings);
+    const double error = computeError(
+        reducedModelResults.displacements,
+        global.fullPatternDisplacements[simulationScenarioIndex],
+        global.reducedToFullInterfaceViMap,
+        global.objectiveNormalizationValues[simulationScenarioIndex]);
+
+    results.objectiveValuePerSimulationScenario[simulationScenarioIndex] =
+        error;
+  }
+  //  if (result.y !=
+  //      std::accumulate(results.objectiveValuePerSimulationScenario.begin(),
+  //                      results.objectiveValuePerSimulationScenario.end(), 0))
+  //                      {
+  //    std::cerr
+  //        << "Sum of per scenario objectives is not equal to result objective"
+  //        << std::endl;
+  //  }
   results.time = elapsed.count() / 1000.0;
   const bool printDebugInfo = false;
   if (printDebugInfo) {
@@ -739,33 +761,6 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
     //  std::cout << result.x << endl;
     std::cout << "Objective value:" << global.minY << endl;
   }
-  //  std::cout << result.y << endl;
-  //  std::cout << minY << endl;
-  //  std::cout << "Time(sec):" << elapsed.count() << std::endl;
-  //  std::cout << "Max function evaluations:" << maxFun << std::endl;
-  //  std::cout << "Initial guess:" << initialGuess << std::endl;
-  //  const size_t maxFun = 200;
-  //  bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
-  //         rhobeg, rhoend, maxFun, w.data());
-  //  std::cout << "Finished first optimization round" << std::endl;
-  //  firstOptimizationRoundResults.resize(6);
-  //  for (int simulationScenarioIndex = SimulationScenario::Axial;
-  //       simulationScenarioIndex !=
-  //       SimulationScenario::NumberOfSimulationScenarios;
-  //       simulationScenarioIndex++) {
-  //    SimulationResults reducedModelResults = simulator.executeSimulation(
-  //        g_reducedPatternSimulationJob[simulationScenarioIndex], false,
-  //        false);
-  //    reducedModelResults.setLabelPrefix("FirstRound");
-  //    firstOptimizationRoundResults[simulationScenarioIndex] =
-  //        std::move(reducedModelResults);
-  //  }
-  //  g_firstRoundIterationIndex = gObjectiveValueHistory.size();
-  //  rhobeg *= 1e1;
-  //  //  rhoend *= 1e2;
-  //  bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
-  //         rhobeg, rhoend, maxFun, w.data());
-  //  std::cout << "Finished second optimization round" << std::endl;
 
   return results;
 }
@@ -784,37 +779,7 @@ ReducedModelOptimizer::createScenarios(
   const double forceMagnitude = 1;
   // Assuming the patterns lays on the x-y plane
   const CoordType patternPlaneNormal(0, 0, 1);
-  // Make the first interface node lay on the x axis
+  // Assumes that the first interface node lay on the x axis
-  //  const size_t fullPatternFirstInterfaceNodeIndex =
-  //      m_fullPatternOppositeInterfaceViMap.begin()->second;
-  //  CoordType fullPatternFirstInterfaceNodePosition =
-  //      m_pFullModelSimulationMesh->vert[fullPatternFirstInterfaceNodeIndex].cP();
-  //  CoordType centerOfMass(0, 0, 0);
-  //  for (size_t vi = 0; vi < pMesh->VN(); vi++) {
-  //    centerOfMass = centerOfMass + pMesh->vert[vi].P();
-  //  }
-  //  centerOfMass /= pMesh->VN();
-  //  vcg::tri::UpdatePosition<SimulationMesh>::Translate(
-  //      *m_pFullModelSimulationMesh, -centerOfMass);
-  //  vcg::tri::UpdatePosition<SimulationMesh>::Translate(
-  //      *m_pReducedPatternSimulationMesh, centerOfMass);
-
-  //  const vcg::Matrix33d R = vcg::RotationMatrix(
-  //      fullPatternFirstInterfaceNodePosition,
-  //      CoordType(fullPatternFirstInterfaceNodePosition.Norm(), 0, 0), false);
-  //  std::for_each(m_pFullModelSimulationMesh->vert.begin(),
-  //                m_pFullModelSimulationMesh->vert.end(), [&](auto &v) {
-  //                  v.P() = R * v.P();
-  //                  v.N() = R * v.N();
-  //                });
-  //  std::for_each(m_pReducedPatternSimulationMesh->vert.begin(),
-  //                m_pReducedPatternSimulationMesh->vert.end(), [&](auto &v) {
-  //                  v.P() = R * v.P();
-  //                  v.N() = R * v.N();
-  //                });
-
-  //  m_pFullModelSimulationMesh->updateEigenEdgeAndVertices();
-  //  m_pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
 
   //// Axial
   SimulationScenario scenarioName = SimulationScenario::Axial;
@@ -945,89 +910,87 @@ ReducedModelOptimizer::createScenarios(
   return scenarios;
 }
 
-// void ReducedModelOptimizer::runBeamOptimization() {
+void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
-//  // load beams
+  assert(global.optimizationSettings.shouldNormalizeObjectiveValue);
-//  VCGEdgeMesh fullBeam;
+  if (global.optimizationSettings.normalizationStrategy ==
-//  fullBeam.loadPly("/home/iason/Models/simple_beam_model_10elem_1m.ply");
+          Settings::NormalizationStrategy::Epsilon ||
-//  VCGEdgeMesh reducedBeam;
+      global.optimizationSettings.normalizationStrategy ==
-//  reducedBeam.loadPly("/home/iason/Models/simple_beam_model_4elem_1m.ply");
+          Settings::NormalizationStrategy::EqualDisplacements) {
-//  fullBeam.registerForDrawing();
-//  reducedBeam.registerForDrawing();
-//  //  polyscope::show();
-//  // maps
-//  std::unordered_map<size_t, size_t> displacementReducedToFullMap;
-//  displacementReducedToFullMap[reducedBeam.VN() / 2] = fullBeam.VN() / 2;
-//  g_reducedToFullViMap = displacementReducedToFullMap;
-//  std::unordered_map<size_t, size_t> jobFullToReducedMap;
-//  jobFullToReducedMap[0] = 0;
-//  jobFullToReducedMap[fullBeam.VN() - 1] = reducedBeam.VN() - 1;
 
-//  // full model simuilation job
+    // Compute the sum of the displacement norms
-//  auto pFullPatternSimulationMesh =
+    std::vector<double> fullPatternDisplacementNormSum(
-//  std::make_shared<SimulationMesh>(fullBeam);
+        NumberOfSimulationScenarios);
-//  pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.02,
+    for (int simulationScenarioIndex : global.simulationScenarioIndices) {
-//  0.02}); pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
+      double displacementNormSum = 0;
-//  std::unordered_map<VertexIndex, std::unordered_set<int>> fixedVertices;
+      for (auto interfaceViPair : global.reducedToFullInterfaceViMap) {
-//  fixedVertices[0] = ::unordered_set<int>({0, 1, 2, 3, 4, 5});
+        const Vector6d &vertexDisplacement =
-//  std::unordered_map<VertexIndex, Vector6d> forces;
+            global.fullPatternDisplacements[simulationScenarioIndex]
-//  forces[fullBeam.VN() - 1] = Vector6d({0, 0, 10, 0, 0, 0});
+                                           [interfaceViPair.second];
-//  const std::string fullBeamSimulationJobLabel = "Pull_Z";
+        displacementNormSum += vertexDisplacement.getTranslation().norm();
-//  std::shared_ptr<SimulationJob> pFullModelSimulationJob =
+      }
-//      make_shared<SimulationJob>(SimulationJob(pFullPatternSimulationMesh,
+      fullPatternDisplacementNormSum[simulationScenarioIndex] =
-//                                               fullBeamSimulationJobLabel,
+          displacementNormSum;
-//                                               fixedVertices, forces));
+    }
-//  auto fullModelResults =
-//  formFinder.executeSimulation(pFullModelSimulationJob);
 
-//  // Optimal reduced model displacements
+    for (int simulationScenarioIndex : global.simulationScenarioIndices) {
-//  const size_t numberOfSimulationScenarios = 1;
+      if (global.optimizationSettings.normalizationStrategy ==
-//  g_optimalReducedModelDisplacements.resize(numberOfSimulationScenarios);
+          Settings::NormalizationStrategy::Epsilon) {
-//  g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1].resize(
+        const double epsilon =
-//      reducedBeam.VN(), 3);
+            global.optimizationSettings.normalizationParameter;
-//  computeDesiredReducedModelDisplacements(
+        if (epsilon > fullPatternDisplacementNormSum[simulationScenarioIndex]) {
-//      fullModelResults, displacementReducedToFullMap,
+          //          std::cout << "Epsilon used in "
-//      g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]);
+          //                    <<
-
+          //                    simulationScenarioStrings[simulationScenarioIndex]
-//  // reduced model simuilation job
+          //                    << std::endl;
-//  auto reducedSimulationMesh = std::make_shared<SimulationMesh>(reducedBeam);
+        }
-//  reducedSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02});
+        global.objectiveNormalizationValues[simulationScenarioIndex] = std::max(
-//  reducedSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
+            fullPatternDisplacementNormSum[simulationScenarioIndex], epsilon);
-//  g_reducedPatternSimulationJob.resize(numberOfSimulationScenarios);
+        //            displacementNormSum;
-//  SimulationJob reducedSimJob;
+      } else { // NormalizationStrategy::EqualDisplacements
-//  computeReducedModelSimulationJob(*pFullModelSimulationJob,
+        global.objectiveNormalizationValues[simulationScenarioIndex] =
-//                                   jobFullToReducedMap, reducedSimJob);
+            fullPatternDisplacementNormSum[simulationScenarioIndex] /
-//  reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] =
+            fullPatternDisplacementNormSum[SimulationScenario::Shear];
-//      reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] * 0.1;
+      }
-//  g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1] =
+    }
-//      make_shared<SimulationJob>(
+  } else if (global.optimizationSettings.normalizationStrategy ==
-//          reducedSimulationMesh, fullBeamSimulationJobLabel,
+             Settings::NormalizationStrategy::MaxDisplacement) {
-//          reducedSimJob.constrainedVertices,
+    for (int simulationScenarioIndex : global.simulationScenarioIndices) {
-//          reducedSimJob.nodalExternalForces);
+      // Compute the quantities for normalizing the obj value
-//  initializeOptimizationParameters(reducedSimulationMesh);
+      const double minB = global.optimizationSettings.xRanges[0].min;
-
+      const double maxRatio = global.optimizationSettings.xRanges[1].max;
-//  //  const std::string simulationJobsPath = "SimulationJobs";
+      const double minE = global.optimizationSettings.xRanges[2].min;
-//  //  std::filesystem::create_directory(simulationJobsPath);
+      const double minHS = 0.3;
-//  //  g_reducedPatternSimulationJob[0].save(simulationJobsPath);
+      std::vector<double> mostFlexibleOptimizationParameters{minB, maxRatio,
-//  //  g_reducedPatternSimulationJob[0].load(
+                                                             minE, minHS};
-//  //      std::filesystem::path(simulationJobsPath)
+      if (global.optimizeInnerHexagonSize) {
-//  //          .append(g_reducedPatternSimulationJob[0].mesh->getLabel() +
+        updateMesh(4, mostFlexibleOptimizationParameters.data());
-//  //                  "_simScenario.json"));
+      } else {
-
+        updateMesh(3, mostFlexibleOptimizationParameters.data());
-//  runOptimization({}, &objective);
+      }
-
+      reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex] =
-//  fullModelResults.registerForDrawing();
+          global.reducedPatternSimulationJobs[simulationScenarioIndex]
-//  SimulationResults reducedModelResults = simulator.executeSimulation(
+              ->getCopy();
-//      g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1]);
+      reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex]
-//  double error = computeError(
+          .pMesh->setLabel("reduced_maxDisplacement");
-//      reducedModelResults,
+      SimulationResults reducedModelResults = simulator.executeSimulation(
-//      g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]);
+          global.reducedPatternSimulationJobs[simulationScenarioIndex]);
-//  reducedModelResults.registerForDrawing();
+      const double errorOfMaxDisplacedReduced = computeRawError(
-//  std::cout << "Error between beams:" << error << endl;
+          reducedModelResults.displacements,
-//  //    registerWorldAxes();
+          global.fullPatternDisplacements[simulationScenarioIndex],
-//  polyscope::show();
+          global.reducedToFullInterfaceViMap);
-//  fullModelResults.unregister();
+      const double errorOfNonDisplacedReduced = computeRawError(
-//  reducedModelResults.unregister();
+          std::vector<Vector6d>(reducedModelResults.displacements.size(),
-//}
+                                Vector6d(0)),
+          global.fullPatternDisplacements[simulationScenarioIndex],
+          global.reducedToFullInterfaceViMap);
+      const double displacementMultiplier =
+          global.optimizationSettings.normalizationParameter;
+      global.objectiveNormalizationValues[simulationScenarioIndex] =
+          displacementMultiplier *
+          std::max(errorOfMaxDisplacedReduced, errorOfNonDisplacedReduced);
+    }
+  }
+}
 
 ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
     const Settings &optimizationSettings,
@@ -1041,10 +1004,10 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
         SimulationScenario::Saddle};
   }
 
-  global.g_optimalReducedModelDisplacements.resize(6);
+  global.g_optimalReducedModelDisplacements.resize(NumberOfSimulationScenarios);
-  global.reducedPatternSimulationJobs.resize(6);
+  global.reducedPatternSimulationJobs.resize(NumberOfSimulationScenarios);
-  global.fullPatternDisplacements.resize(6);
+  global.fullPatternDisplacements.resize(NumberOfSimulationScenarios);
-  global.reducedPatternMaximumDisplacementNormSum.resize(6);
+  global.objectiveNormalizationValues.resize(NumberOfSimulationScenarios);
   global.g_firstRoundIterationIndex = 0;
   global.minY = std::numeric_limits<double>::max();
   global.numOfSimulationCrashes = 0;
@@ -1052,7 +1015,8 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
   global.optimizationSettings = optimizationSettings;
   global.fullPatternSimulationJobs =
       createScenarios(m_pFullPatternSimulationMesh);
-  reducedPatternMaximumDisplacementSimulationJobs.resize(6);
+  reducedPatternMaximumDisplacementSimulationJobs.resize(
+      NumberOfSimulationScenarios);
   //  polyscope::removeAllStructures();
 
   FormFinder::Settings simulationSettings;
@@ -1071,40 +1035,11 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
                                      reducedPatternSimulationJob);
     global.reducedPatternSimulationJobs[simulationScenarioIndex] =
         std::make_shared<SimulationJob>(reducedPatternSimulationJob);
+  }
 
-    if (optimizationSettings.shouldNormalizeObjectiveValue) {
+  if (global.optimizationSettings.normalizationStrategy !=
-      // Compute the quantities for normalizing the obj value
+      Settings::NormalizationStrategy::NonNormalized) {
-      const double minB = optimizationSettings.xRanges[0].min;
+    computeObjectiveValueNormalizationFactors();
-      const double maxRatio = optimizationSettings.xRanges[1].max;
-      const double minE = optimizationSettings.xRanges[2].min;
-      const double minHS = 0.3;
-      std::vector<double> mostFlexibleOptimizationParameters{minB, maxRatio,
-                                                             minE, minHS};
-      if (global.optimizeInnerHexagonSize) {
-        updateMesh(4, mostFlexibleOptimizationParameters.data());
-      } else {
-        updateMesh(3, mostFlexibleOptimizationParameters.data());
-      }
-      reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex] =
-          global.reducedPatternSimulationJobs[simulationScenarioIndex]
-              ->getCopy();
-      reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex]
-          .pMesh->setLabel("reduced_maxDisplacement");
-      SimulationResults reducedModelResults = simulator.executeSimulation(
-          global.reducedPatternSimulationJobs[simulationScenarioIndex],
-          simulationSettings);
-      const double errorOfMaxDisplacedReduced = computeRawError(
-          reducedModelResults.displacements, fullModelResults.displacements,
-          global.reducedToFullInterfaceViMap);
-      const double errorOfNonDisplacedReduced = computeRawError(
-          std::vector<Vector6d>(reducedModelResults.displacements.size(),
-                                Vector6d(0)),
-          fullModelResults.displacements, global.reducedToFullInterfaceViMap);
-      const double displacementMultiplier = 2;
-      global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex] =
-          displacementMultiplier *
-          std::max(errorOfMaxDisplacedReduced, errorOfNonDisplacedReduced);
-    }
   }
   Results optResults = runOptimization(optimizationSettings);
   for (int simulationScenarioIndex : global.simulationScenarioIndices) {
@@ -1113,12 +1048,9 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
     optResults.reducedPatternSimulationJobs.push_back(
         global.reducedPatternSimulationJobs[simulationScenarioIndex]);
   }
-  optResults.draw();
-
 #ifdef POLYSCOPE_DEFINED
-//  visualizeResults(
+  optResults.draw();
-//      global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs,
-//      global.simulationScenarioIndices, global.reducedToFullInterfaceViMap);
 #endif // POLYSCOPE_DEFINED
+
   return optResults;
 }

src/reducedmodeloptimizer.hpp

@@ -8,6 +8,9 @@
 #include "matplot/matplot.h"
 #include <Eigen/Dense>
 
+#ifdef POLYSCOPE_DEFINED
+#include "polyscope/color_management.h"
+#endif // POLYSCOPE_DEFINED
 using FullPatternVertexIndex = VertexIndex;
 using ReducedPatternVertexIndex = VertexIndex;
 
@@ -21,6 +24,21 @@ class ReducedModelOptimizer {
   std::unordered_map<size_t, size_t> nodeToSlot;
   std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
   std::vector<double> initialGuess;
+#ifdef POLYSCOPE_DEFINED
+  struct StaticColors {
+    glm::vec3 fullInitial;
+    glm::vec3 fullDeformed;
+    glm::vec3 reducedInitial;
+    glm::vec3 reducedDeformed;
+    StaticColors() {
+      fullInitial = {0.416666, 0.109804, 0.890196};
+      fullDeformed = {0.583333, 0.890196, 0.109804};
+      reducedInitial = {0.890196, 0.109804, 0.193138};
+      reducedDeformed = {0.109804, 0.890196, 0.806863};
+    }
+  };
+  inline static StaticColors colors;
+#endif // POLYSCOPE_DEFINED
 
 public:
   enum SimulationScenario {
@@ -44,10 +62,19 @@ public:
   struct Results;
 
   struct Settings {
+    enum NormalizationStrategy {
+      NonNormalized,
+      Epsilon,
+      MaxDisplacement,
+      EqualDisplacements
+    };
+    inline static vector<std::string> normalizationStrategyStrings{
+        "NonNormalized", "Epsilon", "MaxDsiplacement", "EqualDisplacements"};
     std::vector<xRange> xRanges;
     int numberOfFunctionCalls{100};
     double solutionAccuracy{1e-2};
-    bool shouldNormalizeObjectiveValue{true};
+    NormalizationStrategy normalizationStrategy{Epsilon};
+    double normalizationParameter{0.003};
 
     std::string toString() const {
       std::string settingsString;
@@ -61,7 +88,7 @@ public:
       settingsString +=
           "FuncCalls=" + std::to_string(numberOfFunctionCalls) +
           " Accuracy=" + std::to_string(solutionAccuracy) +
-          " Norm=" + (shouldNormalizeObjectiveValue ? "yes" : "no");
+          " Norm=" + normalizationStrategyStrings[normalizationStrategy];
 
       return settingsString;
     }
@@ -75,6 +102,7 @@ public:
       }
       os << "Function Calls";
       os << "Solution Accuracy";
+      os << "Normalization strategy";
       //        os << std::endl;
     }
 
@@ -87,11 +115,211 @@ public:
       }
       os << numberOfFunctionCalls;
       os << solutionAccuracy;
+      os << normalizationStrategyStrings[normalizationStrategy] + "_" +
+                std::to_string(normalizationParameter);
+    }
+  };
+  struct Results {
+    double time{-1};
+    int numberOfSimulationCrashes{0};
+    std::vector<double> x;
+    double objectiveValue;
+    double rawObjectiveValue;
+    std::vector<double> objectiveValuePerSimulationScenario;
+    std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
+    std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
+
+    void save(const string &saveToPath) const {
+      assert(std::filesystem::is_directory(saveToPath));
+
+      const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
+      for (int simulationJobIndex = 0;
+           simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
+        const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
+            fullPatternSimulationJobs[simulationJobIndex];
+        std::filesystem::path simulationJobFolderPath(
+            std::filesystem::path(saveToPath)
+                .append(pFullPatternSimulationJob->label));
+        std::filesystem::create_directory(simulationJobFolderPath);
+        const auto fullPatternDirectoryPath =
+            std::filesystem::path(simulationJobFolderPath).append("Full");
+        std::filesystem::create_directory(fullPatternDirectoryPath);
+        pFullPatternSimulationJob->save(fullPatternDirectoryPath.string());
+        const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
+            reducedPatternSimulationJobs[simulationJobIndex];
+        const auto reducedPatternDirectoryPath =
+            std::filesystem::path(simulationJobFolderPath).append("Reduced");
+        if (!std::filesystem::exists(reducedPatternDirectoryPath)) {
+          std::filesystem::create_directory(reducedPatternDirectoryPath);
+        }
+        pReducedPatternSimulationJob->save(
+            reducedPatternDirectoryPath.string());
+      }
+    }
+    void load(const string &loadFromPath) {
+      assert(std::filesystem::is_directory(loadFromPath));
+
+      for (const auto &directoryEntry :
+           filesystem::directory_iterator(loadFromPath)) {
+        const auto simulationScenarioPath = directoryEntry.path();
+        if (!std::filesystem::is_directory(simulationScenarioPath)) {
+          continue;
+        }
+        // Load reduced 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());
+            fullPatternSimulationJobs.push_back(
+                std::make_shared<SimulationJob>(job));
+          }
+        }
+
+        // Load full pattern files
+        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());
+            reducedPatternSimulationJobs.push_back(
+                std::make_shared<SimulationJob>(job));
+          }
+        }
+      }
+    }
+#if POLYSCOPE_DEFINED
+    void draw() const {
+      initPolyscope();
+      FormFinder simulator;
+      assert(fullPatternSimulationJobs.size() ==
+             reducedPatternSimulationJobs.size());
+      fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
+          colors.fullInitial);
+      reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(
+          colors.reducedInitial);
+
+      const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
+      for (int simulationJobIndex = 0;
+           simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
+        // Drawing of full pattern results
+        const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
+            fullPatternSimulationJobs[simulationJobIndex];
+        pFullPatternSimulationJob->registerForDrawing(
+            fullPatternSimulationJobs[0]->pMesh->getLabel());
+        SimulationResults fullModelResults =
+            simulator.executeSimulation(pFullPatternSimulationJob);
+        fullModelResults.registerForDrawing(colors.fullDeformed);
+        // Drawing of reduced pattern results
+        const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
+            reducedPatternSimulationJobs[simulationJobIndex];
+        SimulationResults reducedModelResults =
+            simulator.executeSimulation(pReducedPatternSimulationJob);
+        reducedModelResults.registerForDrawing(colors.reducedDeformed);
+        polyscope::options::programName =
+            fullPatternSimulationJobs[0]->pMesh->getLabel();
+        polyscope::show();
+        // Save a screensh
+        const std::string screenshotFilename =
+            "/home/iason/Coding/Projects/Approximating shapes with flat "
+            "patterns/RodModelOptimizationForPatterns/Results/Images/" +
+            pFullPatternSimulationJob->pMesh->getLabel() + "_" +
+            simulationScenarioStrings[simulationJobIndex];
+        polyscope::screenshot(screenshotFilename, false);
+        fullModelResults.unregister();
+        reducedModelResults.unregister();
+        //                    double error = computeError(
+        //                        reducedModelResults.displacements,fullModelResults.displacements,
+        //                        );
+        //                    std::cout << "Error of simulation scenario "
+        //                              <<
+        //                              simulationScenarioStrings[simulationScenarioIndex]
+        //                              << " is "
+        //                              << error << std::endl;
+      }
+    }
+#endif // POLYSCOPE_DEFINED
+    void saveMeshFiles() const {
+      const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
+      assert(numberOfSimulationJobs != 0 &&
+             fullPatternSimulationJobs.size() ==
+                 reducedPatternSimulationJobs.size());
+
+      fullPatternSimulationJobs[0]->pMesh->savePly(
+          "FullPattern_undeformed.ply");
+      reducedPatternSimulationJobs[0]->pMesh->savePly(
+          "ReducedPattern_undeformed.ply");
+      FormFinder simulator;
+      for (int simulationJobIndex = 0;
+           simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
+        // Drawing of full pattern results
+        const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
+            fullPatternSimulationJobs[simulationJobIndex];
+        SimulationResults fullModelResults =
+            simulator.executeSimulation(pFullPatternSimulationJob);
+        fullModelResults.saveDeformedModel();
+
+        // Drawing of reduced pattern results
+        const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
+            reducedPatternSimulationJobs[simulationJobIndex];
+        SimulationResults reducedModelResults =
+            simulator.executeSimulation(pReducedPatternSimulationJob);
+        reducedModelResults.saveDeformedModel();
+      }
+    }
+
+    void
+    writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization,
+                  csvFile &os) {
+      os << "Total raw Obj value";
+      os << "Total Obj value";
+      for (int simulationScenarioIndex = 0;
+           simulationScenarioIndex <
+           SimulationScenario::NumberOfSimulationScenarios;
+           simulationScenarioIndex++) {
+        os << "Obj value " + simulationScenarioStrings[simulationScenarioIndex];
+      }
+      for (const ReducedModelOptimizer::xRange &range :
+           settings_optimization.xRanges) {
+        os << range.label;
+      }
+      os << "Time(s)";
+      os << "#Crashes";
+    }
+
+    void
+    writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization,
+                   csvFile &os) const {
+      os << rawObjectiveValue;
+      os << objectiveValue;
+      for (double scenarioObjValue : objectiveValuePerSimulationScenario) {
+        os << scenarioObjValue;
+      }
+      for (const double &optimalX : x) {
+        os << optimalX;
+      }
+
+      for (int unusedXVarCounter = 0;
+           unusedXVarCounter < settings_optimization.xRanges.size() - x.size();
+           unusedXVarCounter++) {
+        os << "-";
+      }
+
+      os << time;
+      if (numberOfSimulationCrashes == 0) {
+        os << "-";
+      } else {
+        os << numberOfSimulationCrashes;
+      }
     }
   };
 
   inline static const std::string simulationScenarioStrings[] = {
-      "Axial", "Shear", "Bending", "Double", "Saddle"};
+      "Axial", "Shear", "Bending", "Dome", "Saddle"};
   Results optimize(const Settings &xRanges,
                    const std::vector<SimulationScenario> &simulationScenarios =
                        std::vector<SimulationScenario>());
@@ -118,7 +346,7 @@ public:
                             std::vector<double> &x);
 
   static double objective(double x0, double x1, double x2, double x3);
-  static double objective(double b, double h, double E);
+  static double objective(double b, double r, double E);
 
   static std::vector<std::shared_ptr<SimulationJob>>
   createScenarios(const std::shared_ptr<SimulationMesh> &pMesh,
@@ -138,20 +366,36 @@ public:
           &fullToReducedInterfaceViMap,
       std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
           &fullPatternOppositeInterfaceViMap);
-  static void visualizeResults(
+  static void
-      const std::shared_ptr<SimulationJob> &fullPatternSimulationJobs,
+  visualizeResults(const std::vector<std::shared_ptr<SimulationJob>>
-      const std::shared_ptr<SimulationJob> &reducedPatternSimulationJobs,
+                       &fullPatternSimulationJobs,
+                   const std::vector<std::shared_ptr<SimulationJob>>
+                       &reducedPatternSimulationJobs,
+                   const std::vector<SimulationScenario> &simulationScenarios,
+                   const std::unordered_map<ReducedPatternVertexIndex,
+                                            FullPatternVertexIndex>
+                       &reducedToFullInterfaceViMap);
+  static void registerResultsForDrawing(
+      const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
+      const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
       const std::unordered_map<ReducedPatternVertexIndex,
                                FullPatternVertexIndex>
-          &reducedToFullInterfaceViMap,
+          &reducedToFullInterfaceViMap);
-      const bool &shouldDraw);
 
-  static double computeError(const std::vector<Vector6d> &reducedPatternResults,
+  static double
-                             const std::vector<Vector6d> &fullPatternResults,
+  computeRawError(const std::vector<Vector6d> &reducedPatternDisplacements,
-                             const std::unordered_map<ReducedPatternVertexIndex,
+                  const std::vector<Vector6d> &fullPatternDisplacements,
-                                                      FullPatternVertexIndex>
+                  const std::unordered_map<ReducedPatternVertexIndex,
-                                 &reducedToFullInterfaceViMap,
+                                           FullPatternVertexIndex>
-                             const double &normalizationFactor);
+                      &reducedToFullInterfaceViMap);
+
+  static double
+  computeError(const std::vector<Vector6d> &reducedPatternDisplacements,
+               const std::vector<Vector6d> &fullPatternDisplacements,
+               const std::unordered_map<ReducedPatternVertexIndex,
+                                        FullPatternVertexIndex>
+                   &reducedToFullInterfaceViMap,
+               const double &normalizationFactor);
 
 private:
   static void computeDesiredReducedModelDisplacements(
@@ -168,194 +412,8 @@ private:
   static void
   initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh);
 
-  static double
-  computeError(const SimulationResults &reducedPatternResults,
-               const Eigen::MatrixX3d &optimalReducedPatternDisplacements);
   static double objective(long n, const double *x);
   FormFinder simulator;
-  double static computeRawError(
+  void computeObjectiveValueNormalizationFactors();
-      const std::vector<Vector6d> &reducedPatternDisplacements,
-      const std::vector<Vector6d> &fullPatternDisplacements,
-      const std::unordered_map<ReducedPatternVertexIndex,
-                               FullPatternVertexIndex>
-          &reducedToFullInterfaceViMap);
 };
-struct ReducedModelOptimizer::Results {
-  double time{-1};
-  int numberOfSimulationCrashes{0};
-  std::vector<double> x;
-  double objectiveValue;
-  std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
-  std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
-
-  void save(const string &saveToPath) const {
-    assert(std::filesystem::is_directory(saveToPath));
-
-    const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
-    for (int simulationJobIndex = 0;
-         simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
-      const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
-          fullPatternSimulationJobs[simulationJobIndex];
-      std::filesystem::path simulationJobFolderPath(
-          std::filesystem::path(saveToPath)
-              .append(pFullPatternSimulationJob->label));
-      std::filesystem::create_directory(simulationJobFolderPath);
-      const auto fullPatternDirectoryPath =
-          std::filesystem::path(simulationJobFolderPath).append("Full");
-      std::filesystem::create_directory(fullPatternDirectoryPath);
-      pFullPatternSimulationJob->save(fullPatternDirectoryPath.string());
-      const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
-          reducedPatternSimulationJobs[simulationJobIndex];
-      const auto reducedPatternDirectoryPath =
-          std::filesystem::path(simulationJobFolderPath).append("Reduced");
-      if (!std::filesystem::exists(reducedPatternDirectoryPath)) {
-        std::filesystem::create_directory(reducedPatternDirectoryPath);
-      }
-      pReducedPatternSimulationJob->save(reducedPatternDirectoryPath.string());
-    }
-  }
-  void load(const string &loadFromPath) {
-    assert(std::filesystem::is_directory(loadFromPath));
-
-    for (const auto &directoryEntry :
-         filesystem::directory_iterator(loadFromPath)) {
-      const auto simulationScenarioPath = directoryEntry.path();
-      if (!std::filesystem::is_directory(simulationScenarioPath)) {
-        continue;
-      }
-      // Load reduced 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());
-          fullPatternSimulationJobs.push_back(
-              std::make_shared<SimulationJob>(job));
-        }
-      }
-
-      // Load full pattern files
-      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());
-          reducedPatternSimulationJobs.push_back(
-              std::make_shared<SimulationJob>(job));
-        }
-      }
-    }
-  }
-#if POLYSCOPE_DEFINED
-  void draw() const {
-    initPolyscope();
-    FormFinder simulator;
-    assert(fullPatternSimulationJobs.size() ==
-           reducedPatternSimulationJobs.size());
-    fullPatternSimulationJobs[0]->pMesh->registerForDrawing();
-    reducedPatternSimulationJobs[0]->pMesh->registerForDrawing();
-
-    const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
-    for (int simulationJobIndex = 0;
-         simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
-      // Drawing of full pattern results
-      const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
-          fullPatternSimulationJobs[simulationJobIndex];
-      pFullPatternSimulationJob->registerForDrawing(
-          fullPatternSimulationJobs[0]->pMesh->getLabel());
-      SimulationResults fullModelResults =
-          simulator.executeSimulation(pFullPatternSimulationJob);
-      fullModelResults.registerForDrawing();
-      // Drawing of reduced pattern results
-      const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
-          reducedPatternSimulationJobs[simulationJobIndex];
-      SimulationResults reducedModelResults =
-          simulator.executeSimulation(pReducedPatternSimulationJob);
-      reducedModelResults.registerForDrawing();
-      polyscope::show();
-      // Save a screensh
-      //    const std::string screenshotFilename =
-      //        "/home/iason/Coding/Projects/Approximating shapes with flat "
-      //        "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
-      //        + m_pFullPatternSimulationMesh->getLabel() + "_" +
-      //        simulationScenarioStrings[simulationScenarioIndex];
-      //    polyscope::screenshot(screenshotFilename, false);
-      fullModelResults.unregister();
-      reducedModelResults.unregister();
-      //    double error = computeError(
-      //        reducedModelResults,
-      //        global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
-      //    std::cout << "Error of simulation scenario "
-      //              << simulationScenarioStrings[simulationScenarioIndex] << "
-      //              is "
-      //              << error << std::endl;
-    }
-  }
-#endif // POLYSCOPE_DEFINED
-  void saveMeshFiles() const {
-    const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
-    assert(numberOfSimulationJobs != 0 &&
-           fullPatternSimulationJobs.size() ==
-               reducedPatternSimulationJobs.size());
-
-    fullPatternSimulationJobs[0]->pMesh->savePly("FullPattern_undeformed.ply");
-    reducedPatternSimulationJobs[0]->pMesh->savePly(
-        "ReducedPattern_undeformed.ply");
-    FormFinder simulator;
-    for (int simulationJobIndex = 0;
-         simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
-      // Drawing of full pattern results
-      const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
-          fullPatternSimulationJobs[simulationJobIndex];
-      SimulationResults fullModelResults =
-          simulator.executeSimulation(pFullPatternSimulationJob);
-      fullModelResults.saveDeformedModel();
-
-      // Drawing of reduced pattern results
-      const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
-          reducedPatternSimulationJobs[simulationJobIndex];
-      SimulationResults reducedModelResults =
-          simulator.executeSimulation(pReducedPatternSimulationJob);
-      reducedModelResults.saveDeformedModel();
-    }
-  }
-
-  void
-  writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization,
-                csvFile &os) {
-    os << "Obj value";
-    for (const ReducedModelOptimizer::xRange &range :
-         settings_optimization.xRanges) {
-      os << range.label;
-    }
-    os << "Time(s)";
-    os << "#Crashes";
-  }
-
-  void
-  writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization,
-                 csvFile &os) const {
-    os << objectiveValue;
-    for (const double &optimalX : x) {
-      os << optimalX;
-    }
-
-    for (int unusedXVarCounter = 0;
-         unusedXVarCounter < settings_optimization.xRanges.size() - x.size();
-         unusedXVarCounter++) {
-      os << "-";
-    }
-
-    os << time;
-    if (numberOfSimulationCrashes == 0) {
-      os << "-";
-    } else {
-      os << numberOfSimulationCrashes;
-    }
-  }
-};
-
 #endif // REDUCEDMODELOPTIMIZER_HPP