1 files changed, 125 insertions, 0 deletions
diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp
new file mode 100644
index 0000000..1e799bb
--- /dev/null
+++ b/src/analysis_tools.cpp
@@ -0,0 +1,125 @@
+
+#include "analysis_tools.hpp"
+
+template <class T>
+bool is_shorter(const std::list<T> &l1, const 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].v[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].v[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].v[0]].r.x - n.G.dual_vertices[n.G.dual_edges[edge].v[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].v[0]].r.y - n.G.dual_vertices[n.G.dual_edges[edge].v[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);
+}
+