Added potential energy error term in addition to the displacement error term

src/main.cpp

@@ -153,60 +153,63 @@ int main(int argc, char *argv[])
         csv_resultsLocalFile << endrow;
     }
 
+    ReducedModelEvaluator::evaluateReducedModel(optimizationResults);
 #ifdef POLYSCOPE_DEFINED
-  std::vector<std::string> scenarioLabels(
-      optimizationResults.objectiveValue.perSimulationScenario_total.size());
-  const double colorAxial = 1;
-  const double colorShear = 3;
-  const double colorBending = 5;
-  const double colorDome = 0.1;
-  const double colorSaddle = 0;
-  std::vector<double> colors(optimizationResults.objectiveValue.perSimulationScenario_total.size());
-  for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) {
-      scenarioLabels[scenarioIndex]
-          = optimizationResults.reducedPatternSimulationJobs[scenarioIndex]->getLabel();
-      if (scenarioLabels[scenarioIndex].rfind("Axial", 0) == 0) {
-          colors[scenarioIndex] = colorAxial;
 
-      } else if (scenarioLabels[scenarioIndex].rfind("Shear", 0) == 0) {
+    std::vector<std::string> scenarioLabels(
-          colors[scenarioIndex] = colorShear;
+        optimizationResults.objectiveValue.perSimulationScenario_total.size());
+    const double colorAxial = 1;
+    const double colorShear = 3;
+    const double colorBending = 5;
+    const double colorDome = 0.1;
+    const double colorSaddle = 0;
+    std::vector<double> colors(
+        optimizationResults.objectiveValue.perSimulationScenario_total.size());
+    for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) {
+        scenarioLabels[scenarioIndex]
+            = optimizationResults.reducedPatternSimulationJobs[scenarioIndex]->getLabel();
+        if (scenarioLabels[scenarioIndex].rfind("Axial", 0) == 0) {
+            colors[scenarioIndex] = colorAxial;
 
-      } else if (scenarioLabels[scenarioIndex].rfind("Bending", 0) == 0) {
+        } else if (scenarioLabels[scenarioIndex].rfind("Shear", 0) == 0) {
-          colors[scenarioIndex] = colorBending;
+            colors[scenarioIndex] = colorShear;
 
-      } else if (scenarioLabels[scenarioIndex].rfind("Dome", 0) == 0) {
+        } else if (scenarioLabels[scenarioIndex].rfind("Bending", 0) == 0) {
-          colors[scenarioIndex] = colorDome;
+            colors[scenarioIndex] = colorBending;
-      } else if (scenarioLabels[scenarioIndex].rfind("Saddle", 0) == 0) {
-          colors[scenarioIndex] = colorSaddle;
-      } else {
-          std::cerr << "Label could not be identified" << std::endl;
-      }
-  }
-  std::vector<double> y(optimizationResults.objectiveValue.perSimulationScenario_total.size());
-  for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) {
-      y[scenarioIndex]
-          //          = optimizationResults.objectiveValue.perSimulationScenario_rawTranslational[scenarioIndex]
-          //            + optimizationResults.objectiveValue.perSimulationScenario_rawRotational[scenarioIndex];
-          = optimizationResults.objectiveValue.perSimulationScenario_total_unweighted[scenarioIndex];
-  }
-  std::vector<double> x = matplot::linspace(0, y.size() - 1, y.size());
-  std::vector<double> markerSizes(y.size(), 5);
-  SimulationResultsReporter::createPlot("scenario index",
-                                        "error",
-                                        x,
-                                        y,
-                                        markerSizes,
-                                        colors,
-                                        std::filesystem::path(resultsOutputDir)
-                                            .append("perScenarioObjectiveValues.png"));
-  //    optimizationResults.saveMeshFiles();
-  std::cout << "Saved results to:" << resultsOutputDir << std::endl;
-  //  optimizationResults.draw();
 
-  //    ReducedModelOptimization::Results optResults;
+        } else if (scenarioLabels[scenarioIndex].rfind("Dome", 0) == 0) {
-  //    optResults.load("/home/iason/Desktop/dlib_ensmallen_comparison/TestSets/"
+            colors[scenarioIndex] = colorDome;
-  //                    "singlePattern_dlib_firstSubmission/12@single_reduced(100000_1.20)");
+        } else if (scenarioLabels[scenarioIndex].rfind("Saddle", 0) == 0) {
-  ReducedModelEvaluator::evaluateReducedModel(optimizationResults);
+            colors[scenarioIndex] = colorSaddle;
+        } else {
+            std::cerr << "Label could not be identified" << std::endl;
+        }
+    }
+    std::vector<double> y(optimizationResults.objectiveValue.perSimulationScenario_total.size());
+    for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) {
+        y[scenarioIndex]
+            //          = optimizationResults.objectiveValue.perSimulationScenario_rawTranslational[scenarioIndex]
+            //            + optimizationResults.objectiveValue.perSimulationScenario_rawRotational[scenarioIndex];
+            = optimizationResults.objectiveValue
+                  .perSimulationScenario_total_unweighted[scenarioIndex];
+    }
+    std::vector<double> x = matplot::linspace(0, y.size() - 1, y.size());
+    std::vector<double> markerSizes(y.size(), 5);
+    SimulationResultsReporter::createPlot("scenario index",
+                                          "error",
+                                          x,
+                                          y,
+                                          markerSizes,
+                                          colors,
+                                          std::filesystem::path(resultsOutputDir)
+                                              .append("perScenarioObjectiveValues.png"));
+    //    optimizationResults.saveMeshFiles();
+    std::cout << "Saved results to:" << resultsOutputDir << std::endl;
+    //  optimizationResults.draw();
+
+    //    ReducedModelOptimization::Results optResults;
+    //    optResults.load("/home/iason/Desktop/dlib_ensmallen_comparison/TestSets/"
+    //                    "singlePattern_dlib_firstSubmission/12@single_reduced(100000_1.20)");
 #endif
   return 0;
 }

src/reducedmodeloptimizer.cpp

@@ -108,9 +108,9 @@ double ReducedModelOptimizer::computeRawRotationalError(
 {
     double rawRotationalError = 0;
     for (const auto &reducedToFullInterfaceViPair : reducedToFullInterfaceViMap) {
-        const double vertexRotationalError
+        const double vertexRotationalError = std::abs(
-            = rotatedQuaternion_fullPattern[reducedToFullInterfaceViPair.second].angularDistance(
+            rotatedQuaternion_fullPattern[reducedToFullInterfaceViPair.second].angularDistance(
-                rotatedQuaternion_reducedPattern[reducedToFullInterfaceViPair.first]);
+                rotatedQuaternion_reducedPattern[reducedToFullInterfaceViPair.first]));
         rawRotationalError += vertexRotationalError;
     }
 
@@ -247,7 +247,7 @@ double ReducedModelOptimizer::objective(const std::vector<double> &x)
     //    FormFinder simulator;
 
     std::for_each(
-        std::execution::par_unseq,
+        //        std::execution::par_unseq,
         global.simulationScenarioIndices.begin(),
         global.simulationScenarioIndices.end(),
         [&](const int &simulationScenarioIndex) {
@@ -267,15 +267,11 @@ double ReducedModelOptimizer::objective(const std::vector<double> &x)
                 std::cout << "Failed simulation" << std::endl;
 #endif
             } else {
-//                double simulationScenarioError = 0;
+                const double simulationScenarioError_PE = std::abs(
-//            constexpr bool usePotentialEnergy = false;
+                    (reducedModelResults.internalPotentialEnergy
-//            if (usePotentialEnergy) {
+                     - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy)
-//                simulationScenarioError += std::abs(
+                    / global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy);
-//                    reducedModelResults.internalPotentialEnergy
+//             else {
-//                    - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy
-//                          / global.fullPatternResults[simulationScenarioIndex]
-//                                .internalPotentialEnergy);
-//            } else {
 #ifdef POLYSCOPE_DEFINED
 #ifdef USE_SCENARIO_WEIGHTS
                 const double scenarioWeight = global.scenarioWeights[simulationScenarioIndex];
@@ -285,7 +281,7 @@ double ReducedModelOptimizer::objective(const std::vector<double> &x)
 #else
                 const double scenarioWeight = 1;
 #endif
-                const double simulationScenarioError = computeError(
+                const double simulationScenarioError_displacements = computeError(
                     global.fullPatternResults[simulationScenarioIndex],
                     reducedModelResults,
                     global.reducedToFullInterfaceViMap,
@@ -294,7 +290,8 @@ double ReducedModelOptimizer::objective(const std::vector<double> &x)
                     scenarioWeight,
                     global.objectiveWeights[simulationScenarioIndex]);
 
-                simulationErrorsPerScenario[simulationScenarioIndex] = simulationScenarioError;
+                simulationErrorsPerScenario[simulationScenarioIndex]
+                    = simulationScenarioError_displacements + simulationScenarioError_PE;
                 //            }
                 //                #ifdef POLYSCOPE_DEFINED
                 //                            reducedJob->pMesh->registerForDrawing(Colors::reducedInitial);
@@ -978,7 +975,8 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces(
             .append(m_pFullPatternSimulationMesh->getLabel() + ".json"));
     const bool fullPatternScenarioMagnitudesExist = std::filesystem::exists(
         patternMaxForceMagnitudesFilePath);
-    if (fullPatternScenarioMagnitudesExist) {
+    constexpr bool recomputeMagnitudes = false;
+    if (fullPatternScenarioMagnitudesExist && !recomputeMagnitudes) {
         nlohmann::json json;
         std::ifstream ifs(patternMaxForceMagnitudesFilePath.string());
         ifs >> json;
@@ -1798,9 +1796,9 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors()
                     continue;
                 }
             }
-            angularDistanceSum += global.fullPatternResults[simulationScenarioIndex]
+            angularDistanceSum += std::abs(global.fullPatternResults[simulationScenarioIndex]
-                                      .rotationalDisplacementQuaternion[fullPatternVi]
+                                               .rotationalDisplacementQuaternion[fullPatternVi]
-                                      .angularDistance(Eigen::Quaterniond::Identity());
+                                               .angularDistance(Eigen::Quaterniond::Identity()));
         }
 
         fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex]
@@ -1916,7 +1914,7 @@ void ReducedModelOptimizer::optimize(
                 .append(m_pFullPatternSimulationMesh->getLabel())
                 .append(pFullPatternSimulationJob->getLabel()));
         //                .append(pFullPatternSimulationJob->getLabel() + ".json")
-        constexpr bool recomputeFullPatternResults = false;
+        constexpr bool recomputeFullPatternResults = true;
         SimulationResults fullPatternResults;
         if (!recomputeFullPatternResults && std::filesystem::exists(jobResultsDirectoryPath)) {
             fullPatternResults.load(std::filesystem::path(jobResultsDirectoryPath).append("Results"),

src/reducedmodeloptimizer.hpp

@@ -66,7 +66,7 @@ public:
     //    inline constexpr static ParameterLabels parameterLabels();
 
     inline static std::array<std::string, ReducedModelOptimization::NumberOfOptimizationVariables>
-        parameterLabels = {"R", "A", "I2", "I3", "J", "Theta", "R"};
+        parameterLabels = {"E", "A", "I2", "I3", "J", "Theta", "R"};
     constexpr static std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
         simulationScenariosResolution = {11, 11, 20, 20, 20};
     constexpr static std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>