added measurements of spanning cluster sizes and conducting backbone sizes

3 files changed, 68 insertions, 30 deletions

diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 8f7bd12..dd11342 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -54,6 +54,8 @@ class network {
     network(const network&);
 
     void set_thresholds(double, std::mt19937&);
+    virtual void break_edge(unsigned e, bool unbreak = false) {fuses[e] = !unbreak;};
+    virtual current_info get_current_info() {current_info empty; return empty;};
 };
 
 class fuse_network : public network {
@@ -63,10 +65,11 @@ class fuse_network : public network {
     fuse_network(const graph&, cholmod_common*);
 
     void fracture(hooks&, double cutoff = 1e-13);
-    current_info get_current_info();
 };
 
 class elastic_network : public network {
+  private:
+    double weight;
   public:
     problem hook_x;
     problem hook_y;
@@ -75,6 +78,7 @@ class elastic_network : public network {
     elastic_network(const elastic_network&);
 
     void fracture(hooks&, double weight = 0.5, double cutoff = 1e-10);
-    current_info get_current_info();
+    void break_edge(unsigned, bool unbreak = false) override;
+    current_info get_current_info() override;
 };
 
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index e885fcb..ed51c8e 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -290,36 +290,16 @@ elastic_network::elastic_network(const elastic_network& n) : network(n), hook_x(
 }
 
 void elastic_network::fracture(hooks& m, double weight, double cutoff) {
+  this->weight = weight;
   m.pre_fracture(*this);
 
   while (true) {
-    current_info cx = hook_x.solve(fuses);
-    current_info cy = hook_y.solve(fuses);
+    current_info ctot = this->get_current_info();
 
-    bool done_x = cx.conductivity < 1.0 / G.edges.size();
-    bool done_y = cy.conductivity < 1.0 / G.edges.size();
-
-    if (done_x && done_y) {
+    if (ctot.conductivity == 0) {
       break;
     }
 
-    current_info ctot;
-    ctot.currents.resize(G.edges.size());
-
-    if (done_x) {
-      for (unsigned i = 0; i < G.edges.size(); i++) {
-        ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity;
-      }
-    } else if (done_y) {
-      for (unsigned i = 0; i < G.edges.size(); i++) {
-        ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity;
-      }
-    } else {
-      for (unsigned i = 0; i < G.edges.size(); i++) {
-        ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2));
-      }
-    }
-
     unsigned max_pos = UINT_MAX;
     long double max_val = std::numeric_limits<long double>::lowest();
 
@@ -337,9 +317,7 @@ void elastic_network::fracture(hooks& m, double weight, double cutoff) {
       throw nofuseex;
     }
 
-    fuses[max_pos] = true;
-    hook_x.break_edge(max_pos);
-    hook_y.break_edge(max_pos);
+    this->break_edge(max_pos);
 
     m.bond_broken(*this, ctot, max_pos);
   }
@@ -347,3 +325,41 @@ void elastic_network::fracture(hooks& m, double weight, double cutoff) {
   m.post_fracture(*this);
 }
 
+current_info elastic_network::get_current_info() {
+  current_info cx = hook_x.solve(fuses);
+  current_info cy = hook_y.solve(fuses);
+
+  bool done_x = cx.conductivity < 1.0 / G.edges.size();
+  bool done_y = cy.conductivity < 1.0 / G.edges.size();
+
+  current_info ctot;
+  ctot.currents.resize(G.edges.size());
+
+  if (done_x && done_y) {
+    ctot.conductivity = 0;
+  } else if (done_x) {
+    ctot.conductivity = cy.conductivity;
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity;
+    }
+  } else if (done_y) {
+    ctot.conductivity = cx.conductivity;
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity;
+    }
+  } else {
+    ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); 
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2));
+    }
+  }
+
+  return ctot;
+}
+
+void elastic_network::break_edge(unsigned e, bool unbreak) {
+  fuses[e] = !unbreak;
+  hook_x.break_edge(e, unbreak);
+  hook_y.break_edge(e, unbreak);
+}
+
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 407048d..33056ef 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -148,7 +148,9 @@ ma::ma(unsigned n, double a, unsigned Mx, unsigned My, double beta) :
   sm(2 * n),
   sa(3 * n),
   sl(2 * n),
-  sn(2 * n)
+  sn(2 * n),
+  ss(2 * n),
+  sb(3 * n)
 {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_";
@@ -178,6 +180,8 @@ ma::~ma() {
   update_distribution_file("sa", sa, model_string);
   update_distribution_file("sl", sl, model_string);
   update_distribution_file("sn", sn, model_string);
+  update_distribution_file("ss", ss, model_string);
+  update_distribution_file("sb", sb, model_string);
 }
 
 void ma::pre_fracture(const network&) {
@@ -219,6 +223,8 @@ void ma::post_fracture(network &n) {
   for (unsigned i = 0; i < num; i++) {
     if (i != crack_component) {
       sm[components[i].size() - 1]++;
+    } else {
+      ss[components[i].size() - 1]++;
     }
     sn[components[i].size() - 1]++;
   }
@@ -228,7 +234,7 @@ void ma::post_fracture(network &n) {
   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;
+      n.break_edge(e, true);
     }
 
     auto cracks = find_minimal_crack(G, n);
@@ -252,6 +258,18 @@ void ma::post_fracture(network &n) {
     sc[new_components[i].size() - 1]++;
   }
 
+  current_info ct = n.get_current_info();
+
+  unsigned conducting_backbone_size = 0;
+
+  for (unsigned i = 0; i < n.G.edges.size(); i++) {
+    if (ct.currents[i] > 1.0 / n.G.edges.size()) {
+      conducting_backbone_size++;
+    }
+  }
+
+  sb[conducting_backbone_size - 1]++;
+
   av_it++;
 
   while (av_it != avalanches.rend()) {