From 2faf0e4598c7c046d58107d23145f95db334200c Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 16 Mar 2020 18:19:09 -0400 Subject: Implemented new behavior when w=0 or w=1. --- lib/include/network.hpp | 10 ++-- lib/src/network.cpp | 130 +++++++++++++++++++++++------------------------ src/analysis_tools.cpp | 20 ++++++-- src/animate_fracture.cpp | 28 ++++++++++ src/fracture.cpp | 44 ++++++++++++++++ 5 files changed, 157 insertions(+), 75 deletions(-) diff --git a/lib/include/network.hpp b/lib/include/network.hpp index d562530..37f1912 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -31,9 +31,10 @@ class network { private: problem px; problem py; - std::unordered_map, bool> seen_pairs; + std::unordered_map, bool> seen_pairs_x; + std::unordered_map, bool> seen_pairs_y; - void update_backbone(const std::vector& c); + void update_backbone(const std::vector& c, const std::array& cond); void break_edge(unsigned, bool unbreak = false); void get_cycle_edges_helper(std::set&, std::set&, unsigned, unsigned) const; bool find_cycle_helper(std::array&, const std::set&, unsigned, unsigned, @@ -47,6 +48,7 @@ private: std::list> get_tie_flaps(unsigned) const; public: + bool two_sides; const graph& G; ClusterTree C; @@ -55,7 +57,7 @@ public: std::vector backbone; std::vector thresholds; - network(const graph&, cholmod_common*); + network(const graph&, bool two_sides, cholmod_common*); void set_thresholds(double, std::mt19937&); void fracture(hooks&); @@ -76,7 +78,7 @@ private: double weight; public: - fuse_network(const graph&, cholmod_common*, double weight = 1.0); + fuse_network(const graph&, cholmod_common*); fuse_network(const fuse_network&); current_info get_current_info() override; diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 561796b..8944c15 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -1,14 +1,16 @@ #include "network.hpp" #include +#include class nofuseException : public std::exception { virtual const char* what() const throw() { return "No valid fuse was available to break."; } } nofuseex; -network::network(const graph& G, cholmod_common* 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()) {} +network::network(const graph& G, bool two_sides, cholmod_common* c) + : px(G, 0, c), py(G, 1, c), two_sides(two_sides), 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) { @@ -38,14 +40,16 @@ void network::fracture(hooks& m) { current_info c = this->get_current_info(); - if (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) { + if ((c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) && two_sides) { + break; + } else if (c.conductivity[0] < min_cond && !two_sides) { break; } if (px.sol.conductivity[0] > min_cond) { - this->update_backbone(px.sol.currents); + this->update_backbone(px.sol.currents, c.conductivity); } else { - this->update_backbone(py.sol.currents); + this->update_backbone(py.sol.currents, c.conductivity); } unsigned max_pos = UINT_MAX; @@ -74,8 +78,8 @@ void network::fracture(hooks& m) { } void network::get_cycle_edges_helper(std::set& cycle_edges, - std::set& seen_vertices, unsigned v_prev, - unsigned v_cur) const { + 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]) { @@ -103,7 +107,7 @@ 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 { + unsigned vPrev, unsigned vCur, unsigned vEnd) const { for (unsigned ei : G.dual_vertices[vCur].ne) { if (fuses[ei]) { if (!cycle_edges.contains(ei)) { @@ -134,14 +138,15 @@ bool network::find_cycle_helper(std::array& sig, const std::set network::find_cycle(const std::set& cycle_edges, unsigned v0, - unsigned v1) const { + unsigned v1) const { std::array sig = {0, 0}; this->find_cycle_helper(sig, cycle_edges, v0, v0, v1); 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 { +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]) { if (!cycle_edges.contains(ei)) { @@ -173,8 +178,8 @@ bool network::get_cycle_helper(std::array& sig, std::set& return false; } -std::pair, std::set> network::get_cycle(const std::set& cycle_edges, unsigned v0, - unsigned v1) const { +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); @@ -202,7 +207,7 @@ std::set network::get_cluster_edges(unsigned v0) const { } void network::get_tie_flaps_helper(std::set& added_edges, unsigned v0, - unsigned vCur) const { + unsigned vCur) const { for (unsigned ei : G.vertices[vCur].ne) { if (!backbone[ei]) { if (!added_edges.contains(ei)) { @@ -245,14 +250,25 @@ std::list> network::get_tie_flaps(unsigned v0) const { 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(); +void network::update_backbone(const std::vector& c, const std::array& cc) { + bool done_x = cc[0] < 1.0 / G.edges.size(); + bool done_y = cc[1] < 1.0 / G.edges.size(); + + auto cycle_check_current = [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)); + }; + auto cycle_check_wrap_x = [](std::array c) { + return (c[0] % 2 == 0 && c[1] % 2 == 1); + }; + auto cycle_check_wrap_y = [](std::array c) { + return (c[0] % 2 == 1 && c[1] % 2 == 0); + }; // 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])) { - if (!seen_pairs.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]})) { + if ((!seen_pairs_x.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]}) && !done_x) || (!seen_pairs_y.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]}) && !done_y)) { // 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. @@ -261,7 +277,7 @@ void network::update_backbone(const std::vector& c) { // 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]); + this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[i].v[1]); if (G.dual_edges[i].crossings[0]) base_sig[0]++; if (G.dual_edges[i].crossings[1]) @@ -272,11 +288,11 @@ void network::update_backbone(const std::vector& c) { 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]); + 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]); + 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]}; + new_sig_1[1] + new_sig_2[1]}; if (G.dual_edges[i].crossings[0]) new_sig[0]++; @@ -292,13 +308,10 @@ void network::update_backbone(const std::vector& c) { // 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 (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_current)) { // If none does, we remove it from the backbone! - seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; + seen_pairs_x[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; + seen_pairs_y[{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 @@ -334,14 +347,13 @@ 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)); - })) { + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_x)) { // 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; + seen_pairs_x[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_y)) { + seen_pairs_y[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true; } } } @@ -371,7 +383,7 @@ void network::update_backbone(const std::vector& c) { if (C.same_component(G.vertices[v].nd[i], G.vertices[v].nd[l])) { unsigned il = i < l ? i : l; unsigned ig = i < l ? l : i; - if (!seen_pairs.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]})) { + if ((!seen_pairs_x.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]}) && !done_x) || (!seen_pairs_y.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]}) && !done_y)) { bool any_intervening1 = false; bool any_intervening2 = false; unsigned ie1 = 0; @@ -449,14 +461,10 @@ void network::update_backbone(const std::vector& c) { 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); - })) { + if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_current)) { // one of our pairs qualifies for trimming! - seen_pairs[{d1, d2}] = true; + seen_pairs_x[{d1, d2}] = true; + seen_pairs_y[{d1, d2}] = true; // We separate the flaps of the bowtie (there may be more than // two!) into distinct sets. @@ -490,11 +498,10 @@ 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; + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_x)) { + seen_pairs_x[{d1, d2}] = true; + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_y)) { + seen_pairs_y[{d1, d2}] = true; } } } @@ -520,7 +527,7 @@ std::string network::write() { output += "\"fuses\"->{"; for (unsigned i = 0; i < G.edges.size(); i++) { if (!fuses[i]) { - output += std::to_string(i) +","; + output += std::to_string(i) + ","; } } output.pop_back(); @@ -536,7 +543,8 @@ std::string network::write() { output += std::to_string(t) + ","; } output.pop_back(); - output += "},\"conductivity\"->{" + std::to_string(c.conductivity[0]) + "," + std::to_string(c.conductivity[1]); + output += "},\"conductivity\"->{" + std::to_string(c.conductivity[0]) + "," + + std::to_string(c.conductivity[1]); output += "},\"currents\"->{"; for (const double& t : c.currents) { output += std::to_string(t) + ","; @@ -547,35 +555,23 @@ std::string network::write() { return output; }; +fuse_network::fuse_network(const graph& G, cholmod_common* c) : network(G, false, c), weight(1.0) {} -fuse_network::fuse_network(const graph& G, cholmod_common* c, double weight) - : network(G, c), weight(weight) {} - -fuse_network::fuse_network(const fuse_network& n) : network(n), weight(n.weight) {} +fuse_network::fuse_network(const fuse_network& n) : network(n), weight(1.0) {} current_info fuse_network::get_current_info() { px.solve(fuses); py.solve(fuses); bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size(); - bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size(); current_info ctot; ctot.currents.resize(G.edges.size()); ctot.conductivity = {px.sol.conductivity[0], py.sol.conductivity[1]}; - if (done_x && !done_y) { - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = fabs(weight * py.sol.currents[i] / py.sol.conductivity[1]); - } - } else if (done_y && !done_x) { - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0]); - } - } else if (!done_x && !done_y) { + if (!done_x) { for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0] + - weight * py.sol.currents[i] / py.sol.conductivity[1]); + ctot.currents[i] = fabs(px.sol.currents[i] / px.sol.conductivity[0]); } } @@ -583,7 +579,7 @@ current_info fuse_network::get_current_info() { } elastic_network::elastic_network(const graph& G, cholmod_common* c, double weight) - : network(G, c), weight(weight) {} + : network(G, true, c), weight(weight) {} elastic_network::elastic_network(const elastic_network& n) : network(n), weight(n.weight) {} @@ -617,7 +613,7 @@ current_info elastic_network::get_current_info() { return ctot; } -percolation_network::percolation_network(const graph& G, cholmod_common* c) : network(G, c) {} +percolation_network::percolation_network(const graph& G, cholmod_common* c) : network(G, true, c) {} percolation_network::percolation_network(const percolation_network& n) : network(n) {} diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp index 2965915..8ad8235 100644 --- a/src/analysis_tools.cpp +++ b/src/analysis_tools.cpp @@ -30,10 +30,22 @@ std::pair, std::set> find_minimal_crack(const } if (all_good) { - if (cycles[0].second.size() > cycles[1].second.size()) { - return cycles[1]; + if (n.two_sides) { + if (cycles[0].second.size() > cycles[1].second.size()) { + return cycles[1]; + } else { + return cycles[0]; + } } else { - return cycles[0]; + if (cycles[0].first[0] % 2 == 0) { + return cycles[0]; + } else { + return cycles[1]; + } + } + } else if (!n.two_sides) { + if (cycles[!not_good].first[0] % 2 == 0) { + return cycles[!not_good]; } } @@ -63,7 +75,7 @@ std::pair, std::set> find_minimal_crack(const pos++; } - if (cycles[!not_good].second.size() > new_cycle.size()) { + if (cycles[!not_good].second.size() > new_cycle.size() || !n.two_sides) { return {{sum_sig_0, sum_sig_1}, new_cycle}; } else { return cycles[!not_good]; diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp index e8de8c1..fd1ffe2 100644 --- a/src/animate_fracture.cpp +++ b/src/animate_fracture.cpp @@ -69,6 +69,33 @@ int main(int argc, char* argv[]) { randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; + if (w == 0.0 || w == 1.0) { + if (use_aN) { + animate meas(sqrt(2*n *a), sqrt( 2*n / a), width, argc, argv); + + for (unsigned trial = 0; trial < N; trial++) { + graph G(n, a, rng); + fuse_network net(G, &c); + net.set_thresholds(beta, rng); + net.fracture(meas); + } + } else { + animate meas(Lx, Ly, width, argc, argv); + const graph G(Lx, Ly); + const fuse_network plain_net(G, &c); + + for (unsigned trial = 0; trial < N; trial++) { + fuse_network net = plain_net; + net.set_thresholds(beta, rng); + try { + net.fracture(meas); + } catch (std::exception &e) { + std::cout << e.what() << std::endl; + getchar(); + } + } + } + } else { if (use_aN) { animate meas(sqrt(2*n *a), sqrt( 2*n / a), width, argc, argv); @@ -89,6 +116,7 @@ int main(int argc, char* argv[]) { net.fracture(meas); } } + } CHOL_F(finish)(&c); diff --git a/src/fracture.cpp b/src/fracture.cpp index 96cdac3..22101c8 100644 --- a/src/fracture.cpp +++ b/src/fracture.cpp @@ -76,6 +76,49 @@ int main(int argc, char* argv[]) { randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; + if (w == 0 || w == 1) { + if (use_aN) { + ma meas(n, a, beta, w); + + for (unsigned trial = 0; trial < N; trial++) { + while (true) { + try { + graph G(n, a, rng); + fuse_network n(G, &c); + n.set_thresholds(beta, rng); + n.fracture(meas); + break; + } catch (std::exception &e) { + std::cout << e.what() << '\n'; + } + } + + if (quit.load()) + break; + } + } else { + ma meas(Lx, Ly, beta, w); + + const graph G(Lx, Ly); + const fuse_network n(G, &c); + + for (unsigned trial = 0; trial < N; trial++) { + while (true) { + try { + fuse_network net = n; + net.set_thresholds(beta, rng); + net.fracture(meas); + break; + } catch (std::exception &e) { + std::cout << e.what() << '\n'; + } + } + + if (quit.load()) + break; + } + } + } else { if (use_aN) { ma meas(n, a, beta, w); @@ -117,6 +160,7 @@ int main(int argc, char* argv[]) { break; } } + } CHOL_F(finish)(&c); -- cgit v1.2.3-54-g00ecf