summaryrefslogtreecommitdiff
path: root/lib/src
diff options
context:
space:
mode:
Diffstat (limited to 'lib/src')
-rw-r--r--lib/src/network.cpp80
1 files changed, 80 insertions, 0 deletions
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);
+}
+