From 948f90b6493da83d10e18f30b0fbb8e937e29c0b Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 24 Jun 2019 21:41:34 -0400
Subject: mostly implemented the ability to find dead bonds using topological
 properties

---
 lib/src/graph.cpp | 54 +++++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 53 insertions(+), 1 deletion(-)

(limited to 'lib/src/graph.cpp')

diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 974d3b3..bfdb952 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -33,6 +33,8 @@ graph::graph(unsigned Nx, unsigned Ny) {
   for (unsigned i = 0; i < nv; i++) {
     vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
     vertices[i].r.y = (double)(i / (Nx / 2));
+    signed f = (1 + i / (Nx / 2)) % 2 == 1 ? 1 : -1;
+    vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv};
     vertices[i].polygon = {
       {vertices[i].r.x - 1.0, vertices[i].r.y},
       {vertices[i].r.x, vertices[i].r.y - 1.0},
@@ -42,6 +44,7 @@ graph::graph(unsigned Nx, unsigned Ny) {
 
     dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
     dual_vertices[i].r.y = (double)(i / (Nx / 2));
+    dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx * (i / (Nx / 2)) + (i % (Nx / 2)), (nv + i - Nx / 2) % nv};
     dual_vertices[i].polygon = {
       {dual_vertices[i].r.x - 1.0, vertices[i].r.y},
       {dual_vertices[i].r.x, vertices[i].r.y - 1.0},
@@ -67,6 +70,29 @@ graph::graph(unsigned Nx, unsigned Ny) {
     }
   }
 
+  for (vertex& v : vertices) {
+    v.ne.resize(v.nd.size());
+  }
+
+  for (unsigned i = 0; i < edges.size(); i++) {
+    for (unsigned vi : edges[i].v) {
+      auto it1 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[0]);
+      auto it2 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[1]);
+
+      unsigned d1 = std::distance(vertices[vi].nd.begin(), it1);
+      unsigned d2 = std::distance(vertices[vi].nd.begin(), it2);
+
+      unsigned dv1 = d1 < d2 ? d1 : d2;
+      unsigned dv2 = d1 < d2 ? d2 : d1;
+
+      if (dv2 - dv1 == 1) {
+        vertices[vi].ne[dv1] = i;
+      } else {
+        vertices[vi].ne[dv2] = i;
+      }
+    }
+  }
+
 }
 
 class eulerException: public std::exception
@@ -239,12 +265,15 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
 
         if (it1 == known_vertices.end()) {
           vi1 = dual_vertices.size();
-          dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{{site->p.x, site->p.y}}});
+          dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{i1},{},{{site->p.x, site->p.y}}});
           known_vertices[t1] = vi1;
         } else {
           vi1 = it1->second;
+          dual_vertices[vi1].nd.push_back(i1);
           dual_vertices[vi1].polygon.push_back({site->p.x, site->p.y});
         }
+
+        vertices[i1].nd.push_back(vi1);
         
         if (i1 < i2 || (i1 == i2 && !self_bonded)) {
           if (i1 == i2) {
@@ -359,6 +388,29 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
     }
   }
 
+  for (vertex& v : vertices) {
+    v.ne.resize(v.nd.size());
+  }
+
+  for (unsigned i = 0; i < edges.size(); i++) {
+    for (unsigned vi : edges[i].v) {
+      auto it1 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[0]);
+      auto it2 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[1]);
+
+      unsigned d1 = std::distance(vertices[vi].nd.begin(), it1);
+      unsigned d2 = std::distance(vertices[vi].nd.begin(), it2);
+
+      unsigned dv1 = d1 < d2 ? d1 : d2;
+      unsigned dv2 = d1 < d2 ? d2 : d1;
+
+      if (dv2 - dv1 == 1) {
+        vertices[vi].ne[dv1] = i;
+      } else {
+        vertices[vi].ne[dv2] = i;
+      }
+    }
+  }
+
   jcv_diagram_free(&diagram);
 }
 
-- 
cgit v1.2.3-70-g09d2