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 --- CMakeLists.txt | 4 +- lib/include/graph.hpp | 1 + lib/include/network.hpp | 50 ++++++++-- lib/src/network.cpp | 204 +++++++++++++++++++++++++++------------- 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 +- 15 files changed, 386 insertions(+), 265 deletions(-) create mode 100644 src/fracture_elastic.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 35e8436..e3f70d8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,8 +1,8 @@ cmake_minimum_required(VERSION 3.8) -set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall -stdlib=libc++ -fuse-ld=lld") -set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -stdlib=libc++ -march=native -fuse-ld=lld") +set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall") +set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -march=native") set (CMAKE_CXX_STANDARD 17) set (CMAKE_C_STANDARD 11) diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp index 1714dc2..45cc478 100644 --- a/lib/include/graph.hpp +++ b/lib/include/graph.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include "array_hash.hpp" diff --git a/lib/include/network.hpp b/lib/include/network.hpp index f12c9c7..ba98086 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -26,25 +26,55 @@ #endif -class network { +class problem { private: - cholmod_dense *b; - cholmod_factor *factor; - cholmod_sparse *voltcurmat; - cholmod_common *c; + const graph& G; + cholmod_dense* b; + cholmod_factor* factor; + cholmod_sparse* voltcurmat; + cholmod_common* c; + unsigned axis; + public: + problem(const graph&, unsigned axis, cholmod_common*); + problem(const graph&, unsigned axis, cholmod_sparse*, cholmod_common*); + problem(const problem&); + ~problem(); + current_info solve(std::vector& fuses); + void break_edge(unsigned, bool unbreak = false); +}; + +class network { public: const graph& G; std::vector fuses; std::vector thresholds; - network(const graph&, cholmod_common*); - network(const network &other); - ~network(); + network(const graph&); + network(const network&); void set_thresholds(double, std::mt19937&); - void break_edge(unsigned, bool unbreak = false); - current_info get_current_info(); +}; + +class fuse_network : public network { + public: + problem ohm; + + fuse_network(const graph&, cholmod_common*); + void fracture(hooks&, double cutoff = 1e-13); + current_info get_current_info(); +}; + +class elastic_network : public network { + public: + problem hook_x; + problem hook_y; + + elastic_network(const graph&, cholmod_common*); + elastic_network(const elastic_network&); + + void fracture(hooks&, double weight = 0.5, double cutoff = 1e-13); + current_info get_current_info(); }; diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 1edc973..8b0e8df 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -17,14 +17,19 @@ class nofuseException: public std::exception } } nofuseex; -network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) { +problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), voltcurmat(vcmat), c(c), axis(axis) { b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c); for (unsigned i = 0; i < G.edges.size(); i++) { - double v0y = G.vertices[G.edges[i].v[0]].r.y; - double v1y = G.vertices[G.edges[i].v[1]].r.y; - - if (G.edges[i].crossings[1]) { - bool ind = v0y < v1y ? 0 : 1; + graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r; + graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r; + + if (G.edges[i].crossings[axis]) { + bool ind; + if (axis == 1) { + ind = v0.y < v1.y ? 0 : 1; + } else { + ind = v0.x < v1.x ? 0 : 1; + } ((double *)b->x)[G.edges[i].v[ind]] += 1.0; ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0; @@ -78,8 +83,10 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges. factor = CHOL_F(analyze)(laplacian, c); CHOL_F(factorize)(laplacian, factor, c); CHOL_F(free_sparse)(&laplacian, c); +} - t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c); +problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) { + cholmod_triplet *t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c); t->nnz = 2 * G.edges.size(); @@ -98,39 +105,70 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges. CHOL_F(free_triplet)(&t, c); } -network::network(const network &other) : c(other.c), G(other.G), fuses(other.fuses), thresholds(other.thresholds) { +problem::problem(const problem& other) : G(other.G), c(other.c), axis(other.axis) { b = CHOL_F(copy_dense)(other.b, c); factor = CHOL_F(copy_factor)(other.factor, c); voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c); } -network::~network() { - CHOL_F(free_sparse)(&voltcurmat, c); +problem::~problem() { CHOL_F(free_dense)(&b, c); CHOL_F(free_factor)(&factor, c); + CHOL_F(free_sparse)(&voltcurmat, c); } -void network::set_thresholds(double beta, std::mt19937& rng) { - if (beta == 0.0) { - /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */ - for (long double& threshold : thresholds) { - threshold = 1.0; - } - } else { - std::uniform_real_distribution d(0.0, 1.0); +current_info problem::solve(std::vector& fuses) { + cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); - for (long double& threshold : thresholds) { - threshold = std::numeric_limits::lowest(); + if (((double *)x->x)[0] != ((double *)x->x)[0]) { + throw nanex; + } - while (threshold == std::numeric_limits::lowest()) { - threshold = logl(d(rng)) / (long double)beta; + cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c); + + double alpha[2] = {1, 0}; + double beta[2] = {0, 0}; + CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); + + std::vector currents(G.edges.size()); + + double total_current = 0; + + for (int i = 0; i < G.edges.size(); i++) { + if (fuses[i]) { + currents[i] = 0; + } else { + currents[i] = ((double *)y->x)[i]; + + graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r; + graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r; + + if (G.edges[i].crossings[axis]) { + bool comp; + if (axis == 1) { + comp = v0.y > v1.y; + } else { + comp = v0.x > v1.x; + } + + if (comp) { + currents[i] += 1.0; + total_current += currents[i]; + } else { + currents[i] -= 1.0; + total_current -= currents[i]; + } } } } + + CHOL_F(free_dense)(&x, c); + CHOL_F(free_dense)(&y, c); + + return {total_current, currents}; } -void network::break_edge(unsigned e, bool unbreak) { - fuses[e] = !unbreak; +void problem::break_edge(unsigned e, bool unbreak) { unsigned v0 = G.edges[e].v[0]; unsigned v1 = G.edges[e].v[1]; @@ -167,78 +205,116 @@ void network::break_edge(unsigned e, bool unbreak) { CHOL_F(free_sparse)(&perm_update_mat, c); CHOL_F(free_sparse)(&update_mat, c); - double v0y = G.vertices[v0].r.y; - double v1y = G.vertices[v1].r.y; + graph::coordinate r0 = G.vertices[v0].r; + graph::coordinate r1 = G.vertices[v1].r; - if (G.edges[e].crossings[1]) { - bool ind = v0y < v1y ? unbreak : !unbreak; + if (G.edges[e].crossings[axis]) { + bool ind; + if (axis == 1) { + ind = r0.y < r1.y ? unbreak : !unbreak; + } else { + ind = r0.x < r1.x ? unbreak : !unbreak; + } ((double *)b->x)[G.edges[e].v[ind]] -= 1.0; ((double *)b->x)[G.edges[e].v[!ind]] += 1.0; } } -current_info network::get_current_info() { - cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); +network::network(const graph& G) : G(G), fuses(G.edges.size()), thresholds(G.edges.size()) {} - if (((double *)x->x)[0] != ((double *)x->x)[0]) { - throw nanex; - } +network::network(const network& n) : G(n.G), fuses(n.fuses), thresholds(n.thresholds) {} - cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c); +void network::set_thresholds(double beta, std::mt19937& rng) { + if (beta == 0.0) { + /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */ + for (long double& threshold : thresholds) { + threshold = 1.0; + } + } else { + std::uniform_real_distribution d(0.0, 1.0); - double alpha[2] = {1, 0}; - double beta[2] = {0, 0}; - CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); + for (long double& threshold : thresholds) { + threshold = std::numeric_limits::lowest(); - std::vector currents(G.edges.size()); + while (threshold == std::numeric_limits::lowest()) { + threshold = logl(d(rng)) / (long double)beta; + } + } + } +} - double total_current = 0; +fuse_network::fuse_network(const graph& G, cholmod_common *c) : network(G), ohm(G, 1, c) { +} - for (int i = 0; i < G.edges.size(); i++) { - if (fuses[i]) { - currents[i] = 0; - } else { - currents[i] = ((double *)y->x)[i]; +void fuse_network::fracture(hooks& m, double cutoff) { + m.pre_fracture(*this); - double v0y = G.vertices[G.edges[i].v[0]].r.y; - double v1y = G.vertices[G.edges[i].v[1]].r.y; + while (true) { + current_info ci = ohm.solve(fuses); - if (G.edges[i].crossings[1]) { - if (v0y > v1y) { - currents[i] += 1.0; - total_current += currents[i]; - } else { - currents[i] -= 1.0; - total_current -= currents[i]; + if (ci.conductivity < cutoff * G.vertices.size()) { + break; + } + + unsigned max_pos = UINT_MAX; + long double max_val = std::numeric_limits::lowest(); + + for (unsigned i = 0; i < G.edges.size(); i++) { + if (!fuses[i] && fabs(ci.currents[i]) > cutoff) { + long double val = logl(fabs(ci.currents[i])) - thresholds[i]; + if (val > max_val) { + max_val = val; + max_pos = i; } } } + if (max_pos == UINT_MAX) { + throw nofuseex; + } + + fuses[max_pos] = true; + ohm.break_edge(max_pos); + + m.bond_broken(*this, ci, max_pos); } - CHOL_F(free_dense)(&x, c); - CHOL_F(free_dense)(&y, c); + m.post_fracture(*this); +} - return {total_current, currents}; +elastic_network::elastic_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) { +} + +elastic_network::elastic_network(const elastic_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) { } -void network::fracture(hooks& m, double cutoff) { +void elastic_network::fracture(hooks& m, double weight, double cutoff) { m.pre_fracture(*this); while (true) { - current_info ci = this->get_current_info(); + current_info cx = hook_x.solve(fuses); + current_info cy = hook_y.solve(fuses); - if (ci.conductivity < cutoff * G.vertices.size()) { + current_info ctot; + + ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); + + if (ctot.conductivity < cutoff * G.vertices.size()) { break; } + ctot.currents.resize(G.edges.size()); + for (unsigned i = 0; i < G.edges.size(); i++) { + ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i], 2) + pow(weight * cy.currents[i], 2)); + } + unsigned max_pos = UINT_MAX; long double max_val = std::numeric_limits::lowest(); for (unsigned i = 0; i < G.edges.size(); i++) { - if (!fuses[i] && fabs(ci.currents[i]) > cutoff) { - long double val = logl(fabs(ci.currents[i])) - thresholds[i]; + if (!fuses[i] && fabs(ctot.currents[i]) > cutoff) { + long double val = logl(fabs(ctot.currents[i])) - thresholds[i]; if (val > max_val) { max_val = val; max_pos = i; @@ -250,9 +326,11 @@ void network::fracture(hooks& m, double cutoff) { throw nofuseex; } - this->break_edge(max_pos); + fuses[max_pos] = true; + hook_x.break_edge(max_pos); + hook_y.break_edge(max_pos); - m.bond_broken(*this, ci, max_pos); + m.bond_broken(*this, ctot, max_pos); } m.post_fracture(*this); 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