grid_icosphere: improved performance

src/grid_icosphere.cc

@@ -2,6 +2,7 @@
 // https://schneide.wordpress.com/2016/07/15/generating-an-icosphere-in-c
 
 #include "vector3d.h"
+#include <inttypes.h>
 
 #include <utility>
 #include <limits>
@@ -18,7 +19,8 @@
 #include <map>
 #endif
 
-using Index = size_t;
+// using Index = size_t;
+using Index = uint32_t;
 using Vertex = Vector3d;
 using Triangle = std::array<Index, 3>;
 using TriangleList = std::vector<Triangle>;
@@ -58,7 +60,7 @@ namespace icosahedron
 // 12 vertices
 static VertexList vertices(12);
 
-void
+static void
 init(void)
 {
   constexpr auto pi_5 = M_PI * 0.2;
@@ -111,6 +113,9 @@ subdivide(VertexList &vertices, const TriangleList &triangles)
   Triangle mid;
 
   auto n = triangles.size();
+#ifdef UseUnorderedMap
+  lookup.reserve(1.5 * n);
+#endif
   TriangleList result(4 * n);
   for (size_t i = 0; i < n; ++i)
     {
@@ -126,20 +131,6 @@ subdivide(VertexList &vertices, const TriangleList &triangles)
   return result;
 }
 
-using IndexedMesh = std::pair<VertexList, TriangleList>;
-
-static IndexedMesh
-make_icosphere(int subdivisions)
-{
-  icosahedron::init();
-  auto vertices = icosahedron::vertices;
-  auto triangles = icosahedron::triangles;
-
-  for (int i = 0; i < subdivisions; ++i) triangles = subdivide(vertices, triangles);
-
-  return { vertices, triangles };
-}
-
 static inline void
 d_normalize(Vertex &v)
 {
@@ -185,9 +176,18 @@ size_t
 gen_icosphere_coords(int subdivisions, bool withBounds, std::vector<double> &xvals, std::vector<double> &yvals,
                      std::vector<double> &xbounds, std::vector<double> &ybounds)
 {
-  const auto mesh = make_icosphere(subdivisions);
-  const auto &vertices = mesh.first;
-  const auto &triangles = mesh.second;
+  icosahedron::init();
+  auto triangles = icosahedron::triangles;
+  auto vertices = icosahedron::vertices;
+
+  size_t numTriangles = std::pow(4, subdivisions) * 20;
+  if (numTriangles > (size_t) std::numeric_limits<Index>::max())
+    fprintf(stderr, "Too many grid cells:%zu (limit=%zu)!\n", numTriangles, (size_t) std::numeric_limits<Index>::max());
+  size_t numVerticies = std::pow(4, subdivisions) * 10 + 2;
+  vertices.reserve(numVerticies);
+  // printf("numTriangles %zu, numVerticies %zu\n", numTriangles, numVerticies);
+
+  for (int i = 0; i < subdivisions; ++i) triangles = subdivide(vertices, triangles);
 
   auto numCells = triangles.size();
   xvals.resize(numCells);
@@ -243,13 +243,13 @@ subDivdeAndVertexForEdgeTest(int subdivisions)
   size_t expectedNumVerticies = std::pow(4, subdivisions) * 10 + 2;
   if (triangles.size() != expectedNumTriangles)
     {
-      std::cout << "ERROR: unexpected number of triangles for T: std::map, expected " << expectedNumTriangles << " got "
-                << triangles.size() << std::endl;
+      std::cout << "ERROR: unexpected number of triangles, expected " << expectedNumTriangles << " got " << triangles.size()
+                << std::endl;
     }
   if (vertices.size() != expectedNumVerticies)
     {
-      std::cout << "ERROR: unexpected number of verticies for T: std::map, expected " << expectedNumVerticies << " got "
-                << vertices.size() << std::endl;
+      std::cout << "ERROR: unexpected number of verticies, expected " << expectedNumVerticies << " got " << vertices.size()
+                << std::endl;
     }
 }