diff options
Diffstat (limited to 'lib/src/network.cpp')
| -rw-r--r-- | lib/src/network.cpp | 78 |
1 files changed, 65 insertions, 13 deletions
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<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edge return sig; } +bool network::get_cycle_helper(std::array<unsigned, 2>& sig, std::set<unsigned>& edges, const std::set<unsigned>& 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<unsigned, 2>& 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::array<unsigned, 2>, std::set<unsigned>> network::get_cycle(const std::set<unsigned>& cycle_edges, unsigned v0, + unsigned v1) const { + std::set<unsigned> edges; + std::array<unsigned, 2> 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<unsigned>& 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<double>& 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<double>& c) { } } } + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), + [](std::array<unsigned, 2> 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<double>& c) { } } } + } else if (std::any_of(all_sigs.begin(), all_sigs.end(), + [](std::array<unsigned, 2> c) { + return ((c[0] % 2 == 0 && c[1] % 2 == 1) || (c[0] % 2 == 1 && c[1] % 2 == 0)); + })) { + seen_pairs[{d1, d2}] = true; } } } |