more progress on crack stress

4 files changed, 279 insertions, 221 deletions

diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index c78fe9f..5ad4c8b 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -13,13 +13,44 @@
 
 #include "array_hash.hpp"
 
-class graph {
+class coordinate {
 public:
-  typedef struct coordinate {
-    double x;
-    double y;
-  } coordinate;
+  double x;
+  double y;
+
+  coordinate operator+(const coordinate& r) const {
+    coordinate rp = {0, 0};
+    rp.x = x + r.x;
+    rp.y = y + r.y;
+
+    return rp;
+  }
+
+  template <class T>
+  coordinate operator*(const T& a) const {
+    coordinate rp = *this;
+    rp.x *= a;
+    rp.y *= a;
+    return rp;
+  }
+
+  coordinate operator-(const coordinate& r) const {
+    coordinate rp = *this;
+    return rp + (r * -1);
+  }
+
+  coordinate operator-() const {
+    return *this * -1;
+  }
 
+  void operator+=(const coordinate& r) {
+    x += r.x;
+    y += r.y;
+  }
+};
+
+class graph {
+public:
   typedef struct vertex {
     coordinate r;
     std::vector<unsigned> nd;
@@ -49,4 +80,6 @@ public:
   graph const rotate();
 
   std::string write() const;
+
+  coordinate Δr(unsigned e) const;
 };
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index c75d466..1c770e9 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -445,14 +445,14 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
   jcv_diagram_free(&diagram);
 }
 
-graph::coordinate reverse(const graph::coordinate& x) { return {x.y, x.x}; }
+coordinate reverse(const coordinate& x) { return {x.y, x.x}; }
 
 graph const graph::rotate() {
   graph g2(*this);
 
   for (graph::vertex& v : g2.vertices) {
     v.r = reverse(v.r);
-    for (graph::coordinate& r : v.polygon) {
+    for (coordinate& r : v.polygon) {
       r = reverse(r);
     }
   }
@@ -464,7 +464,7 @@ graph const graph::rotate() {
 
   for (graph::vertex& v : g2.dual_vertices) {
     v.r = reverse(v.r);
-    for (graph::coordinate& r : v.polygon) {
+    for (coordinate& r : v.polygon) {
       r = reverse(r);
     }
   }
@@ -490,7 +490,7 @@ std::string graph::write() const {
   output += "},\"vp\"->{";
   for (const graph::vertex &v : vertices) {
     output += "{";
-    for (const graph::coordinate &r : v.polygon) {
+    for (const coordinate &r : v.polygon) {
       output += "{" + std::to_string(r.x) + "," + std::to_string(r.y) + "},";
     }
     output.pop_back();
@@ -505,7 +505,7 @@ std::string graph::write() const {
   output += "},\"up\"->{";
   for (const graph::vertex &v : dual_vertices) {
     output += "{";
-    for (const graph::coordinate &r : v.polygon) {
+    for (const coordinate &r : v.polygon) {
       output += "{" + std::to_string(r.x) + "," + std::to_string(r.y) + "},";
     }
     output.pop_back();
@@ -536,3 +536,23 @@ std::string graph::write() const {
 
   return output;
 }
+
+coordinate graph::Δr(unsigned e) const {
+  coordinate tmp = dual_vertices[dual_edges[e].v[1]].r - dual_vertices[dual_edges[e].v[0]].r;
+  if (dual_edges[e].crossings[0]) {
+    if (tmp.x > 0) {
+      tmp.x -= L.x;
+    } else {
+      tmp.x += L.x;
+    }
+  }
+  if (dual_edges[e].crossings[1]) {
+    if (tmp.y > 0) {
+      tmp.y -= L.y;
+    } else {
+      tmp.y += L.y;
+    }
+  }
+
+  return tmp;
+}
diff --git a/lib/src/problem.cpp b/lib/src/problem.cpp
index ab87b0c..a029d66 100644
--- a/lib/src/problem.cpp
+++ b/lib/src/problem.cpp
@@ -9,8 +9,8 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse* vcmat, cholmod_c
     : G(G), axis(axis), voltcurmat(vcmat), c(c) {
   b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c);
   for (unsigned i = 0; i < G.edges.size(); i++) {
-    graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
-    graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
+    coordinate v0 = G.vertices[G.edges[i].v[0]].r;
+    coordinate v1 = G.vertices[G.edges[i].v[1]].r;
 
     if (G.edges[i].crossings[axis]) {
       bool ind;
@@ -131,8 +131,8 @@ void problem::solve(std::vector<bool>& fuses) {
     } else {
       sol.currents[i] = ((double*)y->x)[i];
 
-      graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
-      graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
+      coordinate v0 = G.vertices[G.edges[i].v[0]].r;
+      coordinate v1 = G.vertices[G.edges[i].v[1]].r;
 
       if (G.edges[i].crossings[axis]) {
         bool comp;
@@ -196,8 +196,8 @@ void problem::break_edge(unsigned e, bool unbreak) {
   CHOL_F(free_sparse)(&perm_update_mat, c);
   CHOL_F(free_sparse)(&update_mat, c);
 
-  graph::coordinate r0 = G.vertices[v0].r;
-  graph::coordinate r1 = G.vertices[v1].r;
+  coordinate r0 = G.vertices[v0].r;
+  coordinate r1 = G.vertices[v1].r;
 
   if (G.edges[e].crossings[axis]) {
     bool ind;
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 16f1225..e7aa528 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -37,7 +37,8 @@ void update_distribution_file(std::string id, const std::vector<uint64_t>& data,
 }
 
 template <class T>
-void update_field_file(std::string id, uint64_t N, const std::vector<T>& data, std::string model_string) {
+void update_field_file(std::string id, uint64_t N, const std::vector<T>& data,
+                       std::string model_string) {
   std::string filename = model_string + id + ".dat";
   std::ifstream file(filename);
 
@@ -100,9 +101,9 @@ tx1[i][1] * tx2[i][1]; fftw_reverse_in[i][1] = tx1[i][0] * tx2[i][1] - tx1[i][1]
 */
 
 void autocorrelation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data,
-                      const std::list<graph::coordinate>& pos, std::array<unsigned, 2> count) {
-  for (std::list<graph::coordinate>::const_iterator it1 = pos.begin(); it1 != pos.end(); it1++) {
-    for (std::list<graph::coordinate>::const_iterator it2 = it1; it2 != pos.end(); it2++) {
+                      const std::list<coordinate>& pos, std::array<unsigned, 2> count) {
+  for (std::list<coordinate>::const_iterator it1 = pos.begin(); it1 != pos.end(); it1++) {
+    for (std::list<coordinate>::const_iterator it2 = it1; it2 != pos.end(); it2++) {
       double Δx_tmp = fabs(it1->x - it2->x);
       double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp;
 
@@ -110,19 +111,21 @@ void autocorrelation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vecto
       double Δy = Δy_tmp < Ly / 2 ? Δy_tmp : Ly - Δy_tmp;
 
       if (count[0] % 2 == 0) {
-        data[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δx / Lx)) + (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δy / Ly))]++;
+        data[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δx / Lx)) +
+             (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δy / Ly))]++;
       } else {
-        data[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δy / Ly)) + (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δx / Lx))]++;
+        data[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δy / Ly)) +
+             (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δx / Lx))]++;
       }
     }
   }
 }
 
 void correlation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data,
-                  const std::list<graph::coordinate>& pos1,
-                  const std::list<graph::coordinate>& pos2) {
-  for (std::list<graph::coordinate>::const_iterator it1 = pos1.begin(); it1 != pos1.end(); it1++) {
-    for (std::list<graph::coordinate>::const_iterator it2 = pos2.begin(); it2 != pos2.end();
+                  const std::list<coordinate>& pos1,
+                  const std::list<coordinate>& pos2) {
+  for (std::list<coordinate>::const_iterator it1 = pos1.begin(); it1 != pos1.end(); it1++) {
+    for (std::list<coordinate>::const_iterator it2 = pos2.begin(); it2 != pos2.end();
          it2++) {
       double Δx_tmp = fabs(it1->x - it2->x);
       double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp;
@@ -157,16 +160,18 @@ pow(fftw_forward_out[i][1], 2); fftw_reverse_in[i][1] = 0.0;
 }
 */
 
-unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, unsigned My) {
+unsigned edge_r_to_ind(coordinate r, double Lx, double Ly, unsigned Mx, unsigned My) {
   return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);
 }
 
 ma::ma(unsigned n, double a, double beta, double weight)
     : sc(2 * n), sn(2 * n), ss(2 * n), sm(2 * n), sl(2 * n), sb(n + 1), sd(3 * n), sa(3 * n),
-      sA(3 * n), sp(3 * n), si(10000), sI(10000), cc(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cn(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
-      cs(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cm(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cl(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
-      cb(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), ca(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cA(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
-      cq(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cS(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), last_clusters(2 * n) {
+      sA(3 * n), sp(3 * n), si(10000), sI(10000), cc(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
+      cn(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), cs(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
+      cm(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), cl(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
+      cb(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), ca(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
+      cA(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), cq(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
+      cS(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), last_clusters(2 * n) {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" +
                    std::to_string(beta) + "_" + std::to_string(weight) + "_";
@@ -276,7 +281,7 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
 void ma::post_fracture(network& n) {
   auto crack = find_minimal_crack(n, avalanches.back().back());
 
-  std::list<graph::coordinate> cl_cs;
+  std::list<coordinate> cl_cs;
   sl[crack.second.size() - 1]++;
   for (unsigned e : crack.second) {
     cl_cs.push_back(n.G.dual_edges[e].r);
@@ -284,7 +289,7 @@ void ma::post_fracture(network& n) {
   Nl += crack.second.size();
   autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cl, cl_cs, crack.first);
 
-  std::vector<std::list<graph::coordinate>> crit_components(n.G.dual_vertices.size());
+  std::vector<std::list<coordinate>> crit_components(n.G.dual_vertices.size());
 
   for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
     crit_components[last_clusters.findroot(i)].push_back(n.G.dual_vertices[i].r);
@@ -298,244 +303,244 @@ void ma::post_fracture(network& n) {
   }
   Nc += n.G.dual_vertices.size();
 
-    std::vector<std::list<graph::coordinate>> components(n.G.dual_vertices.size());
+  std::vector<std::list<coordinate>> components(n.G.dual_vertices.size());
 
-    for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
-      components[n.C.findroot(i)].push_back(n.G.dual_vertices[i].r);
-    }
+  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
+    components[n.C.findroot(i)].push_back(n.G.dual_vertices[i].r);
+  }
 
-    unsigned crack_component = n.C.findroot(n.G.dual_edges[avalanches.back().back()].v[0]);
+  unsigned crack_component = n.C.findroot(n.G.dual_edges[avalanches.back().back()].v[0]);
 
-    for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
-      if (components[i].size() > 0) {
-        if (i != crack_component) {
-          sm[components[i].size() - 1]++;
-          autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cm, components[i], crack.first);
-        } else {
-          ss[components[i].size() - 1]++;
-          autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cs, components[i], crack.first);
-          Ns += components[i].size();
-        }
-        sn[components[i].size() - 1]++;
-        autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cn, components[i], crack.first);
+  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
+    if (components[i].size() > 0) {
+      if (i != crack_component) {
+        sm[components[i].size() - 1]++;
+        autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cm, components[i], crack.first);
+      } else {
+        ss[components[i].size() - 1]++;
+        autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cs, components[i], crack.first);
+        Ns += components[i].size();
       }
+      sn[components[i].size() - 1]++;
+      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cn, components[i], crack.first);
     }
-    Nn += n.G.dual_vertices.size();
+  }
+  Nn += n.G.dual_vertices.size();
 
-    std::vector<bool> vertex_in(n.G.vertices.size());
+  std::vector<bool> vertex_in(n.G.vertices.size());
 
-    for (unsigned i = 0; i < n.G.edges.size(); i++) {
-      if (!last_backbone[i]) {
-        vertex_in[n.G.edges[i].v[0]] = true;
-        vertex_in[n.G.edges[i].v[1]] = true;
-      }
+  for (unsigned i = 0; i < n.G.edges.size(); i++) {
+    if (!last_backbone[i]) {
+      vertex_in[n.G.edges[i].v[0]] = true;
+      vertex_in[n.G.edges[i].v[1]] = true;
     }
+  }
 
-    unsigned bb_size = 0;
+  unsigned bb_size = 0;
 
-    std::list<graph::coordinate> cb_co;
+  std::list<coordinate> cb_co;
 
-    for (unsigned i = 0; i < n.G.vertices.size(); i++) {
-      if (vertex_in[i]) {
-        bb_size++;
-        cb_co.push_back(n.G.vertices[i].r);
-      }
+  for (unsigned i = 0; i < n.G.vertices.size(); i++) {
+    if (vertex_in[i]) {
+      bb_size++;
+      cb_co.push_back(n.G.vertices[i].r);
     }
+  }
 
-    sb[bb_size]++;
-    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cb, cb_co, crack.first);
-    Nb += bb_size;
+  sb[bb_size]++;
+  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cb, cb_co, crack.first);
+  Nb += bb_size;
 
-    auto av_it = avalanches.rbegin();
-    av_it++;
+  auto av_it = avalanches.rbegin();
+  av_it++;
 
-    while (av_it != avalanches.rend()) {
-      sa[(*av_it).size() - 1]++;
-      Na += (*av_it).size();
-      Nq += (*av_it).size();
-      std::list<graph::coordinate> ca_co;
-      for (unsigned e : (*av_it)) {
-        ca_co.push_back(n.G.edges[e].r);
-      }
-      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ca, ca_co, crack.first);
-      autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cq, ca_co, {0, 1});
-      av_it++;
+  while (av_it != avalanches.rend()) {
+    sa[(*av_it).size() - 1]++;
+    Na += (*av_it).size();
+    Nq += (*av_it).size();
+    std::list<coordinate> ca_co;
+    for (unsigned e : (*av_it)) {
+      ca_co.push_back(n.G.edges[e].r);
     }
+    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ca, ca_co, crack.first);
+    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cq, ca_co, {0, 1});
+    av_it++;
+  }
 
-    sA[avalanches.back().size() - 1]++;
-    std::list<graph::coordinate> cA_co;
-    for (unsigned e : avalanches.back()) {
-      cA_co.push_back(n.G.edges[e].r);
-    }
-    autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cA, cA_co, crack.first);
-    NA += avalanches.back().size();
+  sA[avalanches.back().size() - 1]++;
+  std::list<coordinate> cA_co;
+  for (unsigned e : avalanches.back()) {
+    cA_co.push_back(n.G.edges[e].r);
+  }
+  autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cA, cA_co, crack.first);
+  NA += avalanches.back().size();
 
-    sd[num - 1]++;
+  sd[num - 1]++;
 
+  // Currents relative to tip of largest crack fragment
+  std::map<unsigned, std::set<unsigned>> fragments;
 
-    // Currents relative to tip of largest crack fragment
-    std::map<unsigned, std::list<unsigned>> fragments;
+  for (unsigned e : crack.second) {
+    if (last_backbone[e]) {
+      unsigned root = last_clusters.findroot(n.G.dual_edges[e].v[0]);
+      if (fragments.contains(root)) {
+        fragments[root].insert(e);
+      } else {
+        fragments[root] = {e};
+      }
+    }
+  }
 
-    for (unsigned e : crack.second) {
-      if (last_backbone[e]) {
-        unsigned root = last_clusters.findroot(n.G.dual_edges[e].v[0]);
-        if (fragments.contains(root)) {
-          fragments[root].push_back(e);
-        } else {
-          fragments[root] = {e};
+  if (fragments.size() > 0) {
+    std::set<unsigned> biggest_fragment =
+        std::max_element(fragments.begin(), fragments.end(), [](const auto& l1, const auto& l2) {
+          return l1.second.size() < l2.second.size();
+        })->second;
+
+    sp[biggest_fragment.size()]++;
+
+    // We have all the edges in the biggest fragment, in no meaningful order.
+    // To find its ends, we will start at a random place and traverse the
+    // fragment until we cannot continue. Then, we go back from where we
+    // started and do the same in the other direction. Along the way we will
+    // compute the (unwrapped) displacement between the two ends so that we
+    // will know which is to the left and which is to the right.
+
+    unsigned ce = *biggest_fragment.begin();
+    unsigned cv = n.G.dual_edges[ce].v[0];
+    coordinate Δr = n.G.Δr(ce); // direction is from v0 to v1
+
+    while (true) {
+      signed ne = -1;
+      unsigned nv;
+
+      for (unsigned e : n.G.dual_vertices[cv].ne) {
+        if (e != ce && biggest_fragment.contains(e)) {
+          ne = e;
+          unsigned v0 = n.G.dual_edges[ne].v[0];
+          unsigned v1 = n.G.dual_edges[ne].v[1];
+          Δr += v0 == cv ? n.G.Δr(ne) : -n.G.Δr(ne);
+          nv = v0 == cv ? v1 : v0;
+
+          break;
         }
       }
+
+      if (ne < 0) {
+        break;
+      } else {
+        ce = ne;
+        cv = nv;
+      }
     }
 
-    if (fragments.size() > 0) {
-      std::list<unsigned> biggest_fragment = std::max_element(fragments.begin(), fragments.end(), [] (const auto& l1, const auto& l2) {return l1.second.size() < l2.second.size();})->second;
+    unsigned end_1e = ce;
+    unsigned end_1v = cv;
 
-      sp[biggest_fragment.size()]++;
+    ce = *biggest_fragment.begin();
+    cv = n.G.dual_edges[ce].v[1];
 
-      std::map<unsigned, std::list<unsigned>> end_finder;
+    while (true) {
+      signed ne = -1;
+      unsigned nv;
 
-      for (unsigned e : biggest_fragment) {
-        unsigned v0 = n.G.dual_edges[e].v[0];
-        unsigned v1 = n.G.dual_edges[e].v[1];
-        if (end_finder.contains(v0)) {
-          end_finder[v0].push_back(v0);
-        } else {
-          end_finder[v0] = {v0};
-        }
-        if (end_finder.contains(v1)) {
-          end_finder[v1].push_back(v1);
-        } else {
-          end_finder[v1] = {v1};
+      for (unsigned e : n.G.dual_vertices[cv].ne) {
+        if (e != ce && biggest_fragment.contains(e)) {
+          ne = e;
+          unsigned v0 = n.G.dual_edges[ne].v[0];
+          unsigned v1 = n.G.dual_edges[ne].v[1];
+          Δr += v0 == cv ? -n.G.Δr(ne) : n.G.Δr(ne);
+          nv = v0 == cv ? v1 : v0;
+
+          break;
         }
       }
 
-      unsigned left_end, right_end;
-
-      bool comp;
+      if (ne < 0) {
+        break;
+      } else {
+        ce = ne;
+        cv = nv;
+      }
+    }
 
-      graph::coordinate Δr = {0,0};
-      unsigned e0v0 = n.G.dual_edges[biggest_fragment.front()].v[0];
-      unsigned e0v1 = n.G.dual_edges[biggest_fragment.front()].v[1];
-      unsigned e1v0 = n.G.dual_edges[*(biggest_fragment.begin()++)].v[0];
-      unsigned e1v1 = n.G.dual_edges[*(biggest_fragment.begin()++)].v[1];
+    unsigned end_0e = ce;
+    unsigned end_0v = cv;
 
-      unsigned last_vertex;
+    bool comp;
 
-      if (e0v0 == e1v0 || e0v0 == e1v1) {
-        last_vertex = e0v1;
-      } else {
-        last_vertex = e1v1;
-      }
+    if (crack.first[0] % 2 == 1) {
+      comp = Δr.x > 0;
+    } else {
+      comp = Δr.y > 0;
+    }
 
+    unsigned left_end, right_end;
 
-      for (unsigned e : biggest_fragment) {
-        unsigned v0 = n.G.dual_edges[e].v[0];
-        unsigned v1 = n.G.dual_edges[e].v[1];
+    if (comp) {
+      left_end = end_0e;
+      right_end = end_1e;
+    } else {
+      left_end = end_1e;
+      right_end = end_0e;
+    }
 
-        unsigned v_forward, v_back;
+    std::cout << crack.first[0] << " " << crack.first[1] << "\n";
+    std::cout << left_end << " " << right_end << "\n";
+    for (unsigned e : biggest_fragment) {
+      std::cout << e << " ";
+    }
+    std::cout << "\n";
+    std::cout << Δr.x << " " << Δr.y << "\n";
 
-        if (v0 == last_vertex) {
-          v_back = v0;
-          v_forward = v1;
-        } else {
-          v_back = v1;
-          v_forward = v0;
-        }
+    for (unsigned i = 0; i < n.G.dual_edges.size(); i++) {
+      if (!last_backbone[i]) {
+        double Δx_tmpl, Δy_tmpl, Δx_tmpr, Δy_tmpr, Δxr, Δxl, Δyr, Δyl;
 
-        double Δx = n.G.dual_vertices[v_forward].r.x - n.G.dual_vertices[v_back].r.x;
+        if (crack.first[0] % 2 == 1) {
+          Δx_tmpl = n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x;
+          Δxl = Δx_tmpl < 0 ? Δx_tmpl + n.G.L.x : Δx_tmpl;
 
-        if (n.G.dual_edges[e].crossings[0]) {
-          if (Δx > 0) {
-            Δx -= n.G.L.x;
-          } else {
-            Δx += n.G.L.x;
-          }
-        }
+          Δy_tmpl = fabs(n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y);
+          Δyl = Δy_tmpl > n.G.L.y / 2 ? n.G.L.y - Δy_tmpl : Δy_tmpl;
 
-        double Δy = n.G.dual_vertices[v_forward].r.y - n.G.dual_vertices[v_back].r.y;
+          Δx_tmpr = n.G.dual_edges[i].r.x - n.G.dual_edges[left_end].r.x;
+          Δxr = Δx_tmpr < 0 ? Δx_tmpr + n.G.L.x : Δx_tmpr;
 
-        if (n.G.dual_edges[e].crossings[1]) {
-          if (Δy > 0) {
-            Δy -= n.G.L.y;
-          } else {
-            Δy += n.G.L.y;
-          }
-        }
+          Δy_tmpr = fabs(n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y);
+          Δyr = Δy_tmpr > n.G.L.y / 2 ? n.G.L.y - Δy_tmpr : Δy_tmpr;
+        } else {
+          Δx_tmpl = n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y;
+          Δxl = Δx_tmpl < 0 ? Δx_tmpl + n.G.L.y : Δx_tmpl;
 
-        Δr.x += Δx;
-        Δr.y += Δy;
-      }
+          Δy_tmpl = fabs(n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x);
+          Δyl = Δy_tmpl > n.G.L.x / 2 ? n.G.L.x - Δy_tmpl : Δy_tmpl;
 
-      if (crack.first[0] % 2 == 1) {
-        comp = Δr.x > 0;
-      } else {
-        comp = Δr.y > 0;
-      }
+          Δx_tmpr = n.G.dual_edges[i].r.y - n.G.dual_edges[left_end].r.y;
+          Δxr = Δx_tmpr < 0 ? Δx_tmpr + n.G.L.y : Δx_tmpr;
 
-      if (comp) {
-        left_end = biggest_fragment.back();
-        right_end = biggest_fragment.front();
-      } else {
-        left_end = biggest_fragment.front();
-        right_end = biggest_fragment.back();
-      }
+          Δy_tmpr = fabs(n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x);
+          Δyr = Δy_tmpr > n.G.L.x / 2 ? n.G.L.x - Δy_tmpr : Δy_tmpr;
+        }
 
-      std::cout << crack.first[0] << " " << crack.first[1] << "\n";
-      std::cout << left_end << " " << right_end << "\n";
-      for (unsigned e : biggest_fragment) {
-        std::cout << e << " ";
-      }
-        std::cout << "\n";
-
-      for (unsigned i = 0; i < n.G.dual_edges.size(); i++) {
-        if (!last_backbone[i]) {
-          double Δx_tmpl, Δy_tmpl, Δx_tmpr, Δy_tmpr, Δxr, Δxl, Δyr, Δyl;
-
-          if (crack.first[0] % 2 == 1) {
-            Δx_tmpl = n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x;
-            Δxl = Δx_tmpl < 0 ? Δx_tmpl + n.G.L.x : Δx_tmpl;
-
-            Δy_tmpl = fabs(n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y);
-            Δyl = Δy_tmpl > n.G.L.y / 2 ? n.G.L.y - Δy_tmpl : Δy_tmpl;
-
-            Δx_tmpr = n.G.dual_edges[i].r.x - n.G.dual_edges[left_end].r.x;
-            Δxr = Δx_tmpr < 0 ? Δx_tmpr + n.G.L.x : Δx_tmpr;
-
-            Δy_tmpr = fabs(n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y);
-            Δyr = Δy_tmpr > n.G.L.y / 2 ? n.G.L.y - Δy_tmpr : Δy_tmpr;
-          } else {
-            Δx_tmpl = n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y;
-            Δxl = Δx_tmpl < 0 ? Δx_tmpl + n.G.L.y : Δx_tmpl;
-
-            Δy_tmpl = fabs(n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x);
-            Δyl = Δy_tmpl > n.G.L.x / 2 ? n.G.L.x - Δy_tmpl : Δy_tmpl;
-
-            Δx_tmpr = n.G.dual_edges[i].r.y - n.G.dual_edges[left_end].r.y;
-            Δxr = Δx_tmpr < 0 ? Δx_tmpr + n.G.L.y : Δx_tmpr;
-
-            Δy_tmpr = fabs(n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x);
-            Δyr = Δy_tmpr > n.G.L.x / 2 ? n.G.L.x - Δy_tmpr : Δy_tmpr;
-          }
-
-          if (Δxl < Δxr) {
-            cS[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δyl / n.G.L.y)) +
-               (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δxl / n.G.L.x))] +=
-                last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
-          } else {
-            cS[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δyr / n.G.L.y)) +
-               (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δxr / n.G.L.x))] +=
-                last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
-          }
+        if (Δxl < Δxr) {
+          cS[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δyl / n.G.L.y)) +
+             (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δxl / n.G.L.x))] +=
+              last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
+        } else {
+          cS[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δyr / n.G.L.y)) +
+             (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δxr / n.G.L.x))] +=
+              last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
         }
       }
-
-      NS++;
-
-    } else {
-      sp[0]++;
     }
 
-    aG.push_back(last_currents.conductivity);
+    NS++;
 
+  } else {
+    sp[0]++;
   }
 
+  aG.push_back(last_currents.conductivity);
+}
+