Objective value takes into account the rotational displacemenets between the full and the reduced pattern.

src/main.cpp

@@ -56,7 +56,7 @@ int main(int argc, char *argv[]) {
   settings_optimization.normalizationStrategy
       = ReducedModelOptimization::Settings::NormalizationStrategy::Epsilon;
   settings_optimization.normalizationParameter = 0.0003;
-  settings_optimization.solutionAccuracy = 0.01;
+  settings_optimization.solutionAccuracy = 0.001;
 
   // Optimize pair
   const std::string pairName =
@@ -97,7 +97,7 @@ int main(int argc, char *argv[]) {
   //                            .string(),
   //                        false);
   csv_results << "Name";
-  optimizationResults.writeHeaderTo(settings_optimization, csv_results);
+  optimizationResults.writeHeaderTo(csv_results);
   settings_optimization.writeHeaderTo(csv_results);
   csv_results << endrow;
   csv_results << pairName;

src/reducedmodeloptimizer.cpp

@@ -13,7 +13,8 @@ struct GlobalOptimizationVariables {
     std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
     //    std::vector<std::vector<Vector6d>> fullPatternDisplacements;
     std::vector<SimulationResults> fullPatternResults;
-    std::vector<double> objectiveNormalizationValues;
+    std::vector<double> translationalDisplacementNormalizationValues;
+    std::vector<double> rotationalDisplacementNormalizationValues;
     std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
     std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
     std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
@@ -67,12 +68,6 @@ double ReducedModelOptimizer::computeRawDisplacementError(
             reducedPatternDisplacements[reducedModelVi][0],
             reducedPatternDisplacements[reducedModelVi][1],
             reducedPatternDisplacements[reducedModelVi][2]);
-        if (!std::isfinite(reducedVertexDisplacement[0]) ||
-            !std::isfinite(reducedVertexDisplacement[1]) ||
-            !std::isfinite(reducedVertexDisplacement[2])) {
-            std::cout << "Displacements are not finite" << std::endl;
-            std::terminate();
-        }
         const VertexIndex fullModelVi = reducedFullViPair.second;
         Eigen::Vector3d fullVertexDisplacement(
             fullPatternDisplacements[fullModelVi][0],
@@ -134,11 +129,24 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
                     reducedModelResults.internalPotentialEnergy
                     - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy);
             } else {
-                simulationScenarioError = computeDisplacementError(
+                const double translationalError = computeDisplacementError(
                     reducedModelResults.displacements,
                     global.fullPatternResults[simulationScenarioIndex].displacements,
                     global.reducedToFullInterfaceViMap,
-                    global.objectiveNormalizationValues[simulationScenarioIndex]);
+                    global.translationalDisplacementNormalizationValues[simulationScenarioIndex]);
+                double rotationalError = 0;
+                for (const auto &reducedToFullInterfaceViPair : global.reducedToFullInterfaceViMap) {
+                    const double vertexRotationalError
+                        = global.fullPatternResults[simulationScenarioIndex]
+                              .rotationalDisplacementQuaternion[reducedToFullInterfaceViPair.second]
+                              .angularDistance(reducedModelResults.rotationalDisplacementQuaternion
+                                                   [reducedToFullInterfaceViPair.first]);
+                    rotationalError += vertexRotationalError;
+                }
+                simulationScenarioError += translationalError
+                                           + rotationalError
+                                                 / global.rotationalDisplacementNormalizationValues
+                                                       [simulationScenarioIndex];
             }
             //#ifdef POLYSCOPE_DEFINED
             //            std::cout << "sim error:" << simulationScenarioError << std::endl;
@@ -549,10 +557,9 @@ void ReducedModelOptimizer::visualizeResults(
     std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
         fullPatternSimulationJobs[0]->pMesh;
     pFullPatternSimulationMesh->registerForDrawing();
-    pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() +
+    pFullPatternSimulationMesh->save(pFullPatternSimulationMesh->getLabel() + "_undeformed.ply");
-                                        "_undeformed.ply");
+    reducedPatternSimulationJobs[0]->pMesh->save(reducedPatternSimulationJobs[0]->pMesh->getLabel()
-    reducedPatternSimulationJobs[0]->pMesh->savePly(
+                                                 + "_undeformed.ply");
-        reducedPatternSimulationJobs[0]->pMesh->getLabel() + "_undeformed.ply");
     double totalError = 0;
     for (const int simulationScenarioIndex : simulationScenarios) {
         const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
@@ -570,8 +577,8 @@ void ReducedModelOptimizer::visualizeResults(
         double normalizationFactor = 1;
         if (global.optimizationSettings.normalizationStrategy !=
             ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) {
-            normalizationFactor =
+            normalizationFactor
-                global.objectiveNormalizationValues[simulationScenarioIndex];
+                = global.translationalDisplacementNormalizationValues[simulationScenarioIndex];
         }
         reducedModelResults.saveDeformedModel();
         fullModelResults.saveDeformedModel();
@@ -680,10 +687,10 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
         std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
     ReducedModelOptimization::Results results;
     results.numberOfSimulationCrashes = global.numOfSimulationCrashes;
-    results.optimalXNameValuePairs.reserve(settings.xRanges.size());
+    results.optimalXNameValuePairs.resize(settings.xRanges.size());
     for(int xVariableIndex=0;xVariableIndex<settings.xRanges.size();xVariableIndex++){
-        results.optimalXNameValuePairs[settings.xRanges[xVariableIndex].label]=
+        results.optimalXNameValuePairs[xVariableIndex]
-         global.minX[xVariableIndex];
+            = std::make_pair(settings.xRanges[xVariableIndex].label, global.minX[xVariableIndex]);
     }
     results.objectiveValue = global.minY;
 #ifdef POLYSCOPE_DEFINED
@@ -713,7 +720,7 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
             reducedModelResults.displacements,
             global.fullPatternResults[simulationScenarioIndex].displacements,
             global.reducedToFullInterfaceViMap,
-            global.objectiveNormalizationValues[simulationScenarioIndex]);
+            global.translationalDisplacementNormalizationValues[simulationScenarioIndex]);
         results.rawObjectiveValue += computeRawDisplacementError(
             reducedModelResults.displacements,
             global.fullPatternResults[simulationScenarioIndex].displacements,
@@ -877,27 +884,61 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
         Settings::NormalizationStrategy::Epsilon) {
 
         // Compute the sum of the displacement norms
-        std::vector<double> fullPatternDisplacementNormSum(
+        std::vector<double> fullPatternTranslationalDisplacementNormSum(NumberOfSimulationScenarios);
-            NumberOfSimulationScenarios);
+        std::vector<double> fullPatternAngularDistance(NumberOfSimulationScenarios);
         for (int simulationScenarioIndex : global.simulationScenarioIndices) {
-            double displacementNormSum = 0;
+            double translationalDisplacementNormSum = 0;
+            double angularDistanceSum = 0;
             for (auto interfaceViPair : global.reducedToFullInterfaceViMap) {
+                const int fullPatternVi = interfaceViPair.second;
+                //If the full pattern vertex is translationally constrained dont take it into account
+                if (global.fullPatternSimulationJobs[simulationScenarioIndex]
+                        ->constrainedVertices.contains(fullPatternVi)) {
+                    const std::unordered_set<int> constrainedDof
+                        = global.fullPatternSimulationJobs[simulationScenarioIndex]
+                              ->constrainedVertices.at(fullPatternVi);
+                    if (constrainedDof.contains(0) && constrainedDof.contains(1)
+                        && constrainedDof.contains(2)) {
+                        continue;
+                    }
+                }
+
                 const Vector6d &vertexDisplacement = global
                                                          .fullPatternResults[simulationScenarioIndex]
-                                                         .displacements[interfaceViPair.second];
+                                                         .displacements[fullPatternVi];
-                displacementNormSum += vertexDisplacement.getTranslation().norm();
+                translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm();
+                //If the full pattern vertex is rotationally constrained dont take it into account
+                if (global.fullPatternSimulationJobs[simulationScenarioIndex]
+                        ->constrainedVertices.contains(fullPatternVi)) {
+                    const std::unordered_set<int> constrainedDof
+                        = global.fullPatternSimulationJobs[simulationScenarioIndex]
+                              ->constrainedVertices.at(fullPatternVi);
+                    if (constrainedDof.contains(3) && constrainedDof.contains(5)
+                        && constrainedDof.contains(4)) {
+                        continue;
+                    }
+                }
+                angularDistanceSum += global.fullPatternResults[simulationScenarioIndex]
+                                          .rotationalDisplacementQuaternion[fullPatternVi]
+                                          .angularDistance(Eigen::Quaterniond::Identity());
             }
-            fullPatternDisplacementNormSum[simulationScenarioIndex] =
+            fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex]
-                displacementNormSum;
+                = translationalDisplacementNormSum;
+            fullPatternAngularDistance[simulationScenarioIndex] = angularDistanceSum;
         }
 
         for (int simulationScenarioIndex : global.simulationScenarioIndices) {
             if (global.optimizationSettings.normalizationStrategy ==
                 Settings::NormalizationStrategy::Epsilon) {
-                const double epsilon =
+                const double epsilon_translationalDisplacement = global.optimizationSettings
-                    global.optimizationSettings.normalizationParameter;
+                                                                     .normalizationParameter;
-                global.objectiveNormalizationValues[simulationScenarioIndex] = std::max(
+                global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
-                    fullPatternDisplacementNormSum[simulationScenarioIndex], epsilon);
+                    = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex],
+                               epsilon_translationalDisplacement);
+                const double epsilon_rotationalDisplacement = vcg::math::ToRad(2.0);
+                global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
+                    = std::max(fullPatternAngularDistance[simulationScenarioIndex],
+                               epsilon_rotationalDisplacement);
             }
         }
     }
@@ -919,7 +960,8 @@ Results ReducedModelOptimizer::optimize(
     global.g_optimalReducedModelDisplacements.resize(NumberOfSimulationScenarios);
     global.reducedPatternSimulationJobs.resize(NumberOfSimulationScenarios);
     global.fullPatternResults.resize(NumberOfSimulationScenarios);
-    global.objectiveNormalizationValues.resize(NumberOfSimulationScenarios);
+    global.translationalDisplacementNormalizationValues.resize(NumberOfSimulationScenarios);
+    global.rotationalDisplacementNormalizationValues.resize(NumberOfSimulationScenarios);
     global.minY = std::numeric_limits<double>::max();
     global.numOfSimulationCrashes = 0;
     global.numberOfFunctionCalls = 0;