added animation example, and did many fixes to the voronoi system

1 files changed, 87 insertions, 216 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index c8cc73c..43483be 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -2,128 +2,6 @@
 #include "measurements.hpp"
 #include <iostream>
 
-template <class T>
-bool is_shorter(std::list<T> l1, std::list<T> l2) {
-  return l1.size() < l2.size();
-}
-
-bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) {
-  return true;
-}
-
-std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
-  Graph Gtmp(n.G.vertices.size());
-  std::list<unsigned int> removed_edges;
-
-  class add_tree_edges : public boost::default_dfs_visitor {
-    public:
-      Graph& G;
-      std::list<unsigned int>& E;
-
-      add_tree_edges(Graph& G, std::list<unsigned int>& E) : G(G), E(E) {}
-
-      void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        boost::add_edge(boost::source(e, g), boost::target(e, g), g[e], G);
-      }
-
-      void back_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        if (!(boost::edge(boost::source(e, g), boost::target(e, g), G).second)) {
-          E.push_back(g[e].index);
-        }
-      }
-  };
-
-  add_tree_edges ate(Gtmp, removed_edges);
-  boost::depth_first_search(G, visitor(ate));
-
-  class find_cycle : public boost::default_dfs_visitor {
-    public:
-      std::list<unsigned int>& E;
-      unsigned int end;
-      struct done{};
-
-      find_cycle(std::list<unsigned int>& E, unsigned int end) : E(E), end(end) {}
-
-      void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
-        if (v == end) {
-          throw done{};
-        }
-      }
-
-      void examine_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        E.push_back(g[e].index);
-      }
-
-      void finish_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
-        E.erase(std::find(E.begin(), E.end(), g[e].index));
-      }
-  };
-
-  std::list<std::list<unsigned int>> cycles;
-
-  for (auto edge : removed_edges) {
-    std::list<unsigned int> cycle = {edge};
-    find_cycle vis(cycle, n.G.dual_edges[edge][1]);
-    std::vector<boost::default_color_type> new_color_map(boost::num_vertices(Gtmp));
-    try {
-    boost::depth_first_visit(Gtmp, n.G.dual_edges[edge][0], vis, boost::make_iterator_property_map(new_color_map.begin(), boost::get(boost::vertex_index, Gtmp), new_color_map[0]));
-    } catch(find_cycle::done const&) {
-      cycles.push_back(cycle);
-    }
-  }
-
-  if (cycles.size() > 1) {
-    std::list<std::valarray<uint8_t>> bool_cycles;
-    for (auto cycle : cycles) {
-      std::valarray<uint8_t> bool_cycle(n.G.edges.size());
-      for (auto v : cycle) {
-        bool_cycle[v] = 1;
-      }
-      bool_cycles.push_back(bool_cycle);
-    }
-
-    // generate all possible cycles by taking xor of the edge sets of the known cycles
-    for (auto it1 = bool_cycles.begin(); it1 != std::prev(bool_cycles.end()); it1++) {
-      for (auto it2 = std::next(it1); it2 != bool_cycles.end(); it2++) {
-        std::valarray<uint8_t> new_bool_cycle = (*it1) ^ (*it2);
-        std::list<unsigned int> new_cycle;
-        unsigned int pos = 0;
-        for (uint8_t included : new_bool_cycle) {
-          if (included) {
-            new_cycle.push_back(pos);
-          }
-          pos++;
-        }
-        cycles.push_back(new_cycle);
-      }
-    }
-
-    // find the cycle representing the crack by counting boundary crossings
-    for (auto cycle : cycles) {
-      std::array<unsigned int, 2> crossing_count{0,0};
-
-      for (auto edge : cycle) {
-        double dx = fabs(n.G.dual_vertices[n.G.dual_edges[edge][0]].r.x - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.x);
-        if (dx > n.G.L.x / 2) {
-          crossing_count[0]++;
-        }
-        double dy = fabs(n.G.dual_vertices[n.G.dual_edges[edge][0]].r.y - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.y);
-        if (dy > n.G.L.y / 2) {
-          crossing_count[1]++;
-        }
-      }
-
-      if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
-        return cycle;
-      }
-    }
-  } else {
-    return cycles.front();
-  }
-
-  exit(5);
-}
-
 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);
@@ -153,8 +31,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 Lx, unsigned int Ly, double beta) {
-  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + 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, unsigned int Mx, unsigned int 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";
   std::ifstream file(filename);
 
   uint64_t N_old = 0;
@@ -226,18 +104,22 @@ void autocorrelation(unsigned int Lx, unsigned int Ly, std::vector<T>& out_data,
   }
 }
 
-ma::ma(unsigned int Lx, unsigned int Ly, double beta) : Lx(Lx), Ly(Ly), beta(beta), G(Lx * Ly / 2),
+unsigned int edge_r_to_ind(graph::coordinate r, unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My) {
+  return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);
+}
+
+ma::ma(unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My, double beta) : Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(Lx * Ly),
   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)
+  sb(log2(Mx < My ? Mx : My) + 1, 0),
+  Ccc((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Css((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Cmm((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Caa((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Cdd((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Cll((Mx / 2 + 1) * (My / 2 + 1), 0),
+  Cee((Mx / 2 + 1) * (My / 2 + 1), 0)
 {
   N = 0;
   Nc = 0;
@@ -246,17 +128,13 @@ ma::ma(unsigned int Lx, unsigned int Ly, double beta) : Lx(Lx), Ly(Ly), beta(bet
   // FFTW setup for correlation functions
   fftw_set_timelimit(1);
 
-  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(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);
-
-  std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_broken.dat";
+  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));
 
-  bondfile.open(filename, std::ifstream::app);
+  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() {
@@ -274,18 +152,16 @@ ma::~ma() {
   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);
-
-  bondfile.close();
+  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("Cdd", Cdd, 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("Cee", Cee, N, Lx, Ly, beta, Mx, My);
 }
 
-void ma::pre_fracture(const network &) {
+void ma::pre_fracture(const network&) {
   lv = std::numeric_limits<long double>::lowest();
   avalanches = {{}};
   boost::remove_edge_if(trivial, G);
@@ -297,13 +173,13 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned int i
     sa[avalanches.back().size() - 1]++;
     Na++;
 
-    std::fill_n(fftw_forward_in, net.G.edges.size(), 0.0);
+    std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
     for (auto e : avalanches.back()) {
-      fftw_forward_in[e] = 1.0;
+      fftw_forward_in[edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0;
     }
 
-    autocorrelation(Lx, Ly, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+    autocorrelation(Mx, My, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
     lv = c;
     avalanches.push_back({i});
@@ -311,121 +187,116 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned int i
     avalanches.back().push_back(i);
   }
 
-  boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], {i}, G);
-
-  bondfile << i << " " << c << std::scientific << " ";
+  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) {
-  bondfile << "\n";
   std::vector<unsigned int> component(boost::num_vertices(G));
   unsigned int num = boost::connected_components(G, &component[0]);
 
   std::list<unsigned int> crack = find_minimal_crack(G, n);
 
-  unsigned int crack_component = component[n.G.dual_edges[crack.front()][0]];
+  // crack surface correlations
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
+
+  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;
+  }
+
+  autocorrelation(Mx, My, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+
+  unsigned int crack_component = component[n.G.dual_edges[crack.front()].v[0]];
+
+  std::vector<std::list<unsigned int>> components(num);
+
+  for (unsigned int i = 0; i < n.G.dual_vertices.size(); i++) {
+    components[component[i]].push_back(i);
+  }
 
   // non-spanning clusters
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
   for (unsigned int i = 0; i < num; i++) {
     if (i != crack_component) {
-      unsigned int size = 0;
-
-      for (unsigned int j = 0; j < n.G.edges.size(); j++) {
-        if (component[n.G.dual_edges[j][0]] == i && n.fuses[j]) {
-          size++;
-          fftw_forward_in[j] = 1.0;
-        } else{
-          fftw_forward_in[j] = 0.0;
-        }
+      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)] += 1.0;
       }
 
-      if (size > 0) {
-        sc[size - 1]++;
-        autocorrelation(Lx, Ly, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-        Nc++;
-      }
-    }
-  }
+      sc[components[i].size() - 1]++;
+      autocorrelation(Mx, My, 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(Ly); be++) {
-    unsigned int bin = pow(2, be);
-
-    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 = 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;
-        }
-      }
-
-      if (in_bin) {
-        sb[be]++;
+      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;
       }
     }
   }
 
-  sb[log2(Ly)]++;
+  unsigned int max_factor = log2(Mx < My ? Mx : My) + 1;
+  std::vector<std::valarray<unsigned int>> bins(max_factor);
+  for (unsigned int i = 0; i < max_factor; i++) {
+    bins[i].resize(Mx * My / (unsigned int)pow(2, 2 * i));
+  }
 
-  for (unsigned int i = 0; i < n.G.edges.size(); i++) {
-    if (n.fuses[i] && component[n.G.dual_edges[i][0]] == crack_component)  {
-      fftw_forward_in[i] = 1.0;
-    } else {
-      fftw_forward_in[i] = 0.0;
+  for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
+    for (unsigned int i = 0; i < max_factor; i++) {
+      bins[i][edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx / pow(2, i), My / pow(2, i))] = 1;
     }
   }
 
-  autocorrelation(Lx, Ly, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  for (unsigned int i =0 ; i < max_factor; i++) {
+    sb[i] += bins[i].sum();
+  }
 
-  // crack surface correlations
-  std::fill_n(fftw_forward_in, n.G.edges.size(), 0.0);
+  // spanning cluster
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
-  for (auto edge : crack) {
-    fftw_forward_in[edge] = 1.0;
+  for (unsigned int i = 0; i < n.G.dual_vertices.size(); i++) {
+    if (component[i] == crack_component)  {
+      fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[i].r, Lx, Ly, Mx, My)] += 1.0;
+    }
   }
 
-  autocorrelation(Lx, Ly, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Mx, My, Cmm, 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];
+    return component[n.G.dual_edges[i].v[0]] == crack_component;
   };
 
   for (auto avalanche : avalanches) {
     if (avalanche.end() != std::find_if(avalanche.begin(), avalanche.end(), inCrack)) {
       for (auto edge : avalanche) {
-        fftw_forward_in[edge] = 1.0;
+        fftw_forward_in[edge_r_to_ind(n.G.edges[edge].r, Lx, Ly, Mx, My)] += 1.0;
       }
     }
   }
 
-  autocorrelation(Lx, Ly, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Mx, My, 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);
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
   // rewind the last avalanche
   for (auto e : avalanches.back()) {
-    boost::remove_edge(n.G.dual_edges[e][0], n.G.dual_edges[e][1], G);
+    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[e] = 1;
+    fftw_forward_in[edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0;
   }
 
-  autocorrelation(Lx, Ly, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Mx, My, 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;
 
+  std::fill_n(fftw_forward_in, Mx * My, 0.0);
+
   for (unsigned int i = 0; i < n.G.edges.size(); i++) {
     if (n.fuses[i]) { 
       total_broken++;
-      fftw_forward_in[i] = 1.0;
-    } else {
-      fftw_forward_in[i] = 0.0;
+      fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
     }
   }
 
-  autocorrelation(Lx, Ly, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+  autocorrelation(Mx, My, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
 
   sd[total_broken]++;