Computation of the max force for each base simulation scenario based on a desired displacement.

TestSet/ReducedPatterns/SingleBar_reduced.ply

@@ -1,25 +0,0 @@
-ply
-format ascii 1.0
-comment VCGLIB generated
-element vertex 7
-property double x
-property double y
-property double z
-property uchar red
-property uchar green
-property uchar blue
-property uchar alpha
-element face 0
-property list uchar int vertex_indices
-element edge 1
-property int vertex1
-property int vertex2
-end_header
-0 0 0 255 255 255 255 
--0.1666666666666666 -0.2886751345948129 0 255 255 255 255 
--0.3333333333333333 -0.5773502691896257 0 255 255 255 255 
-0 -0.8660254037844387 0 255 255 255 255 
-0.3333333333333333 -0.5773502691896258 0 255 255 255 255 
-0.1666666666666666 -0.288675134594813 0 255 255 255 255 
-0.1 -0.5773502691896258 0 255 255 255 255 
-0 3 

src/main.cpp

@@ -1,5 +1,5 @@
-#include "drmsimulationmodel.hpp"
 #include "csvfile.hpp"
+#include "drmsimulationmodel.hpp"
 #include "edgemesh.hpp"
 #include "reducedmodeloptimizer.hpp"
 #include "simulationhistoryplotter.hpp"
@@ -57,6 +57,8 @@ int main(int argc, char *argv[]) {
       = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
   settings_optimization.normalizationParameter = 0.0003;
   settings_optimization.solutionAccuracy = 0.001;
+  settings_optimization.objectiveWeights.translational = std::atof(argv[4]);
+  settings_optimization.objectiveWeights.rotational = 2 - std::atof(argv[4]);
 
   // Optimize pair
   const std::string pairName =
@@ -66,26 +68,35 @@ int main(int argc, char *argv[]) {
   assert(interfaceNodeIndex==numberOfNodesPerSlot[0]+numberOfNodesPerSlot[3]);
   ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
   optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges.size());
-  ReducedPatternOptimization::Results optimizationResults = optimizer.optimize(
+  ReducedPatternOptimization::Results optimizationResults
-      settings_optimization);
+      //   = optimizer.optimize(
+      //      settings_optimization);
+      = optimizer.optimize(settings_optimization);
 
   // Export results
-  const bool input_resultDirectoryDefined = argc >= 5;
+  const bool input_resultDirectoryDefined = argc >= 6;
-  std::string optimizationResultsDirectory =
+  std::string optimizationResultsDirectory = input_resultDirectoryDefined
-      input_resultDirectoryDefined ? argv[4] : std::filesystem::current_path().append("OptimizationResults");
+                                                 ? argv[5]
+                                                 : std::filesystem::current_path().append(
+                                                     "OptimizationResults");
   std::string resultsOutputDir;
   if (optimizationResults.numberOfSimulationCrashes != 0) {
-    const auto crashedJobsDirPath =
+      const auto crashedJobsDirPath = std::filesystem::path(
-        std::filesystem::path(optimizationResultsDirectory)
+          optimizationResultsDirectory.append("CrashedJobs")
-            .append("CrashedJobs")
+              .append(
-            .append(pairName);
+                  pairName + "(" + std::to_string(settings_optimization.numberOfFunctionCalls) + "_"
+                  + to_string_with_precision(settings_optimization.objectiveWeights.translational)
+                  + ")"));
       std::filesystem::create_directories(crashedJobsDirPath);
       resultsOutputDir = crashedJobsDirPath.string();
   } else {
       std::filesystem::path convergedJobsDirPath(
           std::filesystem::path(optimizationResultsDirectory)
               .append("ConvergedJobs")
-            .append(pairName));
+              .append(
+                  pairName + "(" + std::to_string(settings_optimization.numberOfFunctionCalls) + "_"
+                  + to_string_with_precision(settings_optimization.objectiveWeights.translational)
+                  + ")"));
       std::filesystem::create_directories(convergedJobsDirPath);
       resultsOutputDir = convergedJobsDirPath.string();
   }

src/reducedmodeloptimizer.cpp

@@ -16,6 +16,8 @@ struct GlobalOptimizationVariables {
     std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
     std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
         reducedToFullInterfaceViMap;
+    std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+        fullPatternInterfaceViPairs;
     matplot::line_handle gPlotHandle;
     std::vector<double> objectiveValueHistory;
     Eigen::VectorXd initialParameters;
@@ -27,6 +29,14 @@ struct GlobalOptimizationVariables {
     int numberOfOptimizationParameters{5};
     ReducedPatternOptimization::Settings optimizationSettings;
     vcg::Triangle3<double> baseTriangle;
+    //Variables for finding the full pattern simulation forces
+    std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh;
+    std::function<void(const double &,
+                       const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>> &,
+                       SimulationJob &)>
+        constructScenarioFunction;
+    FullPatternVertexIndex interfaceViForComputingScenarioError;
+    double desiredMaxDisplacementValue;
 } global;
 
 double ReducedModelOptimizer::computeDisplacementError(
@@ -112,7 +122,8 @@ double ReducedModelOptimizer::computeError(
                                  simulationResults_reducedPattern.rotationalDisplacementQuaternion,
                                  reducedToFullInterfaceViMap,
                                  normalizationFactor_rotationalDisplacement);
-    return 1.5 * translationalError + 0.5 * rotationalError;
+    return global.optimizationSettings.objectiveWeights.translational * translationalError
+           + global.optimizationSettings.objectiveWeights.rotational * rotationalError;
 }
 
 double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3,
@@ -129,6 +140,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
     //        }
     //    std::cout << std::endl;
 
+    //    std::cout << x[n - 2] << " " << x[n - 1] << std::endl;
     //  const Element &e =
     //  global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout <<
     //  e.axialConstFactor << " " << e.torsionConstFactor << " "
@@ -136,6 +148,11 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
     //            e.secondBendingConstFactor
     //            << std::endl;
     updateMesh(n, x);
+    //    global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing();
+    //    global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing();
+    //    polyscope::show();
+    //    global.reducedPatternSimulationJobs[0]->pMesh->unregister();
+
     //  std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
     //            << e.firstBendingConstFactor << " " <<
     //            e.secondBendingConstFactor
@@ -144,10 +161,11 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
     // run simulations
     double totalError = 0;
     LinearSimulationModel simulator;
-//    FormFinder simulator;
+    //    FormFinder simulator;
     for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
-        SimulationResults reducedModelResults = simulator.executeSimulation(
+        const std::shared_ptr<SimulationJob> &reducedJob
-            global.reducedPatternSimulationJobs[simulationScenarioIndex]);
+            = global.reducedPatternSimulationJobs[simulationScenarioIndex];
+        SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob);
         //    std::string filename;
         if (!reducedModelResults.converged) {
             totalError += std::numeric_limits<double>::max();
@@ -218,7 +236,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
     }
 #ifdef POLYSCOPE_DEFINED
     ++global.numberOfFunctionCalls;
-    if (global.numberOfFunctionCalls % 100 == 0) {
+    if (global.numberOfFunctionCalls % 1000 == 0) {
         std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl;
     }
 #endif
@@ -276,7 +294,7 @@ void ReducedModelOptimizer::computeMaps(
     std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
         &fullToReducedInterfaceViMap,
     std::vector<std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex>>
-        &fullPatternOppositeInterfaceViMap)
+        &fullPatternOppositeInterfaceViPairs)
 {
     // Compute the offset between the interface nodes
     const size_t interfaceSlotIndex = 4; // bottom edge
@@ -320,7 +338,11 @@ void ReducedModelOptimizer::computeMaps(
         pu_reducedModel.remap[baseTriangleInterfaceVi];
     const size_t reducedModelInterfaceVertexOffset
         = reducedPattern.VN() /*- 1*/ /*- reducedModelBaseTriangleInterfaceVi*/;
-    reducedPattern.createFan({1}); //TODO: should be an input parameter from main
+    Results::applyOptimizationResults_innerHexagon(initialHexagonSize,
+                                                   0,
+                                                   global.baseTriangle,
+                                                   reducedPattern);
+    reducedPattern.createFan({0}); //TODO: should be an input parameter from main
     for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
         reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex +
                                     reducedModelBaseTriangleInterfaceVi] =
@@ -331,20 +353,21 @@ void ReducedModelOptimizer::computeMaps(
     constructInverseMap(reducedToFullInterfaceViMap, fullToReducedInterfaceViMap);
 
     // Create opposite vertex map
-    fullPatternOppositeInterfaceViMap.clear();
+    fullPatternOppositeInterfaceViPairs.clear();
-    fullPatternOppositeInterfaceViMap.reserve(fanSize / 2);
+    fullPatternOppositeInterfaceViPairs.reserve(fanSize / 2);
     //    for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) {
     for (int fanIndex = 0; fanIndex < fanSize / 2; fanIndex++) {
         const size_t vi0 = fullModelBaseTriangleInterfaceVi
                            + fanIndex * fullPatternInterfaceVertexOffset;
         const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset;
         assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN());
-        fullPatternOppositeInterfaceViMap.emplace_back(std::make_pair(vi0, vi1));
+        fullPatternOppositeInterfaceViPairs.emplace_back(std::make_pair(vi0, vi1));
     }
+    global.fullPatternInterfaceViPairs = fullPatternOppositeInterfaceViPairs;
 
+#if POLYSCOPE_DEFINED
     const bool debugMapping = false;
     if (debugMapping) {
-#if POLYSCOPE_DEFINED
         reducedPattern.registerForDrawing();
 
         //    std::vector<glm::vec3> colors_reducedPatternExcludedEdges(
@@ -361,8 +384,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 :
+        for (const std::pair<size_t, size_t> oppositeVerts : fullPatternOppositeInterfaceViPairs) {
-             fullPatternOppositeInterfaceViMap) {
             auto color = polyscope::getNextUniqueColor();
             nodeColorsOpposite[oppositeVerts.first] = color;
             nodeColorsOpposite[oppositeVerts.second] = color;
@@ -395,8 +417,8 @@ void ReducedModelOptimizer::computeMaps(
                                    nodeColorsReducedToFull_full)
             ->setEnabled(true);
         polyscope::show();
-#endif
     }
+#endif
 }
 
 void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern,
@@ -462,12 +484,14 @@ const double I3=global.initialParameters(4) * x[4];
             e.I3 = I3;
     }
     assert(pReducedPatternSimulationMesh->EN() == 12);
-    assert(n>=2);
+    assert(n >= 2);
 
-    CoordType center_barycentric(1, 0, 0);
+    //    CoordType center_barycentric(1, 0, 0);
-    CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5);
+    //    CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5);
-    CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric)
+    //    CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric)
-                                        * x[n - 2]);
+    //                                        * x[n - 2]);
+
+    CoordType movableVertex_barycentric(1 - x[n - 2], x[n - 2] / 2, x[n - 2] / 2);
     CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0)
                                                       * movableVertex_barycentric[0]
                                                   + global.baseTriangle.cP(1)
@@ -480,7 +504,7 @@ const double I3=global.initialParameters(4) * x[4];
 
     for (int rotationCounter = 0;
          rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {
-        pReducedPatternSimulationMesh->vert[2 * rotationCounter + 1].P()
+        pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
             = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
                                   vcg::math::ToRad(60.0 * rotationCounter))
               * baseTriangleMovableVertexPosition;
@@ -690,26 +714,33 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults(
 
         assert(global.minX[xVariableIndex] == optimizationResult_dlib.x(xVariableIndex));
         optimalX[xVariableIndex] = optimizationResult_dlib.x(xVariableIndex);
+        std::cout << results.optimalXNameValuePairs[xVariableIndex].first << ":"
+                  << optimalX[xVariableIndex] << " ";
     }
+    std::cout << std::endl;
 
     // Compute obj value per simulation scenario and the raw objective value
     updateMesh(optimalX.size(), optimalX.data());
     //    results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios);
+    //TODO:use push_back it will make the code more readable
     LinearSimulationModel simulator;
     results.objectiveValue.totalRaw = 0;
     results.objectiveValue.perSimulationScenario_translational.resize(
-        totalNumberOfSimulationScenarios);
+        global.simulationScenarioIndices.size());
     results.objectiveValue.perSimulationScenario_rawTranslational.resize(
-        totalNumberOfSimulationScenarios);
+        global.simulationScenarioIndices.size());
-    results.objectiveValue.perSimulationScenario_rotational.resize(totalNumberOfSimulationScenarios);
+    results.objectiveValue.perSimulationScenario_rotational.resize(
+        global.simulationScenarioIndices.size());
     results.objectiveValue.perSimulationScenario_rawRotational.resize(
-        totalNumberOfSimulationScenarios);
+        global.simulationScenarioIndices.size());
-    results.objectiveValue.perSimulationScenario_total.resize(totalNumberOfSimulationScenarios);
+    results.objectiveValue.perSimulationScenario_total.resize(
-    for (int simulationScenarioIndex : global.simulationScenarioIndices) {
+        global.simulationScenarioIndices.size());
+    for (int i = 0; i < global.simulationScenarioIndices.size(); i++) {
+        const int simulationScenarioIndex = global.simulationScenarioIndices[i];
         SimulationResults reducedModelResults = simulator.executeSimulation(
             global.reducedPatternSimulationJobs[simulationScenarioIndex]);
 
-        results.objectiveValue.perSimulationScenario_total[simulationScenarioIndex] = computeError(
+        results.objectiveValue.perSimulationScenario_total[i] = computeError(
             global.fullPatternResults[simulationScenarioIndex],
             reducedModelResults,
             global.reducedToFullInterfaceViMap,
@@ -721,15 +752,13 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults(
             global.fullPatternResults[simulationScenarioIndex].displacements,
             reducedModelResults.displacements,
             global.reducedToFullInterfaceViMap);
-        results.objectiveValue.perSimulationScenario_rawTranslational[simulationScenarioIndex]
+        results.objectiveValue.perSimulationScenario_rawTranslational[i] = rawTranslationalError;
-            = rawTranslationalError;
         //Raw rotational
         const double rawRotationalError = computeRawRotationalError(
             global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion,
             reducedModelResults.rotationalDisplacementQuaternion,
             global.reducedToFullInterfaceViMap);
-        results.objectiveValue.perSimulationScenario_rawRotational[simulationScenarioIndex]
+        results.objectiveValue.perSimulationScenario_rawRotational[i] = rawRotationalError;
-            = rawRotationalError;
 
         //Normalized translational
         const double normalizedTranslationalError = computeDisplacementError(
@@ -737,8 +766,7 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults(
             reducedModelResults.displacements,
             global.reducedToFullInterfaceViMap,
             global.translationalDisplacementNormalizationValues[simulationScenarioIndex]);
-        results.objectiveValue.perSimulationScenario_translational[simulationScenarioIndex]
+        results.objectiveValue.perSimulationScenario_translational[i] = normalizedTranslationalError;
-            = normalizedTranslationalError;
         //            const double test_normalizedTranslationError = computeDisplacementError(
         //                global.fullPatternResults[simulationScenarioIndex].displacements,
         //                reducedModelResults.displacements,
@@ -750,8 +778,7 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults(
             reducedModelResults.rotationalDisplacementQuaternion,
             global.reducedToFullInterfaceViMap,
             global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
-        results.objectiveValue.perSimulationScenario_rotational[simulationScenarioIndex]
+        results.objectiveValue.perSimulationScenario_rotational[i] = normalizedRotationalError;
-            = normalizedRotationalError;
         //            const double test_normalizedRotationalError = computeRotationalError(
         //                global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion,
         //                reducedModelResults.rotationalDisplacementQuaternion,
@@ -766,16 +793,15 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults(
         std::cout << "translation normalization value:"
                   << global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
                   << std::endl;
-        std::cout << "Translational error:" << normalizedTranslationalError << std::endl;
         std::cout << "raw Rotational error:" << rawRotationalError << std::endl;
         std::cout << "rotational normalization value:"
                   << global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
                   << std::endl;
+        std::cout << "Translational error:" << normalizedTranslationalError << std::endl;
         std::cout << "Rotational error:" << normalizedRotationalError << std::endl;
         //        results.objectiveValuePerSimulationScenario[simulationScenarioIndex]
         //            = normalizedTranslationalError + normalizedRotationalError;
-        std::cout << "Objective value:"
+        std::cout << "Total Error value:" << results.objectiveValue.perSimulationScenario_total[i]
-                  << results.objectiveValue.perSimulationScenario_total[simulationScenarioIndex]
                   << std::endl;
         results.objectiveValue.totalRaw += rawTranslationalError + rawRotationalError;
         std::cout << std::endl;
@@ -786,6 +812,9 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults(
         std::cout << "Finished optimizing." << endl;
     std::cout << "Total optimal objective value:" << results.objectiveValue.total << std::endl;
     assert(global.minY == optimizationResult_dlib.y);
+    if (global.minY != optimizationResult_dlib.y) {
+        std::cerr << "ERROR in objective value" << std::endl;
+    }
     }
 
     for (int simulationScenarioIndex : global.simulationScenarioIndices) {
@@ -830,19 +859,239 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::runOptimization(const
     return results;
 }
 
-//TODO: create a function that takes as arguments the magnitude range, a lambda that generates the simulation job and the base sim scenario index
+void ReducedModelOptimizer::constructAxialSimulationScenario(
+    const double &forceMagnitude,
+    const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+        &oppositeInterfaceViPairs,
+    SimulationJob &job)
+{
+    for (auto viPairIt = oppositeInterfaceViPairs.begin();
+         viPairIt != oppositeInterfaceViPairs.end();
+         viPairIt++) {
+        if (viPairIt != oppositeInterfaceViPairs.begin()) {
+            CoordType forceDirection(1, 0, 0);
+            job.nodalExternalForces[viPairIt->first]
+                = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
+                  * forceMagnitude;
+            job.constrainedVertices[viPairIt->second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
+        }
+    }
+
+}
+
+void ReducedModelOptimizer::constructShearSimulationScenario(
+    const double &forceMagnitude,
+    const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+        &oppositeInterfaceViPairs,
+    SimulationJob &job)
+{
+    for (auto viPairIt = oppositeInterfaceViPairs.begin();
+         viPairIt != oppositeInterfaceViPairs.end();
+         viPairIt++) {
+        if (viPairIt != oppositeInterfaceViPairs.begin()) {
+            CoordType forceDirection(0, 1, 0);
+            const auto viPair = *viPairIt;
+            job.nodalExternalForces[viPair.first]
+                = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
+                  * forceMagnitude;
+            job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
+        }
+    }
+}
+
+void ReducedModelOptimizer::constructBendingSimulationScenario(
+    const double &forceMagnitude,
+    const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+        &oppositeInterfaceViPairs,
+    SimulationJob &job)
+{
+    for (const auto &viPair : oppositeInterfaceViPairs) {
+        job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 1, 0, 0, 0}) * forceMagnitude;
+        job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
+    }
+}
+
+void ReducedModelOptimizer::constructDomeSimulationScenario(
+    const double &forceMagnitude,
+    const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+        &oppositeInterfaceViPairs,
+    SimulationJob &job)
+{
+    for (auto viPairIt = oppositeInterfaceViPairs.begin();
+         viPairIt != oppositeInterfaceViPairs.end();
+         viPairIt++) {
+        const auto viPair = *viPairIt;
+        CoordType interfaceVector = (job.pMesh->vert[viPair.first].cP()
+                                     - job.pMesh->vert[viPair.second].cP());
+        VectorType momentAxis = vcg::RotationMatrix(VectorType(0, 0, 1), vcg::math::ToRad(90.0))
+                                * interfaceVector.Normalize();
+        if (viPairIt == oppositeInterfaceViPairs.begin()) {
+            job.nodalForcedDisplacements[viPair.first]
+                = Eigen::Vector3d(-interfaceVector[0],
+                                  -interfaceVector[1],
+                                  0)
+                  * 0.001
+                  * std::abs(
+                      forceMagnitude); //NOTE:Should the forced displacement change relatively to the magnitude?
+            //                  * std::abs(forceMagnitude / maxForceMagnitude_dome);
+            job.nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceVector[0],
+                                                                          interfaceVector[1],
+                                                                          0)
+                                                          * 0.001 * std::abs(forceMagnitude);
+            //                  * std::abs(forceMagnitude / maxForceMagnitude_dome);
+            //            CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
+            //                          ^ CoordType(0, 0, -1).Normalize();
+            //            nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude
+            //                                        * 0.0001;
+            //            nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude
+            //                                         * 0.0001;
+        } else {
+            job.nodalExternalForces[viPair.first]
+                = Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]}) * forceMagnitude
+                  / 5;
+            job.nodalExternalForces[viPair.second]
+                = Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]})
+                  * forceMagnitude / 5;
+            job.constrainedVertices[viPair.first] = std::unordered_set<DoFType>{2};
+            job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{2};
+        }
+    }
+}
+
+void ReducedModelOptimizer::constructSaddleSimulationScenario(
+    const double &forceMagnitude,
+    const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+        &oppositeInterfaceViPairs,
+    SimulationJob &job)
+{
+    for (auto viPairIt = oppositeInterfaceViPairs.begin();
+         viPairIt != oppositeInterfaceViPairs.end();
+         viPairIt++) {
+        const auto &viPair = *viPairIt;
+        CoordType v = (job.pMesh->vert[viPair.first].cP() - job.pMesh->vert[viPair.second].cP())
+                      ^ CoordType(0, 0, -1).Normalize();
+        if (viPairIt == oppositeInterfaceViPairs.begin()) {
+            job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0})
+                                                    * forceMagnitude;
+            job.nodalExternalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0})
+                                                     * forceMagnitude;
+        } else {
+            job.constrainedVertices[viPair.first] = std::unordered_set<DoFType>{2};
+            job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
+
+            job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0})
+                                                    * forceMagnitude / 2;
+            job.nodalExternalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0})
+                                                     * forceMagnitude / 2;
+        }
+    }
+}
+
+double fullPatternMaxSimulationForceObjective(const double &forceMagnitude)
+{
+    const double &desiredDisplacementValue = global.desiredMaxDisplacementValue;
+    SimulationJob job;
+    job.pMesh = global.pFullPatternSimulationMesh;
+    global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job);
+    //    ReducedModelOptimizer::axialSimulationScenario(forceMagnitude,
+    //                                                   global.fullPatternInterfaceViPairs,
+    //                                                   job);
+
+    DRMSimulationModel simulator;
+    DRMSimulationModel::Settings settings;
+    settings.totalExternalForcesNormPercentageTermination = 1e-2;
+    //    settings.shouldUseTranslationalKineticEnergyThreshold = true;
+    //    settings.totalTranslationalKineticEnergyThreshold = 1e-7;
+    SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
+                                                            settings);
+    const double error = std::abs(
+        //        results.displacements[global.fullPatternInterfaceViPairs[1].first].getTranslation().norm()
+        results.displacements[global.interfaceViForComputingScenarioError].getTranslation().norm()
+        - desiredDisplacementValue);
+
+    return error;
+}
+
+double ReducedModelOptimizer::getFullPatternMaxSimulationForce(const BaseSimulationScenario &scenario)
+{
+    double forceMagnitude = 1;
+    global.desiredMaxDisplacementValue = 0.1
+                                         * vcg::Distance(global.baseTriangle.cP(0),
+                                                         (global.baseTriangle.cP(1)
+                                                          + global.baseTriangle.cP(2))
+                                                             / 2);
+    const double optimizationEpsilon = global.desiredMaxDisplacementValue * 0.5;
+    switch (scenario) {
+    case Axial:
+        global.constructScenarioFunction = &ReducedModelOptimizer::constructAxialSimulationScenario;
+        global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
+        dlib::find_min_single_variable(&fullPatternMaxSimulationForceObjective,
+                                       forceMagnitude,
+                                       1e-8,
+                                       1e8,
+                                       optimizationEpsilon);
+        break;
+    case Shear:
+        global.constructScenarioFunction = &ReducedModelOptimizer::constructShearSimulationScenario;
+        global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
+        dlib::find_min_single_variable(&fullPatternMaxSimulationForceObjective,
+                                       forceMagnitude,
+                                       1e-8,
+                                       1e8,
+                                       optimizationEpsilon);
+        break;
+    case Bending:
+        global.constructScenarioFunction = &ReducedModelOptimizer::constructBendingSimulationScenario;
+        global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first;
+        dlib::find_min_single_variable(&fullPatternMaxSimulationForceObjective,
+                                       forceMagnitude,
+                                       1e-8,
+                                       1e8,
+                                       optimizationEpsilon,
+                                       200);
+        break;
+    case Dome:
+        global.desiredMaxDisplacementValue *= 2;
+        global.constructScenarioFunction = &ReducedModelOptimizer::constructDomeSimulationScenario;
+        global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
+        dlib::find_min_single_variable(&fullPatternMaxSimulationForceObjective,
+                                       forceMagnitude,
+                                       1e-8,
+                                       1e8,
+                                       optimizationEpsilon,
+                                       200);
+        break;
+    case Saddle:
+        global.desiredMaxDisplacementValue *= 2;
+        global.constructScenarioFunction = &ReducedModelOptimizer::constructSaddleSimulationScenario;
+        global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first;
+        dlib::find_min_single_variable(&fullPatternMaxSimulationForceObjective,
+                                       forceMagnitude,
+                                       1e-8,
+                                       1e8,
+                                       optimizationEpsilon,
+                                       150);
+        break;
+    }
+
+    return forceMagnitude;
+}
+
+//TODO: Make more compact
 std::vector<std::shared_ptr<SimulationJob>>
 ReducedModelOptimizer::createFullPatternSimulationScenarios(
     const std::shared_ptr<SimulationMesh> &pMesh)
 {
     std::vector<std::shared_ptr<SimulationJob>> scenarios;
     scenarios.resize(totalNumberOfSimulationScenarios);
-    std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
+
-    std::unordered_map<VertexIndex, Vector6d> nodalForces;
+    SimulationJob job;
+    job.pMesh = pMesh;
 
     //// Axial
-    const double maxForceMagnitude_axial = 500;
+    const double maxForceMagnitude_axial = getFullPatternMaxSimulationForce(
-    const double minForceMagnitude_axial = -500;
+        BaseSimulationScenario::Axial);
+    const double minForceMagnitude_axial = -maxForceMagnitude_axial;
     const int numberOfSimulationScenarios_axial
         = simulationScenariosResolution[BaseSimulationScenario::Axial];
     const double forceMagnitudeStep_axial = numberOfSimulationScenarios_axial == 1
@@ -853,38 +1102,27 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios(
         = std::accumulate(simulationScenariosResolution.begin(),
                           simulationScenariosResolution.begin() + BaseSimulationScenario::Axial,
                           0);
-
     for (int axialSimulationScenarioIndex = 0;
          axialSimulationScenarioIndex < numberOfSimulationScenarios_axial;
          axialSimulationScenarioIndex++) {
+        job.nodalExternalForces.clear();
+        job.constrainedVertices.clear();
+        job.nodalForcedDisplacements.clear();
+        job.label = baseSimulationScenarioNames[BaseSimulationScenario::Axial] + "_"
+                    + std::to_string(axialSimulationScenarioIndex);
+
         const double forceMagnitude = (forceMagnitudeStep_axial * axialSimulationScenarioIndex
                                        + minForceMagnitude_axial);
-        for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
+        constructAxialSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job);
-             viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
-             viPairIt++) {
-            if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) {
-                CoordType forceDirection(1, 0, 0);
-                const auto viPair = *viPairIt;
-                nodalForces[viPair.first]
-                    = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
-                      * forceMagnitude;
-                fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
-            }
-        }
 
         scenarios[baseSimulationScenarioIndexOffset_axial + axialSimulationScenarioIndex]
-            = std::make_shared<SimulationJob>(
+            = std::make_shared<SimulationJob>(job);
-                SimulationJob(pMesh,
-                              baseSimulationScenarioNames[BaseSimulationScenario::Axial] + "_"
-                                  + std::to_string(axialSimulationScenarioIndex),
-                              fixedVertices,
-                              nodalForces,
-                              {}));
     }
 
     //// Shear
-    const double maxForceMagnitude_shear = 50;
+    const double maxForceMagnitude_shear = getFullPatternMaxSimulationForce(
-    const double minForceMagnitude_shear = -50;
+        BaseSimulationScenario::Shear);
+    const double minForceMagnitude_shear = -maxForceMagnitude_shear;
     const int numberOfSimulationScenarios_shear
         = simulationScenariosResolution[BaseSimulationScenario::Shear];
     const double forceMagnitudeStep_shear = numberOfSimulationScenarios_shear == 1
@@ -898,35 +1136,21 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios(
     for (int shearSimulationScenarioIndex = 0;
          shearSimulationScenarioIndex < numberOfSimulationScenarios_shear;
          shearSimulationScenarioIndex++) {
-        fixedVertices.clear();
+        job.constrainedVertices.clear();
-        nodalForces.clear();
+        job.nodalExternalForces.clear();
+        job.label = baseSimulationScenarioNames[BaseSimulationScenario::Shear] + "_"
+                    + std::to_string(shearSimulationScenarioIndex);
         const double forceMagnitude = (forceMagnitudeStep_shear * shearSimulationScenarioIndex
                                        + minForceMagnitude_shear);
-        for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
+        constructShearSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job);
-             viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
-             viPairIt++) {
-            if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) {
-                CoordType forceDirection(0, 1, 0);
-                const auto viPair = *viPairIt;
-                nodalForces[viPair.first]
-                    = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
-                      * forceMagnitude;
-                fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
-            }
-        }
         scenarios[baseSimulationScenarioIndexOffset_shear + shearSimulationScenarioIndex]
-            = std::make_shared<SimulationJob>(
+            = std::make_shared<SimulationJob>(job);
-                SimulationJob(pMesh,
-                              baseSimulationScenarioNames[BaseSimulationScenario::Shear] + "_"
-                                  + std::to_string(shearSimulationScenarioIndex),
-                              fixedVertices,
-                              nodalForces,
-                              {}));
     }
 
     ////  Bending
-    const double maxForceMagnitude_bending = 0.005;
+    const double maxForceMagnitude_bending = getFullPatternMaxSimulationForce(
-    const double minForceMagnitude_bending = -0.005;
+        BaseSimulationScenario::Bending);
+    const double minForceMagnitude_bending = -maxForceMagnitude_bending;
     const int numberOfSimulationScenarios_bending
         = simulationScenariosResolution[BaseSimulationScenario::Bending];
     const double forceMagnitudeStep_bending = numberOfSimulationScenarios_bending == 1
@@ -941,27 +1165,23 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios(
     for (int bendingSimulationScenarioIndex = 0;
          bendingSimulationScenarioIndex < numberOfSimulationScenarios_bending;
          bendingSimulationScenarioIndex++) {
-        fixedVertices.clear();
+        job.nodalExternalForces.clear();
-        nodalForces.clear();
+        job.constrainedVertices.clear();
+        job.label = baseSimulationScenarioNames[BaseSimulationScenario::Bending] + "_"
+                    + std::to_string(bendingSimulationScenarioIndex);
         const double forceMagnitude = (forceMagnitudeStep_bending * bendingSimulationScenarioIndex
                                        + minForceMagnitude_bending);
-        for (const auto &viPair : m_fullPatternOppositeInterfaceViPairs) {
+        constructBendingSimulationScenario(forceMagnitude,
-            nodalForces[viPair.first] = Vector6d({0, 0, 1, 0, 0, 0}) * forceMagnitude;
+                                           m_fullPatternOppositeInterfaceViPairs,
-            fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
+                                           job);
-        }
         scenarios[baseSimulationScenarioIndexOffset_bending + bendingSimulationScenarioIndex]
-            = std::make_shared<SimulationJob>(
+            = std::make_shared<SimulationJob>(job);
-                SimulationJob(pMesh,
-                              baseSimulationScenarioNames[BaseSimulationScenario::Bending] + "_"
-                                  + std::to_string(bendingSimulationScenarioIndex),
-                              fixedVertices,
-                              nodalForces,
-                              {}));
     }
 
     //// Dome
-    const double maxForceMagnitude_dome = 0.025;
+    const double maxForceMagnitude_dome = getFullPatternMaxSimulationForce(
-    const double minForceMagnitude_dome = -0.025;
+        BaseSimulationScenario::Dome);
+    const double minForceMagnitude_dome = -maxForceMagnitude_dome;
     const int numberOfSimulationScenarios_dome
         = simulationScenariosResolution[BaseSimulationScenario::Dome];
     const double forceMagnitudeStep_dome = numberOfSimulationScenarios_dome == 1
@@ -975,58 +1195,22 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios(
     for (int domeSimulationScenarioIndex = 0;
          domeSimulationScenarioIndex < numberOfSimulationScenarios_dome;
          domeSimulationScenarioIndex++) {
-        fixedVertices.clear();
+        job.constrainedVertices.clear();
-        nodalForces.clear();
+        job.nodalExternalForces.clear();
-        std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
+        job.nodalForcedDisplacements.clear();
+        job.label = baseSimulationScenarioNames[BaseSimulationScenario::Dome] + "_"
+                    + std::to_string(domeSimulationScenarioIndex);
         const double forceMagnitude = (forceMagnitudeStep_dome * domeSimulationScenarioIndex
                                        + minForceMagnitude_dome);
-        for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
+        constructDomeSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job);
-             viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
-             viPairIt++) {
-            const auto viPair = *viPairIt;
-            CoordType interfaceVector = (pMesh->vert[viPair.first].cP()
-                                         - pMesh->vert[viPair.second].cP());
-            VectorType momentAxis = vcg::RotationMatrix(VectorType(0, 0, 1), vcg::math::ToRad(90.0))
-                                    * interfaceVector.Normalize();
-            if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) {
-                nodalForcedDisplacements[viPair.first] = Eigen::Vector3d(-interfaceVector[0],
-                                                                         -interfaceVector[1],
-                                                                         0)
-                                                         * std::abs(forceMagnitude);
-                nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceVector[0],
-                                                                          interfaceVector[1],
-                                                                          0)
-                                                          * std::abs(forceMagnitude);
-                //            CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
-                //                          ^ CoordType(0, 0, -1).Normalize();
-                //            nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude
-                //                                        * 0.0001;
-                //            nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude
-                //                                         * 0.0001;
-            } else {
-                nodalForces[viPair.first]
-                    = Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]})
-                      * forceMagnitude / 5;
-                nodalForces[viPair.second]
-                    = Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]})
-                      * forceMagnitude / 5;
-                fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
-                fixedVertices[viPair.second] = std::unordered_set<DoFType>{2};
-            }
-        }
         scenarios[baseSimulationScenarioIndexOffset_dome + domeSimulationScenarioIndex]
-            = std::make_shared<SimulationJob>(
+            = std::make_shared<SimulationJob>(job);
-                SimulationJob(pMesh,
-                              baseSimulationScenarioNames[BaseSimulationScenario::Dome] + "_"
-                                  + std::to_string(domeSimulationScenarioIndex),
-                              fixedVertices,
-                              nodalForces,
-                              nodalForcedDisplacements));
     }
 
     //// Saddle
-    const double maxForceMagnitude_saddle = 0.005;
+    const double maxForceMagnitude_saddle = getFullPatternMaxSimulationForce(
-    const double minForceMagnitude_saddle = -0.005;
+        BaseSimulationScenario::Saddle);
+    const double minForceMagnitude_saddle = -maxForceMagnitude_saddle;
     const int numberOfSimulationScenarios_saddle
         = simulationScenariosResolution[BaseSimulationScenario::Saddle];
     const double forceMagnitudeStep_saddle = numberOfSimulationScenarios_saddle == 1
@@ -1041,37 +1225,18 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios(
     for (int saddleSimulationScenarioIndex = 0;
          saddleSimulationScenarioIndex < numberOfSimulationScenarios_saddle;
          saddleSimulationScenarioIndex++) {
-        fixedVertices.clear();
+        job.constrainedVertices.clear();
-        nodalForces.clear();
+        job.nodalExternalForces.clear();
+        job.nodalForcedDisplacements.clear();
+        job.label = baseSimulationScenarioNames[BaseSimulationScenario::Saddle] + "_"
+                    + std::to_string(saddleSimulationScenarioIndex);
         const double forceMagnitude = (forceMagnitudeStep_saddle * saddleSimulationScenarioIndex
                                        + minForceMagnitude_saddle);
-        for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
+        constructSaddleSimulationScenario(forceMagnitude,
-             viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
+                                          m_fullPatternOppositeInterfaceViPairs,
-             viPairIt++) {
+                                          job);
-            const auto &viPair = *viPairIt;
-            CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
-                          ^ CoordType(0, 0, -1).Normalize();
-            if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) {
-                nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude;
-                nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude;
-            } else {
-                fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
-                fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
-
-                nodalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude
-                                            / 2;
-                nodalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude
-                                             / 2;
-            }
-        }
         scenarios[baseSimulationScenarioIndexOffset_saddle + saddleSimulationScenarioIndex]
-            = std::make_shared<SimulationJob>(
+            = std::make_shared<SimulationJob>(job);
-                SimulationJob(pMesh,
-                              baseSimulationScenarioNames[BaseSimulationScenario::Saddle] + "_"
-                                  + std::to_string(saddleSimulationScenarioIndex),
-                              fixedVertices,
-                              nodalForces,
-                              {}));
     }
 
     return scenarios;
@@ -1181,32 +1346,38 @@ Results ReducedModelOptimizer::optimize(
     global.numOfSimulationCrashes = 0;
     global.numberOfFunctionCalls = 0;
     global.optimizationSettings = optimizationSettings;
+    global.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh;
     global.fullPatternSimulationJobs = createFullPatternSimulationScenarios(
         m_pFullPatternSimulationMesh);
     //  polyscope::removeAllStructures();
 
     DRMSimulationModel::Settings simulationSettings;
     simulationSettings.shouldDraw = false;
-    //    global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
+    const bool drawFullPatternSimulationResults = false;
-    //        ReducedPatternOptimization::Colors::fullInitial);
+    if (drawFullPatternSimulationResults) {
+        global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
+            ReducedPatternOptimization::Colors::fullInitial);
+    }
     //    LinearSimulationModel linearSimulator;
     for (int simulationScenarioIndex : global.simulationScenarioIndices) {
         const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
             = global.fullPatternSimulationJobs[simulationScenarioIndex];
         SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
                                                                            simulationSettings);
+        if (drawFullPatternSimulationResults) {
             //                SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation(
             //                    pFullPatternSimulationJob);
-        //        fullPatternResults.registerForDrawing(ReducedPatternOptimization::Colors::fullDeformed,
+            fullPatternResults.registerForDrawing(ReducedPatternOptimization::Colors::fullDeformed,
-        //                                              true,
+                                                  true,
-        //                                              true);
+                                                  true);
             //                fullPatternResults_linear.labelPrefix += "_linear";
             //                fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed,
             //                                                             true,
             //                                                             true);
-        //        polyscope::show();
+            polyscope::show();
-        //        fullPatternResults.unregister();
+            fullPatternResults.unregister();
             //                fullPatternResults_linear.unregister();
+        }
         global.fullPatternResults[simulationScenarioIndex] = fullPatternResults;
         SimulationJob reducedPatternSimulationJob;
         reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
@@ -1217,12 +1388,14 @@ Results ReducedModelOptimizer::optimize(
             = std::make_shared<SimulationJob>(reducedPatternSimulationJob);
         //        std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl;
     }
-    //    global.fullPatternSimulationJobs[0]->pMesh->unregister();
+    if (drawFullPatternSimulationResults) {
-
+        global.fullPatternSimulationJobs[0]->pMesh->unregister();
+    }
     //    if (global.optimizationSettings.normalizationStrategy
     //        != Settings::NormalizationStrategy::NonNormalized) {
     computeObjectiveValueNormalizationFactors();
     //    }
+    PolyscopeInterface::deinitPolyscope();
     Results optResults = runOptimization(optimizationSettings);
 
     return optResults;

src/reducedmodeloptimizer.hpp

@@ -30,12 +30,13 @@ class ReducedModelOptimizer
     std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
 
 public:
-    constexpr static std::array<int, 5> simulationScenariosResolution = {2, 2, 6, 6, 6};
+    constexpr static std::array<int, 5> simulationScenariosResolution = {4, 4, 4, 4, 4};
     inline static int totalNumberOfSimulationScenarios
         = std::accumulate(simulationScenariosResolution.begin(),
                           simulationScenariosResolution.end(),
                           0);
     inline static int fanSize{6};
+    inline static double initialHexagonSize{0.3};
     inline static VectorType patternPlaneNormal{0, 0, 1};
     ReducedPatternOptimization::Results optimize(
         const ReducedPatternOptimization::Settings &xRanges,
@@ -132,6 +133,35 @@ public:
             &reducedToFullInterfaceViMap,
         const double &normalizationFactor_translationalDisplacement,
         const double &normalizationFactor_rotationalDisplacement);
+    static void constructAxialSimulationScenario(
+        const double &forceMagnitude,
+        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            &oppositeInterfaceViPairs,
+        SimulationJob &job);
+
+    static void constructShearSimulationScenario(
+        const double &forceMagnitude,
+        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            &oppositeInterfaceViPairs,
+        SimulationJob &job);
+
+    static void constructBendingSimulationScenario(
+        const double &forceMagnitude,
+        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            &oppositeInterfaceViPairs,
+        SimulationJob &job);
+
+    static void constructDomeSimulationScenario(
+        const double &forceMagnitude,
+        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            &oppositeInterfaceViPairs,
+        SimulationJob &job);
+
+    static void constructSaddleSimulationScenario(
+        const double &forceMagnitude,
+        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
+            &oppositeInterfaceViPairs,
+        SimulationJob &job);
 
 private:
     static void computeDesiredReducedModelDisplacements(
@@ -153,6 +183,9 @@ private:
     static ReducedPatternOptimization::Results getResults(
         const dlib::function_evaluation &optimizationResult_dlib,
         const ReducedPatternOptimization::Settings &settings);
+    std::vector<double> getFullPatternMaxSimulationForces();
+    double getFullPatternMaxSimulationForce(
+        const ReducedPatternOptimization::BaseSimulationScenario &scenario);
 };
 void updateMesh(long n, const double *x);
 #endif // REDUCEDMODELOPTIMIZER_HPP