From 7f7bfc7789c7c0d5b8a8daeda82b349ee1ea52a0 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sat, 27 Apr 2019 18:02:58 -0400 Subject: fixed problems with the current cutoff --- lib/include/network.hpp | 2 +- lib/src/network.cpp | 29 +++++++++++++++++------------ src/measurements.cpp | 3 +++ 3 files changed, 21 insertions(+), 13 deletions(-) diff --git a/lib/include/network.hpp b/lib/include/network.hpp index ba98086..d95b3c3 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -74,7 +74,7 @@ class elastic_network : public network { elastic_network(const graph&, cholmod_common*); elastic_network(const elastic_network&); - void fracture(hooks&, double weight = 0.5, double cutoff = 1e-13); + void fracture(hooks&, double weight = 0.5, double cutoff = 1e-11); current_info get_current_info(); }; diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 712bf07..0cacdf9 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -296,21 +296,26 @@ void elastic_network::fracture(hooks& m, double weight, double cutoff) { current_info cx = hook_x.solve(fuses); current_info cy = hook_y.solve(fuses); - current_info ctot; - - ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); + bool done_x = cx.conductivity < cutoff * G.vertices.size() || cx.conductivity != cx.conductivity; + bool done_y = cy.conductivity < cutoff * G.vertices.size() || cy.conductivity != cy.conductivity; - if (ctot.conductivity < cutoff * G.vertices.size()) { + if (done_x && done_y) { break; } + current_info ctot; ctot.currents.resize(G.edges.size()); - for (unsigned i = 0; i < G.edges.size(); i++) { - if (cx.conductivity < cutoff * G.vertices.size()) { - ctot.currents[i] = cy.currents[i] / cy.conductivity; - } else if (cy.conductivity < cutoff * G.vertices.size()) { - ctot.currents[i] = cx.currents[i] / cx.conductivity; - } else { + + 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)); } } @@ -319,8 +324,8 @@ void elastic_network::fracture(hooks& m, double weight, double cutoff) { long double max_val = std::numeric_limits::lowest(); for (unsigned i = 0; i < G.edges.size(); i++) { - if (!fuses[i] && fabs(ctot.currents[i]) > cutoff) { - long double val = logl(fabs(ctot.currents[i])) - thresholds[i]; + if (!fuses[i] && ctot.currents[i] > cutoff) { + long double val = logl(ctot.currents[i]) - thresholds[i]; if (val > max_val) { max_val = val; max_pos = i; diff --git a/src/measurements.cpp b/src/measurements.cpp index 5a79bd7..bf053e9 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -187,6 +187,9 @@ void ma::post_fracture(network &n) { auto post_cracks = find_minimal_crack(G, n); std::vector component(boost::num_vertices(G)); unsigned num = boost::connected_components(G, &component[0]); + if (post_cracks.size() > 2 || post_cracks.size() == 0) { + throw badcycleex; + } unsigned crack_component = component[n.G.dual_edges[post_cracks.front().second.front()].v[0]]; std::vector> components(num); -- cgit v1.2.3-70-g09d2