From 3f7f20f21f583ca2de566bea08a87eac4b17ad29 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 23 Sep 2019 23:19:13 -0400 Subject: successfully implemented the backbone trimming with homespun graph algorithms, measurements have mostly been disabled and need to be migrated --- lib/include/network.hpp | 56 +++-- lib/src/graph.cpp | 22 +- lib/src/network.cpp | 472 ++++++++++++++++------------------------ src/analysis_tools.cpp | 2 + src/analysis_tools.hpp | 4 +- src/animate.cpp | 9 +- src/animate.hpp | 1 - src/animate_fracture_square.cpp | 2 +- src/measurements.cpp | 13 +- src/measurements.hpp | 2 - 10 files changed, 258 insertions(+), 325 deletions(-) diff --git a/lib/include/network.hpp b/lib/include/network.hpp index bf9015c..8748ea9 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -18,34 +18,46 @@ #define CURRENT_CUTOFF 1e-10 +class fuse_network; +class elastic_network; +class percolation_network; + class network { - private: - problem px; - problem py; - std::unordered_map, bool> seen_pairs; + friend class fuse_network; + friend class elastic_network; + friend class percolation_network; - 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; - std::array find_cycle(const std::set&, unsigned, unsigned) const; + private: + problem px; + problem py; + std::unordered_map, bool> seen_pairs; + + 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; + 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; + void get_tie_flaps_helper(std::set&, unsigned, unsigned) const; + std::list> get_tie_flaps(unsigned) const; public: - const graph& G; + const graph& G; - ClusterTree C; + ClusterTree C; - std::vector fuses; - std::vector backbone; - std::vector thresholds; + std::vector fuses; + std::vector backbone; + std::vector thresholds; - network(const graph&, cholmod_common*); - network(const network&); + network(const graph&, cholmod_common*); - void set_thresholds(double, std::mt19937&); - void fracture(hooks&, bool one_axis = true); + void set_thresholds(double, std::mt19937&); + void fracture(hooks&, bool one_axis = true); - virtual current_info get_current_info() const {current_info empty; return empty;}; + virtual current_info get_current_info() {current_info empty; return empty;}; }; class fuse_network : public network { @@ -56,7 +68,7 @@ class fuse_network : public network { fuse_network(const graph&, cholmod_common*, double weight = 1.0); fuse_network(const fuse_network&); - current_info get_current_info() const override; + current_info get_current_info() override; }; class elastic_network : public network { @@ -67,7 +79,7 @@ class elastic_network : public network { elastic_network(const graph&, cholmod_common*, double weight = 0.5); elastic_network(const elastic_network&); - current_info get_current_info() const override; + current_info get_current_info() override; }; class percolation_network : public network { @@ -75,6 +87,6 @@ class percolation_network : public network { percolation_network(const graph&, cholmod_common*); percolation_network(const percolation_network&); - current_info get_current_info() const override; + current_info get_current_info() override; }; diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp index bfdb952..a5063de 100644 --- a/lib/src/graph.cpp +++ b/lib/src/graph.cpp @@ -44,7 +44,7 @@ graph::graph(unsigned Nx, unsigned Ny) { dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))); dual_vertices[i].r.y = (double)(i / (Nx / 2)); - dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx * (i / (Nx / 2)) + (i % (Nx / 2)), (nv + i - Nx / 2) % nv}; + dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv}; dual_vertices[i].polygon = { {dual_vertices[i].r.x - 1.0, vertices[i].r.y}, {dual_vertices[i].r.x, vertices[i].r.y - 1.0}, @@ -93,6 +93,16 @@ graph::graph(unsigned Nx, unsigned Ny) { } } + for (vertex& v : dual_vertices) { + v.ne.reserve(v.nd.size()); + } + + for (unsigned i = 0; i < dual_edges.size(); i++) { + for (unsigned vi : dual_edges[i].v) { + dual_vertices[vi].ne.push_back(i); + } + } + } class eulerException: public std::exception @@ -411,6 +421,16 @@ void graph::helper(unsigned nv, std::mt19937& rng) { } } + for (vertex& v : dual_vertices) { + v.ne.reserve(v.nd.size()); + } + + for (unsigned i = 0; i < dual_edges.size(); i++) { + for (unsigned vi : dual_edges[i].v) { + dual_vertices[vi].ne.push_back(i); + } + } + jcv_diagram_free(&diagram); } diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 2a187c5..0cf50c7 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -77,6 +77,7 @@ void network::fracture(hooks& m, bool one_axis) { void network::get_cycle_edges_helper(std::set& cycle_edges, std::set& seen_vertices, unsigned v_prev, unsigned v_cur) const { + seen_vertices.insert(v_cur); for (unsigned ei : G.dual_vertices[v_cur].ne) { if (fuses[ei]) { const std::array& e = G.dual_edges[ei].v; @@ -84,9 +85,8 @@ void network::get_cycle_edges_helper(std::set& cycle_edges, if (vn != v_prev) { auto it = seen_vertices.find(vn); - if (it == seen_vertices.end()) { - // if (seen_vertices.contains()) { // need to wait for better c++20 support... - cycle_edges.insert(vn); + if (it != seen_vertices.end()) { + cycle_edges.insert(ei); } else { this->get_cycle_edges_helper(cycle_edges, seen_vertices, v_cur, vn); } @@ -107,16 +107,21 @@ std::set network::get_cycle_edges(unsigned v0) const { bool network::find_cycle_helper(std::array& sig, const std::set& cycle_edges, unsigned vPrev, unsigned vCur, unsigned vEnd) const { for (unsigned ei : G.dual_vertices[vCur].ne) { if (fuses[ei]) { - const std::array& e = G.dual_edges[ei].v; - unsigned vn = e[0] == vCur ? e[1] : e[0]; - if (vn != vPrev) { - if (vn == vEnd) { - return true; - } else { - if (find_cycle_helper(sig, cycle_edges, vCur, vn, vEnd)) { + 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) { if (G.dual_edges[ei].crossings[0]) sig[0]++; if (G.dual_edges[ei].crossings[1]) sig[1]++; return true; + } else { + if (this->find_cycle_helper(sig, cycle_edges, vCur, vn, vEnd)) { + if (G.dual_edges[ei].crossings[0]) sig[0]++; + if (G.dual_edges[ei].crossings[1]) sig[1]++; + return true; + } } } } @@ -127,202 +132,169 @@ bool network::find_cycle_helper(std::array& sig, const std::set network::find_cycle(const std::set& cycle_edges, unsigned v0, unsigned v1) const { - std::array sig; - this->find_cycle_helper(sig, cycle_edges, v0, v0, + std::array sig = {0, 0}; + this->find_cycle_helper(sig, cycle_edges, v0, v0, v1); + return sig; } -void network::update_backbone(const std::vector& c) { - bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size(); - bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size(); - - class find_cycle : public boost::default_dfs_visitor { - public: - const graph& G; - std::array& E; - unsigned end; - struct done {}; - - find_cycle(const graph& G, std::array& E, unsigned end) : G(G), E(E), end(end) {} +void network::get_cluster_edges_helper(std::set& seen_edges, unsigned v) const { + for (unsigned ei : G.vertices[v].ne) { + if (!backbone[ei]) { + auto it = seen_edges.find(ei); + if (it == seen_edges.end()) { + const std::array& e = G.edges[ei].v; + unsigned vn = e[0] == v ? e[1] : e[0]; - void discover_vertex(boost::graph_traits::vertex_descriptor v, const Graph& g) { - if (v == end) { - throw done{}; + seen_edges.insert(ei); + this->get_cluster_edges_helper(seen_edges, vn); } } + } +} - void examine_edge(boost::graph_traits::edge_descriptor ed, const Graph& g) { - unsigned e = g[ed].index; - if (G.dual_edges[e].crossings[0]) - E[0]++; - if (G.dual_edges[e].crossings[1]) - E[1]++; - } +std::set network::get_cluster_edges(unsigned v0) const { + std::set cluster_edges; + this->get_cluster_edges_helper(cluster_edges, v0); + return cluster_edges; +} - void finish_edge(boost::graph_traits::edge_descriptor ed, const Graph& g) { - unsigned e = g[ed].index; - if (G.dual_edges[e].crossings[0]) - E[0]--; - if (G.dual_edges[e].crossings[1]) - E[1]--; +void network::get_tie_flaps_helper(std::set& added_edges, unsigned v0, unsigned vCur) const { + for (unsigned ei : G.vertices[vCur].ne) { + if (!backbone[ei]) { + auto it = added_edges.find(ei); + if (it == added_edges.end()) { + const std::array& e = G.edges[ei].v; + unsigned vn = e[0] == vCur ? e[1] : e[0]; + + if (vn != v0) { + added_edges.insert(ei); + this->get_tie_flaps_helper(added_edges, v0, vn); + } + } } - }; + } +} - class get_cluster_edges : public boost::default_dfs_visitor { - public: - std::set& E; +std::list> network::get_tie_flaps(unsigned v0) const { + std::list> tie_flaps; + + for (unsigned ei : G.vertices[v0].ne) { + if (!backbone[ei]) { + bool seen_edge = false; + for (const std::set& flap : tie_flaps) { + auto it = flap.find(ei); + if (it != flap.end()) { + seen_edge = true; + break; + } + } - get_cluster_edges(std::set& E) : E(E) {} + if (!seen_edge) { + tie_flaps.push_back({ei}); - void examine_edge(boost::graph_traits::edge_descriptor ed, const Graph& g) { - unsigned e = g[ed].index; - E.insert(e); + const std::array& e = G.edges[ei].v; + unsigned vn = e[0] == v0 ? e[1] : e[0]; + + this->get_tie_flaps_helper(tie_flaps.back(), v0, vn); + } } - }; + } + + return tie_flaps; +} + +void network::update_backbone(const std::vector& c) { + bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size(); + bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size(); + // First, we will check for lollipops! for (unsigned i = 0; i < G.edges.size(); i++) { if ((!backbone[i]) && C.same_component(G.dual_edges[i].v[0], G.dual_edges[i].v[1])) { auto it_search = seen_pairs.find({G.dual_edges[i].v[0], G.dual_edges[i].v[1]}); if (it_search == seen_pairs.end()) { - std::list cedges = this->get_cycle_edges(G.dual_edges[i].v[0], G.dual_edges[i].v[1]); - - std::array cycle = {0, 0}; + // This is a candidate lollipop stem. First, we will identify any + // cycles in the dual cluster that impinges on the stem and mark them + // by an edge that uniquely severs each. + std::set cedges = this->get_cycle_edges(G.dual_edges[i].v[0]); + + // Now, we create a crossing signature for each cycle. First, we do it + // for the new cycle that would form by removing the stem. + std::array base_sig = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[i].v[1]); if (G.dual_edges[i].crossings[0]) - cycle[0]++; + base_sig[0]++; if (G.dual_edges[i].crossings[1]) - cycle[1]++; - find_cycle vis2(G, cycle, G.dual_edges[i].v[1]); - std::vector cm2(G.dual_vertices.size()); - try { - boost::depth_first_visit(dG, G.dual_edges[i].v[0], vis2, - boost::make_iterator_property_map( - cm2.begin(), boost::get(boost::vertex_index, dG), cm2[0])); - } catch (find_cycle::done const&) { - std::list> other_cycles = {cycle}; - for (unsigned e : cedges) { - std::array tcycle = {0, 0}; - if (G.dual_edges[i].crossings[0]) - tcycle[0]++; - if (G.dual_edges[i].crossings[1]) - tcycle[1]++; - if (G.dual_edges[e].crossings[0]) - tcycle[0]++; - if (G.dual_edges[e].crossings[1]) - tcycle[1]++; - find_cycle vis3(G, tcycle, G.dual_edges[e].v[0]); - std::vector cm3(G.dual_vertices.size()); - try { - boost::depth_first_visit( - dG, G.dual_edges[i].v[0], vis3, - boost::make_iterator_property_map(cm3.begin(), - boost::get(boost::vertex_index, dG), cm3[0])); - } catch (find_cycle::done const&) { - } - find_cycle vis4(G, tcycle, G.dual_edges[e].v[1]); - std::vector cm4(G.dual_vertices.size()); - try { - boost::depth_first_visit( - dG, G.dual_edges[i].v[1], vis4, - boost::make_iterator_property_map(cm4.begin(), - boost::get(boost::vertex_index, dG), cm4[0])); - } catch (find_cycle::done const&) { - } - other_cycles.push_back(tcycle); - } + base_sig[1]++; - if (std::any_of(other_cycles.begin(), other_cycles.end(), - [done_x, done_y](std::array c) { - return (c[0] % 2 == 0 && c[1] % 2 == 0) || - ((c[0] % 2 == 0 && done_x) || (c[1] % 2 == 0 && done_y)); - })) { - seen_pairs[{G.edges[i].v[0], G.edges[i].v[1]}] = true; - backbone[i] = true; - boost::remove_edge(G.edges[i].v[0], G.edges[i].v[1], bG); - - for (unsigned j = 0; j < 2; j++) { - std::set component_edges = {}; - get_cluster_edges visc(component_edges); - std::vector cm2(G.vertices.size()); - boost::depth_first_visit( - bG, G.edges[i].v[j], visc, - boost::make_iterator_property_map(cm2.begin(), - boost::get(boost::vertex_index, bG), cm2[0])); - - std::array total_cur = {0.0, 0.0}; - for (unsigned e : component_edges) { - if (G.edges[e].crossings[0]) { - if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) { - total_cur[0] += c[e]; - } else { - total_cur[0] -= c[e]; - } + // Then, we do it for each of the existing cycles we found in the + // previous step. + std::list> all_sigs = {base_sig}; + for (unsigned e : cedges) { + std::array new_sig_1 = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[e].v[0]); + std::array new_sig_2 = this->find_cycle(cedges, G.dual_edges[i].v[1], G.dual_edges[e].v[1]); + std::array new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]}; + + if (G.dual_edges[i].crossings[0]) + new_sig[0]++; + if (G.dual_edges[i].crossings[1]) + new_sig[1]++; + if (G.dual_edges[e].crossings[0]) + new_sig[0]++; + if (G.dual_edges[e].crossings[1]) + new_sig[1]++; + + all_sigs.push_back(new_sig); + } + + // Now, having found all cycles involving the candidate stem, we check + // if any of them spans the torus and therefore can carry current. + if (std::any_of(all_sigs.begin(), all_sigs.end(), + [done_x, done_y](std::array 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! + seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; + backbone[i] = true; + + // We're not done yet: we've severed the stem but the pop remains! We + // check each side of the lollipop and sum up all of the currents + // that span an edge. Any component without a sufficiently large net + // current must be disconnected from the current-carrying cluster and + // can also be removed from the backbone. + for (unsigned j = 0; j < 2; j++) { + std::set cluster_edges = this->get_cluster_edges(G.edges[i].v[j]); + std::array total_current = {0.0, 0.0}; + + for (unsigned e : cluster_edges) { + if (G.edges[e].crossings[0]) { + if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) { + total_current[0] += c[e]; + } else { + total_current[0] -= c[e]; } - if (G.edges[e].crossings[1]) { - if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) { - total_cur[1] += c[e]; - } else { - total_cur[1] -= c[e]; - } + } + if (G.edges[e].crossings[1]) { + if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) { + total_current[1] += c[e]; + } else { + total_current[1] -= c[e]; } } + } - if (fabs(total_cur[0]) < 1.0 / G.edges.size() && - fabs(total_cur[1]) < 1.0 / G.edges.size()) { - for (unsigned e : component_edges) { - backbone[e] = true; - boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG); - } + if (fabs(total_current[0]) < 1.0 / G.edges.size() && + fabs(total_current[1]) < 1.0 / G.edges.size()) { + for (unsigned e : cluster_edges) { + backbone[e] = true; } } } } - for (unsigned e : cedges) { - boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG); - } } } } - class get_tie_edges : public boost::default_dfs_visitor { - public: - std::set& LE; - std::set& RE; - unsigned start_edge; - bool initialized; - bool done_left; - unsigned disc; - unsigned dr; - - get_tie_edges(std::set& LE, std::set& RE) : LE(LE), RE(RE) { - initialized = false; - done_left = false; - disc = 0; - dr = 0; - } - - void discover_vertex(boost::graph_traits::vertex_descriptor vd, const Graph& g) { - if (!initialized && disc > 0) { - start_edge = vd; - initialized = true; - } - disc++; - } - - void finish_vertex(boost::graph_traits::vertex_descriptor vd, const Graph& g) { - if (vd == start_edge) { - done_left = true; - } - } - - void examine_edge(boost::graph_traits::edge_descriptor ed, const Graph& g) { - unsigned e = g[ed].index; - if (!done_left) { - LE.insert(e); - } else if (LE.find(e) == LE.end()) { - RE.insert(e); - } - } - }; - + // Now, we will check for bow ties! std::set bb_verts; for (unsigned i = 0; i < G.edges.size(); i++) { @@ -385,139 +357,75 @@ void network::update_backbone(const std::vector& c) { d1 = G.vertices[v].nd[il]; d2 = G.vertices[v].nd[ig]; - std::list cedges = {}; - get_cycle_edges vis(cedges); - std::vector cm1(G.dual_vertices.size()); - boost::depth_first_visit( - dG, d1, vis, - boost::make_iterator_property_map(cm1.begin(), - boost::get(boost::vertex_index, dG), cm1[0])); - - for (unsigned e : cedges) { - boost::remove_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], dG); - } + std::set cedges = this->get_cycle_edges(d1); - std::array cycle = {0, 0}; + std::array base_sig = this->find_cycle(cedges, d1, d2); for (unsigned k = p1; k < p2; k++) { if (G.dual_edges[G.vertices[v].ne[k]].crossings[0]) - cycle[0]++; + base_sig[0]++; if (G.dual_edges[G.vertices[v].ne[k]].crossings[1]) - cycle[1]++; + base_sig[1]++; } - find_cycle vis2(G, cycle, d2); - std::vector cm2(G.dual_vertices.size()); - try { - boost::depth_first_visit( - dG, d1, vis2, - boost::make_iterator_property_map(cm2.begin(), - boost::get(boost::vertex_index, dG), cm2[0])); - } catch (find_cycle::done const&) { - std::list> other_cycles = {cycle}; - for (unsigned e : cedges) { - std::array tcycle = {0, 0}; - for (unsigned k = p1; k < p2; k++) { - if (G.dual_edges[G.vertices[v].ne[k]].crossings[0]) - tcycle[0]++; - if (G.dual_edges[G.vertices[v].ne[k]].crossings[1]) - tcycle[1]++; - } - if (G.dual_edges[e].crossings[0]) - tcycle[0]++; - if (G.dual_edges[e].crossings[1]) - tcycle[1]++; - find_cycle vis3(G, tcycle, G.dual_edges[e].v[0]); - std::vector cm3(G.dual_vertices.size()); - try { - boost::depth_first_visit( - dG, d1, vis3, - boost::make_iterator_property_map( - cm3.begin(), boost::get(boost::vertex_index, dG), cm3[0])); - } catch (find_cycle::done const&) { - } - find_cycle vis4(G, tcycle, G.dual_edges[e].v[1]); - std::vector cm4(G.dual_vertices.size()); - try { - boost::depth_first_visit( - dG, d2, vis4, - boost::make_iterator_property_map( - cm4.begin(), boost::get(boost::vertex_index, dG), cm4[0])); - } catch (find_cycle::done const&) { - } - other_cycles.push_back(tcycle); + + // Then, we do it for each of the existing cycles we found in the + // previous step. + std::list> all_sigs = {base_sig}; + for (unsigned e : cedges) { + std::array new_sig_1 = this->find_cycle(cedges, d1, G.dual_edges[e].v[0]); + std::array new_sig_2 = this->find_cycle(cedges, d2, G.dual_edges[e].v[1]); + std::array new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]}; + + for (unsigned k = p1; k < p2; k++) { + if (G.dual_edges[G.vertices[v].ne[k]].crossings[0]) + new_sig[0]++; + if (G.dual_edges[G.vertices[v].ne[k]].crossings[1]) + new_sig[1]++; } + if (G.dual_edges[e].crossings[0]) + new_sig[0]++; + if (G.dual_edges[e].crossings[1]) + new_sig[1]++; - if (std::any_of(other_cycles.begin(), other_cycles.end(), - [done_x, done_y](std::array c) { - return ((c[0] % 2 == 0 && c[1] % 2 == 0) || - (c[0] % 2 == 0 && done_x)) || - (c[1] % 2 == 0 && done_y); - })) { - seen_pairs[{d1, d2}] = true; - - std::set left_edges = {}; - std::set right_edges = {}; - get_tie_edges visc(left_edges, right_edges); - std::vector cm2(G.vertices.size()); - boost::depth_first_visit( - bG, v, visc, - boost::make_iterator_property_map( - cm2.begin(), boost::get(boost::vertex_index, bG), cm2[0])); - - std::array total_cur_left = {0.0, 0.0}; - for (unsigned e : left_edges) { - if (G.edges[e].crossings[0]) { - if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) { - total_cur_left[0] += c[e]; - } else { - total_cur_left[0] -= c[e]; - } - } - if (G.edges[e].crossings[1]) { - if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) { - total_cur_left[1] += c[e]; - } else { - total_cur_left[1] -= c[e]; - } - } - } + all_sigs.push_back(new_sig); + } + + if (std::any_of(all_sigs.begin(), all_sigs.end(), + [done_x, done_y](std::array c) { + return ((c[0] % 2 == 0 && c[1] % 2 == 0) || + (c[0] % 2 == 0 && done_x)) || + (c[1] % 2 == 0 && done_y); + })) { + seen_pairs[{d1, d2}] = true; - std::array total_cur_right = {0.0, 0.0}; - for (unsigned e : right_edges) { + std::list> flaps = this->get_tie_flaps(v); + + for (const std::set& flap : flaps) { + std::array total_current = {0.0, 0.0}; + for (unsigned e : flap) { if (G.edges[e].crossings[0]) { if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) { - total_cur_right[0] += c[e]; + total_current[0] += c[e]; } else { - total_cur_right[0] -= c[e]; + total_current[0] -= c[e]; } } if (G.edges[e].crossings[1]) { if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) { - total_cur_right[1] += c[e]; + total_current[1] += c[e]; } else { - total_cur_right[1] -= c[e]; + total_current[1] -= c[e]; } } } - if (fabs(total_cur_left[0]) < 1.0 / G.edges.size() && - fabs(total_cur_left[1]) < 1.0 / G.edges.size()) { - for (unsigned e : left_edges) { + if (fabs(total_current[0]) < 1.0 / G.edges.size() && + fabs(total_current[1]) < 1.0 / G.edges.size()) { + for (unsigned e : flap) { backbone[e] = true; - boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG); - } - } - if (fabs(total_cur_right[0]) < 1.0 / G.edges.size() && - fabs(total_cur_right[1]) < 1.0 / G.edges.size()) { - for (unsigned e : right_edges) { - backbone[e] = true; - boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG); } } } } - for (unsigned e : cedges) { - boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG); - } } } } diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp index 2ef74c6..5f180ee 100644 --- a/src/analysis_tools.cpp +++ b/src/analysis_tools.cpp @@ -6,6 +6,7 @@ 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; } @@ -120,4 +121,5 @@ std::list, std::list>> find_minimal_ return output; } +*/ diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp index c3bffa2..2a4c5d2 100644 --- a/src/analysis_tools.hpp +++ b/src/analysis_tools.hpp @@ -11,7 +11,7 @@ template bool is_shorter(const std::list &, const std::list &); -bool trivial(boost::detail::edge_desc_impl); +//bool trivial(boost::detail::edge_desc_impl); -std::list, std::list>> find_minimal_crack(const Graph &, const network &); +//std::list, std::list>> find_minimal_crack(const Graph &, const network &); diff --git a/src/animate.cpp b/src/animate.cpp index a7b6173..6c844d2 100644 --- a/src/animate.cpp +++ b/src/animate.cpp @@ -2,7 +2,7 @@ #include "animate.hpp" #include -animate::animate(double Lx, double Ly, unsigned window_size, int argc, char *argv[]) : G(2 * (unsigned)ceil(Lx * Ly / 2)) { +animate::animate(double Lx, double Ly, unsigned window_size, int argc, char *argv[]) { glutInit(&argc, argv); glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); glutInitWindowSize((unsigned)(Lx / Ly * window_size), window_size); @@ -16,7 +16,6 @@ animate::animate(double Lx, double Ly, unsigned window_size, int argc, char *arg void animate::pre_fracture(const network &n) { lv = std::numeric_limits::lowest(); avalanches = {}; - boost::remove_edge_if(trivial, G); seen_guy = false; glClearColor(1.0f, 1.0f, 1.0f, 1.0f ); @@ -59,7 +58,6 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i) avalanches.back().push_back(i); } - boost::add_edge(n.G.dual_edges[i].v[0], n.G.dual_edges[i].v[1], {i}, G); glClearColor(1.0f, 1.0f, 1.0f, 1.0f ); glClear(GL_COLOR_BUFFER_BIT); @@ -110,6 +108,11 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i) draw = true; } + if (cur.currents[j] > 1e-9 && n.backbone[j]) { + weird = true; + } + + if (draw) { graph::coordinate r1 = n.G.vertices[n.G.edges[j].v[0]].r; graph::coordinate r2 = n.G.vertices[n.G.edges[j].v[1]].r; diff --git a/src/animate.hpp b/src/animate.hpp index 865b3bd..d14d755 100644 --- a/src/animate.hpp +++ b/src/animate.hpp @@ -9,7 +9,6 @@ class animate : public hooks { private: - Graph G; bool seen_guy; public: long double lv; diff --git a/src/animate_fracture_square.cpp b/src/animate_fracture_square.cpp index bc4c3a2..12de7da 100644 --- a/src/animate_fracture_square.cpp +++ b/src/animate_fracture_square.cpp @@ -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, false); + tmp_network.fracture(meas, true); if (quit.load()) break; diff --git a/src/measurements.cpp b/src/measurements.cpp index 984b80c..3004d7e 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -144,7 +144,6 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u } ma::ma(unsigned n, double a, double beta, double weight, bool one) : - G(2 * n), sn(2 * n), ss(2 * n), sm(2 * n), @@ -179,7 +178,6 @@ ma::ma(unsigned n, double a, double beta, double weight, bool one) : } ma::ma(unsigned Lx, unsigned Ly, double beta, double weight, bool one) : - G(Lx * Ly / 2), sn(Lx * Ly / 2), ss(Lx * Ly / 2), sm(Lx * Ly / 2), @@ -228,7 +226,6 @@ ma::~ma() { void ma::pre_fracture(const network&) { lv = std::numeric_limits::lowest(); - boost::remove_edge_if(trivial, G); avalanches = {}; num = 0; } @@ -250,14 +247,11 @@ void ma::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); num++; } void ma::post_fracture(network &n) { - auto post_cracks = find_minimal_crack(G, n); - std::vector component(boost::num_vertices(G)); - unsigned num = boost::connected_components(G, &component[0]); +/* auto post_cracks = find_minimal_crack(G, n); if (post_cracks.size() > 2 || post_cracks.size() == 0) { throw badcycleex; } @@ -271,10 +265,6 @@ void ma::post_fracture(network &n) { 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> components(num); - for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) { components[component[i]].push_back(n.G.dual_vertices[i].r); } @@ -336,6 +326,7 @@ void ma::post_fracture(network &n) { } 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 b22c327..879f511 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -19,8 +19,6 @@ class ma : public hooks { // - interface for turning on and off specific measurements // private: - Graph G; - unsigned num; // measurement storage -- cgit v1.2.3-54-g00ecf