From 02dafc73e73a56d0d857a7dc4a4471fdf7abc2d9 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 12 Nov 2019 17:50:36 -0500 Subject: fixed the crack stress, adding new measurements for mean ordinary stress --- src/measurements.cpp | 84 ++++++++++++++++++++++++++++++++++------------------ src/measurements.hpp | 4 ++- 2 files changed, 59 insertions(+), 29 deletions(-) diff --git a/src/measurements.cpp b/src/measurements.cpp index e7aa528..446d29c 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -171,7 +171,7 @@ ma::ma(unsigned n, double a, double beta, double weight) 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) { + cS(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 2), NS(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 2), cT(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 2), NT(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 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) + "_"; @@ -190,7 +190,6 @@ ma::ma(unsigned n, double a, double beta, double weight) Na = 0; NA = 0; Nq = 0; - NS = 0; } ma::ma(unsigned Lx, unsigned Ly, double beta, double weight) @@ -200,7 +199,7 @@ ma::ma(unsigned Lx, unsigned Ly, double beta, double weight) 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)), cS((Lx / 2 + 1) * (Ly / 2 + 1)), last_clusters(Lx * Ly / 2) { + cq((Lx / 2 + 1) * (Ly / 2 + 1)), cS(Lx * Ly / 2), NS(Lx * Ly / 2), cT(Lx * Ly / 2), NT(Lx * Ly / 2), 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) + "_"; @@ -219,7 +218,6 @@ ma::ma(unsigned Lx, unsigned Ly, double beta, double weight) Na = 0; NA = 0; Nq = 0; - NS = 0; } ma::~ma() { @@ -244,7 +242,10 @@ 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); + update_field_file("cS", 0, cS, model_string); + update_field_file("NS", 0, NS, model_string); + update_field_file("cT", 0, cT, model_string); + update_field_file("NT", 0, NT, model_string); } void ma::pre_fracture(const network&) { @@ -485,14 +486,6 @@ void ma::post_fracture(network& n) { right_end = end_0e; } - 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"; - 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; @@ -504,10 +497,10 @@ void ma::post_fracture(network& n) { Δ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; + Δx_tmpr = n.G.dual_edges[i].r.x - n.G.dual_edges[right_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); + Δy_tmpr = fabs(n.G.dual_edges[right_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; @@ -516,31 +509,66 @@ void ma::post_fracture(network& n) { Δ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; + Δx_tmpr = n.G.dual_edges[i].r.y - n.G.dual_edges[right_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); + Δy_tmpr = fabs(n.G.dual_edges[right_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]; - } + cS[(unsigned)floor(Mx * (Δyl / n.G.L.y)) + + (Mx / 2) * (unsigned)floor(My * (Δxl / n.G.L.x))] += + last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2]; + NS[(unsigned)floor(Mx * (Δyl / n.G.L.y)) + + (Mx / 2) * (unsigned)floor(My * (Δxl / n.G.L.x))]++; + NS[(unsigned)floor(Mx * (Δyr / n.G.L.y)) + + (Mx / 2) * (unsigned)floor(My * (Δxr / n.G.L.x))]++; + cS[(unsigned)floor(Mx * (Δyr / n.G.L.y)) + + (Mx / 2) * (unsigned)floor(My * (Δxr / n.G.L.x))] += + last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2]; } } - NS++; - } else { sp[0]++; } + for (unsigned i = 0; i < n.G.dual_edges.size(); i++) { + if (!last_backbone[i]) { + for (unsigned j = 0; j < n.G.dual_edges.size(); j++) { + if (!last_backbone[j]) { + double Δx_tmpl, Δy_tmpl, Δx_tmpr, Δy_tmpr, Δxr, Δxl, Δyr, Δyl; + + if (crack.first[0] % 2 == 1) { + Δx_tmpl = fabs(n.G.dual_edges[j].r.x - n.G.dual_edges[i].r.x); + Δxl = Δx_tmpl > n.G.L.x / 2 ? n.G.L.x - Δx_tmpl : Δx_tmpl; + + Δy_tmpl = fabs(n.G.dual_edges[j].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; + } else { + Δx_tmpl = fabs(n.G.dual_edges[j].r.y - n.G.dual_edges[i].r.y); + Δxl = Δx_tmpl > n.G.L.y / 2 ? n.G.L.y - Δx_tmpl : Δx_tmpl; + + Δy_tmpl = fabs(n.G.dual_edges[j].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; + + } + + cS[(unsigned)floor(Mx * (Δyl / n.G.L.y)) + + (Mx / 2) * (unsigned)floor(My * (Δxl / n.G.L.x))] += + last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2]; + NS[(unsigned)floor(Mx * (Δyl / n.G.L.y)) + + (Mx / 2) * (unsigned)floor(My * (Δxl / n.G.L.x))]++; + NS[(unsigned)floor(Mx * (Δyr / n.G.L.y)) + + (Mx / 2) * (unsigned)floor(My * (Δxr / n.G.L.x))]++; + cS[(unsigned)floor(Mx * (Δyr / n.G.L.y)) + + (Mx / 2) * (unsigned)floor(My * (Δxr / n.G.L.x))] += + last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2]; + } + } + } + } + aG.push_back(last_currents.conductivity); } diff --git a/src/measurements.hpp b/src/measurements.hpp index c0458b5..0786cf8 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -29,7 +29,6 @@ class ma : public hooks { uint64_t Na; uint64_t NA; uint64_t Nq; - uint64_t NS; // measurement storage std::vector sc; @@ -57,6 +56,9 @@ class ma : public hooks { std::list> aG; std::vector cS; + std::vector NS; + std::vector cT; + std::vector NT; std::vector sp; ClusterTree last_clusters; -- cgit v1.2.3-54-g00ecf