Refactoring

beamformfinder.cpp

@@ -15,9 +15,9 @@ void FormFinder::runUnitTests() {
   VCGEdgeMesh beam;
   //  const size_t spanGridSize = 11;
   //  mesh.createSpanGrid(spanGridSize);
-  beam.loadPly(std::filesystem::path(groundOfTruthFolder)
+  beam.load(std::filesystem::path(groundOfTruthFolder)
-                   .append("simpleBeam.ply")
+                .append("simpleBeam.ply")
-                   .string());
+                .string());
   std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
   fixedVertices[0] = std::unordered_set<DoFType>{0, 1, 2, 3};
   fixedVertices[beam.VN() - 1] = std::unordered_set<DoFType>{1, 2};
@@ -65,9 +65,9 @@ void FormFinder::runUnitTests() {
   VCGEdgeMesh shortSpanGrid;
   //  const size_t spanGridSize = 11;
   //  mesh.createSpanGrid(spanGridSize);
-  shortSpanGrid.loadPly(std::filesystem::path(groundOfTruthFolder)
+  shortSpanGrid.load(std::filesystem::path(groundOfTruthFolder)
-                            .append("shortSpanGridshell.ply")
+                         .append("shortSpanGridshell.ply")
-                            .string());
+                         .string());
 
   fixedVertices.clear();
   //// Corner nodes
@@ -123,9 +123,9 @@ void FormFinder::runUnitTests() {
 
   // Third example of the paper
   VCGEdgeMesh longSpanGrid;
-  longSpanGrid.loadPly(std::filesystem::path(groundOfTruthFolder)
+  longSpanGrid.load(std::filesystem::path(groundOfTruthFolder)
-                           .append("longSpanGridshell.ply")
+                        .append("longSpanGridshell.ply")
-                           .string());
+                        .string());
   const size_t spanGridSize = 11;
 
   fixedVertices.clear();
@@ -793,7 +793,7 @@ void FormFinder::updateResidualForcesOnTheFly(
           std::vector<std::pair<int, Vector6d>>(4, {-1, Vector6d()}));
   //  omp_lock_t writelock;
   //  omp_init_lock(&writelock);
-  //#pragma omp parallel for //schedule(static) num_threads(8)
+  //#pragma omp parallel for schedule(static) num_threads(8)
   for (int ei = 0; ei < pMesh->EN(); ei++) {
     const EdgeType &e = pMesh->edge[ei];
     const SimulationMesh::VertexType &ev_j = *e.cV(0);
@@ -1019,10 +1019,23 @@ void FormFinder::updateNodalExternalForces(
     }
     externalMomentsNorm += std::sqrt(pow(nodalExternalForce[3], 2) +
                                      pow(nodalExternalForce[4], 2));
-    //    CoordType v = (mesh->vert[nodeIndex].cP() -
+
-    //    mesh->vert[0].cP()).Normalize(); const double forceMagnitude = 0.1;
+    /*
-    //    nodalExternalForce[3] = v[0] * forceMagnitude;
+     * The external moments are given as a rotation around an axis.
-    //    nodalExternalForce[4] = v[1] * forceMagnitude;
+     * In this implementation we model moments as rotation of the normal vector
+     * and because of that we need to transform the moments.
+     */
+
+    if (externalMomentsNorm != 0) {
+      VectorType momentVector(
+          nodalExternalForce[3], nodalExternalForce[4],
+          nodalExternalForce[5]); // rotation around this vector
+      VectorType transformedVector = momentVector ^ VectorType(0, 0, 1);
+      nodalExternalForce[3] = transformedVector[0];
+      nodalExternalForce[4] = transformedVector[1];
+      nodalExternalForce[5] = transformedVector[5];
+    }
+
     node.force.external = nodalExternalForce;
   }
 }
@@ -1162,26 +1175,26 @@ void FormFinder::updateNodalMasses() {
       assert(rotationalSumSk_I3 != 0);
       assert(rotationalSumSk_J != 0);
     }
-    pMesh->nodes[v].translationalMass =
+    pMesh->nodes[v].mass.translational =
         gamma * pow(mSettings.Dtini, 2) * 2 * translationalSumSk;
-    pMesh->nodes[v].rotationalMass_I2 =
+    pMesh->nodes[v].mass.rotationalI2 =
         gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I2;
-    pMesh->nodes[v].rotationalMass_I3 =
+    pMesh->nodes[v].mass.rotationalI3 =
         gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I3;
-    pMesh->nodes[v].rotationalMass_J =
+    pMesh->nodes[v].mass.rotationalJ =
         gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_J;
 
     assert(std::pow(mSettings.Dtini, 2.0) * translationalSumSk /
-               pMesh->nodes[v].translationalMass <
+               pMesh->nodes[v].mass.translational <
            2);
     assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_I2 /
-               pMesh->nodes[v].rotationalMass_I2 <
+               pMesh->nodes[v].mass.rotationalI2 <
            2);
     assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_I3 /
-               pMesh->nodes[v].rotationalMass_I3 <
+               pMesh->nodes[v].mass.rotationalI3 <
            2);
     assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_J /
-               pMesh->nodes[v].rotationalMass_J <
+               pMesh->nodes[v].mass.rotationalJ <
            2);
   }
 }
@@ -1193,16 +1206,16 @@ void FormFinder::updateNodalAccelerations() {
     for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
       if (dofi == DoF::Ux || dofi == DoF::Uy || dofi == DoF::Uz) {
         node.acceleration[dofi] =
-            node.force.residual[dofi] / node.translationalMass;
+            node.force.residual[dofi] / node.mass.translational;
       } else if (dofi == DoF::Nx) {
         node.acceleration[dofi] =
-            node.force.residual[dofi] / node.rotationalMass_J;
+            node.force.residual[dofi] / node.mass.rotationalJ;
       } else if (dofi == DoF::Ny) {
         node.acceleration[dofi] =
-            node.force.residual[dofi] / node.rotationalMass_I3;
+            node.force.residual[dofi] / node.mass.rotationalI3;
       } else if (dofi == DoF::Nr) {
         node.acceleration[dofi] =
-            node.force.residual[dofi] / node.rotationalMass_I2;
+            node.force.residual[dofi] / node.mass.rotationalI2;
       }
     }
   }
@@ -1239,7 +1252,6 @@ void FormFinder::updateNodePosition(
     const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
         &fixedVertices) {
   Node &node = pMesh->nodes[v];
-  CoordType previousLocation = v.cP();
   const VertexIndex &vi = pMesh->nodes[v].vi;
 
   VectorType displacementVector(0, 0, 0);
@@ -1252,7 +1264,6 @@ void FormFinder::updateNodePosition(
   if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(2)) {
     displacementVector += VectorType(0, 0, node.displacements[2]);
   }
-  node.previousLocation = previousLocation;
   v.P() = node.initialLocation + displacementVector;
   if (shouldApplyInitialDistortion && mCurrentSimulationStep < 40) {
     VectorType desiredInitialDisplacement(0, 0, 0.01);
@@ -1411,16 +1422,16 @@ void FormFinder::updateKineticEnergy() {
         std::pow(node.velocity[0], 2) + std::pow(node.velocity[1], 2) +
         std::pow(node.velocity[2], 2));
     const double nodeTranslationalKineticEnergy =
-        0.5 * node.translationalMass * pow(translationalVelocityNorm, 2);
+        0.5 * node.mass.translational * pow(translationalVelocityNorm, 2);
 
     const double nodeRotationalKineticEnergy =
-        0.5 * (node.rotationalMass_J * pow(node.velocity[3], 2) +
+        0.5 * (node.mass.rotationalJ * pow(node.velocity[3], 2) +
-               +node.rotationalMass_I3 * pow(node.velocity[4], 2) +
+               +node.mass.rotationalI3 * pow(node.velocity[4], 2) +
-               +node.rotationalMass_I2 * pow(node.velocity[5], 2));
+               +node.mass.rotationalI2 * pow(node.velocity[5], 2));
 
     node.kineticEnergy +=
         nodeTranslationalKineticEnergy /*+ nodeRotationalKineticEnergy*/;
-    assert(node.kineticEnergy < 10000000000000000000);
+    assert(node.kineticEnergy < 1e15);
 
     pMesh->currentTotalKineticEnergy += node.kineticEnergy;
     pMesh->currentTotalTranslationalKineticEnergy +=
@@ -1472,13 +1483,13 @@ void FormFinder::updatePositionsOnTheFly(
       //      assert(rotationalSumSk_I3 != 0);
       //      assert(rotationalSumSk_J != 0);
     }
-    pMesh->nodes[v].translationalMass =
+    pMesh->nodes[v].mass.translational =
         gamma * pow(mSettings.Dtini, 2) * 2 * translationalSumSk;
-    pMesh->nodes[v].rotationalMass_I2 =
+    pMesh->nodes[v].mass.rotationalI2 =
         gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I2;
-    pMesh->nodes[v].rotationalMass_I3 =
+    pMesh->nodes[v].mass.rotationalI3 =
         gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I3;
-    pMesh->nodes[v].rotationalMass_J =
+    pMesh->nodes[v].mass.rotationalJ =
         gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_J;
 
     //    assert(std::pow(mSettings.Dtini, 2.0) * translationalSumSk /
@@ -1500,16 +1511,16 @@ void FormFinder::updatePositionsOnTheFly(
     for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
       if (dofi == DoF::Ux || dofi == DoF::Uy || dofi == DoF::Uz) {
         node.acceleration[dofi] =
-            node.force.residual[dofi] / node.translationalMass;
+            node.force.residual[dofi] / node.mass.translational;
       } else if (dofi == DoF::Nx) {
         node.acceleration[dofi] =
-            node.force.residual[dofi] / node.rotationalMass_J;
+            node.force.residual[dofi] / node.mass.rotationalJ;
       } else if (dofi == DoF::Ny) {
         node.acceleration[dofi] =
-            node.force.residual[dofi] / node.rotationalMass_I3;
+            node.force.residual[dofi] / node.mass.rotationalI3;
       } else if (dofi == DoF::Nr) {
         node.acceleration[dofi] =
-            node.force.residual[dofi] / node.rotationalMass_I2;
+            node.force.residual[dofi] / node.mass.rotationalI2;
       }
     }
   }
@@ -1789,4 +1800,4 @@ mesh->currentTotalPotentialEnergykN*/
     }
   }
   return results;
-}
+}

beamformfinder.hpp

@@ -3,8 +3,6 @@
 
 #include "elementalmesh.hpp"
 #include "matplot/matplot.h"
-//#include "polyscope/curve_network.h"
-//#include "polyscope/polyscope.h"
 #include "simulationresult.hpp"
 #include <Eigen/Dense>
 #include <filesystem>

edgemesh.cpp

@@ -1,5 +1,6 @@
 #include "edgemesh.hpp"
 #include "vcg/simplex/face/topology.h"
+#include <wrap/nanoply/include/nanoplyWrapper.hpp>
 
 Eigen::MatrixX2i VCGEdgeMesh::getEigenEdges() const { return eigenEdges; }
 
@@ -45,7 +46,8 @@ void VCGEdgeMesh::GeneratedRegularSquaredPattern(
     for (int k = 0; k <= samplesNo; k++) {
       const double t = double(k) / samplesNo;
       const double a = (1 - t) * angle;
-      // const double r = vcg::math::Sin(t*M_PI_2) /*(1-((1-t)*(1-t)))*/ * 0.5;
+      // const double r = vcg::math::Sin(t*M_PI_2) /*(1-((1-t)*(1-t)))*/ *
+      // 0.5;
       const double r = t * 0.5; // linear
 
       vcg::Point2d p(vcg::math::Cos(a), vcg::math::Sin(a));
@@ -174,7 +176,7 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth,
   return true;
 }
 
-bool VCGEdgeMesh::loadPly(const std::string plyFilename) {
+bool VCGEdgeMesh::load(const string &plyFilename) {
 
   std::string usedPath = plyFilename;
   if (std::filesystem::path(plyFilename).is_relative()) {
@@ -182,7 +184,6 @@ bool VCGEdgeMesh::loadPly(const std::string plyFilename) {
   }
   assert(std::filesystem::exists(usedPath));
   this->Clear();
-  const bool useDefaultImporter = false;
   if (!loadUsingNanoply(usedPath)) {
     std::cerr << "Error: Unable to open " + usedPath << std::endl;
     return false;
@@ -217,7 +218,8 @@ bool VCGEdgeMesh::loadUsingNanoply(const std::string &plyFilename) {
   return true;
 }
 
-// bool VCGEdgeMesh::plyFileHasAllRequiredFields(const std::string &plyFilename)
+// bool VCGEdgeMesh::plyFileHasAllRequiredFields(const std::string
+// &plyFilename)
 // {
 //  const nanoply::Info info(plyFilename);
 //  const std::vector<nanoply::PlyElement>::const_iterator edgeElemIt =
@@ -238,28 +240,6 @@ bool VCGEdgeMesh::loadUsingNanoply(const std::string &plyFilename) {
 //                            plyPropertyBeamMaterialID);
 //}
 
-bool VCGEdgeMesh::hasPlyEdgeProperty(
-    const std::string &plyFilename,
-    const std::vector<nanoply::PlyProperty> &edgeProperties,
-    const std::string &plyEdgePropertyName) {
-  const bool hasEdgeProperty = hasProperty(edgeProperties, plyEdgePropertyName);
-  if (!hasEdgeProperty) {
-    std::cerr << "Ply file " + plyFilename +
-                     " is missing the propertry:" + plyEdgePropertyName
-              << std::endl;
-    return false;
-  }
-  return true;
-}
-
-bool VCGEdgeMesh::hasProperty(const std::vector<nanoply::PlyProperty> &v,
-                              const std::string &propertyName) {
-  return v.end() != std::find_if(v.begin(), v.end(),
-                                 [&](const nanoply::PlyProperty &plyProperty) {
-                                   return plyProperty.name == propertyName;
-                                 });
-}
-
 Eigen::MatrixX3d VCGEdgeMesh::getNormals() const {
   Eigen::MatrixX3d vertexNormals;
   vertexNormals.resize(VN(), 3);
@@ -289,8 +269,10 @@ void VCGEdgeMesh::getEdges(Eigen::MatrixX3d &edgeStartingPoints,
 VCGEdgeMesh::VCGEdgeMesh() {}
 
 void VCGEdgeMesh::updateEigenEdgeAndVertices() {
+#ifdef POLYSCOPE_DEFINED
   getEdges(eigenEdges);
   getVertices(eigenVertices);
+#endif
 }
 
 void VCGEdgeMesh::copy(VCGEdgeMesh &mesh) {
@@ -329,9 +311,15 @@ void VCGEdgeMesh::getEdges(Eigen::MatrixX2i &edges) {
   edges.resize(EN(), 2);
   for (int edgeIndex = 0; edgeIndex < EN(); edgeIndex++) {
     const VCGEdgeMesh::EdgeType &edge = this->edge[edgeIndex];
-    const size_t nodeIndex0 = vcg::tri::Index<VCGEdgeMesh>(*this, edge.cV(0));
+    assert(!edge.IsD());
-    const size_t nodeIndex1 = vcg::tri::Index<VCGEdgeMesh>(*this, edge.cV(1));
+    auto vp0 = edge.cV(0);
-    edges.row(edgeIndex) = Eigen::Vector2i(nodeIndex0, nodeIndex1);
+    auto vp1 = edge.cV(1);
+    assert(vcg::tri::IsValidPointer(*this, vp0) &&
+           vcg::tri::IsValidPointer(*this, vp1));
+    const size_t vi0 = vcg::tri::Index<VCGEdgeMesh>(*this, vp0);
+    const size_t vi1 = vcg::tri::Index<VCGEdgeMesh>(*this, vp1);
+    assert(vi0 != -1 && vi1 != -1);
+    edges.row(edgeIndex) = Eigen::Vector2i(vi0, vi1);
   }
 }
 
@@ -342,13 +330,14 @@ void VCGEdgeMesh::printVertexCoordinates(const size_t &vi) const {
 
 #ifdef POLYSCOPE_DEFINED
 void VCGEdgeMesh::registerForDrawing(
-    const std::optional<glm::vec3> &desiredColor) const {
+    const std::optional<glm::vec3> &desiredColor,
+    const bool &shouldEnable) const {
   initPolyscope();
   const double drawingRadius = 0.002;
-  auto polyscopeHandle_edgeMesh =
+  auto polyscopeHandle_edgeMesh = polyscope::registerCurveNetwork(
-      polyscope::registerCurveNetwork(label, getEigenVertices(),
+      label, getEigenVertices(), getEigenEdges());
-                                      getEigenEdges())
+  polyscopeHandle_edgeMesh->setEnabled(shouldEnable);
-          ->setRadius(drawingRadius, true);
+  polyscopeHandle_edgeMesh->setRadius(drawingRadius, true);
   if (desiredColor.has_value()) {
     polyscopeHandle_edgeMesh->setColor(desiredColor.value());
   }

edgemesh.hpp

@@ -1,14 +1,13 @@
 #ifndef EDGEMESH_HPP
 #define EDGEMESH_HPP
 #include "beam.hpp"
-#include "externvariables.hpp"
+#include "mesh.hpp"
-#include "polymesh.hpp"
+//#include "polymesh.hpp"
 #include "utilities.hpp"
-#include "vcgtrimesh.hpp"
+//#include "vcgtrimesh.hpp"
 #include <vcg/complex/complex.h>
 #include <vector>
 #include <wrap/io_trimesh/import.h>
-#include <wrap/nanoply/include/nanoplyWrapper.hpp>
 
 using EdgeIndex = size_t;
 
@@ -29,7 +28,8 @@ class VCGEdgeMeshEdgeType
                        vcg::edge::VEAdj> {};
 
 class VCGEdgeMesh : public vcg::tri::TriMesh<std::vector<VCGEdgeMeshVertexType>,
-                                             std::vector<VCGEdgeMeshEdgeType>> {
+                                             std::vector<VCGEdgeMeshEdgeType>>,
+                    Mesh {
 
 protected:
   Eigen::MatrixX2i eigenEdges;
@@ -54,19 +54,11 @@ public:
 
   Eigen::MatrixX3d getNormals() const;
 
-  bool hasProperty(const std::vector<nanoply::PlyProperty> &v,
-                   const std::string &propertyName);
-
-  bool
-  hasPlyEdgeProperty(const std::string &plyFilename,
-                     const std::vector<nanoply::PlyProperty> &edgeProperties,
-                     const std::string &plyEdgePropertyName);
-
   bool plyFileHasAllRequiredFields(const std::string &plyFilename);
 
   bool loadUsingNanoply(const std::string &plyFilename);
 
-  virtual bool loadPly(const std::string plyFilename);
+  bool load(const std::string &plyFilename) override;
 
   bool savePly(const std::string plyFilename);
 
@@ -80,7 +72,8 @@ public:
   void printVertexCoordinates(const size_t &vi) const;
 #ifdef POLYSCOPE_DEFINED
   void registerForDrawing(
-      const std::optional<glm::vec3> &desiredColor = std::nullopt) const;
+      const std::optional<glm::vec3> &desiredColor = std::nullopt,
+      const bool &shouldEnable = true) const;
   void unregister() const;
 #endif
   std::string getLabel() const;

elementalmesh.cpp

@@ -1,4 +1,5 @@
 #include "elementalmesh.hpp"
+#include <wrap/nanoply/include/nanoplyWrapper.hpp>
 
 SimulationMesh::SimulationMesh() {
   elements = vcg::tri::Allocator<VCGEdgeMesh>::GetPerEdgeAttribute<Element>(
@@ -45,8 +46,8 @@ SimulationMesh::~SimulationMesh() {
                                                                    nodes);
 }
 
-SimulationMesh::SimulationMesh(FlatPattern &pattern) {
+SimulationMesh::SimulationMesh(PatternGeometry &pattern) {
-  vcg::tri::MeshAssert<FlatPattern>::VertexNormalNormalized(pattern);
+  vcg::tri::MeshAssert<PatternGeometry>::VertexNormalNormalized(pattern);
 
   vcg::tri::Append<VCGEdgeMesh, ConstVCGEdgeMesh>::MeshCopy(*this, pattern);
   elements = vcg::tri::Allocator<VCGEdgeMesh>::GetPerEdgeAttribute<Element>(
@@ -218,6 +219,7 @@ void SimulationMesh::initializeElements() {
     element.derivativeR_j.resize(6);
     element.derivativeR_jplus1.resize(6);
   }
+  updateElementalFrames();
 }
 
 void SimulationMesh::updateElementalLengths() {
@@ -235,6 +237,14 @@ void SimulationMesh::updateElementalLengths() {
   }
 }
 
+void SimulationMesh::updateElementalFrames() {
+  for (const EdgeType &e : edge) {
+    const VectorType elementNormal =
+        (e.cV(0)->cN() + e.cV(1)->cN()).Normalize();
+    elements[e].frame = computeElementFrame(e.cP(0), e.cP(1), elementNormal);
+  }
+}
+
 void SimulationMesh::setBeamCrossSection(
     const CrossSectionType &beamDimensions) {
   for (size_t ei = 0; ei < EN(); ei++) {
@@ -267,7 +277,7 @@ std::vector<ElementMaterial> SimulationMesh::getBeamMaterial() {
   return beamMaterial;
 }
 
-bool SimulationMesh::loadPly(const string &plyFilename) {
+bool SimulationMesh::load(const string &plyFilename) {
   this->Clear();
   //  assert(plyFileHasAllRequiredFields(plyFilename));
   // Load the ply file

elementalmesh.hpp

@@ -3,7 +3,7 @@
 
 #include "Eigen/Dense"
 #include "edgemesh.hpp"
-#include "flatpattern.hpp"
+#include "trianglepatterngeometry.hpp"
 
 struct Element;
 struct Node;
@@ -24,14 +24,14 @@ public:
   PerEdgeAttributeHandle<Element> elements;
   PerVertexAttributeHandle<Node> nodes;
   ~SimulationMesh();
-  SimulationMesh(FlatPattern &pattern);
+  SimulationMesh(PatternGeometry &pattern);
   SimulationMesh(ConstVCGEdgeMesh &edgeMesh);
   SimulationMesh(SimulationMesh &elementalMesh);
   void updateElementalLengths();
 
   std::vector<VCGEdgeMesh::EdgePointer>
   getIncidentElements(const VCGEdgeMesh::VertexType &v);
-  virtual bool loadPly(const string &plyFilename);
+  virtual bool load(const string &plyFilename);
   std::vector<CrossSectionType> getBeamDimensions();
   std::vector<ElementMaterial> getBeamMaterial();
 
@@ -47,6 +47,7 @@ public:
   void setBeamMaterial(const double &pr, const double &ym);
   void reset();
   SimulationMesh();
+  void updateElementalFrames();
 };
 
 struct Element {
@@ -131,14 +132,17 @@ struct Node {
     bool hasExternalForce() const { return external.isZero(); }
   };
 
+  struct Mass {
+    double translational;
+    double rotationalI2;
+    double rotationalI3;
+    double rotationalJ;
+  };
+
+  Mass mass;
   VertexIndex vi;
   CoordType initialLocation;
-  CoordType previousLocation;
   CoordType initialNormal;
-  double translationalMass;
-  double rotationalMass_I2;
-  double rotationalMass_I3;
-  double rotationalMass_J;
   Vector6d acceleration{0};
   Forces force;
   Vector6d velocity{0};

flatpattern.cpp

@@ -1,424 +0,0 @@
-#include "flatpattern.hpp"
-#include "trianglepatterngeometry.hpp"
-#include <filesystem>
-
-FlatPattern::FlatPattern() {}
-
-FlatPattern::FlatPattern(const std::string &filename, bool addNormalsIfAbsent) {
-  if (!std::filesystem::exists(std::filesystem::path(filename))) {
-    assert(false);
-    std::cerr << "No flat pattern with name " << filename << std::endl;
-    return;
-  }
-  if (!loadPly(filename)) {
-    assert(false);
-    std::cerr << "File could not be loaded  " << filename << std::endl;
-    return;
-  }
-  if (addNormalsIfAbsent) {
-    bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001;
-    if (normalsAreAbsent) {
-      for (auto &v : vert) {
-        v.N() = CoordType(0, 0, 1);
-      }
-    }
-  }
-
-  vcg::tri::UpdateTopology<FlatPattern>::VertexEdge(*this);
-  const double baseTriangleCentralEdgeSize =
-      (vert[0].cP() - vert[3].cP()).Norm();
-
-  updateEigenEdgeAndVertices();
-}
-
-FlatPattern::FlatPattern(const std::vector<size_t> &numberOfNodesPerSlot,
-                         const std::vector<vcg::Point2i> &edges) {
-  add(numberOfNodesPerSlot, edges);
-  // add normals
-  bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001;
-  if (normalsAreAbsent) {
-    for (auto &v : vert) {
-      v.N() = CoordType(0, 0, 1);
-    }
-  }
-  baseTriangleHeight = (vert[0].cP() - vert[3].cP()).Norm();
-  updateEigenEdgeAndVertices();
-}
-
-bool FlatPattern::createHoneycombAtom() {
-  VCGEdgeMesh honeycombQuarter;
-  const VCGEdgeMesh::CoordType n(0, 0, 1);
-  const double H = 0.2;
-  const double height = 1.5 * H;
-  const double width = 0.2;
-  const double theta = 70;
-  const double dy = tan(vcg::math::ToRad(90 - theta)) * width / 2;
-  vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(
-      honeycombQuarter, VCGEdgeMesh::CoordType(0, height / 2, 0), n);
-  vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(
-      honeycombQuarter, VCGEdgeMesh::CoordType(0, H / 2 - dy, 0), n);
-  vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(
-      honeycombQuarter, VCGEdgeMesh::CoordType(width / 2, H / 2, 0), n);
-  vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(
-      honeycombQuarter, VCGEdgeMesh::CoordType(width / 2, 0, 0), n);
-  vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(honeycombQuarter, 0, 1);
-  vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(honeycombQuarter, 1, 2);
-  vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(honeycombQuarter, 2, 3);
-
-  VCGEdgeMesh honeycombAtom;
-  // Top right
-  vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::MeshCopy(honeycombAtom,
-                                                       honeycombQuarter);
-  // Bottom right
-  vcg::Matrix44d rotM;
-  rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0));
-  vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(honeycombQuarter, rotM);
-  vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(honeycombAtom,
-                                                   honeycombQuarter);
-  // Bottom left
-  rotM.SetRotateDeg(180, vcg::Point3d(0, 1, 0));
-  vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(honeycombQuarter, rotM);
-  vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(honeycombAtom,
-                                                   honeycombQuarter);
-  // Top left
-  rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0));
-  vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(honeycombQuarter, rotM);
-  vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(honeycombAtom,
-                                                   honeycombQuarter);
-
-  for (VertexType &v : honeycombAtom.vert) {
-    v.P()[2] = 0;
-  }
-
-  return true;
-}
-
-void FlatPattern::copy(FlatPattern &copyFrom) {
-  VCGEdgeMesh::copy(copyFrom);
-  baseTriangleHeight = copyFrom.getBaseTriangleHeight();
-}
-
-void FlatPattern::deleteDanglingEdges() {
-  for (VertexType &v : vert) {
-    std::vector<VCGEdgeMesh::EdgePointer> incidentElements;
-    vcg::edge::VEStarVE(&v, incidentElements);
-    if (incidentElements.size() == 1) {
-      vcg::tri::Allocator<VCGEdgeMesh>::DeleteEdge(*this, *incidentElements[0]);
-    }
-    if (incidentElements.size() == 1) {
-      vcg::tri::Allocator<VCGEdgeMesh>::DeleteVertex(*this, v);
-    }
-  }
-
-  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateVertex(*this);
-  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateEdge(*this);
-  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
-}
-
-void FlatPattern::scale(const double &desiredBaseTriangleCentralEdgeSize) {
-  this->baseTriangleHeight = desiredBaseTriangleCentralEdgeSize;
-  const double baseTriangleCentralEdgeSize =
-      (vert[0].cP() - vert[3].cP()).Norm();
-  const double scaleRatio =
-      desiredBaseTriangleCentralEdgeSize / baseTriangleCentralEdgeSize;
-
-  vcg::tri::UpdatePosition<VCGEdgeMesh>::Scale(*this, scaleRatio);
-}
-
-void FlatPattern::deleteDanglingVertices() {
-  vcg::tri::Allocator<FlatPattern>::PointerUpdater<VertexPointer> pu;
-  deleteDanglingVertices(pu);
-}
-
-void FlatPattern::deleteDanglingVertices(
-    vcg::tri::Allocator<FlatPattern>::PointerUpdater<VertexPointer> &pu) {
-  for (VertexType &v : vert) {
-    std::vector<FlatPattern::EdgePointer> incidentElements;
-    vcg::edge::VEStarVE(&v, incidentElements);
-    if (incidentElements.size() == 0) {
-      vcg::tri::Allocator<FlatPattern>::DeleteVertex(*this, v);
-    }
-  }
-
-  vcg::tri::Allocator<FlatPattern>::CompactVertexVector(*this, pu);
-  vcg::tri::Allocator<FlatPattern>::CompactEdgeVector(*this);
-
-  updateEigenEdgeAndVertices();
-}
-
-void FlatPattern::tilePattern(VCGEdgeMesh &pattern, VCGPolyMesh &tileInto,
-                              const bool &shouldDeleteDanglingEdges) {
-  const size_t middleIndex = 3;
-  double xOffset =
-      vcg::Distance(pattern.vert[0].cP(), pattern.vert[middleIndex].cP()) /
-      std::tan(M_PI / 3);
-  CoordType patternCoord0 = pattern.vert[0].cP();
-  CoordType patternBottomRight =
-      pattern.vert[middleIndex].cP() + CoordType(xOffset, 0, 0);
-  CoordType patternBottomLeft =
-      pattern.vert[middleIndex].cP() - CoordType(xOffset, 0, 0);
-  std::vector<vcg::Point3d> patternTrianglePoints{
-      patternCoord0, patternBottomRight, patternBottomLeft};
-  CoordType pointOnPattern =
-      patternCoord0 + (patternBottomLeft - patternCoord0) ^
-      (patternBottomRight - patternCoord0);
-
-  std::vector<CoordType> faceCenters(FN());
-  VCGTriMesh tileIntoEdgeMesh;
-
-  for (VCGPolyMesh::FaceType &f : tileInto.face) {
-    std::vector<VCGPolyMesh::VertexType *> incidentVertices;
-    vcg::face::VFIterator<PFace> vfi(
-        &f, 0); // initialize the iterator to the first face
-                //    vcg::face::Pos p(vfi.F(), f.cV(0));
-                //    vcg::face::VVOrderedStarFF(p, incidentVertices);
-    size_t numberOfNodes = 0;
-    CoordType centerOfFace(0, 0, 0);
-    for (size_t vi = 0; vi < f.VN(); vi++) {
-      numberOfNodes++;
-      centerOfFace = centerOfFace + f.cP(vi);
-    }
-    centerOfFace /= f.VN();
-    vcg::tri::Allocator<VCGTriMesh>::AddVertex(tileIntoEdgeMesh, centerOfFace,
-                                               vcg::Color4b::Yellow);
-
-    //    const size_t vi = vcg::tri::Index<VCGPolyMesh>(tileInto, v);
-    //    std::cout << "vertex " << vi << " has incident vertices:" <<
-    //    std::endl;
-    for (size_t vi = 0; vi < f.VN(); vi++) {
-      //      size_t f = 0;
-      //                std::cout << vcg::tri::Index<VCGTriMesh>(tileInto,
-      //  / incidentVertices[f]) /                          << std::endl;
-
-      // Compute transformation matrix M
-      //        vcg::Matrix44d M;
-      std::vector<vcg::Point3d> meshTrianglePoints{
-          centerOfFace, f.cP(vi), vi + 1 == f.VN() ? f.cP(0) : f.cP(vi + 1)};
-      CoordType faceNormal = ((meshTrianglePoints[1] - meshTrianglePoints[0]) ^
-                              (meshTrianglePoints[2] - meshTrianglePoints[0]))
-                                 .Normalize();
-      auto fit = vcg::tri::Allocator<VCGTriMesh>::AddFace(
-          tileIntoEdgeMesh, meshTrianglePoints[0], meshTrianglePoints[1],
-          meshTrianglePoints[2]);
-      fit->N() = faceNormal;
-      CoordType pointOnTriMesh =
-          meshTrianglePoints[0] +
-              (meshTrianglePoints[1] - meshTrianglePoints[0]) ^
-          (meshTrianglePoints[2] - meshTrianglePoints[0]);
-      vcg::Matrix44d M;
-      //        vcg::ComputeRigidMatchMatrix(meshTrianglePoints,
-      //        patternTrianglePoints,
-      //                                     M);
-      vcg::Matrix44d A_prime;
-      A_prime[0][0] = meshTrianglePoints[0][0];
-      A_prime[1][0] = meshTrianglePoints[0][1];
-      A_prime[2][0] = meshTrianglePoints[0][2];
-      A_prime[3][0] = 1;
-      A_prime[0][1] = meshTrianglePoints[1][0];
-      A_prime[1][1] = meshTrianglePoints[1][1];
-      A_prime[2][1] = meshTrianglePoints[1][2];
-      A_prime[3][1] = 1;
-      A_prime[0][2] = meshTrianglePoints[2][0];
-      A_prime[1][2] = meshTrianglePoints[2][1];
-      A_prime[2][2] = meshTrianglePoints[2][2];
-      A_prime[3][2] = 1;
-      A_prime[0][3] = pointOnTriMesh[0];
-      A_prime[1][3] = pointOnTriMesh[1];
-      A_prime[2][3] = pointOnTriMesh[2];
-      A_prime[3][3] = 1;
-      vcg::Matrix44d A;
-      A[0][0] = patternTrianglePoints[0][0];
-      A[1][0] = patternTrianglePoints[0][1];
-      A[2][0] = patternTrianglePoints[0][2];
-      A[3][0] = 1;
-      A[0][1] = patternTrianglePoints[1][0];
-      A[1][1] = patternTrianglePoints[1][1];
-      A[2][1] = patternTrianglePoints[1][2];
-      A[3][1] = 1;
-      A[0][2] = patternTrianglePoints[2][0];
-      A[1][2] = patternTrianglePoints[2][1];
-      A[2][2] = patternTrianglePoints[2][2];
-      A[3][2] = 1;
-      A[0][3] = pointOnPattern[0];
-      A[1][3] = pointOnPattern[1];
-      A[2][3] = pointOnPattern[2];
-      A[3][3] = 1;
-      M = A_prime * vcg::Inverse(A);
-
-      VCGEdgeMesh transformedPattern;
-      vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::MeshCopy(transformedPattern,
-                                                           pattern);
-      vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(transformedPattern, M);
-      for (VertexType &v : transformedPattern.vert) {
-        v.N() = faceNormal;
-      }
-
-      vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(*this,
-                                                       transformedPattern);
-    }
-  }
-
-  //  vcg::tri::Clean<VCGEdgeMesh>::MergeCloseVertex(*this, 0.0000000001);
-  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateVertex(*this);
-  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateEdge(*this);
-  //  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
-
-  vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
-
-  vcg::tri::Clean<VCGEdgeMesh>::MergeCloseVertex(*this, 0.0000000001);
-  deleteDanglingVertices();
-  deleteDanglingEdges();
-  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateVertex(*this);
-  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateEdge(*this);
-  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
-  updateEigenEdgeAndVertices();
-  savePly("tiledPattern.ply");
-
-  vcg::tri::Clean<VCGTriMesh>::MergeCloseVertex(tileIntoEdgeMesh, 0.0000000001);
-  vcg::tri::Clean<VCGTriMesh>::RemoveDegenerateVertex(tileIntoEdgeMesh);
-  vcg::tri::Clean<VCGTriMesh>::RemoveDegenerateEdge(tileIntoEdgeMesh);
-  tileIntoEdgeMesh.savePly("tileIntoTriMesh.ply");
-}
-
-void FlatPattern::createFan(const size_t &fanSize) {
-  FlatPattern rotatedPattern;
-  vcg::tri::Append<FlatPattern, FlatPattern>::MeshCopy(rotatedPattern, *this);
-  for (int rotationCounter = 1; rotationCounter < fanSize; rotationCounter++) {
-    vcg::Matrix44d R;
-    auto rotationAxis = vcg::Point3d(0, 0, 1);
-    R.SetRotateDeg(360 / fanSize, rotationAxis);
-    vcg::tri::UpdatePosition<FlatPattern>::Matrix(rotatedPattern, R);
-    vcg::tri::Append<FlatPattern, FlatPattern>::Mesh(*this, rotatedPattern);
-  }
-
-  removeDuplicateVertices();
-  //  const double precision = 1e-4;
-  //    for (size_t vi = 0; vi < VN(); vi++) {
-  //      vert[vi].P()[0] = std::round(vert[vi].P()[0] * (1 / precision)) *
-  //      precision; vert[vi].P()[1] = std::round(vert[vi].P()[1] * (1 /
-  //      precision)) * precision; vert[vi].P()[2] = std::round(vert[vi].P()[2]
-  //      * (1 / precision)) * precision;
-  //    }
-
-  updateEigenEdgeAndVertices();
-}
-
-void FlatPattern::removeDuplicateVertices() {
-  vcg::tri::Clean<VCGEdgeMesh>::MergeCloseVertex(*this, 0.0000000001);
-  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateVertex(*this);
-  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateEdge(*this);
-  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
-  vcg::tri::UpdateTopology<FlatPattern>::VertexEdge(*this);
-}
-
-double FlatPattern::getBaseTriangleHeight() const { return baseTriangleHeight; }
-
-void FlatPattern::tilePattern(VCGEdgeMesh &pattern, VCGTriMesh &tileInto) {
-
-  const size_t middleIndex = 3;
-  double xOffset =
-      vcg::Distance(pattern.vert[0].cP(), pattern.vert[middleIndex].cP()) /
-      std::tan(M_PI / 3);
-  CoordType patternCoord0 = pattern.vert[0].cP();
-  CoordType patternBottomRight =
-      pattern.vert[middleIndex].cP() + CoordType(xOffset, 0, 0);
-  CoordType patternBottomLeft =
-      pattern.vert[middleIndex].cP() - CoordType(xOffset, 0, 0);
-  std::vector<vcg::Point3d> patternTrianglePoints{
-      patternCoord0, patternBottomRight, patternBottomLeft};
-  CoordType pointOnPattern =
-      patternCoord0 + (patternBottomLeft - patternCoord0) ^
-      (patternBottomRight - patternCoord0);
-
-  for (VCGTriMesh::VertexType &v : tileInto.vert) {
-    const auto centralNodeColor = vcg::Color4<unsigned char>(64, 64, 64, 255);
-    const bool isCentralNode = v.cC() == centralNodeColor;
-    if (isCentralNode) {
-      std::vector<VCGTriMesh::VertexType *> incidentVertices;
-      vcg::face::VFIterator<VCGTriMeshFace> vfi(
-          &v); // initialize the iterator tohe first face
-      vcg::face::Pos p(vfi.F(), &v);
-      vcg::face::VVOrderedStarFF(p, incidentVertices);
-
-      const size_t vi = vcg::tri::Index<VCGTriMesh>(tileInto, v);
-      std::cout << "vertex " << vi << " has incident vertices:" << std::endl;
-      for (size_t f = 0; f < incidentVertices.size(); f++) {
-        //      size_t f = 0;
-        std::cout << vcg::tri::Index<VCGTriMesh>(tileInto, incidentVertices[f])
-                  << std::endl;
-
-        // Compute transformation matrix M
-        //        vcg::Matrix44d M;
-        std::vector<vcg::Point3d> meshTrianglePoints{
-            v.cP(), incidentVertices[f]->cP(),
-            f + 1 == incidentVertices.size() ? incidentVertices[0]->cP()
-                                             : incidentVertices[f + 1]->cP()};
-        CoordType faceNormal =
-            ((meshTrianglePoints[1] - meshTrianglePoints[0]) ^
-             (meshTrianglePoints[2] - meshTrianglePoints[0]))
-                .Normalize();
-        CoordType pointOnTriMesh =
-            meshTrianglePoints[0] +
-                (meshTrianglePoints[1] - meshTrianglePoints[0]) ^
-            (meshTrianglePoints[2] - meshTrianglePoints[0]);
-        vcg::Matrix44d M;
-        //        vcg::ComputeRigidMatchMatrix(meshTrianglePoints,
-        //        patternTrianglePoints,
-        //                                     M);
-        vcg::Matrix44d A_prime;
-        A_prime[0][0] = meshTrianglePoints[0][0];
-        A_prime[1][0] = meshTrianglePoints[0][1];
-        A_prime[2][0] = meshTrianglePoints[0][2];
-        A_prime[3][0] = 1;
-        A_prime[0][1] = meshTrianglePoints[1][0];
-        A_prime[1][1] = meshTrianglePoints[1][1];
-        A_prime[2][1] = meshTrianglePoints[1][2];
-        A_prime[3][1] = 1;
-        A_prime[0][2] = meshTrianglePoints[2][0];
-        A_prime[1][2] = meshTrianglePoints[2][1];
-        A_prime[2][2] = meshTrianglePoints[2][2];
-        A_prime[3][2] = 1;
-        A_prime[0][3] = pointOnTriMesh[0];
-        A_prime[1][3] = pointOnTriMesh[1];
-        A_prime[2][3] = pointOnTriMesh[2];
-        A_prime[3][3] = 1;
-        vcg::Matrix44d A;
-        A[0][0] = patternTrianglePoints[0][0];
-        A[1][0] = patternTrianglePoints[0][1];
-        A[2][0] = patternTrianglePoints[0][2];
-        A[3][0] = 1;
-        A[0][1] = patternTrianglePoints[1][0];
-        A[1][1] = patternTrianglePoints[1][1];
-        A[2][1] = patternTrianglePoints[1][2];
-        A[3][1] = 1;
-        A[0][2] = patternTrianglePoints[2][0];
-        A[1][2] = patternTrianglePoints[2][1];
-        A[2][2] = patternTrianglePoints[2][2];
-        A[3][2] = 1;
-        A[0][3] = pointOnPattern[0];
-        A[1][3] = pointOnPattern[1];
-        A[2][3] = pointOnPattern[2];
-        A[3][3] = 1;
-        M = A_prime * vcg::Inverse(A);
-
-        VCGEdgeMesh transformedPattern;
-        vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::MeshCopy(transformedPattern,
-                                                             pattern);
-        vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(transformedPattern, M);
-        for (VertexType &v : transformedPattern.vert) {
-          v.N() = faceNormal;
-        }
-
-        vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(*this,
-                                                         transformedPattern);
-      }
-    }
-  }
-
-  vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
-  deleteDanglingVertices();
-  deleteDanglingEdges();
-  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
-
-  updateEigenEdgeAndVertices();
-}

flatpattern.hpp

@@ -1,35 +0,0 @@
-#ifndef FLATPATTERN_HPP
-#define FLATPATTERN_HPP
-
-#include "trianglepatterngeometry.hpp"
-
-class FlatPattern : public FlatPatternGeometry {
-public:
-  FlatPattern();
-  FlatPattern(const std::string &filename, bool addNormalsIfAbsent = true);
-  FlatPattern(const std::vector<size_t> &numberOfNodesPerSlot,
-              const std::vector<vcg::Point2i> &edges);
-
-  bool createHoneycombAtom();
-  void copy(FlatPattern &copyFrom);
-
-  void tilePattern(VCGEdgeMesh &pattern, VCGTriMesh &tileInto);
-  void tilePattern(VCGEdgeMesh &pattern, VCGPolyMesh &tileInto,
-                   const bool &shouldDeleteDanglingEdges);
-
-  void createFan(const size_t &fanSize = 6);
-
-  void deleteDanglingVertices(
-      vcg::tri::Allocator<FlatPattern>::PointerUpdater<VertexPointer> &pu);
-  void deleteDanglingVertices();
-  void scale(const double &desiredBaseTriangleCentralEdgeSize);
-
-  double getBaseTriangleHeight() const;
-
-private:
-  void deleteDanglingEdges();
-  void removeDuplicateVertices();
-  double baseTriangleHeight;
-};
-
-#endif // FLATPATTERN_HPP

mesh.hpp

@@ -0,0 +1,20 @@
+#ifndef MESH_HPP
+#define MESH_HPP
+
+#include <string>
+
+class Mesh {
+  std::string label;
+
+public:
+  virtual ~Mesh() = default;
+  virtual bool load(const std::string &filePath) { return false; }
+  std::string getLabel() const;
+  void setLabel(const std::string &newLabel);
+};
+
+inline std::string Mesh::getLabel() const { return label; }
+
+inline void Mesh::setLabel(const std::string &newLabel) { label = newLabel; }
+
+#endif // MESH_HPP

patternIO.cpp

@@ -1,116 +1,141 @@
 #include "patternIO.hpp"
+#include "trianglepatterngeometry.hpp"
 #include <filesystem>
 #include <fstream>
 #include <iostream>
 #include <sstream>
 
-PatternIO::PatternIO() {}
+void PatternIO::write(std::ofstream &fileStream, const Pattern &pattern) {
+  fileStream << "p " << pattern.name << " " << pattern.edges.size();
+  for (const vcg::Point2i &edge : pattern.edges) {
+    fileStream << " " << edge[0] << " " << edge[1];
+  }
+}
+void PatternIO::write(std::ofstream &fileStream,
+                      const std::vector<vcg::Point3d> &vertices) {
+  fileStream << "#Nodes" << '\n';
+  for (vcg::Point3d node : vertices) {
+    fileStream << "n " << node.X() << " " << node.Y() << " " << node.Z()
+               << '\n';
+  }
+  fileStream << "#Patterns" << '\n';
+}
 
-void PatternIO::save(const std::string &filepath,
+void PatternIO::exportPLY(const std::string &inputFilePath,
+                          const int &startIndex, const int &endIndex) {
+  assert(std::filesystem::path(inputFilePath).extension() == ".patt");
+  const std::string outputDir =
+      std::filesystem::path(inputFilePath)
+          .parent_path()
+          .append("PLYFiles")
+          .append(std::to_string(startIndex) + "_" + std::to_string(endIndex))
+          .string();
+  exportPLY(inputFilePath, outputDir, startIndex, endIndex);
+}
+
+void PatternIO::exportPLY(const std::string &patternSetFilePath,
+                          const std::string &outputDirectoryPath,
+                          const int &startIndex, const int &endIndex) {
+  assert(std::filesystem::path(patternSetFilePath).extension() == ".patt");
+  std::vector<vcg::Point3d> vertices;
+  std::vector<Pattern> patterns;
+  load(patternSetFilePath, vertices, patterns, startIndex, endIndex);
+  std::filesystem::create_directories(outputDirectoryPath);
+
+  for (int patternIndex = 0; patternIndex < patterns.size(); patternIndex++) {
+    PatternGeometry p;
+    p.add(vertices, patterns[patternIndex].edges);
+    auto tiled = p.createTile(p);
+    p.savePly(
+        std::filesystem::path(outputDirectoryPath)
+            .append(std::to_string(patterns[patternIndex].name) + ".ply"));
+    tiled.savePly(std::filesystem::path(outputDirectoryPath)
+                      .append("tiled_" +
+                              std::to_string(patterns[patternIndex].name) +
+                              ".ply"));
+  }
+}
+
+void PatternIO::save(const std::string &filePath,
                      const PatternSet &patternSet) {
   std::ofstream fileStream;
-  if (std::filesystem::exists(filepath)) {
+  if (!std::filesystem::exists(filePath)) {
-    fileStream.open(filepath, std::ios_base::app);
+    fileStream.open(filePath);
+    write(fileStream, patternSet.nodes);
   } else {
-    fileStream.open(filepath);
+    fileStream.open(filePath, std::ios_base::app);
-    fileStream << "#Nodes" << std::endl;
-    for (vcg::Point2d node : patternSet.nodes) {
-      fileStream << "n " << node.X() << " " << node.Y() << std::endl;
-    }
-    fileStream << "#Patterns" << std::endl;
   }
 
   for (const Pattern &pattern : patternSet.patterns) {
-    fileStream << "p " << pattern.labels.size() << " " << pattern.edges.size()
+    write(fileStream, pattern);
-               << "  ";
+    fileStream << '\n';
-    for (size_t labelIndex = 0; labelIndex < pattern.labels.size() - 1;
-         labelIndex++) {
-      fileStream << pattern.labels[labelIndex] << " ";
-    }
-    fileStream << pattern.labels[pattern.labels.size() - 1] << " ";
-
-    for (const vcg::Point2i &edge : pattern.edges) {
-      fileStream << " " << edge[0] << " " << edge[1] << " ";
-    }
-    fileStream << std::endl;
   }
+
+  fileStream.close();
 }
 
-void PatternIO::save(const std::string &filepath, const Pattern &pattern) {
+void PatternIO::save(std::ofstream &fileStream, const Pattern &pattern) {
+  assert(fileStream.is_open());
+  write(fileStream, pattern);
+  fileStream << '\n';
+}
+
+void PatternIO::save(const std::string &filePath, const Pattern &pattern) {
+  if (!std::filesystem::exists(filePath)) {
+    std::cerr << "File must already exist. The node information must already "
+                 "be in the file."
+              << '\n';
+    std::terminate();
+  }
   std::ofstream fileStream;
-  if (std::filesystem::exists(filepath)) {
+  fileStream.open(filePath, std::ios_base::app);
-    fileStream.open(filepath, std::ios_base::app);
-  } else {
-    fileStream.open(filepath);
-    fileStream << "#Nodes" << std::endl;
-    fileStream << "#Patterns" << std::endl;
-  }
 
-  fileStream << "p " << pattern.labels.size() << " " << pattern.edges.size()
+  write(fileStream, pattern);
-             << "  ";
+  fileStream << '\n';
-  for (size_t labelIndex = 0; labelIndex < pattern.labels.size() - 1;
-       labelIndex++) {
-    fileStream << pattern.labels[labelIndex] << " ";
-  }
-  fileStream << pattern.labels[pattern.labels.size() - 1] << " ";
-
-  for (const vcg::Point2i &edge : pattern.edges) {
-    fileStream << " " << edge[0] << " " << edge[1] << " ";
-  }
-  fileStream << std::endl;
 }
 
-void PatternIO::load(const std::string &filepath, PatternSet &patternSet,
+void PatternIO::load(const std::string &filepath,
-                     const std::vector<PatternLabel> &targetLabels) {
+                     std::vector<vcg::Point3d> &vertices,
+                     std::vector<Pattern> &patterns, const int &startIndex,
+                     const int &endIndex) {
   std::ifstream fileStream(filepath);
   std::string line;
-  std::vector<PatternLabel> sortedTargetPatternLabels(targetLabels);
+  int patternIndex = 0;
-  std::sort(sortedTargetPatternLabels.begin(), sortedTargetPatternLabels.end());
+  while (std::getline(fileStream, line) && patternIndex <= endIndex) {
-  while (std::getline(fileStream, line)) {
+    //    std::cout << line << '\n';
-    std::cout << line << std::endl;
     std::istringstream iss(line);
     char lineID;
     iss >> lineID;
     if (lineID == 'n') {
-      double x, y;
+      double x, y, z = 0;
       iss >> x >> y;
-      std::cout << x << " " << y << std::endl;
+      //      std::cout << x << " " << y << z << '\n';
+      vertices.push_back(vcg::Point3d(x, y, z));
     } else if (lineID == 'p') {
+      patternIndex++;
+      if (startIndex > patternIndex) {
+        continue;
+      }
       Pattern pattern;
-      size_t numberOfLabels, numberOfEdges;
+      iss >> pattern.name;
-      iss >> numberOfLabels >> numberOfEdges;
+      int numberOfEdges;
-      std::cout << "NumberOfLabels:" << numberOfLabels << std::endl;
+      iss >> numberOfEdges;
-      std::cout << "NumberOfEdges:" << numberOfEdges << std::endl;
+      for (int edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
-      std::vector<PatternLabel> patternLabels;
-      for (size_t labelIndex = 0; labelIndex < numberOfLabels; labelIndex++) {
-        size_t patternLabel;
-        iss >> patternLabel;
-        PatternLabel pl = static_cast<PatternLabel>(patternLabel);
-        std::cout << "Pattern label read:" << patternLabel << std::endl;
-        patternLabels.push_back(pl);
-      }
-      if (!targetLabels.empty()) {
-        std::sort(patternLabels.begin(), patternLabels.end());
-        std::vector<PatternLabel> labelIntersection;
-        std::set_intersection(patternLabels.begin(), patternLabels.end(),
-                              sortedTargetPatternLabels.begin(),
-                              sortedTargetPatternLabels.end(),
-                              labelIntersection.begin());
-        if (!labelIntersection.empty()) {
-          pattern.labels = patternLabels;
-        } else {
-          continue;
-        }
-      } else {
-        pattern.labels = patternLabels;
-      }
-      for (size_t edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
         size_t ni0, ni1;
         iss >> ni0 >> ni1;
         vcg::Point2i edge(ni0, ni1);
         pattern.edges.push_back(edge);
-        std::cout << "Node indices read:" << ni0 << "," << ni1 << std::endl;
       }
+      patterns.push_back(pattern);
+
+      //      PatternGeometry fpg;
+      //      fpg.add(vertices, pattern.edges);
+      //      fpg.registerForDrawing(glm::vec3(1, 0, 0));
+      //      PatternGeometry tiledPattern = fpg.createTile(fpg);
+      //      tiledPattern.deleteDanglingVertices();
+      //      tiledPattern.registerForDrawing();
+      //      polyscope::show();
+      //      polyscope::removeAllStructures();
     }
   }
 }

patternIO.hpp

@@ -3,38 +3,62 @@
 
 #include <fstream>
 #include <string>
+#include <unordered_set>
 #include <vcg/space/deprecated_point2.h>
+#include <vcg/space/deprecated_point3.h>
 #include <vector>
 
-enum PatternLabel {
+// enum PatternLabel {
-  Valid = 0,
+//  Valid = 0,
-  MultipleCC,
+//  MultipleCC,
-  DanglingEdge,
+//  DanglingEdge,
-  IntersectingEdges,
+//  IntersectingEdges,
-  ArticulationPoints
+//  ArticulationPoints
-};
+//};
-
+namespace PatternIO {
+using PatternName = int;
 struct Pattern {
   std::vector<vcg::Point2i> edges;
-  std::vector<PatternLabel> labels;
+  PatternName name{-1};
+
+  bool operator==(const Pattern &p) const { return name == p.name; }
+};
+struct PatternHashFunction {
+  std::size_t operator()(const Pattern &p) const {
+    return std::hash<int>()(p.name);
+  }
 };
 
 /*
  * A set of planar patterns using the same nodes
  * */
 struct PatternSet {
-  std::vector<vcg::Point2d> nodes;
+  std::vector<vcg::Point3d> nodes;
+
   std::vector<Pattern> patterns;
 };
 
-class PatternIO {
+void save(const std::string &filePath, const Pattern &pattern);
+void save(const std::string &filePath, const PatternSet &patterns);
+void load(const std::string &filePath, PatternSet &patternSet);
+void save(const std::string &filePath,
+          const std::vector<vcg::Point3d> &vertices);
+void write(std::ofstream &fileStream, const Pattern &pattern);
+void write(std::ofstream &fileStream,
+           const std::vector<vcg::Point3d> &vertices);
+void exportPLY(const std::string &patternSetFilePath,
+               const std::string &outputDirectoryPath, const int &startIndex,
+               const int &endIndex);
+void exportPLY(const std::string &outputDirectoryPath,
+               const std::string &patternSetFilePath,
+               const PatternName &patternToExport);
+void save(std::ofstream &fileStream, const Pattern &pattern);
+void load(const std::string &filepath, std::vector<vcg::Point3d> &vertices,
+          std::vector<Pattern> &patterns, const int &startIndex,
+          const int &endIndex);
+void exportPLY(const std::string &inputFilePath, const int &startIndex,
+               const int &indexStep);
 
-public:
+} // namespace PatternIO
-  PatternIO();
-  static void save(const std::string &filepath, const Pattern &pattern);
-  static void save(const std::string &filepath, const PatternSet &patterns);
-  static void load(const std::string &filepath, PatternSet &patternSet,
-                   const std::vector<PatternLabel> &targetLabels = {});
-};
 
 #endif // PATTERNEXPORTER_HPP

polymesh.hpp

@@ -1,9 +1,9 @@
 #ifndef POLYMESH_HPP
 #define POLYMESH_HPP
+#include "mesh.hpp"
 #include "vcg/complex/complex.h"
 #include <filesystem>
 #include <wrap/io_trimesh/import.h>
-//#include <wrap/nanoply/include/nanoplyWrapper.hpp>
 
 class PFace;
 class PVertex;
@@ -37,13 +37,18 @@ class PFace
           > {};
 
 class VCGPolyMesh
-    : public vcg::tri::TriMesh<std::vector<PVertex>, // the vector of vertices
+    : public vcg::tri::TriMesh<std::vector<PVertex>, std::vector<PFace>>,
-                               std::vector<PFace>    // the vector of faces
+      public Mesh {
-                               > {
 public:
-  void loadFromPlyFile(const std::string &filename) {
+  virtual bool load(const std::string &filename) override {
-    vcg::tri::io::ImporterOBJ<VCGPolyMesh>::Info info;
+    int mask;
-    vcg::tri::io::ImporterOBJ<VCGPolyMesh>::Open(*this, filename.c_str(), info);
+    vcg::tri::io::Importer<VCGPolyMesh>::LoadMask(filename.c_str(), mask);
+    int error = vcg::tri::io::Importer<VCGPolyMesh>::Open(
+        *this, filename.c_str(), mask);
+    if (error != 0) {
+      return false;
+    }
+    //    vcg::tri::io::ImporterOBJ<VCGPolyMesh>::Open();
     //    unsigned int mask = 0;
     //    mask |= nanoply::NanoPlyWrapper<VCGPolyMesh>::IO_VERTCOORD;
     //    mask |= nanoply::NanoPlyWrapper<VCGPolyMesh>::IO_VERTNORMAL;
@@ -58,7 +63,13 @@ public:
     vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
     //    vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
     vcg::tri::UpdateNormal<VCGPolyMesh>::PerVertexNormalized(*this);
+    return true;
   }
+#ifdef POLYSCOPE_DEFINED
+  void registerForDrawing(){
+
+  }
+#endif
 };
 
 #endif // POLYMESH_HPP

simulationresult.hpp

@@ -109,6 +109,7 @@ public:
     return json.dump();
   }
   bool load(const std::string &jsonFilename) {
+    const bool beVerbose = false;
     if (std::filesystem::path(jsonFilename).extension() != ".json") {
       std::cerr << "A json file is expected as input. The given file has the "
                    "following extension:"
@@ -124,7 +125,9 @@ public:
       return false;
     }
 
-    std::cout << "Loading json file:" << jsonFilename << std::endl;
+    if (beVerbose) {
+      std::cout << "Loading json file:" << jsonFilename << std::endl;
+    }
     nlohmann::json json;
     std::ifstream ifs(jsonFilename);
     ifs >> json;
@@ -136,7 +139,7 @@ public:
           std::filesystem::path(
               std::filesystem::path(jsonFilename).parent_path())
               .append(relativeFilepath);
-      pMesh->loadPly(meshFilepath.string());
+      pMesh->load(meshFilepath.string());
       pMesh->setLabel(
           json[jsonLabels.meshLabel]); // FIXME: This should be exported using
                                        // nanoply but nanoply might not be able
@@ -148,9 +151,11 @@ public:
           //      auto conV =
           json[jsonLabels.constrainedVertices]
               .get<std::unordered_map<VertexIndex, std::unordered_set<int>>>();
-      std::cout << "Loaded constrained vertices. Number of constrained "
+      if (beVerbose) {
-                   "vertices found:"
+        std::cout << "Loaded constrained vertices. Number of constrained "
-                << constrainedVertices.size() << std::endl;
+                     "vertices found:"
+                  << constrainedVertices.size() << std::endl;
+      }
     }
 
     if (json.contains(jsonLabels.nodalForces)) {
@@ -159,8 +164,10 @@ public:
       for (const auto &forces : f) {
         nodalExternalForces[forces.first] = Vector6d(forces.second);
       }
-      std::cout << "Loaded forces. Number of forces found:"
+      if (beVerbose) {
-                << nodalExternalForces.size() << std::endl;
+        std::cout << "Loaded forces. Number of forces found:"
+                  << nodalExternalForces.size() << std::endl;
+      }
     }
 
     if (json.contains(jsonLabels.label)) {
@@ -389,7 +396,8 @@ struct SimulationResults {
     polyscope::removeCurveNetwork(getLabel());
   }
   void registerForDrawing(
-      const std::optional<glm::vec3> &desiredColor = std::nullopt) const {
+      const std::optional<glm::vec3> &desiredColor = std::nullopt,
+      const bool &shouldEnable = true) const {
     polyscope::options::groundPlaneEnabled = false;
     polyscope::view::upDir = polyscope::view::UpDir::ZUp;
     const std::string branchName = "Branch:Polyscope";
@@ -405,10 +413,10 @@ struct SimulationResults {
     //        mesh->getEigenEdges());
 
     const std::shared_ptr<SimulationMesh> &mesh = job->pMesh;
-    auto polyscopeHandle_deformedEdmeMesh =
+    auto polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(
-        polyscope::registerCurveNetwork(getLabel(), mesh->getEigenVertices(),
+        getLabel(), mesh->getEigenVertices(), mesh->getEigenEdges());
-                                        mesh->getEigenEdges())
+    polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable);
-            ->setRadius(0.002, true);
+    polyscopeHandle_deformedEdmeMesh->setRadius(0.002, true);
     if (desiredColor.has_value()) {
       polyscopeHandle_deformedEdmeMesh->setColor(desiredColor.value());
     }

topologyenumerator.cpp

@@ -1,609 +0,0 @@
-#include "topologyenumerator.hpp"
-#include <algorithm>
-#include <iostream>
-#include <math.h>
-#include <numeric>
-#include <unordered_set>
-
-const bool debugIsOn{false};
-const bool exportArticulationPointsPatterns{false};
-const bool savePlyFiles{true};
-
-// size_t binomialCoefficient(size_t n, size_t m) {
-//  assert(n > m);
-//  return tgamma(n + 1) / (tgamma(m + 1) * tgamma(n - m + 1));
-//}
-
-// void TopologyEnumerator::createLabelMesh(
-//    const std::vector<vcg::Point3d> vertices,
-//    const std::filesystem::path &savePath) const {
-//  const std::string allOnes(patternTopology.getNumberOfPossibleEdges(), '1');
-//  const std::vector<vcg::Point2i> allEdges =
-//      TrianglePatternTopology::convertToEdges(allOnes, vertices.size());
-//  TrianglePatternGeometry labelMesh;
-//  std::vector<vcg::Point3d> labelVertices(allEdges.size());
-//  for (size_t edgeIndex = 0; edgeIndex < allEdges.size(); edgeIndex++) {
-//    const vcg::Point3d edgeMidpoint =
-//        (vertices[allEdges[edgeIndex][0]] + vertices[allEdges[edgeIndex][1]])
-//        / 2;
-//    labelVertices[edgeIndex] = edgeMidpoint;
-//  }
-//  labelMesh.set(labelVertices);
-//  labelMesh.savePly(std::filesystem::path(savePath)
-//                        .append(std::string("labelMesh.ply"))
-//                        .string());
-//}
-
-size_t TopologyEnumerator::getEdgeIndex(size_t ni0, size_t ni1) const {
-  if (ni1 <= ni0) {
-    std::swap(ni0, ni1);
-  }
-  assert(ni1 > ni0);
-  const size_t &n = numberOfNodes;
-  return (n * (n - 1) / 2) - (n - ni0) * ((n - ni0) - 1) / 2 + ni1 - ni0 - 1;
-}
-
-TopologyEnumerator::TopologyEnumerator() {}
-
-void TopologyEnumerator::computeValidPatterns(
-    const std::vector<size_t> &reducedNumberOfNodesPerSlot) {
-  assert(reducedNumberOfNodesPerSlot.size() == 5);
-  assert(reducedNumberOfNodesPerSlot[0] == 0 ||
-         reducedNumberOfNodesPerSlot[0] == 1);
-  assert(reducedNumberOfNodesPerSlot[1] == 0 ||
-         reducedNumberOfNodesPerSlot[1] == 1);
-  std::vector<size_t> numberOfNodesPerSlot{
-      reducedNumberOfNodesPerSlot[0], reducedNumberOfNodesPerSlot[1],
-      reducedNumberOfNodesPerSlot[1], reducedNumberOfNodesPerSlot[2],
-      reducedNumberOfNodesPerSlot[3], reducedNumberOfNodesPerSlot[2],
-      reducedNumberOfNodesPerSlot[4]};
-  // Generate an edge mesh wih all possible edges
-  numberOfNodes = std::accumulate(numberOfNodesPerSlot.begin(),
-                                  numberOfNodesPerSlot.end(), 0);
-  const size_t numberOfAllPossibleEdges =
-      numberOfNodes * (numberOfNodes - 1) / 2;
-
-  std::vector<vcg::Point2i> allPossibleEdges(numberOfAllPossibleEdges);
-  const int &n = numberOfNodes;
-  for (size_t edgeIndex = 0; edgeIndex < numberOfAllPossibleEdges;
-       edgeIndex++) {
-    const int ni0 =
-        n - 2 -
-        std::floor(std::sqrt(-8 * edgeIndex + 4 * n * (n - 1) - 7) / 2.0 - 0.5);
-    const int ni1 =
-        edgeIndex + ni0 + 1 - n * (n - 1) / 2 + (n - ni0) * ((n - ni0) - 1) / 2;
-    allPossibleEdges[edgeIndex] = vcg::Point2i(ni0, ni1);
-  }
-  FlatPatternGeometry patternGeometryAllEdges;
-  patternGeometryAllEdges.add(numberOfNodesPerSlot, allPossibleEdges);
-  // Create Results path
-  auto resultPath =
-      //      std::filesystem::path("/home/iason/Documents/PhD/Research/Enumerating\\
-      //      "
-      //                            "2d\\ connections\\ of\\ nodes");
-      std::filesystem::current_path()
-          .parent_path()
-          .parent_path()
-          .parent_path()
-          .parent_path();
-  assert(std::filesystem::exists(resultPath));
-
-  auto allResultsPath = resultPath.append("Results");
-  std::filesystem::create_directory(allResultsPath);
-  std::string setupString;
-  //  for (size_t numberOfNodes : reducedNumberOfNodesPerSlot) {
-  for (size_t numberOfNodesPerSlotIndex = 0;
-       numberOfNodesPerSlotIndex < reducedNumberOfNodesPerSlot.size();
-       numberOfNodesPerSlotIndex++) {
-    std::string elemID;
-    if (numberOfNodesPerSlotIndex == 0 || numberOfNodesPerSlotIndex == 1) {
-      elemID = "v";
-    } else if (numberOfNodesPerSlotIndex == 2 ||
-               numberOfNodesPerSlotIndex == 3) {
-      elemID = "e";
-    } else {
-      elemID = "c";
-    }
-    setupString +=
-        std::to_string(reducedNumberOfNodesPerSlot[numberOfNodesPerSlotIndex]) +
-        elemID + "_";
-  }
-  setupString += std::to_string(FlatPatternGeometry().getFanSize()) + "fan";
-  if (debugIsOn) {
-    setupString += "_debug";
-  }
-  auto resultsPath = std::filesystem::path(allResultsPath).append(setupString);
-  //    std::filesystem::remove_all(resultsPath); // delete previous results
-  std::filesystem::create_directory(resultsPath);
-  if (debugIsOn) {
-    patternGeometryAllEdges.savePly(std::filesystem::path(resultsPath)
-                                        .append("allPossibleEdges.ply")
-                                        .string());
-  }
-  //  statistics.numberOfPossibleEdges = numberOfAllPossibleEdges;
-
-  std::vector<vcg::Point2i> validEdges =
-      getValidEdges(numberOfNodesPerSlot, resultsPath, patternGeometryAllEdges,
-                    allPossibleEdges);
-  FlatPatternGeometry patternAllValidEdges;
-  patternAllValidEdges.add(patternGeometryAllEdges.getVertices(), validEdges);
-  if (debugIsOn) {
-    // Export all valid edges in a ply
-    patternAllValidEdges.savePly(
-        std::filesystem::path(resultsPath).append("allValidEdges.ply").string());
-  }
-  //  statistics.numberOfValidEdges = validEdges.size();
-
-  // Find pairs of intersecting edges
-  std::unordered_map<size_t, std::unordered_set<size_t>> intersectingEdges =
-      patternAllValidEdges.getIntersectingEdges(
-          statistics.numberOfIntersectingEdgePairs);
-  if (debugIsOn) {
-    auto intersectingEdgesPath = std::filesystem::path(resultsPath)
-                                     .append("All_intersecting_edge_pairs");
-    std::filesystem::create_directory(intersectingEdgesPath);
-    // Export intersecting pairs in ply files
-    for (auto mapIt = intersectingEdges.begin();
-         mapIt != intersectingEdges.end(); mapIt++) {
-      for (auto setIt = mapIt->second.begin(); setIt != mapIt->second.end();
-           setIt++) {
-        FlatPatternGeometry intersectingEdgePair;
-        const size_t ei0 = mapIt->first;
-        const size_t ei1 = *setIt;
-        vcg::tri::Allocator<FlatPatternGeometry>::AddEdge(
-            intersectingEdgePair,
-            patternGeometryAllEdges.getVertices()[validEdges[ei0][0]],
-            patternGeometryAllEdges.getVertices()[validEdges[ei0][1]]);
-        vcg::tri::Allocator<FlatPatternGeometry>::AddEdge(
-            intersectingEdgePair,
-            patternGeometryAllEdges.getVertices()[validEdges[ei1][0]],
-            patternGeometryAllEdges.getVertices()[validEdges[ei1][1]]);
-        intersectingEdgePair.savePly(
-            std::filesystem::path(intersectingEdgesPath)
-                .append(std::to_string(mapIt->first) + "_" +
-                        std::to_string(*setIt) + ".ply")
-                .string());
-      }
-    }
-  }
-
-  //  assert(validEdges.size() == allPossibleEdges.size() -
-  //  coincideEdges.size() -
-  //                                   duplicateEdges.size());
-
-  PatternSet patternSet;
-  const std::vector<vcg::Point3d> nodes = patternGeometryAllEdges.getVertices();
-  const size_t numberOfNodes = nodes.size();
-  patternSet.nodes.resize(numberOfNodes);
-  for (size_t nodeIndex = 0; nodeIndex < numberOfNodes; nodeIndex++) {
-    patternSet.nodes[nodeIndex] =
-        vcg::Point2d(nodes[nodeIndex][0], nodes[nodeIndex][1]);
-  }
-  if (std::filesystem::exists(std::filesystem::path(resultsPath)
-                                  .append("patterns.patt")
-                                  .string())) {
-    std::filesystem::remove(
-        std::filesystem::path(resultsPath).append("patterns.patt"));
-  }
-  for (size_t numberOfEdges = 2; numberOfEdges < validEdges.size();
-       numberOfEdges++) {
-    //  for (size_t numberOfEdges = 1; numberOfEdges < 3; numberOfEdges++) {
-    std::cout << "Computing " + setupString << " with " << numberOfEdges
-              << " edges." << std::endl;
-    auto perEdgeResultPath = std::filesystem::path(resultsPath)
-                                 .append(std::to_string(numberOfEdges));
-    //        if (std::filesystem::exists(perEdgeResultPath)) {
-    //          continue;
-    //        }
-    std::filesystem::create_directory(perEdgeResultPath);
-    computeValidPatterns(numberOfNodesPerSlot, numberOfEdges, perEdgeResultPath,
-                         patternGeometryAllEdges.getVertices(),
-                         intersectingEdges, validEdges, patternSet);
-    //    statistics.print(setupString, perEdgeResultPath);
-    PatternIO::save(
-        std::filesystem::path(resultsPath).append("patterns.patt").string(),
-        patternSet);
-  }
-}
-
-void TopologyEnumerator::computeEdgeNodes(
-    const std::vector<size_t> &numberOfNodesPerSlot,
-    std::vector<size_t> &nodesEdge0, std::vector<size_t> &nodesEdge1,
-    std::vector<size_t> &nodesEdge2) {
-  // Create vectors holding the node indices of each pattern node of each
-  // triangle edge
-  size_t nodeIndex = 0;
-  if (numberOfNodesPerSlot[0] != 0) {
-    nodesEdge0.push_back(nodeIndex++);
-  }
-  if (numberOfNodesPerSlot[1] != 0)
-    nodesEdge1.push_back(nodeIndex++);
-  if (numberOfNodesPerSlot[2] != 0)
-    nodesEdge2.push_back(nodeIndex++);
-
-  if (numberOfNodesPerSlot[3] != 0) {
-    for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[3];
-         edgeNodeIndex++) {
-      nodesEdge0.push_back(nodeIndex++);
-    }
-  }
-  if (numberOfNodesPerSlot[4] != 0) {
-    for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[4];
-         edgeNodeIndex++) {
-      nodesEdge1.push_back(nodeIndex++);
-    }
-  }
-
-  if (numberOfNodesPerSlot[5] != 0) {
-    for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[5];
-         edgeNodeIndex++) {
-      nodesEdge2.push_back(nodeIndex++);
-    }
-  }
-  if (numberOfNodesPerSlot[1] != 0) {
-    assert(numberOfNodesPerSlot[2]);
-    nodesEdge0.push_back(1);
-    nodesEdge1.push_back(2);
-  }
-
-  if (numberOfNodesPerSlot[0] != 0) {
-    nodesEdge2.push_back(0);
-  }
-}
-
-std::unordered_set<size_t> TopologyEnumerator::computeCoincideEdges(
-    const std::vector<size_t> &numberOfNodesPerSlot) {
-  /*
-   * A coincide edge is defined as an edge connection between two nodes that lay
-   * on a triangle edge and which have another node in between
-   * */
-  std::vector<size_t> nodesEdge0; // left edge
-  std::vector<size_t> nodesEdge1; // bottom edge
-  std::vector<size_t> nodesEdge2; // right edge
-  computeEdgeNodes(numberOfNodesPerSlot, nodesEdge0, nodesEdge1, nodesEdge2);
-
-  std::vector<size_t> coincideEdges0 = getCoincideEdges(nodesEdge0);
-  std::vector<size_t> coincideEdges1 = getCoincideEdges(nodesEdge1);
-  std::vector<size_t> coincideEdges2 = getCoincideEdges(nodesEdge2);
-  std::unordered_set<size_t> coincideEdges{coincideEdges0.begin(),
-                                           coincideEdges0.end()};
-  std::copy(coincideEdges1.begin(), coincideEdges1.end(),
-            std::inserter(coincideEdges, coincideEdges.end()));
-  std::copy(coincideEdges2.begin(), coincideEdges2.end(),
-            std::inserter(coincideEdges, coincideEdges.end()));
-
-  if (numberOfNodesPerSlot[0] && numberOfNodesPerSlot[1]) {
-    coincideEdges.insert(getEdgeIndex(0, 2));
-  }
-
-  if (numberOfNodesPerSlot[0] && numberOfNodesPerSlot[2]) {
-    assert(numberOfNodesPerSlot[1]);
-    coincideEdges.insert(getEdgeIndex(0, 3));
-  }
-
-  return coincideEdges;
-}
-
-std::unordered_set<size_t> TopologyEnumerator::computeDuplicateEdges(
-    const std::vector<size_t> &numberOfNodesPerSlot) {
-  /*
-   * A duplicate edges are all edges the "right" edge since due to rotational
-   * symmetry "left" edge=="right" edge
-   * */
-  std::unordered_set<size_t> duplicateEdges;
-  std::vector<size_t> nodesEdge0; // left edge
-  std::vector<size_t> nodesEdge1; // bottom edge
-  std::vector<size_t> nodesEdge2; // right edge
-  computeEdgeNodes(numberOfNodesPerSlot, nodesEdge0, nodesEdge1, nodesEdge2);
-  if (numberOfNodesPerSlot[5]) {
-    for (size_t edge2NodeIndex = 0; edge2NodeIndex < nodesEdge2.size() - 1;
-         edge2NodeIndex++) {
-      const size_t nodeIndex = nodesEdge2[edge2NodeIndex];
-      const size_t nextNodeIndex = nodesEdge2[edge2NodeIndex + 1];
-      duplicateEdges.insert(getEdgeIndex(nodeIndex, nextNodeIndex));
-    }
-  }
-
-  return duplicateEdges;
-}
-
-std::vector<vcg::Point2i> TopologyEnumerator::getValidEdges(
-    const std::vector<size_t> &numberOfNodesPerSlot,
-    const std::filesystem::path &resultsPath,
-    const FlatPatternGeometry &patternGeometryAllEdges,
-    const std::vector<vcg::Point2i> &allPossibleEdges) {
-
-  std::unordered_set<size_t> coincideEdges =
-      computeCoincideEdges(numberOfNodesPerSlot);
-  // Export each coincide edge into a ply file
-  if (!coincideEdges.empty() && debugIsOn) {
-    auto coincideEdgesPath =
-        std::filesystem::path(resultsPath).append("Coincide_edges");
-    std::filesystem::create_directories(coincideEdgesPath);
-    for (auto coincideEdgeIndex : coincideEdges) {
-      FlatPatternGeometry::EdgeType e =
-          patternGeometryAllEdges.edge[coincideEdgeIndex];
-      FlatPatternGeometry singleEdgeMesh;
-      vcg::Point3d p0 = e.cP(0);
-      vcg::Point3d p1 = e.cP(1);
-      std::vector<vcg::Point3d> edgeVertices;
-      edgeVertices.push_back(p0);
-      edgeVertices.push_back(p1);
-      singleEdgeMesh.add(edgeVertices);
-      singleEdgeMesh.add(std::vector<vcg::Point2i>{vcg::Point2i{0, 1}});
-      singleEdgeMesh.savePly(std::filesystem::path(coincideEdgesPath)
-                                 .append(std::to_string(coincideEdgeIndex))
-                                 .string() +
-                             ".ply");
-    }
-  }
-  statistics.numberOfCoincideEdges = coincideEdges.size();
-
-  // Compute duplicate edges
-  std::unordered_set<size_t> duplicateEdges =
-      computeDuplicateEdges(numberOfNodesPerSlot);
-  if (!duplicateEdges.empty() && debugIsOn) {
-    // Export duplicate edges in a single ply file
-    auto duplicateEdgesPath =
-        std::filesystem::path(resultsPath).append("duplicate");
-    std::filesystem::create_directory(duplicateEdgesPath);
-    FlatPatternGeometry patternDuplicateEdges;
-    for (auto duplicateEdgeIndex : duplicateEdges) {
-      FlatPatternGeometry::EdgeType e =
-          patternGeometryAllEdges.edge[duplicateEdgeIndex];
-      vcg::Point3d p0 = e.cP(0);
-      vcg::Point3d p1 = e.cP(1);
-      vcg::tri::Allocator<FlatPatternGeometry>::AddEdge(
-          patternDuplicateEdges, p0, p1);
-    }
-    patternDuplicateEdges.savePly(
-        std::filesystem::path(duplicateEdgesPath).append("duplicateEdges.ply").string());
-  }
-  statistics.numberOfDuplicateEdges = duplicateEdges.size();
-
-  // Create the set of all possible edges without coincide and duplicate edges
-  std::vector<vcg::Point2i> validEdges;
-  for (size_t edgeIndex = 0; edgeIndex < allPossibleEdges.size(); edgeIndex++) {
-    if (coincideEdges.count(edgeIndex) == 0 &&
-        duplicateEdges.count(edgeIndex) == 0) {
-      validEdges.push_back(allPossibleEdges[edgeIndex]);
-    }
-  }
-
-  return validEdges;
-}
-
-void TopologyEnumerator::computeValidPatterns(
-    const std::vector<size_t> &numberOfNodesPerSlot,
-    const size_t &numberOfDesiredEdges,
-    const std::filesystem::path &resultsPath,
-    const std::vector<vcg::Point3d> &allVertices,
-    const std::unordered_map<size_t, std::unordered_set<size_t>>
-        &intersectingEdges,
-    const std::vector<vcg::Point2i> &validEdges, PatternSet &patternsSet) {
-  assert(numberOfNodesPerSlot.size() == 7);
-
-  // Iterate over all patterns which have numberOfDesiredEdges edges from
-  // from the validEdges Identify patterns that contain dangling edges
-  const bool enoughValidEdgesExist = validEdges.size() >= numberOfDesiredEdges;
-  if (!enoughValidEdgesExist) {
-    std::filesystem::remove_all(resultsPath); // delete previous results folder
-    return;
-  }
-  assert(enoughValidEdgesExist);
-
-  // Create pattern result paths
-  auto validPatternsPath = std::filesystem::path(resultsPath).append("Valid");
-  std::filesystem::create_directory(validPatternsPath);
-
-  const size_t numberOfPatterns = FlatPatternGeometry::binomialCoefficient(
-      validEdges.size(), numberOfDesiredEdges);
-  statistics.numberOfPatterns = numberOfPatterns;
-
-  // Initialize pattern binary representation
-  std::string patternBinaryRepresentation;
-  patternBinaryRepresentation = std::string(numberOfDesiredEdges, '1');
-  patternBinaryRepresentation +=
-      std::string(validEdges.size() - numberOfDesiredEdges, '0');
-  std::sort(patternBinaryRepresentation.begin(),
-            patternBinaryRepresentation.end());
-  size_t patternIndex = 0;
-  do {
-    patternIndex++;
-    const std::string patternName = std::to_string(patternIndex);
-    //    std::cout << "Pattern name:" + patternBinaryRepresentation <<
-    //    std::endl; isValidPattern(patternBinaryRepresentation, validEdges,
-    //                   numberOfDesiredEdges);
-    // Create the geometry of the pattern
-    // Compute the pattern edges from the binary representation
-    std::vector<vcg::Point2i> patternEdges(numberOfDesiredEdges);
-    size_t patternEdgeIndex = 0;
-    for (size_t validEdgeIndex = 0;
-         validEdgeIndex < patternBinaryRepresentation.size();
-         validEdgeIndex++) {
-      if (patternBinaryRepresentation[validEdgeIndex] == '1') {
-        assert(patternEdgeIndex < numberOfDesiredEdges);
-        patternEdges[patternEdgeIndex++] = validEdges[validEdgeIndex];
-      }
-    }
-    Pattern pattern;
-    pattern.edges = patternEdges;
-
-    FlatPatternGeometry patternGeometry;
-    patternGeometry.add(allVertices, patternEdges);
-
-    // Check if pattern contains intersecting edges
-    const bool patternContainsIntersectingEdges =
-        patternGeometry.hasIntersectingEdges(patternBinaryRepresentation,
-                                             intersectingEdges);
-    // Export the tiled ply file if it contains intersecting edges
-    if (patternContainsIntersectingEdges) {
-      // create the tiled geometry of the pattern
-      statistics.numberOfPatternsWithIntersectingEdges++;
-      if (debugIsOn) {
-        if (savePlyFiles) {
-          FlatPatternGeometry tiledPatternGeometry =
-              FlatPatternGeometry::createTile(patternGeometry);
-          auto intersectingPatternsPath =
-              std::filesystem::path(resultsPath).append("Intersecting");
-          std::filesystem::create_directory(intersectingPatternsPath);
-          patternGeometry.savePly(
-              std::filesystem::path(intersectingPatternsPath)
-                  .append(patternName)
-                  .string() +
-              ".ply");
-          tiledPatternGeometry.savePly(
-              std::filesystem::path(intersectingPatternsPath)
-                  .append(patternName + "_tiled")
-                  .string() +
-              ".ply");
-        }
-        pattern.labels.push_back(PatternLabel::IntersectingEdges);
-      } else {
-        continue; // should be uncommented in order to improve performance
-      }
-    }
-
-    // Compute the tiled valence
-    const bool tiledPatternHasDanglingEdges = patternGeometry.hasDanglingEdges(
-        numberOfNodesPerSlot); // marks the nodes with valence>=1
-                               // Create the tiled geometry of the pattern
-    const bool hasFloatingComponents =
-        !patternGeometry.isFullyConnectedWhenTiled();
-    FlatPatternTopology topology(numberOfNodesPerSlot, patternEdges);
-    const bool hasArticulationPoints = topology.containsArticulationPoints();
-    FlatPatternGeometry tiledPatternGeometry =
-        FlatPatternGeometry::createTile(
-            patternGeometry); // the marked nodes of hasDanglingEdges are
-    // duplicated here
-    // Check dangling edges with vcg method
-    //    const bool vcg_tiledPatternHasDangling =
-    //        tiledPatternGeometry.hasUntiledDanglingEdges();
-    if (tiledPatternHasDanglingEdges /*&& !hasFloatingComponents &&
-        !hasArticulationPoints*/) {
-      statistics.numberOfPatternsWithADanglingEdgeOrNode++;
-      if (debugIsOn) {
-        if (savePlyFiles) {
-          auto danglingEdgesPath =
-              std::filesystem::path(resultsPath).append("Dangling");
-          std::filesystem::create_directory(danglingEdgesPath);
-          patternGeometry.savePly(std::filesystem::path(danglingEdgesPath)
-                                      .append(patternName)
-                                      .string() +
-                                  ".ply");
-          tiledPatternGeometry.savePly(std::filesystem::path(danglingEdgesPath)
-                                           .append(patternName + "_tiled")
-                                           .string() +
-                                       ".ply");
-        }
-        pattern.labels.push_back(PatternLabel::DanglingEdge);
-      } else {
-        continue;
-      }
-    }
-
-    if (hasFloatingComponents /*&& !hasArticulationPoints &&
-        !tiledPatternHasDanglingEdges*/) {
-      statistics.numberOfPatternsWithMoreThanASingleCC++;
-      if (debugIsOn) {
-        if (savePlyFiles) {
-          auto moreThanOneCCPath =
-              std::filesystem::path(resultsPath).append("MoreThanOneCC");
-          std::filesystem::create_directory(moreThanOneCCPath);
-          patternGeometry.savePly(std::filesystem::path(moreThanOneCCPath)
-                                      .append(patternName)
-                                      .string() +
-                                  ".ply");
-          tiledPatternGeometry.savePly(std::filesystem::path(moreThanOneCCPath)
-                                           .append(patternName + "_tiled")
-                                           .string() +
-                                       ".ply");
-        }
-        pattern.labels.push_back(PatternLabel::MultipleCC);
-      } else {
-        continue;
-      }
-    }
-
-    if (hasArticulationPoints /*&& !hasFloatingComponents &&
-        !tiledPatternHasDanglingEdges*/) {
-      statistics.numberOfPatternsWithArticulationPoints++;
-      if (exportArticulationPointsPatterns || debugIsOn) {
-        if (savePlyFiles) {
-          auto articulationPointsPath =
-              std::filesystem::path(resultsPath).append("ArticulationPoints");
-          std::filesystem::create_directory(articulationPointsPath);
-          patternGeometry.savePly(std::filesystem::path(articulationPointsPath)
-                                      .append(patternName)
-                                      .string() +
-                                  ".ply");
-          tiledPatternGeometry.savePly(
-              std::filesystem::path(articulationPointsPath)
-                  .append(patternName + "_tiled")
-                  .string() +
-              ".ply");
-
-          //          std::cout << "Pattern:" << patternName << std::endl;
-        }
-        pattern.labels.push_back(PatternLabel::ArticulationPoints);
-      } else {
-        continue;
-      }
-    }
-
-    const bool isValidPattern =
-        !patternContainsIntersectingEdges && !tiledPatternHasDanglingEdges &&
-        !hasFloatingComponents && !hasArticulationPoints;
-    if (isValidPattern) {
-      statistics.numberOfValidPatterns++;
-      if (savePlyFiles) {
-        //        if (numberOfDesiredEdges == 4) {
-        //          std::cout << "Saving:"
-        //                    << std::filesystem::path(validPatternsPath)
-        //                               .append(patternName)
-        //                               .string() +
-        //                           ".ply"
-        //                    << std::endl;
-        //        }
-        patternGeometry.savePly(std::filesystem::path(validPatternsPath)
-                                    .append(patternName)
-                                    .string() +
-                                ".ply");
-        tiledPatternGeometry.savePly(std::filesystem::path(validPatternsPath)
-                                         .append(patternName + "_tiled")
-                                         .string() +
-                                     ".ply");
-      }
-      pattern.labels.push_back(PatternLabel::Valid);
-    }
-
-    assert(!pattern.labels.empty());
-    patternsSet.patterns.push_back(pattern);
-    //    assert(vcg_tiledPatternHasDangling == tiledPatternHasDanglingEdges);
-  } while (std::next_permutation(patternBinaryRepresentation.begin(),
-                                 patternBinaryRepresentation.end()));
-}
-
-std::vector<size_t> TopologyEnumerator::getCoincideEdges(
-    const std::vector<size_t> &edgeNodeIndices) const {
-  std::vector<size_t> coincideEdges;
-  if (edgeNodeIndices.size() < 3)
-    return coincideEdges;
-  for (size_t edgeNodeIndex = 0; edgeNodeIndex < edgeNodeIndices.size() - 2;
-       edgeNodeIndex++) {
-    const size_t &firstNodeIndex = edgeNodeIndices[edgeNodeIndex];
-    for (size_t secondEdgeNodeIndex = edgeNodeIndex + 2;
-         secondEdgeNodeIndex < edgeNodeIndices.size(); secondEdgeNodeIndex++) {
-      const size_t &secondNodeIndex = edgeNodeIndices[secondEdgeNodeIndex];
-      coincideEdges.push_back(getEdgeIndex(firstNodeIndex, secondNodeIndex));
-    }
-  }
-  return coincideEdges;
-}
-
-bool TopologyEnumerator::isValidPattern(
-    const std::string &patternBinaryRepresentation,
-    const std::vector<vcg::Point2i> &validEdges,
-    const size_t &numberOfDesiredEdges) const {
-  return true;
-}

topologyenumerator.hpp

@@ -157,7 +157,7 @@ private:
   std::vector<vcg::Point2i>
   getValidEdges(const std::vector<size_t> &numberOfNodesPerSlot,
                 const std::filesystem::path &resultsPath,
-                const FlatPatternGeometry &patternGeometryAllEdges,
+                const PatternGeometry &patternGeometryAllEdges,
                 const std::vector<vcg::Point2i> &allPossibleEdges);
   std::unordered_set<size_t> computeDuplicateEdges();
   std::unordered_set<size_t>
@@ -175,6 +175,7 @@ private:
       const std::vector<vcg::Point3d> &vertices,
       const std::unordered_map<size_t, std::unordered_set<size_t>>
           &intersectingEdges,
-      const std::vector<vcg::Point2i> &validEdges, PatternSet &patternsSet);
+      const std::vector<vcg::Point2i> &validEdges);
+}ges, PatternSet &patternsSet);
 };
-#endif // TOPOLOGYENUMERATOR_HPP
+#

trianglepatterngeometry.cpp

@@ -8,10 +8,10 @@
 #include <vcg/space/intersection2.h>
 #include <wrap/io_trimesh/export.h>
 
-size_t FlatPatternGeometry::computeTiledValence(
+size_t PatternGeometry::computeTiledValence(
     const size_t &nodeIndex,
     const std::vector<size_t> &numberOfNodesPerSlot) const {
-  std::vector<FlatPatternGeometry::EdgeType *> connectedEdges;
+  std::vector<PatternGeometry::EdgeType *> connectedEdges;
   vcg::edge::VEStarVE(&vert[nodeIndex], connectedEdges);
   const size_t nodeValence = connectedEdges.size();
   assert(nodeSlot.count(nodeIndex) != 0);
@@ -25,8 +25,7 @@ size_t FlatPatternGeometry::computeTiledValence(
     } else {
       correspondingNodeIndex = nodeIndex - 1;
     }
-    std::vector<FlatPatternGeometry::EdgeType *>
+    std::vector<PatternGeometry::EdgeType *> connectedEdgesCorrespondingNode;
-        connectedEdgesCorrespondingNode;
     vcg::edge::VEStarVE(&vert[correspondingNodeIndex],
                         connectedEdgesCorrespondingNode);
     size_t correspondingNodeValence = connectedEdgesCorrespondingNode.size();
@@ -55,8 +54,7 @@ size_t FlatPatternGeometry::computeTiledValence(
                                                0);
     }
     assert(correspondingNodeIndex < vn);
-    std::vector<FlatPatternGeometry::EdgeType *>
+    std::vector<PatternGeometry::EdgeType *> connectedEdgesCorrespondingNode;
-        connectedEdgesCorrespondingNode;
     vcg::edge::VEStarVE(&vert[correspondingNodeIndex],
                         connectedEdgesCorrespondingNode);
     size_t correspondingNodeValence = connectedEdgesCorrespondingNode.size();
@@ -72,15 +70,13 @@ size_t FlatPatternGeometry::computeTiledValence(
   return 0;
 }
 
-size_t FlatPatternGeometry::getFanSize() const { return fanSize; }
+size_t PatternGeometry::getFanSize() const { return fanSize; }
 
-double FlatPatternGeometry::getTriangleEdgeSize() const {
+double PatternGeometry::getTriangleEdgeSize() const { return triangleEdgeSize; }
-  return triangleEdgeSize;
-}
 
-FlatPatternGeometry::FlatPatternGeometry() {}
+PatternGeometry::PatternGeometry() {}
 
-std::vector<vcg::Point3d> FlatPatternGeometry::getVertices() const {
+std::vector<vcg::Point3d> PatternGeometry::getVertices() const {
   std::vector<VCGEdgeMesh::CoordType> verts(VN());
   for (size_t vi = 0; vi < VN(); vi++) {
     verts[vi] = vert[vi].cP();
@@ -88,73 +84,70 @@ std::vector<vcg::Point3d> FlatPatternGeometry::getVertices() const {
   return verts;
 }
 
-FlatPatternGeometry
+PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern) {
-FlatPatternGeometry::createTile(FlatPatternGeometry &pattern) {
 
-  const size_t fanSize = FlatPatternGeometry().getFanSize();
+  const size_t fanSize = PatternGeometry().getFanSize();
-  FlatPatternGeometry fan(createFan(pattern));
+  PatternGeometry fan(createFan(pattern));
-  FlatPatternGeometry tile(fan);
+  PatternGeometry tile(fan);
 
   if (fanSize % 2 == 1) {
     vcg::Matrix44d R;
     auto rotationAxis = vcg::Point3d(0, 0, 1);
     R.SetRotateDeg(180, rotationAxis);
-    vcg::tri::UpdatePosition<FlatPatternGeometry>::Matrix(fan, R);
+    vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, R);
   }
   vcg::Matrix44d T;
   const double centerAngle = 2 * M_PI / fanSize;
-  const double triangleHeight = std::sin((M_PI - centerAngle) / 2) *
+  const double triangleHeight =
-                                FlatPatternGeometry().triangleEdgeSize;
+      std::sin((M_PI - centerAngle) / 2) * PatternGeometry().triangleEdgeSize;
   T.SetTranslate(0, -2 * triangleHeight, 0);
-  vcg::tri::UpdatePosition<FlatPatternGeometry>::Matrix(fan, T);
+  vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, T);
 
-  FlatPatternGeometry fanOfFan = createFan(fan);
+  PatternGeometry fanOfFan = createFan(fan);
-  vcg::tri::Append<FlatPatternGeometry, FlatPatternGeometry>::Mesh(tile,
+  vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(tile, fanOfFan);
-                                                                   fanOfFan);
+  vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(tile, 0.0000005);
-  vcg::tri::Clean<FlatPatternGeometry>::MergeCloseVertex(tile, 0.0000005);
+  vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(tile);
-  vcg::tri::Allocator<FlatPatternGeometry>::CompactEveryVector(tile);
+  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(tile);
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::VertexEdge(tile);
+  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(tile);
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::EdgeEdge(tile);
 
   for (size_t vi = 0; vi < pattern.vn; vi++) {
     tile.vert[vi].C() = vcg::Color4b::Blue;
   }
 
+  tile.updateEigenEdgeAndVertices();
   return tile;
 }
 
-FlatPatternGeometry
+PatternGeometry PatternGeometry::createFan(PatternGeometry &pattern) {
-FlatPatternGeometry::createFan(FlatPatternGeometry &pattern) {
+  const size_t fanSize = PatternGeometry().getFanSize();
-  const size_t fanSize = FlatPatternGeometry().getFanSize();
+  PatternGeometry fan(pattern);
-  FlatPatternGeometry fan(pattern);
+  PatternGeometry rotatedPattern(pattern);
-  FlatPatternGeometry rotatedPattern(pattern);
   for (int rotationCounter = 1; rotationCounter < fanSize; rotationCounter++) {
     vcg::Matrix44d R;
     auto rotationAxis = vcg::Point3d(0, 0, 1);
     R.SetRotateDeg(360 / fanSize, rotationAxis);
-    vcg::tri::UpdatePosition<FlatPatternGeometry>::Matrix(rotatedPattern, R);
+    vcg::tri::UpdatePosition<PatternGeometry>::Matrix(rotatedPattern, R);
-    vcg::tri::Append<FlatPatternGeometry, FlatPatternGeometry>::Mesh(
+    vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(fan,
-        fan, rotatedPattern);
+                                                             rotatedPattern);
   }
   return fan;
 }
 
-FlatPatternGeometry::FlatPatternGeometry(FlatPatternGeometry &other) {
+PatternGeometry::PatternGeometry(PatternGeometry &other) {
-  vcg::tri::Append<FlatPatternGeometry, FlatPatternGeometry>::MeshCopy(*this,
+  vcg::tri::Append<PatternGeometry, PatternGeometry>::MeshCopy(*this, other);
-                                                                       other);
   this->vertices = other.getVertices();
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::VertexEdge(*this);
+  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::EdgeEdge(*this);
+  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
 }
 
-bool FlatPatternGeometry::savePly(const std::string plyFilename) {
+bool PatternGeometry::savePly(const std::string plyFilename) {
 
-  int returnValue = vcg::tri::io::ExporterPLY<FlatPatternGeometry>::Save(
+  int returnValue = vcg::tri::io::ExporterPLY<PatternGeometry>::Save(
       *this, plyFilename.c_str(),
       vcg::tri::io::Mask::IOM_EDGEINDEX | vcg::tri::io::Mask::IOM_VERTCOLOR,
       false);
   if (returnValue != 0) {
-    std::cerr << vcg::tri::io::ExporterPLY<FlatPatternGeometry>::ErrorMsg(
+    std::cerr << vcg::tri::io::ExporterPLY<PatternGeometry>::ErrorMsg(
                      returnValue)
               << std::endl;
   }
@@ -162,44 +155,44 @@ bool FlatPatternGeometry::savePly(const std::string plyFilename) {
   return static_cast<bool>(returnValue);
 }
 
-void FlatPatternGeometry::add(const std::vector<vcg::Point3d> &vertices) {
+void PatternGeometry::add(const std::vector<vcg::Point3d> &vertices) {
   this->vertices = vertices;
   std::for_each(vertices.begin(), vertices.end(), [&](const vcg::Point3d &p) {
-    vcg::tri::Allocator<FlatPatternGeometry>::AddVertex(*this, p);
+    vcg::tri::Allocator<PatternGeometry>::AddVertex(*this, p);
   });
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::VertexEdge(*this);
+  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::EdgeEdge(*this);
+  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
   updateEigenEdgeAndVertices();
 }
 
-void FlatPatternGeometry::add(const std::vector<vcg::Point2i> &edges) {
+void PatternGeometry::add(const std::vector<vcg::Point2i> &edges) {
   std::for_each(edges.begin(), edges.end(), [&](const vcg::Point2i &e) {
-    vcg::tri::Allocator<FlatPatternGeometry>::AddEdge(*this, e[0], e[1]);
+    vcg::tri::Allocator<PatternGeometry>::AddEdge(*this, e[0], e[1]);
   });
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::VertexEdge(*this);
+  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::EdgeEdge(*this);
+  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
   updateEigenEdgeAndVertices();
 }
 
-void FlatPatternGeometry::add(const std::vector<vcg::Point3d> &vertices,
+void PatternGeometry::add(const std::vector<vcg::Point3d> &vertices,
-                              const std::vector<vcg::Point2i> &edges) {
+                          const std::vector<vcg::Point2i> &edges) {
   add(vertices);
   add(edges);
   updateEigenEdgeAndVertices();
 }
 
-void FlatPatternGeometry::add(const std::vector<size_t> &numberOfNodesPerSlot,
+void PatternGeometry::add(const std::vector<size_t> &numberOfNodesPerSlot,
-                              const std::vector<vcg::Point2i> &edges) {
+                          const std::vector<vcg::Point2i> &edges) {
   assert(numberOfNodesPerSlot.size() == 7);
   auto vertices =
       constructVertexVector(numberOfNodesPerSlot, fanSize, triangleEdgeSize);
   add(vertices, edges);
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::VertexEdge(*this);
+  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::EdgeEdge(*this);
+  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
   updateEigenEdgeAndVertices();
 }
 
-std::vector<vcg::Point3d> FlatPatternGeometry::constructVertexVector(
+std::vector<vcg::Point3d> PatternGeometry::constructVertexVector(
     const std::vector<size_t> &numberOfNodesPerSlot, const size_t &fanSize,
     const double &triangleEdgeSize) {
 
@@ -294,7 +287,7 @@ std::vector<vcg::Point3d> FlatPatternGeometry::constructVertexVector(
   return vertices;
 }
 
-void FlatPatternGeometry::constructNodeToTiledValenceMap(
+void PatternGeometry::constructNodeToTiledValenceMap(
     const std::vector<size_t> &numberOfNodesPerSlot) {
   for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
     const size_t tiledValence =
@@ -303,7 +296,7 @@ void FlatPatternGeometry::constructNodeToTiledValenceMap(
   }
 }
 
-bool FlatPatternGeometry::hasDanglingEdges(
+bool PatternGeometry::hasDanglingEdges(
     const std::vector<size_t> &numberOfNodesPerSlot) {
   if (nodeSlot.empty()) {
     FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeSlot);
@@ -325,7 +318,7 @@ bool FlatPatternGeometry::hasDanglingEdges(
   return hasDanglingEdges;
 }
 
-bool FlatPatternGeometry::hasUntiledDanglingEdges() {
+bool PatternGeometry::hasUntiledDanglingEdges() {
   //  vcg::tri::Clean<TrianglePatternGeometry>::MergeCloseVertex(*this,
   //  0.0000005);
   //  vcg::tri::Allocator<TrianglePatternGeometry>::CompactEveryVector(*this);
@@ -333,7 +326,7 @@ bool FlatPatternGeometry::hasUntiledDanglingEdges() {
   //  vcg::tri::UpdateTopology<TrianglePatternGeometry>::EdgeEdge(*this);
   bool hasDanglingEdges = false;
   for (size_t vi = 0; vi < vn; vi++) {
-    std::vector<FlatPatternGeometry::EdgeType *> connectedEdges;
+    std::vector<PatternGeometry::EdgeType *> connectedEdges;
     vcg::edge::VEStarVE(&vert[vi], connectedEdges);
     const size_t nodeValence = connectedEdges.size();
     if (nodeValence == 1) {
@@ -350,7 +343,7 @@ bool FlatPatternGeometry::hasUntiledDanglingEdges() {
 
 // TODO: The function expects that the numberOfNodesPerSlot follows a specific
 // format and that the vertex container was populated in a particular order.
-void FlatPatternGeometry::constructCorresponginNodeMap(
+void PatternGeometry::constructCorresponginNodeMap(
     const std::vector<size_t> &numberOfNodesPerSlot) {
   assert(vn != 0 && !nodeSlot.empty() && correspondingNode.empty() &&
          numberOfNodesPerSlot.size() == 7);
@@ -382,7 +375,7 @@ void FlatPatternGeometry::constructCorresponginNodeMap(
   }
 }
 
-bool FlatPatternGeometry::isFullyConnectedWhenTiled() {
+bool PatternGeometry::isFullyConnectedWhenTiled() {
   assert(vn != 0 /* && !correspondingNode.empty()*/);
   //  TrianglePatternGeometry copyOfPattern(*this);
 
@@ -392,7 +385,7 @@ bool FlatPatternGeometry::isFullyConnectedWhenTiled() {
   for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) {
     if (nodeSlot[nodeIndex] == 1 || nodeSlot[nodeIndex] == 4 ||
         nodeSlot[nodeIndex] == 2) {
-      std::vector<FlatPatternGeometry::EdgeType *> connectedEdges;
+      std::vector<PatternGeometry::EdgeType *> connectedEdges;
       vcg::edge::VEStarVE(&vert[nodeIndex], connectedEdges);
       const size_t nodeValence = connectedEdges.size();
       if (nodeValence != 0) {
@@ -405,15 +398,14 @@ bool FlatPatternGeometry::isFullyConnectedWhenTiled() {
     return false;
   }
 
-  FlatPatternGeometry fanedPattern = createFan(*this);
+  PatternGeometry fanedPattern = createFan(*this);
-  vcg::tri::Clean<FlatPatternGeometry>::MergeCloseVertex(fanedPattern,
+  vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(fanedPattern, 0.000000005);
-                                                         0.000000005);
+  vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(fanedPattern);
-  vcg::tri::Allocator<FlatPatternGeometry>::CompactEveryVector(fanedPattern);
+  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(fanedPattern);
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::VertexEdge(fanedPattern);
+  vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(fanedPattern);
-  vcg::tri::UpdateTopology<FlatPatternGeometry>::EdgeEdge(fanedPattern);
+  std::vector<std::pair<int, PatternGeometry::EdgePointer>> eCC;
-  std::vector<std::pair<int, FlatPatternGeometry::EdgePointer>> eCC;
+  vcg::tri::Clean<PatternGeometry>::edgeMeshConnectedComponents(fanedPattern,
-  vcg::tri::Clean<FlatPatternGeometry>::edgeMeshConnectedComponents(
+                                                                eCC);
-      fanedPattern, eCC);
 
   const bool sideInterfaceIsConnected = 1 == eCC.size();
   if (!sideInterfaceIsConnected) {
@@ -490,7 +482,7 @@ bool FlatPatternGeometry::isFullyConnectedWhenTiled() {
   return true;
 }
 
-bool FlatPatternGeometry::hasIntersectingEdges(
+bool PatternGeometry::hasIntersectingEdges(
     const std::string &patternBinaryRepresentation,
     const std::unordered_map<size_t, std::unordered_set<size_t>>
         &intersectingEdges) {
@@ -520,7 +512,7 @@ bool FlatPatternGeometry::hasIntersectingEdges(
 }
 
 std::unordered_map<size_t, std::unordered_set<size_t>>
-FlatPatternGeometry::getIntersectingEdges(
+PatternGeometry::getIntersectingEdges(
     size_t &numberOfIntersectingEdgePairs) const {
   std::unordered_map<size_t, std::unordered_set<size_t>> intersectingEdges;
   numberOfIntersectingEdgePairs = 0;
@@ -569,3 +561,434 @@ FlatPatternGeometry::getIntersectingEdges(
   }
   return intersectingEdges;
 }
+
+PatternGeometry::PatternGeometry(const std::string &filename,
+                                 bool addNormalsIfAbsent) {
+  if (!std::filesystem::exists(std::filesystem::path(filename))) {
+    assert(false);
+    std::cerr << "No flat pattern with name " << filename << std::endl;
+    return;
+  }
+  if (!load(filename)) {
+    assert(false);
+    std::cerr << "File could not be loaded  " << filename << std::endl;
+    return;
+  }
+  if (addNormalsIfAbsent) {
+    addNormals();
+  }
+
+  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
+  baseTriangleHeight = computeBaseTriangleHeight();
+
+  updateEigenEdgeAndVertices();
+}
+
+double PatternGeometry::computeBaseTriangleHeight() const {
+  return (vert[0].cP() - vert[3].cP()).Norm();
+}
+
+void PatternGeometry::addNormals() {
+  bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001;
+  if (normalsAreAbsent) {
+    for (auto &v : vert) {
+      v.N() = CoordType(0, 0, 1);
+    }
+  }
+}
+
+PatternGeometry::PatternGeometry(
+    const std::vector<size_t> &numberOfNodesPerSlot,
+    const std::vector<vcg::Point2i> &edges) {
+  add(numberOfNodesPerSlot, edges);
+  addNormals();
+  baseTriangleHeight = computeBaseTriangleHeight();
+  updateEigenEdgeAndVertices();
+}
+
+bool PatternGeometry::createHoneycombAtom() {
+  VCGEdgeMesh honeycombQuarter;
+  const VCGEdgeMesh::CoordType n(0, 0, 1);
+  const double H = 0.2;
+  const double height = 1.5 * H;
+  const double width = 0.2;
+  const double theta = 70;
+  const double dy = tan(vcg::math::ToRad(90 - theta)) * width / 2;
+  vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(
+      honeycombQuarter, VCGEdgeMesh::CoordType(0, height / 2, 0), n);
+  vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(
+      honeycombQuarter, VCGEdgeMesh::CoordType(0, H / 2 - dy, 0), n);
+  vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(
+      honeycombQuarter, VCGEdgeMesh::CoordType(width / 2, H / 2, 0), n);
+  vcg::tri::Allocator<VCGEdgeMesh>::AddVertex(
+      honeycombQuarter, VCGEdgeMesh::CoordType(width / 2, 0, 0), n);
+  vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(honeycombQuarter, 0, 1);
+  vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(honeycombQuarter, 1, 2);
+  vcg::tri::Allocator<VCGEdgeMesh>::AddEdge(honeycombQuarter, 2, 3);
+
+  VCGEdgeMesh honeycombAtom;
+  // Top right
+  vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::MeshCopy(honeycombAtom,
+                                                       honeycombQuarter);
+  // Bottom right
+  vcg::Matrix44d rotM;
+  rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0));
+  vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(honeycombQuarter, rotM);
+  vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(honeycombAtom,
+                                                   honeycombQuarter);
+  // Bottom left
+  rotM.SetRotateDeg(180, vcg::Point3d(0, 1, 0));
+  vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(honeycombQuarter, rotM);
+  vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(honeycombAtom,
+                                                   honeycombQuarter);
+  // Top left
+  rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0));
+  vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(honeycombQuarter, rotM);
+  vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(honeycombAtom,
+                                                   honeycombQuarter);
+
+  for (VertexType &v : honeycombAtom.vert) {
+    v.P()[2] = 0;
+  }
+
+  return true;
+}
+
+void PatternGeometry::copy(PatternGeometry &copyFrom) {
+  VCGEdgeMesh::copy(copyFrom);
+  baseTriangleHeight = copyFrom.getBaseTriangleHeight();
+}
+
+void PatternGeometry::deleteDanglingEdges() {
+  for (VertexType &v : vert) {
+    std::vector<VCGEdgeMesh::EdgePointer> incidentElements;
+    vcg::edge::VEStarVE(&v, incidentElements);
+    if (incidentElements.size() == 1) {
+      vcg::tri::Allocator<VCGEdgeMesh>::DeleteEdge(*this, *incidentElements[0]);
+    }
+    if (incidentElements.size() == 1) {
+      vcg::tri::Allocator<VCGEdgeMesh>::DeleteVertex(*this, v);
+    }
+  }
+
+  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateVertex(*this);
+  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateEdge(*this);
+  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
+}
+
+void PatternGeometry::scale(const double &desiredBaseTriangleCentralEdgeSize) {
+  this->baseTriangleHeight = desiredBaseTriangleCentralEdgeSize;
+  const double baseTriangleCentralEdgeSize =
+      (vert[0].cP() - vert[3].cP()).Norm();
+  const double scaleRatio =
+      desiredBaseTriangleCentralEdgeSize / baseTriangleCentralEdgeSize;
+
+  vcg::tri::UpdatePosition<VCGEdgeMesh>::Scale(*this, scaleRatio);
+}
+
+void PatternGeometry::deleteDanglingVertices() {
+  vcg::tri::Allocator<PatternGeometry>::PointerUpdater<VertexPointer> pu;
+  deleteDanglingVertices(pu);
+}
+
+void PatternGeometry::deleteDanglingVertices(
+    vcg::tri::Allocator<PatternGeometry>::PointerUpdater<VertexPointer> &pu) {
+  for (VertexType &v : vert) {
+    std::vector<PatternGeometry::EdgePointer> incidentElements;
+    vcg::edge::VEStarVE(&v, incidentElements);
+    if (incidentElements.size() == 0 && !v.IsD()) {
+      vcg::tri::Allocator<PatternGeometry>::DeleteVertex(*this, v);
+    }
+  }
+
+  vcg::tri::Allocator<PatternGeometry>::CompactVertexVector(*this, pu);
+  vcg::tri::Allocator<PatternGeometry>::CompactEdgeVector(*this);
+
+  updateEigenEdgeAndVertices();
+}
+
+void PatternGeometry::tilePattern(VCGEdgeMesh& pattern, VCGPolyMesh &tileInto,const int& interfaceNodeIndex,
+                                  const bool &shouldDeleteDanglingEdges) {
+  double xOffset =
+      vcg::Distance(pattern.vert[0].cP(), pattern.vert[interfaceNodeIndex].cP()) /
+      std::tan(M_PI / 3);
+  CoordType patternCoord0 = pattern.vert[0].cP();
+  CoordType patternBottomRight =
+      pattern.vert[interfaceNodeIndex].cP() + CoordType(xOffset, 0, 0);
+  CoordType patternBottomLeft =
+      pattern.vert[interfaceNodeIndex].cP() - CoordType(xOffset, 0, 0);
+  std::vector<vcg::Point3d> patternTrianglePoints{
+      patternCoord0, patternBottomRight, patternBottomLeft};
+  CoordType pointOnPattern =
+      patternCoord0 + (patternBottomLeft - patternCoord0) ^
+      (patternBottomRight - patternCoord0);
+
+  std::vector<CoordType> faceCenters(FN());
+  VCGTriMesh tileIntoEdgeMesh;
+
+  for (VCGPolyMesh::FaceType &f : tileInto.face) {
+    std::vector<VCGPolyMesh::VertexType *> incidentVertices;
+    vcg::face::VFIterator<PFace> vfi(
+        &f, 0); // initialize the iterator to the first face
+                //    vcg::face::Pos p(vfi.F(), f.cV(0));
+                //    vcg::face::VVOrderedStarFF(p, incidentVertices);
+    size_t numberOfNodes = 0;
+    CoordType centerOfFace(0, 0, 0);
+    for (size_t vi = 0; vi < f.VN(); vi++) {
+      numberOfNodes++;
+      centerOfFace = centerOfFace + f.cP(vi);
+    }
+    centerOfFace /= f.VN();
+    vcg::tri::Allocator<VCGTriMesh>::AddVertex(tileIntoEdgeMesh, centerOfFace,
+                                               vcg::Color4b::Yellow);
+
+    //    const size_t vi = vcg::tri::Index<VCGPolyMesh>(tileInto, v);
+    //    std::cout << "vertex " << vi << " has incident vertices:" <<
+    //    std::endl;
+    for (size_t vi = 0; vi < f.VN(); vi++) {
+      //      size_t f = 0;
+      //                std::cout << vcg::tri::Index<VCGTriMesh>(tileInto,
+      //  / incidentVertices[f]) /                          << std::endl;
+
+      // Compute transformation matrix M
+      //        vcg::Matrix44d M;
+      std::vector<vcg::Point3d> meshTrianglePoints{
+          centerOfFace, f.cP(vi), vi + 1 == f.VN() ? f.cP(0) : f.cP(vi + 1)};
+      CoordType faceNormal = ((meshTrianglePoints[1] - meshTrianglePoints[0]) ^
+                              (meshTrianglePoints[2] - meshTrianglePoints[0]))
+                                 .Normalize();
+      auto fit = vcg::tri::Allocator<VCGTriMesh>::AddFace(
+          tileIntoEdgeMesh, meshTrianglePoints[0], meshTrianglePoints[1],
+          meshTrianglePoints[2]);
+      fit->N() = faceNormal;
+      CoordType pointOnTriMesh =
+          meshTrianglePoints[0] +
+              (meshTrianglePoints[1] - meshTrianglePoints[0]) ^
+          (meshTrianglePoints[2] - meshTrianglePoints[0]);
+      vcg::Matrix44d M;
+      //        vcg::ComputeRigidMatchMatrix(meshTrianglePoints,
+      //        patternTrianglePoints,
+      //                                     M);
+      vcg::Matrix44d A_prime;
+      A_prime[0][0] = meshTrianglePoints[0][0];
+      A_prime[1][0] = meshTrianglePoints[0][1];
+      A_prime[2][0] = meshTrianglePoints[0][2];
+      A_prime[3][0] = 1;
+      A_prime[0][1] = meshTrianglePoints[1][0];
+      A_prime[1][1] = meshTrianglePoints[1][1];
+      A_prime[2][1] = meshTrianglePoints[1][2];
+      A_prime[3][1] = 1;
+      A_prime[0][2] = meshTrianglePoints[2][0];
+      A_prime[1][2] = meshTrianglePoints[2][1];
+      A_prime[2][2] = meshTrianglePoints[2][2];
+      A_prime[3][2] = 1;
+      A_prime[0][3] = pointOnTriMesh[0];
+      A_prime[1][3] = pointOnTriMesh[1];
+      A_prime[2][3] = pointOnTriMesh[2];
+      A_prime[3][3] = 1;
+      vcg::Matrix44d A;
+      A[0][0] = patternTrianglePoints[0][0];
+      A[1][0] = patternTrianglePoints[0][1];
+      A[2][0] = patternTrianglePoints[0][2];
+      A[3][0] = 1;
+      A[0][1] = patternTrianglePoints[1][0];
+      A[1][1] = patternTrianglePoints[1][1];
+      A[2][1] = patternTrianglePoints[1][2];
+      A[3][1] = 1;
+      A[0][2] = patternTrianglePoints[2][0];
+      A[1][2] = patternTrianglePoints[2][1];
+      A[2][2] = patternTrianglePoints[2][2];
+      A[3][2] = 1;
+      A[0][3] = pointOnPattern[0];
+      A[1][3] = pointOnPattern[1];
+      A[2][3] = pointOnPattern[2];
+      A[3][3] = 1;
+      M = A_prime * vcg::Inverse(A);
+
+      VCGEdgeMesh transformedPattern;
+      vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::MeshCopy(transformedPattern,
+                                                           pattern);
+      vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(transformedPattern, M);
+      for (VertexType &v : transformedPattern.vert) {
+        v.N() = faceNormal;
+      }
+
+      vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(*this,
+                                                       transformedPattern);
+    }
+  }
+
+  //  vcg::tri::Clean<VCGEdgeMesh>::MergeCloseVertex(*this, 0.0000000001);
+  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateVertex(*this);
+  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateEdge(*this);
+  //  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
+
+  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveDuplicateVertex(*this);
+  vcg::tri::Clean<VCGEdgeMesh>::MergeCloseVertex(*this, 0.000061524);
+  vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
+
+  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveUnreferencedVertex(*this);
+  //  vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
+  //    vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateVertex(*this);
+  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateEdge(*this);
+  //  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
+  deleteDanglingVertices();
+  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
+  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateVertex(*this);
+  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveUnreferencedVertex(*this);
+  //  vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
+  //  deleteDanglingEdges();
+  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveUnreferencedVertex(*this);
+  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateVertex(*this);
+  //  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateEdge(*this);
+  //  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
+  updateEigenEdgeAndVertices();
+  savePly("tiledPattern.ply");
+}
+
+void PatternGeometry::createFan(const size_t &fanSize) {
+  PatternGeometry rotatedPattern;
+  vcg::tri::Append<PatternGeometry, PatternGeometry>::MeshCopy(rotatedPattern,
+                                                               *this);
+  for (int rotationCounter = 1; rotationCounter < fanSize; rotationCounter++) {
+    vcg::Matrix44d R;
+    auto rotationAxis = vcg::Point3d(0, 0, 1);
+    R.SetRotateDeg(360 / fanSize, rotationAxis);
+    vcg::tri::UpdatePosition<PatternGeometry>::Matrix(rotatedPattern, R);
+    vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(*this,
+                                                             rotatedPattern);
+  }
+
+  removeDuplicateVertices();
+  //  const double precision = 1e-4;
+  //    for (size_t vi = 0; vi < VN(); vi++) {
+  //      vert[vi].P()[0] = std::round(vert[vi].P()[0] * (1 / precision)) *
+  //      precision; vert[vi].P()[1] = std::round(vert[vi].P()[1] * (1 /
+  //      precision)) * precision; vert[vi].P()[2] = std::round(vert[vi].P()[2]
+  //      * (1 / precision)) * precision;
+  //    }
+
+  updateEigenEdgeAndVertices();
+}
+
+void PatternGeometry::removeDuplicateVertices() {
+  vcg::tri::Clean<VCGEdgeMesh>::MergeCloseVertex(*this, 0.0000000001);
+  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateVertex(*this);
+  vcg::tri::Clean<VCGEdgeMesh>::RemoveDegenerateEdge(*this);
+  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
+  vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
+}
+
+double PatternGeometry::getBaseTriangleHeight() const {
+  return baseTriangleHeight;
+}
+
+void PatternGeometry::tilePattern(PatternGeometry &pattern, VCGTriMesh &tileInto) {
+
+  const size_t middleIndex = 3;
+  double xOffset =
+      vcg::Distance(pattern.vert[0].cP(), pattern.vert[middleIndex].cP()) /
+      std::tan(M_PI / 3);
+  CoordType patternCoord0 = pattern.vert[0].cP();
+  CoordType patternBottomRight =
+      pattern.vert[middleIndex].cP() + CoordType(xOffset, 0, 0);
+  CoordType patternBottomLeft =
+      pattern.vert[middleIndex].cP() - CoordType(xOffset, 0, 0);
+  std::vector<vcg::Point3d> patternTrianglePoints{
+      patternCoord0, patternBottomRight, patternBottomLeft};
+  CoordType pointOnPattern =
+      patternCoord0 + (patternBottomLeft - patternCoord0) ^
+      (patternBottomRight - patternCoord0);
+
+  for (VCGTriMesh::VertexType &v : tileInto.vert) {
+    const auto centralNodeColor = vcg::Color4<unsigned char>(64, 64, 64, 255);
+    const bool isCentralNode = v.cC() == centralNodeColor;
+    if (isCentralNode) {
+      std::vector<VCGTriMesh::VertexType *> incidentVertices;
+      vcg::face::VFIterator<VCGTriMeshFace> vfi(
+          &v); // initialize the iterator tohe first face
+      vcg::face::Pos p(vfi.F(), &v);
+      vcg::face::VVOrderedStarFF(p, incidentVertices);
+
+      const size_t vi = vcg::tri::Index<VCGTriMesh>(tileInto, v);
+      std::cout << "vertex " << vi << " has incident vertices:" << std::endl;
+      for (size_t f = 0; f < incidentVertices.size(); f++) {
+        //      size_t f = 0;
+        std::cout << vcg::tri::Index<VCGTriMesh>(tileInto, incidentVertices[f])
+                  << std::endl;
+
+        // Compute transformation matrix M
+        //        vcg::Matrix44d M;
+        std::vector<vcg::Point3d> meshTrianglePoints{
+            v.cP(), incidentVertices[f]->cP(),
+            f + 1 == incidentVertices.size() ? incidentVertices[0]->cP()
+                                             : incidentVertices[f + 1]->cP()};
+        CoordType faceNormal =
+            ((meshTrianglePoints[1] - meshTrianglePoints[0]) ^
+             (meshTrianglePoints[2] - meshTrianglePoints[0]))
+                .Normalize();
+        CoordType pointOnTriMesh =
+            meshTrianglePoints[0] +
+                (meshTrianglePoints[1] - meshTrianglePoints[0]) ^
+            (meshTrianglePoints[2] - meshTrianglePoints[0]);
+        vcg::Matrix44d M;
+        //        vcg::ComputeRigidMatchMatrix(meshTrianglePoints,
+        //        patternTrianglePoints,
+        //                                     M);
+        vcg::Matrix44d A_prime;
+        A_prime[0][0] = meshTrianglePoints[0][0];
+        A_prime[1][0] = meshTrianglePoints[0][1];
+        A_prime[2][0] = meshTrianglePoints[0][2];
+        A_prime[3][0] = 1;
+        A_prime[0][1] = meshTrianglePoints[1][0];
+        A_prime[1][1] = meshTrianglePoints[1][1];
+        A_prime[2][1] = meshTrianglePoints[1][2];
+        A_prime[3][1] = 1;
+        A_prime[0][2] = meshTrianglePoints[2][0];
+        A_prime[1][2] = meshTrianglePoints[2][1];
+        A_prime[2][2] = meshTrianglePoints[2][2];
+        A_prime[3][2] = 1;
+        A_prime[0][3] = pointOnTriMesh[0];
+        A_prime[1][3] = pointOnTriMesh[1];
+        A_prime[2][3] = pointOnTriMesh[2];
+        A_prime[3][3] = 1;
+        vcg::Matrix44d A;
+        A[0][0] = patternTrianglePoints[0][0];
+        A[1][0] = patternTrianglePoints[0][1];
+        A[2][0] = patternTrianglePoints[0][2];
+        A[3][0] = 1;
+        A[0][1] = patternTrianglePoints[1][0];
+        A[1][1] = patternTrianglePoints[1][1];
+        A[2][1] = patternTrianglePoints[1][2];
+        A[3][1] = 1;
+        A[0][2] = patternTrianglePoints[2][0];
+        A[1][2] = patternTrianglePoints[2][1];
+        A[2][2] = patternTrianglePoints[2][2];
+        A[3][2] = 1;
+        A[0][3] = pointOnPattern[0];
+        A[1][3] = pointOnPattern[1];
+        A[2][3] = pointOnPattern[2];
+        A[3][3] = 1;
+        M = A_prime * vcg::Inverse(A);
+
+        VCGEdgeMesh transformedPattern;
+        vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::MeshCopy(transformedPattern,
+                                                             pattern);
+        vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(transformedPattern, M);
+        for (VertexType &v : transformedPattern.vert) {
+          v.N() = faceNormal;
+        }
+
+        vcg::tri::Append<VCGEdgeMesh, VCGEdgeMesh>::Mesh(*this,
+                                                         transformedPattern);
+      }
+    }
+  }
+
+  vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
+  deleteDanglingVertices();
+  deleteDanglingEdges();
+  vcg::tri::Allocator<VCGEdgeMesh>::CompactEveryVector(*this);
+
+  updateEigenEdgeAndVertices();
+}

trianglepatterngeometry.hpp

@@ -4,8 +4,10 @@
 #include <string>
 #include <unordered_map>
 #include <unordered_set>
+#include "vcgtrimesh.hpp"
+#include "polymesh.hpp"
 
-class FlatPatternGeometry : public VCGEdgeMesh {
+class PatternGeometry : public VCGEdgeMesh {
 private:
   size_t
   computeTiledValence(const size_t &nodeIndex,
@@ -13,6 +15,12 @@ private:
   size_t getNodeValence(const size_t &nodeIndex) const;
   size_t getNodeSlot(const size_t &nodeIndex) const;
 
+  void addNormals();
+  void deleteDanglingEdges();
+  void removeDuplicateVertices();
+  double baseTriangleHeight;
+  double computeBaseTriangleHeight() const;
+
   const size_t fanSize{6};
   std::vector<VCGEdgeMesh::CoordType> vertices;
   const double triangleEdgeSize{1}; // radius edge
@@ -27,12 +35,12 @@ private:
       const std::vector<size_t> &numberOfNodesPerSlot);
 
 public:
-  FlatPatternGeometry();
+  PatternGeometry();
   /*The following function should be a copy constructor with
    * a const argument but this is impossible due to the
    * vcglib interface.
    * */
-  FlatPatternGeometry(FlatPatternGeometry &other);
+  PatternGeometry(PatternGeometry &other);
   bool savePly(const std::string plyFilename);
   void add(const std::vector<vcg::Point3d> &vertices);
   void add(const std::vector<vcg::Point2i> &edges);
@@ -45,8 +53,8 @@ public:
                         const size_t &fanSize, const double &triangleEdgeSize);
   bool hasDanglingEdges(const std::vector<size_t> &numberOfNodesPerSlot);
   std::vector<vcg::Point3d> getVertices() const;
-  static FlatPatternGeometry createFan(FlatPatternGeometry &pattern);
+  static PatternGeometry createFan(PatternGeometry &pattern);
-  static FlatPatternGeometry createTile(FlatPatternGeometry &pattern);
+  static PatternGeometry createTile(PatternGeometry &pattern);
   double getTriangleEdgeSize() const;
   bool hasUntiledDanglingEdges();
   std::unordered_map<size_t, std::unordered_set<size_t>>
@@ -64,5 +72,31 @@ public:
           &intersectingEdges);
   bool isPatternValid();
   size_t getFanSize() const;
+  void add(const std::vector<vcg::Point2d> &vertices,
+           const std::vector<vcg::Point2i> &edges);
+
+  PatternGeometry(const std::vector<size_t> &numberOfNodesPerSlot,
+                  const std::vector<vcg::Point2i> &edges);
+  PatternGeometry(const std::string &filename, bool addNormalsIfAbsent = true);
+
+  bool createHoneycombAtom();
+  void copy(PatternGeometry &copyFrom);
+
+  void tilePattern(PatternGeometry &pattern, VCGTriMesh &tileInto);
+  void tilePattern(VCGEdgeMesh &pattern, VCGPolyMesh &tileInto, const int &interfaceNodeIndex,
+                   const bool &shouldDeleteDanglingEdges);
+
+  void createFan(const size_t &fanSize = 6);
+
+  void deleteDanglingVertices(
+      vcg::tri::Allocator<PatternGeometry>::PointerUpdater<VertexPointer> &pu);
+  void deleteDanglingVertices();
+  void scale(const double &desiredBaseTriangleCentralEdgeSize);
+
+  double getBaseTriangleHeight() const;
+
+  PatternGeometry(const std::vector<vcg::Point2d> &vertices,
+                  const std::vector<vcg::Point2i> &edges);
+
 };
 #endif // FLATPATTERNGEOMETRY_HPP

utilities.hpp

@@ -13,6 +13,11 @@ struct Vector6d : public std::array<double, 6> {
     }
   }
 
+  Vector6d(const std::vector<double> &v) {
+    assert(v.size() == 6);
+    std::copy(v.begin(), v.end(), this->begin());
+  }
+
   Vector6d(const double &d) {
     for (size_t i = 0; i < 6; i++) {
       this->operator[](i) = d;

vcgtrimesh.cpp

@@ -1,4 +1,5 @@
 #include "vcgtrimesh.hpp"
+#include <wrap/nanoply/include/nanoplyWrapper.hpp>
 #include "wrap/io_trimesh/import_obj.h"
 #include "wrap/io_trimesh/import_off.h"
 #include <filesystem>

vcgtrimesh.hpp

@@ -1,7 +1,6 @@
 #ifndef VCGTRIMESH_HPP
 #define VCGTRIMESH_HPP
 #include <vcg/complex/complex.h>
-#include <wrap/nanoply/include/nanoplyWrapper.hpp>
 
 using VertexIndex = size_t;
 class VCGTriMeshVertex;