partially fixed problems that arise in small systems by reworking the way that boundary edges are identified

7 files changed, 87 insertions, 64 deletions

diff --git a/lib/include/array_hash.hpp b/lib/include/array_hash.hpp
new file mode 100644
index 0000000..3b35e4a
--- /dev/null
+++ b/lib/include/array_hash.hpp
@@ -0,0 +1,27 @@
+
+#pragma once
+
+#include <array>
+#include <unordered_map>
+
+namespace std
+{
+    template<typename T, size_t N>
+    struct hash<array<T, N> >
+    {
+        typedef array<T, N> argument_type;
+        typedef size_t result_type;
+
+        result_type operator()(const argument_type& a) const
+        {
+            hash<T> hasher;
+            result_type h = 0;
+            for (result_type i = 0; i < N; ++i)
+            {
+                h = h * 31 + hasher(a[i]);
+            }
+            return h;
+        }
+    };
+}
+
diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index c1ad4f0..ea08214 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -9,6 +9,7 @@
 #include <list>
 #include <exception>
 
+#include "array_hash.hpp"
 
 class graph {
   public:
@@ -25,6 +26,7 @@ class graph {
     typedef struct edge {
       std::array<unsigned int, 2> v;
       coordinate r;
+      std::array<bool, 2> crossings;
     } edge;
 
     coordinate L;
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index d955e43..70173e6 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -12,6 +12,7 @@
 #include "graph.hpp"
 #include "hooks.hpp"
 #include "current_info.hpp"
+#include "array_hash.hpp"
 
 #ifdef FRACTURE_LONGINT
 
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 2fee5fd..7564fa7 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -6,7 +6,6 @@
 #define JCV_REAL_TYPE double
 #define JCV_ATAN2 atan2
 #define JCV_FLT_MAX 1.7976931348623157E+308
-#define NDEBUG
 #include <jc_voronoi.h>
 
 // actual mod for floats
@@ -65,27 +64,6 @@ graph::graph(unsigned int Nx, unsigned int Ny) {
   }
 }
 
-namespace std
-{
-    template<typename T, size_t N>
-    struct hash<array<T, N> >
-    {
-        typedef array<T, N> argument_type;
-        typedef size_t result_type;
-
-        result_type operator()(const argument_type& a) const
-        {
-            hash<T> hasher;
-            result_type h = 0;
-            for (result_type i = 0; i < N; ++i)
-            {
-                h = h * 31 + hasher(a[i]);
-            }
-            return h;
-        }
-    };
-}
-
 class eulerException: public std::exception
 {
   virtual const char* what() const throw()
@@ -248,9 +226,13 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
         }
         
         if (i1 < i2) {
+          unsigned int nn = neighbor->index % 9;
+          bool crossed_x = (nn == 1) || (nn == 2) || (nn > 4);
+          bool crossed_y = nn > 2;
           edges.push_back({{i1, i2},
               {mod((site->p.x + neighbor->p.x) / 2, L.x),
-               mod((site->p.y + neighbor->p.y) / 2, L.y)}
+               mod((site->p.y + neighbor->p.y) / 2, L.y)},
+               {crossed_x, crossed_y}
                });
 
           jcv_graphedge *en;
@@ -276,9 +258,13 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
             vi2 = it2->second;
           }
 
+          bool dcrossed_x = (unsigned int)floor(e->pos[0].x / L.x) != (unsigned int)floor(e->pos[1].x / L.x);
+          bool dcrossed_y = (unsigned int)floor(e->pos[0].y / L.y) != (unsigned int)floor(e->pos[1].y / L.y);
+
           dual_edges.push_back({{vi1, vi2},
               {mod((e->pos[0].x + e->pos[1].x) / 2, L.x),
-               mod((e->pos[0].y + e->pos[1].y) / 2, L.y)}
+               mod((e->pos[0].y + e->pos[1].y) / 2, L.y)},
+               {dcrossed_x, dcrossed_y}
                });
         }
 
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 4d4ed2d..3812f43 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -7,7 +7,7 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
     double v0y = G.vertices[G.edges[i].v[0]].r.y;
     double v1y = G.vertices[G.edges[i].v[1]].r.y;
 
-    if (fabs(v0y - v1y) > G.L.y / 2) {
+    if (G.edges[i].crossings[1]) {
       bool ind = v0y < v1y ? 0 : 1;
 
       ((double *)b->x)[G.edges[i].v[ind]] += 1.0;
@@ -19,8 +19,6 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
 
   cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz,  1, CHOLMOD_REAL, c);
 
-  t->nnz = nnz;
-
   for (unsigned int i = 0; i < G.vertices.size(); i++) {
     ((CHOL_INT *)t->i)[i] = i;
     ((CHOL_INT *)t->j)[i] = i;
@@ -29,26 +27,37 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
 
   unsigned int terms = G.vertices.size();
 
+  std::unordered_map<std::array<unsigned int, 2>, unsigned int> known_edges;
+
   for (unsigned int i = 0; i < G.edges.size(); i++) {
     unsigned int v0 = G.edges[i].v[0];
     unsigned int v1 = G.edges[i].v[1];
 
+    ((double *)t->x)[v0]++;
+    ((double *)t->x)[v1]++;
+
     unsigned int s0 = v0 < v1 ? v0 : v1;
     unsigned int s1 = v0 < v1 ? v1 : v0;
 
-    ((CHOL_INT *)t->i)[terms] = s1;
-    ((CHOL_INT *)t->j)[terms] = s0;
-    ((double *)t->x)[terms] = -1.0;
+    auto it = known_edges.find({s0, s1});
 
-    ((double *)t->x)[v0]++;
-    ((double *)t->x)[v1]++;
+    if (it == known_edges.end()) {
+      ((CHOL_INT *)t->i)[terms] = s1;
+      ((CHOL_INT *)t->j)[terms] = s0;
+      ((double *)t->x)[terms] = -1.0;
 
-    terms++;
+      known_edges[{s0, s1}] = terms;
+      terms++;
+    } else {
+      ((double *)t->x)[it->second] -= 1.0;
+    }
   }
 
   ((double *)t->x)[0]++;
 
-  cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, nnz, c);
+  t->nnz = terms;
+
+  cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, terms, c);
   CHOL_F(free_triplet)(&t, c);
   factor = CHOL_F(analyze)(laplacian, c);
   CHOL_F(factorize)(laplacian, factor, c);
@@ -138,7 +147,7 @@ void network::break_edge(unsigned int e, bool unbreak) {
   double v0y = G.vertices[v0].r.y;
   double v1y = G.vertices[v1].r.y;
 
-  if (fabs(v0y - v1y) > G.L.y / 2) {
+  if (G.edges[e].crossings[1]) {
     bool ind = v0y < v1y ? unbreak : !unbreak;
 
     ((double *)b->x)[G.edges[e].v[ind]] -= 1.0;
@@ -172,7 +181,7 @@ current_info network::get_current_info() {
       double v0y = G.vertices[G.edges[i].v[0]].r.y;
       double v1y = G.vertices[G.edges[i].v[1]].r.y;
 
-      if (fabs(v0y - v1y) > G.L.y / 2) {
+      if (G.edges[i].crossings[1]) {
         if (v0y > v1y) {
           currents[i] += 1.0;
           total_current += currents[i];
diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp
index 1e799bb..34e4ea0 100644
--- a/src/analysis_tools.cpp
+++ b/src/analysis_tools.cpp
@@ -102,12 +102,10 @@ std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
       std::array<unsigned int, 2> crossing_count{0,0};
 
       for (auto edge : cycle) {
-        double dx = fabs(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r.x - n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r.x);
-        if (dx > n.G.L.x / 2) {
+        if (n.G.dual_edges[edge].crossings[0]) {
           crossing_count[0]++;
         }
-        double dy = fabs(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r.y - n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r.y);
-        if (dy > n.G.L.y / 2) {
+        if (n.G.dual_edges[edge].crossings[1]) {
           crossing_count[1]++;
         }
       }
diff --git a/src/animate.cpp b/src/animate.cpp
index 26f4996..a24f0f2 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -25,15 +25,15 @@ void animate::pre_fracture(const network &n) {
     graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
     graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
 
-    if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
-      if (r1.x < n.G.L.x / 2) {
+    if (n.G.edges[i].crossings[0]) {
+      if (r1.x < r2.x) {
         r1.x += n.G.L.x;
       } else {
         r2.x += n.G.L.x;
       }
     }
-    if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
-      if (r1.y < n.G.L.y / 2) {
+    if (n.G.edges[i].crossings[1]) {
+      if (r1.y < r2.y) {
         r1.y += n.G.L.y;
       } else {
         r2.y += n.G.L.y;
@@ -63,16 +63,16 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned in
     graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
     graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
 
-    if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
-      if (r1.x < n.G.L.x / 2) {
-        r2.x -= n.G.L.x;
+    if (n.G.edges[i].crossings[0]) {
+      if (r1.x < r2.x) {
+        r1.x += n.G.L.x;
       } else {
         r2.x += n.G.L.x;
       }
     }
-    if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
-      if (r1.y < n.G.L.y / 2) {
-        r2.y -= n.G.L.y;
+    if (n.G.edges[i].crossings[1]) {
+      if (r1.y < r2.y) {
+        r1.y += n.G.L.y;
       } else {
         r2.y += n.G.L.y;
       }
@@ -111,20 +111,20 @@ void animate::post_fracture(network &n) {
         graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
         graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
 
-        if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
-          if (r1.x < n.G.L.x / 2) {
-            r1.x += n.G.L.x;
-          } else {
-            r2.x += n.G.L.x;
-          }
-        }
-        if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
-          if (r1.y < n.G.L.y / 2) {
-            r1.y += n.G.L.y;
-          } else {
-            r2.y += n.G.L.y;
-          }
-        }
+    if (n.G.edges[i].crossings[0]) {
+      if (r1.x < r2.x) {
+        r1.x += n.G.L.x;
+      } else {
+        r2.x += n.G.L.x;
+      }
+    }
+    if (n.G.edges[i].crossings[1]) {
+      if (r1.y < r2.y) {
+        r1.y += n.G.L.y;
+      } else {
+        r2.y += n.G.L.y;
+      }
+    }
 
         glVertex2d(r1.x, r1.y);
         glVertex2d(r2.x, r2.y);