diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2019-05-17 16:38:24 -0400 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2019-05-17 16:38:24 -0400 |
commit | 1e939e597964fa081b347e40af2be1069867b906 (patch) | |
tree | 395de78337d7ed85971fa3e84860ee4f4e5a5956 /src/perc_meas.cpp | |
parent | 456d748b8cbf141e1fe6e82a7c0106f469e28ba6 (diff) | |
download | fuse_networks-1e939e597964fa081b347e40af2be1069867b906.tar.gz fuse_networks-1e939e597964fa081b347e40af2be1069867b906.tar.bz2 fuse_networks-1e939e597964fa081b347e40af2be1069867b906.zip |
added support for running percolation using the fuse network rules
Diffstat (limited to 'src/perc_meas.cpp')
-rw-r--r-- | src/perc_meas.cpp | 167 |
1 files changed, 167 insertions, 0 deletions
diff --git a/src/perc_meas.cpp b/src/perc_meas.cpp new file mode 100644 index 0000000..4220a93 --- /dev/null +++ b/src/perc_meas.cpp @@ -0,0 +1,167 @@ + +#include "perc_meas.hpp" +#include <iostream> +#include <cstdio> + +class badcycleException: public std::exception +{ + virtual const char* what() const throw() + { + return "Could not find a valid cycle on the broken system."; + } +} badcycleex; + +void update_distribution_file(std::string id, const std::vector<uint64_t>& data, std::string model_string) { + std::string filename = model_string + id + ".dat"; + std::ifstream file(filename); + + std::vector<uint64_t> data_old(data.size(), 0); + + if (file.is_open()) { + for (unsigned i = 0; i < data.size(); i++) { + uint64_t num; + file >> num; + data_old[i] = num; + } + + file.close(); + } + + std::ofstream file_out(filename); + + for (unsigned i = 0; i < data.size(); i++) { + file_out <<std::fixed<< data_old[i] + data[i] << " "; + } + + file_out.close(); +} + +void update_distribution_file(std::string id, const std::vector<std::vector<uint64_t>>& data, std::string model_string) { + std::string filename = model_string + id + ".dat"; + std::ifstream file(filename); + + std::vector<std::vector<uint64_t>> data_old(data.size()); + + for (unsigned i = 0; i < data.size(); i++) { + data_old[i].resize(data[0].size()); + } + + if (file.is_open()) { + for (unsigned i = 0; i < data.size(); i++) { + for (unsigned j = 0; j < data[0].size(); j++) { + uint64_t num; + file >> num; + data_old[i][j] = num; + } + } + + file.close(); + } + + std::ofstream file_out(filename); + + for (unsigned i = 0; i < data.size(); i++) { + for (unsigned j = 0; j < data[0].size(); j++) { + file_out << std::fixed << data_old[i][j] + data[i][j] << " "; + } + file_out << "\n"; + } + + file_out.close(); +} + +pm::pm(unsigned n, double a) : + G(2 * n), + rank(2 * n), + parent(2 * n), + ds(&rank[0], &parent[0]), + sm(2 * n), + sl(2 * n), + sn(3 * n), + ss(2 * n), + sb(3 * n), + sp(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); + } +} + +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("ss", ss, model_string); + update_distribution_file("sb", sb, model_string); + update_distribution_file("sp", sp, 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; +} + +void pm::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); + 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 i = 0; i < components.size(); i++) { + unsigned comp_size = 0; + BOOST_FOREACH(VertexIndex child_index, components[i]) { + comp_size++; + } + + sn[p][comp_size - 1]++; + } + + unsigned bb_size = 0; + + for (unsigned i = 0; i <net.G.edges.size(); i++) { + if (cur.currents[i] > 1.0 / net.G.edges.size()) { + bb_size++; + } + } + + sb[p][bb_size - 1]++; + p++; +} + +void pm::post_fracture(network &n) { + auto post_cracks = find_minimal_crack(G, 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) { + throw badcycleex; + } + for (auto c : post_cracks) { + sl[c.second.size() - 1]++; + } + unsigned crack_component = component[n.G.dual_edges[post_cracks.front().second.front()].v[0]]; + + std::vector<std::list<graph::coordinate>> components(num); + + for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) { + components[component[i]].push_back(n.G.dual_vertices[i].r); + } + + for (unsigned i = 0; i < num; i++) { + if (i != crack_component) { + sm[components[i].size() - 1]++; + } else { + ss[components[i].size() - 1]++; + } + sn[p][components[i].size() - 1]++; + } + + sp[p - 1]++; +} + |