From ffa995abf6d7d4855818c4df23ec580d36460457 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sat, 9 Nov 2019 10:45:44 -0500 Subject: started implementing crack stress measurements, but probably doesn't work --- src/measurements.cpp | 156 +++++++++++++++++++++++++++++++++++++++++++++++++-- src/measurements.hpp | 7 +++ 2 files changed, 159 insertions(+), 4 deletions(-) diff --git a/src/measurements.cpp b/src/measurements.cpp index 2432198..2e45168 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -2,6 +2,7 @@ #include "measurements.hpp" #include #include +#include class badcycleException : public std::exception { virtual const char* what() const throw() { @@ -162,10 +163,10 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u 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), si(10000), sI(10000), cc(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cn(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), + 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)), last_clusters(2 * n) { + 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) + "_"; @@ -188,12 +189,12 @@ ma::ma(unsigned n, double a, double beta, double weight) ma::ma(unsigned Lx, unsigned Ly, double beta, double weight) : sc(Lx * Ly / 2), sn(Lx * Ly / 2), ss(Lx * Ly / 2), sm(Lx * Ly / 2), sl(Lx * Ly / 2), - sb(Lx * Ly / 2 + 1), sd(Lx * Ly), sa(Lx * Ly), sA(Lx * Ly), si(10000), sI(10000), + sb(Lx * Ly / 2 + 1), sd(Lx * Ly), sa(Lx * Ly), sA(Lx * Ly), sp(Lx * Ly), si(10000), sI(10000), cc((Lx / 2 + 1) * (Ly / 2 + 1)), cn((Lx / 2 + 1) * (Ly / 2 + 1)), cs((Lx / 2 + 1) * (Ly / 2 + 1)), cm((Lx / 2 + 1) * (Ly / 2 + 1)), cl((Lx / 2 + 1) * (Ly / 2 + 1)), cb((Lx / 2 + 1) * (Ly / 2 + 1)), ca((Lx / 2 + 1) * (Ly / 2 + 1)), cA((Lx / 2 + 1) * (Ly / 2 + 1)), - cq((Lx / 2 + 1) * (Ly / 2 + 1)), last_clusters(Lx * Ly / 2) { + cq((Lx / 2 + 1) * (Ly / 2 + 1)), cS((Lx / 2 + 1) * (Ly / 2 + 1)), last_clusters(Lx * Ly / 2) { if (beta != 0.0) { model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + std::to_string(weight) + "_"; @@ -212,6 +213,7 @@ ma::ma(unsigned Lx, unsigned Ly, double beta, double weight) Na = 0; NA = 0; Nq = 0; + NS = 0; } ma::~ma() { @@ -224,6 +226,7 @@ ma::~ma() { update_distribution_file("sd", sd, model_string); update_distribution_file("sa", sa, model_string); update_distribution_file("sA", sA, model_string); + update_distribution_file("sp", sp, model_string); update_distribution_file("si", si, model_string); update_distribution_file("sI", sI, model_string); update_field_file("cc", Nc, cc, model_string); @@ -235,6 +238,7 @@ ma::~ma() { update_field_file("ca", Na, ca, model_string); update_field_file("cA", NA, cA, model_string); update_field_file("cq", Nq, cq, model_string); + update_field_file("cS", NS, cS, model_string); } void ma::pre_fracture(const network&) { @@ -247,6 +251,8 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) { long double c = net.thresholds[i] - logl(cur.currents[i]); if (c > lv) { last_clusters = net.C; + last_currents = cur; + last_backbone = net.backbone; lv = c; avalanches.push_back({i}); } else { @@ -364,5 +370,147 @@ void ma::post_fracture(network& n) { NA += avalanches.back().size(); sd[num - 1]++; + + + // Currents relative to tip of largest crack fragment + std::map> fragments; + + for (unsigned e : crack.second) { + 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::list 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()]++; + + unsigned left_end, right_end; + + bool comp; + + 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 last_vertex; + + if (e0v0 == e1v0 || e0v0 == e1v1) { + last_vertex = e0v1; + } else { + last_vertex = e1v1; + } + + + for (unsigned e : biggest_fragment) { + unsigned v0 = n.G.dual_edges[e].v[0]; + unsigned v1 = n.G.dual_edges[e].v[1]; + + unsigned v_forward, v_back; + + if (v0 == last_vertex) { + v_back = v0; + v_forward = v1; + } else { + v_back = v1; + v_forward = v0; + } + + double Δx = n.G.dual_vertices[v_forward].r.x - n.G.dual_vertices[v_back].r.x; + + if (n.G.dual_edges[e].crossings[0]) { + if (Δx > 0) { + Δx -= n.G.L.x; + } else { + Δx += n.G.L.x; + } + } + + double Δy = n.G.dual_vertices[v_forward].r.y - n.G.dual_vertices[v_back].r.y; + + if (n.G.dual_edges[e].crossings[1]) { + if (Δy > 0) { + Δy -= n.G.L.y; + } else { + Δy += n.G.L.y; + } + } + + Δr.x += Δx; + Δr.y += Δy; + } + + if (crack.first[0] % 2 == 1) { + comp = Δr.x > 0; + } else { + comp = Δr.y > 0; + } + + if (comp) { + left_end = biggest_fragment.back(); + right_end = biggest_fragment.front(); + } else { + left_end = biggest_fragment.front(); + right_end = biggest_fragment.back(); + } + + std::cout << n.G.dual_edges[left_end].r.x << " " << n.G.dual_edges[left_end].r.y << " " << n.G.dual_edges[right_end].r.x << " " << n.G.dual_edges[right_end].r.y << " "<< crack.first[0] % 2 << " " << Δr.x << " " << Δr.y << " " << biggest_fragment.size() << "\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 ? Δy_tmpl - n.G.L.y / 2 : Δ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 ? Δy_tmpr - n.G.L.y / 2 : Δ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 ? Δy_tmpl - n.G.L.x / 2 : Δ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 ? Δy_tmpr - n.G.L.x / 2 : Δ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]; + } + } + } + + NS++; + + } else { + sp[0]++; + } + + aG.push_back(last_currents.conductivity); + } diff --git a/src/measurements.hpp b/src/measurements.hpp index ebf79e5..c0458b5 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -29,6 +29,7 @@ class ma : public hooks { uint64_t Na; uint64_t NA; uint64_t Nq; + uint64_t NS; // measurement storage std::vector sc; @@ -54,7 +55,13 @@ class ma : public hooks { std::vector cA; std::vector cq; + std::list> aG; + std::vector cS; + std::vector sp; + ClusterTree last_clusters; + current_info last_currents; + std::vector last_backbone; public: long double lv; -- cgit v1.2.3-54-g00ecf