From c4e85c16d3f43c8b95966f17601e5cc2859a2ac6 Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Mon, 12 Jul 2021 20:54:16 +0200
Subject: [PATCH] full cellshape computed in hybrid geometries

---
 geometry.cpp     |   4 +-
 graph.cpp        | 100 ++++++++++++++++++++++++++++-------------------
 nonisotropic.cpp |   5 +++
 3 files changed, 67 insertions(+), 42 deletions(-)

diff --git a/geometry.cpp b/geometry.cpp
index ca6404cb..d9640127 100644
--- a/geometry.cpp
+++ b/geometry.cpp
@@ -117,7 +117,9 @@ struct subcellshape {
   transmatrix to_cellcenter;
   transmatrix from_cellcenter;
   /** \brief for adjacent directions a,b, next_dir[a][b] is the next direction adjacent to a, in (counter?)clockwise order from b */
-  vector<vector<char>> next_dir;  
+  vector<vector<char>> next_dir;
+  /** useful in product geometries */
+  vector<hyperpoint> walltester;
 
   /** compute all the properties based on `faces`, for the main heptagon cellshape */
   void compute_hept();
diff --git a/graph.cpp b/graph.cpp
index a292996c..3b077c89 100644
--- a/graph.cpp
+++ b/graph.cpp
@@ -3963,6 +3963,57 @@ EX void gridline(const shiftmatrix& V, const hyperpoint h1, const hyperpoint h2,
   gridline(V, h1, V, h2, col, prec);
   }
 
+EX subcellshape& generate_subcellshape_if_needed(cell *c, int id) {
+  if(isize(cgi.subshapes) <= id) cgi.subshapes.resize(id+1);
+  
+  auto& ss = cgi.subshapes[id];
+  if(!ss.faces.empty()) return ss;
+  
+  cell *c1 = hybri ? hybrid::get_where(c).first : c;
+  
+  if(prod || WDIM == 2) for(int i=0; i<c1->type; i++) {
+    hyperpoint w;
+    auto f = [&] { 
+      /* mirror image of C0 in the axis h1-h2 */
+      hyperpoint h1 = get_corner_position(c1, i);
+      hyperpoint h2 = get_corner_position(c1, i+1);
+      transmatrix T = gpushxto0(h1);
+      T = spintox(T * h2) * T;
+      w = T * C0;
+      w[1] = -w[1];
+      w = iso_inverse(T) * w;
+      };
+    if(prod) PIU(f());
+    else f();
+    ss.walltester.push_back(w);
+    }
+
+  for(int i=0; i<c1->type; i++)
+    ss.faces.push_back({hybrid::get_corner(c1, i, 0, -1), hybrid::get_corner(c1, i, 0, +1), hybrid::get_corner(c1, i, 1, +1), hybrid::get_corner(c1, i, 1, -1)});
+
+  for(int a: {0,1}) {
+    vector<hyperpoint> l;
+    int z = a ? 1 : -1;
+    hyperpoint ctr = zpush0(z * cgi.plevel/2);
+    for(int i=0; i<c1->type; i++)
+      if(prod || WDIM == 2)
+        l.push_back(hybrid::get_corner(c1, i, 0, z));
+      else {
+        l.push_back(ctr);
+        l.push_back(hybrid::get_corner(c1, i, 0, z));
+        l.push_back(hybrid::get_corner(c1, i+1, 1, z));
+        l.push_back(ctr);
+        l.push_back(hybrid::get_corner(c1, i, 1, z));
+        l.push_back(hybrid::get_corner(c1, i, 0, z));
+        }
+    if(a == 0) std::reverse(l.begin()+1, l.end());
+    ss.faces.push_back(l);
+    }
+  
+  ss.compute_hept();
+  return ss;
+  }
+
 int hrmap::wall_offset(cell *c) {
   int id = currentmap->full_shvid(c);
 
@@ -3973,51 +4024,18 @@ int hrmap::wall_offset(cell *c) {
   int &wo = wop.first;
   if(!wop.second) wop.second = c;
   if(wo == -1) {
-    cell *c1 = hybri ? hybrid::get_where(c).first : c;
+    auto& ss = generate_subcellshape_if_needed(c, id);
     wo = isize(cgi.shWall3D);
-    int won = wo + c->type + (WDIM == 2 ? 2 : 0);
+
     if(!cgi.wallstart.empty()) cgi.wallstart.pop_back();
-    cgi.reserve_wall3d(won);
+    cgi.reserve_wall3d(wo + isize(ss.faces));
+
     
-    if(prod || WDIM == 2) for(int i=0; i<c1->type; i++) {
-      hyperpoint w;
-      auto f = [&] { 
-        /* mirror image of C0 in the axis h1-h2 */
-        hyperpoint h1 = get_corner_position(c1, i);
-        hyperpoint h2 = get_corner_position(c1, i+1);
-        transmatrix T = gpushxto0(h1);
-        T = spintox(T * h2) * T;
-        w = T * C0;
-        w[1] = -w[1];
-        w = iso_inverse(T) * w;
-        };
-      if(prod) PIU(f());
-      else f();
-      cgi.walltester[wo + i] = w;
-      } 
-
-    for(int i=0; i<c1->type; i++)
-     cgi.make_wall(wo + i, {hybrid::get_corner(c1, i, 0, -1), hybrid::get_corner(c1, i, 0, +1), hybrid::get_corner(c1, i, 1, +1), hybrid::get_corner(c1, i, 1, -1)});
-
-    for(int a: {0,1}) {
-      vector<hyperpoint> l;
-      int z = a ? 1 : -1;
-      hyperpoint ctr = zpush0(z * cgi.plevel/2);
-      for(int i=0; i<c1->type; i++)
-        if(prod || WDIM == 2)
-          l.push_back(hybrid::get_corner(c1, i, 0, z));
-        else {
-          l.push_back(ctr);
-          l.push_back(hybrid::get_corner(c1, i, 0, z));
-          l.push_back(hybrid::get_corner(c1, i+1, 1, z));
-          l.push_back(ctr);
-          l.push_back(hybrid::get_corner(c1, i, 1, z));
-          l.push_back(hybrid::get_corner(c1, i, 0, z));
-          }
-      if(a == 0) std::reverse(l.begin()+1, l.end());
-      cgi.make_wall(won-2+a, l);
+    for(int i=0; i<isize(ss.faces); i++) {
+      cgi.make_wall(wo + i, ss.faces[i]);
+      cgi.walltester[wo + i] = ss.walltester[i];
       }
-
+    
     cgi.wallstart.push_back(isize(cgi.raywall));
     cgi.compute_cornerbonus();
     cgi.extra_vertices();
diff --git a/nonisotropic.cpp b/nonisotropic.cpp
index 78632c9d..a31d7548 100644
--- a/nonisotropic.cpp
+++ b/nonisotropic.cpp
@@ -1306,6 +1306,11 @@ EX namespace hybrid {
   
     virtual transmatrix spin_to(cell *c, int d, ld bonus) override { if(d >= c->type-2) return Id; c = get_where(c).first; return in_underlying([&] { return currentmap->spin_to(c, d, bonus); }); }
     virtual transmatrix spin_from(cell *c, int d, ld bonus) override { if(d >= c->type-2) return Id; c = get_where(c).first; return in_underlying([&] { return currentmap->spin_from(c, d, bonus); }); }
+
+    subcellshape& get_cellshape(cell *c) override {      
+      int id = full_shvid(c);
+      return generate_subcellshape_if_needed(c, id);      
+      }
     };
   
   hrmap_hybrid* hmap() { return (hrmap_hybrid*) currentmap; }