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

1 files changed, 31 insertions, 129 deletions

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++;
 }