The linear simulation model is used for optimizing the reduced model.

src/linearsimulationmodel.cpp

@@ -0,0 +1,21 @@
+#include "linearsimulationmodel.hpp"
+#include <Eigen/Core>
+#include <filesystem>
+//#include <igl/list_to_matrix.h>
+#include <iostream>
+#include <vcg/complex/algorithms/create/platonic.h>
+#include <vcg/complex/algorithms/update/normal.h>
+
+
+
+
+
+
+
+
+
+
+
+
+
+

src/linearsimulationmodel.hpp

@@ -0,0 +1,215 @@
+#ifndef LINEARSIMULATIONMODEL_HPP
+#define LINEARSIMULATIONMODEL_HPP
+
+//#include "beam.hpp"
+#include "simulationresult.hpp"
+#include "threed_beam_fea.h"
+#include <filesystem>
+#include <vector>
+
+// struct BeamSimulationProperties {
+//  float crossArea;
+//  float I2;
+//  float I3;
+//  float polarInertia;
+//  float G;
+//  // Properties used by fea
+//  float EA;
+//  float EIz;
+//  float EIy;
+//  float GJ;
+
+//  BeamSimulationProperties(const BeamDimensions &dimensions,
+//                           const BeamMaterial &material);
+//};
+
+// struct NodalForce {
+//  int index;
+//  int dof;
+//  double magnitude;
+//};
+
+// struct SimulationJob {
+//  Eigen::MatrixX3d nodes;
+//  Eigen::MatrixX2i elements;
+//  Eigen::MatrixX3d elementalNormals;
+//  Eigen::VectorXi fixedNodes;
+//  std::vector<NodalForce> nodalForces;
+//  std::vector<BeamDimensions> beamDimensions;
+//  std::vector<BeamMaterial> beamMaterial;
+//};
+
+// struct SimulationResults {
+//  std::vector<Eigen::VectorXd> edgeForces; ///< Force values per force
+//  component
+//                                           ///< #force components x #edges
+//  Eigen::MatrixXd
+//      nodalDisplacements; ///< The displacement of each node #nodes x 3
+//  SimulationResults(const fea::Summary &feaSummary);
+//  SimulationResults() {}
+//};
+
+class LinearSimulationModel {
+public:
+    LinearSimulationModel(){
+
+    }
+    static std::vector<fea::Elem>
+    getFeaElements(const std::shared_ptr<SimulationJob> &job) {
+        const int numberOfEdges = job->pMesh->EN();
+        std::vector<fea::Elem> elements(numberOfEdges);
+        for (int edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
+                const SimulationMesh::CoordType &evn0 =
+                        job->pMesh->edge[edgeIndex].cV(0)->cN();
+                const SimulationMesh::CoordType &evn1 =
+                        job->pMesh->edge[edgeIndex].cV(1)->cN();
+                const std::vector<double> nAverageVector{(evn0[0] + evn1[0]) / 2,
+                            (evn0[1] + evn1[1]) / 2,
+                            (evn0[2] + evn1[2]) / 2};
+                const Element &element = job->pMesh->elements[edgeIndex];
+                const double E = element.material.youngsModulus;
+                fea::Props feaProperties(E * element.A, E * element.I3, E * element.I2,
+                                         element.G * element.J, nAverageVector);
+                const int vi0 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(0));
+                const int vi1 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(1));
+                elements[edgeIndex] = fea::Elem(vi0, vi1, feaProperties);
+            }
+
+        return elements;
+    }
+    static std::vector<fea::Node>
+    getFeaNodes(const std::shared_ptr<SimulationJob> &job) {
+        const int numberOfNodes = job->pMesh->VN();
+        std::vector<fea::Node> feaNodes(numberOfNodes);
+        for (int vi = 0; vi < numberOfNodes; vi++) {
+                const CoordType &p = job->pMesh->vert[vi].cP();
+                feaNodes[vi] = fea::Node(p[0], p[1], p[2]);
+            }
+
+        return feaNodes;
+    }
+    static std::vector<fea::BC>
+    getFeaFixedNodes(const std::shared_ptr<SimulationJob> &job) {
+        std::vector<fea::BC> boundaryConditions;
+        boundaryConditions.reserve(job->constrainedVertices.size() * 6);
+        for (auto fixedVertex : job->constrainedVertices) {
+                const int vertexIndex = fixedVertex.first;
+                for (int dofIndex : fixedVertex.second) {
+                        boundaryConditions.emplace_back(
+                                    fea::BC(vertexIndex, static_cast<fea::DOF>(dofIndex), 0));
+                    }
+            }
+
+        return boundaryConditions;
+    }
+    static std::vector<fea::Force>
+    getFeaNodalForces(const std::shared_ptr<SimulationJob> &job) {
+        std::vector<fea::Force> nodalForces;
+        nodalForces.reserve(job->nodalExternalForces.size() * 6);
+
+        for (auto nodalForce : job->nodalExternalForces) {
+                for (int dofIndex = 0; dofIndex < 6; dofIndex++) {
+                        if (nodalForce.second[dofIndex] == 0) {
+                                continue;
+                            }
+                        nodalForces.emplace_back(
+                                    fea::Force(nodalForce.first, dofIndex, nodalForce.second[dofIndex]));
+                    }
+            }
+
+        return nodalForces;
+    }
+    static SimulationResults getResults(const fea::Summary &feaSummary) {
+        SimulationResults results;
+
+        results.executionTime = feaSummary.total_time_in_ms * 1000;
+        // displacements
+        results.displacements.resize(feaSummary.num_nodes);
+        for (int vi = 0; vi < feaSummary.num_nodes; vi++) {
+                results.displacements[vi] = Vector6d(feaSummary.nodal_displacements[vi]);
+            }
+
+        //  // Convert forces
+        //  // Convert to vector of eigen matrices of the form force component-> per
+        //  // Edge
+        //  const int numDof = 6;
+        //  const size_t numberOfEdges = feaSummary.element_forces.size();
+        //  edgeForces =
+        //      std::vector<Eigen::VectorXd>(numDof, Eigen::VectorXd(2 *
+        //      numberOfEdges));
+        //  for (gsl::index edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
+        //    for (gsl::index forceComponentIndex = 0; forceComponentIndex < numDof;
+        //         forceComponentIndex++) {
+        //      (edgeForces[forceComponentIndex])(2 * edgeIndex) =
+        //          feaSummary.element_forces[edgeIndex][forceComponentIndex];
+        //      (edgeForces[forceComponentIndex])(2 * edgeIndex + 1) =
+        //          feaSummary.element_forces[edgeIndex][numDof +
+        //          forceComponentIndex];
+        //    }
+        //  }
+        return results;
+    }
+
+  SimulationResults
+  executeSimulation(const std::shared_ptr<SimulationJob> &simulationJob) {
+      assert(simulationJob->pMesh->VN() != 0);
+      fea::Job job(getFeaNodes(simulationJob), getFeaElements(simulationJob));
+      //  printInfo(job);
+
+      // create the default options
+      fea::Options opts;
+      opts.save_elemental_forces = false;
+      opts.save_nodal_displacements = false;
+      opts.save_nodal_forces = false;
+      opts.save_report = false;
+      opts.save_tie_forces = false;
+      //  if (!elementalForcesOutputFilepath.empty()) {
+      //    opts.save_elemental_forces = true;
+      //    opts.elemental_forces_filename = elementalForcesOutputFilepath;
+      //  }
+      //  if (!nodalDisplacementsOutputFilepath.empty()) {
+      //    opts.save_nodal_displacements = true;
+      //    opts.nodal_displacements_filename = nodalDisplacementsOutputFilepath;
+      //  }
+
+      // have the program output status updates
+      opts.verbose = false;
+
+      // form an empty vector of ties
+      std::vector<fea::Tie> ties;
+
+      // also create an empty list of equations
+      std::vector<fea::Equation> equations;
+      auto fixedVertices = getFeaFixedNodes(simulationJob);
+      auto nodalForces = getFeaNodalForces(simulationJob);
+      fea::Summary feaResults =
+              fea::solve(job, fixedVertices, nodalForces, ties, equations, opts);
+
+      SimulationResults results = getResults(feaResults);
+      results.job = simulationJob;
+      return results;
+  }
+  //  SimulationResults getResults() const;
+  //  void setResultsNodalDisplacementCSVFilepath(const std::string
+  //  &outputPath); void setResultsElementalForcesCSVFilepath(const std::string
+  //  &outputPath);
+
+private:
+  //  std::string nodalDisplacementsOutputFilepath{"nodal_displacement.csv"};
+  //  std::string elementalForcesOutputFilepath{"elemental_forces.csv"};
+  //  SimulationResults results;
+
+  static void printInfo(const fea::Job &job) {
+      std::cout << "Details regarding the fea::Job:" << std::endl;
+      std::cout << "Nodes:" << std::endl;
+      for (fea::Node n : job.nodes) {
+              std::cout << n << std::endl;
+          }
+      std::cout << "Elements:" << std::endl;
+      for (Eigen::Vector2i e : job.elems) {
+              std::cout << e << std::endl;
+          }
+  }
+};
+
+#endif // LINEARSIMULATIONMODEL_HPP

src/main.cpp

@@ -1,7 +1,6 @@
 #include "beamformfinder.hpp"
 #include "csvfile.hpp"
 #include "edgemesh.hpp"
-#include "flatpattern.hpp"
 #include "reducedmodeloptimizer.hpp"
 #include "simulationhistoryplotter.hpp"
 #include "trianglepattterntopology.hpp"
@@ -29,13 +28,13 @@ int main(int argc, char *argv[]) {
   // Populate the pattern pair to be optimized
   ////Full pattern
   const std::string filepath_fullPattern = argv[1];
-  FlatPattern fullPattern(filepath_fullPattern);
+  PatternGeometry fullPattern(filepath_fullPattern);
   fullPattern.setLabel(
       std::filesystem::path(filepath_fullPattern).stem().string());
   fullPattern.scale(0.03);
   ////Reduced pattern
   const std::string filepath_reducedPattern = argv[2];
-  FlatPattern reducedPattern(filepath_reducedPattern);
+  PatternGeometry reducedPattern(filepath_reducedPattern);
   reducedPattern.setLabel(
       std::filesystem::path(filepath_reducedPattern).stem().string());
   reducedPattern.scale(0.03);
@@ -44,15 +43,16 @@ int main(int argc, char *argv[]) {
   ReducedModelOptimizer::xRange beamWidth{"B", 0.5, 1.5};
   ReducedModelOptimizer::xRange beamDimensionsRatio{"bOverh", 0.7, 1.3};
   ReducedModelOptimizer::xRange beamE{"E", 0.1, 1.9};
-  ReducedModelOptimizer::xRange innerHexagonSize{"HS", 0.1, 0.9};
+  ReducedModelOptimizer::xRange innerHexagonSize{"HexSize", 0.1, 0.8};
+  ReducedModelOptimizer::xRange innerHexagonAngle{"HexAngle", -29.5, 29.5};
   ReducedModelOptimizer::Settings settings_optimization;
   settings_optimization.xRanges = {beamWidth, beamDimensionsRatio, beamE,
-                                   innerHexagonSize};
+                                   innerHexagonSize, innerHexagonAngle};
   const bool input_numberOfFunctionCallsDefined = argc >= 4;
   settings_optimization.numberOfFunctionCalls =
       input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100;
   settings_optimization.normalizationStrategy =
-      ReducedModelOptimizer::Settings::NormalizationStrategy::NonNormalized;
+      ReducedModelOptimizer::Settings::NormalizationStrategy::Epsilon;
   settings_optimization.normalizationParameter = 0.0003;
   settings_optimization.solutionAccuracy = 0.01;
 

src/reducedmodeloptimizer.cpp

@@ -1,6 +1,7 @@
 #include "reducedmodeloptimizer.hpp"
-#include "flatpattern.hpp"
+#include "linearsimulationmodel.hpp"
 #include "simulationhistoryplotter.hpp"
+#include "trianglepatterngeometry.hpp"
 #include "trianglepattterntopology.hpp"
 #include <chrono>
 #include <dlib/global_optimization.h>
@@ -16,22 +17,17 @@ struct GlobalOptimizationVariables {
         reducedToFullInterfaceViMap;
     matplot::line_handle gPlotHandle;
     std::vector<double> gObjectiveValueHistory;
-  Eigen::Vector2d g_initialX;
-  std::unordered_set<size_t> reducedPatternExludedEdges;
     Eigen::VectorXd g_initialParameters;
     std::vector<ReducedModelOptimizer::SimulationScenario>
         simulationScenarioIndices;
     std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)};
+    double innerHexagonInitialRotationAngle{30};
     double g_innerHexagonInitialPos{0};
-  bool optimizeInnerHexagonSize{false};
-  std::vector<SimulationResults> firstOptimizationRoundResults;
-  int g_firstRoundIterationIndex{0};
     double minY{DBL_MAX};
     std::vector<double> minX;
-  std::vector<std::vector<double>> failedSimulationsXRatio;
     int numOfSimulationCrashes{false};
     int numberOfFunctionCalls{0};
-  int numberOfOptimizationParameters{3};
+    int numberOfOptimizationParameters{5};
     ReducedModelOptimizer::Settings optimizationSettings;
 } global;
 
@@ -131,17 +127,22 @@ void updateMesh(long n, const double *x) {
     //            e.secondBendingConstFactor
     //            << std::endl;
 
-  if (global.optimizeInnerHexagonSize) {
     assert(pReducedPatternSimulationMesh->EN() == 12);
-    for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN();
+    auto R = vcg::RotationMatrix(
-         vi += 2) {
+        ReducedModelOptimizer::patternPlaneNormal,
-      pReducedPatternSimulationMesh->vert[vi].P() =
+        vcg::math::ToRad(x[4] - global.innerHexagonInitialRotationAngle));
-          global.g_innerHexagonVectors[vi / 2] * x[n - 1];
+    //    for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN();
-    }
+    //         vi += 2) {
+    for (int rotationCounter = 0;
+         rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {
+        pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() =
+            R * global.g_innerHexagonVectors[rotationCounter] * x[3];
     }
 
     pReducedPatternSimulationMesh->reset();
+#ifdef POLYSCOPE_DEFINED
     pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
+#endif
 }
 
 double ReducedModelOptimizer::objective(double b, double r, double E) {
@@ -150,8 +151,9 @@ double ReducedModelOptimizer::objective(double b, double r, double E) {
 }
 
 double ReducedModelOptimizer::objective(double b, double h, double E,
-                                        double innerHexagonSize) {
+                                        double innerHexagonSize,
-  std::vector<double> x{b, h, E, innerHexagonSize};
+                                        double innerHexagonRotationAngle) {
+    std::vector<double> x{b, h, E, innerHexagonSize, innerHexagonRotationAngle};
     return ReducedModelOptimizer::objective(x.size(), x.data());
 }
 
@@ -172,24 +174,24 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
 
     // run simulations
     double totalError = 0;
-  FormFinder simulator;
+    LinearSimulationModel simulator;
-  FormFinder::Settings simulationSettings;
-  //    simulationSettings.shouldDraw = true;
     for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
-    //    std::cout << "Executing "
-    //              << simulationScenarioStrings[simulationScenarioIndex]
-    //              << std::endl;
         SimulationResults reducedModelResults = simulator.executeSimulation(
-        global.reducedPatternSimulationJobs[simulationScenarioIndex],
+            global.reducedPatternSimulationJobs[simulationScenarioIndex]);
-        simulationSettings);
+        //#ifdef POLYSCOPE_DEFINED
+        //    global.reducedPatternSimulationJobs[simulationScenarioIndex]
+        //        ->pMesh->registerForDrawing(colors.reducedInitial);
+        //    reducedModelResults.registerForDrawing(colors.reducedDeformed);
+        //    polyscope::show();
+        //    reducedModelResults.unregister();
+        //#endif
         //    std::string filename;
         if (!reducedModelResults.converged) {
             totalError += std::numeric_limits<double>::max();
             global.numOfSimulationCrashes++;
-      bool beVerbose = false; // TODO: use an extern or data member for that
+#ifdef POLYSCOPE_DEFINED
-      if (beVerbose) {
             std::cout << "Failed simulation" << std::endl;
-      }
+#endif
         } else {
             double simulationScenarioError = computeError(
                 reducedModelResults.displacements,
@@ -234,11 +236,12 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
         global.minY = totalError;
         global.minX.assign(x, x + n);
     }
-  // if (++global.numberOfFunctionCalls %100== 0) {
+#ifdef POLYSCOPE_DEFINED
-  // std::cout << "Number of function calls:" << global.numberOfFunctionCalls
+    if (++global.numberOfFunctionCalls % 100 == 0) {
-  //          << std::endl;
+        std::cout << "Number of function calls:" << global.numberOfFunctionCalls
-  //}
+                  << std::endl;
-
+    }
+#endif
     // compute error and return it
     //  global.gObjectiveValueHistory.push_back(error);
     //  auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
@@ -257,7 +260,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
 }
 
 void ReducedModelOptimizer::createSimulationMeshes(
-    FlatPattern &fullModel, FlatPattern &reducedModel,
+    PatternGeometry &fullModel, PatternGeometry &reducedModel,
     std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
     std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
     if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) {
@@ -278,8 +281,8 @@ void ReducedModelOptimizer::createSimulationMeshes(
     pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
 }
 
-void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel,
+void ReducedModelOptimizer::createSimulationMeshes(
-                                                   FlatPattern &reducedModel) {
+    PatternGeometry &fullModel, PatternGeometry &reducedModel) {
     ReducedModelOptimizer::createSimulationMeshes(
         fullModel, reducedModel, m_pFullPatternSimulationMesh,
         m_pReducedPatternSimulationMesh);
@@ -288,7 +291,7 @@ void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel,
 void ReducedModelOptimizer::computeMaps(
     const std::unordered_set<size_t> &reducedModelExcludedEdges,
     const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
-    FlatPattern &fullPattern, FlatPattern &reducedPattern,
+    PatternGeometry &fullPattern, PatternGeometry &reducedPattern,
     std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
         &reducedToFullInterfaceViMap,
     std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
@@ -305,7 +308,8 @@ void ReducedModelOptimizer::computeMaps(
     const size_t baseTriangleInterfaceVi =
         *(slotToNode.find(interfaceSlotIndex)->second.begin());
 
-  vcg::tri::Allocator<FlatPattern>::PointerUpdater<FlatPattern::VertexPointer>
+    vcg::tri::Allocator<PatternGeometry>::PointerUpdater<
+        PatternGeometry::VertexPointer>
         pu_fullModel;
     fullPattern.deleteDanglingVertices(pu_fullModel);
     const size_t fullModelBaseTriangleInterfaceVi =
@@ -322,19 +326,19 @@ void ReducedModelOptimizer::computeMaps(
 
     // Save excluded edges TODO:this changes the global object. Should this be a
     // function parameter?
-  global.reducedPatternExludedEdges.clear();
+    //  global.reducedPatternExludedEdges.clear();
-  const size_t fanSize = 6;
+    //  const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN();
-  const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN();
+    //  for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
-  for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
+    //    for (const size_t ei : reducedModelExcludedEdges) {
-    for (const size_t ei : reducedModelExcludedEdges) {
+    //      global.reducedPatternExludedEdges.insert(
-      global.reducedPatternExludedEdges.insert(
+    //          fanIndex * reducedBaseTriangleNumberOfEdges + ei);
-          fanIndex * reducedBaseTriangleNumberOfEdges + ei);
+    //    }
-    }
+    //  }
-  }
 
     // Construct reduced->full and full->reduced interface vi map
     reducedToFullInterfaceViMap.clear();
-  vcg::tri::Allocator<FlatPattern>::PointerUpdater<FlatPattern::VertexPointer>
+    vcg::tri::Allocator<PatternGeometry>::PointerUpdater<
+        PatternGeometry::VertexPointer>
         pu_reducedModel;
     reducedPattern.deleteDanglingVertices(pu_reducedModel);
     const size_t reducedModelBaseTriangleInterfaceVi =
@@ -365,17 +369,18 @@ void ReducedModelOptimizer::computeMaps(
     if (debugMapping) {
 #if POLYSCOPE_DEFINED
         reducedPattern.registerForDrawing();
-    std::vector<glm::vec3> colors_reducedPatternExcludedEdges(
+
-        reducedPattern.EN(), glm::vec3(0, 0, 0));
+        //    std::vector<glm::vec3> colors_reducedPatternExcludedEdges(
-    for (const size_t ei : global.reducedPatternExludedEdges) {
+        //        reducedPattern.EN(), glm::vec3(0, 0, 0));
-      colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0);
+        //    for (const size_t ei : global.reducedPatternExludedEdges) {
-    }
+        //      colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0);
-    const std::string label = reducedPattern.getLabel();
+        //    }
-    polyscope::getCurveNetwork(label)
+        //    const std::string label = reducedPattern.getLabel();
-        ->addEdgeColorQuantity("Excluded edges",
+        //    polyscope::getCurveNetwork(label)
-                               colors_reducedPatternExcludedEdges)
+        //        ->addEdgeColorQuantity("Excluded edges",
-        ->setEnabled(true);
+        //                               colors_reducedPatternExcludedEdges)
-    polyscope::show();
+        //        ->setEnabled(true);
+        //    polyscope::show();
 
         std::vector<glm::vec3> nodeColorsOpposite(fullPattern.VN(),
                                                   glm::vec3(0, 0, 0));
@@ -418,7 +423,7 @@ void ReducedModelOptimizer::computeMaps(
 }
 
 void ReducedModelOptimizer::computeMaps(
-    FlatPattern &fullPattern, FlatPattern &reducedPattern,
+    PatternGeometry &fullPattern, PatternGeometry &reducedPattern,
     const std::unordered_set<size_t> &reducedModelExcludedEdges) {
     ReducedModelOptimizer::computeMaps(
         reducedModelExcludedEdges, slotToNode, fullPattern, reducedPattern,
@@ -433,7 +438,7 @@ ReducedModelOptimizer::ReducedModelOptimizer(
 }
 
 void ReducedModelOptimizer::initializePatterns(
-    FlatPattern &fullPattern, FlatPattern &reducedPattern,
+    PatternGeometry &fullPattern, PatternGeometry &reducedPattern,
     const std::unordered_set<size_t> &reducedModelExcludedEdges) {
     // fullPattern.setLabel("full_pattern_" + fullPattern.getLabel());
     // reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel());
@@ -443,23 +448,22 @@ void ReducedModelOptimizer::initializePatterns(
     polyscope::removeAllStructures();
 #endif
     // Create copies of the input models
-  FlatPattern copyFullPattern;
+    PatternGeometry copyFullPattern;
-  FlatPattern copyReducedPattern;
+    PatternGeometry copyReducedPattern;
     copyFullPattern.copy(fullPattern);
     copyReducedPattern.copy(reducedPattern);
-  global.optimizeInnerHexagonSize = copyReducedPattern.EN() == 2;
+    /*
-  if (global.optimizeInnerHexagonSize) {
+   * Here we create the vector that connects the central node with the bottom
+   * 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++) {
-    const int fanSize = 6;
-    const CoordType rotationAxis(0, 0, 1);
-    for (int rotationCounter = 0; rotationCounter < fanSize;
-         rotationCounter++) {
         VectorType rotatedVector =
-          vcg::RotationMatrix(rotationAxis,
+            vcg::RotationMatrix(patternPlaneNormal,
                                 vcg::math::ToRad(rotationCounter * 60.0)) *
             baseTriangle_leftBottomNode;
         global.g_innerHexagonVectors[rotationCounter] = rotatedVector;
@@ -469,7 +473,9 @@ void ReducedModelOptimizer::initializePatterns(
     const double innerHexagonInitialPos_y =
         copyReducedPattern.vert[0].cP()[1] / global.g_innerHexagonVectors[0][1];
     global.g_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);
@@ -477,10 +483,8 @@ void ReducedModelOptimizer::initializePatterns(
 
 void ReducedModelOptimizer::initializeOptimizationParameters(
     const std::shared_ptr<SimulationMesh> &mesh) {
-  global.numberOfOptimizationParameters = 3;
+    global.numberOfOptimizationParameters = 5;
-  global.g_initialParameters.resize(
+    global.g_initialParameters.resize(global.numberOfOptimizationParameters);
-      global.optimizeInnerHexagonSize ? ++global.numberOfOptimizationParameters
-                                      : global.numberOfOptimizationParameters);
     // Save save the beam stiffnesses
     //  for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) {
     //    Element &e = pReducedModelElementalMesh->elements[ei];
@@ -493,13 +497,12 @@ void ReducedModelOptimizer::initializeOptimizationParameters(
     //    }
     const double initialB = std::sqrt(mesh->elements[0].A);
     const double initialRatio = 1;
-  ;
     global.g_initialParameters(0) = initialB;
     global.g_initialParameters(1) = initialRatio;
     global.g_initialParameters(2) = mesh->elements[0].material.youngsModulus;
-  if (global.optimizeInnerHexagonSize) {
     global.g_initialParameters(3) = global.g_innerHexagonInitialPos;
-  }
+    global.innerHexagonInitialRotationAngle = 30;
+    global.g_initialParameters(4) = global.innerHexagonInitialRotationAngle;
     //  g_initialParameters =
     //      m_pReducedPatternSimulationMesh->elements[0].properties.E;
     //  for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) {
@@ -680,19 +683,11 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
 
     auto start = std::chrono::system_clock::now();
     dlib::function_evaluation result;
-  if (global.optimizeInnerHexagonSize) {
+    double (*objF)(double, double, double, double, double) = &objective;
-    double (*objF)(double, double, double, double) = &objective;
     result = dlib::find_min_global(
         objF, xMin, xMax,
         dlib::max_function_calls(settings.numberOfFunctionCalls),
         std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
-  } else {
-    double (*objF)(double, double, double) = &objective;
-    result = dlib::find_min_global(
-        objF, xMin, xMax,
-        dlib::max_function_calls(settings.numberOfFunctionCalls),
-        std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
-  }
     auto end = std::chrono::system_clock::now();
     auto elapsed =
         std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
@@ -711,13 +706,7 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
         auto temp = global.optimizationSettings.normalizationStrategy;
         global.optimizationSettings.normalizationStrategy =
             Settings::NormalizationStrategy::NonNormalized;
-    if (global.optimizeInnerHexagonSize) {
+        results.rawObjectiveValue = objective(results.x.size(), results.x.data());
-      results.rawObjectiveValue = objective(4, results.x.data());
-
-    } else {
-      results.rawObjectiveValue =
-          objective(results.x[0], results.x[1], results.x[2]);
-    }
         global.optimizationSettings.normalizationStrategy = temp;
 
     } else {
@@ -765,10 +754,6 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
     return results;
 }
 
-void ReducedModelOptimizer::setInitialGuess(std::vector<double> v) {
-  initialGuess = v;
-}
-
 std::vector<std::shared_ptr<SimulationJob>>
 ReducedModelOptimizer::createScenarios(
     const std::shared_ptr<SimulationMesh> &pMesh) {
@@ -777,9 +762,6 @@ ReducedModelOptimizer::createScenarios(
     std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
     std::unordered_map<VertexIndex, Vector6d> nodalForces;
     const double forceMagnitude = 1;
-  // Assuming the patterns lays on the x-y plane
-  const CoordType patternPlaneNormal(0, 0, 1);
-  // Assumes that the first interface node lays on the y axis
 
     //// Axial
     SimulationScenario scenarioName = SimulationScenario::Axial;
@@ -899,12 +881,12 @@ ReducedModelOptimizer::createScenarios(
             fixedVertices[viPair.second] = std::unordered_set<DoFType>{2};
         }
         CoordType v =
-        (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
+            (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) ^
-            .Normalize();
+            CoordType(0, 0, -1).Normalize();
         nodalForces[viPair.first] =
-        Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude * 0.1;
+            Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude * 1;
         nodalForces[viPair.second] =
-        Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 0.1;
+            Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 1;
     }
     scenarios[scenarioName] = std::make_shared<SimulationJob>(
         SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
@@ -918,21 +900,21 @@ ReducedModelOptimizer::createScenarios(
          viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) {
         const auto &viPair = *viPairIt;
         CoordType v =
-        (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
+            (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) ^
-            .Normalize();
+            CoordType(0, 0, -1).Normalize();
         if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
             nodalForces[viPair.first] =
-          Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.02 * forceMagnitude;
+                Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.2 * forceMagnitude;
             nodalForces[viPair.second] =
-          Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.02 * forceMagnitude;
+                Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.2 * 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}) * 0.01 * forceMagnitude;
+                Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.1 * forceMagnitude;
             nodalForces[viPair.second] =
-          Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.01 * forceMagnitude;
+                Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.1 * forceMagnitude;
         }
     }
     scenarios[scenarioName] = std::make_shared<SimulationJob>(
@@ -943,7 +925,6 @@ ReducedModelOptimizer::createScenarios(
 }
 
 void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
-  assert(global.optimizationSettings.shouldNormalizeObjectiveValue);
     if (global.optimizationSettings.normalizationStrategy ==
         Settings::NormalizationStrategy::Epsilon) {
 
@@ -997,7 +978,6 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
     global.reducedPatternSimulationJobs.resize(NumberOfSimulationScenarios);
     global.fullPatternDisplacements.resize(NumberOfSimulationScenarios);
     global.objectiveNormalizationValues.resize(NumberOfSimulationScenarios);
-  global.g_firstRoundIterationIndex = 0;
     global.minY = std::numeric_limits<double>::max();
     global.numOfSimulationCrashes = 0;
     global.numberOfFunctionCalls = 0;

src/reducedmodeloptimizer.hpp

@@ -5,6 +5,7 @@
 #include "csvfile.hpp"
 #include "edgemesh.hpp"
 #include "elementalmesh.hpp"
+#include "linearsimulationmodel.hpp"
 #include "matplot/matplot.h"
 #include <Eigen/Dense>
 
@@ -23,7 +24,6 @@ class ReducedModelOptimizer {
       m_fullPatternOppositeInterfaceViMap;
   std::unordered_map<size_t, size_t> nodeToSlot;
   std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
-  std::vector<double> initialGuess;
 #ifdef POLYSCOPE_DEFINED
   struct StaticColors {
     glm::vec3 fullInitial;
@@ -41,6 +41,8 @@ class ReducedModelOptimizer {
 #endif // POLYSCOPE_DEFINED
 
 public:
+  inline static int fanSize{6};
+  inline static VectorType patternPlaneNormal{0, 0, 1};
   enum SimulationScenario {
     Axial,
     Shear,
@@ -49,7 +51,6 @@ public:
     Saddle,
     NumberOfSimulationScenarios
   };
-
       struct xRange{
       std::string label;
     double min;
@@ -62,9 +63,14 @@ public:
   struct Results;
 
   struct Settings {
-    enum NormalizationStrategy { NonNormalized, Epsilon };
+    enum NormalizationStrategy {
+      NonNormalized,
+      Epsilon,
+      MaxDisplacement,
+      EqualDisplacements
+    };
     inline static vector<std::string> normalizationStrategyStrings{
-        "NonNormalized", "Epsilon"};
+        "NonNormalized", "Epsilon", "MaxDsiplacement", "EqualDisplacements"};
     std::vector<xRange> xRanges;
     int numberOfFunctionCalls{100};
     double solutionAccuracy{1e-2};
@@ -110,11 +116,8 @@ public:
       }
       os << numberOfFunctionCalls;
       os << solutionAccuracy;
-      if (normalizationStrategy == Epsilon) {
+      os << normalizationStrategyStrings[normalizationStrategy] + "_" +
-        os << "Epsilon_" + std::to_string(normalizationParameter);
+                std::to_string(normalizationParameter);
-      } else {
-        os << "NonNormalized";
-      }
     }
   };
   struct Results {
@@ -194,12 +197,13 @@ public:
     void draw() const {
       initPolyscope();
       FormFinder simulator;
+      LinearSimulationModel linearSimulator;
       assert(fullPatternSimulationJobs.size() ==
              reducedPatternSimulationJobs.size());
       fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
           colors.fullInitial);
       reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(
-          colors.reducedInitial);
+          colors.reducedInitial, false);
 
       const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
       for (int simulationJobIndex = 0;
@@ -212,14 +216,24 @@ public:
         SimulationResults fullModelResults =
             simulator.executeSimulation(pFullPatternSimulationJob);
         fullModelResults.registerForDrawing(colors.fullDeformed);
+        SimulationResults fullModelLinearResults =
+            linearSimulator.executeSimulation(pFullPatternSimulationJob);
+        fullModelLinearResults.setLabelPrefix("linear");
+        fullModelLinearResults.registerForDrawing(colors.fullDeformed, false);
         // Drawing of reduced pattern results
         const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
             reducedPatternSimulationJobs[simulationJobIndex];
         SimulationResults reducedModelResults =
             simulator.executeSimulation(pReducedPatternSimulationJob);
         reducedModelResults.registerForDrawing(colors.reducedDeformed);
+        SimulationResults reducedModelLinearResults =
+            linearSimulator.executeSimulation(pReducedPatternSimulationJob);
+        reducedModelLinearResults.setLabelPrefix("linear");
+        reducedModelLinearResults.registerForDrawing(colors.reducedDeformed,
+                                                     false);
         polyscope::options::programName =
             fullPatternSimulationJobs[0]->pMesh->getLabel();
+        polyscope::view::resetCameraToHomeView();
         polyscope::show();
         // Save a screensh
         const std::string screenshotFilename =
@@ -230,12 +244,14 @@ public:
         polyscope::screenshot(screenshotFilename, false);
         fullModelResults.unregister();
         reducedModelResults.unregister();
+        reducedModelLinearResults.unregister();
+        fullModelLinearResults.unregister();
         //                    double error = computeError(
         //                        reducedModelResults.displacements,fullModelResults.displacements,
         //                        );
         //                    std::cout << "Error of simulation scenario "
         //                              <<
-        //                              simulationScenarioStrings[simulationScenarioIndex]
+        //                              simula         simulationScenarioStrings[simulationScenarioIndex]
         //                              << " is "
         //                              << error << std::endl;
       }
@@ -333,31 +349,29 @@ public:
   getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
 
   void initializePatterns(
-      FlatPattern &fullPattern, FlatPattern &reducedPatterm,
+      PatternGeometry &fullPattern, PatternGeometry &reducedPatterm,
       const std::unordered_set<size_t> &reducedModelExcludedEges);
 
-  void setInitialGuess(std::vector<double> v);
-
-  static void runBeamOptimization();
-
   static void runSimulation(const std::string &filename,
                             std::vector<double> &x);
 
-  static double objective(double x0, double x1, double x2, double x3);
+  static double objective(double x0, double x1, double x2, double x3,
+                          double innerHexagonRotationAngle);
   static double objective(double b, double r, double E);
 
   static std::vector<std::shared_ptr<SimulationJob>>
     createScenarios(const std::shared_ptr<SimulationMesh> &pMesh,
                     const std::unordered_map<size_t, size_t>
                         &fullPatternOppositeInterfaceViMap);
+
   static void createSimulationMeshes(
-      FlatPattern &fullModel, FlatPattern &reducedModel,
+      PatternGeometry &fullModel, PatternGeometry &reducedModel,
       std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
       std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh);
   static void computeMaps(
       const std::unordered_set<size_t> &reducedModelExcludedEdges,
       const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
-      FlatPattern &fullPattern, FlatPattern &reducedPattern,
+      PatternGeometry &fullPattern, PatternGeometry &reducedPattern,
       std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
           &reducedToFullInterfaceViMap,
       std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
@@ -403,10 +417,10 @@ private:
   static Results runOptimization(const Settings &settings);
   std::vector<std::shared_ptr<SimulationJob>>
   createScenarios(const std::shared_ptr<SimulationMesh> &pMesh);
-  void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern,
+  void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern,
                    const std::unordered_set<size_t> &reducedModelExcludedEges);
-  void createSimulationMeshes(FlatPattern &fullModel,
+  void createSimulationMeshes(PatternGeometry &fullModel,
-                              FlatPattern &reducedModel);
+                              PatternGeometry &reducedModel);
   static void
   initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh);