readded spatial measuremnets only for voronoi, and only for square systems

4 files changed, 61 insertions, 18 deletions

diff --git a/src/animate.cpp b/src/animate.cpp
index b1031dc..a7b6173 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -17,6 +17,7 @@ void animate::pre_fracture(const network &n) {
   lv = std::numeric_limits<long double>::lowest();
   avalanches = {};
   boost::remove_edge_if(trivial, G);
+  seen_guy = false;
 
   glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
   glLineWidth(5);
@@ -135,7 +136,7 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
       }
  glEnd();
  glFlush();
- if (nw > 2) {std::cout << "\n"; getchar();}
+ if (nw > 2 && !seen_guy) {seen_guy = true; getchar();}
 }
 
 void animate::post_fracture(network &n) {
diff --git a/src/animate.hpp b/src/animate.hpp
index fcb89fc..865b3bd 100644
--- a/src/animate.hpp
+++ b/src/animate.hpp
@@ -10,6 +10,7 @@
 class animate : public hooks {
   private:
     Graph G;
+    bool seen_guy;
   public:
     long double lv;
     std::list<std::list<unsigned>> avalanches;
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 96b4ccb..aca9cc6 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -38,15 +38,13 @@ 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, std::string model_string, unsigned Mx, unsigned My) {
-  std::string filename = model_string + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat";
+void update_field_file(std::string id, const std::vector<T>& data, std::string model_string) {
+  std::string filename = model_string + id + ".dat";
   std::ifstream file(filename);
 
-  uint64_t N_old = 0;
   std::vector<T> data_old(data.size(), 0);
 
   if (file.is_open()) {
-    file >> N_old;
     for (unsigned j = 0; j < data.size(); j++) {
       file >> data_old[j];
     }
@@ -55,7 +53,6 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned N, s
 
   std::ofstream file_out(filename);
 
-  file_out <<std::fixed<< N_old + N <<  "\n";
   for (unsigned j = 0; j < data.size(); j++) {
     file_out << data_old[j] + data[j] << " ";
   }
@@ -92,7 +89,7 @@ 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) {
+void autocorrelation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data, const std::list<graph::coordinate>& pos, std::array<unsigned, 2> count) {
   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);
@@ -101,7 +98,11 @@ 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[(unsigned)(Mx * (Δx / Lx)) + (Mx / 2) * (unsigned)(My * (Δy / Ly))]++;
+      if (count[0] % 2 == 0) {
+        data[(unsigned)(Mx * (Δx / Lx)) + (Mx / 2) * (unsigned)(My * (Δy / Ly))]++;
+      } else {
+        data[(unsigned)(Mx * (Δy / Ly)) + (Mx / 2) * (unsigned)(My * (Δx / Lx))]++;
+      }
     }
   }
 }
@@ -153,7 +154,14 @@ ma::ma(unsigned n, double a, double beta, double weight, bool one) :
   sa(3 * n),
   sA(3 * n),
   si(10000),
-  sI(10000)
+  sI(10000),
+  cn(pow((unsigned)sqrt(n), 2)),
+  cs(pow((unsigned)sqrt(n), 2)),
+  cm(pow((unsigned)sqrt(n), 2)),
+  cl(pow((unsigned)sqrt(n), 2)),
+  cb(pow((unsigned)sqrt(n), 2)),
+  ca(pow((unsigned)sqrt(n), 2)),
+  cA(pow((unsigned)sqrt(n), 2))
 {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_" + std::to_string(weight) + "_";
@@ -205,6 +213,13 @@ ma::~ma() {
   update_distribution_file("sA", sA, model_string);
   update_distribution_file("si", si, model_string);
   update_distribution_file("sI", sI, model_string);
+  update_field_file("cl", cl, model_string);
+  update_field_file("cm", cm, model_string);
+  update_field_file("cs", cs, model_string);
+  update_field_file("cn", cn, model_string);
+  update_field_file("cb", cb, model_string);
+  update_field_file("ca", ca, model_string);
+  update_field_file("cA", cA, model_string);
 }
 
 void ma::pre_fracture(const network&) {
@@ -240,18 +255,18 @@ void ma::post_fracture(network &n) {
   std::vector<unsigned> component(boost::num_vertices(G));
   unsigned num = boost::connected_components(G, &component[0]);
   if (post_cracks.size() > 2 || post_cracks.size() == 0) {
-    for (auto c : post_cracks) {
-      for (unsigned e : c.second) {
-        std::cout << e << " ";
-      }
-      std::cout << "\n";
-    }
-      getchar();
     throw badcycleex;
   }
+
   for (auto c : post_cracks) {
+    std::list<graph::coordinate> cl_cs;
     sl[c.second.size() - 1]++;
+    for (unsigned e : c.second) {
+      cl_cs.push_back(n.G.dual_edges[e].r);
+    }
+    autocorrelation2(n.G.L.x, n.G.L.y, 2 * sqrt(cl.size()), 2 * sqrt(cl.size()), cl, cl_cs, c.first);
   }
+
   unsigned crack_component = component[n.G.dual_edges[post_cracks.front().second.front()].v[0]];
 
   std::vector<std::list<graph::coordinate>> components(num);
@@ -263,10 +278,13 @@ void ma::post_fracture(network &n) {
   for (unsigned i = 0; i < num; i++) {
     if (i != crack_component) {
       sm[components[i].size() - 1]++;
+      autocorrelation2(n.G.L.x, n.G.L.y, 2 * sqrt(cm.size()), 2 * sqrt(cm.size()), cm, components[i], post_cracks.front().first);
     } else {
       ss[components[i].size() - 1]++;
+      autocorrelation2(n.G.L.x, n.G.L.y, 2 * sqrt(cs.size()), 2 * sqrt(cs.size()), cs, components[i], post_cracks.front().first);
     }
     sn[components[i].size() - 1]++;
+    autocorrelation2(n.G.L.x, n.G.L.y, 2 * sqrt(cn.size()), 2 * sqrt(cn.size()), cn, components[i], post_cracks.front().first);
   }
 
   std::vector<bool> vertex_in(n.G.vertices.size());
@@ -280,21 +298,37 @@ void ma::post_fracture(network &n) {
 
   unsigned bb_size = 0;
 
+  std::list<graph::coordinate> cb_co;
+
   for (unsigned i = 0; i < n.G.vertices.size(); i++) {
-    if (vertex_in[i]) bb_size++;
+    if (vertex_in[i]) {
+      bb_size++;
+      cb_co.push_back(n.G.vertices[i].r);
+    }
   }
 
   sb[bb_size]++;
+  autocorrelation2(n.G.L.x, n.G.L.y, 2 * sqrt(cb.size()), 2 * sqrt(cb.size()), cb, cb_co, post_cracks.front().first);
 
   auto av_it = avalanches.rbegin();
   av_it++;
 
   while (av_it != avalanches.rend()) {
     sa[(*av_it).size() - 1]++;
+    std::list<graph::coordinate> ca_co;
+    for (unsigned e : (*av_it)) {
+      ca_co.push_back(n.G.edges[e].r);
+    }
+    autocorrelation2(n.G.L.x, n.G.L.y, 2 * sqrt(ca.size()), 2 * sqrt(ca.size()), ca, ca_co, post_cracks.front().first);
     av_it++;
   }
 
   sA[avalanches.back().size() - 1]++;
+  std::list<graph::coordinate> cA_co;
+  for (unsigned e : avalanches.back()) {
+    cA_co.push_back(n.G.edges[e].r);
+  }
+  autocorrelation2(n.G.L.x, n.G.L.y, 2 * sqrt(cA.size()), 2 * sqrt(cA.size()), cA, cA_co, post_cracks.front().first);
 
   sd[num - 1]++;
 }
diff --git a/src/measurements.hpp b/src/measurements.hpp
index c211974..c69404b 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -30,13 +30,20 @@ class ma : public hooks {
     std::vector<uint64_t> sl; // final avalanche size distribution
     std::vector<uint64_t> sb; // final avalanche size distribution
     std::vector<uint64_t> sd; // final avalanche size distribution
-    std::vector<uint64_t> sc; // non-spanning cluster size distribution
     std::vector<uint64_t> sa; // non-final avalanche size distribution
     std::vector<uint64_t> sA; // non-final avalanche size distribution
 
     std::vector<uint64_t> si;
     std::vector<uint64_t> sI;
 
+    std::vector<uint64_t> cn;
+    std::vector<uint64_t> cs;
+    std::vector<uint64_t> cm;
+    std::vector<uint64_t> cl;
+    std::vector<uint64_t> cb;
+    std::vector<uint64_t> ca;
+    std::vector<uint64_t> cA;
+
   public:
     long double lv;