From 6e6d0a7839ea6fd76048634c036fa824a7f5c4bf Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Mon, 12 Jul 2021 11:07:22 +0200
Subject: [PATCH] the main (heptagon) cellshape now also uses subcellshape;
 removed dirdist and get_vertices and introduced get_cellshape instead

---
 asonov.cpp        |  31 ++++---
 barriers.cpp      |   2 +-
 binary-tiling.cpp |  11 ++-
 cell.cpp          |  14 +--
 celldrawer.cpp    |  20 +++--
 commandline.cpp   |   3 +-
 crystal.cpp       |   2 +-
 euclid.cpp        |  32 ++-----
 fake.cpp          |  19 ++--
 geometry.cpp      |  21 +++--
 geometry2.cpp     |   6 --
 polygons.cpp      |   5 +-
 raycaster.cpp     |   2 +-
 reg3.cpp          | 220 +++++++++++++++++++++-------------------------
 14 files changed, 182 insertions(+), 206 deletions(-)

diff --git a/asonov.cpp b/asonov.cpp
index b4564682..63727159 100644
--- a/asonov.cpp
+++ b/asonov.cpp
@@ -110,24 +110,27 @@ EX void prepare() {
 
 EX void prepare_walls() {
 
-  cgi.cellshape.clear();
+  auto& hsh = cgi.heptshape;
+  hsh = unique_ptr<subcellshape>(new subcellshape);
+  auto& cs = hsh->faces;
+  cs.clear();
   
   auto pt = [&] (int x, int y, int z) { return asonov::tx*x/2 + asonov::ty*y/2 + asonov::tz*z/2 + C0; };
   
-  cgi.cellshape.push_back({pt(-1,-1,+1), pt(00,+1,+1), pt(+1,+1,+1)});
-  cgi.cellshape.push_back({pt(00,-1,+1), pt(+1,+1,+1), pt(+1,-1,+1)});
-  cgi.cellshape.push_back({pt(-1,+1,+1), pt(00,+1,+1), pt(-1,-1,+1)});
-  cgi.cellshape.push_back({pt(-1,-1,+1), pt(+1,+1,+1), pt(00,-1,+1)});
-  cgi.cellshape.push_back({pt(+1,-1,-1), pt(+1,00,-1), pt(+1,+1,-1), pt(+1,+1,+1), pt(+1,-1,+1)});
-  cgi.cellshape.push_back({pt(-1,+1,-1), pt(-1,+1,+1), pt(00,+1,+1), pt(+1,+1,+1), pt(+1,+1,-1)});
-  cgi.cellshape.push_back({pt(-1,-1,-1), pt(-1,00,-1), pt(+1,-1,-1)});
-  cgi.cellshape.push_back({pt(-1,00,-1), pt(-1,+1,-1), pt(+1,-1,-1)});
-  cgi.cellshape.push_back({pt(-1,+1,-1), pt(+1,00,-1), pt(+1,-1,-1)});
-  cgi.cellshape.push_back({pt(-1,+1,-1), pt(+1,+1,-1), pt(+1,00,-1)});
-  cgi.cellshape.push_back({pt(-1,+1,-1), pt(-1,00,-1), pt(-1,-1,-1), pt(-1,-1,+1), pt(-1,+1,+1)});
-  cgi.cellshape.push_back({pt(+1,-1,-1), pt(+1,-1,+1), pt(00,-1,+1), pt(-1,-1,+1), pt(-1,-1,-1)});
+  cs.push_back({pt(-1,-1,+1), pt(00,+1,+1), pt(+1,+1,+1)});
+  cs.push_back({pt(00,-1,+1), pt(+1,+1,+1), pt(+1,-1,+1)});
+  cs.push_back({pt(-1,+1,+1), pt(00,+1,+1), pt(-1,-1,+1)});
+  cs.push_back({pt(-1,-1,+1), pt(+1,+1,+1), pt(00,-1,+1)});
+  cs.push_back({pt(+1,-1,-1), pt(+1,00,-1), pt(+1,+1,-1), pt(+1,+1,+1), pt(+1,-1,+1)});
+  cs.push_back({pt(-1,+1,-1), pt(-1,+1,+1), pt(00,+1,+1), pt(+1,+1,+1), pt(+1,+1,-1)});
+  cs.push_back({pt(-1,-1,-1), pt(-1,00,-1), pt(+1,-1,-1)});
+  cs.push_back({pt(-1,00,-1), pt(-1,+1,-1), pt(+1,-1,-1)});
+  cs.push_back({pt(-1,+1,-1), pt(+1,00,-1), pt(+1,-1,-1)});
+  cs.push_back({pt(-1,+1,-1), pt(+1,+1,-1), pt(+1,00,-1)});
+  cs.push_back({pt(-1,+1,-1), pt(-1,00,-1), pt(-1,-1,-1), pt(-1,-1,+1), pt(-1,+1,+1)});
+  cs.push_back({pt(+1,-1,-1), pt(+1,-1,+1), pt(00,-1,+1), pt(-1,-1,+1), pt(-1,-1,-1)});
   
-  reg3::make_vertices_only();
+  hsh->compute_hept();
   }
 
 transmatrix coord_to_matrix(coord c, coord zero) {
diff --git a/barriers.cpp b/barriers.cpp
index fef196ac..ec8dcb4b 100644
--- a/barriers.cpp
+++ b/barriers.cpp
@@ -753,7 +753,7 @@ EX bool valid_dir(const vector<char>& ad, int j, cell *c) {
   if(!bch) return ad[j] == 1;
 
   if(ad[j] != 2) return false;
-  auto ad1 = currentmap->dirdist(c, j);
+  auto ad1 = currentmap->get_cellshape(c).dirdist[j];
   int a = 0;
   for(auto& dd: ad1) if(dd == 1) a++;
 
diff --git a/binary-tiling.cpp b/binary-tiling.cpp
index dbff0e17..9d163863 100644
--- a/binary-tiling.cpp
+++ b/binary-tiling.cpp
@@ -472,8 +472,12 @@ EX namespace bt {
       return gm * where;
       }
 
-    vector<hyperpoint> get_vertices(cell* c) override {
-      vector<hyperpoint> res;
+    subcellshape& get_cellshape(cell *c) override {
+      if(cgi.heptshape) 
+        return *cgi.heptshape;
+
+      cgi.heptshape = (std::unique_ptr<subcellshape>) (new subcellshape);
+      vector<hyperpoint>& res = cgi.heptshape->vertices_only;
       ld yy = log(2) / 2;
       auto add = [&] (hyperpoint h) { 
         res.push_back(bt::parabolic3(h[0], h[1]) * xpush0(yy*h[2]));
@@ -512,7 +516,8 @@ EX namespace bt {
           }
         default: ;
         }
-      return res;
+      cgi.heptshape->vertices_only_local = res;
+      return *cgi.heptshape;
       }
     
     ld spin_angle(cell *c, int d) override {
diff --git a/cell.cpp b/cell.cpp
index 4032a4ca..20e92765 100644
--- a/cell.cpp
+++ b/cell.cpp
@@ -45,7 +45,6 @@ struct hrmap {
     }
   virtual void draw_all();
   virtual void draw_at(cell *at, const shiftmatrix& where);
-  virtual vector<hyperpoint> get_vertices(cell*);
 
   virtual void virtualRebase(heptagon*& base, transmatrix& at) {
     printf("virtualRebase called unexpectedly\n"); 
@@ -76,13 +75,13 @@ struct hrmap {
     return currentmap->iadj(c, i);
     }
 
-  /** \brief in 3D honeycombs, returns a vector<bool> v, where v[j] iff faces i and j are adjacent */
-  virtual const vector<char>& dirdist(cell *c, int i) { throw hr_exception("dirdist called unexpectedly"); }
+  virtual subcellshape& get_cellshape(cell *c) { throw hr_exception("get_cellshape called unexpectedly"); }
 
   /** \brief in 3D honeycombs, returns a cellwalker res at cw->move(j) such that the face pointed at by cw and res share an edge */
   virtual cellwalker strafe(cellwalker cw, int j) { throw hr_exception("strafe called unexpectedly"); }
 
-  const vector<char>& dirdist(cellwalker cw) { return dirdist(cw.at, cw.spin); }
+  /** \brief in 3D honeycombs, returns a vector<bool> v, where v[j] iff faces i and j are adjacent */
+  const vector<char>& dirdist(cellwalker cw) { return get_cellshape(cw.at).dirdist[cw.spin]; }
   };
 
 /** hrmaps which are based on regular non-Euclidean 2D tilings, possibly quotient  
@@ -1287,7 +1286,8 @@ EX vector<cell*> adj_minefield_cells(cell *c) {
     }
   else if(adj_memo.count(c)) return adj_memo[c];
   else {
-    const vector<hyperpoint> vertices = currentmap->get_vertices(c);
+    auto& ss = currentmap->get_cellshape(c);
+    const vector<hyperpoint>& vertices = ss.vertices_only_local;
     manual_celllister cl;
     cl.add(c);
     for(int i=0; i<isize(cl.lst); i++) {
@@ -1295,7 +1295,9 @@ EX vector<cell*> adj_minefield_cells(cell *c) {
       bool shares = false;
       if(c != c1) {
         transmatrix T = currentmap->relative_matrix(c1->master, c->master, C0);
-        for(hyperpoint h: vertices) for(hyperpoint h2: vertices)
+        auto& ss1 = currentmap->get_cellshape(c1);
+        auto& vertices1 = ss1.vertices_only_local;
+        for(hyperpoint h: vertices) for(hyperpoint h2: vertices1)
           if(hdist(h, T * h2) < 1e-6) shares = true;
         if(shares) res.push_back(c1);
         }
diff --git a/celldrawer.cpp b/celldrawer.cpp
index 46865f2b..4dab9017 100644
--- a/celldrawer.cpp
+++ b/celldrawer.cpp
@@ -864,12 +864,14 @@ void celldrawer::draw_grid() {
       if(!bt::in() && c->move(t) < c) continue;
       dynamicval<color_t> g(poly_outline, gridcolor(c, c->move(t)));          
       if(fat_edges && reg3::in()) {
-        for(int i=0; i<S7; i++) if(c < c->move(i) || fake::split()) {
-          for(int j=0; j<cgi.face; j++) {
-            int jj = j == cgi.face-1 ? 0 : j+1;
-            int jjj = jj == cgi.face-1 ? 0 : jj+1;
-            hyperpoint a = cgi.cellshape[i][j];
-            hyperpoint b = cgi.cellshape[i][jj];
+        auto& ss = currentmap->get_cellshape(c);
+        for(int i=0; i<c->type; i++) if(c < c->move(i) || fake::split()) {
+          int face = isize(ss.faces[i]);
+          for(int j=0; j<face; j++) {
+            int jj = j == face-1 ? 0 : j+1;
+            int jjj = jj == face-1 ? 0 : jj+1;
+            hyperpoint a = ss.faces[i][j];
+            hyperpoint b = ss.faces[i][jj];
             if(cgflags & qIDEAL) {
               ld mm = cgi.ultra_mirror_part;
               if((cgflags & qULTRA) && !reg3::ultra_mirror_in())
@@ -879,9 +881,9 @@ void celldrawer::draw_grid() {
             gridline(V, a, b, gridcolor(c, c->move(t)), prec);            
 
             if(reg3::ultra_mirror_in()) {
-              hyperpoint a = cgi.cellshape[i][j];
-              hyperpoint b = cgi.cellshape[i][jj];
-              hyperpoint d = cgi.cellshape[i][jjj];
+              hyperpoint a = ss.faces[i][j];
+              hyperpoint b = ss.faces[i][jj];
+              hyperpoint d = ss.faces[i][jjj];
               auto& mm = cgi.ultra_mirror_part;
               tie(a, d) = make_pair(normalize(lerp(a, b, mm)), normalize(lerp(d, b, mm)));
               gridline(V, a, d, stdgridcolor, prec);
diff --git a/commandline.cpp b/commandline.cpp
index c801ba87..6f619906 100644
--- a/commandline.cpp
+++ b/commandline.cpp
@@ -250,7 +250,8 @@ int arg::readCommon() {
     }
   else if(argis("-face-vertex")) {
     PHASE(3);  start_game();
-    View = cspin(0, 2, M_PI/2) * spintox(cgi.vertices_only[0]);
+    auto &ss = cgi.get_cellshape(cwt.at);
+    View = cspin(0, 2, M_PI/2) * spintox(ss.vertices_only_local[0]);
     }
   else if(argis("-face-face")) {
     PHASE(3);  start_game();
diff --git a/crystal.cpp b/crystal.cpp
index 1cc745c6..d0d3c01e 100644
--- a/crystal.cpp
+++ b/crystal.cpp
@@ -647,7 +647,7 @@ struct hrmap_crystal : hrmap_standard {
       if(go > 1e-2) continue;
 
       for(int s=0; s<S7; s++) 
-        if(cgi.dirs_adjacent[d][s])
+        if(cgi.heptshape->dirdist[d][s] == 1)
           for(int t=0; t<S7; t++) 
             if(st1[t] == st[s]) {
               if(hdist(U * tC0(cgi.adjmoves[t]), tC0(cgi.adjmoves[s])) > gdist + .1)
diff --git a/euclid.cpp b/euclid.cpp
index d7e71d31..ad43108e 100644
--- a/euclid.cpp
+++ b/euclid.cpp
@@ -292,28 +292,8 @@ EX namespace euc {
       return eumove(d);
       }
     
-    vector<hyperpoint> get_vertices(cell* c) override {
-      vector<hyperpoint> res;
-      if(S7 < 14)
-        for(ld a: {-.5,.5}) for(ld b: {-.5,.5}) for(ld c: {-.5, .5}) res.push_back(hpxy3(a,b,c));
-      if(S7 == 12) {
-        res.push_back(hpxy3(1,0,0));
-        res.push_back(hpxy3(-1,0,0));
-        res.push_back(hpxy3(0,1,0));
-        res.push_back(hpxy3(0,-1,0));
-        res.push_back(hpxy3(0,0,1));
-        res.push_back(hpxy3(0,0,-1));
-        }
-      if(S7 == 14) {
-        for(ld a: {-1.,-.5,0.,.5,1.})
-        for(ld b: {-1.,-.5,0.,.5,1.})
-        for(ld c: {-1.,-.5,0.,.5,1.})
-          if(a == 0 || b == 0 || c == 0)
-          if(a == .5 || a == -.5 || b == .5 || b == -.5 || c == .5 || c == -.5)
-          if(a == 1 || a == -1 || b == 1 || b == -1 || c == 1 || c == -1)
-            res.push_back(hpxy3(a,b,c));
-        }
-      return res;
+    subcellshape& get_cellshape(cell* c) override {
+      return *cgi.heptshape;
       }
     };
   
@@ -1300,7 +1280,10 @@ EX void generate() {
 
   auto v = euc::get_shifttable();
   
-  auto& cs = cgi.cellshape;
+  auto& hsh = cgi.heptshape;
+  hsh = unique_ptr<subcellshape>(new subcellshape);
+
+  auto& cs = hsh->faces;
 
   cgi.loop = 4;
   cgi.schmid = 3;
@@ -1365,9 +1348,8 @@ EX void generate() {
       }
     }
   
-  reg3::make_vertices_only();
+  hsh->compute_hept();
   #endif
-
   }
 
 /** @brief returns true if the current geometry is based on this module 
diff --git a/fake.cpp b/fake.cpp
index 32957929..2e8eef73 100644
--- a/fake.cpp
+++ b/fake.cpp
@@ -332,11 +332,11 @@ EX vector<vector<hyperpoint>> befake(const vector<vector<hyperpoint>>& v) {
 EX ld compute_around(bool setup) {
   auto &ucgi = *underlying_cgip;
   
-  auto fcs = befake(ucgi.cellshape);
+  auto fcs = befake(ucgi.heptshape->faces);
   
   if(setup) {
-    cgi.cellshape = fcs;
-    cgi.vertices_only = befake(ucgi.vertices_only);
+    cgi.heptshape->faces = fcs;
+    cgi.heptshape->compute_hept();
     }
   
   hyperpoint h = Hypc;
@@ -394,11 +394,10 @@ EX void generate() {
   cgi.face = ucgi.face;
   cgi.schmid = ucgi.schmid;
 
-  for(int a=0; a<16; a++)
-  for(int b=0; b<16; b++) {
-    cgi.dirs_adjacent[a][b] = ucgi.dirs_adjacent[a][b];
-    cgi.next_dir[a][b] = ucgi.next_dir[a][b];
-    }
+  auto& hsh = cgi.heptshape;
+  hsh = unique_ptr<subcellshape>(new subcellshape);
+  
+  *hsh = *ucgi.heptshape;
 
   for(int b=0; b<12; b++)
     cgi.spins[b] = ucgi.spins[b];
@@ -487,7 +486,7 @@ EX void compute_scale() {
     }
   else if(euclid) scale = 1;
   else if(have_ideal) {
-    hyperpoint h0 = underlying_cgip->cellshape[0][0];
+    hyperpoint h0 = underlying_cgip->heptshape->faces[0][0];
     auto s = kleinize(h0);
     ld d = hypot_d(LDIM, s);
     scale = 1/d;
@@ -517,7 +516,7 @@ EX void compute_scale() {
     
     /* ultra a bit earlier */    
     if(underlying == gRhombic3 || underlying == gBitrunc3) {
-      auto fcs = befake(underlying_cgip->cellshape[0][0]);
+      auto fcs = befake(underlying_cgip->heptshape->faces[0][0]);
       set_flag(ginf[gFake].flags, qULTRA, material(fcs) < 0);
       }
     }
diff --git a/geometry.cpp b/geometry.cpp
index e2fd6fbd..7fdaba9e 100644
--- a/geometry.cpp
+++ b/geometry.cpp
@@ -98,7 +98,7 @@ struct gi_extension {
   virtual ~gi_extension() {}
   };
 
-/** for subdivided 3D cells */
+/** both for 'heptagon' 3D cells and subdivided 3D cells */
 struct subcellshape {
   vector<vector<hyperpoint>> faces;
   vector<vector<hyperpoint>> faces_local;
@@ -109,6 +109,17 @@ struct subcellshape {
   hyperpoint cellcenter;
   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;  
+
+  /** compute all the properties based on `faces`, for the main heptagon cellshape */
+  void compute_hept();
+
+  /** compute all the properties based on `faces`, for subcells */
+  void compute_sub();
+
+  /** common part of compute_hept and compute_sub */
+  void compute_common();
   };
 
 /** basic geometry parameters */
@@ -143,24 +154,18 @@ struct geometry_information {
 
   int loop, face, schmid;
 
-  vector<vector<hyperpoint>> cellshape;
-  vector<hyperpoint> vertices_only;
-  
   transmatrix spins[32], adjmoves[32];
   
+  unique_ptr<struct subcellshape> heptshape;  
   vector<struct subcellshape> subshapes;
 
   ld adjcheck;
   ld strafedist;
-  vector<vector<char>> dirdist;
 
   ld ultra_mirror_dist, ultra_material_part, ultra_mirror_part;
   
   vector<transmatrix> ultra_mirrors;  
 
-  /** \brief for adjacent directions a,b, next_dir[a][b] is the next direction adjacent to a, in (counter?)clockwise order from b */
-  int next_dir[32][32];
-
   vector<pair<string, string> > rels;
   int xp_order, r_order, rx_order;
 
diff --git a/geometry2.cpp b/geometry2.cpp
index 8b324b13..49b3f4e0 100644
--- a/geometry2.cpp
+++ b/geometry2.cpp
@@ -734,12 +734,6 @@ EX hyperpoint get_warp_corner(cell *c, int cid) {
   return ddspin(c,cid,M_PI/S7) * xpush0(cgi.tessf/2);
   }
 
-vector<hyperpoint> hrmap::get_vertices(cell* c) {
-  vector<hyperpoint> res;
-  for(int i=0; i<c->type; i++) res.push_back(get_corner_position(c, i, 3));
-  return res;
-  }
-
 EX map<cell*, map<cell*, vector<transmatrix>>> brm_structure;
 
 EX void generate_brm(cell *c1) {
diff --git a/polygons.cpp b/polygons.cpp
index 587ac6f7..1e69d080 100644
--- a/polygons.cpp
+++ b/polygons.cpp
@@ -1047,8 +1047,9 @@ void geometry_information::create_wall3d() {
     }
 
   if(euc::in() || reg3::in() || asonov::in()) {
-    for(int w=0; w<isize(cgi.cellshape); w++)
-      make_wall(w, cgi.cellshape[w]);
+    auto& faces = cgi.heptshape->faces;
+    for(int w=0; w<isize(faces); w++)
+      make_wall(w, faces[w]);
     }
 
   if(geometry == gSol) {
diff --git a/raycaster.cpp b/raycaster.cpp
index debfb0c5..ae99ae8e 100644
--- a/raycaster.cpp
+++ b/raycaster.cpp
@@ -701,7 +701,7 @@ void enable_raycaster() {
             "}\n";
 
       if(hyperbolic && reg3::ultra_mirror_in()) {
-        fmain += "for(int i="+its(S7*2)+"; i<"+its(S7*2+isize(cgi.vertices_only))+"; i++) {\n";
+        fmain += "for(int i="+its(S7*2)+"; i<"+its(S7*2+isize(cgi.heptshape->vertices_only))+"; i++) {\n";
         fmain += "mat4 uMi = " + getM("i") + ";";
         fmain +=
           "    mediump float v = ((position - uMi * position)[3] / (uMi * tangent - tangent)[3]);\n"
diff --git a/reg3.cpp b/reg3.cpp
index 2e0753b4..d296324e 100644
--- a/reg3.cpp
+++ b/reg3.cpp
@@ -13,6 +13,65 @@
 namespace hr {
 #if MAXMDIM >= 4
 
+void subcellshape::compute_common() {
+  reg3::make_vertices_only(vertices_only, faces);
+
+  faces_local = faces;
+  for(auto& face: faces_local) for(auto& v: face) v = from_cellcenter * v;
+
+  vertices_only_local = vertices_only;
+  for(auto& v: vertices_only_local) v = from_cellcenter * v;
+
+  int N = isize(faces);
+
+  dirdist.resize(N);
+  for(int i=0; i<N; i++) {
+    auto& da = dirdist[i];
+    da.resize(N, false);
+    set<unsigned> cface;
+    for(auto& v: faces[i]) cface.insert(bucketer(v));
+    for(int j=0; j<N; j++) {
+      int mutual = 0;
+      for(auto& w: faces[j]) if(cface.count(bucketer(w))) mutual++;
+      da[j] = i == j ? 0 : mutual == 2 ? 1 : INFD;
+      }
+    }
+  floyd_warshall(dirdist);
+
+  next_dir.resize(N);
+  for(int a=0; a<N; a++) next_dir[a].resize(N);
+  
+  for(int a=0; a<N; a++)
+  for(int b=0; b<N; b++)
+    if(dirdist[a][b] == 1) 
+      for(int c=0; c<N; c++)
+        if(dirdist[a][c] == 1 && dirdist[b][c] == 1) {
+          transmatrix t = build_matrix(tC0(cgi.adjmoves[a]), tC0(cgi.adjmoves[b]), tC0(cgi.adjmoves[c]), C0);
+          if(det(t) > 0) next_dir[a][b] = c;
+          }  
+  }
+
+void subcellshape::compute_hept() {
+  cellcenter = C0;
+  to_cellcenter = Id;
+  from_cellcenter = Id;
+  compute_common();
+  }
+
+void subcellshape::compute_sub() {
+  hyperpoint gres = Hypc;
+  for(auto& face: faces) {
+    hyperpoint res = Hypc;
+    for(auto& vertex: face) 
+      res += vertex;
+    face_centers.push_back(normalize(res));
+    gres += res;
+    }
+  cellcenter = normalize(gres);
+  to_cellcenter = rgpushxto0(cellcenter);
+  from_cellcenter = gpushxto0(cellcenter);
+  }
+
 /** \brief regular three-dimensional tessellations */
 EX namespace reg3 {
 
@@ -47,14 +106,15 @@ EX namespace reg3 {
 
     if(cgflags & qULTRA) {
     
-      for(auto& v: cgi.vertices_only) {
+      for(auto& v: cgi.heptshape->vertices_only) {
       
         hyperpoint nei;
+        auto& faces = cgi.heptshape->faces;
       
-        for(int i=0; i<isize(cgi.cellshape); i++)
-        for(int j=0; j<isize(cgi.cellshape[i]); j++)
-          if(sqhypot_d(WDIM, cgi.cellshape[i][j]-v) < 1e-6)
-            nei = cgi.cellshape[i][j?j-1:j+1]; 
+        for(int i=0; i<isize(faces); i++)
+        for(int j=0; j<isize(faces[i]); j++)
+          if(sqhypot_d(WDIM, faces[i][j]-v) < 1e-6)
+            nei = faces[i][j?j-1:j+1]; 
             
         transmatrix T = spintox(v);
         hyperpoint a = T * v;
@@ -106,10 +166,6 @@ EX namespace reg3 {
       }
     }
 
-  EX void make_vertices_only() {
-    make_vertices_only(cgi.vertices_only, cgi.cellshape);
-    }
-
   EX void generate() {
 
     if(fake::in()) {
@@ -117,12 +173,14 @@ EX namespace reg3 {
       return;
       }
   
+    auto& hsh = cgi.heptshape;
+    hsh = unique_ptr<subcellshape>(new subcellshape);
+    
     int& loop = cgi.loop;
     int& face = cgi.face;
     auto& spins = cgi.spins;
-    auto& cellshape = cgi.cellshape;
+    auto& cellshape = hsh->faces;
     auto& adjcheck = cgi.adjcheck;
-    auto& dirdist = cgi.dirdist;
   
     int& mid = cgi.schmid;
     mid = 3;
@@ -255,21 +313,6 @@ EX namespace reg3 {
 
     adjcheck = hdist(tC0(cgi.adjmoves[0]), tC0(cgi.adjmoves[1])) * 1.0001;
 
-    int numedges = 0;
-    dirdist.resize(S7);
-    for(int a=0; a<S7; a++) {
-      dirdist[a].resize(S7);
-      for(int b=0; b<S7; b++) {
-        dirdist[a][b] =
-          a == b ? 0 :
-          hdist(tC0(cgi.adjmoves[a]), tC0(cgi.adjmoves[b])) < adjcheck ? 1 :
-          INFD;
-        if(dirdist[a][b] == 1) numedges++;
-        }
-      }
-    floyd_warshall(dirdist);
-    DEBB(DF_GEOM, ("numedges = ", numedges));
-    
     if(loop == 4) cgi.strafedist = adjcheck;
     else cgi.strafedist = hdist(cgi.adjmoves[0] * C0, cgi.adjmoves[1] * C0);
     
@@ -280,18 +323,9 @@ EX namespace reg3 {
       for(auto& vv: cellshape) for(auto& v: vv) v = T * v;
       }
 
-    make_vertices_only();
+    hsh->compute_hept();
     compute_ultra();
     
-    for(int a=0; a<S7; a++)
-    for(int b=0; b<S7; b++)
-      if(cgi.dirdist[a][b] == 1) 
-        for(int c=0; c<S7; c++)
-          if(cgi.dirdist[a][c] == 1 && cgi.dirdist[b][c] == 1) {
-            transmatrix t = build_matrix(tC0(cgi.adjmoves[a]), tC0(cgi.adjmoves[b]), tC0(cgi.adjmoves[c]), C0);
-            if(det(t) > 1e-3) cgi.next_dir[a][b] = c;
-            }  
-    
     generate_subcells();
     }
 
@@ -299,14 +333,14 @@ EX namespace reg3 {
     if(S7 != 6) throw hr_exception("generate_plain_subcubes but no cubes");
     auto& ssh = cgi.subshapes;
     const int sub = subcube_count;
-    auto vx = abs(cgi.cellshape[0][0][0]);
-    auto vz = abs(cgi.cellshape[0][0][3]);
+    auto vx = abs(cgi.heptshape->faces[0][0][0]);
+    auto vz = abs(cgi.heptshape->faces[0][0][3]);
     for(int x=1-sub; x<sub; x+=2)
     for(int y=1-sub; y<sub; y+=2)
     for(int z=1-sub; z<sub; z+=2) {
       ssh.emplace_back();
       auto &ss = ssh.back();
-      ss.faces = cgi.cellshape;
+      ss.faces = cgi.heptshape->faces;
       for(auto& face: ss.faces) for(auto& v: face) {
         v[0] += vx * x;
         v[1] += vx * y;
@@ -318,7 +352,7 @@ EX namespace reg3 {
   
   EX void generate_coxeter(flagtype f) {
     auto& ssh = cgi.subshapes;
-    for(auto& fac: cgi.cellshape) {
+    for(auto& fac: cgi.heptshape->faces) {
       hyperpoint facectr = Hypc;
       vector<hyperpoint> ring;
       hyperpoint last = fac.back();
@@ -384,8 +418,8 @@ EX namespace reg3 {
     if(S7 != 6) throw hr_exception("generate_plain_subcubes but no cubes");
     const int sub = subcube_count;
     if(1) {
-      auto vx = abs(cgi.cellshape[0][0][0]);
-      auto vz = abs(cgi.cellshape[0][0][3]);
+      auto vx = abs(cgi.heptshape->faces[0][0][0]);
+      auto vz = abs(cgi.heptshape->faces[0][0][3]);
       auto step = hdist0(tC0(cgi.adjmoves[0]));
       array<int, 3> co;
       int s = bch ? 1 : 2;
@@ -596,7 +630,7 @@ EX namespace reg3 {
     
       case eVariation::pure: {
         cgi.subshapes.emplace_back();
-        cgi.subshapes[0].faces = cgi.cellshape;
+        cgi.subshapes[0].faces = cgi.heptshape->faces;
         break;
         }
     
@@ -604,39 +638,7 @@ EX namespace reg3 {
         throw hr_exception("unknown variation in generate_subcells");
       }
 
-    for(auto& ss: cgi.subshapes) {
-      make_vertices_only(ss.vertices_only, ss.faces);
-      hyperpoint gres = Hypc;
-      for(auto& face: ss.faces) {
-        hyperpoint res = Hypc;
-        for(auto& vertex: face) 
-          res += vertex;
-        ss.face_centers.push_back(normalize(res));
-        gres += res;
-        }
-      ss.cellcenter = normalize(gres);
-      ss.to_cellcenter = rgpushxto0(ss.cellcenter);
-      ss.from_cellcenter = gpushxto0(ss.cellcenter);
-      ss.faces_local = ss.faces;
-      for(auto& face: ss.faces_local) for(auto& v: face) v = ss.from_cellcenter * v;
-      ss.vertices_only_local = ss.vertices_only;
-      for(auto& v: ss.vertices_only_local) v = ss.from_cellcenter * v;
-      
-      int N = isize(ss.faces);
-      ss.dirdist.resize(N);
-      for(int i=0; i<N; i++) {
-        auto& da = ss.dirdist[i];
-        da.resize(N, false);
-        set<unsigned> cface;
-        for(auto& v: ss.faces[i]) cface.insert(bucketer(v));
-        for(int j=0; j<N; j++) {
-          int mutual = 0;
-          for(auto& w: ss.faces[j]) if(cface.count(bucketer(w))) mutual++;
-          da[j] = i == j ? 0 : mutual == 2 ? 1 : INFD;
-          }
-        }
-      floyd_warshall(ss.dirdist);
-      }
+    for(auto& ss: cgi.subshapes) ss.compute_sub();
 
     println(hlog, "subcells generated = ", isize(cgi.subshapes));
     }
@@ -680,10 +682,10 @@ EX namespace reg3 {
     void initialize(int cell_count);
     vector<cell*>& allcells() override { return acells; }
 
-    vector<hyperpoint> get_vertices(cell* c) override { 
-      if(PURE) return cgi.vertices_only; 
+    subcellshape& get_cellshape(cell *c) override {
+      if(PURE) return *cgi.heptshape ;
       int id = local_id.at(c).second;
-      return cgi.subshapes[id].vertices_only_local;
+      return cgi.subshapes[id];
       }
 
     transmatrix master_relative(cell *c, bool get_inverse) override {
@@ -699,11 +701,6 @@ EX namespace reg3 {
     
     transmatrix ray_iadj(cell *c, int i) override;
 
-    const vector<char>& dirdist(cell *c, int i) override { 
-      int id = local_id.at(c).second;
-      return cgi.subshapes[id].dirdist[i];
-      }
-
     cellwalker strafe(cellwalker cw, int j) override {
       int id = local_id.at(cw.at).first;
       return cellwalker(cw.at->cmove(j), strafe_data[id][j][cw.spin]);
@@ -760,7 +757,7 @@ EX namespace reg3 {
       va.emplace_back(vertex_adjacency_info{a, Id, {}});
       
       set<unsigned> buckets;
-      for(auto& v: cgi.vertices_only) buckets.insert(bucketer(v));
+      for(auto& v: cgi.heptshape->vertices_only) buckets.insert(bucketer(v));
 
       if(cgflags & qIDEAL) {
         for(int d=0; d<S7; d++) {
@@ -777,7 +774,7 @@ EX namespace reg3 {
           if(found) continue;
 
           bool found_va = false;
-          for(auto& w: cgi.vertices_only)
+          for(auto& w: cgi.heptshape->vertices_only)
             if(buckets.count(bucketer(T*w)))
               found_va = true;
           if(!found_va) continue;
@@ -1031,8 +1028,9 @@ EX namespace reg3 {
     void make_plane(cellwalker cw) {
       if(plane.count(cw)) return;
       plane.insert(cw);
-      for(int i=0; i<S7; i++)
-        if(cgi.dirdist[i][cw.spin] == 1)
+      auto& ss = get_cellshape(cw.at);
+      for(int i=0; i<cw.at->type; i++)
+        if(ss.dirdist[i][cw.spin] == 1)
           make_plane(strafe(cw, i));
       }
     
@@ -1131,19 +1129,12 @@ EX namespace reg3 {
   
     void build_reps() {
       // start_game();
-      for(int a=0; a<12; a++)
-      for(int b=0; b<12; b++)
-        if(cgi.dirdist[a][b] == 1) 
-          for(int c=0; c<12; c++)
-            if(cgi.dirdist[a][c] == 1 && cgi.dirdist[b][c] == 1) {
-              transmatrix t = build_matrix(tC0(cgi.adjmoves[a]), tC0(cgi.adjmoves[b]), tC0(cgi.adjmoves[c]), C0);
-              if(det(t) > 0) cgi.next_dir[a][b] = c;
-              }  
+      auto& hsh = cgi.heptshape;
       
       set<coord> boundaries;
       
       for(int a=0; a<12; a++)
-      for(int b=0; b<12; b++) if(cgi.dirdist[a][b] == 1) {
+      for(int b=0; b<12; b++) if(hsh->dirdist[a][b] == 1) {
         coord res = crystal::c0;
         int sa = a, sb = b;
         do {
@@ -1152,7 +1143,7 @@ EX namespace reg3 {
           sa = flip(sa);
           sb = flip(sb);
           swap(sa, sb);
-          sb = cgi.next_dir[sa][sb];
+          sb = hsh->next_dir[sa][sb];
           // sb = next_dirsa][sb];
           }
         while(a != sa || b != sb);
@@ -1518,12 +1509,8 @@ EX namespace reg3 {
       return T;
       }
     
-    vector<hyperpoint> get_vertices(cell* c) override {
-      return cgi.vertices_only;
-      }
-
-    const vector<char>& dirdist(cell *c, int i) override { 
-      return cgi.dirdist[i];
+    subcellshape& get_cellshape(cell *c) override {
+      return *cgi.heptshape;
       }
 
     cellwalker strafe(cellwalker cw, int j) override {
@@ -1964,10 +1951,10 @@ EX namespace reg3 {
       return quotient_map->tmatrices_cell[aid][d];
       }     
     
-    vector<hyperpoint> get_vertices(cell* c) override {
-      if(PURE) return cgi.vertices_only;
+    subcellshape& get_cellshape(cell *c) override {
+      if(PURE) return *cgi.heptshape;
       int aid = cell_id.at(c);
-      return quotient_map->get_vertices(quotient_map->acells[aid]);
+      return quotient_map->get_cellshape(quotient_map->acells[aid]);
       }
     
     map<cell*, int> cell_id;
@@ -2006,12 +1993,6 @@ EX namespace reg3 {
       return quotient_map->ray_iadj(quotient_map->acells[aid], i);
       }
 
-    const vector<char>& dirdist(cell *c, int i) override { 
-      if(PURE) return cgi.dirdist[i];
-      int aid = cell_id.at(c);
-      return quotient_map->dirdist(quotient_map->acells[aid], i);
-      }
-
     cellwalker strafe(cellwalker cw, int j) override {
 
       hyperpoint hfront = tC0(cgi.adjmoves[cw.spin]);
@@ -2398,7 +2379,8 @@ EX void construct_relations() {
   formulas.push_back("");
 
   all.push_back(Id);
-  hyperpoint v = cgi.cellshape[0][0];
+  auto& faces = cgi.heptshape->faces;
+  hyperpoint v = faces[0][0];
   auto add = [&] (transmatrix T) {
     for(int i=0; i<isize(all); i++) if(eqmatrix(all[i], T)) return i;
     int S = isize(all);
@@ -2406,14 +2388,14 @@ EX void construct_relations() {
     return S;
     };
   
-  println(hlog, cgi.cellshape);
+  println(hlog, faces);
 
-  println(hlog, "cellshape = ", isize(cgi.cellshape));
+  println(hlog, "cellshape = ", isize(faces));
   bool ok = true;
   int last_i = -1;
-  for(auto& v: cgi.cellshape) for(hyperpoint h: v) {
+  for(auto& v: faces) for(hyperpoint h: v) {
     int i = 0, j = 0;
-    for(auto& uv: cgi.cellshape) for(hyperpoint u: uv) {
+    for(auto& uv: faces) for(hyperpoint u: uv) {
       if(hdist(h, cgi.full_X*u) < 5e-2) i++;
       if(hdist(h, cgi.full_R*u) < 5e-2) j++;
       }
@@ -2427,7 +2409,7 @@ EX void construct_relations() {
   
   auto work = [&] (transmatrix T, int p, char c) {
     if(hdist0(tC0(T)) > 5) return;
-    for(auto& hv: cgi.cellshape) for(hyperpoint h: hv) if(hdist(T * h, v) < 1e-4) goto ok;
+    for(auto& hv: faces) for(hyperpoint h: hv) if(hdist(T * h, v) < 1e-4) goto ok;
     return;
     ok:
     int id = add(T);