fixed problems with the current cutoff

2 files changed, 18 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<long double>::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;