From f861195c1416220c4039bda4d1cbf8c3aab07528 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sun, 24 Feb 2019 17:19:52 -0500
Subject: added support for specifiying voronoi graphs with number of points
and aspect ratio, and added support for "infinite" beta (everything has the
some threshold)
---
lib/include/graph.hpp | 3 ++
lib/src/graph.cpp | 21 ++++++++---
lib/src/network.cpp | 17 ++++++---
src/fracture.cpp | 74 +++++++++++++++++++++++++++----------
src/measurements.cpp | 100 ++++++++++++++++++++++++++++++++++++--------------
src/measurements.hpp | 6 +--
6 files changed, 161 insertions(+), 60 deletions(-)
diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index 1089925..032a4d1 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -39,5 +39,8 @@ class graph {
graph(unsigned Nx, unsigned Ny);
graph(double Lx, double Ly, std::mt19937& rng);
+ graph(unsigned n, double a, std::mt19937& rng);
+
+ void helper(unsigned n, std::mt19937& rng);
};
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 526ffeb..8aa8877 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -141,27 +141,38 @@ unsigned get_triangle_signature(unsigned j1, unsigned j2, unsigned j3) {
}
graph::graph(double Lx, double Ly, std::mt19937& rng) {
- L = {Lx, Ly};
-
// randomly choose N to be floor(Lx * Ly / 2) or ceil(Lx * Ly / 2) with
// probability proportional to the distance from each
std::uniform_real_distribution<double> d(0.0, 1.0);
unsigned N = round(Lx * Ly / 2 + d(rng) - 0.5);
- unsigned nv = N;
+ L = {Lx, Ly};
+
+ this->helper(N, rng);
+}
+
+graph::graph(unsigned n, double a, std::mt19937& rng) {
+ L = {2 * n * a, 2 * n / a};
+
+ this->helper(n, rng);
+}
+
+void graph::helper(unsigned n, std::mt19937& rng) {
+ std::uniform_real_distribution<double> d(0.0, 1.0);
+ unsigned nv = n;
vertices.resize(nv);
// the coordinates of the lattice, from which a delaunay triangulation
// and corresponding voronoi tessellation will be built. Everyone is in the
// rectangle (0, 0) (Lx, Ly)
for (vertex &v : vertices) {
- v = {{Lx * d(rng), Ly * d(rng)}};
+ v = {{L.x * d(rng), L.y * d(rng)}};
}
// set up the voronoi objects
jcv_diagram diagram;
memset(&diagram, 0, sizeof(jcv_diagram));
- jcv_rect bounds = {{-Lx, -Ly}, {2 * Lx, 2 * Ly}};
+ jcv_rect bounds = {{-L.x, -L.y}, {2 * L.x, 2 * L.y}};
std::vector<jcv_point> points(9 * nv);
for (unsigned i = 0; i < nv; i++) {
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 027946a..43b0fab 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -95,13 +95,20 @@ network::~network() {
}
void network::set_thresholds(double beta, std::mt19937& rng) {
- std::uniform_real_distribution<long double> d(0.0, 1.0);
+ if (beta == 0.0) {
+ /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */
+ for (long double& threshold : thresholds) {
+ threshold = 1.0;
+ }
+ } else {
+ std::uniform_real_distribution<long double> d(0.0, 1.0);
- for (long double& threshold : thresholds) {
- threshold = std::numeric_limits<long double>::lowest();
+ for (long double& threshold : thresholds) {
+ threshold = std::numeric_limits<long double>::lowest();
- while (threshold == std::numeric_limits<long double>::lowest()) {
- threshold = logl(d(rng)) / (long double)beta;
+ while (threshold == std::numeric_limits<long double>::lowest()) {
+ threshold = logl(d(rng)) / (long double)beta;
+ }
}
}
}
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 53a11a5..b1bde45 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -35,7 +35,11 @@ int main(int argc, char* argv[]) {
double Ly = 16;
double beta = 0.5;
- while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
+ unsigned n = 128;
+ double a = 1.0;
+ bool use_aN = false;
+
+ while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:")) != -1) {
switch (opt) {
case 'N':
N = (unsigned)atof(optarg);
@@ -49,6 +53,14 @@ int main(int argc, char* argv[]) {
case 'b':
beta = atof(optarg);
break;
+ case 'n':
+ n = (unsigned)atof(optarg);
+ use_aN = true;
+ break;
+ case 'a':
+ a = atof(optarg);
+ use_aN = true;
+ break;
default:
exit(1);
}
@@ -57,30 +69,54 @@ int main(int argc, char* argv[]) {
cholmod_common c;
CHOL_F(start)(&c);
- // fourier transforms of powers of two are faster
- unsigned Mx = pow(2, ceil(log2(4*Lx)));
- unsigned My = pow(2, ceil(log2(4*Ly)));
-
- ma meas(Lx, Ly, Mx, My, beta);
-
randutils::auto_seed_128 seeds;
std::mt19937 rng{seeds};
- for (unsigned trial = 0; trial < N; trial++) {
- while (true) {
- try {
- graph G(Lx, Ly, rng);
- network network(G, &c);
- network.set_thresholds(beta, rng);
- network.fracture(meas);
- break;
- } catch (std::exception &e) {
- std::cout << e.what() << '\n';
+ if (use_aN) {
+ unsigned Mx = pow(2, ceil(log2(4*2*n*a)));
+ unsigned My = pow(2, ceil(log2(4*2*n/a)));
+
+ ma meas(n, a, Mx, My, beta);
+
+ for (unsigned trial = 0; trial < N; trial++) {
+ while (true) {
+ try {
+ graph G(n, a, rng);
+ network network(G, &c);
+ network.set_thresholds(beta, rng);
+ network.fracture(meas);
+ break;
+ } catch (std::exception &e) {
+ std::cout << e.what() << '\n';
+ }
}
+
+ if (quit.load())
+ break;
}
+ } else {
+ // fourier transforms of powers of two are faster
+ unsigned Mx = pow(2, ceil(log2(4*Lx)));
+ unsigned My = pow(2, ceil(log2(4*Ly)));
+
+ ma meas(Lx, Ly, Mx, My, beta);
+
+ for (unsigned trial = 0; trial < N; trial++) {
+ while (true) {
+ try {
+ graph G(Lx, Ly, rng);
+ network network(G, &c);
+ network.set_thresholds(beta, rng);
+ network.fracture(meas);
+ break;
+ } catch (std::exception &e) {
+ std::cout << e.what() << '\n';
+ }
+ }
- if (quit.load())
- break;
+ if (quit.load())
+ break;
+ }
}
CHOL_F(finish)(&c);
diff --git a/src/measurements.cpp b/src/measurements.cpp
index ed4b8eb..0b570cd 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -3,8 +3,8 @@
#include <iostream>
#include <cstdio>
-void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned N, double Lx, double Ly, double beta) {
- std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat";
+void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned N, std::string model_string) {
+ std::string filename = model_string + id + ".dat";
std::ifstream file(filename);
uint64_t N_old = 0;
@@ -32,8 +32,8 @@ void update_distribution_file(std::string id, const std::vector<uint64_t>& data,
}
template <class T>
-void update_field_file(std::string id, const std::vector<T>& data, unsigned N, double Lx, double Ly, double beta, unsigned Mx, unsigned My) {
- std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat";
+void update_field_file(std::string id, const std::vector<T>& data, unsigned N, std::string model_string, unsigned Mx, unsigned My) {
+ std::string filename = model_string + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat";
std::ifstream file(filename);
uint64_t N_old = 0;
@@ -95,7 +95,7 @@ void autocorrelation2(double Lx, double Ly, unsigned Mx, unsigned My, std::vecto
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)]++;
+ data[(unsigned)(Mx * (Δx / Lx)) + (Mx / 2) * (unsigned)(My * (Δy / Ly))]++;
}
}
}
@@ -136,8 +136,48 @@ 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(unsigned n, double a, unsigned Mx, unsigned My, double beta) :
+ Mx(Mx), My(My), G(2 * n),
+ sc(3 * n, 0),
+ ss(3 * n, 0),
+ sm(3 * n, 0),
+ sa(3 * n, 0),
+ sl(3 * n, 0),
+ sD(3 * n, 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),
+ cDD((Mx / 2) * (My / 2), 0),
+ csD((Mx / 2) * (My / 2), 0)
+{
+ N = 0;
+ Nc = 0;
+ Na = 0;
+
+ if (beta != 0.0) {
+ model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_";
+ } else {
+ model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_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(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)),
+ 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),
@@ -156,6 +196,12 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta) :
Nc = 0;
Na = 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);
@@ -182,20 +228,20 @@ ma::~ma() {
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_field_file("ccc", ccc, Nc, Lx, Ly, beta, Mx, My);
- update_field_file("css", css, N, Lx, Ly, beta, Mx, My);
- 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("csD", csD, N, Lx, Ly, beta, Mx, My);
+ 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("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("cDD", cDD, N, model_string, Mx, My);
+ update_field_file("csD", csD, N, model_string, Mx, My);
//stress_file.close();
}
@@ -212,7 +258,7 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
sa[avalanches.back().size() - 1]++;
Na++;
- autocorrelation2(Lx, Ly, Mx, My, caa, avalanches.back());
+ autocorrelation2(net.G.L.x, net.G.L.y, Mx, My, caa, avalanches.back());
lv = c;
avalanches.push_back({net.G.edges[i].r});
@@ -236,7 +282,7 @@ void ma::post_fracture(network &n) {
sr.push_back(n.G.dual_edges[edge].r);
}
- autocorrelation2(Lx, Ly, Mx, My, css, sr);
+ 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]];
@@ -251,7 +297,7 @@ void ma::post_fracture(network &n) {
sc[components[i].size() - 1]++;
Nc++;
- autocorrelation2(Lx, Ly, Mx, My, ccc, components[i]);
+ autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, ccc, components[i]);
}
}
@@ -259,7 +305,7 @@ void ma::post_fracture(network &n) {
sm[components[crack_component].size() - 1]++;
- autocorrelation2(Lx, Ly, Mx, My, cmm, components[crack_component]);
+ autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cmm, components[crack_component]);
/// damage at end
std::list<graph::coordinate> Dr;
@@ -272,15 +318,15 @@ void ma::post_fracture(network &n) {
sD[Dr.size() - 1]++;
- autocorrelation2(Lx, Ly, Mx, My, cDD, Dr);
- correlation2(Lx, Ly, Mx, My, csD, sr, Dr);
+ 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);
// ********************** LAST AVALANCHE *********************
// rewind the last avalanche
sl[avalanches.back().size() - 1]++;
- autocorrelation2(Lx, Ly, Mx, My, cll, avalanches.back());
+ autocorrelation2(n.G.L.x, n.G.L.y, Mx, My, cll, avalanches.back());
N++;
}
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 0fd4644..600c6ac 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -28,11 +28,8 @@ class ma : public hooks {
fftw_plan reverse_plan;
*/
unsigned N;
- double Lx;
- double Ly;
unsigned Mx;
unsigned My;
- double beta;
Graph G;
// std::ofstream stress_file;
@@ -56,10 +53,11 @@ class ma : public hooks {
public:
long double lv;
-
std::list<std::list<graph::coordinate>> avalanches;
+ std::string model_string;
ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta);
+ ma(unsigned n, double a, unsigned Mx, unsigned My, double beta);
~ma();
void pre_fracture(const network &) override;
--
cgit v1.2.3-70-g09d2