summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/include/network.hpp8
-rw-r--r--lib/src/network.cpp68
-rw-r--r--src/measurements.cpp22
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<long double>::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()) {