Refactoring

src/main.cpp

@@ -65,7 +65,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, reducedPattern, {},settings_optimization.xRanges.size());
+  optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges.size());
   ReducedModelOptimization::Results optimizationResults =
       optimizer.optimize(settings_optimization);
 
@@ -106,7 +106,7 @@ int main(int argc, char *argv[]) {
   csv_results << endrow;
 
 #ifdef POLYSCOPE_DEFINED
-    optimizationResults.draw();
+  optimizationResults.draw();
 #endif
   return 0;
 }

src/reducedmodeloptimizer.cpp

@@ -39,14 +39,14 @@ struct GlobalOptimizationVariables {
 std::vector<SimulationJob> reducedPatternMaximumDisplacementSimulationJobs;
 
 double ReducedModelOptimizer::computeDisplacementError(
-    const std::vector<Vector6d> &reducedPatternDisplacements,
     const std::vector<Vector6d> &fullPatternDisplacements,
+    const std::vector<Vector6d> &reducedPatternDisplacements,
     const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
         &reducedToFullInterfaceViMap,
     const double &normalizationFactor)
 {
-    const double rawError = computeRawDisplacementError(reducedPatternDisplacements,
+    const double rawError = computeRawDisplacementError(fullPatternDisplacements,
-                                                        fullPatternDisplacements,
+                                                        reducedPatternDisplacements,
                                                         reducedToFullInterfaceViMap);
     if (global.optimizationSettings.normalizationStrategy !=
         Settings::NormalizationStrategy::NonNormalized) {
@@ -56,8 +56,8 @@ double ReducedModelOptimizer::computeDisplacementError(
 }
 
 double ReducedModelOptimizer::computeRawDisplacementError(
-    const std::vector<Vector6d> &reducedPatternDisplacements,
     const std::vector<Vector6d> &fullPatternDisplacements,
+    const std::vector<Vector6d> &reducedPatternDisplacements,
     const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
         &reducedToFullInterfaceViMap)
 {
@@ -81,6 +81,58 @@ double ReducedModelOptimizer::computeRawDisplacementError(
 
     return error;
 }
+
+double ReducedModelOptimizer::computeRawRotationalError(
+    const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
+    const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
+    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+        &reducedToFullInterfaceViMap)
+{
+    double rawRotationalError = 0;
+    for (const auto &reducedToFullInterfaceViPair : reducedToFullInterfaceViMap) {
+        const double vertexRotationalError
+            = rotatedQuaternion_fullPattern[reducedToFullInterfaceViPair.second].angularDistance(
+                rotatedQuaternion_reducedPattern[reducedToFullInterfaceViPair.first]);
+        rawRotationalError += vertexRotationalError;
+    }
+
+    return rawRotationalError;
+}
+
+double ReducedModelOptimizer::computeRotationalError(
+    const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
+    const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
+    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+        &reducedToFullInterfaceViMap,
+    const double &normalizationFactor)
+{
+    const double rawRotationalError = computeRawRotationalError(rotatedQuaternion_fullPattern,
+                                                                rotatedQuaternion_reducedPattern,
+                                                                reducedToFullInterfaceViMap);
+    return rawRotationalError / normalizationFactor;
+}
+
+double ReducedModelOptimizer::computeError(
+    const SimulationResults &simulationResults_fullPattern,
+    const SimulationResults &simulationResults_reducedPattern,
+    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+        &reducedToFullInterfaceViMap,
+    const double &normalizationFactor_translationalDisplacement,
+    const double &normalizationFactor_rotationalDisplacement)
+{
+    const double translationalError
+        = computeDisplacementError(simulationResults_fullPattern.displacements,
+                                   simulationResults_reducedPattern.displacements,
+                                   reducedToFullInterfaceViMap,
+                                   normalizationFactor_translationalDisplacement);
+    const double rotationalError
+        = computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion,
+                                 simulationResults_reducedPattern.rotationalDisplacementQuaternion,
+                                 reducedToFullInterfaceViMap,
+                                 normalizationFactor_rotationalDisplacement);
+    return translationalError + rotationalError;
+}
+
 double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3,
                                         double innerHexagonSize,
                                         double innerHexagonRotationAngle) {
@@ -129,24 +181,12 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
                     reducedModelResults.internalPotentialEnergy
                     - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy);
             } else {
-                const double translationalError = computeDisplacementError(
+                simulationScenarioError += computeError(
-                    reducedModelResults.displacements,
+                    global.fullPatternResults[simulationScenarioIndex],
-                    global.fullPatternResults[simulationScenarioIndex].displacements,
+                    reducedModelResults,
                     global.reducedToFullInterfaceViMap,
-                    global.translationalDisplacementNormalizationValues[simulationScenarioIndex]);
+                    global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
-                double rotationalError = 0;
+                    global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
-                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;
@@ -246,7 +286,6 @@ void ReducedModelOptimizer::createSimulationMeshes(
 }
 
 void ReducedModelOptimizer::computeMaps(
-    const std::unordered_set<size_t> &reducedModelExcludedEdges,
     const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
     PatternGeometry &fullPattern, PatternGeometry &reducedPattern,
     std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
@@ -375,13 +414,15 @@ void ReducedModelOptimizer::computeMaps(
     }
 }
 
-void ReducedModelOptimizer::computeMaps(
+void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern,
-    PatternGeometry &fullPattern, PatternGeometry &reducedPattern,
+                                        PatternGeometry &reducedPattern)
-    const std::unordered_set<size_t> &reducedModelExcludedEdges) {
+{
-    ReducedModelOptimizer::computeMaps(
+    ReducedModelOptimizer::computeMaps(slotToNode,
-        reducedModelExcludedEdges, slotToNode, fullPattern, reducedPattern,
+                                       fullPattern,
-        global.reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap,
+                                       reducedPattern,
-        m_fullPatternOppositeInterfaceViMap);
+                                       global.reducedToFullInterfaceViMap,
+                                       m_fullToReducedInterfaceViMap,
+                                       m_fullPatternOppositeInterfaceViMap);
 }
 
 ReducedModelOptimizer::ReducedModelOptimizer(
@@ -390,11 +431,10 @@ ReducedModelOptimizer::ReducedModelOptimizer(
     FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode);
 }
 
-void ReducedModelOptimizer::initializePatterns(
+void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
-    PatternGeometry &fullPattern, PatternGeometry &reducedPattern,
+                                               PatternGeometry &reducedPattern,
-    const std::unordered_set<size_t> &reducedModelExcludedEdges,const int& optimizationParameters) {
+                                               const int &optimizationParameters)
-    // fullPattern.setLabel("full_pattern_" + fullPattern.getLabel());
+{
-    // reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel());
     assert(fullPattern.VN() == reducedPattern.VN() &&
            fullPattern.EN() >= reducedPattern.EN());
 #if POLYSCOPE_DEFINED
@@ -406,31 +446,8 @@ void ReducedModelOptimizer::initializePatterns(
     copyFullPattern.copy(fullPattern);
     copyReducedPattern.copy(reducedPattern);
     global.baseTriangle = copyReducedPattern.getBaseTriangle();
-    /*
+
-   * Here we create the vector that connects the central node with the bottom
+    computeMaps(copyFullPattern, copyReducedPattern);
-   * left node of the base triangle. During the optimization the vi%2==0 nodes
-   * move on these vectors.
-   * */
-    //    const double h = copyReducedPattern.getBaseTriangleHeight();
-    //    double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0));
-    //    VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h,
-    //                                           0);
-    //    for (int rotationCounter = 0; rotationCounter < fanSize; rotationCounter++) {
-    //        VectorType rotatedVector =
-    //            vcg::RotationMatrix(patternPlaneNormal,
-    //                                vcg::math::ToRad(rotationCounter * 60.0)) *
-    //            baseTriangle_leftBottomNode;
-    //        global.g_innerHexagonVectors[rotationCounter] = rotatedVector;
-    //    }
-    //    const double innerHexagonInitialPos_x =
-    //        copyReducedPattern.vert[0].cP()[0] / global.g_innerHexagonVectors[0][0];
-    //    const double innerHexagonInitialPos_y =
-    //        copyReducedPattern.vert[0].cP()[1] / global.g_innerHexagonVectors[0][1];
-    //    global.innerHexagonInitialPos = innerHexagonInitialPos_x;
-    //    global.innerHexagonInitialRotationAngle =
-    //        30; /* NOTE: Here I assume that the CW reduced pattern is given as input.
-    //             This is not very generic */
-    computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges);
     createSimulationMeshes(copyFullPattern, copyReducedPattern);
     initializeOptimizationParameters(m_pReducedPatternSimulationMesh,optimizationParameters);
 }
@@ -716,17 +733,31 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
     for (int simulationScenarioIndex:global.simulationScenarioIndices) {
         SimulationResults reducedModelResults = simulator.executeSimulation(
             global.reducedPatternSimulationJobs[simulationScenarioIndex]);
-        const double error = computeDisplacementError(
+
-            reducedModelResults.displacements,
+        const double rawTranslationalError = computeRawDisplacementError(
             global.fullPatternResults[simulationScenarioIndex].displacements,
-            global.reducedToFullInterfaceViMap,
-            global.translationalDisplacementNormalizationValues[simulationScenarioIndex]);
-        results.rawObjectiveValue += computeRawDisplacementError(
             reducedModelResults.displacements,
-            global.fullPatternResults[simulationScenarioIndex].displacements,
             global.reducedToFullInterfaceViMap);
-        results.objectiveValuePerSimulationScenario[simulationScenarioIndex] =
+        const double rawRotationalError = computeRawRotationalError(
-            error;
+            global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion,
+            reducedModelResults.rotationalDisplacementQuaternion,
+            global.reducedToFullInterfaceViMap);
+
+        results.rawObjectiveValue += rawTranslationalError + rawRotationalError;
+
+        const double normalizedTranslationalError
+            = rawTranslationalError
+              / global.translationalDisplacementNormalizationValues[simulationScenarioIndex];
+        const double normalizedRotationalError
+            = rawRotationalError
+              / global.rotationalDisplacementNormalizationValues[simulationScenarioIndex];
+        std::cout << "Simulation scenario:" << simulationScenarioStrings[simulationScenarioIndex]
+                  << std::endl;
+        std::cout << "Translational error:" << normalizedTranslationalError << std::endl;
+        std::cout << "Rotational error:" << normalizedRotationalError << std::endl;
+        std::cout << std::endl;
+        results.objectiveValuePerSimulationScenario[simulationScenarioIndex]
+            = normalizedTranslationalError + normalizedRotationalError;
     }
     results.time = elapsed.count() / 1000.0;
     const bool printDebugInfo = false;
@@ -880,8 +911,6 @@ ReducedModelOptimizer::createScenarios(
 }
 
 void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
-    if (global.optimizationSettings.normalizationStrategy ==
-        Settings::NormalizationStrategy::Epsilon) {
 
         // Compute the sum of the displacement norms
         std::vector<double> fullPatternTranslationalDisplacementNormSum(NumberOfSimulationScenarios);
@@ -939,10 +968,12 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
                 global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
                     = std::max(fullPatternAngularDistance[simulationScenarioIndex],
                                epsilon_rotationalDisplacement);
+            } else {
+                global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = 1;
+                global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1;
             }
         }
     }
-}
 
 Results ReducedModelOptimizer::optimize(
     const Settings &optimizationSettings,

src/reducedmodeloptimizer.hpp

@@ -47,8 +47,9 @@ public:
   SimulationJob
   getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
 
-  void initializePatterns(PatternGeometry &fullPattern, PatternGeometry &reducedPatterm,
+  void initializePatterns(PatternGeometry &fullPattern,
-      const std::unordered_set<size_t> &reducedModelExcludedEges, const int &optimizationParameters);
+                          PatternGeometry &reducedPatterm,
+                          const int &optimizationParameters);
 
   static void runSimulation(const std::string &filename,
                             std::vector<double> &x);
@@ -66,16 +67,15 @@ public:
       PatternGeometry &fullModel, PatternGeometry &reducedModel,
       std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
       std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh);
-  static void computeMaps(
+  static void computeMaps(const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
-      const std::unordered_set<size_t> &reducedModelExcludedEdges,
+                          PatternGeometry &fullPattern,
-      const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
+                          PatternGeometry &reducedPattern,
-      PatternGeometry &fullPattern, PatternGeometry &reducedPattern,
+                          std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
-      std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+                              &reducedToFullInterfaceViMap,
-          &reducedToFullInterfaceViMap,
+                          std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
-      std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
+                              &fullToReducedInterfaceViMap,
-          &fullToReducedInterfaceViMap,
+                          std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
-      std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
+                              &fullPatternOppositeInterfaceViMap);
-          &fullPatternOppositeInterfaceViMap);
   static void visualizeResults(
       const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs,
       const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs,
@@ -90,38 +90,59 @@ public:
           &reducedToFullInterfaceViMap);
 
   static double computeRawDisplacementError(
-      const std::vector<Vector6d> &reducedPatternDisplacements,
       const std::vector<Vector6d> &fullPatternDisplacements,
+      const std::vector<Vector6d> &reducedPatternDisplacements,
       const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
           &reducedToFullInterfaceViMap);
 
   static double computeDisplacementError(
-      const std::vector<Vector6d> &reducedPatternDisplacements,
       const std::vector<Vector6d> &fullPatternDisplacements,
+      const std::vector<Vector6d> &reducedPatternDisplacements,
       const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
           &reducedToFullInterfaceViMap,
       const double &normalizationFactor);
 
-  static double objective(double E, double A, double innerHexagonSize, double innerHexagonRotationAngle);
+  static double objective(double E,
-  private:
+                          double A,
-      static void computeDesiredReducedModelDisplacements(
+                          double innerHexagonSize,
-          const SimulationResults &fullModelResults,
+                          double innerHexagonRotationAngle);
-          const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
+  static double computeRawRotationalError(
-          Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
+      const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
-      static ReducedModelOptimization::Results runOptimization(
+      const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
-          const ReducedModelOptimization::Settings &settings);
+      const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
-      std::vector<std::shared_ptr<SimulationJob>> createScenarios(
+          &reducedToFullInterfaceViMap);
-          const std::shared_ptr<SimulationMesh> &pMesh);
-      void computeMaps(PatternGeometry &fullModel,
-                       PatternGeometry &reducedPattern,
-                       const std::unordered_set<size_t> &reducedModelExcludedEges);
-      void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
-      static void initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh,
-                                                   const int &optimizationParamters);
 
-      static double objective(long n, const double *x);
+  static double computeRotationalError(
-      DRMSimulationModel simulator;
+      const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
-      void computeObjectiveValueNormalizationFactors();
+      const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
+      const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+          &reducedToFullInterfaceViMap,
+      const double &normalizationFactor);
+  static double computeError(
+      const SimulationResults &simulationResults_fullPattern,
+      const SimulationResults &simulationResults_reducedPattern,
+      const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+          &reducedToFullInterfaceViMap,
+      const double &normalizationFactor_translationalDisplacement,
+      const double &normalizationFactor_rotationalDisplacement);
+
+  private:
+  static void computeDesiredReducedModelDisplacements(
+      const SimulationResults &fullModelResults,
+      const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
+      Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
+  static ReducedModelOptimization::Results runOptimization(
+      const ReducedModelOptimization::Settings &settings);
+  std::vector<std::shared_ptr<SimulationJob>> createScenarios(
+      const std::shared_ptr<SimulationMesh> &pMesh);
+  void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern);
+  void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
+  static void initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh,
+                                               const int &optimizationParamters);
+
+  static double objective(long n, const double *x);
+  DRMSimulationModel simulator;
+  void computeObjectiveValueNormalizationFactors();
 };
 void updateMesh(long n, const double *x);
 #endif // REDUCEDMODELOPTIMIZER_HPP