From 0a79391dc49a9662100969e53b55fef4d367df77 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 12 Nov 2019 11:38:29 -0500 Subject: more progress on crack stress --- lib/include/graph.hpp | 43 +++++- lib/src/graph.cpp | 30 +++- lib/src/problem.cpp | 12 +- src/measurements.cpp | 415 +++++++++++++++++++++++++------------------------- 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 + 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 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& 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& data, } template -void update_field_file(std::string id, uint64_t N, const std::vector& data, std::string model_string) { +void update_field_file(std::string id, uint64_t N, const std::vector& 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& data, - const std::list& pos, std::array count) { - for (std::list::const_iterator it1 = pos.begin(); it1 != pos.end(); it1++) { - for (std::list::const_iterator it2 = it1; it2 != pos.end(); it2++) { + const std::list& pos, std::array count) { + for (std::list::const_iterator it1 = pos.begin(); it1 != pos.end(); it1++) { + for (std::list::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& data, - const std::list& pos1, - const std::list& pos2) { - for (std::list::const_iterator it1 = pos1.begin(); it1 != pos1.end(); it1++) { - for (std::list::const_iterator it2 = pos2.begin(); it2 != pos2.end(); + const std::list& pos1, + const std::list& pos2) { + for (std::list::const_iterator it1 = pos1.begin(); it1 != pos1.end(); it1++) { + for (std::list::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 cl_cs; + std::list 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> crit_components(n.G.dual_vertices.size()); + std::vector> 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> components(n.G.dual_vertices.size()); + std::vector> 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 vertex_in(n.G.vertices.size()); + std::vector 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 cb_co; + std::list 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 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 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 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 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> fragments; - // Currents relative to tip of largest crack fragment - std::map> 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 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 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> 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); +} + -- cgit v1.2.3-54-g00ecf