Modular Variable comparison

src/main.cpp

@@ -38,20 +38,21 @@ int main(int argc, char *argv[]) {
   reducedPattern.scale(0.03, interfaceNodeIndex);
 
   // Set the optization settings
-  ReducedPatternOptimization::xRange beamE{"E", 0.001, 100};
+  ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000};
   ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000};
+  ReducedPatternOptimization::xRange beamI{"I", 0.001, 1000};
   ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000};
   ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000};
   ReducedPatternOptimization::xRange beamJ{"J", 0.001, 1000};
   ReducedPatternOptimization::xRange innerHexagonSize{"HexSize", 0.05, 0.95};
   ReducedPatternOptimization::xRange innerHexagonAngle{"HexAngle", -30.0, 30.0};
   ReducedPatternOptimization::Settings settings_optimization;
-  //  settings_optimization.xRanges
+  settings_optimization.xRanges
-  //      = {beamE, beamA, beamJ, beamI2, beamI3, innerHexagonSize, innerHexagonAngle};
+      = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle};
-  settings_optimization.xRanges = {beamE, beamA, beamI2, innerHexagonSize, innerHexagonAngle};
   const bool input_numberOfFunctionCallsDefined = argc >= 4;
-  settings_optimization.numberOfFunctionCalls =
+  settings_optimization.numberOfFunctionCalls = input_numberOfFunctionCallsDefined
-      input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100;
+                                                    ? std::atoi(argv[3])
+                                                    : 100;
   settings_optimization.normalizationStrategy
       = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
 
@@ -116,9 +117,7 @@ int main(int argc, char *argv[]) {
       const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
       assert(interfaceNodeIndex == numberOfNodesPerSlot[0] + numberOfNodesPerSlot[3]);
       ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
-      optimizer.initializePatterns(fullPattern,
+      optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges);
-                                   reducedPattern,
-                                   settings_optimization.xRanges.size());
       optimizer.optimize(settings_optimization, optimizationResults);
       optimizationResults.label = optimizationName;
       optimizationResults.baseTriangleFullPattern.copy(fullPattern);

src/reducedmodeloptimizer.cpp

@@ -20,7 +20,9 @@ struct GlobalOptimizationVariables {
     std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
         fullPatternInterfaceViPairs;
     matplot::line_handle gPlotHandle;
-    std::vector<double> objectiveValueHistory;
+    std::vector<double> objectiveValueHistoryY;
+    std::vector<double> objectiveValueHistoryX;
+    std::vector<double> plotColors;
     Eigen::VectorXd initialParameters;
     std::vector<int> simulationScenarioIndices;
     double minY{DBL_MAX};
@@ -40,6 +42,10 @@ struct GlobalOptimizationVariables {
     double desiredMaxDisplacementValue;
     double desiredMaxRotationAngle;
     std::string currentScenarioName;
+    std::vector<std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                                   const double &newValue)>>
+        updateReducedPatternFunctions_material;
+
 } global;
 
 double ReducedModelOptimizer::computeDisplacementError(
@@ -165,7 +171,7 @@ double ReducedModelOptimizer::objective(const dlib::matrix<double, 0, 1> &x)
     //    polyscope::show();
     //    global.reducedPatternSimulationJobs[0]->pMesh->unregister();
 
-    //  std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
+    //  std::cout << e.axialConstFact|A,I2,I3,J,r,thetaor << " " << e.torsionConstFactor << " "
     //            << e.firstBendingConstFactor << " " <<
     //            e.secondBendingConstFactor
     //            << std::endl;
@@ -249,28 +255,36 @@ double ReducedModelOptimizer::objective(const dlib::matrix<double, 0, 1> &x)
     //  std::cout << error << std::endl;
     if (totalError < global.minY) {
         global.minY = totalError;
-        global.minX.assign(x.begin(), x.begin() + n);
+        std::cout << "New best:" << totalError << std::endl;
+        //        global.minX.assign(x.begin(), x.begin() + n);
+        std::cout.precision(17);
+        for (int i = 0; i < x.size(); i++) {
+            std::cout << x(i) << " ";
+        }
+        std::cout << std::endl;
+        global.objectiveValueHistoryY.push_back(std::log(totalError));
+        global.objectiveValueHistoryX.push_back(global.numberOfFunctionCalls);
+        global.plotColors.push_back(0.1);
+        //        auto xPlot = matplot::linspace(0,
+        //                                       global.objectiveValueHistoryY.size(),
+        //                                       global.objectiveValueHistoryY.size());
+        global.gPlotHandle = matplot::scatter(global.objectiveValueHistoryX,
+                                              global.objectiveValueHistoryY,
+                                              4,
+                                              global.plotColors);
+        //        matplot::show();
+        //                SimulationResultsReporter::createPlot("Number of Steps",
+        //                                              "Objective value",
+        //                                              global.objectiveValueHistoryY);
     }
+//    global.objectiveValueHistory.push_back(totalError);
+//    global.plotColors.push_back(10);
 #ifdef POLYSCOPE_DEFINED
     ++global.numberOfFunctionCalls;
-    global.objectiveValueHistory.push_back(
-        std::sqrt(totalError / global.simulationScenarioIndices.size()));
     if (global.optimizationSettings.numberOfFunctionCalls >= 100
         && global.numberOfFunctionCalls % (global.optimizationSettings.numberOfFunctionCalls / 100)
                == 0) {
         std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl;
-        auto xPlot = matplot::linspace(0,
-                                       global.objectiveValueHistory.size(),
-                                       global.objectiveValueHistory.size());
-        //      std::vector<double> colors(global.gObjectiveValueHistory.size(), 2);
-        //      if (global.g_firstRoundIterationIndex != 0) {
-        //        for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(),
-        //                 [](double &c) { c = 0.7; });
-        //      }
-        global.gPlotHandle = matplot::scatter(xPlot, global.objectiveValueHistory);
-        SimulationResultsReporter::createPlot("Number of Steps",
-                                              "Objective value",
-                                              global.objectiveValueHistory);
     }
 #endif
     //     compute error and return it
@@ -476,7 +490,7 @@ ReducedModelOptimizer::ReducedModelOptimizer(const std::vector<size_t> &numberOf
 
 void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
                                                PatternGeometry &reducedPattern,
-                                               const int &optimizationParameters)
+                                               const std::vector<xRange> &optimizationParameters)
 {
     assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN());
     fullPatternNumberOfEdges = fullPattern.EN();
@@ -492,113 +506,12 @@ void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
 
     computeMaps(copyFullPattern, copyReducedPattern);
     createSimulationMeshes(copyFullPattern, copyReducedPattern);
+    initializeUpdateReducedPatternFunctions();
     initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters);
 }
 
-void ReducedModelOptimizer::initializeOptimizationParameters(
+void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions()
-    const std::shared_ptr<SimulationMesh> &mesh, const int &optimizationParamters)
 {
-    global.numberOfOptimizationParameters = optimizationParamters;
-    global.initialParameters.resize(global.numberOfOptimizationParameters);
-
-    switch (optimizationParamters) {
-    case 2:
-        function_updateReducedPattern_material =
-            [](const dlib::matrix<double, 0, 1> &x,
-               std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {};
-        break;
-
-    case 3:
-        global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
-        function_updateReducedPattern_material =
-            [&](const dlib::matrix<double, 0, 1> &x,
-                std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
-                const double E = global.initialParameters(0) * x(0);
-                for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
-                    Element &e = pReducedPatternSimulationMesh->elements[ei];
-                    e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
-                }
-            };
-        break;
-    case 4:
-        global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
-        global.initialParameters(1) = mesh->elements[0].A;
-        function_updateReducedPattern_material =
-            [&](const dlib::matrix<double, 0, 1> &x,
-                std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
-                const double E = global.initialParameters(0) * x(0);
-                const double A = global.initialParameters(1) * x(1);
-                const double beamWidth = std::sqrt(A);
-                const double beamHeight = beamWidth;
-                for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
-                    Element &e = pReducedPatternSimulationMesh->elements[ei];
-                    e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight));
-                    e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
-                }
-            };
-        break;
-    case 5:
-        global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
-        global.initialParameters(1) = mesh->elements[0].A;
-        global.initialParameters(2) = mesh->elements[0].inertia.I2;
-        function_updateReducedPattern_material =
-            [&](const dlib::matrix<double, 0, 1> &x,
-                std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
-                const double E = global.initialParameters(0) * x(0);
-                const double A = global.initialParameters(1) * x(1);
-                const double beamWidth = std::sqrt(A);
-                const double beamHeight = beamWidth;
-                const double I = global.initialParameters(2) * x(2);
-                for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
-                    Element &e = pReducedPatternSimulationMesh->elements[ei];
-                    e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight));
-                    e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
-                    e.inertia.I2 = I;
-                    e.inertia.I3 = I;
-                }
-            };
-        break;
-    case 7:
-        global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
-        global.initialParameters(1) = mesh->elements[0].A;
-        global.initialParameters(2) = mesh->elements[0].inertia.J;
-        global.initialParameters(3) = mesh->elements[0].inertia.I2;
-        global.initialParameters(4) = mesh->elements[0].inertia.I3;
-        function_updateReducedPattern_material = [&](const dlib::matrix<double, 0, 1> &x,
-                                                     std::shared_ptr<SimulationMesh>
-                                                         &pReducedPatternSimulationMesh) {
-            const double E = global.initialParameters(0) * x(0);
-            const double A = global.initialParameters(1) * x(1);
-            const double beamWidth = std::sqrt(A);
-            const double beamHeight = beamWidth;
-
-            const double J = global.initialParameters(2) * x(2);
-            const double I2 = global.initialParameters(3) * x(3);
-            const double I3 = global.initialParameters(4) * x(4);
-            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
-                Element &e = pReducedPatternSimulationMesh->elements[ei];
-                e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight));
-                e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
-                e.inertia.J = J;
-                e.inertia.I2 = I2;
-                e.inertia.I3 = I3;
-            }
-            //#ifdef POLYSCOPE_DEFINED
-            //    pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
-            //    pReducedPatternSimulationMesh->registerForDrawing();
-            //    std::cout << "Angle:" + std::to_string(x[n - 1]) + " size:" + std::to_string(x[n - 2])
-            //              << std::endl;
-            //    std::cout << "Verts:" << pReducedPatternSimulationMesh->VN() << std::endl;
-            //    polyscope::show();
-            //#endif
-        };
-        break;
-    default:
-        std::cerr << "Wrong number of parameters:" << optimizationParamters << std::endl;
-        std::terminate();
-        break;
-    }
-
     function_updateReducedPattern_geometry = [&](const dlib::matrix<double, 0, 1> &x,
                                                  std::shared_ptr<SimulationMesh>
                                                      &pReducedPatternSimulationMesh) {
@@ -630,6 +543,120 @@ void ReducedModelOptimizer::initializeOptimizationParameters(
                   * baseTriangleMovableVertexPosition;
         }
     };
+
+    function_updateReducedPattern_material_E =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newE) {
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.setMaterial(ElementMaterial(e.material.poissonsRatio, newE));
+            }
+        };
+    function_updateReducedPattern_material_A =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newA) {
+            const double beamWidth = std::sqrt(newA);
+            const double beamHeight = beamWidth;
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight));
+            }
+        };
+    function_updateReducedPattern_material_I =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newI) {
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.inertia.I2 = newI;
+                e.inertia.I3 = newI;
+            }
+        };
+    function_updateReducedPattern_material_I2 =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newI2) {
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.inertia.I2 = newI2;
+            }
+        };
+    function_updateReducedPattern_material_I3 =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newI3) {
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.inertia.I3 = newI3;
+            }
+        };
+    function_updateReducedPattern_material_J =
+        [&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newJ) {
+            for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
+                Element &e = pReducedPatternSimulationMesh->elements[ei];
+                e.inertia.J = newJ;
+            }
+        };
+}
+
+void ReducedModelOptimizer::initializeOptimizationParameters(
+    const std::shared_ptr<SimulationMesh> &mesh, const std::vector<xRange> &optimizationParameters)
+{
+    global.numberOfOptimizationParameters = optimizationParameters.size();
+    global.initialParameters.resize(global.numberOfOptimizationParameters);
+
+    for (int optimizationParameterIndex = 0;
+         optimizationParameterIndex < optimizationParameters.size();
+         optimizationParameterIndex++) {
+        const xRange &xOptimizationParameter = optimizationParameters[optimizationParameterIndex];
+        if (xOptimizationParameter.label == "E") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0]
+                                                                       .material.youngsModulus;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_E);
+        } else if (xOptimizationParameter.label == "A") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0].A;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_A);
+        } else if (xOptimizationParameter.label == "I") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I2;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_I);
+        } else if (xOptimizationParameter.label == "I2") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I2;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_I2);
+        } else if (xOptimizationParameter.label == "I3") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I3;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_I3);
+        } else if (xOptimizationParameter.label == "J") {
+            global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.J;
+            global.updateReducedPatternFunctions_material.push_back(
+                function_updateReducedPattern_material_J);
+        }
+        //NOTE:Assuming that I2,I3 and J are be passed after E and A because the update of E and A changes I2,I3,J.
+        //        case "HexSize":
+        //            updateReducedPatternFunctions_material.push_back(
+        //                function_updateReducedPattern_material_E);
+        //            break;
+        //        case "HexAngle":
+        //            updateReducedPatternFunctions_material.push_back(
+        //                function_updateReducedPattern_material_E);
+        //            break;
+    }
+
+    if (global.updateReducedPatternFunctions_material.size() != 0) {
+        function_updateReducedPattern_material =
+            [&](const dlib::matrix<double, 0, 1> &x,
+                std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+                for (int optimizationParameterIndex = 0;
+                     optimizationParameterIndex
+                     < global.updateReducedPatternFunctions_material.size();
+                     optimizationParameterIndex++) {
+                    global.updateReducedPatternFunctions_material[optimizationParameterIndex](
+                        pReducedPatternSimulationMesh,
+                        global.initialParameters(optimizationParameterIndex)
+                            * x(optimizationParameterIndex));
+                }
+            };
+    } else {
+        function_updateReducedPattern_material =
+            [](const dlib::matrix<double, 0, 1> &x,
+               std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {};
+    }
 }
 
 void ReducedModelOptimizer::computeReducedModelSimulationJob(
@@ -795,16 +822,15 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
         if (xVariableIndex < settings.xRanges.size() - 2) {
             results.optimalXNameValuePairs[xVariableIndex]
                 = std::make_pair(settings.xRanges[xVariableIndex].label,
-                                 global.minX[xVariableIndex]
+                                 optimizationResult_dlib.x(xVariableIndex)
                                      * global.initialParameters(xVariableIndex));
         } else {
             //Hex size and angle are pure values (not multipliers of the initial values)
             results.optimalXNameValuePairs[xVariableIndex]
                 = std::make_pair(settings.xRanges[xVariableIndex].label,
-                                 global.minX[xVariableIndex]);
+                                 optimizationResult_dlib.x(xVariableIndex));
         }
 
-        assert(global.minX[xVariableIndex] == optimizationResult_dlib.x(xVariableIndex));
         optimalX[xVariableIndex] = optimizationResult_dlib.x(xVariableIndex);
 #ifdef POLYSCOPE_DEFINED
         std::cout << results.optimalXNameValuePairs[xVariableIndex].first << ":"
@@ -834,9 +860,9 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza
         global.simulationScenarioIndices.size());
     results.objectiveValue.perSimulationScenario_total.resize(
         global.simulationScenarioIndices.size());
-#ifdef POLYSCOPE_DEFINED
+    //#ifdef POLYSCOPE_DEFINED
-    global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed);
+    //    global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed);
-#endif
+    //#endif
 
     results.perScenario_fullPatternPotentialEnergy.resize(global.simulationScenarioIndices.size());
     for (int i = 0; i < global.simulationScenarioIndices.size(); i++) {
@@ -1000,98 +1026,180 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces(
 void ReducedModelOptimizer::runOptimization(const Settings &settings,
                                             ReducedPatternOptimization::Results &results)
 {
-    global.objectiveValueHistory.clear();
+#if POLYSCOPE_DEFINED
+    global.objectiveValueHistoryY.clear();
+    global.objectiveValueHistoryY.reserve(settings.numberOfFunctionCalls);
+    global.objectiveValueHistoryX.reserve(settings.numberOfFunctionCalls);
+    global.plotColors.reserve(settings.numberOfFunctionCalls);
+#endif
     double (*objF)(const dlib::matrix<double, 0, 1> &) = &objective;
 
-    //Geometry optimization of the reduced pattern
+    dlib::function_evaluation optimalResult;
-    constexpr int numberOfGeometryOptimizationParameters = 2;
+    constexpr bool shouldJointlyOptimizeGeometryAndMaterial = false;
-    dlib::matrix<double, 0, 1> xGeometryMin(numberOfGeometryOptimizationParameters);
+    const auto hexAngleParameterIt = std::find_if(settings.xRanges.begin(),
-    dlib::matrix<double, 0, 1> xGeometryMax(numberOfGeometryOptimizationParameters);
+                                                  settings.xRanges.end(),
-    const int geometryParametersOffset
+                                                  [](const xRange &x) {
-        = global.numberOfOptimizationParameters
+                                                      return x.label == "HexAngle";
-          - numberOfGeometryOptimizationParameters; //I assume the geometry parameters come at the end
+                                                  });
-    int paramIndex = 0;
+    const bool hasHexAngleParameter = hexAngleParameterIt != settings.xRanges.end();
-    for (int i = geometryParametersOffset; i < global.numberOfOptimizationParameters;
+    const auto hexSizeParameterIt = std::find_if(settings.xRanges.begin(),
-         i++, paramIndex++) {
+                                                 settings.xRanges.end(),
-        xGeometryMin(paramIndex) = settings.xRanges[i].min;
+                                                 [](const xRange &x) {
-        xGeometryMax(paramIndex) = settings.xRanges[i].max;
+                                                     return x.label == "HexSize";
-    }
+                                                 });
-
+    const bool hasHexSizeParameter = hexSizeParameterIt != settings.xRanges.end();
-    //Set reduced pattern update functions to be used during optimization
+    const bool hasBothGeometricalParameters = hasHexAngleParameter && hasHexSizeParameter;
-    function_updateReducedPattern =
+    assert(hasBothGeometricalParameters || (!hasHexAngleParameter && !hasHexSizeParameter));
-        [&](const dlib::matrix<double, 0, 1> &x,
+    const int numberOfGeometryOptimizationParameters = hasBothGeometricalParameters ? 2 : 0;
-            std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
-            function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
-            pReducedPatternSimulationMesh->reset();
-        };
-
-    dlib::function_evaluation result_dlib_geometry
-        = dlib::find_min_global(objF,
-                                xGeometryMin,
-                                xGeometryMax,
-                                dlib::max_function_calls(settings.numberOfFunctionCalls),
-                                std::chrono::hours(24 * 365 * 290),
-                                settings.solverAccuracy);
-    //Material optimization of the reduced pattern
     const int numberOfMaterialOptimizationParameters = global.numberOfOptimizationParameters
                                                        - numberOfGeometryOptimizationParameters;
-    std::cout << "opt size:" << result_dlib_geometry.x(0) << std::endl;
+    const bool hasGeometryAndMaterialParameters = hasBothGeometricalParameters
-    std::cout << "opt size:" << global.minX[0] << std::endl;
+                                                  && numberOfMaterialOptimizationParameters != 0;
-    dlib::function_evaluation result_dlib_material;
+    if (!shouldJointlyOptimizeGeometryAndMaterial && hasGeometryAndMaterialParameters) {
-    if (numberOfMaterialOptimizationParameters != 0) {
+        //Geometry optimization of the reduced pattern
-        dlib::matrix<double, 0, 1> xMaterialMin(numberOfMaterialOptimizationParameters);
+        dlib::matrix<double, 0, 1> xGeometryMin(numberOfGeometryOptimizationParameters);
-        dlib::matrix<double, 0, 1> xMaterialMax(numberOfMaterialOptimizationParameters);
+        dlib::matrix<double, 0, 1> xGeometryMax(numberOfGeometryOptimizationParameters);
-        for (int i = 0; i < numberOfMaterialOptimizationParameters; i++) {
+        //        const int hexAngleParameterIndex = std::distance(settings.xRanges.begin(),
-            xMaterialMin(i) = settings.xRanges[i].min;
+        //                                                         hexAngleParameterIt);
-            xMaterialMax(i) = settings.xRanges[i].max;
+        //        const int hexSizeParameterIndex = std::distance(settings.xRanges.begin(),
+        //                                                        hexSizeParameterIt);
+        xGeometryMin(0) = hexSizeParameterIt->min;
+        xGeometryMax(0) = hexSizeParameterIt->max;
+        xGeometryMin(1) = hexAngleParameterIt->min;
+        xGeometryMax(1) = hexAngleParameterIt->max;
+
+        //Set reduced pattern update functions to be used during optimization
+        function_updateReducedPattern =
+            [&](const dlib::matrix<double, 0, 1> &x,
+                std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+                function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
+                pReducedPatternSimulationMesh->reset();
+            };
+        const int optimizationFunctionCalls_geometry
+            = global.numberOfOptimizationParameters == numberOfGeometryOptimizationParameters
+                  ? settings.numberOfFunctionCalls
+                  : static_cast<int>(settings.numberOfFunctionCalls / 3.0);
+        dlib::function_evaluation result_dlib_geometry
+            = dlib::find_min_global(objF,
+                                    xGeometryMin,
+                                    xGeometryMax,
+                                    dlib::max_function_calls(optimizationFunctionCalls_geometry),
+                                    std::chrono::hours(24 * 365 * 290),
+                                    settings.solverAccuracy);
+        //Material optimization of the reduced pattern
+        //        std::cout << "opt size:" << result_dlib_geometry.x(0) << std::endl;
+        //        std::cout << "opt size:" << global.minX[0] << std::endl;
+        dlib::function_evaluation result_dlib_material;
+        if (numberOfMaterialOptimizationParameters != 0) {
+            dlib::matrix<double, 0, 1> xMaterialMin(numberOfMaterialOptimizationParameters);
+            dlib::matrix<double, 0, 1> xMaterialMax(numberOfMaterialOptimizationParameters);
+            for (int i = 0; i < numberOfMaterialOptimizationParameters; i++) {
+                xMaterialMin(i) = settings.xRanges[i].min;
+                xMaterialMax(i) = settings.xRanges[i].max;
+            }
+            function_updateReducedPattern =
+                [&](const dlib::matrix<double, 0, 1> &x,
+                    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+                    function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
+                    function_updateReducedPattern_geometry(
+                        result_dlib_geometry.x,
+                        pReducedPatternSimulationMesh); //FIXME: I shouldn't need to call this
+                    pReducedPatternSimulationMesh->reset();
+                };
+            result_dlib_material = dlib::find_min_global(objF,
+                                                         xMaterialMin,
+                                                         xMaterialMax,
+                                                         dlib::max_function_calls(static_cast<int>(
+                                                             2.0 * settings.numberOfFunctionCalls
+                                                             / 3.0)),
+                                                         std::chrono::hours(24 * 365 * 290),
+                                                         settings.solverAccuracy);
         }
+        //    constexpr bool useBOBYQA = false;
+        //    if (useBOBYQA) {
+        //        const size_t npt = 2 * global.numberOfOptimizationParameters;
+        //        //      ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
+        //        const double rhobeg = 0.1;
+        //        //  const double rhobeg = 10;
+        //        const double rhoend = rhobeg * 1e-6;
+        //        //  const size_t maxFun = 10 * (x.size() ^ 2);
+        //        const size_t bobyqa_maxFunctionCalls = 200;
+        //        dlib::find_min_bobyqa(objF,
+        //                              result_dlib.x,
+        //                              npt,
+        //                              xMaterialMin,
+        //                              xMaterialMax,
+        //                              rhobeg,
+        //                              rhoend,
+        //                              bobyqa_maxFunctionCalls);
+        //    }
+
+        optimalResult.x.set_size(global.numberOfOptimizationParameters);
+        std::copy(result_dlib_material.x.begin(),
+                  result_dlib_material.x.end(),
+                  optimalResult.x.begin());
+        //                  debug_x.begin());
+        std::copy(result_dlib_geometry.x.begin(),
+                  result_dlib_geometry.x.end(),
+                  optimalResult.x.begin() + numberOfMaterialOptimizationParameters);
+        std::cout << "opt x:";
+        for (const auto optx : optimalResult.x) {
+            std::cout << optx << " ";
+        }
+        std::cout << std::endl;
         function_updateReducedPattern =
             [&](const dlib::matrix<double, 0, 1> &x,
                 std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
                 function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
+                function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
                 pReducedPatternSimulationMesh->reset();
             };
-        result_dlib_material = dlib::find_min_global(objF,
+        //        const auto meshLabel = global.reducedPatternSimulationJobs[0]->pMesh->getLabel();
-                                                     xMaterialMin,
+        //        global.reducedPatternSimulationJobs[0]->pMesh->setLabel("pre");
-                                                     xMaterialMax,
+        //        global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::reducedInitial);
-                                                     dlib::max_function_calls(
+        //        polyscope::show();
-                                                         settings.numberOfFunctionCalls),
+        //        global.reducedPatternSimulationJobs[0]->pMesh->unregister();
-                                                     std::chrono::hours(24 * 365 * 290),
+        optimalResult.y = objective(optimalResult.x);
-                                                     settings.solverAccuracy);
+    } else {
-    }
+        dlib::matrix<double, 0, 1> xMin(global.numberOfOptimizationParameters);
-    //    constexpr bool useBOBYQA = false;
+        dlib::matrix<double, 0, 1> xMax(global.numberOfOptimizationParameters);
-    //    if (useBOBYQA) {
+        for (int i = 0; i < global.numberOfOptimizationParameters; i++) {
-    //        const size_t npt = 2 * global.numberOfOptimizationParameters;
+            xMin(i) = settings.xRanges[i].min;
-    //        //      ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
+            xMax(i) = settings.xRanges[i].max;
-    //        const double rhobeg = 0.1;
+        }
-    //        //  const double rhobeg = 10;
+        if (numberOfGeometryOptimizationParameters != 0
-    //        const double rhoend = rhobeg * 1e-6;
+            && numberOfMaterialOptimizationParameters != 0) {
-    //        //  const size_t maxFun = 10 * (x.size() ^ 2);
+            function_updateReducedPattern =
-    //        const size_t bobyqa_maxFunctionCalls = 200;
+                [&](const dlib::matrix<double, 0, 1> &x,
-    //        dlib::find_min_bobyqa(objF,
+                    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
-    //                              result_dlib.x,
+                    function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
-    //                              npt,
+                    function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
-    //                              xMaterialMin,
+                    pReducedPatternSimulationMesh->reset();
-    //                              xMaterialMax,
+                };
-    //                              rhobeg,
+        } else if (numberOfGeometryOptimizationParameters == 0) {
-    //                              rhoend,
+            function_updateReducedPattern =
-    //                              bobyqa_maxFunctionCalls);
+                [&](const dlib::matrix<double, 0, 1> &x,
-    //    }
+                    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+                    function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
+                    pReducedPatternSimulationMesh->reset();
+                };
 
-    dlib::function_evaluation optimalResult;
+        } else {
-    optimalResult.x.set_size(global.numberOfOptimizationParameters);
+            function_updateReducedPattern =
-    std::copy(result_dlib_material.x.begin(), result_dlib_material.x.end(), optimalResult.x.begin());
+                [&](const dlib::matrix<double, 0, 1> &x,
-    std::copy(result_dlib_geometry.x.begin(),
+                    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
-              result_dlib_geometry.x.end(),
+                    function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
-              optimalResult.x.begin() + geometryParametersOffset);
+                    pReducedPatternSimulationMesh->reset();
-    function_updateReducedPattern =
+                };
-        [&](const dlib::matrix<double, 0, 1> &x,
+        }
-            std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+        optimalResult = dlib::find_min_global(objF,
-            function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
+                                              xMin,
-            function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
+                                              xMax,
-            pReducedPatternSimulationMesh->reset();
+                                              dlib::max_function_calls(
-        };
+                                                  settings.numberOfFunctionCalls),
-    optimalResult.y = objective(optimalResult.x);
+                                              std::chrono::hours(24 * 365 * 290),
+                                              settings.solverAccuracy);
+    }
     getResults(optimalResult, settings, results);
 }
 

src/reducedmodeloptimizer.hpp

@@ -67,9 +67,10 @@ public:
 
     SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
 
-    void initializePatterns(PatternGeometry &fullPattern,
+    void initializePatterns(
-                            PatternGeometry &reducedPatterm,
+        PatternGeometry &fullPattern,
-                            const int &optimizationParameters);
+        PatternGeometry &reducedPatterm,
+        const std::vector<ReducedPatternOptimization::xRange> &optimizationParameters);
 
     static void runSimulation(const std::string &filename, std::vector<double> &x);
 
@@ -183,9 +184,24 @@ public:
     static std::function<void(const dlib::matrix<double, 0, 1> &x,
                               std::shared_ptr<SimulationMesh> &m)>
         function_updateReducedPattern_material;
-    static std::function<void(const dlib::matrix<double, 0, 1> &x,
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
-                              std::shared_ptr<SimulationMesh> &m)>
+                              const double &newE)>
         function_updateReducedPattern_material_E;
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                              const double &newA)>
+        function_updateReducedPattern_material_A;
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                              const double &newI)>
+        function_updateReducedPattern_material_I;
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                              const double &newI2)>
+        function_updateReducedPattern_material_I2;
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                              const double &newI3)>
+        function_updateReducedPattern_material_I3;
+    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                              const double &newJ)>
+        function_updateReducedPattern_material_J;
     static std::function<void(const dlib::matrix<double, 0, 1> &x,
                               std::shared_ptr<SimulationMesh> &m)>
         function_updateReducedPattern_geometry;
@@ -203,8 +219,9 @@ private:
             &maxForceMagnitudes);
     void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern);
     void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
-    static void initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh,
+    static void initializeOptimizationParameters(
-                                                 const int &optimizationParamters);
+        const std::shared_ptr<SimulationMesh> &mesh,
+        const std::vector<ReducedPatternOptimization::xRange> &optimizationParamters);
 
     static double objective(long n, const double *x);
     DRMSimulationModel simulator;
@@ -224,11 +241,30 @@ private:
         const std::vector<ReducedPatternOptimization::BaseSimulationScenario>
             &desiredBaseSimulationScenarioIndices);
     static double objective(const dlib::matrix<double, 0, 1> &x);
+    static void initializeUpdateReducedPatternFunctions();
 };
 inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
     ReducedModelOptimizer::function_updateReducedPattern;
 inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
     ReducedModelOptimizer::function_updateReducedPattern_material;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newE)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_E;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newA)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_A;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newI)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_I;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newI2)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_I2;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newI3)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_I3;
+inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
+                          const double &newJ)>
+    ReducedModelOptimizer::function_updateReducedPattern_material_J;
 inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
     ReducedModelOptimizer::function_updateReducedPattern_geometry;