From cb1b2e6822bdd1d1644ff2dad2d6157858e105b0 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 24 Apr 2019 23:31:40 -0400 Subject: many changes to introduce two-component, elastic-like fracture --- src/analysis_tools.cpp | 74 +++++++++---------- src/analysis_tools.hpp | 2 +- src/animate.cpp | 32 ++++++-- src/animate_fracture.cpp | 6 +- src/animate_fracture_square.cpp | 4 +- src/fracture.cpp | 12 +-- src/fracture_elastic.cpp | 87 ++++++++++++++++++++++ src/fracture_square.cpp | 6 +- src/measurements.cpp | 160 ++++++++-------------------------------- src/measurements.hpp | 5 +- src/sample_fracture_square.cpp | 4 +- 11 files changed, 202 insertions(+), 190 deletions(-) create mode 100644 src/fracture_elastic.cpp (limited to 'src') diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp index bc8095e..dea20f0 100644 --- a/src/analysis_tools.cpp +++ b/src/analysis_tools.cpp @@ -18,7 +18,7 @@ bool trivial(boost::detail::edge_desc_impl) return true; } -std::list find_minimal_crack(const Graph& G, const network& n) { +std::list, std::list>> find_minimal_crack(const Graph& G, const network& n) { Graph Gtmp(n.G.vertices.size()); std::list removed_edges; @@ -79,53 +79,53 @@ std::list find_minimal_crack(const Graph& G, const network& n) { } } - 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); + 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 pos = 0; - for (uint8_t included : new_bool_cycle) { - if (included) { - new_cycle.push_back(pos); - } - pos++; + // 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 pos = 0; + for (uint8_t included : new_bool_cycle) { + if (included) { + new_cycle.push_back(pos); } - cycles.push_back(new_cycle); + 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}; + std::list, std::list>> output; - for (auto edge : cycle) { - if (n.G.dual_edges[edge].crossings[0]) { - crossing_count[0]++; - } - if (n.G.dual_edges[edge].crossings[1]) { - crossing_count[1]++; - } - } + // find the cycle representing the crack by counting boundary crossings + for (auto cycle : cycles) { + std::array crossing_count{0,0}; - if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) { - return cycle; + for (auto edge : cycle) { + if (n.G.dual_edges[edge].crossings[0]) { + crossing_count[0]++; + } + if (n.G.dual_edges[edge].crossings[1]) { + crossing_count[1]++; } } - } else { - return cycles.front(); + + if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) { + output.push_back({{1, 0}, cycle}); + } else if (crossing_count[0] % 2 == 0 && crossing_count[1] % 2 == 1) { + output.push_back({{0, 1}, cycle}); + } } - throw badcycleex; + return output; } diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp index 4f3f285..3beb1a8 100644 --- a/src/analysis_tools.hpp +++ b/src/analysis_tools.hpp @@ -25,5 +25,5 @@ bool is_shorter(const std::list &, const std::list &); bool trivial(boost::detail::edge_desc_impl); -std::list find_minimal_crack(const Graph &, const network &); +std::list, std::list>> find_minimal_crack(const Graph &, const network &); diff --git a/src/animate.cpp b/src/animate.cpp index 5bae15e..be0beef 100644 --- a/src/animate.cpp +++ b/src/animate.cpp @@ -1,5 +1,6 @@ #include "animate.hpp" +#include animate::animate(double Lx, double Ly, unsigned window_size, int argc, char *argv[]) : G(2 * (unsigned)ceil(Lx * Ly / 2)) { glutInit(&argc, argv); @@ -85,20 +86,41 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i) } void animate::post_fracture(network &n) { + + std::list crack; +// unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]]; + unsigned crack_component = 10000; + + std::default_random_engine gen; + std::uniform_real_distribution dis(0.0,1.0); + + bool cycle_present = true; + auto av_it = avalanches.rbegin(); + + while (cycle_present) { + for (unsigned e : *av_it) { + boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G); + n.fuses[e] = false; + } + + auto cracks = find_minimal_crack(G, n); + std::cout << cracks.size() << "\n"; + if (cracks.size() == 0) { + cycle_present = false; + } + + av_it++; + } + std::vector component(boost::num_vertices(G)); unsigned num = boost::connected_components(G, &component[0]); - std::list crack = find_minimal_crack(G, n); - unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]]; - std::vector> components(num); for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) { components[component[i]].push_back(i); } - std::default_random_engine gen; - std::uniform_real_distribution dis(0.0,1.0); char key; while ((key = getchar()) != 'n') { diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp index 66afc83..c3b4a69 100644 --- a/src/animate_fracture.cpp +++ b/src/animate_fracture.cpp @@ -64,9 +64,9 @@ int main(int argc, char* argv[]) { for (unsigned trial = 0; trial < N; trial++) { graph G(Lx, Ly, rng); - network network(G, &c); - network.set_thresholds(beta, rng); - network.fracture(meas); + elastic_network elastic_network(G, &c); + elastic_network.set_thresholds(beta, rng); + elastic_network.fracture(meas); if (quit.load()) break; diff --git a/src/animate_fracture_square.cpp b/src/animate_fracture_square.cpp index d2efde9..84cc8a5 100644 --- a/src/animate_fracture_square.cpp +++ b/src/animate_fracture_square.cpp @@ -63,10 +63,10 @@ int main(int argc, char* argv[]) { std::mt19937 rng{seeds}; graph G(Lx, Ly); - network perm_network(G, &c); + elastic_network perm_network(G, &c); for (unsigned trial = 0; trial < N; trial++) { - network tmp_network(perm_network); + elastic_network tmp_network(perm_network); tmp_network.set_thresholds(beta, rng); tmp_network.fracture(meas); diff --git a/src/fracture.cpp b/src/fracture.cpp index f697617..fb51d3e 100644 --- a/src/fracture.cpp +++ b/src/fracture.cpp @@ -82,9 +82,9 @@ int main(int argc, char* argv[]) { while (true) { try { graph G(n, a, rng); - network network(G, &c); - network.set_thresholds(beta, rng); - network.fracture(meas); + elastic_network fuse_network(G, &c); + fuse_network.set_thresholds(beta, rng); + fuse_network.fracture(meas); break; } catch (std::exception &e) { std::cout << e.what() << '\n'; @@ -105,9 +105,9 @@ int main(int argc, char* argv[]) { while (true) { try { graph G(Lx, Ly, rng); - network network(G, &c); - network.set_thresholds(beta, rng); - network.fracture(meas); + elastic_network fuse_network(G, &c); + fuse_network.set_thresholds(beta, rng); + fuse_network.fracture(meas); break; } catch (std::exception &e) { std::cout << e.what() << '\n'; diff --git a/src/fracture_elastic.cpp b/src/fracture_elastic.cpp new file mode 100644 index 0000000..dec4a2c --- /dev/null +++ b/src/fracture_elastic.cpp @@ -0,0 +1,87 @@ + +#include +#include + +#include + +#include "randutils/randutils.hpp" + +#include +#include +#include +#include "measurements.hpp" + +#include +#include +#include + +std::atomic quit(false); // signal flag + +void got_signal(int) { + quit.store(true); +} + +int main(int argc, char* argv[]) { + struct sigaction sa; + memset( &sa, 0, sizeof(sa) ); + sa.sa_handler = got_signal; + sigfillset(&sa.sa_mask); + sigaction(SIGINT, &sa, NULL); + + int opt; + + unsigned N = 1; + unsigned Lx = 16; + unsigned Ly = 16; + double beta = 0.5; + + while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) { + switch (opt) { + case 'N': + N = (unsigned)atof(optarg); + break; + case 'X': + Lx = atoi(optarg); + break; + case 'Y': + Ly = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + default: + exit(1); + } + } + + cholmod_common c; + CHOL_F(start)(&c); + + ma meas(Lx, Ly, 2*Lx, 2*Ly, beta); + graph G(Lx, Ly); + network perm_network(G, &c); + + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + for (unsigned trial = 0; trial < N; trial++) { + while (true) { + try { + network tmp_network(perm_network); + tmp_network.set_thresholds(beta, rng); + tmp_network.fracture(meas); + break; + } catch (std::exception &e) { + std::cout << e.what() << '\n'; + } + } + + if (quit.load()) + break; + } + + CHOL_F(finish)(&c); + + return 0; +} + diff --git a/src/fracture_square.cpp b/src/fracture_square.cpp index dec4a2c..96e6030 100644 --- a/src/fracture_square.cpp +++ b/src/fracture_square.cpp @@ -57,9 +57,9 @@ int main(int argc, char* argv[]) { cholmod_common c; CHOL_F(start)(&c); - ma meas(Lx, Ly, 2*Lx, 2*Ly, beta); + ma meas(Lx, Ly, beta); graph G(Lx, Ly); - network perm_network(G, &c); + elastic_network perm_network(G, &c); randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; @@ -67,7 +67,7 @@ int main(int argc, char* argv[]) { for (unsigned trial = 0; trial < N; trial++) { while (true) { try { - network tmp_network(perm_network); + elastic_network tmp_network(perm_network); tmp_network.set_thresholds(beta, rng); tmp_network.fracture(meas); break; diff --git a/src/measurements.cpp b/src/measurements.cpp index 2af1ba3..f2598ee 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -3,15 +3,13 @@ #include #include -void update_distribution_file(std::string id, const std::vector& data, unsigned N, std::string model_string) { +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); - uint64_t N_old = 0; std::vector data_old(data.size(), 0); if (file.is_open()) { - file >> N_old; for (unsigned i = 0; i < data.size(); i++) { uint64_t num; file >> num; @@ -23,7 +21,6 @@ void update_distribution_file(std::string id, const std::vector& data, std::ofstream file_out(filename); - file_out <::lowest(); - avalanches = {{}}; boost::remove_edge_if(trivial, G); + avalanches = {}; } void ma::bond_broken(const network& net, const current_info& cur, unsigned i) { long double c = logl(cur.conductivity / fabs(cur.currents[i])) + net.thresholds[i]; - if (c > lv && avalanches.back().size() > 0) { - sa[avalanches.back().size() - 1]++; - Na++; - - autocorrelation2(net.G.L.x, net.G.L.y, Mx, My, caa, avalanches.back()); - + if (c > lv) { + if (avalanches.size() > 0) { + sa[avalanches.back().size() - 1]++; + } lv = c; - avalanches.push_back({net.G.edges[i].r}); - last_avalanche = {i}; + avalanches.push_back({i}); } else { - avalanches.back().push_back(net.G.edges[i].r); - last_avalanche.push_back(i); + avalanches.back().push_back(i); } 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) { + auto post_cracks = find_minimal_crack(G, n); std::vector component(boost::num_vertices(G)); unsigned num = boost::connected_components(G, &component[0]); - - std::list crack = find_minimal_crack(G, n); - - ss[crack.size() - 1]++; - std::list sr; - - for (auto edge : crack) { - sr.push_back(n.G.dual_edges[edge].r); - } - - autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, css, sr); - - unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]]; + unsigned crack_component = component[n.G.dual_edges[post_cracks.front().second.front()].v[0]]; std::vector> components(num); @@ -303,43 +229,26 @@ void ma::post_fracture(network &n) { } for (unsigned i = 0; i < num; i++) { - if (i != crack_component && components[i].size() > 0) { - sc[components[i].size() - 1]++; - Nc++; - - autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ccc, components[i]); + if (i != crack_component) { + sm[components[i].size() - 1]++; } } - // spanning cluster - - sm[components[crack_component].size() - 1]++; - - autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cmm, components[crack_component]); + auto av_it = avalanches.rbegin(); - /// damage at end - std::list Dr; - - for (unsigned i = 0; i < n.G.edges.size(); i++) { - if (n.fuses[i]) { - Dr.push_back(n.G.edges[i].r); + while (true) { + for (unsigned e : *av_it) { + boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G); + n.fuses[e] = false; } - } - - sD[Dr.size() - 1]++; - autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cDD, Dr); - correlation2(n.G.L.x, n.G.L.y, Mx, My, csD, sr, Dr); + auto cracks = find_minimal_crack(G, n); - // ********************** LAST AVALANCHE ********************* - - // rewind the last avalanche - sl[avalanches.back().size() - 1]++; - - autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cll, avalanches.back()); + if (cracks.size() == 0) { + break; + } - for (unsigned e : last_avalanche) { - boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G); + av_it++; } num = boost::connected_components(G, &component[0]); @@ -351,14 +260,7 @@ void ma::post_fracture(network &n) { } for (unsigned i = 0; i < num; i++) { - if (new_components[i].size() > 0) { - sb[new_components[i].size() - 1]++; - Nb++; - - autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cbb, new_components[i]); - } + sc[new_components[i].size() - 1]++; } - - N++; } diff --git a/src/measurements.hpp b/src/measurements.hpp index 610c266..7c9d49c 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -35,6 +35,7 @@ class ma : public hooks { // measurement storage std::vector sc; // non-spanning cluster size distribution + std::vector sn; // non-spanning cluster size distribution std::vector ss; // minimal spanning cluster size distribution std::vector sm; // spanning cluster size distribution std::vector sa; // non-final avalanche size distribution @@ -56,11 +57,11 @@ class ma : public hooks { public: long double lv; - std::list> avalanches; + std::list> avalanches; std::list last_avalanche; std::string model_string; - ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta); + ma(unsigned Lx, unsigned Ly, double beta); ma(unsigned n, double a, unsigned Mx, unsigned My, double beta); ~ma(); diff --git a/src/sample_fracture_square.cpp b/src/sample_fracture_square.cpp index fb2eb33..d4bef2c 100644 --- a/src/sample_fracture_square.cpp +++ b/src/sample_fracture_square.cpp @@ -52,10 +52,10 @@ int main(int argc, char* argv[]) { std::mt19937 rng{seeds}; graph G(Lx, Ly); - network perm_network(G, &c); + elastic_network perm_network(G, &c); for (unsigned trial = 0; trial < N; trial++) { - network tmp_network(perm_network); + elastic_network tmp_network(perm_network); tmp_network.set_thresholds(beta, rng); tmp_network.fracture(meas); } -- cgit v1.2.3-70-g09d2