Using POLYSCOPE_DEFINED macro for enabling and disabling visualization related code.

src/main.cpp

@@ -1,11 +1,7 @@
 #include "beamformfinder.hpp"
 #include "csvfile.hpp"
 #include "edgemesh.hpp"
-#include "externvariables.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"
@@ -17,30 +13,23 @@
 #include <string>
 #include <vcg/complex/algorithms/update/position.h>
 
-bool printDebugInfo;
-
 int main(int argc, char *argv[]) {
-  if (argc < 5) {
+  if (argc < 3) {
-    std::cerr << "Specify at least if D(ebug) or R(elease) and the pattern "
+    std::cerr << "Specify at least the two pattern filepaths to be "
-                 "pair to be optimized."
+                 "optimized.Exiting.."
               << std::endl;
     std::terminate();
   }
 
-  if (argv[1] == "D") {
-    printDebugInfo = true;
-  } else {
-    printDebugInfo = false;
-  }
   // Populate the pattern pair to be optimized
   ////Full pattern
-  const std::string filepath_fullPattern = argv[2];
+  const std::string filepath_fullPattern = argv[1];
   FlatPattern fullPattern(filepath_fullPattern);
   fullPattern.setLabel(
       std::filesystem::path(filepath_fullPattern).stem().string());
   fullPattern.scale(0.03);
   ////Reduced pattern
-  const std::string filepath_reducedPattern = argv[3];
+  const std::string filepath_reducedPattern = argv[2];
   FlatPattern reducedPattern(filepath_reducedPattern);
   reducedPattern.setLabel(
       std::filesystem::path(filepath_reducedPattern).stem().string());
@@ -56,13 +45,11 @@ int main(int argc, char *argv[]) {
                                    innerHexagonSize};
   const bool input_numberOfFunctionCallsDefined = argc >= 5;
   settings_optimization.numberOfFunctionCalls =
-      input_numberOfFunctionCallsDefined ? std::atoi(argv[4]) : 100;
+      input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100;
   // Optimize pair
   const std::string pairName =
       fullPattern.getLabel() + "@" + reducedPattern.getLabel();
-  if (printDebugInfo) {
+  std::cout << "Optimizing " << pairName << std::endl;
-    std::cout << "Optimizing " << pairName << std::endl;
-  }
 
   const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
   ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
@@ -72,7 +59,7 @@ int main(int argc, char *argv[]) {
   // Export results
   const bool input_resultDirectoryDefined = argc >= 6;
   std::string optimiziationResultsDirectory =
-      input_resultDirectoryDefined ? argv[5] : "OptimizationResults";
+      input_resultDirectoryDefined ? argv[4] : "OptimizationResults";
   std::filesystem::path dirPath_thisOptimization(
       std::filesystem::path(optimiziationResultsDirectory).append(pairName));
   std::filesystem::create_directories(dirPath_thisOptimization);

src/reducedmodeloptimizer.cpp

@@ -10,7 +10,7 @@ struct GlobalOptimizationVariables {
   std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
   std::vector<std::vector<Vector6d>> fullPatternDisplacements;
   std::vector<double> fullPatternDisplacementNormSum;
-  std::vector<SimulationJob> g_fullPatternSimulationJob;
+  std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
   std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
   std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
       reducedToFullInterfaceViMap;
@@ -33,94 +33,7 @@ struct GlobalOptimizationVariables {
   int numberOfFunctionCalls{0};
   int numberOfOptimizationParameters{3};
   ReducedModelOptimizer::Settings optimizationSettings;
-};
+} global;
-
-// static GlobalOptimizationVariables global;
-
-const static int MAX_THREAD = 64;
-#if defined(_MSC_VER)
-__declspec(align(64)) GlobalOptimizationVariables tls[MAX_THREAD];
-#elif defined(__GNUC__)
-GlobalOptimizationVariables tls[MAX_THREAD] __attribute__((aligned(64)));
-#endif
-//#pragma omp threadprivate(global)
-
-// struct OptimizationCallback {
-//  double operator()(const size_t &iterations, const Eigen::VectorXd &x,
-//                    const double &fval, Eigen::VectorXd &gradient) const {
-//    // run simulation
-//    //    SimulationResults reducedModelResults =
-//    //        simulator.executeSimulation(reducedModelSimulationJob);
-//    //    reducedModelResults.draw(reducedModelSimulationJob);
-//    gObjectiveValueHistory.push_back(fval);
-//    auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
-//                                   gObjectiveValueHistory.size());
-//    gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory);
-//    //    const std::string plotImageFilename = "objectivePlot.png";
-//    //    matplot::save(plotImageFilename);
-//    //    if (numberOfOptimizationRounds % 30 == 0) {
-//    //      std::filesystem::copy_file(
-//    //          std::filesystem::path(plotImageFilename),
-//    //          std::filesystem::path("objectivePlot_copy.png"));
-//    //    }
-//    //    std::stringstream ss;
-//    //    ss << x;
-//    //    reducedModelResults.simulationLabel = ss.str();
-//    //    SimulationResultsReporter resultsReporter;
-//    //    resultsReporter.reportResults(
-//    //        {reducedModelResults},
-//    //        std::filesystem::current_path().append("Results"));
-
-//    return true;
-//  }
-//};
-
-// struct Objective {
-//  double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const {
-//    assert(x.rows() == 4);
-
-//    //      drawSimulationJob(simulationJob);
-//    // Set mesh from x
-//    std::shared_ptr<SimulationMesh> reducedModel =
-//        g_reducedPatternSimulationJob.mesh;
-//    for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) {
-//      if (g_reducedPatternExludedEdges.contains(ei)) {
-//        continue;
-//      }
-//      Element &e = reducedModel->elements[ei];
-//      e.axialConstFactor = g_initialStiffnessFactors(ei, 0) * x(0);
-//      e.torsionConstFactor = g_initialStiffnessFactors(ei, 1) * x(1);
-//      e.firstBendingConstFactor = g_initialStiffnessFactors(ei, 2) * x(2);
-//      e.secondBendingConstFactor = g_initialStiffnessFactors(ei, 3) * x(3);
-//    }
-//    // run simulation
-//    SimulationResults reducedModelResults =
-//        simulator.executeSimulation(g_reducedPatternSimulationJob);
-//    //    std::stringstream ss;
-//    //    ss << x;
-//    //    reducedModelResults.simulationLabel = ss.str();
-//    //    SimulationResultsReporter resultsReporter;
-//    //    resultsReporter.reportResults(
-//    //        {reducedModelResults},
-//    //        std::filesystem::current_path().append("Results"));
-//    // compute error and return it
-//    double error = 0;
-//    for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) {
-//      VertexIndex reducedModelVi = reducedFullViPair.first;
-//      Eigen::Vector3d vertexDisplacement(
-//          reducedModelResults.displacements[reducedModelVi][0],
-//          reducedModelResults.displacements[reducedModelVi][1],
-//          reducedModelResults.displacements[reducedModelVi][2]);
-//      Eigen::Vector3d errorVector =
-//          Eigen::Vector3d(
-//              g_optimalReducedModelDisplacements.row(reducedModelVi)) -
-//          vertexDisplacement;
-//      error += errorVector.norm();
-//    }
-
-//    return error;
-//  }
-//};
 
 double ReducedModelOptimizer::computeError(
     const std::vector<Vector6d> &reducedPatternDisplacements,
@@ -128,7 +41,6 @@ double ReducedModelOptimizer::computeError(
     const double &interfaceDisplacementsNormSum,
     const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
         &reducedToFullInterfaceViMap) {
-  auto &global = tls[omp_get_thread_num()];
   double error = 0;
   for (const auto reducedFullViPair : reducedToFullInterfaceViMap) {
     VertexIndex reducedModelVi = reducedFullViPair.first;
@@ -162,7 +74,6 @@ double ReducedModelOptimizer::computeError(
 }
 
 void updateMesh(long n, const double *x) {
-  auto &global = tls[omp_get_thread_num()];
   std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
       global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
           ->pMesh;
@@ -215,12 +126,10 @@ void updateMesh(long n, const double *x) {
       pReducedPatternSimulationMesh->vert[vi].P() =
           global.g_innerHexagonVectors[vi / 2] * x[n - 1];
     }
-
-    pReducedPatternSimulationMesh->reset();
-    pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
-    //    pReducedPatternSimulationMesh->registerForDrawing("Optimized
-    //    hexagon"); polyscope::show();
   }
+
+  pReducedPatternSimulationMesh->reset();
+  pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
 }
 
 double ReducedModelOptimizer::objective(double b, double h, double E) {
@@ -235,7 +144,6 @@ double ReducedModelOptimizer::objective(double b, double h, double E,
 }
 
 double ReducedModelOptimizer::objective(long n, const double *x) {
-  auto &global = tls[omp_get_thread_num()];
   //  std::cout.precision(17);
 
   //  for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
@@ -257,40 +165,15 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
   double error = 0;
   FormFinder simulator;
   FormFinder::Settings simulationSettings;
-  //    simulationSettings.shouldDraw = true;
+  //      simulationSettings.shouldDraw = true;
   for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
     SimulationResults reducedModelResults = simulator.executeSimulation(
         global.reducedPatternSimulationJobs[simulationScenarioIndex],
         simulationSettings);
     std::string filename;
-    if (!reducedModelResults.converged /*&&
+    if (!reducedModelResults.converged) {
-        g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]
-                    ->pMesh->elements[0]
-                    .A > 1e-8 &
-            x[0] / x[1] < 60*/) {
       std::cout << "Failed simulation" << std::endl;
       error += std::numeric_limits<double>::max();
-      filename = "/home/iason/Coding/Projects/Approximating shapes with flat "
-                 "patterns/RodModelOptimizationForPatterns/build/"
-                 "ProblematicSimulationJobs/nonConv_dimensions.txt";
-      //      if (failedSimulationsXRatio.empty()) {
-      //        failedSimulationsXRatio.resize(2);
-      //      }
-      //      failedSimulationsXRatio[0].push_back(std::log(x[0] / x[1]));
-      //      failedSimulationsXRatio[1].push_back(
-      //          std::log(g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]
-      //                       ->pMesh->elements[0]
-      //                       .A));
-
-      //      SimulationResultsReporter::createPlot(
-      //          "log(b/h)", "log(A)", failedSimulationsXRatio[0],
-      //          failedSimulationsXRatio[1], "ratioToAPlot.png");
-      //      std::cout << "Failed simulation" << std::endl;
-      //      simulationSettings.shouldDraw = true;
-      //      simulationSettings.debugMessages = true;
-      //      simulator.executeSimulation(
-      //          g_reducedPatternSimulationJob[simulationScenarioIndex],
-      //          simulationSettings);
       global.numOfSimulationCrashes++;
     } else {
       error += computeError(
@@ -298,9 +181,6 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
           global.fullPatternDisplacements[simulationScenarioIndex],
           global.fullPatternDisplacementNormSum[simulationScenarioIndex],
           global.reducedToFullInterfaceViMap);
-      filename = "/home/iason/Coding/Projects/Approximating shapes with flat "
-                 "patterns/RodModelOptimizationForPatterns/build/"
-                 "ProblematicSimulationJobs/conv_dimensions.txt";
     }
     std::ofstream out(filename, std::ios_base::app);
     auto pMesh =
@@ -323,6 +203,13 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
         << " \nscenario:" + simulationScenarioStrings[simulationScenarioIndex] +
                "\n\n";
     out.close();
+#ifdef POLYSCOPE_DEFINED
+    ReducedModelOptimizer::visualizeResults(
+        global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs,
+        std::vector<SimulationScenario>{
+            static_cast<SimulationScenario>(simulationScenarioIndex)},
+        global.reducedToFullInterfaceViMap);
+#endif // POLYSCOPE_DEFINED
   }
   //  std::cout << error << std::endl;
   if (error < global.minY) {
@@ -391,7 +278,6 @@ void ReducedModelOptimizer::computeMaps(
     std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
         &fullPatternOppositeInterfaceViMap) {
 
-  auto &global = tls[omp_get_thread_num()];
   // Compute the offset between the interface nodes
   const size_t interfaceSlotIndex = 4; // bottom edge
   assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() &&
@@ -459,6 +345,7 @@ void ReducedModelOptimizer::computeMaps(
 
   const bool debugMapping = false;
   if (debugMapping) {
+#if POLYSCOPE_DEFINED
     reducedPattern.registerForDrawing();
     std::vector<glm::vec3> colors_reducedPatternExcludedEdges(
         reducedPattern.EN(), glm::vec3(0, 0, 0));
@@ -508,13 +395,13 @@ void ReducedModelOptimizer::computeMaps(
                                nodeColorsReducedToFull_full)
         ->setEnabled(true);
     polyscope::show();
+#endif
   }
 }
 
 void ReducedModelOptimizer::computeMaps(
     FlatPattern &fullPattern, FlatPattern &reducedPattern,
     const std::unordered_set<size_t> &reducedModelExcludedEdges) {
-  auto &global = tls[omp_get_thread_num()];
   ReducedModelOptimizer::computeMaps(
       reducedModelExcludedEdges, slotToNode, fullPattern, reducedPattern,
       global.reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap,
@@ -534,13 +421,14 @@ void ReducedModelOptimizer::initializePatterns(
   // reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel());
   assert(fullPattern.VN() == reducedPattern.VN() &&
          fullPattern.EN() >= reducedPattern.EN());
+#if POLYSCOPE_DEFINED
   polyscope::removeAllStructures();
+#endif
   // Create copies of the input models
   FlatPattern copyFullPattern;
   FlatPattern copyReducedPattern;
   copyFullPattern.copy(fullPattern);
   copyReducedPattern.copy(reducedPattern);
-  auto &global = tls[omp_get_thread_num()];
   global.optimizeInnerHexagonSize = copyReducedPattern.EN() == 2;
   if (global.optimizeInnerHexagonSize) {
     const double h = copyReducedPattern.getBaseTriangleHeight();
@@ -571,7 +459,6 @@ void ReducedModelOptimizer::initializePatterns(
 
 void ReducedModelOptimizer::initializeOptimizationParameters(
     const std::shared_ptr<SimulationMesh> &mesh) {
-  auto &global = tls[omp_get_thread_num()];
   global.numberOfOptimizationParameters = 3;
   global.g_initialParameters.resize(
       global.optimizeInnerHexagonSize ? ++global.numberOfOptimizationParameters
@@ -642,6 +529,74 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob(
   //      reducedModelNodalForcedNormals;
 }
 
+#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 interfaceDisplacementNormSum = 0;
+    for (const auto &interfaceViPair : reducedToFullInterfaceViMap) {
+      const int fullPatternInterfaceIndex = interfaceViPair.second;
+      Eigen::Vector3d fullPatternDisplacementVector(
+          fullModelResults.displacements[fullPatternInterfaceIndex][0],
+          fullModelResults.displacements[fullPatternInterfaceIndex][1],
+          fullModelResults.displacements[fullPatternInterfaceIndex][2]);
+      interfaceDisplacementNormSum += fullPatternDisplacementVector.norm();
+    }
+    reducedModelResults.saveDeformedModel();
+    fullModelResults.saveDeformedModel();
+    double error = computeError(
+        reducedModelResults.displacements, fullModelResults.displacements,
+        interfaceDisplacementNormSum, reducedToFullInterfaceViMap);
+    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;
+}
+#endif // POLYSCOPE_DEFINED
+
 void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
     const SimulationResults &fullModelResults,
     const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
@@ -666,7 +621,6 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
 
 ReducedModelOptimizer::Results
 ReducedModelOptimizer::runOptimization(const Settings &settings) {
-  auto &global = tls[omp_get_thread_num()];
 
   global.gObjectiveValueHistory.clear();
 
@@ -739,12 +693,10 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
   results.x = global.minX;
   results.objectiveValue = global.minY;
   results.time = elapsed.count() / 1000.0;
-  if (printDebugInfo) {
+  std::cout << "Finished optimizing." << endl;
-    std::cout << "Finished optimizing." << endl;
+  //  std::cout << "Solution x:" << endl;
-    //  std::cout << "Solution x:" << endl;
+  //  std::cout << result.x << endl;
-    //  std::cout << result.x << endl;
+  std::cout << "Objective value:" << global.minY << 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;
@@ -1035,71 +987,9 @@ ReducedModelOptimizer::createScenarios(
 //  reducedModelResults.unregister();
 //}
 
-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();
-  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 interfaceDisplacementNormSum = 0;
-    for (const auto &interfaceViPair : reducedToFullInterfaceViMap) {
-      const int fullPatternInterfaceIndex = interfaceViPair.second;
-      Eigen::Vector3d fullPatternDisplacementVector(
-          fullModelResults.displacements[fullPatternInterfaceIndex][0],
-          fullModelResults.displacements[fullPatternInterfaceIndex][1],
-          fullModelResults.displacements[fullPatternInterfaceIndex][2]);
-      interfaceDisplacementNormSum += fullPatternDisplacementVector.norm();
-    }
-    double error = computeError(
-        reducedModelResults.displacements, fullModelResults.displacements,
-        interfaceDisplacementNormSum, reducedToFullInterfaceViMap);
-    std::cout << "Error of simulation scenario "
-              << simulationScenarioStrings[simulationScenarioIndex] << " is "
-              << error << std::endl;
-    totalError += error;
-    reducedModelResults.registerForDrawing();
-    //    firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
-    //    reducedModelResults.saveDeformedModel();
-    //    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;
-}
-
 ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
     const Settings &optimizationSettings,
     const std::vector<SimulationScenario> &simulationScenarios) {
-  auto &global = tls[omp_get_thread_num()];
 
   global.simulationScenarioIndices = simulationScenarios;
   if (global.simulationScenarioIndices.empty()) {
@@ -1109,8 +999,6 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
         SimulationScenario::Saddle};
   }
 
-  std::vector<std::shared_ptr<SimulationJob>> simulationJobs =
-      createScenarios(m_pFullPatternSimulationMesh);
   global.g_optimalReducedModelDisplacements.resize(6);
   global.reducedPatternSimulationJobs.resize(6);
   global.fullPatternDisplacements.resize(6);
@@ -1120,13 +1008,15 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
   global.numOfSimulationCrashes = 0;
   global.numberOfFunctionCalls = 0;
   global.optimizationSettings = optimizationSettings;
+  global.fullPatternSimulationJobs =
+      createScenarios(m_pFullPatternSimulationMesh);
   //  polyscope::removeAllStructures();
 
   FormFinder::Settings settings;
   //  settings.shouldDraw = true;
   for (int simulationScenarioIndex : global.simulationScenarioIndices) {
     const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
-        simulationJobs[simulationScenarioIndex];
+        global.fullPatternSimulationJobs[simulationScenarioIndex];
     SimulationResults fullModelResults =
         simulator.executeSimulation(pFullPatternSimulationJob, settings);
     global.fullPatternDisplacements[simulationScenarioIndex] =
@@ -1159,7 +1049,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
   Results optResults = runOptimization(optimizationSettings);
   for (int simulationScenarioIndex : global.simulationScenarioIndices) {
     optResults.fullPatternSimulationJobs.push_back(
-        simulationJobs[simulationScenarioIndex]);
+        global.fullPatternSimulationJobs[simulationScenarioIndex]);
     optResults.reducedPatternSimulationJobs.push_back(
         global.reducedPatternSimulationJobs[simulationScenarioIndex]);
   }

src/reducedmodeloptimizer.hpp

@@ -191,50 +191,6 @@ struct ReducedModelOptimizer::Results {
   std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
   std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
 
-  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;
-    }
-  }
   void save(const string &saveToPath) const {
     assert(std::filesystem::is_directory(saveToPath));
 
@@ -296,7 +252,52 @@ struct ReducedModelOptimizer::Results {
       }
     }
   }
+#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 &&