Comparison over the full pattern test set for the three reduced models and exporting of the results of each optimization round.

CMakeLists.txt

@@ -24,8 +24,8 @@ download_project(PROJ               MATPLOTPLUSPLUS
                  )
 add_subdirectory(${MATPLOTPLUSPLUS_SOURCE_DIR})
 
-if (EXISTS ${MYSOURCES_SOURCE_DIR})
 set(MYSOURCES_SOURCE_DIR "/home/iason/Coding/Libraries/MySources")
+if (EXISTS ${MYSOURCES_SOURCE_DIR})
 else()
 ##MySources
 download_project(PROJ               MYSOURCES

src/main.cpp

@@ -30,13 +30,13 @@ int main(int argc, char *argv[]) {
   std::vector<vcg::Point2i> CWreducedModelEdges{vcg::Point2i(1, 5),
                                                 vcg::Point2i(3, 1)};
   FlatPattern CWReducedModel(numberOfNodesPerSlot, CWreducedModelEdges);
-  CWReducedModel.setLabel("CW_reduced");
+  CWReducedModel.setLabel("CCW_reduced");
   CWReducedModel.scale(0.03);
 
   std::vector<vcg::Point2i> CCWreducedModelEdges{vcg::Point2i(1, 5),
                                                  vcg::Point2i(3, 5)};
   FlatPattern CCWReducedModel(numberOfNodesPerSlot, CCWreducedModelEdges);
-  CCWReducedModel.setLabel("CCW_reduced");
+  CCWReducedModel.setLabel("CW_reduced");
   CCWReducedModel.scale(0.03);
 
   std::vector<FlatPattern *> reducedModels{&singleBarReducedModel,
@@ -48,7 +48,7 @@ int main(int argc, char *argv[]) {
   ReducedModelOptimizer::xRange beamE{"E", 0.1, 1.9};
   ReducedModelOptimizer::xRange innerHexagonSize{"HexagonSize", 0.1, 0.9};
   // Test set of full patterns
-  std::string fullPatternsTestSetDirectory = "../../TestSet";
+  std::string fullPatternsTestSetDirectory = "../TestSet";
   if (!std::filesystem::exists(
           std::filesystem::path(fullPatternsTestSetDirectory))) {
     std::cerr << "Full pattern directory does not exist: "
@@ -77,6 +77,7 @@ int main(int argc, char *argv[]) {
       pFullPattern->scale(0.03);
       FlatPattern *pReducedPattern = new FlatPattern();
       pReducedPattern->copy(*reducedModels[reducedPatternIndex]);
+      //      pReducedPattern->copy(*reducedModels[2]);
       patternPairs.push_back(std::make_pair(pFullPattern, pReducedPattern));
     }
   }
@@ -90,7 +91,7 @@ int main(int argc, char *argv[]) {
   //  for (settings_optimization.numberOfFunctionCalls = 100;
   //       settings_optimization.numberOfFunctionCalls < 5000;
   //       settings_optimization.numberOfFunctionCalls += 100) {
-  settings_optimization.numberOfFunctionCalls = 1000;
+  settings_optimization.numberOfFunctionCalls = 10;
   const std::string optimizationSettingsString =
       settings_optimization.toString();
   std::filesystem::path thisOptimizationDirectory(
@@ -129,47 +130,60 @@ int main(int argc, char *argv[]) {
   auto runtime_ms =
       std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
 
-  //  for (int patternPairIndex = 0; patternPairIndex < patternPairs.size();
+  for (int patternPairIndex = 0; patternPairIndex < patternPairs.size();
-  //       patternPairIndex++) {
+       patternPairIndex++) {
-  //    std::filesystem::path fullPatternPath(
+    std::filesystem::path saveToPath(
-  //        std::filesystem::path(thisOptimizationDirectory)
+        std::filesystem::path(thisOptimizationDirectory)
-  //            .append(patternPairs[patternPairIndex].first->getLabel()));
+            .append(patternPairs[patternPairIndex].first->getLabel() + "@" +
-  //    std::filesystem::create_directories(std::filesystem::path(fullPatternPath));
+                    patternPairs[patternPairIndex].second->getLabel()));
+    std::filesystem::create_directories(std::filesystem::path(saveToPath));
 
-  //    std::filesystem::path reducedPatternPath(fullPatternPath.append(
+    optimizationResults_testSet[patternPairIndex].save(saveToPath);
-  //        patternPairs[patternPairIndex].second->getLabel()));
+    //    optimizationResults_testSet[patternPairIndex].draw();
-  //    std::filesystem::create_directories(
+  }
-  //        std::filesystem::path(reducedPatternPath));
-
-  //    optimizationResults_testSet[patternPairIndex].save(reducedPatternPath);
-  //  }
   csvFile statistics(std::filesystem::path(thisOptimizationDirectory)
                          .append("statistics.csv")
                          .string(),
                      false);
 
-  for (const auto &range : settings_optimization.xRanges) {
+  // Write header to csv
-    statistics << range.min << range.max;
+  statistics << "FullPattern@ReducedPattern"
-  }
+             << "Obj value";
-  statistics << settings_optimization.numberOfFunctionCalls;
+  for (const ReducedModelOptimizer::xRange &range :
-  if (totalNumberOfSimulationCrashes == 0) {
+       settings_optimization.xRanges) {
-    statistics << "No crashes";
+    statistics << range.label;
-  } else {
-    statistics << totalNumberOfSimulationCrashes;
   }
+  statistics << "Time(s)";
+  statistics << "#Crashes";
+  statistics << endrow;
 
-  statistics << totalError;
+  // Write data
   for (int patternPairIndex = 0; patternPairIndex < patternPairs.size();
        patternPairIndex++) {
-    statistics << patternPairs[patternPairIndex].first->getLabel();
+    statistics << patternPairs[patternPairIndex].first->getLabel() + "@" +
+                      patternPairs[patternPairIndex].second->getLabel();
     statistics << optimizationResults_testSet[patternPairIndex].objectiveValue;
     for (const double &optimalX :
          optimizationResults_testSet[patternPairIndex].x) {
       statistics << optimalX;
     }
-  }
-  statistics << runtime_ms.count() / 1000.0;
 
+    for (int unusedXVarCounter = 0;
+         unusedXVarCounter <
+         settings_optimization.xRanges.size() -
+             optimizationResults_testSet[patternPairIndex].x.size();
+         unusedXVarCounter++) {
+      statistics << "-";
+    }
+
+    statistics << optimizationResults_testSet[patternPairIndex].time;
+    if (totalNumberOfSimulationCrashes == 0) {
+      statistics << "No crashes";
+    } else {
+      statistics << totalNumberOfSimulationCrashes;
+    }
+    statistics << endrow;
+  }
   statistics << endrow;
   //  }
 

src/reducedmodeloptimizer.cpp

@@ -10,7 +10,7 @@
 struct GlobalOptimizationVariables {
   std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
   std::vector<SimulationJob> g_fullPatternSimulationJob;
-  std::vector<std::shared_ptr<SimulationJob>> g_reducedPatternSimulationJob;
+  std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
   std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
       reducedToFullInterfaceViMap;
   std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
@@ -18,7 +18,7 @@ struct GlobalOptimizationVariables {
   matplot::line_handle gPlotHandle;
   std::vector<double> gObjectiveValueHistory;
   Eigen::Vector2d g_initialX;
-  std::unordered_set<size_t> g_reducedPatternExludedEdges;
+  std::unordered_set<size_t> reducedPatternExludedEdges;
   Eigen::VectorXd g_initialParameters;
   std::vector<ReducedModelOptimizer::SimulationScenario>
       simulationScenarioIndices;
@@ -127,7 +127,7 @@ double ReducedModelOptimizer::computeError(
     const Eigen::MatrixX3d &optimalReducedPatternDisplacements) {
   double error = 0;
   auto &global = tls[omp_get_thread_num()];
-  for (const auto reducedFullViPair : global.g_reducedToFullViMap) {
+  for (const auto reducedFullViPair : global.reducedToFullInterfaceViMap) {
     VertexIndex reducedModelVi = reducedFullViPair.first;
     //    const auto pos =
     //        g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP();
@@ -152,10 +152,43 @@ double ReducedModelOptimizer::computeError(
   return error;
 }
 
+double ReducedModelOptimizer::computeError(
+    const std::vector<Vector6d> &reducedPatternDisplacements,
+    const std::vector<Vector6d> &fullPatternDisplacements,
+    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+        &reducedToFullInterfaceViMap) {
+  double error = 0;
+  for (const auto reducedFullViPair : reducedToFullInterfaceViMap) {
+    VertexIndex reducedModelVi = reducedFullViPair.first;
+    //    const auto pos =
+    //        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[reducedModelVi][0],
+        reducedPatternDisplacements[reducedModelVi][1],
+        reducedPatternDisplacements[reducedModelVi][2]);
+    if (!std::isfinite(reducedVertexDisplacement[0]) ||
+        !std::isfinite(reducedVertexDisplacement[1]) ||
+        !std::isfinite(reducedVertexDisplacement[2])) {
+      std::terminate();
+    }
+    Eigen::Vector3d fullVertexDisplacement(
+        fullPatternDisplacements[reducedFullViPair.second][0],
+        fullPatternDisplacements[reducedFullViPair.second][1],
+        fullPatternDisplacements[reducedFullViPair.second][2]);
+    Eigen::Vector3d errorVector =
+        fullVertexDisplacement - reducedVertexDisplacement;
+    //    error += errorVector.squaredNorm();
+    error += errorVector.norm();
+  }
+  return error;
+}
+
 void updateMesh(long n, const double *x) {
   auto &global = tls[omp_get_thread_num()];
   std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
-      global.g_reducedPatternSimulationJob[global.simulationScenarioIndices[0]]
+      global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
           ->pMesh;
   //  const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0];
   //  std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << "
@@ -233,8 +266,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
   //    std::cout << "x[" + std::to_string(parameterIndex) + "]="
   //              << x[parameterIndex] << std::endl;
   //  }
-  const Element &e =
+  const Element &e = global.reducedPatternSimulationJobs[0]->pMesh->elements[0];
-      global.g_reducedPatternSimulationJob[0]->pMesh->elements[0];
   //  std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
   //            << e.firstBendingConstFactor << " " <<
   //            e.secondBendingConstFactor
@@ -252,7 +284,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
   //    simulationSettings.shouldDraw = true;
   for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
     SimulationResults reducedModelResults = simulator.executeSimulation(
-        global.g_reducedPatternSimulationJob[simulationScenarioIndex],
+        global.reducedPatternSimulationJobs[simulationScenarioIndex],
         simulationSettings);
     std::string filename;
     if (!reducedModelResults.converged /*&&
@@ -295,7 +327,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
     std::ofstream out(filename, std::ios_base::app);
     auto pMesh =
         global
-            .g_reducedPatternSimulationJob[global.simulationScenarioIndices[0]]
+            .reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
             ->pMesh;
 
     for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
@@ -314,6 +346,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
                "\n\n";
     out.close();
   }
+  //  std::cout << error << std::endl;
   if (error < global.minY) {
     global.minY = error;
     global.minX.assign(x, x + n);
@@ -341,9 +374,47 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
   return error;
 }
 
+void ReducedModelOptimizer::createSimulationMeshes(
+    FlatPattern &fullModel, FlatPattern &reducedModel,
+    std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
+    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+  if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) {
+    std::cerr << "Error: A rectangular cross section is expected." << std::endl;
+    terminate();
+  }
+  // Full pattern
+  pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
+  pFullPatternSimulationMesh->setBeamCrossSection(
+      CrossSectionType{0.002, 0.002});
+  pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
+
+  // Reduced pattern
+  pReducedPatternSimulationMesh =
+      std::make_shared<SimulationMesh>(reducedModel);
+  pReducedPatternSimulationMesh->setBeamCrossSection(
+      CrossSectionType{0.002, 0.002});
+  pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
+}
+
+void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel,
+                                                   FlatPattern &reducedModel) {
+  ReducedModelOptimizer::createSimulationMeshes(
+      fullModel, reducedModel, m_pFullPatternSimulationMesh,
+      m_pReducedPatternSimulationMesh);
+}
+
 void ReducedModelOptimizer::computeMaps(
+    const std::unordered_set<size_t> &reducedModelExcludedEdges,
+    const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
     FlatPattern &fullPattern, FlatPattern &reducedPattern,
-    const std::unordered_set<size_t> &reducedModelExcludedEges) {
+    std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+        &reducedToFullInterfaceViMap,
+    std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
+        &fullToReducedInterfaceViMap,
+    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() &&
@@ -368,20 +439,20 @@ void ReducedModelOptimizer::computeMaps(
   //  std::cout << "Dups in fan pair:" << duplicateVerticesPerFanPair <<
   //  std::endl;
 
-  // Save excluded edges
+  // Save excluded edges TODO:this changes the global object. Should this be a
-  auto &global = tls[omp_get_thread_num()];
+  // function parameter?
-  global.g_reducedPatternExludedEdges.clear();
+  global.reducedPatternExludedEdges.clear();
   const size_t fanSize = 6;
   const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN();
   for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
-    for (const size_t ei : reducedModelExcludedEges) {
+    for (const size_t ei : reducedModelExcludedEdges) {
-      global.g_reducedPatternExludedEdges.insert(
+      global.reducedPatternExludedEdges.insert(
           fanIndex * reducedBaseTriangleNumberOfEdges + ei);
     }
   }
 
   // Construct reduced->full and full->reduced interface vi map
-  global.reducedToFullInterfaceViMap.clear();
+  reducedToFullInterfaceViMap.clear();
   vcg::tri::Allocator<FlatPattern>::PointerUpdater<FlatPattern::VertexPointer>
       pu_reducedModel;
   reducedPattern.deleteDanglingVertices(pu_reducedModel);
@@ -391,36 +462,30 @@ void ReducedModelOptimizer::computeMaps(
       reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/;
   reducedPattern.createFan();
   for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
-    global.reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset *
+    reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex +
-                                           fanIndex +
+                                reducedModelBaseTriangleInterfaceVi] =
-                                       reducedModelBaseTriangleInterfaceVi] =
         fullModelBaseTriangleInterfaceVi +
         fanIndex * fullPatternInterfaceVertexOffset;
   }
-  m_fullToReducedInterfaceViMap.clear();
+  fullToReducedInterfaceViMap.clear();
-  constructInverseMap(global.reducedToFullInterfaceViMap,
+  constructInverseMap(reducedToFullInterfaceViMap, fullToReducedInterfaceViMap);
-                      m_fullToReducedInterfaceViMap);
 
-  //  fullPattern.setLabel("FullPattern");
-  //  reducedPattern.setLabel("ReducedPattern");
   // Create opposite vertex map
-  m_fullPatternOppositeInterfaceViMap.clear();
+  fullPatternOppositeInterfaceViMap.clear();
   for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) {
     const size_t vi0 = fullModelBaseTriangleInterfaceVi +
                        fanIndex * fullPatternInterfaceVertexOffset;
     const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset;
     assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN());
-    m_fullPatternOppositeInterfaceViMap[vi0] = vi1;
+    fullPatternOppositeInterfaceViMap[vi0] = vi1;
   }
 
-  global.g_reducedToFullViMap = global.reducedToFullInterfaceViMap;
-
   const bool debugMapping = false;
   if (debugMapping) {
     reducedPattern.registerForDrawing();
     std::vector<glm::vec3> colors_reducedPatternExcludedEdges(
         reducedPattern.EN(), glm::vec3(0, 0, 0));
-    for (const size_t ei : global.g_reducedPatternExludedEdges) {
+    for (const size_t ei : global.reducedPatternExludedEdges) {
       colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0);
     }
     const std::string label = reducedPattern.getLabel();
@@ -433,7 +498,7 @@ void ReducedModelOptimizer::computeMaps(
     std::vector<glm::vec3> nodeColorsOpposite(fullPattern.VN(),
                                               glm::vec3(0, 0, 0));
     for (const std::pair<size_t, size_t> oppositeVerts :
-         m_fullPatternOppositeInterfaceViMap) {
+         fullPatternOppositeInterfaceViMap) {
       auto color = polyscope::getNextUniqueColor();
       nodeColorsOpposite[oppositeVerts.first] = color;
       nodeColorsOpposite[oppositeVerts.second] = color;
@@ -469,21 +534,14 @@ void ReducedModelOptimizer::computeMaps(
   }
 }
 
-void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel,
+void ReducedModelOptimizer::computeMaps(
-                                                   FlatPattern &reducedModel) {
+    FlatPattern &fullPattern, FlatPattern &reducedPattern,
-  if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) {
+    const std::unordered_set<size_t> &reducedModelExcludedEdges) {
-    std::cerr << "Error: A rectangular cross section is expected." << std::endl;
+  auto &global = tls[omp_get_thread_num()];
-    terminate();
+  ReducedModelOptimizer::computeMaps(
-  }
+      reducedModelExcludedEdges, slotToNode, fullPattern, reducedPattern,
-  m_pReducedPatternSimulationMesh =
+      global.reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap,
-      std::make_shared<SimulationMesh>(reducedModel);
+      m_fullPatternOppositeInterfaceViMap);
-  m_pReducedPatternSimulationMesh->setBeamCrossSection(
-      CrossSectionType{0.002, 0.002});
-  m_pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
-  m_pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
-  m_pFullPatternSimulationMesh->setBeamCrossSection(
-      CrossSectionType{0.002, 0.002});
-  m_pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
 }
 
 ReducedModelOptimizer::ReducedModelOptimizer(
@@ -636,8 +694,6 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
 
   global.gObjectiveValueHistory.clear();
 
-  const size_t n = global.g_initialParameters.rows();
-  assert(n == 3);
   // g_optimizeInnerHexagonSize ? 5: 4;
   //  const size_t npt = (n + 1) * (n + 2) / 2;
   //  //      ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
@@ -701,8 +757,13 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
         std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
   }
   auto end = std::chrono::system_clock::now();
-  auto elapsed = std::chrono::duration_cast<std::chrono::seconds>(end - start);
+  auto elapsed =
-  Results results{global.numOfSimulationCrashes, global.minX, global.minY};
+      std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
+  Results results;
+  results.numberOfSimulationCrashes = global.numOfSimulationCrashes;
+  results.x = global.minX;
+  results.objectiveValue = global.minY;
+  results.time = elapsed.count() / 1000.0;
   std::cout << "Finished optimizing." << endl;
   //  std::cout << "Solution x:" << endl;
   //  std::cout << result.x << endl;
@@ -1000,28 +1061,36 @@ ReducedModelOptimizer::createScenarios(
 void ReducedModelOptimizer::visualizeResults(
     const std::vector<std::shared_ptr<SimulationJob>>
         &fullPatternSimulationJobs,
-    const std::vector<SimulationScenario> &simulationScenarios) {
+    const std::vector<std::shared_ptr<SimulationJob>>
-  m_pFullPatternSimulationMesh->registerForDrawing();
+        &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(
-        m_pFullPatternSimulationMesh->getLabel());
+        pFullPatternSimulationMesh->getLabel());
     SimulationResults fullModelResults =
         simulator.executeSimulation(pFullPatternSimulationJob);
     fullModelResults.registerForDrawing();
     fullModelResults.saveDeformedModel();
-    auto &global = tls[omp_get_thread_num()];
     const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
-        global.g_reducedPatternSimulationJob[simulationScenarioIndex];
+        reducedPatternSimulationJobs[simulationScenarioIndex];
     SimulationResults reducedModelResults =
         simulator.executeSimulation(pReducedPatternSimulationJob);
-    double error = computeError(
+    double error = computeError(reducedModelResults.displacements,
-        reducedModelResults,
+                                fullModelResults.displacements,
-        global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
+                                reducedToFullInterfaceViMap);
     std::cout << "Error of simulation scenario "
               << simulationScenarioStrings[simulationScenarioIndex] << " is "
               << error << std::endl;
+    totalError += error;
     reducedModelResults.registerForDrawing();
     //    firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
     //    reducedModelResults.saveDeformedModel();
@@ -1029,8 +1098,8 @@ void ReducedModelOptimizer::visualizeResults(
     const std::string screenshotFilename =
         "/home/iason/Coding/Projects/Approximating shapes with flat "
         "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
-        "Images" +
+        "Images/" +
-        m_pFullPatternSimulationMesh->getLabel() + "_" +
+        pFullPatternSimulationMesh->getLabel() + "_" +
         simulationScenarioStrings[simulationScenarioIndex];
     polyscope::show();
     polyscope::screenshot(screenshotFilename, false);
@@ -1038,6 +1107,7 @@ void ReducedModelOptimizer::visualizeResults(
     reducedModelResults.unregister();
     //    firstOptimizationRoundResults[simulationScenarioIndex].unregister();
   }
+  std::cout << "Total error:" << totalError << std::endl;
 }
 
 ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
@@ -1056,7 +1126,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
   std::vector<std::shared_ptr<SimulationJob>> simulationJobs =
       createScenarios(m_pFullPatternSimulationMesh);
   global.g_optimalReducedModelDisplacements.resize(6);
-  global.g_reducedPatternSimulationJob.resize(6);
+  global.reducedPatternSimulationJobs.resize(6);
   global.g_firstRoundIterationIndex = 0;
   global.minY = std::numeric_limits<double>::max();
   global.numOfSimulationCrashes = 0;
@@ -1072,18 +1142,28 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
     global.g_optimalReducedModelDisplacements[simulationScenarioIndex].resize(
         m_pReducedPatternSimulationMesh->VN(), 3);
     computeDesiredReducedModelDisplacements(
-        fullModelResults, global.g_reducedToFullViMap,
+        fullModelResults, global.reducedToFullInterfaceViMap,
         global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
     SimulationJob reducedPatternSimulationJob;
     reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
     computeReducedModelSimulationJob(*pFullPatternSimulationJob,
                                      m_fullToReducedInterfaceViMap,
                                      reducedPatternSimulationJob);
-    global.g_reducedPatternSimulationJob[simulationScenarioIndex] =
+    global.reducedPatternSimulationJobs[simulationScenarioIndex] =
         std::make_shared<SimulationJob>(reducedPatternSimulationJob);
   }
   Results optResults = runOptimization(xRanges, &objective);
-  updateMesh(optResults.x.size(), optResults.x.data());
+  for (int simulationScenarioIndex : global.simulationScenarioIndices) {
-  visualizeResults(simulationJobs, global.simulationScenarioIndices);
+    optResults.fullPatternSimulationJobs.push_back(
+        simulationJobs[simulationScenarioIndex]);
+    optResults.reducedPatternSimulationJobs.push_back(
+        global.reducedPatternSimulationJobs[simulationScenarioIndex]);
+  }
+  //  updateMesh(optResults.x.size(), optResults.x.data());
+  //  optResults.draw();
+
+  //  visualizeResults(simulationJobs, global.simulationScenarioIndices);
+  //  visualizeResults(simulationJobs, global.g_reducedPatternSimulationJob,
+  //                   global.g_simulationScenarioIndices);
   return optResults;
 }

src/reducedmodeloptimizer.hpp

@@ -31,11 +31,7 @@ public:
     Saddle,
     NumberOfSimulationScenarios
   };
-  struct Results {
+
-    int numberOfSimulationCrashes{0};
-    std::vector<double> x;
-    double objectiveValue;
-  };
   struct xRange {
     std::string label;
     double min;
@@ -45,6 +41,7 @@ public:
              "]";
     }
   };
+  struct Results;
 
   struct Settings {
     std::vector<xRange> xRanges;
@@ -120,6 +117,40 @@ public:
   static double objective(double x0, double x1, double x2, double x3);
   static double objective(double b, double h, double E);
 
+  static std::vector<std::shared_ptr<SimulationJob>>
+  createScenarios(const std::shared_ptr<SimulationMesh> &pMesh,
+                  const std::unordered_map<size_t, size_t>
+                      &fullPatternOppositeInterfaceViMap);
+  static void createSimulationMeshes(
+      FlatPattern &fullModel, FlatPattern &reducedModel,
+      std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
+      std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh);
+  static void computeMaps(
+      const std::unordered_set<size_t> &reducedModelExcludedEdges,
+      const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
+      FlatPattern &fullPattern, FlatPattern &reducedPattern,
+      std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+          &reducedToFullInterfaceViMap,
+      std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
+          &fullToReducedInterfaceViMap,
+      std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
+          &fullPatternOppositeInterfaceViMap);
+  static void
+  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);
+
+  static double computeError(const std::vector<Vector6d> &reducedPatternResults,
+                             const std::vector<Vector6d> &fullPatternResults,
+                             const std::unordered_map<ReducedPatternVertexIndex,
+                                                      FullPatternVertexIndex>
+                                 &reducedToFullInterfaceViMap);
+
 private:
   void
   visualizeResults(const std::vector<std::shared_ptr<SimulationJob>>
@@ -147,5 +178,116 @@ private:
   static double objective(long n, const double *x);
   FormFinder simulator;
 };
+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 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));
+
+    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("FullPattern");
+      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("ReducedPattern");
+      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();
+      // Load reduced pattern files
+      for (const auto &fileEntry : filesystem::directory_iterator(
+               std::filesystem::path(simulationScenarioPath)
+                   .append("FullPattern"))) {
+        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("ReducedPattern"))) {
+        const auto filepath = fileEntry.path();
+        if (filepath.extension() == ".json") {
+          SimulationJob job;
+          job.load(filepath.string());
+          reducedPatternSimulationJobs.push_back(
+              std::make_shared<SimulationJob>(job));
+        }
+      }
+    }
+  }
+};
 
 #endif // REDUCEDMODELOPTIMIZER_HPP

src/reducedmodeloptimizer.hpp.autosave

@@ -1,151 +0,0 @@
-#ifndef REDUCEDMODELOPTIMIZER_HPP
-#define REDUCEDMODELOPTIMIZER_HPP
-
-#include "beamformfinder.hpp"
-#include "csvfile.hpp"
-#include "edgemesh.hpp"
-#include "elementalmesh.hpp"
-#include "matplot/matplot.h"
-#include <Eigen/Dense>
-
-using FullPatternVertexIndex = VertexIndex;
-using ReducedPatternVertexIndex = VertexIndex;
-
-class ReducedModelOptimizer {
-  std::shared_ptr<SimulationMesh> m_pReducedPatternSimulationMesh;
-  std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh;
-  std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
-      m_fullToReducedInterfaceViMap;
-  std::unordered_map<FullPatternVertexIndex, FullPatternVertexIndex>
-      m_fullPatternOppositeInterfaceViMap;
-  std::unordered_map<size_t, size_t> nodeToSlot;
-  std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
-  std::vector<double> initialGuess;
-
-public:
-  enum SimulationScenario {
-    Axial,
-    Shear,
-    Bending,
-    Dome,
-    Saddle,
-    NumberOfSimulationScenarios
-  };
-  struct Results {
-    int numberOfSimulationCrashes{0};
-    std::vector<double> x;
-    double objectiveValue;
-  };
-  struct xRange {
-    std::string label;
-    double min;
-    double max;
-    std::string toString() const {
-      return label + "=[" + std::to_string(min) + "," + std::to_string(max) +
-             "]";
-    }
-  };
-
-  struct Settings {
-    std::vector<xRange> xRanges;
-    int numberOfFunctionCalls{100};
-    double solutionAccuracy{1e-5};
-
-    std::string toString() const {
-      std::string settingsString;
-      if (!xRanges.empty()) {
-        std::string xRangesString;
-        for (const xRange &range : xRanges) {
-          xRangesString += range.toString() + " ";
-        }
-        settingsString += xRangesString;
-      }
-      settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls) +
-                        " Accuracy=" + std::to_string(solutionAccuracy);
-
-      return settingsString;
-    }
-
-    void writeTo(csvFile &csv) const {
-      // Create settings csv header
-      if (!xRanges.empty()) {
-        for (const xRange &range : xRanges) {
-          csv << range.label + " min";
-          csv << range.label + " max";
-        }
-      }
-      csv << "Function Calls";
-      csv << "Solution Accuracy";
-      csv << endrow;
-
-      if (!xRanges.empty()) {
-        for (const xRange &range : xRanges) {
-          csv << range.min;
-          csv << range.max;
-        }
-      }
-      csv << numberOfFunctionCalls;
-      csv << solutionAccuracy;
-      csv << endrow;
-    }
-  };
-
-  inline static const std::string simulationScenarioStrings[] = {
-      "Axial", "Shear", "Bending", "Double", "Saddle"};
-  Results optimize(const Settings &xRanges,
-                   const std::vector<SimulationScenario> &simulationScenarios =
-                       std::vector<SimulationScenario>());
-  double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;
-
-  ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot);
-  static void computeReducedModelSimulationJob(
-      const SimulationJob &simulationJobOfFullModel,
-      const std::unordered_map<size_t, size_t> &simulationJobFullToReducedMap,
-      SimulationJob &simulationJobOfReducedModel);
-
-  SimulationJob
-  getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
-
-  void initializePatterns(
-      FlatPattern &fullPattern, FlatPattern &reducedPatterm,
-      const std::unordered_set<size_t> &reducedModelExcludedEges);
-
-  void setInitialGuess(std::vector<double> v);
-
-  static void runBeamOptimization();
-
-  static void runSimulation(const std::string &filename,
-                            std::vector<double> &x);
-
-  static double objective(double x0, double x1, double x2, double x3);
-  static double objective(double b, double h, double E);
-
-private:
-  void
-  visualizeResults(const std::vector<std::shared_ptr<SimulationJob>>
-                       &fullPatternSimulationJobs,
-                   const std::vector<SimulationScenario> &simulationScenarios);
-  static void computeDesiredReducedModelDisplacements(
-      const SimulationResults &fullModelResults,
-      const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
-      Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
-  static Results runOptimization(const Settings &settings,
-                                 double (*pObjectiveFunction)(long,
-                                                              const double *));
-  std::vector<std::shared_ptr<SimulationJob>>
-  createScenarios(const std::shared_ptr<SimulationMesh> &pMesh);
-  void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern,
-                   const std::unordered_set<size_t> &reducedModelExcludedEges);
-  void createSimulationMeshes(FlatPattern &fullModel,
-                              FlatPattern &reducedModel);
-  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;
-};
-
-#endif // REDUCEDMODELOPTIMIZER_HPP