1 files changed, 167 insertions, 0 deletions
diff --git a/src/perc_meas.cpp b/src/perc_meas.cpp
new file mode 100644
index 0000000..4220a93
--- /dev/null
+++ b/src/perc_meas.cpp
@@ -0,0 +1,167 @@
+
+#include "perc_meas.hpp"
+#include <iostream>
+#include <cstdio>
+
+class badcycleException: public std::exception
+{
+  virtual const char* what() const throw()
+  {
+    return "Could not find a valid cycle on the broken system.";
+  }
+} badcycleex;
+
+void update_distribution_file(std::string id, const std::vector<uint64_t>& data, std::string model_string) {
+  std::string filename = model_string + id + ".dat";
+  std::ifstream file(filename);
+
+  std::vector<uint64_t> data_old(data.size(), 0);
+
+  if (file.is_open()) {
+    for (unsigned i = 0; i < data.size(); i++) {
+      uint64_t num;
+      file >> num;
+      data_old[i] = num;
+    }
+
+    file.close();
+  }
+
+  std::ofstream file_out(filename);
+
+  for (unsigned i = 0; i < data.size(); i++) {
+    file_out <<std::fixed<< data_old[i] + data[i] << " ";
+  }
+
+  file_out.close();
+}
+
+void update_distribution_file(std::string id, const std::vector<std::vector<uint64_t>>& data, std::string model_string) {
+  std::string filename = model_string + id + ".dat";
+  std::ifstream file(filename);
+
+  std::vector<std::vector<uint64_t>> data_old(data.size());
+
+  for (unsigned i = 0; i < data.size(); i++) {
+    data_old[i].resize(data[0].size());
+  }
+
+  if (file.is_open()) {
+    for (unsigned i = 0; i < data.size(); i++) {
+      for (unsigned j = 0; j < data[0].size(); j++) {
+        uint64_t num;
+        file >> num;
+        data_old[i][j] = num;
+      }
+    }
+
+    file.close();
+  }
+
+  std::ofstream file_out(filename);
+
+  for (unsigned i = 0; i < data.size(); i++) {
+    for (unsigned j = 0; j < data[0].size(); j++) {
+      file_out << std::fixed << data_old[i][j] + data[i][j] << " ";
+    }
+    file_out << "\n";
+  }
+
+  file_out.close();
+}
+
+pm::pm(unsigned n, double a) :
+  G(2 * n),
+  rank(2 * n),
+  parent(2 * n),
+  ds(&rank[0], &parent[0]),
+  sm(2 * n),
+  sl(2 * n),
+  sn(3 * n),
+  ss(2 * n),
+  sb(3 * n),
+  sp(3 * n)
+{
+  model_string = "percolation_" + std::to_string(n) + "_" + std::to_string(a) + "_";
+  for (std::vector<uint64_t> &x : sn) {
+    x.resize(2 * n);
+  }
+  for (std::vector<uint64_t> &x : sb) {
+    x.resize(3 * n);
+  }
+}
+
+pm::~pm() {
+  update_distribution_file("sm", sm, model_string);
+  update_distribution_file("sl", sl, model_string);
+  update_distribution_file("sn", sn, model_string);
+  update_distribution_file("ss", ss, model_string);
+  update_distribution_file("sb", sb, model_string);
+  update_distribution_file("sp", sp, model_string);
+}
+
+void pm::pre_fracture(const network&) {
+  boost::remove_edge_if(trivial, G);
+  initialize_incremental_components(G, ds);
+  incremental_components(G, ds);
+  p = 0;
+}
+
+void pm::bond_broken(const network& net, const current_info& cur, unsigned i) {
+  boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
+  ds.union_set(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1]);
+
+  boost::component_index<VertexIndex> components(parent.begin(), parent.end());
+  std::vector<unsigned> counts(components.size());
+
+  for (unsigned i = 0; i < components.size(); i++) {
+    unsigned comp_size = 0;
+    BOOST_FOREACH(VertexIndex child_index, components[i]) {
+      comp_size++;
+    }
+
+    sn[p][comp_size - 1]++;
+  }
+
+  unsigned bb_size = 0;
+
+  for (unsigned i = 0; i <net.G.edges.size(); i++) {
+    if (cur.currents[i] > 1.0 / net.G.edges.size()) {
+      bb_size++;
+    }
+  }
+
+  sb[p][bb_size - 1]++;
+  p++;
+}
+
+void pm::post_fracture(network &n) {
+  auto post_cracks = find_minimal_crack(G, 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) {
+    throw badcycleex;
+  }
+  for (auto c : post_cracks) {
+    sl[c.second.size() - 1]++;
+  }
+  unsigned crack_component = component[n.G.dual_edges[post_cracks.front().second.front()].v[0]];
+
+  std::vector<std::list<graph::coordinate>> components(num);
+
+  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
+    components[component[i]].push_back(n.G.dual_vertices[i].r);
+  }
+
+  for (unsigned i = 0; i < num; i++) {
+    if (i != crack_component) {
+      sm[components[i].size() - 1]++;
+    } else {
+      ss[components[i].size() - 1]++;
+    }
+    sn[p][components[i].size() - 1]++;
+  }
+
+  sp[p - 1]++;
+}
+