#include "measurements.hpp" #include #include class badcycleException : public std::exception { virtual const char* what() const throw() { return "Could not find a valid cycle on the broken system."; } } badcycleex; void update_distribution_file(std::string id, const std::vector& data, std::string model_string) { std::string filename = model_string + id + ".dat"; std::ifstream file(filename); std::vector data_old(data.size(), 0); if (file.is_open()) { for (unsigned i = 0; i < data.size(); i++) { uint64_t num; file >> num; data_old[i] = num; } file.close(); } std::ofstream file_out(filename); for (unsigned i = 0; i < data.size(); i++) { file_out << std::fixed << data_old[i] + data[i] << " "; } file_out.close(); } template 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); std::vector data_old(data.size(), 0); uint64_t old_N = 0; if (file.is_open()) { file >> old_N; for (unsigned j = 0; j < data.size(); j++) { file >> data_old[j]; } file.close(); } std::ofstream file_out(filename); file_out << old_N + N << "\n"; for (unsigned j = 0; j < data.size(); j++) { file_out << data_old[j] + data[j] << " "; } file_out.close(); } /* fftw_complex* data_transform(unsigned Mx, unsigned My, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) { fftw_execute(forward_plan); fftw_complex* output = (fftw_complex*)malloc(Mx * (My / 2 + 1) * sizeof(fftw_complex)); for (unsigned i = 0; i < Mx * (My / 2 + 1); i++) { output[i][0] = fftw_forward_out[i][0]; output[i][1] = fftw_forward_out[i][1]; } return output; } */ /* template void correlation(unsigned Mx, unsigned My, std::vector>& data, const fftw_complex* tx1, const fftw_complex* tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) { for (unsigned i = 0; i < Mx * (My / 2 + 1); i++) { fftw_reverse_in[i][0] = tx1[i][0] * tx2[i][0] + tx1[i][1] * tx2[i][1]; fftw_reverse_in[i][1] = tx1[i][0] * tx2[i][1] - tx1[i][1] * tx2[i][0]; } fftw_execute(reverse_plan); for (unsigned j = 0; j < data.size(); j++) { for (unsigned i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) { data[j][i] += (T)pow(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My), j + 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++) { double Δx_tmp = fabs(it1->x - it2->x); double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp; double Δy_tmp = fabs(it1->y - it2->y); 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))]++; } else { 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(); it2++) { double Δx_tmp = fabs(it1->x - it2->x); double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp; double Δy_tmp = fabs(it1->y - it2->y); double Δy = Δy_tmp < Ly / 2 ? Δy_tmp : Ly - Δy_tmp; data[floor((Mx * Δx) / Lx) + (Mx / 2) * floor((My * Δy) / Ly)]++; } } } /* template void autocorrelation(unsigned Mx, unsigned My, std::vector>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) { fftw_execute(forward_plan); for (unsigned i = 0; i < My * (Mx / 2 + 1); i++) { fftw_reverse_in[i][0] = pow(fftw_forward_out[i][0], 2) + pow(fftw_forward_out[i][1], 2); fftw_reverse_in[i][1] = 0.0; } fftw_execute(reverse_plan); for (unsigned j = 0; j < out_data.size(); j++) { for (unsigned i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) { out_data[j][i] += (T)pow(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My), j + 1); } } } */ unsigned edge_r_to_ind(graph::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), sC(2 * n), sn(2 * n), ss(2 * n), sm(2 * n), sl(2 * n), sb(n + 1), sd(3 * n), sD(3 * n), sf(3 * n), sa(3 * n), sA(3 * n), se(3 * n), sE(3 * n), si(10000), sI(10000), sy(10000), sY(10000), sz(10000), sZ(10000), cc(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), 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)), ce(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cE(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cq(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), last_clusters(2 * n), last_clusters_v(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) + "_"; } else { model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_INF_" + std::to_string(weight) + "_"; } Mx = 4 * ceil(sqrt(n)) * a; My = 4 * ceil(sqrt(n)) / a; Nc = 0; Nl = 0; Nm = 0; Ns = 0; Nn = 0; Nb = 0; Na = 0; NA = 0; Nq = 0; Ne = 0; NE = 0; NC = 0; } ma::ma(unsigned Lx, unsigned Ly, double beta, double weight) : sc(Lx * Ly / 2), 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), sD(Lx * Ly), sf(Lx * Ly), sa(Lx * Ly), sA(Lx * Ly), se(Lx * Ly), sE(Lx * Ly), si(10000), sI(10000), sy(10000), sY(10000), sz(10000), sZ(10000), cc((Lx / 2 + 1) * (Ly / 2 + 1)), 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)), ce((Lx / 2 + 1) * (Ly / 2 + 1)), cE((Lx / 2 + 1) * (Ly / 2 + 1)), cq((Lx / 2 + 1) * (Ly / 2 + 1)), last_clusters(Lx * Ly / 2), last_clusters_v(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) + "_"; } else { model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_INF_" + std::to_string(weight) + "_"; } Mx = Lx; My = Ly; Nc = 0; NC = 0; Nl = 0; Nm = 0; Ns = 0; Nn = 0; Nb = 0; Na = 0; NA = 0; Ne = 0; NE = 0; Nq = 0; } ma::~ma() { update_distribution_file("sc", sc, model_string); update_distribution_file("sC", sC, model_string); update_distribution_file("sn", sn, model_string); update_distribution_file("ss", ss, model_string); update_distribution_file("sm", sm, model_string); update_distribution_file("sl", sl, model_string); update_distribution_file("sb", sb, model_string); update_distribution_file("sd", sd, model_string); update_distribution_file("sD", sD, model_string); update_distribution_file("sf", sf, model_string); update_distribution_file("sa", sa, model_string); update_distribution_file("sA", sA, model_string); update_distribution_file("se", se, model_string); update_distribution_file("sE", sE, model_string); update_distribution_file("si", si, model_string); update_distribution_file("sI", sI, model_string); update_distribution_file("sy", sy, model_string); update_distribution_file("sY", sY, model_string); update_distribution_file("sz", sz, model_string); update_distribution_file("sZ", sZ, model_string); update_field_file("cc", Nc, cc, model_string); update_field_file("cC", NC, cC, model_string); update_field_file("cl", Nl, cl, model_string); update_field_file("cm", Nm, cm, model_string); update_field_file("cs", Ns, cs, model_string); update_field_file("cn", Nn, cn, model_string); update_field_file("cb", Nb, cb, model_string); update_field_file("ca", Na, ca, model_string); update_field_file("cA", NA, cA, model_string); update_field_file("ce", Ne, ce, model_string); update_field_file("cE", NE, cE, model_string); update_field_file("cq", Nq, cq, model_string); } void ma::pre_fracture(const network&) { lc = std::numeric_limits::lowest(); lv = std::numeric_limits::lowest(); moduliA = {}; moduliE = {}; avalanches = {}; evalanches = {}; num = 0; } void ma::bond_broken(const network& net, const current_info& cur, unsigned i) { long double c = net.thresholds[i] - logl(cur.currents[i]); long double v = net.thresholds[i] - logl(cur.currents[i] * cur.conductivity[0]); if (c > lc) { last_clusters = net.C; lc = c; moduliA.push_back(cur.conductivity[0]); avalanches.push_back({i}); } else { avalanches.back().push_back(i); } if (v > lv) { last_clusters_v = net.C; lv = v; moduliE.push_back(cur.conductivity[0]); evalanches.push_back({i}); } else { evalanches.back().push_back(i); } for (unsigned j = 0; j < cur.currents.size(); j++) { if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0 && (!net.backbone[j] || j == i)) { si[(unsigned)(10000 * (logl(cur.currents[j]) + 100) / 100)]++; } else if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0 && net.backbone[j] && j != i) { sI[(unsigned)(10000 * (logl(cur.currents[j]) + 100) / 100)]++; } } num++; } void ma::post_fracture(network& n) { auto crack = find_minimal_crack(n, avalanches.back().back()); std::list cl_cs; sl[crack.second.size() - 1]++; for (unsigned e : crack.second) { cl_cs.push_back(n.G.dual_edges[e].r); } 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()); 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); } for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) { if (crit_components[i].size() > 0) { sc[crit_components[i].size() - 1]++; } autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cc, crit_components[i], crack.first); } Nc += n.G.dual_vertices.size(); std::vector> crit_components_v(n.G.dual_vertices.size()); for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) { crit_components_v[last_clusters_v.findroot(i)].push_back(n.G.dual_vertices[i].r); } for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) { if (crit_components_v[i].size() > 0) { sC[crit_components_v[i].size() - 1]++; } autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cC, crit_components_v[i], crack.first); } NC += 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); } 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); } } Nn += n.G.dual_vertices.size(); std::vector vertex_in(n.G.vertices.size()); for (unsigned i = 0; i < n.G.edges.size(); i++) { if (!n.backbone[i]) { vertex_in[n.G.edges[i].v[0]] = true; vertex_in[n.G.edges[i].v[1]] = true; } } unsigned bb_size = 0; 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); } } 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++; 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++; } double mAFinal = log(moduliA.back()); if (mAFinal < 2 && mAFinal > -18) { sZ[(unsigned)(10000 * (mAFinal + 18) / 20)]++; } auto am_it_curr = moduliA.rbegin(); auto am_it_next = std::next(am_it_curr); while (am_it_next != moduliA.rend()) { double mm = log(*am_it_next - *am_it_curr); if (mm < 2 && mm > -18) { sz[(unsigned)(10000 * (mm + 18) / 20)]++; } am_it_curr++; am_it_next++; } double mEFinal = log(moduliE.back()); if (mEFinal < 2 && mEFinal > -18) { sY[(unsigned)(10000 * (mEFinal + 18) / 20)]++; } auto em_it_curr = moduliE.rbegin(); auto em_it_next = std::next(em_it_curr); while (em_it_next != moduliE.rend()) { double mm = log(*em_it_next - *em_it_curr); if (mm < 2 && mm > -18) { sy[(unsigned)(10000 * (mm + 18) / 20)]++; } em_it_curr++; em_it_next++; } auto ev_it = evalanches.rbegin(); ev_it++; while (ev_it != evalanches.rend()) { se[(*ev_it).size() - 1]++; Ne += (*ev_it).size(); std::list ce_co; for (unsigned e : (*ev_it)) { ce_co.push_back(n.G.edges[e].r); } autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ce, ce_co, crack.first); ev_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(); sE[evalanches.back().size() - 1]++; std::list cE_co; for (unsigned e : evalanches.back()) { cE_co.push_back(n.G.edges[e].r); } autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cE, cE_co, crack.first); NE += evalanches.back().size(); sd[num - 1]++; sD[num - avalanches.back().size()]++; sf[num - evalanches.back().size()]++; }