added animation example, and did many fixes to the voronoi system

14 files changed, 849 insertions, 315 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index a804ac5..35e8436 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 cmake_minimum_required(VERSION 3.8)
 
-set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall -stdlib=libc++")
-set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -stdlib=libc++")
+set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall -stdlib=libc++ -fuse-ld=lld")
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -stdlib=libc++ -march=native -fuse-ld=lld")
 
 set (CMAKE_CXX_STANDARD 17)
 set (CMAKE_C_STANDARD 11)
diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index 75d54f4..d746ddd 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -4,8 +4,11 @@
 #include <cmath>
 #include <vector>
 #include <array>
-#include <map>
+#include <unordered_map>
 #include <random>
+#include <list>
+#include <exception>
+
 
 class graph {
   public:
@@ -16,9 +19,13 @@ class graph {
 
     typedef struct vertex {
       coordinate r;
+      std::vector<coordinate> polygon;
     } vertex;
 
-    typedef std::array<unsigned int, 2> edge;
+    typedef struct edge {
+      std::array<unsigned int, 2> v;
+      coordinate r;
+    } edge;
 
     coordinate L;
 
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index f25cf87..10f62b7 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -1,7 +1,6 @@
 
 #include "graph.hpp"
 #include <cstring>
-#include <iostream>
 
 #define JC_VORONOI_IMPLEMENTATION
 #define JCV_REAL_TYPE double
@@ -24,9 +23,21 @@ graph::graph(unsigned int Nx, unsigned int Ny) {
   for (unsigned int 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));
+    vertices[i].polygon = {
+      {vertices[i].r.x - 0.5, vertices[i].r.y},
+      {vertices[i].r.x, vertices[i].r.y - 0.5},
+      {vertices[i].r.x + 0.5, vertices[i].r.y},
+      {vertices[i].r.x, vertices[i].r.y + 0.5}
+    };
 
     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].polygon = {
+      {dual_vertices[i].r.x - 0.5, vertices[i].r.y},
+      {dual_vertices[i].r.x, vertices[i].r.y - 0.5},
+      {dual_vertices[i].r.x + 0.5, vertices[i].r.y},
+      {dual_vertices[i].r.x, vertices[i].r.y + 0.5}
+    };
   }
 
   for (unsigned int x = 0; x < Ny; x++) {
@@ -34,40 +45,70 @@ graph::graph(unsigned int Nx, unsigned int Ny) {
       unsigned int v1 = (Nx * x) / 2 + ((y + x % 2) / 2) % (Nx / 2);
       unsigned int v2 = ((Nx * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2)) % nv;
 
-      edges.push_back({v1, v2});
+      edges.push_back({{v1, v2}, {0.5 + (double)y, 0.5 + (double)x}});
 
       unsigned int dv1 = (Nx * x) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2);
       unsigned int dv2 = ((Nx * (x + 1)) / 2 + ((y + x % 2) / 2) % (Nx / 2)) % nv;
 
-      dual_edges.push_back({dv1, dv2});
+      dual_edges.push_back({{dv1, dv2}, {0.5 + (double)y, 0.5 + (double)x}});
     }
   }
 }
 
+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()
+  {
+    return "The voronoi tessellation generated has the incorrect Euler characteristic for a torus and is malformed.";
+  }
+} eulerex;
+
+
 graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) {
   L = {(double)Nx, (double)Ny};
 
-  unsigned int dnv = Nx * Ny / 2;
+  unsigned int nv = Nx * Ny / 2;
 
-  dual_vertices.resize(dnv);
+  vertices.resize(nv);
   
   std::normal_distribution<double> d(0.0, spread);
 
-  // the coordinates of the dual lattice, from which a delaunay triangulation
-  // and corresponding voronoi tesselation will be built
-  for (unsigned int i = 0; i < dnv; i++) {
+  // the coordinates of the lattice, from which a delaunay triangulation
+  // and corresponding voronoi tessellation will be built
+  for (unsigned int i = 0; i < nv; i++) {
     double rx = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))) + d(rng);
     double ry = (double)(i / (Nx / 2)) + d(rng);
-    dual_vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}};
+    vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}};
   }
 
-  // to make the resulting tesselation periodic, we tile eight (!) copies of
-  // the original dual points for a total of nine. note that the index of each
+  // to make the resulting tessellation periodic, we tile eight (!) copies of
+  // the original points for a total of nine. note that the index of each
   // point quotient 9 is equal to the original index (we'll use this to
   // translate back)
   std::vector<jcv_point> points;
-  points.reserve(9 * dnv);
-  for (const vertex &v : dual_vertices) {
+  points.reserve(9 * nv);
+  for (const vertex &v : vertices) {
     points.push_back({v.r.x, v.r.y});
     points.push_back({v.r.x + L.x, v.r.y + 0.0});
     points.push_back({v.r.x - L.x, v.r.y + 0.0});
@@ -82,81 +123,151 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread)
   jcv_diagram diagram;
   memset(&diagram, 0, sizeof(jcv_diagram));
 
-  jcv_diagram_generate(9 * dnv, points.data(), NULL, &diagram);
+  jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram);
 
   const jcv_site* sites = jcv_diagram_get_sites(&diagram);
 
-  struct coorcomp {
-    bool operator() (const coordinate& lhs, const coordinate& rhs) const
-    {return lhs.x<rhs.x;}
+  struct paircomp {
+    bool operator() (const std::array<unsigned int, 2>& lhs, const std::array<unsigned int, 2>& rhs) const
+    {if (lhs[0] != lhs[1]) return lhs[0] < lhs[1];
+      else return rhs[0] < rhs[1];
+    }
   };
 
-  std::map<coordinate, unsigned int, coorcomp> known_vertices;
+  std::unordered_map<std::array<unsigned int, 3>, unsigned int> known_vertices;
 
   for (int i = 0; i < diagram.numsites; i++) {
     const jcv_site* site = &sites[i];
 
     // we only care about processing the cells of our original, central sites
     if (site->index % 9 == 0) {
+      unsigned int i1 = (unsigned int)(site->index / 9);
       const jcv_graphedge* e = site->edges;
+      const jcv_graphedge* ep = site->edges;
+      while (ep->next) {
+        ep = ep->next;
+      }
       // for each edge on the site's cell
       while(e) {
-        // assess whether the dual edge needs to be added. only neighboring
-        // dual sites whose index is larger than ours are given edges, so we
+        // assess whether the edge needs to be added. only neighboring
+        // sites whose index is larger than ours are given edges, so we
         // don't double up later
         const jcv_site* neighbor = e->neighbor;
-        unsigned int index = (unsigned int)(site->index / 9);
-        unsigned int real_index = (unsigned int)(neighbor->index / 9);
-        
-        if (index < real_index) {
-          dual_edges.push_back({index, real_index});
-          coordinate r1 = {1e-10 * round(1e10 * fmod(L.x + e->pos[0].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[0].y, L.y))};
-          coordinate r2 = {1e-10 * round(1e10 * fmod(L.x + e->pos[1].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[1].y, L.y))};
+        unsigned int i2 = (unsigned int)(neighbor->index / 9);
+
+        vertices[i1].polygon.push_back({e->pos[0].x, e->pos[0].y});
+
+        unsigned int i3p = (unsigned int)(ep->neighbor->index / 9);
+        std::array<unsigned int, 3> t1 = {i1, i2, i3p};
+        std::sort(t1.begin(), t1.end());
 
-          std::map<coordinate, unsigned int>::iterator it1 = known_vertices.find(r1);
-          std::map<coordinate, unsigned int>::iterator it2 = known_vertices.find(r2);
+        auto it1 = known_vertices.find(t1);
 
-          unsigned int vi1, vi2;
+        unsigned int vi1;
 
-          if (it1 == known_vertices.end()) {
-            vi1 = vertices.size();
-            vertices.push_back({r1});
-            known_vertices[r1] = vi1;
+        if (it1 == known_vertices.end()) {
+          vi1 = dual_vertices.size();
+          dual_vertices.push_back({{fmod(L.x + e->pos[0].x, L.x), fmod(L.y + e->pos[0].y, L.y)},{{site->p.x, site->p.y}}});
+          known_vertices[t1] = vi1;
+        } else {
+          vi1 = it1->second;
+          dual_vertices[vi1].polygon.push_back({site->p.x, site->p.y});
+        }
+        
+        if (i1 < i2) {
+          edges.push_back({{i1, i2},
+              {fmod(L.x + (site->p.x + neighbor->p.x) / 2, L.x),
+               fmod(L.y + (site->p.y + neighbor->p.y) / 2, L.y)}
+               });
+
+          jcv_graphedge *en;
+          if (e->next == NULL) {
+            en = site->edges;
           } else {
-            vi1 = it1->second;
+            en = e->next;
           }
+
+          unsigned int i3n = (unsigned int)(en->neighbor->index / 9);
+          std::array<unsigned int, 3> t2 = {i1, i2, i3n};
+          std::sort(t2.begin(), t2.end());
+
+          auto it2 = known_vertices.find(t2);
+
+          unsigned int vi2;
+
           if (it2 == known_vertices.end()) {
-            vi2 = vertices.size();
-            vertices.push_back({r2});
-            known_vertices[r2] = vi2;
+            vi2 = dual_vertices.size();
+            dual_vertices.push_back({{fmod(L.x + e->pos[1].x, L.x), fmod(L.y + e->pos[1].y, L.y)},{}});
+            known_vertices[t2] = vi2;
           } else {
             vi2 = it2->second;
           }
 
-          edges.push_back({vi1, vi2});
+          dual_edges.push_back({{vi1, vi2},
+              {fmod(L.x + (e->pos[0].x + e->pos[1].x) / 2, L.x),
+               fmod(L.y + (e->pos[0].y + e->pos[1].y) / 2, L.y)}
+               });
         }
 
+        ep = e;
         e = e->next;
       }
     }
   }
 
-  for (edge &e : edges) { 
-    std::cout << e[0] << " " << e[1] << " ";
-  }
-  std::cout << "\n";
-  for (vertex &v : vertices) { 
-    std::cout << v.r.x << " " << v.r.y << " ";
+  if (edges.size() != vertices.size() + dual_vertices.size()) {
+    throw eulerex;
   }
-  std::cout << "\n";
-  for (edge &e : dual_edges) { 
-    std::cout << e[0] << " " << e[1] << " ";
-  }
-  std::cout << "\n";
-  for (vertex &v : dual_vertices) { 
-    std::cout << v.r.x << " " << v.r.y << " ";
+
+  for (vertex &v : dual_vertices) {
+    if (fabs(v.polygon[0].x - v.polygon[1].x) > L.x / 2) {
+      if (v.polygon[0].x < L.x / 2) {
+        v.polygon[0].x += L.x;
+      } else {
+        v.polygon[1].x += L.x;
+      }
+    }
+
+    if (fabs(v.polygon[0].y - v.polygon[1].y) > L.y / 2) {
+      if (v.polygon[0].y < L.y / 2) {
+        v.polygon[0].y += L.y;
+      } else {
+        v.polygon[1].y += L.y;
+      }
+    }
+
+    if (fabs(v.polygon[2].x - v.polygon[1].x) > L.x / 2) {
+      if (v.polygon[2].x < L.x / 2) {
+        v.polygon[2].x += L.x;
+      } else {
+        v.polygon[1].x += L.x;
+      }
+    }
+
+    if (fabs(v.polygon[2].y - v.polygon[1].y) > L.y / 2) {
+      if (v.polygon[2].y < L.y / 2) {
+        v.polygon[2].y += L.y;
+      } else {
+        v.polygon[1].y += L.y;
+      }
+    }
+
+    if (fabs(v.polygon[2].x - v.polygon[0].x) > L.x / 2) {
+      if (v.polygon[2].x < L.x / 2) {
+        v.polygon[2].x += L.x;
+      } else {
+        v.polygon[0].x += L.x;
+      }
+    }
+
+    if (fabs(v.polygon[2].y - v.polygon[0].y) > L.y / 2) {
+      if (v.polygon[2].y < L.y / 2) {
+        v.polygon[2].y += L.y;
+      } else {
+        v.polygon[0].y += L.y;
+      }
+    }
   }
-  std::cout << "\n";
 
   jcv_diagram_free(&diagram);
 }
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 9cb1007..4d4ed2d 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -4,14 +4,14 @@
 network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) {
   b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c);
   for (unsigned int i = 0; i < G.edges.size(); i++) {
-    double v0y = G.vertices[G.edges[i][0]].r.y;
-    double v1y = G.vertices[G.edges[i][1]].r.y;
+    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) {
       bool ind = v0y < v1y ? 0 : 1;
 
-      ((double *)b->x)[G.edges[i][ind]] += 1.0;
-      ((double *)b->x)[G.edges[i][!ind]] -= 1.0;
+      ((double *)b->x)[G.edges[i].v[ind]] += 1.0;
+      ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0;
     }
   }
 
@@ -30,8 +30,8 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
   unsigned int terms = G.vertices.size();
 
   for (unsigned int i = 0; i < G.edges.size(); i++) {
-    unsigned int v0 = G.edges[i][0];
-    unsigned int v1 = G.edges[i][1];
+    unsigned int v0 = G.edges[i].v[0];
+    unsigned int v1 = G.edges[i].v[1];
 
     unsigned int s0 = v0 < v1 ? v0 : v1;
     unsigned int s1 = v0 < v1 ? v1 : v0;
@@ -60,11 +60,11 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
 
   for (unsigned int i = 0; i < G.edges.size(); i++) {
     ((CHOL_INT *)t->i)[2 * i] = i;
-    ((CHOL_INT *)t->j)[2 * i] = G.edges[i][0];
+    ((CHOL_INT *)t->j)[2 * i] = G.edges[i].v[0];
     ((double *)t->x)[2 * i] = 1.0;
 
     ((CHOL_INT *)t->i)[2 * i + 1] = i;
-    ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i][1];
+    ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1];
     ((double *)t->x)[2 * i + 1] = -1.0;
   }
 
@@ -99,8 +99,8 @@ void network::set_thresholds(double beta, std::mt19937& rng) {
 
 void network::break_edge(unsigned int e, bool unbreak) {
   fuses[e] = !unbreak;
-  unsigned int v0 = G.edges[e][0];
-  unsigned int v1 = G.edges[e][1];
+  unsigned int v0 = G.edges[e].v[0];
+  unsigned int v1 = G.edges[e].v[1];
 
   unsigned int n = factor->n;
 
@@ -141,8 +141,8 @@ void network::break_edge(unsigned int e, bool unbreak) {
   if (fabs(v0y - v1y) > G.L.y / 2) {
     bool ind = v0y < v1y ? unbreak : !unbreak;
 
-    ((double *)b->x)[G.edges[e][ind]] -= 1.0;
-    ((double *)b->x)[G.edges[e][!ind]] += 1.0;
+    ((double *)b->x)[G.edges[e].v[ind]] -= 1.0;
+    ((double *)b->x)[G.edges[e].v[!ind]] += 1.0;
   }
 }
 
@@ -169,8 +169,8 @@ current_info network::get_current_info() {
     } else {
       currents[i] = ((double *)y->x)[i];
 
-      double v0y = G.vertices[G.edges[i][0]].r.y;
-      double v1y = G.vertices[G.edges[i][1]].r.y;
+      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 (v0y > v1y) {
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index d93783c..283eb10 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,7 +1,18 @@
 
+add_library(analysis_tools analysis_tools.cpp)
+target_link_libraries(analysis_tools frac)
+
 add_library(measurements measurements.cpp)
-target_link_libraries(measurements frac)
+target_link_libraries(measurements frac analysis_tools)
+
+add_library(animate animate.cpp)
+target_link_libraries(animate frac analysis_tools GL GLU glut)
 
 add_executable(fracture fracture.cpp)
-target_link_libraries(fracture frac measurements fftw3 cholmod profiler)
+target_link_libraries(fracture frac measurements fftw3 cholmod profiler GL GLU glut)
+
+add_executable(animate_fracture animate_fracture.cpp)
+target_link_libraries(animate_fracture frac animate fftw3 cholmod profiler GL GLU glut)
 
+add_executable(test test.cpp)
+target_link_libraries(test frac measurements fftw3 cholmod profiler)
diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp
new file mode 100644
index 0000000..1e799bb
--- /dev/null
+++ b/src/analysis_tools.cpp
@@ -0,0 +1,125 @@
+
+#include "analysis_tools.hpp"
+
+template <class T>
+bool is_shorter(const std::list<T> &l1, const std::list<T> &l2) {
+  return l1.size() < l2.size();
+}
+
+bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) {
+  return true;
+}
+
+std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
+  Graph Gtmp(n.G.vertices.size());
+  std::list<unsigned int> removed_edges;
+
+  class add_tree_edges : public boost::default_dfs_visitor {
+    public:
+      Graph& G;
+      std::list<unsigned int>& E;
+
+      add_tree_edges(Graph& G, std::list<unsigned int>& E) : G(G), E(E) {}
+
+      void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
+        boost::add_edge(boost::source(e, g), boost::target(e, g), g[e], G);
+      }
+
+      void back_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
+        if (!(boost::edge(boost::source(e, g), boost::target(e, g), G).second)) {
+          E.push_back(g[e].index);
+        }
+      }
+  };
+
+  add_tree_edges ate(Gtmp, removed_edges);
+  boost::depth_first_search(G, visitor(ate));
+
+  class find_cycle : public boost::default_dfs_visitor {
+    public:
+      std::list<unsigned int>& E;
+      unsigned int end;
+      struct done{};
+
+      find_cycle(std::list<unsigned int>& E, unsigned int end) : E(E), end(end) {}
+
+      void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
+        if (v == end) {
+          throw done{};
+        }
+      }
+
+      void examine_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
+        E.push_back(g[e].index);
+      }
+
+      void finish_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
+        E.erase(std::find(E.begin(), E.end(), g[e].index));
+      }
+  };
+
+  std::list<std::list<unsigned int>> cycles;
+
+  for (auto edge : removed_edges) {
+    std::list<unsigned int> cycle = {edge};
+    find_cycle vis(cycle, n.G.dual_edges[edge].v[1]);
+    std::vector<boost::default_color_type> new_color_map(boost::num_vertices(Gtmp));
+    try {
+    boost::depth_first_visit(Gtmp, n.G.dual_edges[edge].v[0], vis, boost::make_iterator_property_map(new_color_map.begin(), boost::get(boost::vertex_index, Gtmp), new_color_map[0]));
+    } catch(find_cycle::done const&) {
+      cycles.push_back(cycle);
+    }
+  }
+
+  if (cycles.size() > 1) {
+    std::list<std::valarray<uint8_t>> bool_cycles;
+    for (auto cycle : cycles) {
+      std::valarray<uint8_t> bool_cycle(n.G.edges.size());
+      for (auto v : cycle) {
+        bool_cycle[v] = 1;
+      }
+      bool_cycles.push_back(bool_cycle);
+    }
+
+    // generate all possible cycles by taking xor of the edge sets of the known cycles
+    for (auto it1 = bool_cycles.begin(); it1 != std::prev(bool_cycles.end()); it1++) {
+      for (auto it2 = std::next(it1); it2 != bool_cycles.end(); it2++) {
+        std::valarray<uint8_t> new_bool_cycle = (*it1) ^ (*it2);
+        std::list<unsigned int> new_cycle;
+        unsigned int pos = 0;
+        for (uint8_t included : new_bool_cycle) {
+          if (included) {
+            new_cycle.push_back(pos);
+          }
+          pos++;
+        }
+        cycles.push_back(new_cycle);
+      }
+    }
+
+    // find the cycle representing the crack by counting boundary crossings
+    for (auto cycle : cycles) {
+      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) {
+          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) {
+          crossing_count[1]++;
+        }
+      }
+
+      if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
+        return cycle;
+      }
+    }
+  } else {
+    return cycles.front();
+  }
+
+  exit(5);
+}
+
diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp
new file mode 100644
index 0000000..34ed967
--- /dev/null
+++ b/src/analysis_tools.hpp
@@ -0,0 +1,29 @@
+
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/connected_components.hpp>
+#include <boost/graph/depth_first_search.hpp>
+#include <boost/range/combine.hpp>
+#include <boost/foreach.hpp>
+
+#include <vector>
+#include <algorithm>
+#include <cmath>
+#include <list>
+#include <valarray>
+
+#include <network.hpp>
+
+struct EdgeProperties {
+  unsigned int index;
+};
+
+typedef boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, boost::no_property, EdgeProperties> Graph;
+
+
+template<class T>
+bool is_shorter(const std::list<T> &, const std::list<T> &);
+
+bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>);
+
+std::list<unsigned int> find_minimal_crack(const Graph &, const network &);
+
diff --git a/src/animate.cpp b/src/animate.cpp
new file mode 100644
index 0000000..8c2d12f
--- /dev/null
+++ b/src/animate.cpp
@@ -0,0 +1,205 @@
+
+#include "animate.hpp"
+
+animate::animate(unsigned int Lx, unsigned int Ly, unsigned int window_size, int argc, char *argv[]) : Lx(Lx), Ly(Ly), G(Lx * Ly) {
+  glutInit(&argc, argv);
+  glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
+  glutInitWindowSize(window_size, (unsigned int )(window_size * ((double)Ly / (double)Lx)));
+  glutCreateWindow("wolff");
+  glClearColor(0.0,0.0,0.0,0.0);
+  glMatrixMode(GL_PROJECTION);
+  glLoadIdentity();
+  gluOrtho2D(-1.0, Lx + 3.0, -1.0 , Ly + 3.0);
+}
+
+void animate::pre_fracture(const network &n) {
+  lv = std::numeric_limits<long double>::lowest();
+  avalanches = {{}};
+  boost::remove_edge_if(trivial, G);
+
+  glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
+  glClear(GL_COLOR_BUFFER_BIT);
+  glBegin(GL_LINES);
+  glColor3f(0.0f, 0.0f, 0.0f);
+  for (unsigned int i = 0; i < n.G.edges.size(); i++) {
+    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) > Lx / 2) {
+      if (r1.x < Lx / 2) {
+        r1.x += Lx;
+      } else {
+        r2.x += Lx;
+      }
+    }
+    if (fabs(r1.y - r2.y) > Ly / 2) {
+      if (r1.y < Ly / 2) {
+        r1.y += Ly;
+      } else {
+        r2.y += Ly;
+      }
+    }
+
+    glVertex2d(r1.x, r1.y);
+    glVertex2d(r2.x, r2.y);
+  }
+  glEnd();
+  glFlush();
+}
+
+void animate::bond_broken(const network& net, const current_info& cur, unsigned int i) {
+  long double c = logl(cur.conductivity / fabs(cur.currents[i])) + net.thresholds[i];
+  if (c > lv && avalanches.back().size() > 0) {
+    lv = c;
+    avalanches.push_back({i});
+  } else {
+    avalanches.back().push_back(i);
+  }
+
+  boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
+
+   glBegin(GL_LINES);          // Each set of 3 vertices form a triangle
+    glColor3f(1.0f, 1.0f, 1.0f); // Blue
+    graph::coordinate r1 = net.G.vertices[net.G.edges[i].v[0]].r;
+    graph::coordinate r2 = net.G.vertices[net.G.edges[i].v[1]].r;
+
+    if (fabs(r1.x - r2.x) > Lx / 2) {
+      if (r1.x < Lx / 2) {
+        r2.x -= Lx;
+      } else {
+        r2.x += Lx;
+      }
+    }
+    if (fabs(r1.y - r2.y) > Ly / 2) {
+      if (r1.y < Ly / 2) {
+        r2.y -= Ly;
+      } else {
+        r2.y += Ly;
+      }
+    }
+
+    glVertex2d(r1.x, r1.y);
+    glVertex2d(r2.x, r2.y);
+ glEnd();
+ glFlush();
+}
+
+void animate::post_fracture(network &n) {
+  std::vector<unsigned int> component(boost::num_vertices(G));
+  unsigned int num = boost::connected_components(G, &component[0]);
+
+  std::list<unsigned int> crack = find_minimal_crack(G, n);
+  unsigned int crack_component = component[n.G.dual_edges[crack.front()].v[0]];
+
+  std::vector<std::list<unsigned int>> components(num);
+
+  for (unsigned int i = 0; i < n.G.dual_vertices.size(); i++) {
+    components[component[i]].push_back(i);
+  }
+
+  std::default_random_engine gen;
+  std::uniform_real_distribution<double> dis(0.0,1.0);
+
+  char key;
+  while ((key = getchar()) != 'n') {
+    glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
+    glClear(GL_COLOR_BUFFER_BIT);
+    glBegin(GL_LINES);
+    glColor3f(0.0f, 0.0f, 0.0f);
+    for (unsigned int i = 0; i < n.G.edges.size(); i++) {
+      if (!n.fuses[i]) {
+        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) > Lx / 2) {
+          if (r1.x < Lx / 2) {
+            r1.x += Lx;
+          } else {
+            r2.x += Lx;
+          }
+        }
+        if (fabs(r1.y - r2.y) > Ly / 2) {
+          if (r1.y < Ly / 2) {
+            r1.y += Ly;
+          } else {
+            r2.y += Ly;
+          }
+        }
+
+        glVertex2d(r1.x, r1.y);
+        glVertex2d(r2.x, r2.y);
+      }
+    }
+    glEnd();
+
+    switch (key) {
+      case 's' :
+        for (auto edge : crack) {
+         glBegin(GL_TRIANGLES);
+          glColor3f(1.0f, 0.0f, 0.0f);
+          glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[0].x, n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[0].y);
+          glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[1].x, n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[1].y);
+          glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[2].x, n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[2].y);
+
+          glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[0].x, n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[0].y);
+          glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[1].x, n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[1].y);
+          glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[2].x, n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[2].y);
+          glEnd();
+        }
+        glFlush();
+        break;
+
+      case 'c' :
+        glBegin(GL_TRIANGLES);          // Each set of 3 vertices form a triangle
+        for (unsigned int i = 0; i < num; i++) {
+          if (i == crack_component) {
+            glColor3d(1.0, 0.0, 0.0);
+          } else {
+            glColor3d(dis(gen), dis(gen), dis(gen));
+          }
+
+          for (auto it = components[i].begin(); it != components[i].end(); it++) {
+            glVertex2d(n.G.dual_vertices[*it].polygon[0].x, n.G.dual_vertices[*it].polygon[0].y);
+            glVertex2d(n.G.dual_vertices[*it].polygon[1].x, n.G.dual_vertices[*it].polygon[1].y);
+            glVertex2d(n.G.dual_vertices[*it].polygon[2].x, n.G.dual_vertices[*it].polygon[2].y);
+          }
+        }
+        glEnd();
+        glFlush();
+        break;
+
+      case 'C' :
+        glBegin(GL_TRIANGLES);          // Each set of 3 vertices form a triangle
+        for (unsigned int i = 0; i < num; i++) {
+          if (components[i].size() > 1) {
+            if (i == crack_component) {
+              glColor3d(1.0, 0.0, 0.0);
+            } else {
+              glColor3d(dis(gen), dis(gen), dis(gen));
+            }
+
+            for (auto it = components[i].begin(); it != components[i].end(); it++) {
+              glVertex2d(n.G.dual_vertices[*it].polygon[0].x, n.G.dual_vertices[*it].polygon[0].y);
+              glVertex2d(n.G.dual_vertices[*it].polygon[1].x, n.G.dual_vertices[*it].polygon[1].y);
+              glVertex2d(n.G.dual_vertices[*it].polygon[2].x, n.G.dual_vertices[*it].polygon[2].y);
+            }
+          }
+        }
+        glEnd();
+        glFlush();
+        break;
+
+      case 'm' :
+        glBegin(GL_TRIANGLES);
+        glColor3d(1.0, 0.0, 0.0);
+        for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
+          glVertex2d(n.G.dual_vertices[*it].polygon[0].x, n.G.dual_vertices[*it].polygon[0].y);
+          glVertex2d(n.G.dual_vertices[*it].polygon[1].x, n.G.dual_vertices[*it].polygon[1].y);
+          glVertex2d(n.G.dual_vertices[*it].polygon[2].x, n.G.dual_vertices[*it].polygon[2].y);
+        }
+        glEnd();
+        glFlush();
+    }
+  }
+}
+
diff --git a/src/animate.hpp b/src/animate.hpp
new file mode 100644
index 0000000..793771d
--- /dev/null
+++ b/src/animate.hpp
@@ -0,0 +1,24 @@
+
+
+#include <functional>
+#include <hooks.hpp>
+
+#include "analysis_tools.hpp"
+
+#include <GL/glut.h>
+
+class animate : public hooks {
+  private:
+    unsigned int Lx;
+    unsigned int Ly;
+    Graph G;
+  public:
+    long double lv;
+    std::list<std::list<unsigned int>> avalanches;
+
+    animate(unsigned int Lx, unsigned int Ly, unsigned int window_size, int argc, char *argv[]);
+
+    void pre_fracture(const network &) override;
+    void bond_broken(const network& net, const current_info& cur, unsigned int i) override;
+    void post_fracture(network &n) override;
+};
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
new file mode 100644
index 0000000..10b63d0
--- /dev/null
+++ b/src/animate_fracture.cpp
@@ -0,0 +1,79 @@
+
+#include <random>
+#include <iostream>
+
+#include <cholmod.h>
+
+#include "randutils/randutils.hpp"
+
+#include <graph.hpp>
+#include <network.hpp>
+#include <hooks.hpp>
+#include "animate.hpp"
+
+#include <csignal>
+#include <cstring>
+#include <atomic>
+
+std::atomic<bool> quit(false);    // signal flag
+
+void got_signal(int) {
+    quit.store(true);
+}
+
+int main(int argc, char* argv[]) {
+  struct sigaction sa;
+  memset( &sa, 0, sizeof(sa) );
+  sa.sa_handler = got_signal;
+  sigfillset(&sa.sa_mask);
+  sigaction(SIGINT, &sa, NULL);
+
+  int opt;
+
+  unsigned int N = 1;
+  unsigned int Lx = 16;
+  unsigned int Ly = 16;
+  double  beta = 0.5;
+
+  while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
+    switch (opt) {
+      case 'N':
+        N = (unsigned int)atof(optarg);
+        break;
+      case 'X':
+        Lx = atoi(optarg);
+        break;
+      case 'Y':
+        Ly = atoi(optarg);
+        break;
+      case 'b':
+        beta = atof(optarg);
+        break;
+      default:
+        exit(1);
+    }
+  }
+
+  cholmod_common c;
+  CHOL_F(start)(&c);
+
+  animate meas(Lx, Ly, 512, argc, argv);
+
+  randutils::auto_seed_128 seeds;
+  std::mt19937 rng{seeds};
+
+  for (unsigned int trial = 0; trial < N; trial++) {
+    graph G(Lx, Ly, rng, 0.4);
+    network network(G, &c);
+    network.set_thresholds(beta, rng);
+    network.fracture(meas);
+
+    if (quit.load())
+      break;
+  }
+
+  CHOL_F(finish)(&c);
+
+  return 0;
+}
+
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 317d91b..4b6fb5b 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -57,17 +57,22 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  graph G(Lx, Ly);
-  network base_network(G, &c);
-  ma meas(Lx, Ly, beta);
+  ma meas(Lx, Ly, Lx, Ly, beta);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
 
   for (unsigned int trial = 0; trial < N; trial++) {
-    network tmp_network(base_network);
-    tmp_network.set_thresholds(beta, rng);
-    tmp_network.fracture(meas);
+    while (true) {
+      try {
+        graph G(Lx, Ly, rng, 0.4);
+        network network(G, &c);
+        network.set_thresholds(beta, rng);
+        network.fracture(meas);
+        break;
+      } catch (std::exception &e) {
+      }
+    }
 
     if (quit.load())
       break;
diff --git a/src/measurements.cpp b/src/measurements.cpp
index c8cc73c..43483be 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -2,128 +2,6 @@
 #include "measurements.hpp"
 #include <iostream>
 
-template <class T>
-bool is_shorter(std::list<T> l1, std::list<T> l2) {
-  return l1.size() < l2.size();
-}
-
-bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) {
-  return true;
-}
-
-std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
-  Graph Gtmp(n.G.vertices.size());
-  std::list<unsigned int> removed_edges;
-
-  class add_tree_edges : public boost::default_dfs_visitor {
-    public:
-      Graph& G;
-      std::list<unsigned int>& E;
-
-      add_tree_edges(Graph& G, std::list<unsigned int>& E) : G(G), E(E) {}
-
-      void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        boost::add_edge(boost::source(e, g), boost::target(e, g), g[e], G);
-      }
-
-      void back_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        if (!(boost::edge(boost::source(e, g), boost::target(e, g), G).second)) {
-          E.push_back(g[e].index);
-        }
-      }
-  };
-
-  add_tree_edges ate(Gtmp, removed_edges);
-  boost::depth_first_search(G, visitor(ate));
-
-  class find_cycle : public boost::default_dfs_visitor {
-    public:
-      std::list<unsigned int>& E;
-      unsigned int end;
-      struct done{};
-
-      find_cycle(std::list<unsigned int>& E, unsigned int end) : E(E), end(end) {}
-
-      void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
-        if (v == end) {
-          throw done{};
-        }
-      }
-
-      void examine_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        E.push_back(g[e].index);
-      }
-
-      void finish_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        E.erase(std::find(E.begin(), E.end(), g[e].index));
-      }
-  };
-
-  std::list<std::list<unsigned int>> cycles;
-
-  for (auto edge : removed_edges) {
-    std::list<unsigned int> cycle = {edge};
-    find_cycle vis(cycle, n.G.dual_edges[edge][1]);
-    std::vector<boost::default_color_type> new_color_map(boost::num_vertices(Gtmp));
-    try {
-    boost::depth_first_visit(Gtmp, n.G.dual_edges[edge][0], vis, boost::make_iterator_property_map(new_color_map.begin(), boost::get(boost::vertex_index, Gtmp), new_color_map[0]));
-    } catch(find_cycle::done const&) {
-      cycles.push_back(cycle);
-    }
-  }
-
-  if (cycles.size() > 1) {
-    std::list<std::valarray<uint8_t>> bool_cycles;
-    for (auto cycle : cycles) {
-      std::valarray<uint8_t> bool_cycle(n.G.edges.size());
-      for (auto v : cycle) {
-        bool_cycle[v] = 1;
-      }
-      bool_cycles.push_back(bool_cycle);
-    }
-
-    // generate all possible cycles by taking xor of the edge sets of the known cycles
-    for (auto it1 = bool_cycles.begin(); it1 != std::prev(bool_cycles.end()); it1++) {
-      for (auto it2 = std::next(it1); it2 != bool_cycles.end(); it2++) {
-        std::valarray<uint8_t> new_bool_cycle = (*it1) ^ (*it2);
-        std::list<unsigned int> new_cycle;
-        unsigned int pos = 0;
-        for (uint8_t included : new_bool_cycle) {
-          if (included) {
-            new_cycle.push_back(pos);
-          }
-          pos++;
-        }
-        cycles.push_back(new_cycle);
-      }
-    }
-
-    // find the cycle representing the crack by counting boundary crossings
-    for (auto cycle : cycles) {
-      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][0]].r.x - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.x);
-        if (dx > n.G.L.x / 2) {
-          crossing_count[0]++;
-        }
-        double dy = fabs(n.G.dual_vertices[n.G.dual_edges[edge][0]].r.y - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.y);
-        if (dy > n.G.L.y / 2) {
-          crossing_count[1]++;
-        }
-      }
-
-      if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
-        return cycle;
-      }
-    }
-  } else {
-    return cycles.front();
-  }
-
-  exit(5);
-}
-
 void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta) {
   std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat";
   std::ifstream file(filename);
@@ -153,8 +31,8 @@ void update_distribution_file(std::string id, const std::vector<uint64_t>& data,
 }
 
 template <class T>
-void update_field_file(std::string id, const std::vector<T>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta) {
-  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat";
+void update_field_file(std::string id, const std::vector<T>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta, unsigned int Mx, unsigned int My) {
+  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat";
   std::ifstream file(filename);
 
   uint64_t N_old = 0;
@@ -226,18 +104,22 @@ void autocorrelation(unsigned int Lx, unsigned int Ly, std::vector<T>& out_data,
   }
 }
 
-ma::ma(unsigned int Lx, unsigned int Ly, double beta) : Lx(Lx), Ly(Ly), beta(beta), G(Lx * Ly / 2),
+unsigned int edge_r_to_ind(graph::coordinate r, unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My) {
+  return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);
+}
+
+ma::ma(unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My, double beta) : Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(Lx * Ly),
   sc(Lx * Ly, 0),
   sa(Lx * Ly, 0),
   sd(Lx * Ly, 0),
-  sb(log2(Ly) + 1, 0),
-  Ccc((Lx / 2 + 1) * (Ly / 2 + 1), 0),
-  Css((Lx / 2 + 1) * (Ly / 2 + 1), 0),
-  Cmm((Lx / 2 + 1) * (Ly / 2 + 1), 0),
-  Caa((Lx / 2 + 1) * (Ly / 2 + 1), 0),
-  Cdd((Lx / 2 + 1) * (Ly / 2 + 1), 0),
-  Cll((Lx / 2 + 1) * (Ly / 2 + 1), 0),
-  Cee((Lx / 2 + 1) * (Ly / 2 + 1), 0)
+  sb(log2(Mx < My ? Mx : My) + 1, 0),
+  Ccc((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Css((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Cmm((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Caa((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Cdd((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Cll((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Cee((Mx / 2 + 1) * (My / 2 + 1), 0)
 {
   N = 0;
   Nc = 0;
@@ -246,17 +128,13 @@ ma::ma(unsigned int Lx, unsigned int Ly, double beta) : Lx(Lx), Ly(Ly), beta(bet
   // FFTW setup for correlation functions
   fftw_set_timelimit(1);
 
-  fftw_forward_in = (double *)fftw_malloc(Lx * Ly * sizeof(double));
-  fftw_forward_out = (fftw_complex *)fftw_malloc(Lx * Ly * sizeof(fftw_complex));
-  fftw_reverse_in = (fftw_complex *)fftw_malloc(Lx * Ly * sizeof(fftw_complex));
-  fftw_reverse_out = (double *)fftw_malloc(Lx * Ly * sizeof(double));
-
-  forward_plan = fftw_plan_dft_r2c_2d(Ly, Lx, fftw_forward_in, fftw_forward_out, 0);
-  reverse_plan = fftw_plan_dft_c2r_2d(Ly, Lx, fftw_reverse_in, fftw_reverse_out, 0);
-
-  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_broken.dat";
+  fftw_forward_in = (double *)fftw_malloc(Mx * My * sizeof(double));
+  fftw_forward_out = (fftw_complex *)fftw_malloc(Mx * My * sizeof(fftw_complex));
+  fftw_reverse_in = (fftw_complex *)fftw_malloc(Mx * My * sizeof(fftw_complex));
+  fftw_reverse_out = (double *)fftw_malloc(Mx * My * sizeof(double));
 
-  bondfile.open(filename, std::ifstream::app);
+  forward_plan = fftw_plan_dft_r2c_2d(My, Mx, fftw_forward_in, fftw_forward_out, 0);
+  reverse_plan = fftw_plan_dft_c2r_2d(My, Mx, fftw_reverse_in, fftw_reverse_out, 0);
 }
 
 ma::~ma() {
@@ -274,18 +152,16 @@ ma::~ma() {
   update_distribution_file("sd", sd, N, Lx, Ly, beta);
   update_distribution_file("sb", sb, N, Lx, Ly, beta);
 
-  update_field_file("Ccc", Ccc, Nc, Lx, Ly, beta);
-  update_field_file("Css", Css, N, Lx, Ly, beta);
-  update_field_file("Cmm", Cmm, N, Lx, Ly, beta);
-  update_field_file("Cdd", Cdd, N, Lx, Ly, beta);
-  update_field_file("Caa", Caa, Na, Lx, Ly, beta);
-  update_field_file("Cll", Cll, N, Lx, Ly, beta);
-  update_field_file("Cee", Cee, N, Lx, Ly, beta);
-
-  bondfile.close();
+  update_field_file("Ccc", Ccc, Nc, Lx, Ly, beta, Mx, My);
+  update_field_file("Css", Css, N, Lx, Ly, beta, Mx, My);
+  update_field_file("Cmm", Cmm, N, Lx, Ly, beta, Mx, My);
+  update_field_file("Cdd", Cdd, N, Lx, Ly, beta, Mx, My);
+  update_field_file("Caa", Caa, Na, Lx, Ly, beta, Mx, My);
+  update_field_file("Cll", Cll, N, Lx, Ly, beta, Mx, My);
+  update_field_file("Cee", Cee, N, Lx, Ly, beta, Mx, My);
 }
 
-void ma::pre_fracture(const network &) {
+void ma::pre_fracture(const network&) {
   lv = std::numeric_limits<long double>::lowest();
   avalanches = {{}};
   boost::remove_edge_if(trivial, G);
@@ -297,13 +173,13 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned int i
     sa[avalanches.back().size() - 1]++;
     Na++;
 
-    std::fill_n(fftw_forward_in, net.G.edges.size(), 0.0);
+    std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
     for (auto e : avalanches.back()) {
-      fftw_forward_in[e] = 1.0;
+      fftw_forward_in[edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0;
     }
 
-    autocorrelation(Lx, Ly, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+    autocorrelation(Mx, My, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
     lv = c;
     avalanches.push_back({i});
@@ -311,121 +187,116 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned int i
     avalanches.back().push_back(i);
   }
 
-  boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], {i}, G);
-
-  bondfile << i << " " << c << std::scientific << " ";
+  boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
 }
 
 void ma::post_fracture(network &n) {
-  bondfile << "\n";
   std::vector<unsigned int> component(boost::num_vertices(G));
   unsigned int num = boost::connected_components(G, &component[0]);
 
   std::list<unsigned int> crack = find_minimal_crack(G, n);
 
-  unsigned int crack_component = component[n.G.dual_edges[crack.front()][0]];
+  // crack surface correlations
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
+
+  for (auto edge : crack) {
+    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My)] = 0.5;
+    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My)] = 0.5;
+  }
+
+  autocorrelation(Mx, My, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+
+  unsigned int crack_component = component[n.G.dual_edges[crack.front()].v[0]];
+
+  std::vector<std::list<unsigned int>> components(num);
+
+  for (unsigned int i = 0; i < n.G.dual_vertices.size(); i++) {
+    components[component[i]].push_back(i);
+  }
 
   // non-spanning clusters
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
   for (unsigned int i = 0; i < num; i++) {
     if (i != crack_component) {
-      unsigned int size = 0;
-
-      for (unsigned int j = 0; j < n.G.edges.size(); j++) {
-        if (component[n.G.dual_edges[j][0]] == i && n.fuses[j]) {
-          size++;
-          fftw_forward_in[j] = 1.0;
-        } else{
-          fftw_forward_in[j] = 0.0;
-        }
+      for (auto it = components[i].begin(); it != components[i].end(); it++) {
+        fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] += 1.0;
       }
 
-      if (size > 0) {
-        sc[size - 1]++;
-        autocorrelation(Lx, Ly, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-        Nc++;
-      }
-    }
-  }
+      sc[components[i].size() - 1]++;
+      autocorrelation(Mx, My, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+      Nc++;
 
-  // bin counting
-  for (unsigned int be = 0; be < log2(Ly); be++) {
-    unsigned int bin = pow(2, be);
-
-    for (unsigned int i = 0; i < Lx * Ly / pow(bin, 2); i++) {
-      bool in_bin = false;
-      for (unsigned int j = 0; j < pow(bin, 2); j++) {
-        unsigned int edge = Lx * (bin * ((i * bin) / Lx) + j / bin) + (i * bin) % Lx + j % bin;
-        if (!n.fuses[edge] && crack_component == component[n.G.dual_edges[edge][0]]) {
-          in_bin = true;
-          break;
-        }
-      }
-
-      if (in_bin) {
-        sb[be]++;
+      for (auto it = components[i].begin(); it != components[i].end(); it++) {
+        fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] = 0.0;
       }
     }
   }
 
-  sb[log2(Ly)]++;
+  unsigned int max_factor = log2(Mx < My ? Mx : My) + 1;
+  std::vector<std::valarray<unsigned int>> bins(max_factor);
+  for (unsigned int i = 0; i < max_factor; i++) {
+    bins[i].resize(Mx * My / (unsigned int)pow(2, 2 * i));
+  }
 
-  for (unsigned int i = 0; i < n.G.edges.size(); i++) {
-    if (n.fuses[i] && component[n.G.dual_edges[i][0]] == crack_component)  {
-      fftw_forward_in[i] = 1.0;
-    } else {
-      fftw_forward_in[i] = 0.0;
+  for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
+    for (unsigned int i = 0; i < max_factor; i++) {
+      bins[i][edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx / pow(2, i), My / pow(2, i))] = 1;
     }
   }
 
-  autocorrelation(Lx, Ly, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  for (unsigned int i =0 ; i < max_factor; i++) {
+    sb[i] += bins[i].sum();
+  }
 
-  // crack surface correlations
-  std::fill_n(fftw_forward_in, n.G.edges.size(), 0.0);
+  // spanning cluster
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
-  for (auto edge : crack) {
-    fftw_forward_in[edge] = 1.0;
+  for (unsigned int i = 0; i < n.G.dual_vertices.size(); i++) {
+    if (component[i] == crack_component)  {
+      fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[i].r, Lx, Ly, Mx, My)] += 1.0;
+    }
   }
 
-  autocorrelation(Lx, Ly, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Mx, My, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   std::function<bool(unsigned int)> inCrack = [&](unsigned int i) -> bool {
-    return (bool)fftw_forward_in[i];
+    return component[n.G.dual_edges[i].v[0]] == crack_component;
   };
 
   for (auto avalanche : avalanches) {
     if (avalanche.end() != std::find_if(avalanche.begin(), avalanche.end(), inCrack)) {
       for (auto edge : avalanche) {
-        fftw_forward_in[edge] = 1.0;
+        fftw_forward_in[edge_r_to_ind(n.G.edges[edge].r, Lx, Ly, Mx, My)] += 1.0;
       }
     }
   }
 
-  autocorrelation(Lx, Ly, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Mx, My, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
-  std::fill_n(fftw_forward_in, n.G.edges.size(), 0.0);
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
   // rewind the last avalanche
   for (auto e : avalanches.back()) {
-    boost::remove_edge(n.G.dual_edges[e][0], n.G.dual_edges[e][1], G);
+    boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
     n.break_edge(e, true);
-    fftw_forward_in[e] = 1;
+    fftw_forward_in[edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0;
   }
 
-  autocorrelation(Lx, Ly, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Mx, My, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   // damage size distribution
   unsigned int total_broken = 0;
 
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
+
   for (unsigned int i = 0; i < n.G.edges.size(); i++) {
     if (n.fuses[i]) { 
       total_broken++;
-      fftw_forward_in[i] = 1.0;
-    } else {
-      fftw_forward_in[i] = 0.0;
+      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
     }
   }
 
-  autocorrelation(Lx, Ly, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Mx, My, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   sd[total_broken]++;
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 2bc5bae..685806b 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -1,34 +1,18 @@
 
-#include <vector>
-#include <algorithm>
-#include <cmath>
 #include <cstring>
-#include <list>
 #include <fstream>
 #include <string>
 #include <cinttypes>
 #include <sstream>
 #include <functional>
 #include <iostream>
-#include <valarray>
 #include <array>
 
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/connected_components.hpp>
-#include <boost/graph/depth_first_search.hpp>
-#include <boost/range/combine.hpp>
-#include <boost/foreach.hpp>
-
 #include <fftw3.h>
 
-#include <network.hpp>
 #include <hooks.hpp>
 
-struct EdgeProperties {
-  unsigned int index;
-};
-
-typedef boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, boost::no_property, EdgeProperties> Graph;
+#include "analysis_tools.hpp"
 
 class ma : public hooks {
   // need:
@@ -44,6 +28,8 @@ class ma : public hooks {
     unsigned int N;
     unsigned int Lx;
     unsigned int Ly;
+    unsigned int Mx;
+    unsigned int My;
     double beta;
     Graph G;
     
@@ -62,15 +48,13 @@ class ma : public hooks {
     uint64_t Nc;
     uint64_t Na;
 
-    std::ofstream bondfile;
-
   public:
     long double lv;
 
 
     std::list<std::list<unsigned int>> avalanches;
 
-    ma(unsigned int Lx, unsigned int Ly, double beta);
+    ma(unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My, double beta);
     ~ma();
 
     void pre_fracture(const network &) override;
diff --git a/src/test.cpp b/src/test.cpp
new file mode 100644
index 0000000..2ebf3ac
--- /dev/null
+++ b/src/test.cpp
@@ -0,0 +1,83 @@
+
+#include <random>
+#include <iostream>
+
+#include <cholmod.h>
+
+#include "randutils/randutils.hpp"
+
+#include <graph.hpp>
+#include <network.hpp>
+#include <hooks.hpp>
+#include "measurements.hpp"
+
+#include <csignal>
+#include <cstring>
+#include <atomic>
+
+std::atomic<bool> quit(false);    // signal flag
+
+void got_signal(int) {
+    quit.store(true);
+}
+
+int main(int argc, char* argv[]) {
+  struct sigaction sa;
+  memset( &sa, 0, sizeof(sa) );
+  sa.sa_handler = got_signal;
+  sigfillset(&sa.sa_mask);
+  sigaction(SIGINT, &sa, NULL);
+
+  int opt;
+
+  unsigned int N = 1;
+  unsigned int Lx = 16;
+  unsigned int Ly = 16;
+  double  beta = 0.5;
+
+  while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
+    switch (opt) {
+      case 'N':
+        N = (unsigned int)atof(optarg);
+        break;
+      case 'X':
+        Lx = atoi(optarg);
+        break;
+      case 'Y':
+        Ly = atoi(optarg);
+        break;
+      case 'b':
+        beta = atof(optarg);
+        break;
+      default:
+        exit(1);
+    }
+  }
+
+  cholmod_common c;
+  CHOL_F(start)(&c);
+
+  randutils::auto_seed_128 seeds;
+  std::mt19937 rng{seeds};
+
+  graph G(Lx, Ly, rng);
+  return 0;
+  /*
+  network base_network(G, &c);
+  ma meas(Lx, Ly, beta);
+
+  for (unsigned int trial = 0; trial < N; trial++) {
+    network tmp_network(base_network);
+    tmp_network.set_thresholds(beta, rng);
+    tmp_network.fracture(meas);
+
+    if (quit.load())
+      break;
+  }
+
+  CHOL_F(finish)(&c);
+
+  return 0;
+  */
+}
+