Exporting simulation scenario when number of simulations is greater than some treshold.

src/main.cpp

@@ -38,10 +38,19 @@
 //}
 
 int main(int argc, char *argv[]) {
+  //  ReducedModelOptimizer::runBeamOptimization();
   //  FlatPattern pattern("/home/iason/Models/valid_6777.ply");
   //  FlatPattern pattern("/home/iason/Models/simple_beam_paper_example.ply");
   //  pattern.savePly("fannedValid.ply");
 
+  SimulationJob job;
+  job.load("../"
+           "ProblematicSimulationJobs/12:24_4_1_2021/"
+           "reduced_pattern_Single bar reduced model_simScenario.json");
+  FormFinder ff;
+  ff.executeSimulation(std::make_shared<SimulationJob>(job), false, false,
+                       true);
+
   // Create reduced models
   //  FormFinder::runUnitTests();
   const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
@@ -86,6 +95,9 @@ int main(int argc, char *argv[]) {
         //        std::filesystem::path(
         //            "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/"
         //            "1v_0v_2e_1e_1c_6fan/3/Valid/222.ply");
+        std::filesystem::path(
+            "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/"
+            "1v_0v_2e_1e_1c_6fan/3/Valid/624.ply");
     //        std::filesystem::path(
     //            "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/"
     //            "1v_0v_2e_1e_1c_6fan/2/Valid/33.ply");
@@ -94,7 +106,7 @@ int main(int argc, char *argv[]) {
     //            "1v_0v_2e_1e_1c_6fan/3/Valid/431.ply");
     //        std::filesystem::path(
     //            "/home/iason/Models/TestSet_validPatterns/59.ply");
-        entry.path();
+    //        entry.path();
     const auto filepathString = filepath.string();
     // Use only the base triangle version
     const std::string tiledSuffix = "_tiled.ply";
@@ -119,9 +131,10 @@ int main(int argc, char *argv[]) {
     //      optimizationExcludedEi.insert(0);
     //    }
     std::cout << "Full pattern:" << filepathString << std::endl;
-    optimizer.initializePatterns(pattern, *reducedModels[1],
+    FlatPattern cp;
-                                 optimizationExcludedEi);
+    cp.copy(*reducedModels[0]);
-    optimizer.optimize();
+    optimizer.initializePatterns(cp, *reducedModels[0], optimizationExcludedEi);
+    optimizer.optimize({ReducedModelOptimizer::SimulationScenario::Bending});
     //    }
   }
 

src/reducedmodeloptimizer.cpp

@@ -9,18 +9,25 @@ const bool gShouldDraw = true;
 
 FormFinder simulator;
 std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
-std::vector<SimulationJob> g_reducedPatternSimulationJob;
+std::vector<SimulationJob> g_fullPatternSimulationJob;
+std::vector<std::shared_ptr<SimulationJob>> g_reducedPatternSimulationJob;
 std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
-    g_reducedToFullInterfaceViMap;
+    reducedToFullInterfaceViMap;
+std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+    g_reducedToFullViMap;
 matplot::line_handle gPlotHandle;
 std::vector<double> gObjectiveValueHistory;
 Eigen::Vector4d g_initialX;
 std::unordered_set<size_t> g_reducedPatternExludedEdges;
 // double g_initialParameters;
 Eigen::VectorXd g_initialParameters;
-ReducedModelOptimizer::SimulationScenario g_chosenSimulationScenarioName;
+std::vector<ReducedModelOptimizer::SimulationScenario>
+    g_simulationScenarioIndices;
 std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)};
 double g_innerHexagonInitialPos = 0;
+bool g_optimizeInnerHexagonSize{false};
+std::vector<SimulationResults> firstOptimizationRoundResults;
+int g_firstRoundIterationIndex{0};
 
 // struct OptimizationCallback {
 //  double operator()(const size_t &iterations, const Eigen::VectorXd &x,
@@ -103,7 +110,7 @@ double ReducedModelOptimizer::computeError(
     const SimulationResults &reducedPatternResults,
     const Eigen::MatrixX3d &optimalReducedPatternDisplacements) {
   double error = 0;
-  for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) {
+  for (const auto reducedFullViPair : g_reducedToFullViMap) {
     VertexIndex reducedModelVi = reducedFullViPair.first;
     //    const auto pos =
     //        g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP();
@@ -124,15 +131,20 @@ double ReducedModelOptimizer::computeError(
 }
 
 void updateMesh(long n, const double *x) {
-  std::shared_ptr<SimulationMesh> pReducedPatternSimulationMesh =
+  std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
-      g_reducedPatternSimulationJob[0].mesh;
+      g_reducedPatternSimulationJob[0]->pMesh;
+  //  const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0];
+  //  std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << "
+  //  "
+  //            << elem.firstBendingConstFactor << " "
+  //            << elem.secondBendingConstFactor << std::endl;
   for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
     Element &e = pReducedPatternSimulationMesh->elements[ei];
     //    if (g_reducedPatternExludedEdges.contains(ei)) {
     //      continue;
     //    }
     //    e.properties.E = g_initialParameters * x[ei];
-    //    e.properties.E = g_initialParameters(ei, 0) * x[ei];
+    //    e.properties.E = g_initialParameters(0) * x[0];
     //    e.properties.G = g_initialParameters(1) * x[1];
     e.properties.A = g_initialParameters(0) * x[0];
     e.properties.J = g_initialParameters(1) * x[1];
@@ -147,36 +159,51 @@ void updateMesh(long n, const double *x) {
         2 * e.properties.E * e.properties.I3 / e.initialLength;
   }
 
-  if (n == 5) {
+  //  std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << "
+  //  "
+  //            << elem.firstBendingConstFactor << " "
+  //            << elem.secondBendingConstFactor << std::endl;
+  //  const Element &e = pReducedPatternSimulationMesh->elements[0];
+  //  std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
+  //            << e.firstBendingConstFactor << " " <<
+  //            e.secondBendingConstFactor
+  //            << std::endl;
+
+  if (g_optimizeInnerHexagonSize) {
     assert(pReducedPatternSimulationMesh->EN() == 12);
     for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN();
          vi += 2) {
       pReducedPatternSimulationMesh->vert[vi].P() =
-          g_innerHexagonVectors[vi / 2] * x[4];
+          g_innerHexagonVectors[vi / 2] * x[n - 1];
     }
 
-    pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
     pReducedPatternSimulationMesh->reset();
+    pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
     //    pReducedPatternSimulationMesh->registerForDrawing("Optimized
     //    hexagon"); polyscope::show();
   }
 }
 
-double ReducedModelOptimizer::objectiveAllScenarios(long n, const double *x) {
+double ReducedModelOptimizer::objective(long n, const double *x) {
   std::cout.precision(17);
 
   for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
     std::cout << "x[" + std::to_string(parameterIndex) + "]="
               << x[parameterIndex] << std::endl;
   }
+  const Element &e = g_reducedPatternSimulationJob[0]->pMesh->elements[0];
+  //  std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
+  //            << e.firstBendingConstFactor << " " <<
+  //            e.secondBendingConstFactor
+  //            << std::endl;
   updateMesh(n, x);
+  std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
+            << e.firstBendingConstFactor << " " << e.secondBendingConstFactor
+            << std::endl;
 
   // run simulations
   double error = 0;
-  for (int simulationScenarioIndex = SimulationScenario::Axial;
+  for (const int simulationScenarioIndex : g_simulationScenarioIndices) {
-       simulationScenarioIndex !=
-       SimulationScenario::NumberOfSimulationScenarios;
-       simulationScenarioIndex++) {
     SimulationResults reducedModelResults = simulator.executeSimulation(
         g_reducedPatternSimulationJob[simulationScenarioIndex], false, false);
     error += computeError(
@@ -188,32 +215,13 @@ double ReducedModelOptimizer::objectiveAllScenarios(long n, const double *x) {
   gObjectiveValueHistory.push_back(error);
   auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
                                  gObjectiveValueHistory.size());
-  gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory);
+  std::vector<double> colors(gObjectiveValueHistory.size(), 2);
-
+  if (g_firstRoundIterationIndex != 0) {
-  return error;
+    for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(),
-}
+             [](double &c) { c = 0.7; });
-
-double ReducedModelOptimizer::objectiveSingleScenario(long n, const double *x) {
-  std::cout.precision(17);
-
-  for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
-    std::cout << "x[" + std::to_string(parameterIndex) + "]="
-              << x[parameterIndex] << std::endl;
   }
-
+  gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors);
-  updateMesh(n, x);
+  std::cout << std::endl;
-
-  SimulationResults reducedModelResults = simulator.executeSimulation(
-      g_reducedPatternSimulationJob[g_chosenSimulationScenarioName], false,
-      false);
-  const double error = computeError(
-      reducedModelResults,
-      g_optimalReducedModelDisplacements[g_chosenSimulationScenarioName]);
-
-  gObjectiveValueHistory.push_back(error);
-  auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
-                                 gObjectiveValueHistory.size());
-  gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory);
 
   return error;
 }
@@ -257,7 +265,7 @@ void ReducedModelOptimizer::computeMaps(
   }
 
   // Construct reduced->full and full->reduced interface vi map
-  g_reducedToFullInterfaceViMap.clear();
+  reducedToFullInterfaceViMap.clear();
   vcg::tri::Allocator<FlatPattern>::PointerUpdater<FlatPattern::VertexPointer>
       pu_reducedModel;
   reducedPattern.deleteDanglingVertices(pu_reducedModel);
@@ -267,13 +275,13 @@ void ReducedModelOptimizer::computeMaps(
       reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/;
   reducedPattern.createFan();
   for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
-    g_reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex +
+    reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex +
                                 reducedModelBaseTriangleInterfaceVi] =
         fullModelBaseTriangleInterfaceVi +
         fanIndex * fullPatternInterfaceVertexOffset;
   }
   m_fullToReducedInterfaceViMap.clear();
-  constructInverseMap(g_reducedToFullInterfaceViMap,
+  constructInverseMap(reducedToFullInterfaceViMap,
                       m_fullToReducedInterfaceViMap);
 
   //  fullPattern.setLabel("FullPattern");
@@ -288,6 +296,8 @@ void ReducedModelOptimizer::computeMaps(
     m_fullPatternOppositeInterfaceViMap[vi0] = vi1;
   }
 
+  g_reducedToFullViMap = reducedToFullInterfaceViMap;
+
   const bool debugMapping = false;
   if (debugMapping) {
     reducedPattern.registerForDrawing();
@@ -322,11 +332,11 @@ void ReducedModelOptimizer::computeMaps(
     std::vector<glm::vec3> nodeColorsReducedToFull_full(fullPattern.VN(),
                                                         glm::vec3(0, 0, 0));
     for (size_t vi = 0; vi < reducedPattern.VN(); vi++) {
-      if (g_reducedToFullInterfaceViMap.contains(vi)) {
+      if (reducedToFullInterfaceViMap.contains(vi)) {
 
         auto color = polyscope::getNextUniqueColor();
         nodeColorsReducedToFull_reduced[vi] = color;
-        nodeColorsReducedToFull_full[g_reducedToFullInterfaceViMap[vi]] = color;
+        nodeColorsReducedToFull_full[reducedToFullInterfaceViMap[vi]] = color;
       }
     }
     polyscope::getCurveNetwork(reducedPattern.getLabel())
@@ -352,10 +362,10 @@ void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel,
   m_pReducedPatternSimulationMesh->setBeamCrossSection(
       CrossSectionType{0.002, 0.002});
   m_pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1);
-  m_pFullModelSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
+  m_pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
-  m_pFullModelSimulationMesh->setBeamCrossSection(
+  m_pFullPatternSimulationMesh->setBeamCrossSection(
       CrossSectionType{0.002, 0.002});
-  m_pFullModelSimulationMesh->setBeamMaterial(0.3, 1);
+  m_pFullPatternSimulationMesh->setBeamMaterial(0.3, 1);
 }
 
 ReducedModelOptimizer::ReducedModelOptimizer(
@@ -367,6 +377,8 @@ ReducedModelOptimizer::ReducedModelOptimizer(
 void ReducedModelOptimizer::initializePatterns(
     FlatPattern &fullPattern, FlatPattern &reducedPattern,
     const std::unordered_set<size_t> &reducedModelExcludedEdges) {
+  fullPattern.setLabel("full_pattern_" + fullPattern.getLabel());
+  reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel());
   assert(fullPattern.VN() == reducedPattern.VN() &&
          fullPattern.EN() >= reducedPattern.EN());
   polyscope::removeAllStructures();
@@ -375,7 +387,8 @@ void ReducedModelOptimizer::initializePatterns(
   FlatPattern copyReducedPattern;
   copyFullPattern.copy(fullPattern);
   copyReducedPattern.copy(reducedPattern);
-  if (copyReducedPattern.EN() == 2) {
+  g_optimizeInnerHexagonSize = copyReducedPattern.EN() == 2;
+  if (g_optimizeInnerHexagonSize) {
     const double h = copyReducedPattern.getBaseTriangleHeight();
     double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0));
     VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h,
@@ -391,19 +404,20 @@ void ReducedModelOptimizer::initializePatterns(
           baseTriangle_leftBottomNode;
       g_innerHexagonVectors[rotationCounter] = rotatedVector;
     }
-  }
-  computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges);
     const double innerHexagonInitialPos_x =
         copyReducedPattern.vert[0].cP()[0] / g_innerHexagonVectors[0][0];
     const double innerHexagonInitialPos_y =
         copyReducedPattern.vert[0].cP()[1] / g_innerHexagonVectors[0][1];
     g_innerHexagonInitialPos = innerHexagonInitialPos_x;
+  }
+  computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges);
   createSimulationMeshes(copyFullPattern, copyReducedPattern);
-  initializeStiffnesses();
+  initializeStiffnesses(m_pReducedPatternSimulationMesh);
 }
 
-void ReducedModelOptimizer::initializeStiffnesses() {
+void ReducedModelOptimizer::initializeStiffnesses(
-  g_initialParameters.resize(4);
+    const std::shared_ptr<SimulationMesh> &mesh) {
+  g_initialParameters.resize(g_optimizeInnerHexagonSize ? 5 : 4);
   // Save save the beam stiffnesses
   //  for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) {
   //    Element &e = pReducedModelElementalMesh->elements[ei];
@@ -417,73 +431,62 @@ void ReducedModelOptimizer::initializeStiffnesses() {
   //  g_initialParameters =
   //      m_pReducedPatternSimulationMesh->elements[0].properties.E;
   //  for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) {
-  //    g_initialParameters(ei, 0) =
-  //        m_pReducedPatternSimulationMesh->elements[ei].properties.E;
   //  }
-  //    g_initialParameters(1) =
+  //  g_initialParameters(0) = mesh->elements[0].properties.E;
-  //    pReducedModelElementalMesh->elements[0].properties.G;
+  //  g_initialParameters(1) = mesh->elements[0].properties.G;
-  g_initialParameters(0) =
+  g_initialParameters(0) = mesh->elements[0].properties.A;
-      m_pReducedPatternSimulationMesh->elements[0].properties.A;
+  g_initialParameters(1) = mesh->elements[0].properties.J;
-  g_initialParameters(1) =
+  g_initialParameters(2) = mesh->elements[0].properties.I2;
-      m_pReducedPatternSimulationMesh->elements[0].properties.J;
+  g_initialParameters(3) = mesh->elements[0].properties.I3;
-  g_initialParameters(2) =
-      m_pReducedPatternSimulationMesh->elements[0].properties.I2;
-  g_initialParameters(3) =
-      m_pReducedPatternSimulationMesh->elements[0].properties.I3;
   //  }
 }
 
 void ReducedModelOptimizer::computeReducedModelSimulationJob(
     const SimulationJob &simulationJobOfFullModel,
+    const std::unordered_map<size_t, size_t> &simulationJobFullToReducedMap,
     SimulationJob &simulationJobOfReducedModel) {
+  assert(simulationJobOfReducedModel.pMesh->VN() != 0);
   std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
       reducedModelFixedVertices;
-  for (auto fullModelFixedVertex : simulationJobOfFullModel.fixedVertices) {
+  for (auto fullModelFixedVertex :
-    reducedModelFixedVertices[m_fullToReducedInterfaceViMap.at(
+       simulationJobOfFullModel.constrainedVertices) {
+    reducedModelFixedVertices[simulationJobFullToReducedMap.at(
         fullModelFixedVertex.first)] = fullModelFixedVertex.second;
   }
 
   std::unordered_map<VertexIndex, Vector6d> reducedModelNodalForces;
   for (auto fullModelNodalForce :
        simulationJobOfFullModel.nodalExternalForces) {
-    reducedModelNodalForces[m_fullToReducedInterfaceViMap.at(
+    reducedModelNodalForces[simulationJobFullToReducedMap.at(
         fullModelNodalForce.first)] = fullModelNodalForce.second;
   }
 
-  std::unordered_map<VertexIndex, VectorType> reducedModelNodalForcedNormals;
+  //  std::unordered_map<VertexIndex, VectorType>
-  for (auto fullModelNodalForcedRotation :
+  //  reducedModelNodalForcedNormals; for (auto fullModelNodalForcedRotation :
-       simulationJobOfFullModel.nodalForcedNormals) {
+  //       simulationJobOfFullModel.nodalForcedNormals) {
-    reducedModelNodalForcedNormals[m_fullToReducedInterfaceViMap.at(
+  //    reducedModelNodalForcedNormals[simulationJobFullToReducedMap.at(
-        fullModelNodalForcedRotation.first)] =
+  //        fullModelNodalForcedRotation.first)] =
-        fullModelNodalForcedRotation.second;
+  //        fullModelNodalForcedRotation.second;
-  }
+  //  }
-
+  simulationJobOfReducedModel.constrainedVertices = reducedModelFixedVertices;
-  simulationJobOfReducedModel = SimulationJob{m_pReducedPatternSimulationMesh,
+  simulationJobOfReducedModel.nodalExternalForces = reducedModelNodalForces;
-                                              reducedModelFixedVertices,
+  simulationJobOfReducedModel.label = simulationJobOfFullModel.getLabel();
-                                              reducedModelNodalForces,
+  //  simulationJobOfReducedModel.nodalForcedNormals =
-                                              {},
+  //      reducedModelNodalForcedNormals;
-                                              reducedModelNodalForcedNormals};
-}
-
-SimulationJob ReducedModelOptimizer::getReducedSimulationJob(
-    const SimulationJob &fullModelSimulationJob) {
-  SimulationJob reducedModelSimulationJob;
-  computeReducedModelSimulationJob(fullModelSimulationJob,
-                                   reducedModelSimulationJob);
-  return reducedModelSimulationJob;
 }
 
 void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
     const SimulationResults &fullModelResults,
+    const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
     Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) {
 
-  optimalDisplacementsOfReducedModel.resize(
+  assert(optimalDisplacementsOfReducedModel.rows() != 0 &&
-      m_pReducedPatternSimulationMesh->VN(), 3);
+         optimalDisplacementsOfReducedModel.cols() == 3);
   optimalDisplacementsOfReducedModel.setZero(
       optimalDisplacementsOfReducedModel.rows(),
       optimalDisplacementsOfReducedModel.cols());
 
-  for (auto reducedFullViPair : g_reducedToFullInterfaceViMap) {
+  for (auto reducedFullViPair : displacementsReducedToFullMap) {
     const VertexIndex fullModelVi = reducedFullViPair.second;
     const Vector6d fullModelViDisplacements =
         fullModelResults.displacements[fullModelVi];
@@ -499,7 +502,7 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization(
 
   gObjectiveValueHistory.clear();
 
-  const size_t n = m_pReducedPatternSimulationMesh->EN() == 12 ? 5 : 4;
+  const size_t n = g_optimizeInnerHexagonSize ? 5 : 4;
   const size_t npt = 2 * n;
   //      ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
   assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2);
@@ -507,39 +510,58 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization(
                              "conditions is not recommended.");
   // Set initial guess of solution
 
-  std::vector<double> x(n, 100);
+  std::vector<double> x(n, 50);
-  x[4] = g_innerHexagonInitialPos;
+  if (g_optimizeInnerHexagonSize) {
-  if (!initialGuess.empty()) {
+    x[n - 1] = g_innerHexagonInitialPos;
-    x = initialGuess;
+  }
-  } //      {0.10000000000000	001, 2, 1.9999999971613847, 6.9560343643347764};
+  /*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)};
   const double xMin = 1e-2;
-  const double xMax = 5000;
+  const double xMax = 500;
   //  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 (n == 5) {
+  if (g_optimizeInnerHexagonSize) {
-    xLow[4] = 0.1;
+    xLow[n - 1] = 0.1;
-    xUpper[4] = 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);
-  const double rhobeg = 0.01;
+  const double rhobeg = 1;
   //  const double rhobeg = 10;
   const double rhoend = rhobeg * 1e-6;
   const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2;
   std::vector<double> w(wSize);
   //  const size_t maxFun = 10 * (x.size() ^ 2);
-  const size_t maxFun = 120;
+  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();
+  bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
+         rhobeg * 1e-3, rhoend, maxFun, w.data());
+  std::cout << "Finished second optimization round" << std::endl;
 
   Eigen::VectorXd eigenX(x.size(), 1);
   for (size_t xi = 0; xi < x.size(); xi++) {
@@ -552,12 +574,13 @@ void ReducedModelOptimizer::setInitialGuess(std::vector<double> v) {
   initialGuess = v;
 }
 
-std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
+std::vector<std::shared_ptr<SimulationJob>>
+ReducedModelOptimizer::createScenarios(
     const std::shared_ptr<SimulationMesh> &pMesh) {
-  std::vector<SimulationJob> scenarios;
+  std::vector<std::shared_ptr<SimulationJob>> scenarios;
+  scenarios.resize(SimulationScenario::NumberOfSimulationScenarios);
   std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
   std::unordered_map<VertexIndex, Vector6d> nodalForces;
-  // NOTE: Assuming that the first interface node lays on the y axis
   const double forceMagnitude = 1;
   // Assuming the patterns lays on the x-y plane
   const CoordType patternPlaneNormal(0, 0, 1);
@@ -594,6 +617,7 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
   //  m_pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
 
   //// Axial
+  SimulationScenario scenarioName = SimulationScenario::Axial;
   for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
     CoordType forceDirection =
         (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
@@ -604,9 +628,12 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
     fixedVertices[viPair.second] =
         std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
   }
-  scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
+  scenarios[scenarioName] = std::make_shared<SimulationJob>(
+      SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
+                    fixedVertices, nodalForces, {}));
 
-  //// In-plane Bending
+  //// Shear
+  scenarioName = SimulationScenario::Shear;
   fixedVertices.clear();
   nodalForces.clear();
   for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
@@ -620,7 +647,9 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
     fixedVertices[viPair.second] =
         std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
   }
-  scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
+  scenarios[scenarioName] = std::make_shared<SimulationJob>(
+      SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
+                    fixedVertices, nodalForces, {}));
 
   //  //// Torsion
   //  fixedVertices.clear();
@@ -645,7 +674,8 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
   //  }
   //  scenarios.push_back({pMesh, fixedVertices, nodalForces});
 
-  ////  Out - of - plane bending.Pull towards Z
+  ////  Bending
+  scenarioName = SimulationScenario::Bending;
   fixedVertices.clear();
   nodalForces.clear();
   for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
@@ -653,9 +683,12 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
     fixedVertices[viPair.second] =
         std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
   }
-  scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
+  scenarios[scenarioName] = std::make_shared<SimulationJob>(
+      SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
+                    fixedVertices, nodalForces, {}));
 
   //// Double using moments
+  scenarioName = SimulationScenario::Dome;
   fixedVertices.clear();
   nodalForces.clear();
   for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
@@ -676,9 +709,12 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
     nodalForces[viPair.second] =
         Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 0.1;
   }
-  scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
+  scenarios[scenarioName] = std::make_shared<SimulationJob>(
+      SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
+                    fixedVertices, nodalForces, {}));
 
   //// Saddle
+  scenarioName = SimulationScenario::Saddle;
   fixedVertices.clear();
   nodalForces.clear();
   for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
@@ -702,76 +738,177 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
           Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.01 * forceMagnitude;
     }
   }
-  scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
+  scenarios[scenarioName] = std::make_shared<SimulationJob>(
+      SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
+                    fixedVertices, nodalForces, {}));
 
   return scenarios;
 }
 
-void ReducedModelOptimizer::optimize(const int &simulationScenario) {
+void ReducedModelOptimizer::runBeamOptimization() {
+  // load beams
+  VCGEdgeMesh fullBeam;
+  fullBeam.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply");
+  VCGEdgeMesh reducedBeam;
+  reducedBeam.loadFromPly("/home/iason/Models/simple_beam_model_4elem_1m.ply");
+  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;
 
-  std::vector<SimulationJob> simulationJobs =
+  // full model simuilation job
-      createScenarios(m_pFullModelSimulationMesh);
+  auto pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullBeam);
+  pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02});
+  pFullPatternSimulationMesh->setBeamMaterial(0.3, 1);
+  std::unordered_map<VertexIndex, std::unordered_set<int>> fixedVertices;
+  fixedVertices[0] = ::unordered_set<int>({0, 1, 2, 3, 4, 5});
+  std::unordered_map<VertexIndex, Vector6d> forces;
+  forces[fullBeam.VN() - 1] = Vector6d({0, 0, 10, 0, 0, 0});
+  const std::string fullBeamSimulationJobLabel = "Pull_Z";
+  std::shared_ptr<SimulationJob> pFullModelSimulationJob =
+      make_shared<SimulationJob>(SimulationJob(pFullPatternSimulationMesh,
+                                               fullBeamSimulationJobLabel,
+                                               fixedVertices, forces));
+  FormFinder formFinder;
+  auto fullModelResults = formFinder.executeSimulation(pFullModelSimulationJob);
+
+  // Optimal reduced model displacements
+  const size_t numberOfSimulationScenarios = 1;
+  g_optimalReducedModelDisplacements.resize(numberOfSimulationScenarios);
+  g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1].resize(
+      reducedBeam.VN(), 3);
+  computeDesiredReducedModelDisplacements(
+      fullModelResults, displacementReducedToFullMap,
+      g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]);
+
+  // reduced model simuilation job
+  auto reducedSimulationMesh = std::make_shared<SimulationMesh>(reducedBeam);
+  reducedSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02});
+  reducedSimulationMesh->setBeamMaterial(0.3, 1);
+  g_reducedPatternSimulationJob.resize(numberOfSimulationScenarios);
+  SimulationJob reducedSimJob;
+  computeReducedModelSimulationJob(*pFullModelSimulationJob,
+                                   jobFullToReducedMap, reducedSimJob);
+  reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] =
+      reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] * 0.1;
+  g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1] =
+      make_shared<SimulationJob>(
+          reducedSimulationMesh, fullBeamSimulationJobLabel,
+          reducedSimJob.constrainedVertices, reducedSimJob.nodalExternalForces);
+  initializeStiffnesses(reducedSimulationMesh);
+
+  //  const std::string simulationJobsPath = "SimulationJobs";
+  //  std::filesystem::create_directory(simulationJobsPath);
+  //  g_reducedPatternSimulationJob[0].save(simulationJobsPath);
+  //  g_reducedPatternSimulationJob[0].load(
+  //      std::filesystem::path(simulationJobsPath)
+  //          .append(g_reducedPatternSimulationJob[0].mesh->getLabel() +
+  //                  "_simScenario.json"));
+
+  runOptimization(&objective);
+
+  fullModelResults.registerForDrawing();
+  SimulationResults reducedModelResults = simulator.executeSimulation(
+      g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1], false,
+      false);
+  double error = computeError(
+      reducedModelResults,
+      g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]);
+  reducedModelResults.registerForDrawing();
+  std::cout << "Error between beams:" << error << endl;
+  //    registerWorldAxes();
+  polyscope::show();
+  fullModelResults.unregister();
+  reducedModelResults.unregister();
+}
+
+void ReducedModelOptimizer::visualizeResults(
+    const std::vector<std::shared_ptr<SimulationJob>>
+        &fullPatternSimulationJobs,
+    const std::vector<SimulationScenario> &simulationScenarios) {
+  m_pFullPatternSimulationMesh->registerForDrawing();
+  double error = 0;
+  for (const int simulationScenarioIndex : simulationScenarios) {
+    const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
+        fullPatternSimulationJobs[simulationScenarioIndex];
+    pFullPatternSimulationJob->registerForDrawing(
+        m_pFullPatternSimulationMesh->getLabel());
+    SimulationResults fullModelResults =
+        simulator.executeSimulation(pFullPatternSimulationJob, false);
+    fullModelResults.registerForDrawing();
+    fullModelResults.saveDeformedModel();
+    const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
+        g_reducedPatternSimulationJob[simulationScenarioIndex];
+    SimulationResults reducedModelResults =
+        simulator.executeSimulation(pReducedPatternSimulationJob, false, false);
+    error += computeError(
+        reducedModelResults,
+        g_optimalReducedModelDisplacements[simulationScenarioIndex]);
+    std::cout << "Error of simulation scenario "
+              << simulationScenarioStrings[simulationScenarioIndex] << " is "
+              << error << std::endl;
+    reducedModelResults.registerForDrawing();
+    firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
+    reducedModelResults.saveDeformedModel();
+    //    registerWorldAxes();
+    polyscope::show();
+    fullModelResults.unregister();
+    reducedModelResults.unregister();
+    firstOptimizationRoundResults[simulationScenarioIndex].unregister();
+  }
+}
+
+void ReducedModelOptimizer::optimize(
+    const std::vector<SimulationScenario> &simulationScenariosIndices) {
+
+  std::vector<std::shared_ptr<SimulationJob>> simulationJobs =
+      createScenarios(m_pFullPatternSimulationMesh);
   g_optimalReducedModelDisplacements.resize(6);
   g_reducedPatternSimulationJob.resize(6);
+  g_firstRoundIterationIndex = 0;
   //  polyscope::removeAllStructures();
 
   for (int simulationScenarioIndex = SimulationScenario::Axial;
        simulationScenarioIndex !=
        SimulationScenario::NumberOfSimulationScenarios;
        simulationScenarioIndex++) {
-    const SimulationJob &fullPatternSimulationJob =
+    const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
         simulationJobs[simulationScenarioIndex];
     //  fullPatternSimulationJob.mesh->savePly(
     //      "Fanned_" + m_pFullModelSimulationMesh->getLabel() + ".ply");
     SimulationResults fullModelResults =
-        simulator.executeSimulation(fullPatternSimulationJob, false);
+        simulator.executeSimulation(pFullPatternSimulationJob);
     //    fullModelResults.label =
     //        "fullModel_" + SimulationScenarioStrings[simulationScenarioIndex];
     //    fullModelResults.registerForDrawing(fullPatternSimulationJob);
+    g_optimalReducedModelDisplacements[simulationScenarioIndex].resize(
+        m_pReducedPatternSimulationMesh->VN(), 3);
     computeDesiredReducedModelDisplacements(
-        fullModelResults,
+        fullModelResults, g_reducedToFullViMap,
         g_optimalReducedModelDisplacements[simulationScenarioIndex]);
-    computeReducedModelSimulationJob(
+    SimulationJob reducedPatternSimulationJob;
-        fullPatternSimulationJob,
+    reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
-        g_reducedPatternSimulationJob[simulationScenarioIndex]);
+    computeReducedModelSimulationJob(*pFullPatternSimulationJob,
+                                     m_fullToReducedInterfaceViMap,
+                                     reducedPatternSimulationJob);
+    g_reducedPatternSimulationJob[simulationScenarioIndex] =
+        std::make_shared<SimulationJob>(reducedPatternSimulationJob);
   }
 
-  if (simulationScenario == -1) {
+  g_simulationScenarioIndices = simulationScenariosIndices;
-    double (*pObjectiveFunction)(long, const double *) = &objectiveAllScenarios;
+  if (simulationScenariosIndices.empty()) {
-    runOptimization(pObjectiveFunction);
+    g_simulationScenarioIndices = {
-  } else { // run chosen
+        SimulationScenario::Axial, SimulationScenario::Shear,
-    g_chosenSimulationScenarioName = SimulationScenario(simulationScenario);
+        SimulationScenario::Bending, SimulationScenario::Dome,
-    double (*pObjectiveFunction)(long, const double *) =
+        SimulationScenario::Saddle};
-        &objectiveSingleScenario;
-    runOptimization(pObjectiveFunction);
   }
 
-  m_pFullModelSimulationMesh->registerForDrawing();
+  runOptimization(&objective);
-  double error = 0;
+  visualizeResults(simulationJobs, g_simulationScenarioIndices);
-  for (int simulationScenarioIndex = SimulationScenario::Axial;
-       simulationScenarioIndex !=
-       SimulationScenario::NumberOfSimulationScenarios;
-       simulationScenarioIndex++) {
-    const SimulationJob &fullPatternSimulationJob =
-        simulationJobs[simulationScenarioIndex];
-    SimulationResults fullModelResults =
-        simulator.executeSimulation(fullPatternSimulationJob, false);
-    fullModelResults.label =
-        "fullModel_" + SimulationScenarioStrings[simulationScenarioIndex];
-    fullModelResults.registerForDrawing(fullPatternSimulationJob);
-    const SimulationJob &reducedPatternSimulationJob =
-        g_reducedPatternSimulationJob[simulationScenarioIndex];
-    SimulationResults reducedModelResults =
-        simulator.executeSimulation(reducedPatternSimulationJob, false, false);
-    error += computeError(
-        reducedModelResults,
-        g_optimalReducedModelDisplacements[simulationScenarioIndex]);
-    reducedModelResults.label =
-        "reducedModel_" + SimulationScenarioStrings[simulationScenarioIndex];
-    reducedModelResults.registerForDrawing(reducedPatternSimulationJob);
-    //    registerWorldAxes();
-    polyscope::show();
-    fullModelResults.unregister();
-    reducedModelResults.unregister();
-  }
 }

src/reducedmodeloptimizer.hpp

@@ -12,7 +12,7 @@ using ReducedPatternVertexIndex = VertexIndex;
 
 class ReducedModelOptimizer {
   std::shared_ptr<SimulationMesh> m_pReducedPatternSimulationMesh;
-  std::shared_ptr<SimulationMesh> m_pFullModelSimulationMesh;
+  std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh;
   std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
       m_fullToReducedInterfaceViMap;
   std::unordered_map<FullPatternVertexIndex, FullPatternVertexIndex>
@@ -26,18 +26,19 @@ public:
     Axial,
     Shear,
     Bending,
-    Double,
+    Dome,
     Saddle,
     NumberOfSimulationScenarios
   };
-  inline static const std::string SimulationScenarioStrings[] = {
+  inline static const std::string simulationScenarioStrings[] = {
       "Axial", "Shear", "Bending", "Double", "Saddle"};
-  void optimize(const int &simulationScenario = -1);
+  void optimize(const std::vector<SimulationScenario> &simulationScenarios);
   double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;
 
   ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot);
-  void computeReducedModelSimulationJob(
+  static void computeReducedModelSimulationJob(
       const SimulationJob &simulationJobOfFullModel,
+      const std::unordered_map<size_t, size_t> &simulationJobFullToReducedMap,
       SimulationJob &simulationJobOfReducedModel);
 
   SimulationJob
@@ -49,25 +50,35 @@ public:
 
   void setInitialGuess(std::vector<double> v);
 
+  static void runBeamOptimization();
+
+  static void runSimulation(const std::string &filename,
+                            std::vector<double> &x);
+
 private:
-  void computeDesiredReducedModelDisplacements(
+  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);
-  Eigen::VectorXd runOptimization(double (*pObjectiveFunction)(long,
+  static Eigen::VectorXd
-                                                               const double *));
+  runOptimization(double (*pObjectiveFunction)(long, const double *));
-  std::vector<SimulationJob>
+  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);
-  void initializeStiffnesses();
+  static void
+  initializeStiffnesses(const std::shared_ptr<SimulationMesh> &mesh);
 
   static double
   computeError(const SimulationResults &reducedPatternResults,
                const Eigen::MatrixX3d &optimalReducedPatternDisplacements);
-  static double objectiveSingleScenario(long n, const double *x);
+  static double objective(long n, const double *x);
-  static double objectiveAllScenarios(long n, const double *x);
 };
 
 #endif // REDUCEDMODELOPTIMIZER_HPP