simplified and sped up measurementsn substantially, correlation functions and some distributions again incompatible

1 files changed, 61 insertions, 112 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index ff217b2..ed4b8eb 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -37,17 +37,12 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned N, d
   std::ifstream file(filename);
 
   uint64_t N_old = 0;
-  std::vector<std::vector<uint64_t>> data_old(data.size());
-  for (unsigned j = 0; j < data.size(); j++) {
-    data_old[j].resize(data[j].size(), 0);
-  }
+  std::vector<T> data_old(data.size(), 0);
 
   if (file.is_open()) {
     file >> N_old;
     for (unsigned j = 0; j < data.size(); j++) {
-      for (unsigned i = 0; i < data[j].size(); i++) {
-        file >> data_old[j][i];
-      }
+      file >> data_old[j];
     }
     file.close();
   }
@@ -56,10 +51,7 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned N, d
 
   file_out <<std::fixed<< N_old + N <<  "\n";
   for (unsigned j = 0; j < data.size(); j++) {
-    for (unsigned i = 0; i < data[j].size(); i++) {
-      file_out << data_old[j][i] + data[j][i] << " ";
-    }
-    file_out << "\n";
+    file_out << data_old[j] + data[j] << " ";
   }
 
   file_out.close();
@@ -94,6 +86,34 @@ void correlation(unsigned Mx, unsigned My, std::vector<std::vector<T>>& data, co
   }
 }
 
+void autocorrelation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data, const std::list<graph::coordinate>& pos) {
+  for (std::list<graph::coordinate>::const_iterator it1 = pos.begin(); it1 != pos.end(); it1++) {
+    for (std::list<graph::coordinate>::const_iterator it2 = it1; it2 != pos.end(); it2++) {
+      double Δx_tmp = fabs(it1->x - it2->x);
+      double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp;
+
+      double Δy_tmp = fabs(it1->y - it2->y);
+      double Δy = Δy_tmp < Ly / 2 ? Δy_tmp : Ly - Δy_tmp;
+
+      data[floor((Mx * Δx) / Lx) + (Mx / 2) * floor((My * Δy) / Ly)]++;
+    }
+  }
+}
+
+void correlation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data, const std::list<graph::coordinate>& pos1, const std::list<graph::coordinate>& pos2) {
+  for (std::list<graph::coordinate>::const_iterator it1 = pos1.begin(); it1 != pos1.end(); it1++) {
+    for (std::list<graph::coordinate>::const_iterator it2 = pos2.begin(); it2 != pos2.end(); it2++) {
+      double Δx_tmp = fabs(it1->x - it2->x);
+      double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp;
+
+      double Δy_tmp = fabs(it1->y - it2->y);
+      double Δy = Δy_tmp < Ly / 2 ? Δy_tmp : Ly - Δy_tmp;
+
+      data[floor((Mx * Δx) / Lx) + (Mx / 2) * floor((My * Δy) / Ly)]++;
+    }
+  }
+}
+
 template <class T>
 void autocorrelation(unsigned Mx, unsigned My, std::vector<std::vector<T>>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
   fftw_execute(forward_plan);
@@ -116,39 +136,28 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u
   return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);
 }
 
-ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncum) :
+ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta) :
   Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), 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),
-  sd(3 * (unsigned)ceil(Lx * Ly / 2), 0),
   sD(3 * (unsigned)ceil(Lx * Ly / 2), 0),
-  ccc(Ncum),
-  css(Ncum),
-  cmm(Ncum),
-  caa(Ncum),
-  cll(Ncum),
-  cdd(Ncum),
-  cDD(Ncum),
-  csD(Ncum)
+  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),
+  cDD((Mx / 2) * (My / 2), 0),
+  csD((Mx / 2) * (My / 2), 0)
 {
   N = 0;
   Nc = 0;
   Na = 0;
-  for (unsigned i = 0; i < Ncum; i++) {
-    ccc[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
-    css[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
-    cmm[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
-    caa[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
-    cll[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
-    cdd[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
-    cDD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
-    csD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
-  }
 
   // FFTW setup for correlation functions
+  /*
   fftw_set_timelimit(1);
 
   fftw_forward_in = (double *)fftw_malloc(Mx * My * sizeof(double));
@@ -158,12 +167,12 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncu
 
   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);
-  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_current.dat";
-  //stress_file.open(filename, std::ifstream::app);
+  */
 }
 
 ma::~ma() {
   // clean up FFTW objects
+  /*
   fftw_free(fftw_forward_in);
   fftw_free(fftw_forward_out);
   fftw_free(fftw_reverse_in);
@@ -171,13 +180,13 @@ ma::~ma() {
   fftw_destroy_plan(forward_plan);
   fftw_destroy_plan(reverse_plan);
   fftw_cleanup();
+  */
 
   update_distribution_file("sc", sc, Nc, Lx, Ly, beta);
   update_distribution_file("ss", ss, N, Lx, Ly, beta);
   update_distribution_file("sm", sm, N, Lx, Ly, beta);
   update_distribution_file("sa", sa, Na, Lx, Ly, beta);
   update_distribution_file("sl", sl, N, Lx, Ly, beta);
-  update_distribution_file("sd", sd, N, Lx, Ly, beta);
   update_distribution_file("sD", sD, N, Lx, Ly, beta);
 
   update_field_file("ccc", ccc, Nc, Lx, Ly, beta, Mx, My);
@@ -185,7 +194,6 @@ ma::~ma() {
   update_field_file("cmm", cmm, N, Lx, Ly, beta, Mx, My);
   update_field_file("caa", caa, Na, Lx, Ly, beta, Mx, My);
   update_field_file("cll", cll, N, Lx, Ly, beta, Mx, My);
-  update_field_file("cdd", cdd, N, Lx, Ly, beta, Mx, My);
   update_field_file("cDD", cDD, N, Lx, Ly, beta, Mx, My);
   update_field_file("csD", csD, N, Lx, Ly, beta, Mx, My);
 
@@ -204,19 +212,12 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
     sa[avalanches.back().size() - 1]++;
     Na++;
 
-    std::fill_n(fftw_forward_in, Mx * My, 0.0);
-
-    for (auto e : avalanches.back()) {
-      unsigned ind = edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My);
-      fftw_forward_in[ind] += 1.0;
-    }
-
-    autocorrelation(Mx, My, caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+    autocorrelation2(Lx, Ly, Mx, My, caa, avalanches.back());
 
     lv = c;
-    avalanches.push_back({i});
+    avalanches.push_back({net.G.edges[i].r});
   } else {
-    avalanches.back().push_back(i);
+    avalanches.back().push_back(net.G.edges[i].r);
   }
 
   boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
@@ -228,110 +229,58 @@ void ma::post_fracture(network &n) {
 
   std::list<unsigned> crack = find_minimal_crack(G, n);
 
-  // crack surface correlations
-  std::fill_n(fftw_forward_in, Mx * My, 0.0);
-
-  ss[crack.size()]++;
+  ss[crack.size() - 1]++;
+  std::list<graph::coordinate> sr;
 
   for (auto edge : crack) {
-    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My)] += 0.5;
-    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My)] += 0.5;
+    sr.push_back(n.G.dual_edges[edge].r);
   }
 
-  fftw_complex *tss = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);
-
-  correlation(Mx, My, css, tss, tss, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation2(Lx, Ly, Mx, My, css, sr);
 
   unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
 
-  std::vector<std::list<unsigned>> components(num);
+  std::vector<std::list<graph::coordinate>> components(num);
 
   for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
-    components[component[i]].push_back(i);
+    components[component[i]].push_back(n.G.dual_vertices[i].r);
   }
 
-  // non-spanning clusters
-  std::fill_n(fftw_forward_in, Mx * My, 0.0);
-
   for (unsigned i = 0; i < num; i++) {
     if (i != crack_component && components[i].size() > 0) {
       sc[components[i].size() - 1]++;
       Nc++;
 
-      for (auto it = components[i].begin(); it != components[i].end(); it++) {
-        unsigned ind = edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My);
-        fftw_forward_in[ind] += 1.0;
-      }
-
-      autocorrelation(Mx, My, ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-
-      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;
-      }
+      autocorrelation2(Lx, Ly, Mx, My, ccc, components[i]);
     }
   }
 
   // spanning cluster
-  // std::fill_n(fftw_forward_in, Mx * My, 0.0); we already reset in the last loop
 
   sm[components[crack_component].size() - 1]++;
 
-  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);
+  autocorrelation2(Lx, Ly, Mx, My, cmm, components[crack_component]);
 
   /// damage at end
-  std::fill_n(fftw_forward_in, Mx * My, 0.0);
-
-  unsigned final_broken = 0;
+  std::list<graph::coordinate> Dr;
 
   for (unsigned i = 0; i < n.G.edges.size(); i++) {
     if (n.fuses[i]) { 
-      final_broken++;
-      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
+      Dr.push_back(n.G.edges[i].r);
     }
   }
 
-  sD[final_broken]++;
-
-  fftw_complex *tDD = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);
-
-  correlation(Mx, My, cDD, tDD, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-  correlation(Mx, My, csD, tss, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  sD[Dr.size() - 1]++;
 
-  free(tss);
-  free(tDD);
+  autocorrelation2(Lx, Ly, Mx, My, cDD, Dr);
+  correlation2(Lx, Ly, Mx, My, csD, sr, Dr);
 
-  std::fill_n(fftw_forward_in, Mx * My, 0.0);
+  // ********************** LAST AVALANCHE *********************
 
   // rewind the last avalanche
   sl[avalanches.back().size() - 1]++;
 
-  // rewind the last avalanche
-  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);
-    fftw_forward_in[edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0;
-  }
-
-  autocorrelation(Mx, My, cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-
-  // damage size distribution
-  unsigned total_broken = 0;
-
-  std::fill_n(fftw_forward_in, Mx * My, 0.0);
-
-  for (unsigned i = 0; i < n.G.edges.size(); i++) {
-    if (n.fuses[i]) { 
-      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
-    }
-  }
-
-  autocorrelation(Mx, My, cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-
-  sd[total_broken]++;
+  autocorrelation2(Lx, Ly, Mx, My, cll, avalanches.back());
 
   N++;
 }