From bad5c476749c99030767d1f67ae5cbe9698c99f5 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sat, 22 Dec 2018 14:30:54 -0500 Subject: system size now can take on non-integer values --- src/animate.cpp | 64 ++++++++++++++++++++++++------------------------ src/animate.hpp | 4 +-- src/animate_fracture.cpp | 16 +++++++----- src/fracture.cpp | 12 ++++----- src/measurements.cpp | 41 ++++++++++++++++--------------- src/measurements.hpp | 6 ++--- 6 files changed, 73 insertions(+), 70 deletions(-) (limited to 'src') diff --git a/src/animate.cpp b/src/animate.cpp index 8c2d12f..26f4996 100644 --- a/src/animate.cpp +++ b/src/animate.cpp @@ -1,15 +1,15 @@ #include "animate.hpp" -animate::animate(unsigned int Lx, unsigned int Ly, unsigned int window_size, int argc, char *argv[]) : Lx(Lx), Ly(Ly), G(Lx * Ly) { +animate::animate(double Lx, double Ly, unsigned int window_size, int argc, char *argv[]) : G(2 * (unsigned int)ceil(Lx * Ly / 2)) { glutInit(&argc, argv); glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); - glutInitWindowSize(window_size, (unsigned int )(window_size * ((double)Ly / (double)Lx))); + glutInitWindowSize((unsigned int)(Lx / Ly * window_size), window_size); glutCreateWindow("wolff"); glClearColor(0.0,0.0,0.0,0.0); glMatrixMode(GL_PROJECTION); glLoadIdentity(); - gluOrtho2D(-1.0, Lx + 3.0, -1.0 , Ly + 3.0); + gluOrtho2D(-1, Lx + 3, -1 , Ly + 3); } void animate::pre_fracture(const network &n) { @@ -25,18 +25,18 @@ void animate::pre_fracture(const network &n) { graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r; graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r; - if (fabs(r1.x - r2.x) > Lx / 2) { - if (r1.x < Lx / 2) { - r1.x += Lx; + if (fabs(r1.x - r2.x) > n.G.L.x / 2) { + if (r1.x < n.G.L.x / 2) { + r1.x += n.G.L.x; } else { - r2.x += Lx; + r2.x += n.G.L.x; } } - if (fabs(r1.y - r2.y) > Ly / 2) { - if (r1.y < Ly / 2) { - r1.y += Ly; + if (fabs(r1.y - r2.y) > n.G.L.y / 2) { + if (r1.y < n.G.L.y / 2) { + r1.y += n.G.L.y; } else { - r2.y += Ly; + r2.y += n.G.L.y; } } @@ -47,8 +47,8 @@ void animate::pre_fracture(const network &n) { glFlush(); } -void animate::bond_broken(const network& net, const current_info& cur, unsigned int i) { - long double c = logl(cur.conductivity / fabs(cur.currents[i])) + net.thresholds[i]; +void animate::bond_broken(const network& n, const current_info& cur, unsigned int i) { + long double c = logl(cur.conductivity / fabs(cur.currents[i])) + n.thresholds[i]; if (c > lv && avalanches.back().size() > 0) { lv = c; avalanches.push_back({i}); @@ -56,25 +56,25 @@ void animate::bond_broken(const network& net, const current_info& cur, unsigned avalanches.back().push_back(i); } - boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G); + boost::add_edge(n.G.dual_edges[i].v[0], n.G.dual_edges[i].v[1], {i}, G); glBegin(GL_LINES); // Each set of 3 vertices form a triangle glColor3f(1.0f, 1.0f, 1.0f); // Blue - graph::coordinate r1 = net.G.vertices[net.G.edges[i].v[0]].r; - graph::coordinate r2 = net.G.vertices[net.G.edges[i].v[1]].r; + graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r; + graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r; - if (fabs(r1.x - r2.x) > Lx / 2) { - if (r1.x < Lx / 2) { - r2.x -= Lx; + if (fabs(r1.x - r2.x) > n.G.L.x / 2) { + if (r1.x < n.G.L.x / 2) { + r2.x -= n.G.L.x; } else { - r2.x += Lx; + r2.x += n.G.L.x; } } - if (fabs(r1.y - r2.y) > Ly / 2) { - if (r1.y < Ly / 2) { - r2.y -= Ly; + if (fabs(r1.y - r2.y) > n.G.L.y / 2) { + if (r1.y < n.G.L.y / 2) { + r2.y -= n.G.L.y; } else { - r2.y += Ly; + r2.y += n.G.L.y; } } @@ -111,18 +111,18 @@ void animate::post_fracture(network &n) { graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r; graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r; - if (fabs(r1.x - r2.x) > Lx / 2) { - if (r1.x < Lx / 2) { - r1.x += Lx; + if (fabs(r1.x - r2.x) > n.G.L.x / 2) { + if (r1.x < n.G.L.x / 2) { + r1.x += n.G.L.x; } else { - r2.x += Lx; + r2.x += n.G.L.x; } } - if (fabs(r1.y - r2.y) > Ly / 2) { - if (r1.y < Ly / 2) { - r1.y += Ly; + if (fabs(r1.y - r2.y) > n.G.L.y / 2) { + if (r1.y < n.G.L.y / 2) { + r1.y += n.G.L.y; } else { - r2.y += Ly; + r2.y += n.G.L.y; } } diff --git a/src/animate.hpp b/src/animate.hpp index 793771d..423eba7 100644 --- a/src/animate.hpp +++ b/src/animate.hpp @@ -9,14 +9,12 @@ class animate : public hooks { private: - unsigned int Lx; - unsigned int Ly; Graph G; public: long double lv; std::list> avalanches; - animate(unsigned int Lx, unsigned int Ly, unsigned int window_size, int argc, char *argv[]); + animate(double Lx, double Ly, unsigned int window_size, int argc, char *argv[]); void pre_fracture(const network &) override; void bond_broken(const network& net, const current_info& cur, unsigned int i) override; diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp index 10b63d0..62c5b1e 100644 --- a/src/animate_fracture.cpp +++ b/src/animate_fracture.cpp @@ -31,24 +31,28 @@ int main(int argc, char* argv[]) { int opt; unsigned int N = 1; - unsigned int Lx = 16; - unsigned int Ly = 16; + double Lx = 16.0; + double Ly = 16.0; double beta = 0.5; + double w = 0.01; - while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) { + while ((opt = getopt(argc, argv, "X:Y:N:b:w:")) != -1) { switch (opt) { case 'N': N = (unsigned int)atof(optarg); break; case 'X': - Lx = atoi(optarg); + Lx = atof(optarg); break; case 'Y': - Ly = atoi(optarg); + Ly = atof(optarg); break; case 'b': beta = atof(optarg); break; + case 'w': + w = atof(optarg); + break; default: exit(1); } @@ -63,7 +67,7 @@ int main(int argc, char* argv[]) { std::mt19937 rng{seeds}; for (unsigned int trial = 0; trial < N; trial++) { - graph G(Lx, Ly, rng, 0.4); + graph G(Lx, Ly, rng, w); network network(G, &c); network.set_thresholds(beta, rng); network.fracture(meas); diff --git a/src/fracture.cpp b/src/fracture.cpp index 4b6fb5b..742a3f9 100644 --- a/src/fracture.cpp +++ b/src/fracture.cpp @@ -31,8 +31,8 @@ int main(int argc, char* argv[]) { int opt; unsigned int N = 1; - unsigned int Lx = 16; - unsigned int Ly = 16; + double Lx = 16; + double Ly = 16; double beta = 0.5; while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) { @@ -41,10 +41,10 @@ int main(int argc, char* argv[]) { N = (unsigned int)atof(optarg); break; case 'X': - Lx = atoi(optarg); + Lx = atof(optarg); break; case 'Y': - Ly = atoi(optarg); + Ly = atof(optarg); break; case 'b': beta = atof(optarg); @@ -57,7 +57,7 @@ int main(int argc, char* argv[]) { cholmod_common c; CHOL_F(start)(&c); - ma meas(Lx, Ly, Lx, Ly, beta); + ma meas(Lx, Ly, ceil(Lx), ceil(Ly), beta); randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; @@ -65,7 +65,7 @@ int main(int argc, char* argv[]) { for (unsigned int trial = 0; trial < N; trial++) { while (true) { try { - graph G(Lx, Ly, rng, 0.4); + graph G(Lx, Ly, rng); network network(G, &c); network.set_thresholds(beta, rng); network.fracture(meas); diff --git a/src/measurements.cpp b/src/measurements.cpp index 43483be..27f239d 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -2,7 +2,7 @@ #include "measurements.hpp" #include -void update_distribution_file(std::string id, const std::vector& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta) { +void update_distribution_file(std::string id, const std::vector& data, unsigned int 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"; std::ifstream file(filename); @@ -31,7 +31,7 @@ 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 Lx, unsigned int Ly, double beta, unsigned int Mx, unsigned int My) { +void update_field_file(std::string id, const std::vector& data, unsigned int N, double Lx, double Ly, double beta, unsigned int Mx, unsigned int 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"; std::ifstream file(filename); @@ -57,16 +57,16 @@ void update_field_file(std::string id, const std::vector& data, unsigned int } template -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++) { +std::vector data_transform(unsigned int Mx, unsigned int My, const std::vector& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) { + for (unsigned int i = 0; i < Mx * My; i++) { fftw_forward_in[i] = (double)data[i]; } fftw_execute(forward_plan); - std::vector output(Lx * (Ly / 2 + 1)); + std::vector output(Mx * (My / 2 + 1)); - for (unsigned int i = 0; i < Lx * (Ly / 2 + 1); i++) { + for (unsigned int i = 0; i < Mx * (My / 2 + 1); i++) { output[i][0] = fftw_forward_out[i][0]; output[i][1] = fftw_forward_out[i][1]; } @@ -75,43 +75,44 @@ std::vector data_transform(unsigned int Lx, unsigned int Ly, const } template -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++) { +void correlation(unsigned int Mx, unsigned int My, 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 < Mx * (My / 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 < (Lx / 2 + 1) * (Ly / 2 + 1); i++) { - data[i] += (T)(fftw_reverse_out[Lx * (i / (Lx / 2 + 1)) + i % (Lx / 2 + 1)] / (Lx * Ly)); + for (unsigned int i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) { + data[i] += (T)(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My)); } } template -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) { +void autocorrelation(unsigned int Mx, unsigned int My, 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 < Ly * (Lx / 2 + 1); i++) { + for (unsigned int i = 0; i < My * (Mx / 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 < (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)); + for (unsigned int i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) { + out_data[i] += (T)(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My)); } } -unsigned int edge_r_to_ind(graph::coordinate r, unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My) { +unsigned int edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned int Mx, unsigned int My) { return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly); } -ma::ma(unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My, double beta) : Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(Lx * Ly), - sc(Lx * Ly, 0), - sa(Lx * Ly, 0), - sd(Lx * Ly, 0), +ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta) : + Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(2 * (unsigned int)ceil(Lx * Ly / 2)), + sc(2 * (unsigned int)ceil(Lx * Ly / 2), 0), + sa(2 * (unsigned int)ceil(Lx * Ly / 2), 0), + sd(2 * (unsigned int)ceil(Lx * Ly / 2), 0), sb(log2(Mx < My ? Mx : My) + 1, 0), Ccc((Mx / 2 + 1) * (My / 2 + 1), 0), Css((Mx / 2 + 1) * (My / 2 + 1), 0), @@ -217,7 +218,7 @@ void ma::post_fracture(network &n) { // non-spanning clusters std::fill_n(fftw_forward_in, Mx * My, 0.0); for (unsigned int i = 0; i < num; i++) { - if (i != crack_component) { + if (i != crack_component && components[i].size() > 0) { 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)] += 1.0; } diff --git a/src/measurements.hpp b/src/measurements.hpp index 685806b..bd93cfc 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -26,8 +26,8 @@ class ma : public hooks { fftw_plan forward_plan; fftw_plan reverse_plan; unsigned int N; - unsigned int Lx; - unsigned int Ly; + double Lx; + double Ly; unsigned int Mx; unsigned int My; double beta; @@ -54,7 +54,7 @@ class ma : public hooks { std::list> avalanches; - ma(unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My, double beta); + ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta); ~ma(); void pre_fracture(const network &) override; -- cgit v1.2.3-54-g00ecf