diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/measurements.cpp | 53 | ||||
| -rw-r--r-- | src/measurements.hpp | 32 | ||||
| -rw-r--r-- | src/perc_meas.cpp | 110 | ||||
| -rw-r--r-- | src/perc_meas.hpp | 10 | ||||
| -rw-r--r-- | src/percolation.cpp | 2 | ||||
| -rw-r--r-- | src/sample.cpp | 4 | ||||
| -rw-r--r-- | src/sample_fracture.cpp | 8 |
7 files changed, 97 insertions, 122 deletions
diff --git a/src/measurements.cpp b/src/measurements.cpp index 33056ef..aa0617c 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -144,13 +144,15 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u ma::ma(unsigned n, double a, unsigned Mx, unsigned My, double beta) : G(2 * n), - sc(2 * n), - sm(2 * n), - sa(3 * n), - sl(2 * n), sn(2 * n), ss(2 * n), - sb(3 * n) + sm(2 * n), + sl(2 * n), + sb(n), + sd(3 * n), + sc(2 * n), + sa(3 * n), + sA(3 * n) { if (beta != 0.0) { model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_"; @@ -161,11 +163,15 @@ ma::ma(unsigned n, double a, unsigned Mx, unsigned My, double beta) : ma::ma(unsigned Lx, unsigned Ly, double beta) : G(Lx * Ly / 2), - sc(Lx * Ly / 2), + sn(Lx * Ly / 2), + ss(Lx * Ly / 2), sm(Lx * Ly / 2), - sa(Lx * Ly), sl(Lx * Ly / 2), - sn(Lx * Ly / 2) + sb(Lx * Ly / 2), + sd(Lx * Ly), + sc(Lx * Ly / 2), + sa(Lx * Ly), + sA(Lx * Ly) { if (beta != 0.0) { model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_"; @@ -175,19 +181,22 @@ ma::ma(unsigned Lx, unsigned Ly, double beta) : } ma::~ma() { - update_distribution_file("sc", sc, model_string); - update_distribution_file("sm", sm, model_string); - update_distribution_file("sa", sa, model_string); - update_distribution_file("sl", sl, 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("sc", sc, model_string); + update_distribution_file("sa", sa, model_string); + update_distribution_file("sA", sA, model_string); } void ma::pre_fracture(const network&) { lv = std::numeric_limits<long double>::lowest(); boost::remove_edge_if(trivial, G); avalanches = {}; + num = 0; } void ma::bond_broken(const network& net, const current_info& cur, unsigned i) { @@ -200,6 +209,7 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) { } boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G); + num++; } void ma::post_fracture(network &n) { @@ -207,6 +217,7 @@ void ma::post_fracture(network &n) { std::vector<unsigned> component(boost::num_vertices(G)); unsigned num = boost::connected_components(G, &component[0]); if (post_cracks.size() > 2 || post_cracks.size() == 0) { + std::cout << post_cracks.size() << "\n"; throw badcycleex; } for (auto c : post_cracks) { @@ -260,15 +271,23 @@ void ma::post_fracture(network &n) { current_info ct = n.get_current_info(); - unsigned conducting_backbone_size = 0; + + std::vector<bool> vertex_in(n.G.vertices.size()); for (unsigned i = 0; i < n.G.edges.size(); i++) { if (ct.currents[i] > 1.0 / n.G.edges.size()) { - conducting_backbone_size++; + vertex_in[n.G.edges[i].v[0]] = true; + vertex_in[n.G.edges[i].v[1]] = true; } } - sb[conducting_backbone_size - 1]++; + unsigned bb_size = 0; + + for (unsigned i = 0; i < n.G.vertices.size(); i++) { + if (vertex_in[i]) bb_size++; + } + + sb[bb_size - 1]++; av_it++; @@ -276,5 +295,9 @@ void ma::post_fracture(network &n) { sa[(*av_it).size() - 1]++; av_it++; } + + sA[avalanches.back().size() - 1]++; + + sd[num - 1]++; } diff --git a/src/measurements.hpp b/src/measurements.hpp index 7c9d49c..51f7080 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -19,40 +19,20 @@ class ma : public hooks { // - interface for turning on and off specific measurements // private: - /* - double *fftw_forward_in; - fftw_complex *fftw_forward_out; - fftw_complex *fftw_reverse_in; - double *fftw_reverse_out; - fftw_plan forward_plan; - fftw_plan reverse_plan; - */ - unsigned N; - unsigned Mx; - unsigned My; Graph G; -// std::ofstream stress_file; + + unsigned num; // measurement storage - std::vector<uint64_t> sc; // non-spanning cluster size distribution std::vector<uint64_t> sn; // non-spanning cluster size distribution std::vector<uint64_t> ss; // minimal spanning cluster size distribution std::vector<uint64_t> sm; // spanning cluster size distribution - std::vector<uint64_t> sa; // non-final avalanche size distribution std::vector<uint64_t> sl; // final avalanche size distribution std::vector<uint64_t> sb; // final avalanche size distribution - std::vector<uint64_t> sD; // post-fracture damage distribution - std::vector<uint64_t> ccc; // cluster-cluster correlations - std::vector<uint64_t> css; // surface-surface correlations - std::vector<uint64_t> cmm; // surface-surface correlations - std::vector<uint64_t> caa; // avalanche-avalanche correlations - std::vector<uint64_t> cll; // damage-damage distribution - std::vector<uint64_t> cbb; // damage-damage distribution - std::vector<uint64_t> cDD; // damage-damage distribution - std::vector<uint64_t> csD; // damage-damage distribution - uint64_t Nc; - uint64_t Na; - uint64_t Nb; + std::vector<uint64_t> sd; // final avalanche size distribution + std::vector<uint64_t> sc; // non-spanning cluster size distribution + std::vector<uint64_t> sa; // non-final avalanche size distribution + std::vector<uint64_t> sA; // non-final avalanche size distribution public: long double lv; diff --git a/src/perc_meas.cpp b/src/perc_meas.cpp index 8270e1f..23b46b3 100644 --- a/src/perc_meas.cpp +++ b/src/perc_meas.cpp @@ -75,32 +75,19 @@ pm::pm(unsigned n, double a) : rank(2 * n), parent(2 * n), ds(&rank[0], &parent[0]), - sm(2 * n), - sl(2 * n), sn(3 * n), - sN(3 * n), ss(2 * n), + sm(2 * n), + sl(2 * n), sb(3 * n), - sB(3 * n), - sp(3 * n), - sr(3 * n), - sf(3 * n) + sd(3 * n) { model_string = "percolation_" + std::to_string(n) + "_" + std::to_string(a) + "_"; for (std::vector<uint64_t> &x : sn) { x.resize(2 * n); } for (std::vector<uint64_t> &x : sb) { - x.resize(3 * n); - } - for (std::vector<uint64_t> &x : sN) { - x.resize(2 * n); - } - for (std::vector<uint64_t> &x : sB) { - x.resize(3 * n); - } - for (std::vector<uint64_t> &x : sf) { - x.resize(3 * n); + x.resize(n); } } @@ -109,100 +96,73 @@ pm::pm(unsigned Lx, unsigned Ly) : rank(Lx * Ly / 2), parent(Lx * Ly / 2), ds(&rank[0], &parent[0]), - sm(Lx * Ly / 2), - sl(Lx * Ly / 2), sn(Lx * Ly), - sN(Lx * Ly), ss(Lx * Ly / 2), + sm(Lx * Ly / 2), + sl(Lx * Ly / 2), sb(Lx * Ly), - sB(Lx * Ly), - sp(Lx * Ly), - sr(Lx * Ly), - sf(Lx * Ly) + sd(Lx * Ly) { model_string = "percolation_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_"; for (std::vector<uint64_t> &x : sn) { x.resize(Lx * Ly / 2); } for (std::vector<uint64_t> &x : sb) { - x.resize(Lx * Ly); - } - for (std::vector<uint64_t> &x : sN) { x.resize(Lx * Ly / 2); } - for (std::vector<uint64_t> &x : sB) { - x.resize(Lx * Ly); - } - for (std::vector<uint64_t> &x : sf) { - x.resize(Lx * Ly); - } } pm::~pm() { - update_distribution_file("sm", sm, model_string); - update_distribution_file("sl", sl, model_string); update_distribution_file("sn", sn, 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("sB", sB, model_string); - update_distribution_file("sp", sp, model_string); - update_distribution_file("sr", sr, model_string); - update_distribution_file("sf", sf, model_string); + update_distribution_file("sd", sd, model_string); } void pm::pre_fracture(const network&) { boost::remove_edge_if(trivial, G); initialize_incremental_components(G, ds); incremental_components(G, ds); - p = 0; r = 0; - last_thres = std::numeric_limits<long double>::lowest(); } void pm::bond_broken(const network& net, const current_info& cur, unsigned i) { - unsigned p_old = p; - for (unsigned j = 0; j < net.G.edges.size(); j++) { - if (!net.fuses[j] && net.thresholds[j] < net.thresholds[i] && net.thresholds[j] > last_thres) { - p++; - } - } - last_thres = net.thresholds[i]; boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G); ds.union_set(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1]); boost::component_index<VertexIndex> components(parent.begin(), parent.end()); std::vector<unsigned> counts(components.size()); - for (unsigned j = 0; j < components.size(); j++) { + BOOST_FOREACH(VertexIndex current_index, components) { unsigned comp_size = 0; - BOOST_FOREACH(VertexIndex child_index, components[j]) { + BOOST_FOREACH(VertexIndex child_index, components[current_index]) { comp_size++; } sn[r][comp_size - 1]++; - for (unsigned k = p_old; k <= p; k++) { - sN[k][comp_size - 1]++; - } } - unsigned bb_size = 0; + std::vector<bool> vertex_in(net.G.vertices.size()); - for (unsigned j = 0; j < net.G.edges.size(); j++) { - if (cur.currents[j] > 1.0 / net.G.edges.size()) { - bb_size++; + for (unsigned i = 0; i < net.G.edges.size(); i++) { + if (cur.currents[i] > 1.0 / pow(net.G.edges.size(), 2)) { + vertex_in[net.G.edges[i].v[0]] = true; + vertex_in[net.G.edges[i].v[1]] = true; } } - sb[r][bb_size - 1]++; - for (unsigned k = p_old; k <= p; k++) { - sB[k][bb_size - 1]++; + unsigned bb_size = 0; + + for (unsigned i = 0; i < net.G.vertices.size(); i++) { + if (vertex_in[i]) bb_size++; } - sf[r][p]++; + sb[r][bb_size - 1]++; - p++; r++; + last_cur = cur; } void pm::post_fracture(network &n) { @@ -232,12 +192,28 @@ void pm::post_fracture(network &n) { for (unsigned j = r; j < sn.size(); j++) { sn[j][components[i].size() - 1]++; } - for (unsigned j = p; j < sn.size(); j++) { - sN[j][components[i].size() - 1]++; + } + + + std::vector<bool> vertex_in(n.G.vertices.size()); + + for (unsigned i = 0; i < n.G.edges.size(); i++) { + if (last_cur.currents[i] > 1.0 / n.G.edges.size()) { + vertex_in[n.G.edges[i].v[0]] = true; + vertex_in[n.G.edges[i].v[1]] = true; } } - sp[p - 1]++; - sr[r - 1]++; + unsigned bb_size = 0; + + for (unsigned i = 0; i < n.G.vertices.size(); i++) { + if (vertex_in[i]) bb_size++; + } + + for (unsigned i = r; i < sb.size(); i++) { + sb[i][bb_size - 1]++; + } + + sd[r - 1]++; } diff --git a/src/perc_meas.hpp b/src/perc_meas.hpp index e310189..e2357db 100644 --- a/src/perc_meas.hpp +++ b/src/perc_meas.hpp @@ -23,26 +23,22 @@ typedef boost::graph_traits<Graph>::vertices_size_type VertexIndex; class pm : public hooks { private: - unsigned p; unsigned r; - long double last_thres; + current_info last_cur; Graph G; std::vector<VertexIndex> rank; std::vector<Vertex> parent; boost::disjoint_sets<VertexIndex*, Vertex*> ds; // measurement storage - std::vector<uint64_t> sc; // non-spanning cluster size distribution std::vector<std::vector<uint64_t>> sn; // non-spanning cluster size distribution - std::vector<std::vector<uint64_t>> sN; // non-spanning cluster size distribution std::vector<uint64_t> ss; // minimal spanning cluster size distribution std::vector<uint64_t> sm; // spanning cluster size distribution std::vector<uint64_t> sl; // final avalanche size distribution std::vector<std::vector<uint64_t>> sb; // final avalanche size distribution - std::vector<std::vector<uint64_t>> sB; // final avalanche size distribution - std::vector<uint64_t> sp; + std::vector<uint64_t> sd; std::vector<uint64_t> sr; - std::vector<std::vector<uint64_t>> sf; // final avalanche size distribution + std::vector<unsigned> sb_tmp; public: std::string model_string; diff --git a/src/percolation.cpp b/src/percolation.cpp index 383d5e6..76731f4 100644 --- a/src/percolation.cpp +++ b/src/percolation.cpp @@ -33,7 +33,7 @@ int main(int argc, char* argv[]) { unsigned N = 1; unsigned Lx = 16; unsigned Ly = 16; - double beta = 0.0001; + double beta = 0.1; unsigned n = 128; double a = 1.0; diff --git a/src/sample.cpp b/src/sample.cpp index 488d20d..febff19 100644 --- a/src/sample.cpp +++ b/src/sample.cpp @@ -58,8 +58,8 @@ void sample::pre_fracture(const network& n) { } void sample::bond_broken(const network& net, const current_info& cur, unsigned i) { - long double c = logl(cur.conductivity / fabs(cur.currents[i])) + net.thresholds[i]; - sample_file << "{" << i << "," << c << "," << cur.conductivity << "},"; + long double c = logl(cur.conductivity[0] / fabs(cur.currents[i])) + net.thresholds[i]; + sample_file << "{" << i << "," << c << "," << cur.conductivity[0] << "},"; } void sample::post_fracture(network &n) { diff --git a/src/sample_fracture.cpp b/src/sample_fracture.cpp index 782b3e4..417e7ec 100644 --- a/src/sample_fracture.cpp +++ b/src/sample_fracture.cpp @@ -20,7 +20,7 @@ int main(int argc, char* argv[]) { int opt; unsigned N = 1; - double Lx = 16; + unsigned Lx = 16; double Ly = 16; double beta = 0.5; @@ -53,10 +53,10 @@ int main(int argc, char* argv[]) { std::mt19937 rng{seeds}; for (unsigned trial = 0; trial < N; trial++) { - graph G(Lx, Ly, rng); - network network(G, &c); + graph G(Lx, 1.0, rng); + elastic_network network(G, &c); network.set_thresholds(beta, rng); - network.fracture(meas); + network.fracture(meas, true); /*graph G2 = G.rotate(); class network network2(G2, &c); |