many changes to introduce two-component, elastic-like fracture

15 files changed, 386 insertions, 265 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 35e8436..e3f70d8 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++ -fuse-ld=lld")
-set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -stdlib=libc++ -march=native -fuse-ld=lld")
+set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall")
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -march=native")
 
 set (CMAKE_CXX_STANDARD 17)
 set (CMAKE_C_STANDARD 11)
diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index 1714dc2..45cc478 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -8,6 +8,7 @@
 #include <random>
 #include <list>
 #include <exception>
+#include <algorithm>
 
 #include "array_hash.hpp"
 
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index f12c9c7..ba98086 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -26,25 +26,55 @@
 
 #endif
 
-class network {
+class problem {
   private:
-    cholmod_dense *b;
-    cholmod_factor *factor;
-    cholmod_sparse *voltcurmat;
-    cholmod_common *c;
+    const graph& G;
+    cholmod_dense* b;
+    cholmod_factor* factor;
+    cholmod_sparse* voltcurmat;
+    cholmod_common* c;
+    unsigned axis;
 
   public:
+    problem(const graph&, unsigned axis, cholmod_common*);
+    problem(const graph&, unsigned axis, cholmod_sparse*, cholmod_common*);
+    problem(const problem&);
+    ~problem();
+    current_info solve(std::vector<bool>& fuses);
+    void break_edge(unsigned, bool unbreak = false);
+};
+
+class network {
+  public:
     const graph& G;
     std::vector<bool> fuses;
     std::vector<long double> thresholds;
 
-    network(const graph&, cholmod_common*);
-    network(const network &other);
-    ~network();
+    network(const graph&);
+    network(const network&);
 
     void set_thresholds(double, std::mt19937&);
-    void break_edge(unsigned, bool unbreak = false);
-    current_info get_current_info();
+};
+
+class fuse_network : public network {
+  public:
+    problem ohm;
+
+    fuse_network(const graph&, cholmod_common*);
+
     void fracture(hooks&, double cutoff = 1e-13);
+    current_info get_current_info();
+};
+
+class elastic_network : public network {
+  public:
+    problem hook_x;
+    problem hook_y;
+
+    elastic_network(const graph&, cholmod_common*);
+    elastic_network(const elastic_network&);
+
+    void fracture(hooks&, double weight = 0.5, double cutoff = 1e-13);
+    current_info get_current_info();
 };
 
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 1edc973..8b0e8df 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -17,14 +17,19 @@ class nofuseException: public std::exception
   }
 } nofuseex;
 
-network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) {
+problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), voltcurmat(vcmat), c(c), axis(axis) {
   b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c);
   for (unsigned i = 0; i < G.edges.size(); i++) {
-    double v0y = G.vertices[G.edges[i].v[0]].r.y;
-    double v1y = G.vertices[G.edges[i].v[1]].r.y;
-
-    if (G.edges[i].crossings[1]) {
-      bool ind = v0y < v1y ? 0 : 1;
+    graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
+    graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
+
+    if (G.edges[i].crossings[axis]) {
+      bool ind;
+      if (axis == 1) {
+        ind = v0.y < v1.y ? 0 : 1;
+      } else {
+        ind = v0.x < v1.x ? 0 : 1;
+      }
 
       ((double *)b->x)[G.edges[i].v[ind]] += 1.0;
       ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0;
@@ -78,8 +83,10 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
   factor = CHOL_F(analyze)(laplacian, c);
   CHOL_F(factorize)(laplacian, factor, c);
   CHOL_F(free_sparse)(&laplacian, c);
+}
 
-  t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c);
+problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) {
+  cholmod_triplet *t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c);
 
   t->nnz = 2 * G.edges.size();
 
@@ -98,39 +105,70 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
   CHOL_F(free_triplet)(&t, c);
 }
 
-network::network(const network &other) : c(other.c), G(other.G), fuses(other.fuses), thresholds(other.thresholds) {
+problem::problem(const problem& other) : G(other.G), c(other.c), axis(other.axis) {
   b = CHOL_F(copy_dense)(other.b, c);
   factor = CHOL_F(copy_factor)(other.factor, c);
   voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c);
 } 
 
-network::~network() {
-  CHOL_F(free_sparse)(&voltcurmat, c);
+problem::~problem() {
   CHOL_F(free_dense)(&b, c);
   CHOL_F(free_factor)(&factor, c);
+  CHOL_F(free_sparse)(&voltcurmat, c);
 }
 
-void network::set_thresholds(double beta, std::mt19937& rng) {
-  if (beta == 0.0) {
-    /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */
-    for (long double& threshold : thresholds) {
-      threshold = 1.0;
-    }
-  } else {
-    std::uniform_real_distribution<long double> d(0.0, 1.0);
+current_info problem::solve(std::vector<bool>& fuses) {
+  cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
 
-    for (long double& threshold : thresholds) {
-      threshold = std::numeric_limits<long double>::lowest();
+  if (((double *)x->x)[0] != ((double *)x->x)[0]) {
+    throw nanex;
+  }
 
-      while (threshold == std::numeric_limits<long double>::lowest()) {
-        threshold = logl(d(rng)) / (long double)beta;
+  cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c);
+
+  double alpha[2] = {1, 0};
+  double beta[2] = {0, 0};
+  CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+
+  std::vector<double> currents(G.edges.size());
+
+  double total_current = 0;
+
+  for (int i = 0; i < G.edges.size(); i++) {
+    if (fuses[i]) {
+      currents[i] = 0;
+    } else {
+      currents[i] = ((double *)y->x)[i];
+
+      graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
+      graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
+
+      if (G.edges[i].crossings[axis]) {
+        bool comp;
+        if (axis == 1) {
+          comp = v0.y > v1.y;
+        } else {
+          comp = v0.x > v1.x;
+        }
+
+        if (comp) {
+          currents[i] += 1.0;
+          total_current += currents[i];
+        } else {
+          currents[i] -= 1.0;
+          total_current -= currents[i];
+        }
       }
     }
   }
+
+  CHOL_F(free_dense)(&x, c);
+  CHOL_F(free_dense)(&y, c);
+
+  return {total_current, currents};
 }
 
-void network::break_edge(unsigned e, bool unbreak) {
-  fuses[e] = !unbreak;
+void problem::break_edge(unsigned e, bool unbreak) {
   unsigned v0 = G.edges[e].v[0];
   unsigned v1 = G.edges[e].v[1];
 
@@ -167,78 +205,116 @@ void network::break_edge(unsigned e, bool unbreak) {
   CHOL_F(free_sparse)(&perm_update_mat, c);
   CHOL_F(free_sparse)(&update_mat, c);
 
-  double v0y = G.vertices[v0].r.y;
-  double v1y = G.vertices[v1].r.y;
+  graph::coordinate r0 = G.vertices[v0].r;
+  graph::coordinate r1 = G.vertices[v1].r;
 
-  if (G.edges[e].crossings[1]) {
-    bool ind = v0y < v1y ? unbreak : !unbreak;
+  if (G.edges[e].crossings[axis]) {
+    bool ind;
+    if (axis == 1) {
+      ind = r0.y < r1.y ? unbreak : !unbreak;
+    } else {
+      ind = r0.x < r1.x ? unbreak : !unbreak;
+    }
 
     ((double *)b->x)[G.edges[e].v[ind]] -= 1.0;
     ((double *)b->x)[G.edges[e].v[!ind]] += 1.0;
   }
 }
 
-current_info network::get_current_info() {
-  cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
+network::network(const graph& G) : G(G), fuses(G.edges.size()), thresholds(G.edges.size()) {}
 
-  if (((double *)x->x)[0] != ((double *)x->x)[0]) {
-    throw nanex;
-  }
+network::network(const network& n) : G(n.G), fuses(n.fuses), thresholds(n.thresholds) {}
 
-  cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c);
+void network::set_thresholds(double beta, std::mt19937& rng) {
+  if (beta == 0.0) {
+    /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */
+    for (long double& threshold : thresholds) {
+      threshold = 1.0;
+    }
+  } else {
+    std::uniform_real_distribution<long double> d(0.0, 1.0);
 
-  double alpha[2] = {1, 0};
-  double beta[2] = {0, 0};
-  CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+    for (long double& threshold : thresholds) {
+      threshold = std::numeric_limits<long double>::lowest();
 
-  std::vector<double> currents(G.edges.size());
+      while (threshold == std::numeric_limits<long double>::lowest()) {
+        threshold = logl(d(rng)) / (long double)beta;
+      }
+    }
+  }
+}
 
-  double total_current = 0;
+fuse_network::fuse_network(const graph& G, cholmod_common *c) : network(G), ohm(G, 1, c) {
+}
 
-  for (int i = 0; i < G.edges.size(); i++) {
-    if (fuses[i]) {
-      currents[i] = 0;
-    } else {
-      currents[i] = ((double *)y->x)[i];
+void fuse_network::fracture(hooks& m, double cutoff) {
+  m.pre_fracture(*this);
 
-      double v0y = G.vertices[G.edges[i].v[0]].r.y;
-      double v1y = G.vertices[G.edges[i].v[1]].r.y;
+  while (true) {
+    current_info ci = ohm.solve(fuses);
 
-      if (G.edges[i].crossings[1]) {
-        if (v0y > v1y) {
-          currents[i] += 1.0;
-          total_current += currents[i];
-        } else {
-          currents[i] -= 1.0;
-          total_current -= currents[i];
+    if (ci.conductivity < cutoff * G.vertices.size()) {
+      break;
+    }
+
+    unsigned max_pos = UINT_MAX;
+    long double max_val = std::numeric_limits<long double>::lowest();
+
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      if (!fuses[i] && fabs(ci.currents[i]) > cutoff) {
+        long double val = logl(fabs(ci.currents[i])) - thresholds[i];
+        if (val > max_val) {
+          max_val = val;
+          max_pos = i;
         }
       }
     }
 
+    if (max_pos == UINT_MAX)  {
+      throw nofuseex;
+    }
+
+    fuses[max_pos] = true;
+    ohm.break_edge(max_pos);
+
+    m.bond_broken(*this, ci, max_pos);
   }
 
-  CHOL_F(free_dense)(&x, c);
-  CHOL_F(free_dense)(&y, c);
+  m.post_fracture(*this);
+}
 
-  return {total_current, currents};
+elastic_network::elastic_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) {
+}
+
+elastic_network::elastic_network(const elastic_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) {
 }
 
-void network::fracture(hooks& m, double cutoff) {
+void elastic_network::fracture(hooks& m, double weight, double cutoff) {
   m.pre_fracture(*this);
 
   while (true) {
-    current_info ci = this->get_current_info();
+    current_info cx = hook_x.solve(fuses);
+    current_info cy = hook_y.solve(fuses);
 
-    if (ci.conductivity < cutoff * G.vertices.size()) {
+    current_info ctot;
+
+    ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2));
+
+    if (ctot.conductivity < cutoff * G.vertices.size()) {
       break;
     }
 
+    ctot.currents.resize(G.edges.size());
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i], 2) + pow(weight * cy.currents[i], 2));
+    }
+
     unsigned max_pos = UINT_MAX;
     long double max_val = std::numeric_limits<long double>::lowest();
 
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      if (!fuses[i] && fabs(ci.currents[i]) > cutoff) {
-        long double val = logl(fabs(ci.currents[i])) - thresholds[i];
+      if (!fuses[i] && fabs(ctot.currents[i]) > cutoff) {
+        long double val = logl(fabs(ctot.currents[i])) - thresholds[i];
         if (val > max_val) {
           max_val = val;
           max_pos = i;
@@ -250,9 +326,11 @@ void network::fracture(hooks& m, double cutoff) {
       throw nofuseex;
     }
 
-    this->break_edge(max_pos);
+    fuses[max_pos] = true;
+    hook_x.break_edge(max_pos);
+    hook_y.break_edge(max_pos);
 
-    m.bond_broken(*this, ci, max_pos);
+    m.bond_broken(*this, ctot, max_pos);
   }
 
   m.post_fracture(*this);
diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp
index bc8095e..dea20f0 100644
--- a/src/analysis_tools.cpp
+++ b/src/analysis_tools.cpp
@@ -18,7 +18,7 @@ bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>)
   return true;
 }
 
-std::list<unsigned> find_minimal_crack(const Graph& G, const network& n) {
+std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> find_minimal_crack(const Graph& G, const network& n) {
   Graph Gtmp(n.G.vertices.size());
   std::list<unsigned> removed_edges;
 
@@ -79,53 +79,53 @@ std::list<unsigned> find_minimal_crack(const Graph& G, const network& n) {
     }
   }
 
-  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);
+  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> new_cycle;
-        unsigned pos = 0;
-        for (uint8_t included : new_bool_cycle) {
-          if (included) {
-            new_cycle.push_back(pos);
-          }
-          pos++;
+  // 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> new_cycle;
+      unsigned pos = 0;
+      for (uint8_t included : new_bool_cycle) {
+        if (included) {
+          new_cycle.push_back(pos);
         }
-        cycles.push_back(new_cycle);
+        pos++;
       }
+      cycles.push_back(new_cycle);
     }
+  }
 
-    // find the cycle representing the crack by counting boundary crossings
-    for (auto cycle : cycles) {
-      std::array<unsigned, 2> crossing_count{0,0};
+  std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> output;
 
-      for (auto edge : cycle) {
-        if (n.G.dual_edges[edge].crossings[0]) {
-          crossing_count[0]++;
-        }
-        if (n.G.dual_edges[edge].crossings[1]) {
-          crossing_count[1]++;
-        }
-      }
+  // find the cycle representing the crack by counting boundary crossings
+  for (auto cycle : cycles) {
+    std::array<unsigned, 2> crossing_count{0,0};
 
-      if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
-        return cycle;
+    for (auto edge : cycle) {
+      if (n.G.dual_edges[edge].crossings[0]) {
+        crossing_count[0]++;
+      }
+      if (n.G.dual_edges[edge].crossings[1]) {
+        crossing_count[1]++;
       }
     }
-  } else {
-    return cycles.front();
+
+    if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
+      output.push_back({{1, 0}, cycle});
+    } else if (crossing_count[0] % 2 == 0 && crossing_count[1] % 2 == 1) {
+      output.push_back({{0, 1}, cycle});
+    }
   }
 
-  throw badcycleex;
+  return output;
 }
 
diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp
index 4f3f285..3beb1a8 100644
--- a/src/analysis_tools.hpp
+++ b/src/analysis_tools.hpp
@@ -25,5 +25,5 @@ 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> find_minimal_crack(const Graph &, const network &);
+std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> find_minimal_crack(const Graph &, const network &);
 
diff --git a/src/animate.cpp b/src/animate.cpp
index 5bae15e..be0beef 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -1,5 +1,6 @@
 
 #include "animate.hpp"
+#include <iostream>
 
 animate::animate(double Lx, double Ly, unsigned window_size, int argc, char *argv[]) : G(2 * (unsigned)ceil(Lx * Ly / 2)) {
   glutInit(&argc, argv);
@@ -85,20 +86,41 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
 }
 
 void animate::post_fracture(network &n) {
+
+  std::list<unsigned> crack;
+//  unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
+  unsigned crack_component = 10000;
+
+  std::default_random_engine gen;
+  std::uniform_real_distribution<double> dis(0.0,1.0);
+
+  bool cycle_present = true;
+  auto av_it = avalanches.rbegin();
+
+  while (cycle_present) {
+    for (unsigned e : *av_it) {
+      boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
+      n.fuses[e] = false;
+    }
+
+    auto cracks = find_minimal_crack(G, n);
+    std::cout << cracks.size() << "\n";
+    if (cracks.size() == 0) {
+      cycle_present = false;
+    }
+
+    av_it++;
+  }
+
   std::vector<unsigned> component(boost::num_vertices(G));
   unsigned num = boost::connected_components(G, &component[0]);
 
-  std::list<unsigned> crack = find_minimal_crack(G, n);
-  unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
-
   std::vector<std::list<unsigned>> components(num);
 
   for (unsigned 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') {
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
index 66afc83..c3b4a69 100644
--- a/src/animate_fracture.cpp
+++ b/src/animate_fracture.cpp
@@ -64,9 +64,9 @@ int main(int argc, char* argv[]) {
 
   for (unsigned trial = 0; trial < N; trial++) {
     graph G(Lx, Ly, rng);
-    network network(G, &c);
-    network.set_thresholds(beta, rng);
-    network.fracture(meas);
+    elastic_network elastic_network(G, &c);
+    elastic_network.set_thresholds(beta, rng);
+    elastic_network.fracture(meas);
 
     if (quit.load())
       break;
diff --git a/src/animate_fracture_square.cpp b/src/animate_fracture_square.cpp
index d2efde9..84cc8a5 100644
--- a/src/animate_fracture_square.cpp
+++ b/src/animate_fracture_square.cpp
@@ -63,10 +63,10 @@ int main(int argc, char* argv[]) {
   std::mt19937 rng{seeds};
 
   graph G(Lx, Ly);
-  network perm_network(G, &c);
+  elastic_network perm_network(G, &c);
 
   for (unsigned trial = 0; trial < N; trial++) {
-    network tmp_network(perm_network);
+    elastic_network tmp_network(perm_network);
     tmp_network.set_thresholds(beta, rng);
     tmp_network.fracture(meas);
 
diff --git a/src/fracture.cpp b/src/fracture.cpp
index f697617..fb51d3e 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -82,9 +82,9 @@ int main(int argc, char* argv[]) {
       while (true) {
         try {
           graph G(n, a, rng);
-          network network(G, &c);
-          network.set_thresholds(beta, rng);
-          network.fracture(meas);
+          elastic_network fuse_network(G, &c);
+          fuse_network.set_thresholds(beta, rng);
+          fuse_network.fracture(meas);
           break;
         } catch (std::exception &e) {
           std::cout << e.what() << '\n';
@@ -105,9 +105,9 @@ int main(int argc, char* argv[]) {
       while (true) {
         try {
           graph G(Lx, Ly, rng);
-          network network(G, &c);
-          network.set_thresholds(beta, rng);
-          network.fracture(meas);
+          elastic_network fuse_network(G, &c);
+          fuse_network.set_thresholds(beta, rng);
+          fuse_network.fracture(meas);
           break;
         } catch (std::exception &e) {
           std::cout << e.what() << '\n';
diff --git a/src/fracture_elastic.cpp b/src/fracture_elastic.cpp
new file mode 100644
index 0000000..dec4a2c
--- /dev/null
+++ b/src/fracture_elastic.cpp
@@ -0,0 +1,87 @@
+
+#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 N = 1;
+  unsigned Lx = 16;
+  unsigned Ly = 16;
+  double  beta = 0.5;
+
+  while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
+    switch (opt) {
+      case 'N':
+        N = (unsigned)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);
+
+  ma meas(Lx, Ly, 2*Lx, 2*Ly, beta);
+  graph G(Lx, Ly);
+  network perm_network(G, &c);
+
+  randutils::auto_seed_128 seeds;
+  std::mt19937 rng{seeds};
+
+  for (unsigned trial = 0; trial < N; trial++) {
+    while (true) {
+      try {
+        network tmp_network(perm_network);
+        tmp_network.set_thresholds(beta, rng);
+        tmp_network.fracture(meas);
+        break;
+      } catch (std::exception &e) {
+        std::cout << e.what() << '\n';
+      }
+    }
+
+    if (quit.load())
+      break;
+  }
+
+  CHOL_F(finish)(&c);
+
+  return 0;
+}
+
diff --git a/src/fracture_square.cpp b/src/fracture_square.cpp
index dec4a2c..96e6030 100644
--- a/src/fracture_square.cpp
+++ b/src/fracture_square.cpp
@@ -57,9 +57,9 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  ma meas(Lx, Ly, 2*Lx, 2*Ly, beta);
+  ma meas(Lx, Ly, beta);
   graph G(Lx, Ly);
-  network perm_network(G, &c);
+  elastic_network perm_network(G, &c);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
@@ -67,7 +67,7 @@ int main(int argc, char* argv[]) {
   for (unsigned trial = 0; trial < N; trial++) {
     while (true) {
       try {
-        network tmp_network(perm_network);
+        elastic_network tmp_network(perm_network);
         tmp_network.set_thresholds(beta, rng);
         tmp_network.fracture(meas);
         break;
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 2af1ba3..f2598ee 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -3,15 +3,13 @@
 #include <iostream>
 #include <cstdio>
 
-void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned N, std::string model_string) {
+void update_distribution_file(std::string id, const std::vector<uint64_t>& data, std::string model_string) {
   std::string filename = model_string + id + ".dat";
   std::ifstream file(filename);
 
-  uint64_t N_old = 0;
   std::vector<uint64_t> data_old(data.size(), 0);
 
   if (file.is_open()) {
-    file >> N_old;
     for (unsigned i = 0; i < data.size(); i++) {
       uint64_t num;
       file >> num;
@@ -23,7 +21,6 @@ void update_distribution_file(std::string id, const std::vector<uint64_t>& data,
 
   std::ofstream file_out(filename);
 
-  file_out <<std::fixed<< N_old + N << "\n";
   for (unsigned i = 0; i < data.size(); i++) {
     file_out <<std::fixed<< data_old[i] + data[i] << " ";
   }
@@ -179,122 +176,51 @@ ma::ma(unsigned n, double a, unsigned Mx, unsigned My, double beta) :
   */
 }
 
-ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta) :
-  Mx(Mx), My(My), G(2 * (unsigned)ceil(Lx * Ly / 2)),
-  sc(3 * (unsigned)ceil(Lx * Ly / 2), 0),
-  ss(3 * (unsigned)ceil(Lx * Ly / 2), 0),
-  sm(3 * (unsigned)ceil(Lx * Ly / 2), 0),
-  sa(3 * (unsigned)ceil(Lx * Ly / 2), 0),
-  sl(3 * (unsigned)ceil(Lx * Ly / 2), 0),
-  sb(3 * (unsigned)ceil(Lx * Ly / 2), 0),
-  sD(3 * (unsigned)ceil(Lx * Ly / 2), 0),
-  ccc((Mx / 2) * (My / 2), 0),
-  css((Mx / 2) * (My / 2), 0),
-  cmm((Mx / 2) * (My / 2), 0),
-  caa((Mx / 2) * (My / 2), 0),
-  cll((Mx / 2) * (My / 2), 0),
-  cbb((Mx / 2) * (My / 2), 0),
-  cDD((Mx / 2) * (My / 2), 0),
-  csD((Mx / 2) * (My / 2), 0)
+ma::ma(unsigned Lx, unsigned Ly, double beta) :
+  G(Lx * Ly / 2),
+  sc(Lx * Ly / 2, 0),
+  sm(Lx * Ly / 2, 0),
+  sa(Lx * Ly, 0)
 {
-  N = 0;
-  Nc = 0;
-  Na = 0;
-  Nb = 0;
-
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_";
   } else {
     model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_INF_";
   }
-
-  // FFTW setup for correlation functions
-  /*
-  fftw_set_timelimit(1);
-
-  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));
-
-  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() {
-  // clean up FFTW objects
-  /*
-  fftw_free(fftw_forward_in);
-  fftw_free(fftw_forward_out);
-  fftw_free(fftw_reverse_in);
-  fftw_free(fftw_reverse_out);
-  fftw_destroy_plan(forward_plan);
-  fftw_destroy_plan(reverse_plan);
-  fftw_cleanup();
-  */
-
-  update_distribution_file("sc", sc, Nc, model_string);
-  update_distribution_file("ss", ss, N, model_string);
-  update_distribution_file("sm", sm, N, model_string);
-  update_distribution_file("sa", sa, Na, model_string);
-  update_distribution_file("sl", sl, N, model_string);
-  update_distribution_file("sb", sb, Nb, model_string);
-  update_distribution_file("sD", sD, N, model_string);
-
-  update_field_file("ccc", ccc, Nc, model_string, Mx, My);
-  update_field_file("css", css, N, model_string, Mx, My);
-  update_field_file("cmm", cmm, N, model_string, Mx, My);
-  update_field_file("caa", caa, Na, model_string, Mx, My);
-  update_field_file("cll", cll, N, model_string, Mx, My);
-  update_field_file("cbb", cbb, Nb, model_string, Mx, My);
-  update_field_file("cDD", cDD, N, model_string, Mx, My);
-  update_field_file("csD", csD, N, model_string, Mx, My);
-
-  //stress_file.close();
+  update_distribution_file("sc", sc, model_string);
+  update_distribution_file("sm", sm, model_string);
+  update_distribution_file("sa", sa, model_string);
 }
 
 void ma::pre_fracture(const network&) {
   lv = std::numeric_limits<long double>::lowest();
-  avalanches = {{}};
   boost::remove_edge_if(trivial, G);
+  avalanches = {};
 }
 
 void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
   long double c = logl(cur.conductivity / fabs(cur.currents[i])) + net.thresholds[i];
-  if (c > lv && avalanches.back().size() > 0) {
-    sa[avalanches.back().size() - 1]++;
-    Na++;
-
-    autocorrelation2(net.G.L.x, net.G.L.y, Mx, My, caa, avalanches.back());
-
+  if (c > lv) {
+    if (avalanches.size() > 0) {
+      sa[avalanches.back().size() - 1]++;
+    }
     lv = c;
-    avalanches.push_back({net.G.edges[i].r});
-    last_avalanche = {i};
+    avalanches.push_back({i});
   } else {
-    avalanches.back().push_back(net.G.edges[i].r);
-    last_avalanche.push_back(i);
+    avalanches.back().push_back(i);
   }
 
   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) {
+  auto post_cracks = find_minimal_crack(G, n);
   std::vector<unsigned> component(boost::num_vertices(G));
   unsigned num = boost::connected_components(G, &component[0]);
-
-  std::list<unsigned> crack = find_minimal_crack(G, n);
-
-  ss[crack.size() - 1]++;
-  std::list<graph::coordinate> sr;
-
-  for (auto edge : crack) {
-    sr.push_back(n.G.dual_edges[edge].r);
-  }
-
-  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, css, sr);
-
-  unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
+  unsigned crack_component = component[n.G.dual_edges[post_cracks.front().second.front()].v[0]];
 
   std::vector<std::list<graph::coordinate>> components(num);
 
@@ -303,43 +229,26 @@ void ma::post_fracture(network &n) {
   }
 
   for (unsigned i = 0; i < num; i++) {
-    if (i != crack_component && components[i].size() > 0) {
-      sc[components[i].size() - 1]++;
-      Nc++;
-
-      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ccc, components[i]);
+    if (i != crack_component) {
+      sm[components[i].size() - 1]++;
     }
   }
 
-  // spanning cluster
-
-  sm[components[crack_component].size() - 1]++;
-
-  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cmm, components[crack_component]);
+  auto av_it = avalanches.rbegin();
 
-  /// damage at end
-  std::list<graph::coordinate> Dr;
-
-  for (unsigned i = 0; i < n.G.edges.size(); i++) {
-    if (n.fuses[i]) { 
-      Dr.push_back(n.G.edges[i].r);
+  while (true) {
+    for (unsigned e : *av_it) {
+      boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
+      n.fuses[e] = false;
     }
-  }
-
-  sD[Dr.size() - 1]++;
 
-  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cDD, Dr);
-  correlation2(n.G.L.x, n.G.L.y, Mx, My, csD, sr, Dr);
+    auto cracks = find_minimal_crack(G, n);
 
-  // ********************** LAST AVALANCHE *********************
-
-  // rewind the last avalanche
-  sl[avalanches.back().size() - 1]++;
-
-  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cll, avalanches.back());
+    if (cracks.size() == 0) {
+      break;
+    }
 
-  for (unsigned e : last_avalanche) {
-    boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
+    av_it++;
   }
 
   num = boost::connected_components(G, &component[0]);
@@ -351,14 +260,7 @@ void ma::post_fracture(network &n) {
   }
 
   for (unsigned i = 0; i < num; i++) {
-    if (new_components[i].size() > 0) {
-      sb[new_components[i].size() - 1]++;
-      Nb++;
-
-      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cbb, new_components[i]);
-    }
+    sc[new_components[i].size() - 1]++;
   }
-
-  N++;
 }
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 610c266..7c9d49c 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -35,6 +35,7 @@ class ma : public hooks {
     
     // measurement storage
     std::vector<uint64_t> sc; // non-spanning cluster size distribution
+    std::vector<uint64_t> sn; // non-spanning cluster size distribution
     std::vector<uint64_t> ss; // minimal spanning cluster size distribution
     std::vector<uint64_t> sm; // spanning cluster size distribution
     std::vector<uint64_t> sa; // non-final avalanche size distribution
@@ -56,11 +57,11 @@ class ma : public hooks {
   public:
     long double lv;
 
-    std::list<std::list<graph::coordinate>> avalanches;
+    std::list<std::list<unsigned>> avalanches;
     std::list<unsigned> last_avalanche;
     std::string model_string;
 
-    ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta);
+    ma(unsigned Lx, unsigned Ly, double beta);
     ma(unsigned n, double a, unsigned Mx, unsigned My, double beta);
     ~ma();
 
diff --git a/src/sample_fracture_square.cpp b/src/sample_fracture_square.cpp
index fb2eb33..d4bef2c 100644
--- a/src/sample_fracture_square.cpp
+++ b/src/sample_fracture_square.cpp
@@ -52,10 +52,10 @@ int main(int argc, char* argv[]) {
   std::mt19937 rng{seeds};
 
   graph G(Lx, Ly);
-  network perm_network(G, &c);
+  elastic_network perm_network(G, &c);
 
   for (unsigned trial = 0; trial < N; trial++) {
-    network tmp_network(perm_network);
+    elastic_network tmp_network(perm_network);
     tmp_network.set_thresholds(beta, rng);
     tmp_network.fracture(meas);
   }