made backbone cleaning more efficient by restricting to one-side-broken only, updated the crack finder to this new paradigm

5 files changed, 73 insertions, 239 deletions

diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp
index 5f180ee..2965915 100644
--- a/src/analysis_tools.cpp
+++ b/src/analysis_tools.cpp
@@ -1,125 +1,72 @@
 
 #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<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> find_minimal_crack(const Graph& G, const network& n) {
-  Graph Gtmp(n.G.vertices.size());
-  std::list<unsigned> removed_edges;
-
-  class add_tree_edges : public boost::default_dfs_visitor {
-    public:
-      Graph& G;
-      std::list<unsigned>& E;
-
-      add_tree_edges(Graph& G, std::list<unsigned>& 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>& E;
-      unsigned end;
-      struct done{};
-
-      find_cycle(std::list<unsigned>& E, unsigned 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>> cycles;
+std::pair<std::array<unsigned, 2>, std::set<unsigned>> find_minimal_crack(const network& n, unsigned i0) {
+  std::set<unsigned> removed_edges = n.get_cycle_edges(n.G.dual_edges[i0].v[0]);
+
+  std::vector<std::pair<std::array<unsigned, 2>, std::set<unsigned>>> cycles;
+  cycles.reserve(removed_edges.size());
+
+  for (unsigned edge : removed_edges) {
+    cycles.push_back(n.get_cycle(removed_edges, n.G.dual_edges[edge].v[0], n.G.dual_edges[edge].v[1])); 
+    cycles.back().second.insert(edge);
+    if (n.G.dual_edges[edge].crossings[0])
+      cycles.back().first[0]++;
+    if (n.G.dual_edges[edge].crossings[1])
+      cycles.back().first[1]++;
+  }
 
-  for (auto edge : removed_edges) {
-    std::list<unsigned> 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) {
+    return cycles.front();
   }
 
-  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 all_good = true;
+  unsigned not_good;
+  for (unsigned i = 0; i < 2; i++) {
+    if (!((cycles[i].first[0] % 2 == 0 && cycles[i].first[1] % 2 == 1) || (cycles[i].first[0] % 2 == 1 && cycles[i].first[1] % 2 == 0))) {
+      all_good = false;
+      not_good = i;
     }
-    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> new_cycle;
-      unsigned pos = 0;
-      for (uint8_t included : new_bool_cycle) {
-        if (included) {
-          new_cycle.push_back(pos);
-        }
-        pos++;
-      }
-      cycles.push_back(new_cycle);
+  if (all_good) {
+    if (cycles[0].second.size() > cycles[1].second.size()) {
+      return cycles[1];
+    } else {
+      return cycles[0];
     }
   }
 
-  std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> output;
+  const auto& cycle1 = cycles[0];
+  const auto& cycle2 = cycles[1];
 
-  // find the cycle representing the crack by counting boundary crossings
-  for (auto cycle : cycles) {
-    std::array<unsigned, 2> crossing_count{0,0};
+  unsigned sum_sig_0 = cycle1.first[0] + cycle2.first[0];
+  unsigned sum_sig_1 = cycle1.first[1] + cycle2.first[1];
 
-    for (auto edge : cycle) {
-      if (n.G.dual_edges[edge].crossings[0]) {
-        crossing_count[0]++;
-      }
-      if (n.G.dual_edges[edge].crossings[1]) {
-        crossing_count[1]++;
-      }
-    }
+  std::valarray<uint8_t> bool_cycle1(n.G.edges.size());
+  std::valarray<uint8_t> bool_cycle2(n.G.edges.size());
+  for (auto v : cycle1.second) {
+    bool_cycle1[v] = 1;
+  }
+  for (auto v : cycle2.second) {
+    bool_cycle2[v] = 1;
+  }
 
-    if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
-      output.push_back({{1, 0}, cycle});
-    } else if (crossing_count[0] % 2 == 0 && crossing_count[1] % 2 == 1) {
-      output.push_back({{0, 1}, cycle});
+  std::valarray<uint8_t> new_bool_cycle = bool_cycle1 ^ bool_cycle2;
+
+  std::set<unsigned> new_cycle;
+  unsigned pos = 0;
+  for (uint8_t included : new_bool_cycle) {
+    if (included) {
+      new_cycle.insert(pos);
     }
+    pos++;
   }
 
-  return output;
+  if (cycles[!not_good].second.size() > new_cycle.size()) {
+    return {{sum_sig_0, sum_sig_1}, new_cycle};
+  } else {
+    return cycles[!not_good];
+  }
 }
-*/
 
diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp
index 2a4c5d2..72d1e82 100644
--- a/src/analysis_tools.hpp
+++ b/src/analysis_tools.hpp
@@ -7,11 +7,5 @@
 
 #include <network.hpp>
 
-
-template<class T>
-bool is_shorter(const std::list<T> &, const std::list<T> &);
-
-//bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>);
-
-//std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> find_minimal_crack(const Graph &, const network &);
+std::pair<std::array<unsigned, 2>, std::set<unsigned>> find_minimal_crack(const network &, unsigned i0);
 
diff --git a/src/animate.cpp b/src/animate.cpp
index 6c844d2..81bf4ac 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -1,5 +1,6 @@
 
 #include "animate.hpp"
+#include <unistd.h>
 #include <iostream>
 
 animate::animate(double Lx, double Ly, unsigned window_size, int argc, char *argv[]) {
@@ -139,65 +140,27 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
       }
  glEnd();
  glFlush();
- if (nw > 2 && !seen_guy) {seen_guy = true; getchar();}
+// if (nw > 2 && !seen_guy) {seen_guy = true; getchar();}
 }
 
 void animate::post_fracture(network &n) {
 
-  /*
-  std::list<unsigned> crack;
-//  unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
-  unsigned crack_component = 10000;
-
-  std::default_random_engine gen;
-  std::uniform_real_distribution<double> dis(0.0,1.0);
-
-  auto av_it = avalanches.rbegin();
-
-  while (true) {
-    for (unsigned e : *av_it) {
-      boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
-      n.fuses[e] = false;
-    }
-
-    auto cracks = find_minimal_crack(G, n);
-
-    if (cracks.size() == 0) {
-      break;
-    }
+  auto crack = find_minimal_crack(n, avalanches.back().back());
 
-    av_it++;
-  }
-
-  std::vector<unsigned> component(boost::num_vertices(G));
-  unsigned num = boost::connected_components(G, &component[0]);
-
-  std::vector<std::list<unsigned>> components(num);
-
-  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
-    components[component[i]].push_back(i);
-  }
-
-
-  char key;
-  while ((key = getchar()) != 'n') {
-    glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
-    glClear(GL_COLOR_BUFFER_BIT);
-    glBegin(GL_LINES);
-    glColor3f(0.0f, 0.0f, 0.0f);
-    for (unsigned i = 0; i < n.G.edges.size(); i++) {
-      if (!n.fuses[i]) {
-        graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
-        graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
+  for (unsigned j : crack.second) {
+  glBegin(GL_LINES);
+  glColor3d(1.0, 0.0, 0.0);
+    graph::coordinate r1 = n.G.dual_vertices[n.G.dual_edges[j].v[0]].r;
+    graph::coordinate r2 = n.G.dual_vertices[n.G.dual_edges[j].v[1]].r;
 
-    if (n.G.edges[i].crossings[0]) {
+    if (n.G.dual_edges[j].crossings[0]) {
       if (r1.x < r2.x) {
         r1.x += n.G.L.x;
       } else {
         r2.x += n.G.L.x;
       }
     }
-    if (n.G.edges[i].crossings[1]) {
+    if (n.G.dual_edges[j].crossings[1]) {
       if (r1.y < r2.y) {
         r1.y += n.G.L.y;
       } else {
@@ -205,82 +168,12 @@ void animate::post_fracture(network &n) {
       }
     }
 
-        glVertex2d(r1.x, r1.y);
-        glVertex2d(r2.x, r2.y);
-      }
-    }
-    glEnd();
-
-    switch (key) {
-      case 's' :
-        for (auto edge : crack) {
-          glColor3f(1.0f, 0.0f, 0.0f);
-         glBegin(GL_POLYGON);
-          for (const graph::coordinate &r : n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon) {
-            glVertex2d(r.x, r.y);
-          }
-          glEnd();
-         glBegin(GL_POLYGON);
-          for (const graph::coordinate &r : n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon) {
-            glVertex2d(r.x, r.y);
-          }
-          glEnd();
-        }
-        glFlush();
-        break;
-
-      case 'c' :
-        for (unsigned i = 0; i < num; i++) {
-          if (i == crack_component) {
-            glColor3d(1.0, 0.0, 0.0);
-          } else {
-            glColor3d(dis(gen), dis(gen), dis(gen));
-          }
-
-          for (auto it = components[i].begin(); it != components[i].end(); it++) {
-        glBegin(GL_POLYGON);          // Each set of 3 vertices form a triangle
-            for (const graph::coordinate &r: n.G.dual_vertices[*it].polygon) {
-              glVertex2d(r.x, r.y);
-            }
-        glEnd();
-          }
-        }
-        glFlush();
-        break;
-
-      case 'C' :
-        for (unsigned i = 0; i < num; i++) {
-          if (components[i].size() > 1) {
-            if (i == crack_component) {
-              glColor3d(1.0, 0.0, 0.0);
-            } else {
-              glColor3d(dis(gen), dis(gen), dis(gen));
-            }
-
-            for (auto it = components[i].begin(); it != components[i].end(); it++) {
-        glBegin(GL_POLYGON);          // Each set of 3 vertices form a triangle
-              for (const graph::coordinate &r :n.G.dual_vertices[*it].polygon) {
-                glVertex2d(r.x, r.y);
-              }
-        glEnd();
-            }
-          }
-        }
-        glFlush();
-        break;
-
-      case 'm' :
-        for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
-        glBegin(GL_POLYGON);
-        glColor3d(1.0, 0.0, 0.0);
-          for (const graph::coordinate &r :n.G.dual_vertices[*it].polygon) {
-            glVertex2d(r.x, r.y);
-          }
-        glEnd();
-        }
-        glFlush();
-    }
+    glVertex2d(r1.x, r1.y);
+    glVertex2d(r2.x, r2.y);
+  glEnd();
   }
-*/
+  glFlush();
+
+  getchar();
 }
 
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
index 9d875e6..2b5e55c 100644
--- a/src/animate_fracture.cpp
+++ b/src/animate_fracture.cpp
@@ -57,7 +57,7 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  animate meas(sqrt(2*n *a), sqrt( 2*n / a), 700, argc, argv);
+  animate meas(sqrt(2*n *a), sqrt( 2*n / a), 1000, argc, argv);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
@@ -66,7 +66,7 @@ int main(int argc, char* argv[]) {
     graph G(n, a, rng);
     elastic_network elastic_network(G, &c);
     elastic_network.set_thresholds(beta, rng);
-    elastic_network.fracture(meas, true);
+    elastic_network.fracture(meas);
 
     if (quit.load())
       break;
diff --git a/src/animate_fracture_square.cpp b/src/animate_fracture_square.cpp
index 12de7da..ad43254 100644
--- a/src/animate_fracture_square.cpp
+++ b/src/animate_fracture_square.cpp
@@ -57,7 +57,7 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  animate meas(Lx, Ly, 512, argc, argv);
+  animate meas(Lx, Ly, 1000, argc, argv);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
@@ -68,7 +68,7 @@ int main(int argc, char* argv[]) {
   for (unsigned trial = 0; trial < N; trial++) {
     elastic_network tmp_network(perm_network);
     tmp_network.set_thresholds(beta, rng);
-    tmp_network.fracture(meas, true);
+    tmp_network.fracture(meas);
 
     if (quit.load())
       break;