From 7984b8dbde1d44fc18b6026693c387d2d433872d Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 12 Nov 2018 10:52:17 -0500 Subject: changes --- lib/src/network.cpp | 8 +- src/CMakeLists.txt | 2 +- src/measurements.cpp | 344 ++++++++++++++++++++++++++------------------------- src/measurements.hpp | 6 +- 4 files changed, 188 insertions(+), 172 deletions(-) diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 9addfdd..a27ca39 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -1,7 +1,7 @@ #include "network.hpp" -network::network(const graph& G, cholmod_common *c) : G(G), c(c), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) { +network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) { b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c); for (unsigned int i = 0; i < G.edges.size(); i++) { double v0y = G.vertices[G.edges[i][0]].r.y; @@ -36,8 +36,8 @@ network::network(const graph& G, cholmod_common *c) : G(G), c(c), fuses(G.edges. unsigned int s0 = v0 < v1 ? v0 : v1; unsigned int s1 = v0 < v1 ? v1 : v0; - ((CHOL_INT *)t->i)[terms] = v1; - ((CHOL_INT *)t->j)[terms] = v0; + ((CHOL_INT *)t->i)[terms] = s1; + ((CHOL_INT *)t->j)[terms] = s0; ((double *)t->x)[terms] = -1.0; ((double *)t->x)[v0]++; @@ -74,7 +74,7 @@ network::network(const graph& G, cholmod_common *c) : G(G), c(c), fuses(G.edges. CHOL_F(free_triplet)(&t, c); } -network::network(const network &other) : G(other.G), thresholds(other.thresholds), fuses(other.fuses), c(other.c) { +network::network(const network &other) : c(other.c), G(other.G), fuses(other.fuses), thresholds(other.thresholds) { b = CHOL_F(copy_dense)(other.b, c); factor = CHOL_F(copy_factor)(other.factor, c); voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 57a817b..d93783c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -3,5 +3,5 @@ add_library(measurements measurements.cpp) target_link_libraries(measurements frac) add_executable(fracture fracture.cpp) -target_link_libraries(fracture frac measurements fftw3 cholmod) +target_link_libraries(fracture frac measurements fftw3 cholmod profiler) diff --git a/src/measurements.cpp b/src/measurements.cpp index 578c70a..0bd72c4 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -11,6 +11,119 @@ bool trivial(boost::detail::edge_desc_impl) return true; } +std::list find_minimal_crack(const Graph& G, const network& n) { + Graph Gtmp(n.G.vertices.size()); + std::list removed_edges; + + class add_tree_edges : public boost::default_dfs_visitor { + public: + Graph& G; + std::list& E; + + add_tree_edges(Graph& G, std::list& E) : G(G), E(E) {} + + void tree_edge(boost::graph_traits::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::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: + unsigned int end; + std::list& E; + struct done{}; + + find_cycle(std::list& E, unsigned int end) : E(E), end(end) {} + + void discover_vertex(boost::graph_traits::vertex_descriptor v, const Graph& g) { + if (v == end) { + throw done{}; + } + } + + void examine_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { + E.push_back(g[e].index); + } + + void finish_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { + E.erase(std::find(E.begin(), E.end(), g[e].index)); + } + }; + + std::list> cycles; + + for (auto edge : removed_edges) { + std::list cycle = {edge}; + find_cycle vis(cycle, n.G.dual_edges[edge][1]); + std::vector 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> bool_cycles; + for (auto cycle : cycles) { + std::valarray 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 new_bool_cycle = (*it1) ^ (*it2); + std::list 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 crossing_count{0,0}; + + for (auto edge : cycle) { + double dx = fabs(n.G.vertices[n.G.dual_edges[edge][0]].r.x - n.G.vertices[n.G.dual_edges[edge][1]].r.x); + if (dx > 0.5) { + crossing_count[0]++; + } + double dy = fabs(n.G.vertices[n.G.dual_edges[edge][0]].r.y - n.G.vertices[n.G.dual_edges[edge][1]].r.y); + if (dy > 0.5) { + 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& data, unsigned int N, unsigned int L, double beta) { std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat"; std::ifstream file(filename); @@ -97,12 +210,8 @@ void correlation(unsigned int L, std::vector& data, const std::vector -void autocorrelation(unsigned int L, std::vector& out_data, const std::vector& 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) { - for (unsigned int i = 0; i < pow(L, 2); i++) { - fftw_forward_in[i] = (double)data[i]; - } - +template +void autocorrelation(unsigned int L, 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 int i = 0; i < L * (L / 2 + 1); i++) { @@ -124,8 +233,11 @@ ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2)) sb(log2(L) + 1, 0), Ccc(pow(L / 2 + 1, 2), 0), Css(pow(L / 2 + 1, 2), 0), + Cmm(pow(L / 2 + 1, 2), 0), + Caa(pow(L / 2 + 1, 2), 0), Cdd(pow(L / 2 + 1, 2), 0), - Caa(pow(L / 2 + 1, 2), 0) + Cll(pow(L / 2 + 1, 2), 0), + Cee(pow(L / 2 + 1, 2), 0) { Nc = 0; Na = 0; @@ -159,34 +271,36 @@ ma::~ma() { update_field_file("Ccc", Ccc, Nc, L, beta); update_field_file("Css", Css, N, L, beta); + update_field_file("Cmm", Cmm, N, L, beta); update_field_file("Cdd", Cdd, N, L, beta); update_field_file("Caa", Caa, Na, L, beta); - + update_field_file("Cll", Cll, N, L, beta); + update_field_file("Cee", Cee, N, L, beta); } void ma::pre_fracture(const network &) { lv = 0; - avalanche.clear(); + avalanches = {{}}; boost::remove_edge_if(trivial, G); } void ma::bond_broken(const network& net, const std::pair>& cur, unsigned int i) { if (cur.first / fabs(cur.second[i]) * net.thresholds[i] > lv) { - sa[avalanche.size()]++; + sa[avalanches.back().size()]++; Na++; - std::vectoravalanche_damage(pow(L, 2), 0); + memset(fftw_forward_in, 0.0, net.G.edges.size()); - for (auto it = avalanche.begin(); it != avalanche.end(); it++) { - avalanche_damage[*it] = 1; + for (auto e : avalanches.back()) { + fftw_forward_in[e] = 1.0; } - autocorrelation(L, Caa, avalanche_damage, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(L, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); lv = cur.first / fabs(cur.second[i]) * net.thresholds[i]; - avalanche.clear(); + avalanches.push_back({i}); } else { - avalanche.push_back(i); + avalanches.back().push_back(i); } boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], {i}, G); @@ -196,135 +310,9 @@ void ma::post_fracture(network &n) { std::vector component(boost::num_vertices(G)); unsigned int num = boost::connected_components(G, &component[0]); - std::vector comp_sizes(num, 0); - - for (unsigned int c : component) { - comp_sizes[c]++; - } + std::list crack = find_minimal_crack(G, n); - unsigned int max_i = 0; - unsigned int max_size = 0; - - for (unsigned int i = 0; i < num; i++) { - if (comp_sizes[i] > max_size) { - max_i = i; - max_size = comp_sizes[i]; - } - } - - unsigned int vertex_in_largest = 0; - - while (component[vertex_in_largest] != max_i) { - vertex_in_largest++; - } - - Graph Gtmp(2 * pow(L / 2, 2)); - std::list removed_edges; - class add_tree_edges : public boost::default_dfs_visitor { - public: - Graph& G; - std::list& E; - - add_tree_edges(Graph& G, std::list& E) : G(G), E(E) {} - - void tree_edge(boost::graph_traits::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::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: - unsigned int end; - std::list& E; - struct done{}; - - find_cycle(std::list& E, unsigned int end) : E(E), end(end) {} - - void discover_vertex(boost::graph_traits::vertex_descriptor v, const Graph& g) { - if (v == end) { - throw done{}; - } - } - - void examine_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { - E.push_back(g[e].index); - } - - void finish_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { - E.erase(std::find(E.begin(), E.end(), g[e].index)); - } - }; - - std::list> cycles; - - for (auto edge : removed_edges) { - std::list cycle = {edge}; - find_cycle vis(cycle, n.G.dual_edges[edge][1]); - std::vector 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); - } - } - - std::list crack; - - if (cycles.size() > 1) { - std::list> bool_cycles; - for (auto cycle : cycles) { - std::valarray bool_cycle(n.G.edges.size()); - for (auto v : cycle) { - bool_cycle[v] = 1; - } - bool_cycles.push_back(bool_cycle); - } - - 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 new_bool_cycle = (*it1) ^ (*it2); - std::list 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); - } - } - - for (auto cycle : cycles) { - std::array crossing_count{0,0}; - - for (auto edge : cycle) { - double dx = fabs(n.G.vertices[n.G.edges[edge][0]].r.x - n.G.vertices[n.G.edges[edge][1]].r.x); - if (dx > 0.5) { - crossing_count[0]++; - } - double dy = fabs(n.G.vertices[n.G.edges[edge][0]].r.y - n.G.vertices[n.G.edges[edge][1]].r.y); - if (dy > 0.5) { - crossing_count[1]++; - } - } - - if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) { - crack = cycle; - break; - } - } - } else { - crack = cycles.front(); - } + unsigned int crack_component = component[n.G.dual_edges[crack.front()][0]]; for (unsigned int be = 0; be < log2(L); be++) { unsigned int bin = pow(2, be); @@ -333,7 +321,7 @@ void ma::post_fracture(network &n) { bool in_bin = false; for (unsigned int j = 0; j < pow(bin, 2); j++) { unsigned int edge = L * (bin * ((i * bin) / L) + j / bin) + (i * bin) % L + j % bin; - if (!n.fuses[edge] && max_i == component[n.G.dual_edges[edge][0]]) { + if (!n.fuses[edge] && crack_component == component[n.G.dual_edges[edge][0]]) { in_bin = true; break; } @@ -347,62 +335,86 @@ void ma::post_fracture(network &n) { sb[log2(L)]++; + 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; + } + } + + autocorrelation(L, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + // crack surface correlations - std::vector crack_damage(pow(L, 2), 0.0); + memset(fftw_forward_in, 0.0, n.G.edges.size()); for (auto edge : crack) { - crack_damage[edge] = 1.0; + fftw_forward_in[edge] = 1.0; + } + + autocorrelation(L, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + + std::function inCrack = [&](unsigned int i) -> bool { + return (bool)fftw_forward_in[i]; + }; + + 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; + } + } } - std::vector t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out); - correlation(L, Css, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(L, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); - std::vector last_crack(pow(L, 2), 0); + memset(fftw_forward_in, 0.0, n.G.edges.size()); // rewind the last avalanche - for (auto e = avalanche.rbegin(); e != avalanche.rend(); e++) { - boost::remove_edge(n.G.dual_edges[*e][0], n.G.dual_edges[*e][1], G); - n.add_edge(*e); - last_crack[*e] = 1; + for (auto e : avalanches.back()) { + boost::remove_edge(n.G.dual_edges[e][0], n.G.dual_edges[e][1], G); + n.add_edge(e); + fftw_forward_in[e] = 1; } + autocorrelation(L, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); // cluster size distribution and cluster-cluster correlation num = boost::connected_components(G, &component[0]); for (unsigned int i = 0; i < num; i++) { - double size = 0; - std::fill(crack_damage.begin(), crack_damage.end(), 0.0); + unsigned int size = 0; - for (unsigned int j = 0; j < pow(L, 2); j++) { - if (component[n.G.edges[j][0]] == i && n.fuses[j]) { + for (unsigned int j = 0; j < n.G.edges.size(); j++) { + if (component[n.G.dual_edges[j][0]] == i && n.fuses[j]) { size++; - crack_damage[j] = 1.0; + fftw_forward_in[j] = 1.0; + } else{ + fftw_forward_in[j] = 0.0; } } if (size > 0) { sc[size - 1]++; - t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out); - correlation(L, Ccc, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(L, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); Nc++; } } - // damage-damage correlations - std::vector t_damage = data_transform(L, n.fuses, forward_plan, fftw_forward_in, fftw_forward_out); - correlation(L, Cdd, t_damage, t_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); - - // damage size distribution unsigned int total_broken = 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; } } + autocorrelation(L, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + sd[total_broken]++; } diff --git a/src/measurements.hpp b/src/measurements.hpp index 4c05d93..621c425 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include #include @@ -49,8 +50,11 @@ class ma : public hooks { std::vector sb; // bin size distribution std::vector Ccc; // cluster-cluster correlations std::vector Css; // surface-surface correlations + std::vector Cmm; // surface-surface correlations std::vector Caa; // avalanche-avalanche correlations std::vector Cdd; // damage-damage distribution + std::vector Cll; // damage-damage distribution + std::vector Cee; // damage-damage distribution uint64_t Nc; uint64_t Na; @@ -61,7 +65,7 @@ class ma : public hooks { double beta; - std::list avalanche; + std::list> avalanches; ma(unsigned int N, unsigned int L, double beta); ~ma(); -- cgit v1.2.3-54-g00ecf