From d9d3b0518ce5e0a52b9a0bae55fa5d8ca5b3c515 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 24 Sep 2019 17:53:08 -0400 Subject: made backbone cleaning more efficient by restricting to one-side-broken only, updated the crack finder to this new paradigm --- lib/include/network.hpp | 8 +- lib/src/network.cpp | 78 ++++++++++++++++---- src/analysis_tools.cpp | 157 +++++++++++++--------------------------- src/analysis_tools.hpp | 8 +- src/animate.cpp | 139 ++++------------------------------- src/animate_fracture.cpp | 4 +- src/animate_fracture_square.cpp | 4 +- 7 files changed, 144 insertions(+), 254 deletions(-) diff --git a/lib/include/network.hpp b/lib/include/network.hpp index 29bf55a..5b474bb 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -35,9 +35,10 @@ private: void update_backbone(const std::vector& c); void break_edge(unsigned, bool unbreak = false); void get_cycle_edges_helper(std::set&, std::set&, unsigned, unsigned) const; - std::set get_cycle_edges(unsigned) const; bool find_cycle_helper(std::array&, const std::set&, unsigned, unsigned, unsigned) const; + bool get_cycle_helper(std::array&, std::set&, const std::set&, unsigned, unsigned, + unsigned) const; std::array find_cycle(const std::set&, unsigned, unsigned) const; void get_cluster_edges_helper(std::set&, unsigned) const; std::set get_cluster_edges(unsigned) const; @@ -56,7 +57,10 @@ public: network(const graph&, cholmod_common*); void set_thresholds(double, std::mt19937&); - void fracture(hooks&, bool one_axis = true); + void fracture(hooks&); + + std::set get_cycle_edges(unsigned) const; + std::pair, std::set> get_cycle(const std::set&, unsigned, unsigned) const; virtual current_info get_current_info() { current_info empty; diff --git a/lib/src/network.cpp b/lib/src/network.cpp index eef7e9a..3e460a9 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -6,8 +6,8 @@ class nofuseException : public std::exception { } nofuseex; network::network(const graph& G, cholmod_common* c) - : G(G), C(G.dual_vertices.size()), fuses(G.edges.size()), backbone(G.edges.size()), - thresholds(G.edges.size()), px(G, 0, c), py(G, 1, c) {} + : px(G, 0, c), py(G, 1, c), G(G), C(G.dual_vertices.size()), fuses(G.edges.size()), + backbone(G.edges.size()), thresholds(G.edges.size()) {} void network::set_thresholds(double beta, std::mt19937& rng) { if (beta == 0.0) { @@ -29,25 +29,22 @@ void network::set_thresholds(double beta, std::mt19937& rng) { } } -void network::fracture(hooks& m, bool one_axis) { +void network::fracture(hooks& m) { m.pre_fracture(*this); while (true) { - double min_cond = 1.0 / G.edges.size(); + const double min_cond = 1.0 / G.edges.size(); current_info c = this->get_current_info(); - if (px.sol.conductivity[0] > min_cond) { - this->update_backbone(px.sol.currents); - } else { - this->update_backbone(py.sol.currents); - } - if (c.conductivity[0] < min_cond && c.conductivity[1] < min_cond) { + if (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) { break; } - if (one_axis && (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond)) { - break; + if (px.sol.conductivity[0] > min_cond) { + this->update_backbone(px.sol.currents); + } else { + this->update_backbone(py.sol.currents); } unsigned max_pos = UINT_MAX; @@ -144,6 +141,48 @@ std::array network::find_cycle(const std::set& cycle_edge return sig; } +bool network::get_cycle_helper(std::array& sig, std::set& edges, const std::set& cycle_edges, + unsigned vPrev, unsigned vCur, unsigned vEnd) const { + for (unsigned ei : G.dual_vertices[vCur].ne) { + if (fuses[ei]) { + auto it = cycle_edges.find(ei); + if (it == cycle_edges.end()) { + const std::array& e = G.dual_edges[ei].v; + unsigned vn = e[0] == vCur ? e[1] : e[0]; + if (vn != vPrev) { + if (vn == vEnd) { + edges.insert(ei); + if (G.dual_edges[ei].crossings[0]) + sig[0]++; + if (G.dual_edges[ei].crossings[1]) + sig[1]++; + return true; + } else { + if (this->get_cycle_helper(sig, edges, cycle_edges, vCur, vn, vEnd)) { + edges.insert(ei); + if (G.dual_edges[ei].crossings[0]) + sig[0]++; + if (G.dual_edges[ei].crossings[1]) + sig[1]++; + return true; + } + } + } + } + } + } + + return false; +} + +std::pair, std::set> network::get_cycle(const std::set& cycle_edges, unsigned v0, + unsigned v1) const { + std::set edges; + std::array sig = {0, 0}; + this->get_cycle_helper(sig, edges, cycle_edges, v0, v0, v1); + return {sig, edges}; +} + void network::get_cluster_edges_helper(std::set& seen_edges, unsigned v) const { for (unsigned ei : G.vertices[v].ne) { if (!backbone[ei]) { @@ -264,7 +303,7 @@ void network::update_backbone(const std::vector& c) { return (c[0] % 2 == 0 && c[1] % 2 == 0) || ((c[0] % 2 == 0 && done_x) || (c[1] % 2 == 0 && done_y)); })) { - // If it does, we remove it from the backbone! + // If none does, we remove it from the backbone! seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; backbone[i] = true; @@ -301,6 +340,14 @@ void network::update_backbone(const std::vector& c) { } } } + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), + [](std::array c) { + return ((c[0] % 2 == 0 && c[1] % 2 == 1) || (c[0] % 2 == 1 && c[1] % 2 == 0)); + })) { + // If the bond separates a wrapping path, breaking it would end the + // fracture and therefore it will never be removed from the backbone + // in this manner. We can skip it in the future. + seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; } } } @@ -450,6 +497,11 @@ void network::update_backbone(const std::vector& c) { } } } + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), + [](std::array c) { + return ((c[0] % 2 == 0 && c[1] % 2 == 1) || (c[0] % 2 == 1 && c[1] % 2 == 0)); + })) { + seen_pairs[{d1, d2}] = true; } } } 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 -bool is_shorter(const std::list &l1, const std::list &l2) { - return l1.size() < l2.size(); -} - -/* -bool trivial(boost::detail::edge_desc_impl) { - return true; -} - -std::list, std::list>> find_minimal_crack(const Graph& G, const network& n) { - Graph Gtmp(n.G.vertices.size()); - std::list removed_edges; - - class add_tree_edges : public boost::default_dfs_visitor { - public: - Graph& G; - std::list& E; - - add_tree_edges(Graph& G, std::list& E) : G(G), E(E) {} - - void tree_edge(boost::graph_traits::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::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& E; - unsigned end; - struct done{}; - - find_cycle(std::list& E, unsigned end) : E(E), end(end) {} - - void discover_vertex(boost::graph_traits::vertex_descriptor v, const Graph& g) { - if (v == end) { - throw done{}; - } - } - - void examine_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { - E.push_back(g[e].index); - } - - void finish_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { - E.erase(std::find(E.begin(), E.end(), g[e].index)); - } - }; - - std::list> cycles; +std::pair, std::set> find_minimal_crack(const network& n, unsigned i0) { + std::set removed_edges = n.get_cycle_edges(n.G.dual_edges[i0].v[0]); + + std::vector, std::set>> 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 cycle = {edge}; - find_cycle vis(cycle, n.G.dual_edges[edge].v[1]); - std::vector 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> bool_cycles; - for (auto cycle : cycles) { - std::valarray 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 new_bool_cycle = (*it1) ^ (*it2); - std::list 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::list>> 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 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 bool_cycle1(n.G.edges.size()); + std::valarray 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 new_bool_cycle = bool_cycle1 ^ bool_cycle2; + + std::set 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 - -template -bool is_shorter(const std::list &, const std::list &); - -//bool trivial(boost::detail::edge_desc_impl); - -//std::list, std::list>> find_minimal_crack(const Graph &, const network &); +std::pair, std::set> 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 #include 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 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 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 component(boost::num_vertices(G)); - unsigned num = boost::connected_components(G, &component[0]); - - std::vector> 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; -- cgit v1.2.3-54-g00ecf