Added support for measuring avalanches in terms of voltage boundaries, not current boundaries.

2 files changed, 88 insertions, 10 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index 2432198..5132358 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -161,11 +161,11 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u
 }
 
 ma::ma(unsigned n, double a, double beta, double weight)
-    : sc(2 * n), sn(2 * n), ss(2 * n), sm(2 * n), sl(2 * n), sb(n + 1), sd(3 * n), sa(3 * n),
-      sA(3 * n), si(10000), sI(10000), cc(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cn(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
+    : sc(2 * n), sC(2 * n), sn(2 * n), ss(2 * n), sm(2 * n), sl(2 * n), sb(n + 1), sd(3 * n), sD(3 * n), sf(3 * n), sa(3 * n),
+      sA(3 * n), se(3 * n), sE(3 * n), si(10000), sI(10000), cc(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cC(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cn(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
       cs(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cm(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cl(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
-      cb(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), ca(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cA(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
-      cq(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), last_clusters(2 * n) {
+      cb(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), ca(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cA(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), ce(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cE(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
+      cq(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), last_clusters(2 * n), last_clusters_v(2 * n) {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" +
                    std::to_string(beta) + "_" + std::to_string(weight) + "_";
@@ -184,16 +184,20 @@ ma::ma(unsigned n, double a, double beta, double weight)
   Na = 0;
   NA = 0;
   Nq = 0;
+  Ne = 0;
+  NE = 0;
+  NC = 0;
 }
 
 ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
-    : sc(Lx * Ly / 2), sn(Lx * Ly / 2), ss(Lx * Ly / 2), sm(Lx * Ly / 2), sl(Lx * Ly / 2),
-      sb(Lx * Ly / 2 + 1), sd(Lx * Ly), sa(Lx * Ly), sA(Lx * Ly), si(10000), sI(10000),
-      cc((Lx / 2 + 1) * (Ly / 2 + 1)), cn((Lx / 2 + 1) * (Ly / 2 + 1)),
+    : sc(Lx * Ly / 2), sC(Lx * Ly / 2), sn(Lx * Ly / 2), ss(Lx * Ly / 2), sm(Lx * Ly / 2), sl(Lx * Ly / 2),
+      sb(Lx * Ly / 2 + 1), sd(Lx * Ly), sD(Lx * Ly), sf(Lx * Ly), sa(Lx * Ly), sA(Lx * Ly), se(Lx * Ly), sE(Lx * Ly), si(10000), sI(10000),
+      cc((Lx / 2 + 1) * (Ly / 2 + 1)), cC((Lx / 2 + 1) * (Ly / 2 + 1)), cn((Lx / 2 + 1) * (Ly / 2 + 1)),
       cs((Lx / 2 + 1) * (Ly / 2 + 1)), cm((Lx / 2 + 1) * (Ly / 2 + 1)),
       cl((Lx / 2 + 1) * (Ly / 2 + 1)), cb((Lx / 2 + 1) * (Ly / 2 + 1)),
       ca((Lx / 2 + 1) * (Ly / 2 + 1)), cA((Lx / 2 + 1) * (Ly / 2 + 1)),
-      cq((Lx / 2 + 1) * (Ly / 2 + 1)), last_clusters(Lx * Ly / 2) {
+      ce((Lx / 2 + 1) * (Ly / 2 + 1)), cE((Lx / 2 + 1) * (Ly / 2 + 1)),
+      cq((Lx / 2 + 1) * (Ly / 2 + 1)), last_clusters(Lx * Ly / 2), last_clusters_v(Lx * Ly / 2) {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" +
                    std::to_string(beta) + "_" + std::to_string(weight) + "_";
@@ -204,6 +208,7 @@ ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
   Mx = Lx;
   My = Ly;
   Nc = 0;
+  NC = 0;
   Nl = 0;
   Nm = 0;
   Ns = 0;
@@ -211,22 +216,30 @@ ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
   Nb = 0;
   Na = 0;
   NA = 0;
+  Ne = 0;
+  NE = 0;
   Nq = 0;
 }
 
 ma::~ma() {
   update_distribution_file("sc", sc, model_string);
+  update_distribution_file("sC", sC, model_string);
   update_distribution_file("sn", sn, model_string);
   update_distribution_file("ss", ss, model_string);
   update_distribution_file("sm", sm, model_string);
   update_distribution_file("sl", sl, model_string);
   update_distribution_file("sb", sb, model_string);
   update_distribution_file("sd", sd, model_string);
+  update_distribution_file("sD", sD, model_string);
+  update_distribution_file("sf", sf, model_string);
   update_distribution_file("sa", sa, model_string);
   update_distribution_file("sA", sA, model_string);
+  update_distribution_file("se", se, model_string);
+  update_distribution_file("sE", sE, model_string);
   update_distribution_file("si", si, model_string);
   update_distribution_file("sI", sI, model_string);
   update_field_file("cc", Nc, cc, model_string);
+  update_field_file("cC", NC, cC, model_string);
   update_field_file("cl", Nl, cl, model_string);
   update_field_file("cm", Nm, cm, model_string);
   update_field_file("cs", Ns, cs, model_string);
@@ -234,24 +247,36 @@ ma::~ma() {
   update_field_file("cb", Nb, cb, model_string);
   update_field_file("ca", Na, ca, model_string);
   update_field_file("cA", NA, cA, model_string);
+  update_field_file("ce", Ne, ce, model_string);
+  update_field_file("cE", NE, cE, model_string);
   update_field_file("cq", Nq, cq, model_string);
 }
 
 void ma::pre_fracture(const network&) {
+  lc = std::numeric_limits<long double>::lowest();
   lv = std::numeric_limits<long double>::lowest();
   avalanches = {};
+  evalanches = {};
   num = 0;
 }
 
 void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
   long double c = net.thresholds[i] - logl(cur.currents[i]);
-  if (c > lv) {
+  long double v = net.thresholds[i] - logl(cur.currents[i] * cur.conductivity[0]);
+  if (c > lc) {
     last_clusters = net.C;
-    lv = c;
+    lc = c;
     avalanches.push_back({i});
   } else {
     avalanches.back().push_back(i);
   }
+  if (v > lv) {
+    last_clusters_v = net.C;
+    lv = v;
+    evalanches.push_back({i});
+  } else {
+    evalanches.back().push_back(i);
+  }
 
   for (unsigned j = 0; j < cur.currents.size(); j++) {
     if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0 &&
@@ -291,6 +316,20 @@ void ma::post_fracture(network& n) {
   }
   Nc += n.G.dual_vertices.size();
 
+  std::vector<std::list<graph::coordinate>> crit_components_v(n.G.dual_vertices.size());
+
+  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
+    crit_components_v[last_clusters_v.findroot(i)].push_back(n.G.dual_vertices[i].r);
+  }
+
+  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
+    if (crit_components_v[i].size() > 0) {
+      sC[crit_components_v[i].size() - 1]++;
+    }
+    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cC, crit_components_v[i], crack.first);
+  }
+  NC += n.G.dual_vertices.size();
+
     std::vector<std::list<graph::coordinate>> components(n.G.dual_vertices.size());
 
     for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
@@ -355,6 +394,20 @@ void ma::post_fracture(network& n) {
       av_it++;
     }
 
+    auto ev_it = evalanches.rbegin();
+    ev_it++;
+
+    while (ev_it != evalanches.rend()) {
+      se[(*ev_it).size() - 1]++;
+      Ne += (*ev_it).size();
+      std::list<graph::coordinate> ce_co;
+      for (unsigned e : (*ev_it)) {
+        ce_co.push_back(n.G.edges[e].r);
+      }
+      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ce, ce_co, crack.first);
+      ev_it++;
+    }
+
     sA[avalanches.back().size() - 1]++;
     std::list<graph::coordinate> cA_co;
     for (unsigned e : avalanches.back()) {
@@ -363,6 +416,16 @@ void ma::post_fracture(network& n) {
     autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cA, cA_co, crack.first);
     NA += avalanches.back().size();
 
+    sE[evalanches.back().size() - 1]++;
+    std::list<graph::coordinate> cE_co;
+    for (unsigned e : evalanches.back()) {
+      cE_co.push_back(n.G.edges[e].r);
+    }
+    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cE, cE_co, crack.first);
+    NE += evalanches.back().size();
+
     sd[num - 1]++;
+    sD[num - avalanches.back().size()]++;
+    sf[num - evalanches.back().size()]++;
   }
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index ebf79e5..14e6580 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -21,46 +21,61 @@ class ma : public hooks {
     unsigned My;
     unsigned num;
     uint64_t Nc;
+    uint64_t NC;
     uint64_t Nl;
     uint64_t Nm;
     uint64_t Ns;
     uint64_t Nn;
     uint64_t Nb;
     uint64_t Na;
+    uint64_t Ne;
     uint64_t NA;
+    uint64_t NE;
     uint64_t Nq;
     
     // measurement storage
     std::vector<uint64_t> sc;
+    std::vector<uint64_t> sC;
     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> sl; // final avalanche size distribution
     std::vector<uint64_t> sb; // final avalanche size distribution
     std::vector<uint64_t> sd; // final avalanche size distribution
+    std::vector<uint64_t> sD; // final avalanche size distribution
+    std::vector<uint64_t> sf; // final avalanche size distribution
     std::vector<uint64_t> sa; // non-final avalanche size distribution
     std::vector<uint64_t> sA; // non-final avalanche size distribution
+    std::vector<uint64_t> se; // non-final avalanche size distribution
+    std::vector<uint64_t> sE; // non-final avalanche size distribution
 
     std::vector<uint64_t> si;
     std::vector<uint64_t> sI;
 
     std::vector<uint64_t> cc;
+    std::vector<uint64_t> cC;
     std::vector<uint64_t> cn;
     std::vector<uint64_t> cs;
     std::vector<uint64_t> cm;
     std::vector<uint64_t> cl;
     std::vector<uint64_t> cb;
     std::vector<uint64_t> ca;
+    std::vector<uint64_t> ce;
     std::vector<uint64_t> cA;
+    std::vector<uint64_t> cE;
     std::vector<uint64_t> cq;
 
     ClusterTree last_clusters;
+    ClusterTree last_clusters_v;
 
   public:
+    long double lc;
     long double lv;
 
     std::list<std::list<unsigned>> avalanches;
+    std::list<std::list<unsigned>> evalanches;
     std::list<unsigned> last_avalanche;
+    std::list<unsigned> last_evalanche;
     std::string model_string;
 
     ma(unsigned Lx, unsigned Ly, double beta, double w);