diff options
Diffstat (limited to 'src/measurements.cpp')
-rw-r--r-- | src/measurements.cpp | 303 |
1 files changed, 87 insertions, 216 deletions
diff --git a/src/measurements.cpp b/src/measurements.cpp index c8cc73c..43483be 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -2,128 +2,6 @@ #include "measurements.hpp" #include <iostream> -template <class T> -bool is_shorter(std::list<T> l1, std::list<T> l2) { - return l1.size() < l2.size(); -} - -bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) { - return true; -} - -std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) { - Graph Gtmp(n.G.vertices.size()); - std::list<unsigned int> removed_edges; - - class add_tree_edges : public boost::default_dfs_visitor { - public: - Graph& G; - std::list<unsigned int>& E; - - add_tree_edges(Graph& G, std::list<unsigned int>& E) : G(G), E(E) {} - - void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { - boost::add_edge(boost::source(e, g), boost::target(e, g), g[e], G); - } - - void back_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { - if (!(boost::edge(boost::source(e, g), boost::target(e, g), G).second)) { - E.push_back(g[e].index); - } - } - }; - - add_tree_edges ate(Gtmp, removed_edges); - boost::depth_first_search(G, visitor(ate)); - - class find_cycle : public boost::default_dfs_visitor { - public: - std::list<unsigned int>& E; - unsigned int end; - struct done{}; - - find_cycle(std::list<unsigned int>& E, unsigned int end) : E(E), end(end) {} - - void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) { - if (v == end) { - throw done{}; - } - } - - void examine_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { - E.push_back(g[e].index); - } - - void finish_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { - E.erase(std::find(E.begin(), E.end(), g[e].index)); - } - }; - - std::list<std::list<unsigned int>> cycles; - - for (auto edge : removed_edges) { - std::list<unsigned int> cycle = {edge}; - find_cycle vis(cycle, n.G.dual_edges[edge][1]); - std::vector<boost::default_color_type> new_color_map(boost::num_vertices(Gtmp)); - try { - boost::depth_first_visit(Gtmp, n.G.dual_edges[edge][0], vis, boost::make_iterator_property_map(new_color_map.begin(), boost::get(boost::vertex_index, Gtmp), new_color_map[0])); - } catch(find_cycle::done const&) { - cycles.push_back(cycle); - } - } - - if (cycles.size() > 1) { - std::list<std::valarray<uint8_t>> bool_cycles; - for (auto cycle : cycles) { - std::valarray<uint8_t> bool_cycle(n.G.edges.size()); - for (auto v : cycle) { - bool_cycle[v] = 1; - } - bool_cycles.push_back(bool_cycle); - } - - // generate all possible cycles by taking xor of the edge sets of the known cycles - for (auto it1 = bool_cycles.begin(); it1 != std::prev(bool_cycles.end()); it1++) { - for (auto it2 = std::next(it1); it2 != bool_cycles.end(); it2++) { - std::valarray<uint8_t> new_bool_cycle = (*it1) ^ (*it2); - std::list<unsigned int> new_cycle; - unsigned int pos = 0; - for (uint8_t included : new_bool_cycle) { - if (included) { - new_cycle.push_back(pos); - } - pos++; - } - cycles.push_back(new_cycle); - } - } - - // find the cycle representing the crack by counting boundary crossings - for (auto cycle : cycles) { - std::array<unsigned int, 2> crossing_count{0,0}; - - for (auto edge : cycle) { - double dx = fabs(n.G.dual_vertices[n.G.dual_edges[edge][0]].r.x - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.x); - if (dx > n.G.L.x / 2) { - crossing_count[0]++; - } - double dy = fabs(n.G.dual_vertices[n.G.dual_edges[edge][0]].r.y - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.y); - if (dy > n.G.L.y / 2) { - crossing_count[1]++; - } - } - - if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) { - return cycle; - } - } - } else { - return cycles.front(); - } - - exit(5); -} - void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta) { std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat"; std::ifstream file(filename); @@ -153,8 +31,8 @@ void update_distribution_file(std::string id, const std::vector<uint64_t>& data, } template <class T> -void update_field_file(std::string id, const std::vector<T>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta) { - std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat"; +void update_field_file(std::string id, const std::vector<T>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta, unsigned int Mx, unsigned int My) { + std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat"; std::ifstream file(filename); uint64_t N_old = 0; @@ -226,18 +104,22 @@ void autocorrelation(unsigned int Lx, unsigned int Ly, std::vector<T>& out_data, } } -ma::ma(unsigned int Lx, unsigned int Ly, double beta) : Lx(Lx), Ly(Ly), beta(beta), G(Lx * Ly / 2), +unsigned int edge_r_to_ind(graph::coordinate r, unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My) { + return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly); +} + +ma::ma(unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My, double beta) : Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(Lx * Ly), sc(Lx * Ly, 0), sa(Lx * Ly, 0), sd(Lx * Ly, 0), - sb(log2(Ly) + 1, 0), - Ccc((Lx / 2 + 1) * (Ly / 2 + 1), 0), - Css((Lx / 2 + 1) * (Ly / 2 + 1), 0), - Cmm((Lx / 2 + 1) * (Ly / 2 + 1), 0), - Caa((Lx / 2 + 1) * (Ly / 2 + 1), 0), - Cdd((Lx / 2 + 1) * (Ly / 2 + 1), 0), - Cll((Lx / 2 + 1) * (Ly / 2 + 1), 0), - Cee((Lx / 2 + 1) * (Ly / 2 + 1), 0) + sb(log2(Mx < My ? Mx : My) + 1, 0), + Ccc((Mx / 2 + 1) * (My / 2 + 1), 0), + Css((Mx / 2 + 1) * (My / 2 + 1), 0), + Cmm((Mx / 2 + 1) * (My / 2 + 1), 0), + Caa((Mx / 2 + 1) * (My / 2 + 1), 0), + Cdd((Mx / 2 + 1) * (My / 2 + 1), 0), + Cll((Mx / 2 + 1) * (My / 2 + 1), 0), + Cee((Mx / 2 + 1) * (My / 2 + 1), 0) { N = 0; Nc = 0; @@ -246,17 +128,13 @@ ma::ma(unsigned int Lx, unsigned int Ly, double beta) : Lx(Lx), Ly(Ly), beta(bet // FFTW setup for correlation functions fftw_set_timelimit(1); - fftw_forward_in = (double *)fftw_malloc(Lx * Ly * sizeof(double)); - fftw_forward_out = (fftw_complex *)fftw_malloc(Lx * Ly * sizeof(fftw_complex)); - fftw_reverse_in = (fftw_complex *)fftw_malloc(Lx * Ly * sizeof(fftw_complex)); - fftw_reverse_out = (double *)fftw_malloc(Lx * Ly * sizeof(double)); - - forward_plan = fftw_plan_dft_r2c_2d(Ly, Lx, fftw_forward_in, fftw_forward_out, 0); - reverse_plan = fftw_plan_dft_c2r_2d(Ly, Lx, fftw_reverse_in, fftw_reverse_out, 0); - - std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_broken.dat"; + fftw_forward_in = (double *)fftw_malloc(Mx * My * sizeof(double)); + fftw_forward_out = (fftw_complex *)fftw_malloc(Mx * My * sizeof(fftw_complex)); + fftw_reverse_in = (fftw_complex *)fftw_malloc(Mx * My * sizeof(fftw_complex)); + fftw_reverse_out = (double *)fftw_malloc(Mx * My * sizeof(double)); - bondfile.open(filename, std::ifstream::app); + forward_plan = fftw_plan_dft_r2c_2d(My, Mx, fftw_forward_in, fftw_forward_out, 0); + reverse_plan = fftw_plan_dft_c2r_2d(My, Mx, fftw_reverse_in, fftw_reverse_out, 0); } ma::~ma() { @@ -274,18 +152,16 @@ ma::~ma() { update_distribution_file("sd", sd, N, Lx, Ly, beta); update_distribution_file("sb", sb, N, Lx, Ly, beta); - update_field_file("Ccc", Ccc, Nc, Lx, Ly, beta); - update_field_file("Css", Css, N, Lx, Ly, beta); - update_field_file("Cmm", Cmm, N, Lx, Ly, beta); - update_field_file("Cdd", Cdd, N, Lx, Ly, beta); - update_field_file("Caa", Caa, Na, Lx, Ly, beta); - update_field_file("Cll", Cll, N, Lx, Ly, beta); - update_field_file("Cee", Cee, N, Lx, Ly, beta); - - bondfile.close(); + update_field_file("Ccc", Ccc, Nc, Lx, Ly, beta, Mx, My); + update_field_file("Css", Css, N, Lx, Ly, beta, Mx, My); + update_field_file("Cmm", Cmm, N, Lx, Ly, beta, Mx, My); + update_field_file("Cdd", Cdd, N, Lx, Ly, beta, Mx, My); + update_field_file("Caa", Caa, Na, Lx, Ly, beta, Mx, My); + update_field_file("Cll", Cll, N, Lx, Ly, beta, Mx, My); + update_field_file("Cee", Cee, N, Lx, Ly, beta, Mx, My); } -void ma::pre_fracture(const network &) { +void ma::pre_fracture(const network&) { lv = std::numeric_limits<long double>::lowest(); avalanches = {{}}; boost::remove_edge_if(trivial, G); @@ -297,13 +173,13 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned int i sa[avalanches.back().size() - 1]++; Na++; - std::fill_n(fftw_forward_in, net.G.edges.size(), 0.0); + std::fill_n(fftw_forward_in, Mx * My, 0.0); for (auto e : avalanches.back()) { - fftw_forward_in[e] = 1.0; + fftw_forward_in[edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0; } - autocorrelation(Lx, Ly, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(Mx, My, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); lv = c; avalanches.push_back({i}); @@ -311,121 +187,116 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned int i avalanches.back().push_back(i); } - boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], {i}, G); - - bondfile << i << " " << c << std::scientific << " "; + boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G); } void ma::post_fracture(network &n) { - bondfile << "\n"; std::vector<unsigned int> component(boost::num_vertices(G)); unsigned int num = boost::connected_components(G, &component[0]); std::list<unsigned int> crack = find_minimal_crack(G, n); - unsigned int crack_component = component[n.G.dual_edges[crack.front()][0]]; + // crack surface correlations + std::fill_n(fftw_forward_in, Mx * My, 0.0); + + for (auto edge : crack) { + fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My)] = 0.5; + fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My)] = 0.5; + } + + autocorrelation(Mx, My, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + + unsigned int crack_component = component[n.G.dual_edges[crack.front()].v[0]]; + + std::vector<std::list<unsigned int>> components(num); + + for (unsigned int i = 0; i < n.G.dual_vertices.size(); i++) { + components[component[i]].push_back(i); + } // non-spanning clusters + std::fill_n(fftw_forward_in, Mx * My, 0.0); for (unsigned int i = 0; i < num; i++) { if (i != crack_component) { - unsigned int size = 0; - - for (unsigned int j = 0; j < n.G.edges.size(); j++) { - if (component[n.G.dual_edges[j][0]] == i && n.fuses[j]) { - size++; - fftw_forward_in[j] = 1.0; - } else{ - fftw_forward_in[j] = 0.0; - } + for (auto it = components[i].begin(); it != components[i].end(); it++) { + fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] += 1.0; } - if (size > 0) { - sc[size - 1]++; - autocorrelation(Lx, Ly, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); - Nc++; - } - } - } + sc[components[i].size() - 1]++; + autocorrelation(Mx, My, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + Nc++; - // bin counting - for (unsigned int be = 0; be < log2(Ly); be++) { - unsigned int bin = pow(2, be); - - for (unsigned int i = 0; i < Lx * Ly / pow(bin, 2); i++) { - bool in_bin = false; - for (unsigned int j = 0; j < pow(bin, 2); j++) { - unsigned int edge = Lx * (bin * ((i * bin) / Lx) + j / bin) + (i * bin) % Lx + j % bin; - if (!n.fuses[edge] && crack_component == component[n.G.dual_edges[edge][0]]) { - in_bin = true; - break; - } - } - - if (in_bin) { - sb[be]++; + for (auto it = components[i].begin(); it != components[i].end(); it++) { + fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] = 0.0; } } } - sb[log2(Ly)]++; + unsigned int max_factor = log2(Mx < My ? Mx : My) + 1; + std::vector<std::valarray<unsigned int>> bins(max_factor); + for (unsigned int i = 0; i < max_factor; i++) { + bins[i].resize(Mx * My / (unsigned int)pow(2, 2 * i)); + } - for (unsigned int i = 0; i < n.G.edges.size(); i++) { - if (n.fuses[i] && component[n.G.dual_edges[i][0]] == crack_component) { - fftw_forward_in[i] = 1.0; - } else { - fftw_forward_in[i] = 0.0; + for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) { + for (unsigned int i = 0; i < max_factor; i++) { + bins[i][edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx / pow(2, i), My / pow(2, i))] = 1; } } - autocorrelation(Lx, Ly, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + for (unsigned int i =0 ; i < max_factor; i++) { + sb[i] += bins[i].sum(); + } - // crack surface correlations - std::fill_n(fftw_forward_in, n.G.edges.size(), 0.0); + // spanning cluster + std::fill_n(fftw_forward_in, Mx * My, 0.0); - for (auto edge : crack) { - fftw_forward_in[edge] = 1.0; + for (unsigned int i = 0; i < n.G.dual_vertices.size(); i++) { + if (component[i] == crack_component) { + fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[i].r, Lx, Ly, Mx, My)] += 1.0; + } } - autocorrelation(Lx, Ly, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(Mx, My, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); std::function<bool(unsigned int)> inCrack = [&](unsigned int i) -> bool { - return (bool)fftw_forward_in[i]; + return component[n.G.dual_edges[i].v[0]] == crack_component; }; for (auto avalanche : avalanches) { if (avalanche.end() != std::find_if(avalanche.begin(), avalanche.end(), inCrack)) { for (auto edge : avalanche) { - fftw_forward_in[edge] = 1.0; + fftw_forward_in[edge_r_to_ind(n.G.edges[edge].r, Lx, Ly, Mx, My)] += 1.0; } } } - autocorrelation(Lx, Ly, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(Mx, My, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); - std::fill_n(fftw_forward_in, n.G.edges.size(), 0.0); + std::fill_n(fftw_forward_in, Mx * My, 0.0); // rewind the last avalanche for (auto e : avalanches.back()) { - boost::remove_edge(n.G.dual_edges[e][0], n.G.dual_edges[e][1], G); + boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G); n.break_edge(e, true); - fftw_forward_in[e] = 1; + fftw_forward_in[edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0; } - autocorrelation(Lx, Ly, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(Mx, My, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); // damage size distribution unsigned int total_broken = 0; + std::fill_n(fftw_forward_in, Mx * My, 0.0); + for (unsigned int i = 0; i < n.G.edges.size(); i++) { if (n.fuses[i]) { total_broken++; - fftw_forward_in[i] = 1.0; - } else { - fftw_forward_in[i] = 0.0; + fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0; } } - autocorrelation(Lx, Ly, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(Mx, My, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); sd[total_broken]++; |