summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/fracture.cpp16
-rw-r--r--src/measurements.cpp118
-rw-r--r--src/measurements.hpp5
3 files changed, 72 insertions, 67 deletions
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 0aad56e..317d91b 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -31,16 +31,20 @@ int main(int argc, char* argv[]) {
int opt;
unsigned int N = 1;
- unsigned int L = 16;
+ unsigned int Lx = 16;
+ unsigned int Ly = 16;
double beta = 0.5;
- while ((opt = getopt(argc, argv, "N:L:b:")) != -1) {
+ while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
switch (opt) {
case 'N':
N = (unsigned int)atof(optarg);
break;
- case 'L':
- L = atoi(optarg);
+ case 'X':
+ Lx = atoi(optarg);
+ break;
+ case 'Y':
+ Ly = atoi(optarg);
break;
case 'b':
beta = atof(optarg);
@@ -53,9 +57,9 @@ int main(int argc, char* argv[]) {
cholmod_common c;
CHOL_F(start)(&c);
- graph G(L);
+ graph G(Lx, Ly);
network base_network(G, &c);
- ma meas(L, beta);
+ ma meas(Lx, Ly, beta);
randutils::auto_seed_128 seeds;
std::mt19937 rng{seeds};
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 57f85f5..e32f07f 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -104,11 +104,11 @@ std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
for (auto edge : cycle) {
double dx = fabs(n.G.vertices[n.G.dual_edges[edge][0]].r.x - n.G.vertices[n.G.dual_edges[edge][1]].r.x);
- if (dx > 0.5) {
+ if (dx > n.G.L.x / 2) {
crossing_count[0]++;
}
double dy = fabs(n.G.vertices[n.G.dual_edges[edge][0]].r.y - n.G.vertices[n.G.dual_edges[edge][1]].r.y);
- if (dy > 0.5) {
+ if (dy > n.G.L.y / 2) {
crossing_count[1]++;
}
}
@@ -124,8 +124,8 @@ std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
exit(5);
}
-void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, unsigned int L, double beta) {
- std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat";
+void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta) {
+ std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat";
std::ifstream file(filename);
uint64_t N_old = 0;
@@ -153,8 +153,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 int N, unsigned int L, double beta) {
- std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat";
+void update_field_file(std::string id, const std::vector<T>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta) {
+ std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat";
std::ifstream file(filename);
uint64_t N_old = 0;
@@ -179,16 +179,16 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned int
}
template <class T>
-std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<T>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) {
- for (unsigned int i = 0; i < pow(L, 2); i++) {
+std::vector<fftw_complex> data_transform(unsigned int Lx, unsigned int Ly, const std::vector<T>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) {
+ for (unsigned int i = 0; i < Lx * Ly; i++) {
fftw_forward_in[i] = (double)data[i];
}
fftw_execute(forward_plan);
- std::vector<fftw_complex> output(L * (L / 2 + 1));
+ std::vector<fftw_complex> output(Lx * (Ly / 2 + 1));
- for (unsigned int i = 0; i < L * (L / 2 + 1); i++) {
+ for (unsigned int i = 0; i < Lx * (Ly / 2 + 1); i++) {
output[i][0] = fftw_forward_out[i][0];
output[i][1] = fftw_forward_out[i][1];
}
@@ -197,47 +197,47 @@ std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<T>& d
}
template <class T>
-void correlation(unsigned int L, std::vector<T>& data, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
- for (unsigned int i = 0; i < L * (L / 2 + 1); i++) {
+void correlation(unsigned int Lx, unsigned int Ly, std::vector<T>& data, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
+ for (unsigned int i = 0; i < Lx * (Ly / 2 + 1); i++) {
fftw_reverse_in[i][0] = tx1[i][0] * tx2[i][0] + tx1[i][1] * tx2[i][1];
fftw_reverse_in[i][1] = tx1[i][0] * tx2[i][1] - tx1[i][1] * tx2[i][0];
}
fftw_execute(reverse_plan);
- for (unsigned int i = 0; i < pow(L / 2 + 1, 2); i++) {
- data[i] += (T)(fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2));
+ for (unsigned int i = 0; i < (Lx / 2 + 1) * (Ly / 2 + 1); i++) {
+ data[i] += (T)(fftw_reverse_out[Lx * (i / (Lx / 2 + 1)) + i % (Lx / 2 + 1)] / (Lx * Ly));
}
}
template <class T>
-void autocorrelation(unsigned int L, 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) {
+void autocorrelation(unsigned int Lx, unsigned int Ly, 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);
- for (unsigned int i = 0; i < L * (L / 2 + 1); i++) {
+ for (unsigned int i = 0; i < Ly * (Lx / 2 + 1); i++) {
fftw_reverse_in[i][0] = pow(fftw_forward_out[i][0], 2) + pow(fftw_forward_out[i][1], 2);
fftw_reverse_in[i][1] = 0.0;
}
fftw_execute(reverse_plan);
- for (unsigned int i = 0; i < pow(L / 2 + 1, 2); i++) {
- out_data[i] += (T)(fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2));
+ for (unsigned int i = 0; i < (Lx / 2 + 1) * (Ly / 2 + 1); i++) {
+ out_data[i] += (T)(fftw_reverse_out[Lx * (i / (Lx / 2 + 1)) + i % (Lx / 2 + 1)] / (Lx * Ly));
}
}
-ma::ma(unsigned int L, double beta) : L(L), beta(beta), G(2 * pow(L / 2, 2)),
- sc(pow(L, 2), 0),
- sa(pow(L, 2), 0),
- sd(pow(L, 2), 0),
- sb(log2(L) + 1, 0),
- Ccc(pow(L / 2 + 1, 2), 0),
- Css(pow(L / 2 + 1, 2), 0),
- Cmm(pow(L / 2 + 1, 2), 0),
- Caa(pow(L / 2 + 1, 2), 0),
- Cdd(pow(L / 2 + 1, 2), 0),
- Cll(pow(L / 2 + 1, 2), 0),
- Cee(pow(L / 2 + 1, 2), 0)
+ma::ma(unsigned int Lx, unsigned int Ly, double beta) : Lx(Lx), Ly(Ly), beta(beta), G(Lx * Ly / 2),
+ sc(Lx * Ly, 0),
+ sa(Lx * Ly, 0),
+ sd(Lx * Ly, 0),
+ sb(log2(Ly) + 1, 0),
+ Ccc((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+ Css((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+ Cmm((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+ Caa((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+ Cdd((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+ Cll((Lx / 2 + 1) * (Ly / 2 + 1), 0),
+ Cee((Lx / 2 + 1) * (Ly / 2 + 1), 0)
{
N = 0;
Nc = 0;
@@ -246,13 +246,13 @@ ma::ma(unsigned int L, double beta) : L(L), beta(beta), G(2 * pow(L / 2, 2)),
// FFTW setup for correlation functions
fftw_set_timelimit(1);
- fftw_forward_in = (double *)fftw_malloc(pow(L, 2) * sizeof(double));
- fftw_forward_out = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex));
- fftw_reverse_in = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex));
- fftw_reverse_out = (double *)fftw_malloc(pow(L, 2) * sizeof(double));
+ fftw_forward_in = (double *)fftw_malloc(Lx * Ly * sizeof(double));
+ fftw_forward_out = (fftw_complex *)fftw_malloc(Lx * Ly * sizeof(fftw_complex));
+ fftw_reverse_in = (fftw_complex *)fftw_malloc(Lx * Ly * sizeof(fftw_complex));
+ fftw_reverse_out = (double *)fftw_malloc(Lx * Ly * sizeof(double));
- forward_plan = fftw_plan_dft_r2c_2d(L, L, fftw_forward_in, fftw_forward_out, 0);
- reverse_plan = fftw_plan_dft_c2r_2d(L, L, fftw_reverse_in, fftw_reverse_out, 0);
+ forward_plan = fftw_plan_dft_r2c_2d(Ly, Lx, fftw_forward_in, fftw_forward_out, 0);
+ reverse_plan = fftw_plan_dft_c2r_2d(Ly, Lx, fftw_reverse_in, fftw_reverse_out, 0);
}
ma::~ma() {
@@ -265,18 +265,18 @@ ma::~ma() {
fftw_destroy_plan(reverse_plan);
fftw_cleanup();
- update_distribution_file("sa", sa, Na, L, beta);
- update_distribution_file("sc", sc, Nc, L, beta);
- update_distribution_file("sd", sd, N, L, beta);
- update_distribution_file("sb", sb, N, L, beta);
-
- update_field_file("Ccc", Ccc, Nc, L, beta);
- update_field_file("Css", Css, N, L, beta);
- update_field_file("Cmm", Cmm, N, L, beta);
- update_field_file("Cdd", Cdd, N, L, beta);
- update_field_file("Caa", Caa, Na, L, beta);
- update_field_file("Cll", Cll, N, L, beta);
- update_field_file("Cee", Cee, N, L, beta);
+ update_distribution_file("sa", sa, Na, Lx, Ly, beta);
+ update_distribution_file("sc", sc, Nc, Lx, Ly, beta);
+ update_distribution_file("sd", sd, N, Lx, Ly, beta);
+ update_distribution_file("sb", sb, N, Lx, Ly, beta);
+
+ update_field_file("Ccc", Ccc, Nc, Lx, Ly, beta);
+ update_field_file("Css", Css, N, Lx, Ly, beta);
+ update_field_file("Cmm", Cmm, N, Lx, Ly, beta);
+ update_field_file("Cdd", Cdd, N, Lx, Ly, beta);
+ update_field_file("Caa", Caa, Na, Lx, Ly, beta);
+ update_field_file("Cll", Cll, N, Lx, Ly, beta);
+ update_field_file("Cee", Cee, N, Lx, Ly, beta);
}
void ma::pre_fracture(const network &) {
@@ -297,7 +297,7 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned int i
fftw_forward_in[e] = 1.0;
}
- autocorrelation(L, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation(Lx, Ly, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
lv = c;
avalanches.push_back({i});
@@ -332,20 +332,20 @@ void ma::post_fracture(network &n) {
if (size > 0) {
sc[size - 1]++;
- autocorrelation(L, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation(Lx, Ly, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
Nc++;
}
}
}
// bin counting
- for (unsigned int be = 0; be < log2(L); be++) {
+ for (unsigned int be = 0; be < log2(Ly); be++) {
unsigned int bin = pow(2, be);
- for (unsigned int i = 0; i < pow(L / bin, 2); i++) {
+ for (unsigned int i = 0; i < Lx * Ly / pow(bin, 2); i++) {
bool in_bin = false;
for (unsigned int j = 0; j < pow(bin, 2); j++) {
- unsigned int edge = L * (bin * ((i * bin) / L) + j / bin) + (i * bin) % L + j % bin;
+ unsigned int edge = Lx * (bin * ((i * bin) / Lx) + j / bin) + (i * bin) % Lx + j % bin;
if (!n.fuses[edge] && crack_component == component[n.G.dual_edges[edge][0]]) {
in_bin = true;
break;
@@ -358,7 +358,7 @@ void ma::post_fracture(network &n) {
}
}
- sb[log2(L)]++;
+ sb[log2(Ly)]++;
for (unsigned int i = 0; i < n.G.edges.size(); i++) {
if (n.fuses[i] && component[n.G.dual_edges[i][0]] == crack_component) {
@@ -368,7 +368,7 @@ void ma::post_fracture(network &n) {
}
}
- autocorrelation(L, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation(Lx, Ly, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
// crack surface correlations
std::fill_n(fftw_forward_in, n.G.edges.size(), 0.0);
@@ -377,7 +377,7 @@ void ma::post_fracture(network &n) {
fftw_forward_in[edge] = 1.0;
}
- autocorrelation(L, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation(Lx, Ly, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
std::function<bool(unsigned int)> inCrack = [&](unsigned int i) -> bool {
return (bool)fftw_forward_in[i];
@@ -391,7 +391,7 @@ void ma::post_fracture(network &n) {
}
}
- autocorrelation(L, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation(Lx, Ly, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
std::fill_n(fftw_forward_in, n.G.edges.size(), 0.0);
@@ -402,7 +402,7 @@ void ma::post_fracture(network &n) {
fftw_forward_in[e] = 1;
}
- autocorrelation(L, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation(Lx, Ly, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
// damage size distribution
unsigned int total_broken = 0;
@@ -416,7 +416,7 @@ void ma::post_fracture(network &n) {
}
}
- autocorrelation(L, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation(Lx, Ly, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
sd[total_broken]++;
diff --git a/src/measurements.hpp b/src/measurements.hpp
index ab31ccd..703dccd 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -42,7 +42,8 @@ class ma : public hooks {
fftw_plan forward_plan;
fftw_plan reverse_plan;
unsigned int N;
- unsigned int L;
+ unsigned int Lx;
+ unsigned int Ly;
double beta;
Graph G;
@@ -67,7 +68,7 @@ class ma : public hooks {
std::list<std::list<unsigned int>> avalanches;
- ma(unsigned int L, double beta);
+ ma(unsigned int Lx, unsigned int Ly, double beta);
~ma();
void pre_fracture(const network &) override;