diff options
Diffstat (limited to 'lib')
-rw-r--r-- | lib/include/network.hpp | 15 | ||||
-rw-r--r-- | lib/src/network.cpp | 80 |
2 files changed, 95 insertions, 0 deletions
diff --git a/lib/include/network.hpp b/lib/include/network.hpp index dd11342..15cb8ea 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -82,3 +82,18 @@ class elastic_network : public network { current_info get_current_info() override; }; +class percolation_network : public network { + private: + double weight; + public: + problem hook_x; + problem hook_y; + + percolation_network(const graph&, cholmod_common*); + percolation_network(const percolation_network&); + + void fracture(hooks&, double weight = 0.5, double cutoff = 1e-10); + void break_edge(unsigned, bool unbreak = false) override; + current_info get_current_info() override; +}; + diff --git a/lib/src/network.cpp b/lib/src/network.cpp index ed51c8e..1d9b3f9 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -363,3 +363,83 @@ void elastic_network::break_edge(unsigned e, bool unbreak) { hook_y.break_edge(e, unbreak); } +percolation_network::percolation_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) { +} + +percolation_network::percolation_network(const percolation_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) { +} + +void percolation_network::fracture(hooks& m, double weight, double cutoff) { + this->weight = weight; + m.pre_fracture(*this); + + while (true) { + current_info ctot = this->get_current_info(); + + if (ctot.conductivity == 0) { + break; + } + + unsigned max_pos = UINT_MAX; + long double max_val = std::numeric_limits<long double>::lowest(); + + for (unsigned i = 0; i < G.edges.size(); i++) { + if (!fuses[i] && ctot.currents[i] > 1.0 / G.edges.size()) { + long double val = - thresholds[i]; + if (val > max_val) { + max_val = val; + max_pos = i; + } + } + } + + if (max_pos == UINT_MAX) { + throw nofuseex; + } + + this->break_edge(max_pos); + + m.bond_broken(*this, ctot, max_pos); + } + + m.post_fracture(*this); +} + +current_info percolation_network::get_current_info() { + current_info cx = hook_x.solve(fuses); + current_info cy = hook_y.solve(fuses); + + bool done_x = cx.conductivity < 1.0 / G.edges.size(); + bool done_y = cy.conductivity < 1.0 / G.edges.size(); + + current_info ctot; + ctot.currents.resize(G.edges.size()); + + if (done_x && done_y) { + ctot.conductivity = 0; + } else if (done_x) { + ctot.conductivity = cy.conductivity; + for (unsigned i = 0; i < G.edges.size(); i++) { + ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity; + } + } else if (done_y) { + ctot.conductivity = cx.conductivity; + for (unsigned i = 0; i < G.edges.size(); i++) { + ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity; + } + } else { + ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); + for (unsigned i = 0; i < G.edges.size(); i++) { + ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2)); + } + } + + return ctot; +} + +void percolation_network::break_edge(unsigned e, bool unbreak) { + fuses[e] = !unbreak; + hook_x.break_edge(e, unbreak); + hook_y.break_edge(e, unbreak); +} + |