summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/CMakeLists.txt9
-rw-r--r--src/fracture.cpp6
-rw-r--r--src/fracture_square.cpp2
-rw-r--r--src/measurements.cpp173
-rw-r--r--src/measurements.hpp22
-rw-r--r--src/sample.cpp69
-rw-r--r--src/sample.hpp30
-rw-r--r--src/sample_fracture.cpp65
-rw-r--r--src/sample_fracture_square.cpp67
9 files changed, 316 insertions, 127 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 4d4099f..eecbf72 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -8,6 +8,9 @@ target_link_libraries(measurements frac analysis_tools)
add_library(animate animate.cpp)
target_link_libraries(animate frac analysis_tools GL GLU glut)
+add_library(sample sample.cpp)
+target_link_libraries(sample frac)
+
add_executable(fracture fracture.cpp)
target_link_libraries(fracture frac measurements fftw3 cholmod profiler GL GLU glut)
@@ -20,3 +23,9 @@ target_link_libraries(animate_fracture frac animate fftw3 cholmod profiler GL GL
add_executable(animate_fracture_square animate_fracture_square.cpp)
target_link_libraries(animate_fracture_square frac animate fftw3 cholmod profiler GL GLU glut)
+add_executable(sample_fracture sample_fracture.cpp)
+target_link_libraries(sample_fracture frac sample cholmod)
+
+add_executable(sample_fracture_square sample_fracture_square.cpp)
+target_link_libraries(sample_fracture_square frac sample cholmod)
+
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 595cf1d..ebacd28 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -61,10 +61,10 @@ int main(int argc, char* argv[]) {
unsigned Mx = pow(2, ceil(log2(4*Lx)));
unsigned My = pow(2, ceil(log2(4*Ly)));
- ma meas(Lx, Ly, Mx, My, beta, 4);
+ ma meas(Lx, Ly, Mx, My, beta);
- randutils::auto_seed_128 seeds;
- std::mt19937 rng{seeds};
+ //randutils::auto_seed_128 seeds;
+ std::mt19937 rng{683129629};
for (unsigned trial = 0; trial < N; trial++) {
while (true) {
diff --git a/src/fracture_square.cpp b/src/fracture_square.cpp
index 47d7d67..dec4a2c 100644
--- a/src/fracture_square.cpp
+++ b/src/fracture_square.cpp
@@ -57,7 +57,7 @@ int main(int argc, char* argv[]) {
cholmod_common c;
CHOL_F(start)(&c);
- ma meas(Lx, Ly, 2*Lx, 2*Ly, beta, 10);
+ ma meas(Lx, Ly, 2*Lx, 2*Ly, beta);
graph G(Lx, Ly);
network perm_network(G, &c);
diff --git a/src/measurements.cpp b/src/measurements.cpp
index ff217b2..ed4b8eb 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -37,17 +37,12 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned N, d
std::ifstream file(filename);
uint64_t N_old = 0;
- std::vector<std::vector<uint64_t>> data_old(data.size());
- for (unsigned j = 0; j < data.size(); j++) {
- data_old[j].resize(data[j].size(), 0);
- }
+ std::vector<T> data_old(data.size(), 0);
if (file.is_open()) {
file >> N_old;
for (unsigned j = 0; j < data.size(); j++) {
- for (unsigned i = 0; i < data[j].size(); i++) {
- file >> data_old[j][i];
- }
+ file >> data_old[j];
}
file.close();
}
@@ -56,10 +51,7 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned N, d
file_out <<std::fixed<< N_old + N << "\n";
for (unsigned j = 0; j < data.size(); j++) {
- for (unsigned i = 0; i < data[j].size(); i++) {
- file_out << data_old[j][i] + data[j][i] << " ";
- }
- file_out << "\n";
+ file_out << data_old[j] + data[j] << " ";
}
file_out.close();
@@ -94,6 +86,34 @@ void correlation(unsigned Mx, unsigned My, std::vector<std::vector<T>>& data, co
}
}
+void autocorrelation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data, const std::list<graph::coordinate>& pos) {
+ for (std::list<graph::coordinate>::const_iterator it1 = pos.begin(); it1 != pos.end(); it1++) {
+ for (std::list<graph::coordinate>::const_iterator it2 = it1; it2 != pos.end(); it2++) {
+ double Δx_tmp = fabs(it1->x - it2->x);
+ double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp;
+
+ double Δy_tmp = fabs(it1->y - it2->y);
+ double Δy = Δy_tmp < Ly / 2 ? Δy_tmp : Ly - Δy_tmp;
+
+ data[floor((Mx * Δx) / Lx) + (Mx / 2) * floor((My * Δy) / Ly)]++;
+ }
+ }
+}
+
+void correlation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vector<uint64_t>& data, const std::list<graph::coordinate>& pos1, const std::list<graph::coordinate>& pos2) {
+ for (std::list<graph::coordinate>::const_iterator it1 = pos1.begin(); it1 != pos1.end(); it1++) {
+ for (std::list<graph::coordinate>::const_iterator it2 = pos2.begin(); it2 != pos2.end(); it2++) {
+ double Δx_tmp = fabs(it1->x - it2->x);
+ double Δx = Δx_tmp < Lx / 2 ? Δx_tmp : Lx - Δx_tmp;
+
+ double Δy_tmp = fabs(it1->y - it2->y);
+ double Δy = Δy_tmp < Ly / 2 ? Δy_tmp : Ly - Δy_tmp;
+
+ data[floor((Mx * Δx) / Lx) + (Mx / 2) * floor((My * Δy) / Ly)]++;
+ }
+ }
+}
+
template <class T>
void autocorrelation(unsigned Mx, unsigned My, std::vector<std::vector<T>>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
fftw_execute(forward_plan);
@@ -116,39 +136,28 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u
return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);
}
-ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncum) :
+ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta) :
Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), 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),
- sd(3 * (unsigned)ceil(Lx * Ly / 2), 0),
sD(3 * (unsigned)ceil(Lx * Ly / 2), 0),
- ccc(Ncum),
- css(Ncum),
- cmm(Ncum),
- caa(Ncum),
- cll(Ncum),
- cdd(Ncum),
- cDD(Ncum),
- csD(Ncum)
+ 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),
+ cDD((Mx / 2) * (My / 2), 0),
+ csD((Mx / 2) * (My / 2), 0)
{
N = 0;
Nc = 0;
Na = 0;
- for (unsigned i = 0; i < Ncum; i++) {
- ccc[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
- css[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
- cmm[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
- caa[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
- cll[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
- cdd[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
- cDD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
- csD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0);
- }
// FFTW setup for correlation functions
+ /*
fftw_set_timelimit(1);
fftw_forward_in = (double *)fftw_malloc(Mx * My * sizeof(double));
@@ -158,12 +167,12 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncu
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);
- std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_current.dat";
- //stress_file.open(filename, std::ifstream::app);
+ */
}
ma::~ma() {
// clean up FFTW objects
+ /*
fftw_free(fftw_forward_in);
fftw_free(fftw_forward_out);
fftw_free(fftw_reverse_in);
@@ -171,13 +180,13 @@ ma::~ma() {
fftw_destroy_plan(forward_plan);
fftw_destroy_plan(reverse_plan);
fftw_cleanup();
+ */
update_distribution_file("sc", sc, Nc, Lx, Ly, beta);
update_distribution_file("ss", ss, N, Lx, Ly, beta);
update_distribution_file("sm", sm, N, Lx, Ly, beta);
update_distribution_file("sa", sa, Na, Lx, Ly, beta);
update_distribution_file("sl", sl, N, Lx, Ly, beta);
- update_distribution_file("sd", sd, N, Lx, Ly, beta);
update_distribution_file("sD", sD, N, Lx, Ly, beta);
update_field_file("ccc", ccc, Nc, Lx, Ly, beta, Mx, My);
@@ -185,7 +194,6 @@ ma::~ma() {
update_field_file("cmm", cmm, N, Lx, Ly, beta, Mx, My);
update_field_file("caa", caa, Na, Lx, Ly, beta, Mx, My);
update_field_file("cll", cll, N, Lx, Ly, beta, Mx, My);
- update_field_file("cdd", cdd, N, Lx, Ly, beta, Mx, My);
update_field_file("cDD", cDD, N, Lx, Ly, beta, Mx, My);
update_field_file("csD", csD, N, Lx, Ly, beta, Mx, My);
@@ -204,19 +212,12 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
sa[avalanches.back().size() - 1]++;
Na++;
- std::fill_n(fftw_forward_in, Mx * My, 0.0);
-
- for (auto e : avalanches.back()) {
- unsigned ind = edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My);
- fftw_forward_in[ind] += 1.0;
- }
-
- autocorrelation(Mx, My, caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation2(Lx, Ly, Mx, My, caa, avalanches.back());
lv = c;
- avalanches.push_back({i});
+ avalanches.push_back({net.G.edges[i].r});
} else {
- avalanches.back().push_back(i);
+ avalanches.back().push_back(net.G.edges[i].r);
}
boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
@@ -228,110 +229,58 @@ void ma::post_fracture(network &n) {
std::list<unsigned> crack = find_minimal_crack(G, n);
- // crack surface correlations
- std::fill_n(fftw_forward_in, Mx * My, 0.0);
-
- ss[crack.size()]++;
+ ss[crack.size() - 1]++;
+ std::list<graph::coordinate> sr;
for (auto edge : crack) {
- fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My)] += 0.5;
- fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My)] += 0.5;
+ sr.push_back(n.G.dual_edges[edge].r);
}
- fftw_complex *tss = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);
-
- correlation(Mx, My, css, tss, tss, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation2(Lx, Ly, Mx, My, css, sr);
unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
- std::vector<std::list<unsigned>> components(num);
+ std::vector<std::list<graph::coordinate>> components(num);
for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
- components[component[i]].push_back(i);
+ components[component[i]].push_back(n.G.dual_vertices[i].r);
}
- // non-spanning clusters
- std::fill_n(fftw_forward_in, Mx * My, 0.0);
-
for (unsigned i = 0; i < num; i++) {
if (i != crack_component && components[i].size() > 0) {
sc[components[i].size() - 1]++;
Nc++;
- for (auto it = components[i].begin(); it != components[i].end(); it++) {
- unsigned ind = edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My);
- fftw_forward_in[ind] += 1.0;
- }
-
- autocorrelation(Mx, My, ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-
- for (auto it = components[i].begin(); it != components[i].end(); it++) {
- fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] = 0.0;
- }
+ autocorrelation2(Lx, Ly, Mx, My, ccc, components[i]);
}
}
// spanning cluster
- // std::fill_n(fftw_forward_in, Mx * My, 0.0); we already reset in the last loop
sm[components[crack_component].size() - 1]++;
- for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
- fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] += 1.0;
- }
-
- autocorrelation(Mx, My, cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation2(Lx, Ly, Mx, My, cmm, components[crack_component]);
/// damage at end
- std::fill_n(fftw_forward_in, Mx * My, 0.0);
-
- unsigned final_broken = 0;
+ std::list<graph::coordinate> Dr;
for (unsigned i = 0; i < n.G.edges.size(); i++) {
if (n.fuses[i]) {
- final_broken++;
- fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
+ Dr.push_back(n.G.edges[i].r);
}
}
- sD[final_broken]++;
-
- fftw_complex *tDD = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out);
-
- correlation(Mx, My, cDD, tDD, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
- correlation(Mx, My, csD, tss, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ sD[Dr.size() - 1]++;
- free(tss);
- free(tDD);
+ autocorrelation2(Lx, Ly, Mx, My, cDD, Dr);
+ correlation2(Lx, Ly, Mx, My, csD, sr, Dr);
- std::fill_n(fftw_forward_in, Mx * My, 0.0);
+ // ********************** LAST AVALANCHE *********************
// rewind the last avalanche
sl[avalanches.back().size() - 1]++;
- // rewind the last avalanche
- for (auto e : avalanches.back()) {
- boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
- n.break_edge(e, true);
- fftw_forward_in[edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0;
- }
-
- autocorrelation(Mx, My, cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-
- // damage size distribution
- unsigned total_broken = 0;
-
- std::fill_n(fftw_forward_in, Mx * My, 0.0);
-
- for (unsigned i = 0; i < n.G.edges.size(); i++) {
- if (n.fuses[i]) {
- fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0;
- }
- }
-
- autocorrelation(Mx, My, cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-
- sd[total_broken]++;
+ autocorrelation2(Lx, Ly, Mx, My, cll, avalanches.back());
N++;
}
diff --git a/src/measurements.hpp b/src/measurements.hpp
index e267064..0fd4644 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -19,12 +19,14 @@ 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;
double Lx;
double Ly;
@@ -40,16 +42,14 @@ class ma : public hooks {
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> sd; // pre-fracture damage distribution
std::vector<uint64_t> sD; // post-fracture damage distribution
- std::vector<std::vector<uint64_t>> ccc; // cluster-cluster correlations
- std::vector<std::vector<uint64_t>> css; // surface-surface correlations
- std::vector<std::vector<uint64_t>> cmm; // surface-surface correlations
- std::vector<std::vector<uint64_t>> caa; // avalanche-avalanche correlations
- std::vector<std::vector<uint64_t>> cll; // damage-damage distribution
- std::vector<std::vector<uint64_t>> cdd; // damage-damage distribution
- std::vector<std::vector<uint64_t>> cDD; // damage-damage distribution
- std::vector<std::vector<uint64_t>> csD; // damage-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> cDD; // damage-damage distribution
+ std::vector<uint64_t> csD; // damage-damage distribution
uint64_t Nc;
uint64_t Na;
@@ -57,9 +57,9 @@ class ma : public hooks {
long double lv;
- std::list<std::list<unsigned>> avalanches;
+ std::list<std::list<graph::coordinate>> avalanches;
- ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncum);
+ ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta);
~ma();
void pre_fracture(const network &) override;
diff --git a/src/sample.cpp b/src/sample.cpp
new file mode 100644
index 0000000..b8af937
--- /dev/null
+++ b/src/sample.cpp
@@ -0,0 +1,69 @@
+
+#include "sample.hpp"
+#include <iostream>
+#include <cstdio>
+
+sample::sample(double Lx, double Ly, double beta) {
+ std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_sample.dat";
+ sample_file.open(filename, std::ifstream::app);
+}
+
+sample::~sample() {
+ sample_file.close();
+}
+
+void sample::pre_fracture(const network& n) {
+ sample_file << "<|";
+ sample_file << "\"vr\"->{";
+ for (const graph::vertex &v : n.G.vertices) {
+ sample_file << "{" << v.r.x << "," << v.r.y << "},";
+ }
+ sample_file << "},\"vp\"->{";
+ for (const graph::vertex &v : n.G.vertices) {
+ sample_file << "{";
+ for (const graph::coordinate &r : v.polygon) {
+ sample_file << "{" << r.x << "," << r.y << "},";
+ }
+ sample_file << "},";
+ }
+ sample_file << "},\"ur\"->{";
+ for (const graph::vertex &v : n.G.dual_vertices) {
+ sample_file << "{" << v.r.x << "," << v.r.y << "},";
+ }
+ sample_file << "},\"up\"->{";
+ for (const graph::vertex &v : n.G.dual_vertices) {
+ sample_file << "{";
+ for (const graph::coordinate &r : v.polygon) {
+ sample_file << "{" << r.x << "," << r.y << "},";
+ }
+ sample_file << "},";
+ }
+ sample_file << "},\"e\"->{";
+ for (const graph::edge &e : n.G.edges) {
+ sample_file << "{" << e.v[0] << "," << e.v[1] << "},";
+ }
+ sample_file << "},\"ec\"->{";
+ for (const graph::edge &e : n.G.edges) {
+ sample_file << "{" << e.crossings[0] << "," << e.crossings[1] << "},";
+ }
+ sample_file << "},\"d\"->{";
+ for (const graph::edge &e : n.G.dual_edges) {
+ sample_file << "{" << e.v[0] << "," << e.v[1] << "},";
+ }
+ sample_file << "},\"dc\"->{";
+ for (const graph::edge &e : n.G.dual_edges) {
+ sample_file << "{" << e.crossings[0] << "," << e.crossings[1] << "},";
+ }
+ sample_file << "},\"data\"->{";
+}
+
+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 << "},";
+}
+
+void sample::post_fracture(network &n) {
+ sample_file << "}";
+ sample_file << "|>\n";
+}
+
diff --git a/src/sample.hpp b/src/sample.hpp
new file mode 100644
index 0000000..d5d6fe5
--- /dev/null
+++ b/src/sample.hpp
@@ -0,0 +1,30 @@
+
+#include <cstring>
+#include <fstream>
+#include <string>
+#include <cinttypes>
+#include <sstream>
+#include <functional>
+#include <iostream>
+#include <array>
+
+#include <hooks.hpp>
+#include <graph.hpp>
+#include <network.hpp>
+
+class sample : public hooks {
+ // need:
+ // - interface for turning on and off specific measurements
+ //
+ private:
+ std::ofstream sample_file;
+
+ public:
+ sample(double, double, double);
+ ~sample();
+
+ void pre_fracture(const network &) override;
+ void bond_broken(const network& net, const current_info& cur, unsigned i) override;
+ void post_fracture(network &n) override;
+};
+
diff --git a/src/sample_fracture.cpp b/src/sample_fracture.cpp
new file mode 100644
index 0000000..42a6a5a
--- /dev/null
+++ b/src/sample_fracture.cpp
@@ -0,0 +1,65 @@
+
+#include <random>
+#include <iostream>
+
+#include <cholmod.h>
+
+#include "randutils/randutils.hpp"
+
+#include <graph.hpp>
+#include <network.hpp>
+#include <hooks.hpp>
+#include "sample.hpp"
+
+#include <csignal>
+#include <cstring>
+#include <atomic>
+
+int main(int argc, char* argv[]) {
+
+ int opt;
+
+ unsigned N = 1;
+ double Lx = 16;
+ double 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 = atof(optarg);
+ break;
+ case 'Y':
+ Ly = atof(optarg);
+ break;
+ case 'b':
+ beta = atof(optarg);
+ break;
+ default:
+ exit(1);
+ }
+ }
+
+ cholmod_common c;
+ CHOL_F(start)(&c);
+
+ sample meas(Lx, Ly, beta);
+
+ randutils::auto_seed_128 seeds;
+ std::mt19937 rng{seeds};
+
+ for (unsigned trial = 0; trial < N; trial++) {
+ graph G(Lx, Ly, rng);
+ network network(G, &c);
+ network.set_thresholds(beta, rng);
+ network.fracture(meas);
+ }
+
+ CHOL_F(finish)(&c);
+
+ return 0;
+}
+
diff --git a/src/sample_fracture_square.cpp b/src/sample_fracture_square.cpp
new file mode 100644
index 0000000..fb2eb33
--- /dev/null
+++ b/src/sample_fracture_square.cpp
@@ -0,0 +1,67 @@
+
+#include <random>
+#include <iostream>
+
+#include <cholmod.h>
+
+#include "randutils/randutils.hpp"
+
+#include <graph.hpp>
+#include <network.hpp>
+#include <hooks.hpp>
+#include "sample.hpp"
+
+#include <csignal>
+#include <cstring>
+#include <atomic>
+
+int main(int argc, char* argv[]) {
+
+ 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);
+
+ sample meas(Lx, Ly, beta);
+
+ randutils::auto_seed_128 seeds;
+ std::mt19937 rng{seeds};
+
+ graph G(Lx, Ly);
+ network perm_network(G, &c);
+
+ for (unsigned trial = 0; trial < N; trial++) {
+ network tmp_network(perm_network);
+ tmp_network.set_thresholds(beta, rng);
+ tmp_network.fracture(meas);
+ }
+
+ CHOL_F(finish)(&c);
+
+ return 0;
+}
+