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);
+}
+