From 4f4cf365eae07e04298459bf8f9e27ad0cfcc834 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 5 Dec 2018 19:16:28 -0500 Subject: now takes Lx and Ly for anisotropic fracture --- src/fracture.cpp | 16 ++++--- src/measurements.cpp | 118 +++++++++++++++++++++++++-------------------------- src/measurements.hpp | 5 ++- 3 files changed, 72 insertions(+), 67 deletions(-) (limited to 'src') 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 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 find_minimal_crack(const Graph& G, const network& n) { exit(5); } -void update_distribution_file(std::string id, const std::vector& 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& 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& data, } template -void update_field_file(std::string id, const std::vector& 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& 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& data, unsigned int } template -std::vector data_transform(unsigned int L, const std::vector& 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 data_transform(unsigned int Lx, unsigned int Ly, const std::vector& 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 output(L * (L / 2 + 1)); + std::vector 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 data_transform(unsigned int L, const std::vector& d } template -void correlation(unsigned int L, std::vector& data, const std::vector& tx1, const std::vector& 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& data, const std::vector& tx1, const std::vector& 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 -void autocorrelation(unsigned int L, std::vector& 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& 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 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> avalanches; - ma(unsigned int L, double beta); + ma(unsigned int Lx, unsigned int Ly, double beta); ~ma(); void pre_fracture(const network &) override; -- cgit v1.2.3-70-g09d2