summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/analysis_tools.cpp74
-rw-r--r--src/analysis_tools.hpp2
-rw-r--r--src/animate.cpp32
-rw-r--r--src/animate_fracture.cpp6
-rw-r--r--src/animate_fracture_square.cpp4
-rw-r--r--src/fracture.cpp12
-rw-r--r--src/fracture_elastic.cpp87
-rw-r--r--src/fracture_square.cpp6
-rw-r--r--src/measurements.cpp160
-rw-r--r--src/measurements.hpp5
-rw-r--r--src/sample_fracture_square.cpp4
11 files changed, 202 insertions, 190 deletions
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<boost::undirected_tag,unsigned long>)
return true;
}
-std::list<unsigned> find_minimal_crack(const Graph& G, const network& n) {
+std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> find_minimal_crack(const Graph& G, const network& n) {
Graph Gtmp(n.G.vertices.size());
std::list<unsigned> removed_edges;
@@ -79,53 +79,53 @@ std::list<unsigned> find_minimal_crack(const Graph& G, const network& n) {
}
}
- if (cycles.size() > 1) {
- std::list<std::valarray<uint8_t>> bool_cycles;
- for (auto cycle : cycles) {
- std::valarray<uint8_t> bool_cycle(n.G.edges.size());
- for (auto v : cycle) {
- bool_cycle[v] = 1;
- }
- bool_cycles.push_back(bool_cycle);
+ std::list<std::valarray<uint8_t>> bool_cycles;
+ for (auto cycle : cycles) {
+ std::valarray<uint8_t> 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<uint8_t> new_bool_cycle = (*it1) ^ (*it2);
- std::list<unsigned> 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<uint8_t> new_bool_cycle = (*it1) ^ (*it2);
+ std::list<unsigned> 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<unsigned, 2> crossing_count{0,0};
+ std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> 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<unsigned, 2> 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<T> &, const std::list<T> &);
bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>);
-std::list<unsigned> find_minimal_crack(const Graph &, const network &);
+std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> 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 <iostream>
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<unsigned> 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<double> 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<unsigned> component(boost::num_vertices(G));
unsigned num = boost::connected_components(G, &component[0]);
- std::list<unsigned> crack = find_minimal_crack(G, n);
- unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
-
std::vector<std::list<unsigned>> 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<double> 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 <random>
+#include <iostream>
+
+#include <cholmod.h>
+
+#include "randutils/randutils.hpp"
+
+#include <graph.hpp>
+#include <network.hpp>
+#include <hooks.hpp>
+#include "measurements.hpp"
+
+#include <csignal>
+#include <cstring>
+#include <atomic>
+
+std::atomic<bool> 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 <iostream>
#include <cstdio>
-void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned N, std::string model_string) {
+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);
- uint64_t N_old = 0;
std::vector<uint64_t> 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<uint64_t>& data,
std::ofstream file_out(filename);
- file_out <<std::fixed<< N_old + N << "\n";
for (unsigned i = 0; i < data.size(); i++) {
file_out <<std::fixed<< data_old[i] + data[i] << " ";
}
@@ -179,122 +176,51 @@ ma::ma(unsigned n, double a, unsigned Mx, unsigned My, double beta) :
*/
}
-ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta) :
- Mx(Mx), My(My), G(2 * (unsigned)ceil(Lx * Ly / 2)),
- sc(3 * (unsigned)ceil(Lx * Ly / 2), 0),
- ss(3 * (unsigned)ceil(Lx * Ly / 2), 0),
- sm(3 * (unsigned)ceil(Lx * Ly / 2), 0),
- sa(3 * (unsigned)ceil(Lx * Ly / 2), 0),
- sl(3 * (unsigned)ceil(Lx * Ly / 2), 0),
- sb(3 * (unsigned)ceil(Lx * Ly / 2), 0),
- sD(3 * (unsigned)ceil(Lx * Ly / 2), 0),
- ccc((Mx / 2) * (My / 2), 0),
- css((Mx / 2) * (My / 2), 0),
- cmm((Mx / 2) * (My / 2), 0),
- caa((Mx / 2) * (My / 2), 0),
- cll((Mx / 2) * (My / 2), 0),
- cbb((Mx / 2) * (My / 2), 0),
- cDD((Mx / 2) * (My / 2), 0),
- csD((Mx / 2) * (My / 2), 0)
+ma::ma(unsigned Lx, unsigned Ly, double beta) :
+ G(Lx * Ly / 2),
+ sc(Lx * Ly / 2, 0),
+ sm(Lx * Ly / 2, 0),
+ sa(Lx * Ly, 0)
{
- N = 0;
- Nc = 0;
- Na = 0;
- Nb = 0;
-
if (beta != 0.0) {
model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_";
} else {
model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_INF_";
}
-
- // FFTW setup for correlation functions
- /*
- fftw_set_timelimit(1);
-
- fftw_forward_in = (double *)fftw_malloc(Mx * My * sizeof(double));
- fftw_forward_out = (fftw_complex *)fftw_malloc(Mx * My * sizeof(fftw_complex));
- fftw_reverse_in = (fftw_complex *)fftw_malloc(Mx * My * sizeof(fftw_complex));
- fftw_reverse_out = (double *)fftw_malloc(Mx * My * sizeof(double));
-
- forward_plan = fftw_plan_dft_r2c_2d(My, Mx, fftw_forward_in, fftw_forward_out, 0);
- reverse_plan = fftw_plan_dft_c2r_2d(My, Mx, fftw_reverse_in, fftw_reverse_out, 0);
- */
}
ma::~ma() {
- // clean up FFTW objects
- /*
- fftw_free(fftw_forward_in);
- fftw_free(fftw_forward_out);
- fftw_free(fftw_reverse_in);
- fftw_free(fftw_reverse_out);
- fftw_destroy_plan(forward_plan);
- fftw_destroy_plan(reverse_plan);
- fftw_cleanup();
- */
-
- update_distribution_file("sc", sc, Nc, model_string);
- update_distribution_file("ss", ss, N, model_string);
- update_distribution_file("sm", sm, N, model_string);
- update_distribution_file("sa", sa, Na, model_string);
- update_distribution_file("sl", sl, N, model_string);
- update_distribution_file("sb", sb, Nb, model_string);
- update_distribution_file("sD", sD, N, model_string);
-
- update_field_file("ccc", ccc, Nc, model_string, Mx, My);
- update_field_file("css", css, N, model_string, Mx, My);
- update_field_file("cmm", cmm, N, model_string, Mx, My);
- update_field_file("caa", caa, Na, model_string, Mx, My);
- update_field_file("cll", cll, N, model_string, Mx, My);
- update_field_file("cbb", cbb, Nb, model_string, Mx, My);
- update_field_file("cDD", cDD, N, model_string, Mx, My);
- update_field_file("csD", csD, N, model_string, Mx, My);
-
- //stress_file.close();
+ update_distribution_file("sc", sc, model_string);
+ update_distribution_file("sm", sm, model_string);
+ update_distribution_file("sa", sa, model_string);
}
void ma::pre_fracture(const network&) {
lv = std::numeric_limits<long double>::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<unsigned> component(boost::num_vertices(G));
unsigned num = boost::connected_components(G, &component[0]);
-
- std::list<unsigned> crack = find_minimal_crack(G, n);
-
- ss[crack.size() - 1]++;
- std::list<graph::coordinate> 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<std::list<graph::coordinate>> 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<graph::coordinate> 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<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
@@ -56,11 +57,11 @@ class ma : public hooks {
public:
long double lv;
- std::list<std::list<graph::coordinate>> avalanches;
+ std::list<std::list<unsigned>> avalanches;
std::list<unsigned> 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);
}