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

6 files changed, 161 insertions, 60 deletions

diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index 1089925..032a4d1 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -39,5 +39,8 @@ class graph {
 
     graph(unsigned Nx, unsigned Ny);
     graph(double Lx, double Ly, std::mt19937& rng);
+    graph(unsigned n, double a, std::mt19937& rng);
+
+    void helper(unsigned n, std::mt19937& rng);
 };
 
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 526ffeb..8aa8877 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -141,27 +141,38 @@ unsigned get_triangle_signature(unsigned j1, unsigned j2, unsigned j3) {
 }
 
 graph::graph(double Lx, double Ly, std::mt19937& rng) {
-  L = {Lx, Ly};
-
   // randomly choose N to be floor(Lx * Ly / 2) or ceil(Lx * Ly / 2) with
   // probability proportional to the distance from each
   std::uniform_real_distribution<double> d(0.0, 1.0);
   unsigned N = round(Lx * Ly / 2 + d(rng) - 0.5);
 
-  unsigned nv = N;
+  L = {Lx, Ly};
+
+  this->helper(N, rng);
+}
+
+graph::graph(unsigned n, double a, std::mt19937& rng) {
+  L = {2 * n * a, 2 * n / a};
+
+  this->helper(n, rng);
+}
+
+void graph::helper(unsigned n, std::mt19937& rng) {
+  std::uniform_real_distribution<double> d(0.0, 1.0);
+  unsigned nv = n;
   vertices.resize(nv);
 
   // the coordinates of the lattice, from which a delaunay triangulation
   // and corresponding voronoi tessellation will be built. Everyone is in the
   // rectangle (0, 0) (Lx, Ly)
   for (vertex &v : vertices) {
-    v = {{Lx * d(rng), Ly * d(rng)}};
+    v = {{L.x * d(rng), L.y * d(rng)}};
   }
 
   // set up the voronoi objects
   jcv_diagram diagram;
   memset(&diagram, 0, sizeof(jcv_diagram));
-  jcv_rect bounds = {{-Lx, -Ly}, {2 * Lx, 2 * Ly}};
+  jcv_rect bounds = {{-L.x, -L.y}, {2 * L.x, 2 * L.y}};
   std::vector<jcv_point> points(9 * nv);
 
   for (unsigned i = 0; i < nv; i++) {
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 027946a..43b0fab 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -95,13 +95,20 @@ network::~network() {
 }
 
 void network::set_thresholds(double beta, std::mt19937& rng) {
-  std::uniform_real_distribution<long double> d(0.0, 1.0);
+  if (beta == 0.0) {
+    /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */
+    for (long double& threshold : thresholds) {
+      threshold = 1.0;
+    }
+  } else {
+    std::uniform_real_distribution<long double> d(0.0, 1.0);
 
-  for (long double& threshold : thresholds) {
-    threshold = std::numeric_limits<long double>::lowest();
+    for (long double& threshold : thresholds) {
+      threshold = std::numeric_limits<long double>::lowest();
 
-    while (threshold == std::numeric_limits<long double>::lowest()) {
-      threshold = logl(d(rng)) / (long double)beta;
+      while (threshold == std::numeric_limits<long double>::lowest()) {
+        threshold = logl(d(rng)) / (long double)beta;
+      }
     }
   }
 }
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 53a11a5..b1bde45 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -35,7 +35,11 @@ int main(int argc, char* argv[]) {
   double Ly = 16;
   double  beta = 0.5;
 
-  while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
+  unsigned n = 128;
+  double a = 1.0;
+  bool use_aN = false;
+
+  while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:")) != -1) {
     switch (opt) {
       case 'N':
         N = (unsigned)atof(optarg);
@@ -49,6 +53,14 @@ int main(int argc, char* argv[]) {
       case 'b':
         beta = atof(optarg);
         break;
+      case 'n':
+        n = (unsigned)atof(optarg);
+        use_aN = true;
+        break;
+      case 'a':
+        a = atof(optarg);
+        use_aN = true;
+        break;
       default:
         exit(1);
     }
@@ -57,30 +69,54 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  // fourier transforms of powers of two are faster
-  unsigned Mx = pow(2, ceil(log2(4*Lx)));
-  unsigned My = pow(2, ceil(log2(4*Ly)));
-
-  ma meas(Lx, Ly, Mx, My, beta);
-
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
 
-  for (unsigned trial = 0; trial < N; trial++) {
-    while (true) {
-      try {
-        graph G(Lx, Ly, rng);
-        network network(G, &c);
-        network.set_thresholds(beta, rng);
-        network.fracture(meas);
-        break;
-      } catch (std::exception &e) {
-        std::cout << e.what() << '\n';
+  if (use_aN) {
+    unsigned Mx = pow(2, ceil(log2(4*2*n*a)));
+    unsigned My = pow(2, ceil(log2(4*2*n/a)));
+
+    ma meas(n, a, Mx, My, beta);
+
+    for (unsigned trial = 0; trial < N; trial++) {
+      while (true) {
+        try {
+          graph G(n, a, rng);
+          network network(G, &c);
+          network.set_thresholds(beta, rng);
+          network.fracture(meas);
+          break;
+        } catch (std::exception &e) {
+          std::cout << e.what() << '\n';
+        }
       }
+
+      if (quit.load())
+        break;
     }
+  } else {
+    // fourier transforms of powers of two are faster
+    unsigned Mx = pow(2, ceil(log2(4*Lx)));
+    unsigned My = pow(2, ceil(log2(4*Ly)));
+
+    ma meas(Lx, Ly, Mx, My, beta);
+
+    for (unsigned trial = 0; trial < N; trial++) {
+      while (true) {
+        try {
+          graph G(Lx, Ly, rng);
+          network network(G, &c);
+          network.set_thresholds(beta, rng);
+          network.fracture(meas);
+          break;
+        } catch (std::exception &e) {
+          std::cout << e.what() << '\n';
+        }
+      }
 
-    if (quit.load())
-      break;
+      if (quit.load())
+        break;
+    }
   }
 
   CHOL_F(finish)(&c);
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++;
 }
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 0fd4644..600c6ac 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -28,11 +28,8 @@ class ma : public hooks {
     fftw_plan reverse_plan;
     */
     unsigned N;
-    double Lx;
-    double Ly;
     unsigned Mx;
     unsigned My;
-    double beta;
     Graph G;
 //    std::ofstream stress_file;
     
@@ -56,10 +53,11 @@ class ma : public hooks {
   public:
     long double lv;
 
-
     std::list<std::list<graph::coordinate>> avalanches;
+    std::string model_string;
 
     ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta);
+    ma(unsigned n, double a, unsigned Mx, unsigned My, double beta);
     ~ma();
 
     void pre_fracture(const network &) override;