summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/measurements.cpp53
-rw-r--r--src/measurements.hpp32
-rw-r--r--src/perc_meas.cpp110
-rw-r--r--src/perc_meas.hpp10
-rw-r--r--src/percolation.cpp2
-rw-r--r--src/sample.cpp4
-rw-r--r--src/sample_fracture.cpp8
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);