added support for specifiying voronoi graphs with number of points and aspect ratio, and added support for "infinite" beta (everything has the some threshold)

1 files changed, 73 insertions, 27 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index ed4b8eb..0b570cd 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -3,8 +3,8 @@
 #include <iostream>
 #include <cstdio>
 
-void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned N, double Lx, double Ly, double beta) {
-  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat";
+void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned N, std::string model_string) {
+  std::string filename = model_string + id + ".dat";
   std::ifstream file(filename);
 
   uint64_t N_old = 0;
@@ -32,8 +32,8 @@ void update_distribution_file(std::string id, const std::vector<uint64_t>& data,
 }
 
 template <class T>
-void update_field_file(std::string id, const std::vector<T>& data, unsigned N, double Lx, double Ly, double beta, unsigned Mx, unsigned My) {
-  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat";
+void update_field_file(std::string id, const std::vector<T>& data, unsigned N, std::string model_string, unsigned Mx, unsigned My) {
+  std::string filename = model_string + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat";
   std::ifstream file(filename);
 
   uint64_t N_old = 0;
@@ -95,7 +95,7 @@ void autocorrelation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vecto
       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)]++;
+      data[(unsigned)(Mx * (Δx / Lx)) + (Mx / 2) * (unsigned)(My * (Δy / Ly))]++;
     }
   }
 }
@@ -136,8 +136,48 @@ 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(unsigned n, double a, unsigned Mx, unsigned My, double beta) :
+  Mx(Mx), My(My), G(2 * n),
+  sc(3 * n, 0),
+  ss(3 * n, 0),
+  sm(3 * n, 0),
+  sa(3 * n, 0),
+  sl(3 * n, 0),
+  sD(3 * n, 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),
+  cDD((Mx / 2) * (My / 2), 0),
+  csD((Mx / 2) * (My / 2), 0)
+{
+  N = 0;
+  Nc = 0;
+  Na = 0;
+
+  if (beta != 0.0) {
+    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_";
+  } else {
+    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_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(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)),
+  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),
@@ -156,6 +196,12 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta) :
   Nc = 0;
   Na = 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);
@@ -182,20 +228,20 @@ ma::~ma() {
   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_field_file("ccc", ccc, Nc, Lx, Ly, beta, Mx, My);
-  update_field_file("css", css, N, Lx, Ly, beta, Mx, My);
-  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("csD", csD, N, Lx, Ly, beta, Mx, My);
+  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("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("cDD", cDD, N, model_string, Mx, My);
+  update_field_file("csD", csD, N, model_string, Mx, My);
 
   //stress_file.close();
 }
@@ -212,7 +258,7 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
     sa[avalanches.back().size() - 1]++;
     Na++;
 
-    autocorrelation2(Lx, Ly, Mx, My, caa, avalanches.back());
+    autocorrelation2(net.G.L.x, net.G.L.y, Mx, My, caa, avalanches.back());
 
     lv = c;
     avalanches.push_back({net.G.edges[i].r});
@@ -236,7 +282,7 @@ void ma::post_fracture(network &n) {
     sr.push_back(n.G.dual_edges[edge].r);
   }
 
-  autocorrelation2(Lx, Ly, Mx, My, css, sr);
+  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]];
 
@@ -251,7 +297,7 @@ void ma::post_fracture(network &n) {
       sc[components[i].size() - 1]++;
       Nc++;
 
-      autocorrelation2(Lx, Ly, Mx, My, ccc, components[i]);
+      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ccc, components[i]);
     }
   }
 
@@ -259,7 +305,7 @@ void ma::post_fracture(network &n) {
 
   sm[components[crack_component].size() - 1]++;
 
-  autocorrelation2(Lx, Ly, Mx, My, cmm, components[crack_component]);
+  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cmm, components[crack_component]);
 
   /// damage at end
   std::list<graph::coordinate> Dr;
@@ -272,15 +318,15 @@ void ma::post_fracture(network &n) {
 
   sD[Dr.size() - 1]++;
 
-  autocorrelation2(Lx, Ly, Mx, My, cDD, Dr);
-  correlation2(Lx, Ly, Mx, My, csD, sr, Dr);
+  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);
 
   // ********************** LAST AVALANCHE *********************
 
   // rewind the last avalanche
   sl[avalanches.back().size() - 1]++;
 
-  autocorrelation2(Lx, Ly, Mx, My, cll, avalanches.back());
+  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cll, avalanches.back());
 
   N++;
 }