cleaned up graph a bit, redefined assert so that certain run errors in jc_voronoi can be caught

1 files changed, 16 insertions, 4 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index 27f239d..15d37d1 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -112,6 +112,8 @@ ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta) :
   Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(2 * (unsigned int)ceil(Lx * Ly / 2)),
   sc(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
   sa(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
+  sC(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
+  sA(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
   sd(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
   sb(log2(Mx < My ? Mx : My) + 1, 0),
   Ccc((Mx / 2 + 1) * (My / 2 + 1), 0),
@@ -125,6 +127,8 @@ ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta) :
   N = 0;
   Nc = 0;
   Na = 0;
+  NC = 0;
+  NA = 0;
 
   // FFTW setup for correlation functions
   fftw_set_timelimit(1);
@@ -150,6 +154,8 @@ ma::~ma() {
 
   update_distribution_file("sa", sa, Na, Lx, Ly, beta);
   update_distribution_file("sc", sc, Nc, Lx, Ly, beta);
+  update_distribution_file("sA", sA, NA, Lx, Ly, beta);
+  update_distribution_file("sC", sC, NC, Lx, Ly, beta);
   update_distribution_file("sd", sd, N, Lx, Ly, beta);
   update_distribution_file("sb", sb, N, Lx, Ly, beta);
 
@@ -172,7 +178,9 @@ void ma::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) {
     sa[avalanches.back().size() - 1]++;
+    sA[avalanches.back().size() - 1]++;
     Na++;
+    NA++;
 
     std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
@@ -224,8 +232,10 @@ void ma::post_fracture(network &n) {
       }
 
       sc[components[i].size() - 1]++;
+      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++;
+      NC++;
 
       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;
@@ -252,10 +262,10 @@ void ma::post_fracture(network &n) {
   // spanning cluster
   std::fill_n(fftw_forward_in, Mx * My, 0.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;
-    }
+  sC[components[crack_component].size() - 1]++;
+  NC++;
+  for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
+    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] += 1.0;
   }
 
   autocorrelation(Mx, My, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
@@ -277,6 +287,8 @@ void ma::post_fracture(network &n) {
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
   // rewind the last avalanche
+  sA[avalanches.back().size() - 1]++;
+  NA++;
   for (auto e : avalanches.back()) {
     boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
     n.break_edge(e, true);