started implementing crack stress measurements, but probably doesn't work

2 files changed, 159 insertions, 4 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index 2432198..2e45168 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -2,6 +2,7 @@
 #include "measurements.hpp"
 #include <cstdio>
 #include <iostream>
+#include <map>
 
 class badcycleException : public std::exception {
   virtual const char* what() const throw() {
@@ -162,10 +163,10 @@ 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)),
+      sA(3 * n), sp(3 * n), si(10000), sI(10000), 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) {
+      cq(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cS(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), last_clusters(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) + "_";
@@ -188,12 +189,12 @@ ma::ma(unsigned n, double a, double beta, double weight)
 
 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),
+      sb(Lx * Ly / 2 + 1), sd(Lx * Ly), sa(Lx * Ly), sA(Lx * Ly), sp(Lx * Ly), si(10000), sI(10000),
       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) {
+      cq((Lx / 2 + 1) * (Ly / 2 + 1)), cS((Lx / 2 + 1) * (Ly / 2 + 1)), last_clusters(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) + "_";
@@ -212,6 +213,7 @@ ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
   Na = 0;
   NA = 0;
   Nq = 0;
+  NS = 0;
 }
 
 ma::~ma() {
@@ -224,6 +226,7 @@ ma::~ma() {
   update_distribution_file("sd", sd, model_string);
   update_distribution_file("sa", sa, model_string);
   update_distribution_file("sA", sA, model_string);
+  update_distribution_file("sp", sp, model_string);
   update_distribution_file("si", si, model_string);
   update_distribution_file("sI", sI, model_string);
   update_field_file("cc", Nc, cc, model_string);
@@ -235,6 +238,7 @@ ma::~ma() {
   update_field_file("ca", Na, ca, model_string);
   update_field_file("cA", NA, cA, model_string);
   update_field_file("cq", Nq, cq, model_string);
+  update_field_file("cS", NS, cS, model_string);
 }
 
 void ma::pre_fracture(const network&) {
@@ -247,6 +251,8 @@ 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) {
     last_clusters = net.C;
+    last_currents = cur;
+    last_backbone = net.backbone;
     lv = c;
     avalanches.push_back({i});
   } else {
@@ -364,5 +370,147 @@ void ma::post_fracture(network& n) {
     NA += avalanches.back().size();
 
     sd[num - 1]++;
+
+
+    // Currents relative to tip of largest crack fragment
+    std::map<unsigned, std::list<unsigned>> fragments;
+
+    for (unsigned e : crack.second) {
+      unsigned root = last_clusters.findroot(n.G.dual_edges[e].v[0]);
+      if (fragments.contains(root)) {
+        fragments[root].push_back(e);
+      } else {
+        fragments[root] = {e};
+      }
+    }
+
+    if (fragments.size() > 0) {
+      std::list<unsigned> biggest_fragment = std::max_element(fragments.begin(), fragments.end(), [] (const auto& l1, const auto& l2) {return l1.second.size() < l2.second.size();})->second;
+
+      sp[biggest_fragment.size()]++;
+
+      unsigned left_end, right_end;
+
+      bool comp;
+
+      graph::coordinate Δr = {0,0};
+      unsigned e0v0 = n.G.dual_edges[biggest_fragment.front()].v[0];
+      unsigned e0v1 = n.G.dual_edges[biggest_fragment.front()].v[1];
+      unsigned e1v0 = n.G.dual_edges[*(biggest_fragment.begin()++)].v[0];
+      unsigned e1v1 = n.G.dual_edges[*(biggest_fragment.begin()++)].v[1];
+
+      unsigned last_vertex;
+
+      if (e0v0 == e1v0 || e0v0 == e1v1) {
+        last_vertex = e0v1;
+      } else {
+        last_vertex = e1v1;
+      }
+
+
+      for (unsigned e : biggest_fragment) {
+        unsigned v0 = n.G.dual_edges[e].v[0];
+        unsigned v1 = n.G.dual_edges[e].v[1];
+
+        unsigned v_forward, v_back;
+
+        if (v0 == last_vertex) {
+          v_back = v0;
+          v_forward = v1;
+        } else {
+          v_back = v1;
+          v_forward = v0;
+        }
+
+        double Δx = n.G.dual_vertices[v_forward].r.x - n.G.dual_vertices[v_back].r.x;
+
+        if (n.G.dual_edges[e].crossings[0]) {
+          if (Δx > 0) {
+            Δx -= n.G.L.x;
+          } else {
+            Δx += n.G.L.x;
+          }
+        }
+
+        double Δy = n.G.dual_vertices[v_forward].r.y - n.G.dual_vertices[v_back].r.y;
+
+        if (n.G.dual_edges[e].crossings[1]) {
+          if (Δy > 0) {
+            Δy -= n.G.L.y;
+          } else {
+            Δy += n.G.L.y;
+          }
+        }
+
+        Δr.x += Δx;
+        Δr.y += Δy;
+      }
+
+      if (crack.first[0] % 2 == 1) {
+        comp = Δr.x > 0;
+      } else {
+        comp = Δr.y > 0;
+      }
+
+      if (comp) {
+        left_end = biggest_fragment.back();
+        right_end = biggest_fragment.front();
+      } else {
+        left_end = biggest_fragment.front();
+        right_end = biggest_fragment.back();
+      }
+
+      std::cout << n.G.dual_edges[left_end].r.x << " " << n.G.dual_edges[left_end].r.y  << "    " << n.G.dual_edges[right_end].r.x << "   " << n.G.dual_edges[right_end].r.y << "   "<< crack.first[0] % 2 << "     " << Δr.x << " " << Δr.y << "    " << biggest_fragment.size() << "\n";
+
+      for (unsigned i = 0; i < n.G.dual_edges.size(); i++) {
+        if (!last_backbone[i]) {
+          double Δx_tmpl, Δy_tmpl, Δx_tmpr, Δy_tmpr, Δxr, Δxl, Δyr, Δyl;
+
+          if (crack.first[0] % 2 == 1) {
+            Δx_tmpl = n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x;
+            Δxl = Δx_tmpl < 0 ? Δx_tmpl + n.G.L.x : Δx_tmpl;
+
+            Δy_tmpl = fabs(n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y);
+            Δyl = Δy_tmpl > n.G.L.y / 2 ? Δy_tmpl - n.G.L.y / 2 : Δy_tmpl;
+
+            Δx_tmpr = n.G.dual_edges[i].r.x - n.G.dual_edges[left_end].r.x;
+            Δxr = Δx_tmpr < 0 ? Δx_tmpr + n.G.L.x : Δx_tmpr;
+
+            Δy_tmpr = fabs(n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y);
+            Δyr = Δy_tmpr > n.G.L.y / 2 ? Δy_tmpr - n.G.L.y / 2 : Δy_tmpr;
+          } else {
+            Δx_tmpl = n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y;
+            Δxl = Δx_tmpl < 0 ? Δx_tmpl + n.G.L.y : Δx_tmpl;
+
+            Δy_tmpl = fabs(n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x);
+            Δyl = Δy_tmpl > n.G.L.x / 2 ? Δy_tmpl - n.G.L.x / 2 : Δy_tmpl;
+
+            Δx_tmpr = n.G.dual_edges[i].r.y - n.G.dual_edges[left_end].r.y;
+            Δxr = Δx_tmpr < 0 ? Δx_tmpr + n.G.L.y : Δx_tmpr;
+
+            Δy_tmpr = fabs(n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x);
+            Δyr = Δy_tmpr > n.G.L.x / 2 ? Δy_tmpr - n.G.L.x / 2 : Δy_tmpr;
+          }
+
+          if (Δxl < Δxr) {
+            cS[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δyl / n.G.L.y)) +
+               (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δxl / n.G.L.x))] +=
+                last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
+          } else {
+            cS[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δyr / n.G.L.y)) +
+               (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δxr / n.G.L.x))] +=
+                last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
+          }
+        }
+      }
+
+      NS++;
+
+    } else {
+      sp[0]++;
+    }
+
+    aG.push_back(last_currents.conductivity);
+
   }
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index ebf79e5..c0458b5 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -29,6 +29,7 @@ class ma : public hooks {
     uint64_t Na;
     uint64_t NA;
     uint64_t Nq;
+    uint64_t NS;
     
     // measurement storage
     std::vector<uint64_t> sc;
@@ -54,7 +55,13 @@ class ma : public hooks {
     std::vector<uint64_t> cA;
     std::vector<uint64_t> cq;
 
+    std::list<std::array<double, 2>> aG;
+    std::vector<double> cS;
+    std::vector<uint64_t> sp;
+
     ClusterTree last_clusters;
+    current_info last_currents;
+    std::vector<bool> last_backbone;
 
   public:
     long double lv;