summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2019-05-17 16:38:24 -0400
committerJaron Kent-Dobias <jaron@kent-dobias.com>2019-05-17 16:38:24 -0400
commit1e939e597964fa081b347e40af2be1069867b906 (patch)
tree395de78337d7ed85971fa3e84860ee4f4e5a5956 /lib
parent456d748b8cbf141e1fe6e82a7c0106f469e28ba6 (diff)
downloadfuse_networks-1e939e597964fa081b347e40af2be1069867b906.tar.gz
fuse_networks-1e939e597964fa081b347e40af2be1069867b906.tar.bz2
fuse_networks-1e939e597964fa081b347e40af2be1069867b906.zip
added support for running percolation using the fuse network rules
Diffstat (limited to 'lib')
-rw-r--r--lib/include/network.hpp15
-rw-r--r--lib/src/network.cpp80
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);
+}
+