now takes Lx and Ly for anisotropic fracture

6 files changed, 93 insertions, 85 deletions

diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index 80cdd66..9a630ff 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -18,12 +18,14 @@ class graph {
 
     typedef std::array<unsigned int, 2> edge;
 
+    coordinate L;
+
     std::vector<vertex> vertices;
     std::vector<edge> edges;
 
     std::vector<vertex> dual_vertices;
     std::vector<edge> dual_edges;
 
-    graph(unsigned int L);
+    graph(unsigned int Nx, unsigned int Ny);
 };
 
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index e5f374e..f2e8065 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -1,10 +1,11 @@
 
 #include "graph.hpp"
 
-graph::graph(unsigned int L) {
-  double dx = 1.0 / L;
-  unsigned int nv = 2 * pow(L / 2, 2);
-  unsigned int ne = pow(L, 2);
+graph::graph(unsigned int Nx, unsigned int Ny) {
+  L = {(double)Nx, (double)Ny};
+
+  unsigned int ne = Nx * Ny;
+  unsigned int nv = ne / 2;
 
   vertices.resize(nv);
   edges.reserve(ne);
@@ -13,22 +14,22 @@ graph::graph(unsigned int L) {
   dual_edges.reserve(ne);
 
   for (unsigned int i = 0; i < nv; i++) {
-    vertices[i].r.x = dx * ((1 + i / (L / 2)) % 2 + 2 * (i % (L / 2)));
-    vertices[i].r.y = dx * (i / (L / 2));
+    vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
+    vertices[i].r.y = (double)(i / (Nx / 2));
 
-    dual_vertices[i].r.x = dx * ((i / (L / 2)) % 2 + 2 * (i % (L / 2)));
-    dual_vertices[i].r.y = dx * (i / (L / 2));
+    dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
+    dual_vertices[i].r.y = (double)(i / (Nx / 2));
   }
 
-  for (unsigned int x = 0; x < L; x++) {
-    for (unsigned int y = 0; y < L; y++) {
-      unsigned int v1 = (L * x) / 2 + ((y + x % 2) / 2) % (L / 2);
-      unsigned int v2 = ((L * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2)) % nv;
+  for (unsigned int x = 0; x < Ny; x++) {
+    for (unsigned int y = 0; y < Nx; y++) {
+      unsigned int v1 = (Nx * x) / 2 + ((y + x % 2) / 2) % (Nx / 2);
+      unsigned int v2 = ((Nx * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2)) % nv;
 
       edges.push_back({v1, v2});
 
-      unsigned int dv1 = (L * x) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2);
-      unsigned int dv2 = ((L * (x + 1)) / 2 + ((y + x % 2) / 2) % (L / 2)) % nv;
+      unsigned int dv1 = (Nx * x) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2);
+      unsigned int dv2 = ((Nx * (x + 1)) / 2 + ((y + x % 2) / 2) % (Nx / 2)) % nv;
 
       dual_edges.push_back({dv1, dv2});
     }
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 9f6cb0d..bc7a0c6 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -7,7 +7,7 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
     double v0y = G.vertices[G.edges[i][0]].r.y;
     double v1y = G.vertices[G.edges[i][1]].r.y;
 
-    if (fabs(v0y - v1y) > 0.5) {
+    if (fabs(v0y - v1y) > G.L.y / 2) {
       bool ind = v0y < v1y ? 0 : 1;
 
       ((double *)b->x)[G.edges[i][ind]] += 1.0;
@@ -138,7 +138,7 @@ void network::break_edge(unsigned int e, bool unbreak) {
   double v0y = G.vertices[v0].r.y;
   double v1y = G.vertices[v1].r.y;
 
-  if (fabs(v0y - v1y) > 0.5) {
+  if (fabs(v0y - v1y) > G.L.y / 2) {
     bool ind = v0y < v1y ? unbreak : !unbreak;
 
     ((double *)b->x)[G.edges[e][ind]] -= 1.0;
@@ -172,7 +172,7 @@ current_info network::get_current_info() {
       double v0y = G.vertices[G.edges[i][0]].r.y;
       double v1y = G.vertices[G.edges[i][1]].r.y;
 
-      if (fabs(v0y - v1y) > 0.5) {
+      if (fabs(v0y - v1y) > G.L.y / 2) {
         if (v0y > v1y) {
           currents[i] += 1.0;
           total_current += currents[i];
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 0aad56e..317d91b 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -31,16 +31,20 @@ int main(int argc, char* argv[]) {
   int opt;
 
   unsigned int N = 1;
-  unsigned int L = 16;
+  unsigned int Lx = 16;
+  unsigned int Ly = 16;
   double  beta = 0.5;
 
-  while ((opt = getopt(argc, argv, "N:L:b:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
     switch (opt) {
       case 'N':
         N = (unsigned int)atof(optarg);
         break;
-      case 'L':
-        L = atoi(optarg);
+      case 'X':
+        Lx = atoi(optarg);
+        break;
+      case 'Y':
+        Ly = atoi(optarg);
         break;
       case 'b':
         beta = atof(optarg);
@@ -53,9 +57,9 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  graph G(L);
+  graph G(Lx, Ly);
   network base_network(G, &c);
-  ma meas(L, beta);
+  ma meas(Lx, Ly, beta);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 57f85f5..e32f07f 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -104,11 +104,11 @@ std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
 
       for (auto edge : cycle) {
         double dx = fabs(n.G.vertices[n.G.dual_edges[edge][0]].r.x - n.G.vertices[n.G.dual_edges[edge][1]].r.x);
-        if (dx > 0.5) {
+        if (dx > n.G.L.x / 2) {
           crossing_count[0]++;
         }
         double dy = fabs(n.G.vertices[n.G.dual_edges[edge][0]].r.y - n.G.vertices[n.G.dual_edges[edge][1]].r.y);
-        if (dy > 0.5) {
+        if (dy > n.G.L.y / 2) {
           crossing_count[1]++;
         }
       }
@@ -124,8 +124,8 @@ std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
   exit(5);
 }
 
-void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, unsigned int L, double beta) {
-  std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat";
+void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta) {
+  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat";
   std::ifstream file(filename);
 
   uint64_t N_old = 0;
@@ -153,8 +153,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 int N, unsigned int L, double beta) {
-  std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat";
+void update_field_file(std::string id, const std::vector<T>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta) {
+  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat";
   std::ifstream file(filename);
 
   uint64_t N_old = 0;
@@ -179,16 +179,16 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned int
 }
 
 template <class T>
-std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<T>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) {
-  for (unsigned int i = 0; i < pow(L, 2); i++) {
+std::vector<fftw_complex> data_transform(unsigned int Lx, unsigned int Ly, const std::vector<T>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) {
+  for (unsigned int i = 0; i < Lx * Ly; i++) {
     fftw_forward_in[i] = (double)data[i];
   }
 
   fftw_execute(forward_plan);
 
-  std::vector<fftw_complex> output(L * (L / 2 + 1));
+  std::vector<fftw_complex> output(Lx * (Ly / 2 + 1));
 
-  for (unsigned int i = 0; i < L * (L / 2 + 1); i++) {
+  for (unsigned int i = 0; i < Lx * (Ly / 2 + 1); i++) {
     output[i][0] = fftw_forward_out[i][0];
     output[i][1] = fftw_forward_out[i][1];
   }
@@ -197,47 +197,47 @@ std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<T>& d
 }
 
 template <class T>
-void correlation(unsigned int L, std::vector<T>& data, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
-  for (unsigned int i = 0; i < L * (L / 2 + 1); i++) {
+void correlation(unsigned int Lx, unsigned int Ly, std::vector<T>& data, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
+  for (unsigned int i = 0; i < Lx * (Ly / 2 + 1); i++) {
     fftw_reverse_in[i][0] = tx1[i][0] * tx2[i][0] + tx1[i][1] * tx2[i][1];
     fftw_reverse_in[i][1] = tx1[i][0] * tx2[i][1] - tx1[i][1] * tx2[i][0];
   }
 
   fftw_execute(reverse_plan);
 
-  for (unsigned int i = 0; i < pow(L / 2 + 1, 2); i++) {
-    data[i] += (T)(fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2));
+  for (unsigned int i = 0; i < (Lx / 2 + 1) * (Ly / 2 + 1); i++) {
+    data[i] += (T)(fftw_reverse_out[Lx * (i / (Lx / 2 + 1)) + i % (Lx / 2 + 1)] / (Lx * Ly));
   }
 }
 
 template <class T>
-void autocorrelation(unsigned int L, 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) {
+void autocorrelation(unsigned int Lx, unsigned int Ly, 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);
 
-  for (unsigned int i = 0; i < L * (L / 2 + 1); i++) {
+  for (unsigned int i = 0; i < Ly * (Lx / 2 + 1); i++) {
     fftw_reverse_in[i][0] = pow(fftw_forward_out[i][0], 2) + pow(fftw_forward_out[i][1], 2);
     fftw_reverse_in[i][1] = 0.0;
   }
 
   fftw_execute(reverse_plan);
 
-  for (unsigned int i = 0; i < pow(L / 2 + 1, 2); i++) {
-    out_data[i] += (T)(fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2));
+  for (unsigned int i = 0; i < (Lx / 2 + 1) * (Ly / 2 + 1); i++) {
+    out_data[i] += (T)(fftw_reverse_out[Lx * (i / (Lx / 2 + 1)) + i % (Lx / 2 + 1)] / (Lx * Ly));
   }
 }
 
-ma::ma(unsigned int L, double beta) : L(L), beta(beta), G(2 * pow(L / 2, 2)),
-  sc(pow(L, 2), 0),
-  sa(pow(L, 2), 0),
-  sd(pow(L, 2), 0),
-  sb(log2(L) + 1, 0),
-  Ccc(pow(L / 2 + 1, 2), 0),
-  Css(pow(L / 2 + 1, 2), 0),
-  Cmm(pow(L / 2 + 1, 2), 0),
-  Caa(pow(L / 2 + 1, 2), 0),
-  Cdd(pow(L / 2 + 1, 2), 0),
-  Cll(pow(L / 2 + 1, 2), 0),
-  Cee(pow(L / 2 + 1, 2), 0)
+ma::ma(unsigned int Lx, unsigned int Ly, double beta) : Lx(Lx), Ly(Ly), beta(beta), G(Lx * Ly / 2),
+  sc(Lx * Ly, 0),
+  sa(Lx * Ly, 0),
+  sd(Lx * Ly, 0),
+  sb(log2(Ly) + 1, 0),
+  Ccc((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+  Css((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+  Cmm((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+  Caa((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+  Cdd((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+  Cll((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+  Cee((Lx / 2 + 1) * (Ly / 2 + 1), 0)
 {
   N = 0;
   Nc = 0;
@@ -246,13 +246,13 @@ ma::ma(unsigned int L, double beta) : L(L), beta(beta), G(2 * pow(L / 2, 2)),
   // FFTW setup for correlation functions
   fftw_set_timelimit(1);
 
-  fftw_forward_in = (double *)fftw_malloc(pow(L, 2) * sizeof(double));
-  fftw_forward_out = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex));
-  fftw_reverse_in = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex));
-  fftw_reverse_out = (double *)fftw_malloc(pow(L, 2) * sizeof(double));
+  fftw_forward_in = (double *)fftw_malloc(Lx * Ly * sizeof(double));
+  fftw_forward_out = (fftw_complex *)fftw_malloc(Lx * Ly * sizeof(fftw_complex));
+  fftw_reverse_in = (fftw_complex *)fftw_malloc(Lx * Ly * sizeof(fftw_complex));
+  fftw_reverse_out = (double *)fftw_malloc(Lx * Ly * sizeof(double));
 
-  forward_plan = fftw_plan_dft_r2c_2d(L, L, fftw_forward_in, fftw_forward_out, 0);
-  reverse_plan = fftw_plan_dft_c2r_2d(L, L, fftw_reverse_in, fftw_reverse_out, 0);
+  forward_plan = fftw_plan_dft_r2c_2d(Ly, Lx, fftw_forward_in, fftw_forward_out, 0);
+  reverse_plan = fftw_plan_dft_c2r_2d(Ly, Lx, fftw_reverse_in, fftw_reverse_out, 0);
 }
 
 ma::~ma() {
@@ -265,18 +265,18 @@ ma::~ma() {
   fftw_destroy_plan(reverse_plan);
   fftw_cleanup();
 
-  update_distribution_file("sa", sa, Na, L, beta);
-  update_distribution_file("sc", sc, Nc, L, beta);
-  update_distribution_file("sd", sd, N, L, beta);
-  update_distribution_file("sb", sb, N, L, beta);
-
-  update_field_file("Ccc", Ccc, Nc, L, beta);
-  update_field_file("Css", Css, N, L, beta);
-  update_field_file("Cmm", Cmm, N, L, beta);
-  update_field_file("Cdd", Cdd, N, L, beta);
-  update_field_file("Caa", Caa, Na, L, beta);
-  update_field_file("Cll", Cll, N, L, beta);
-  update_field_file("Cee", Cee, N, L, 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);
+
+  update_field_file("Ccc", Ccc, Nc, Lx, Ly, beta);
+  update_field_file("Css", Css, N, Lx, Ly, beta);
+  update_field_file("Cmm", Cmm, N, Lx, Ly, beta);
+  update_field_file("Cdd", Cdd, N, Lx, Ly, beta);
+  update_field_file("Caa", Caa, Na, Lx, Ly, beta);
+  update_field_file("Cll", Cll, N, Lx, Ly, beta);
+  update_field_file("Cee", Cee, N, Lx, Ly, beta);
 }
 
 void ma::pre_fracture(const network &) {
@@ -297,7 +297,7 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned int i
       fftw_forward_in[e] = 1.0;
     }
 
-    autocorrelation(L, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+    autocorrelation(Lx, Ly, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
     lv = c;
     avalanches.push_back({i});
@@ -332,20 +332,20 @@ void ma::post_fracture(network &n) {
 
       if (size > 0) {
         sc[size - 1]++;
-        autocorrelation(L, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+        autocorrelation(Lx, Ly, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
         Nc++;
       }
     }
   }
 
   // bin counting
-  for (unsigned int be = 0; be < log2(L); be++) {
+  for (unsigned int be = 0; be < log2(Ly); be++) {
     unsigned int bin = pow(2, be);
 
-    for (unsigned int i = 0; i < pow(L / bin, 2); i++) {
+    for (unsigned int i = 0; i < Lx * Ly / pow(bin, 2); i++) {
       bool in_bin = false;
       for (unsigned int j = 0; j < pow(bin, 2); j++) {
-        unsigned int edge = L * (bin * ((i * bin) / L) + j / bin) + (i * bin) % L + j % bin;
+        unsigned int edge = Lx * (bin * ((i * bin) / Lx) + j / bin) + (i * bin) % Lx + j % bin;
         if (!n.fuses[edge] && crack_component == component[n.G.dual_edges[edge][0]]) {
           in_bin = true;
           break;
@@ -358,7 +358,7 @@ void ma::post_fracture(network &n) {
     }
   }
 
-  sb[log2(L)]++;
+  sb[log2(Ly)]++;
 
   for (unsigned int i = 0; i < n.G.edges.size(); i++) {
     if (n.fuses[i] && component[n.G.dual_edges[i][0]] == crack_component)  {
@@ -368,7 +368,7 @@ void ma::post_fracture(network &n) {
     }
   }
 
-  autocorrelation(L, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Lx, Ly, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   // crack surface correlations
   std::fill_n(fftw_forward_in, n.G.edges.size(), 0.0);
@@ -377,7 +377,7 @@ void ma::post_fracture(network &n) {
     fftw_forward_in[edge] = 1.0;
   }
 
-  autocorrelation(L, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Lx, Ly, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   std::function<bool(unsigned int)> inCrack = [&](unsigned int i) -> bool {
     return (bool)fftw_forward_in[i];
@@ -391,7 +391,7 @@ void ma::post_fracture(network &n) {
     }
   }
 
-  autocorrelation(L, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Lx, Ly, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   std::fill_n(fftw_forward_in, n.G.edges.size(), 0.0);
 
@@ -402,7 +402,7 @@ void ma::post_fracture(network &n) {
     fftw_forward_in[e] = 1;
   }
 
-  autocorrelation(L, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Lx, Ly, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   // damage size distribution
   unsigned int total_broken = 0;
@@ -416,7 +416,7 @@ void ma::post_fracture(network &n) {
     }
   }
 
-  autocorrelation(L, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Lx, Ly, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   sd[total_broken]++;
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index ab31ccd..703dccd 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -42,7 +42,8 @@ class ma : public hooks {
     fftw_plan forward_plan;
     fftw_plan reverse_plan;
     unsigned int N;
-    unsigned int L;
+    unsigned int Lx;
+    unsigned int Ly;
     double beta;
     Graph G;
     
@@ -67,7 +68,7 @@ class ma : public hooks {
 
     std::list<std::list<unsigned int>> avalanches;
 
-    ma(unsigned int L, double beta);
+    ma(unsigned int Lx, unsigned int Ly, double beta);
     ~ma();
 
     void pre_fracture(const network &) override;