From 456d748b8cbf141e1fe6e82a7c0106f469e28ba6 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 7 May 2019 14:05:32 -0400 Subject: added measurements of spanning cluster sizes and conducting backbone sizes --- lib/include/network.hpp | 8 ++++-- lib/src/network.cpp | 68 ++++++++++++++++++++++++++++++------------------- src/measurements.cpp | 22 ++++++++++++++-- 3 files changed, 68 insertions(+), 30 deletions(-) diff --git a/lib/include/network.hpp b/lib/include/network.hpp index 8f7bd12..dd11342 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -54,6 +54,8 @@ class network { network(const network&); void set_thresholds(double, std::mt19937&); + virtual void break_edge(unsigned e, bool unbreak = false) {fuses[e] = !unbreak;}; + virtual current_info get_current_info() {current_info empty; return empty;}; }; class fuse_network : public network { @@ -63,10 +65,11 @@ class fuse_network : public network { fuse_network(const graph&, cholmod_common*); void fracture(hooks&, double cutoff = 1e-13); - current_info get_current_info(); }; class elastic_network : public network { + private: + double weight; public: problem hook_x; problem hook_y; @@ -75,6 +78,7 @@ class elastic_network : public network { elastic_network(const elastic_network&); void fracture(hooks&, double weight = 0.5, double cutoff = 1e-10); - current_info get_current_info(); + 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 e885fcb..ed51c8e 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -290,36 +290,16 @@ elastic_network::elastic_network(const elastic_network& n) : network(n), hook_x( } void elastic_network::fracture(hooks& m, double weight, double cutoff) { + this->weight = weight; m.pre_fracture(*this); while (true) { - current_info cx = hook_x.solve(fuses); - current_info cy = hook_y.solve(fuses); + current_info ctot = this->get_current_info(); - bool done_x = cx.conductivity < 1.0 / G.edges.size(); - bool done_y = cy.conductivity < 1.0 / G.edges.size(); - - if (done_x && done_y) { + if (ctot.conductivity == 0) { break; } - current_info ctot; - ctot.currents.resize(G.edges.size()); - - if (done_x) { - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity; - } - } else if (done_y) { - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity; - } - } else { - 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)); - } - } - unsigned max_pos = UINT_MAX; long double max_val = std::numeric_limits::lowest(); @@ -337,9 +317,7 @@ void elastic_network::fracture(hooks& m, double weight, double cutoff) { throw nofuseex; } - fuses[max_pos] = true; - hook_x.break_edge(max_pos); - hook_y.break_edge(max_pos); + this->break_edge(max_pos); m.bond_broken(*this, ctot, max_pos); } @@ -347,3 +325,41 @@ void elastic_network::fracture(hooks& m, double weight, double cutoff) { m.post_fracture(*this); } +current_info elastic_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 elastic_network::break_edge(unsigned e, bool unbreak) { + fuses[e] = !unbreak; + hook_x.break_edge(e, unbreak); + hook_y.break_edge(e, unbreak); +} + diff --git a/src/measurements.cpp b/src/measurements.cpp index 407048d..33056ef 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -148,7 +148,9 @@ ma::ma(unsigned n, double a, unsigned Mx, unsigned My, double beta) : sm(2 * n), sa(3 * n), sl(2 * n), - sn(2 * n) + sn(2 * n), + ss(2 * n), + sb(3 * n) { if (beta != 0.0) { model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_"; @@ -178,6 +180,8 @@ ma::~ma() { update_distribution_file("sa", sa, model_string); update_distribution_file("sl", sl, model_string); update_distribution_file("sn", sn, model_string); + update_distribution_file("ss", ss, model_string); + update_distribution_file("sb", sb, model_string); } void ma::pre_fracture(const network&) { @@ -219,6 +223,8 @@ void ma::post_fracture(network &n) { for (unsigned i = 0; i < num; i++) { if (i != crack_component) { sm[components[i].size() - 1]++; + } else { + ss[components[i].size() - 1]++; } sn[components[i].size() - 1]++; } @@ -228,7 +234,7 @@ void ma::post_fracture(network &n) { while (true) { for (unsigned e : *av_it) { boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G); - n.fuses[e] = false; + n.break_edge(e, true); } auto cracks = find_minimal_crack(G, n); @@ -252,6 +258,18 @@ void ma::post_fracture(network &n) { sc[new_components[i].size() - 1]++; } + current_info ct = n.get_current_info(); + + unsigned conducting_backbone_size = 0; + + for (unsigned i = 0; i < n.G.edges.size(); i++) { + if (ct.currents[i] > 1.0 / n.G.edges.size()) { + conducting_backbone_size++; + } + } + + sb[conducting_backbone_size - 1]++; + av_it++; while (av_it != avalanches.rend()) { -- cgit v1.2.3-54-g00ecf