From 21203053a99c1bc6ce5d502131dc9d964d30d842 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sun, 23 Dec 2018 11:33:02 -0500 Subject: cleaned up graph a bit, redefined assert so that certain run errors in jc_voronoi can be caught --- lib/src/graph.cpp | 39 +++++++++++++++++++++++++++------------ src/fracture.cpp | 9 +++++++-- src/measurements.cpp | 20 ++++++++++++++++---- src/measurements.hpp | 4 ++++ 4 files changed, 54 insertions(+), 18 deletions(-) diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp index 3d5647f..b9d4f03 100644 --- a/lib/src/graph.cpp +++ b/lib/src/graph.cpp @@ -6,8 +6,12 @@ #define JCV_REAL_TYPE double #define JCV_ATAN2 atan2 #define JCV_FLT_MAX 1.7976931348623157E+308 +#undef assert +#define assert_error_report_helper(cond) "assertion failed: " +#define assert(cond) {if(!(cond)) { std::cerr << assert_error_report_helper(cond) "\n"; throw assert_error_report_helper(cond); } } #include +// actual mod for floats double mod(double a, double b) { if (a >= 0) { return fmod(a, b); @@ -95,26 +99,30 @@ class eulerException: public std::exception graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) { 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 d(0.0, 1.0); unsigned int N = round(Lx * Ly / 2 + d(rng) - 0.5); unsigned int nv = N; - vertices.resize(nv); // the coordinates of the lattice, from which a delaunay triangulation - // and corresponding voronoi tessellation will be built + // and corresponding voronoi tessellation will be built. Everyone is in the + // rectangle (0, 0) (Lx, Ly) for (unsigned int i = 0; i < nv; i++) { vertices[i] = {{L.x * d(rng), L.y * d(rng)}}; } - double max_difference = std::numeric_limits::max(); + // set up the voronoi objects jcv_diagram diagram; memset(&diagram, 0, sizeof(jcv_diagram)); + jcv_rect bounds = {{-Lx, -Ly}, {2 * Lx, 2 * Ly}}; std::vector points(9 * nv); double rstep = sqrt(step); + double max_difference = std::numeric_limits::max(); while (max_difference > pow(relax, 2) * 2 * N) { double cur_diff = 0; // to make the resulting tessellation periodic, we tile eight (!) copies of @@ -123,6 +131,7 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) // translate back) for (unsigned int i = 0; i < nv; i++) { const vertex& v = vertices[i]; + points[9 * i + 0] = {v.r.x + 0.0, v.r.y + 0.0}; points[9 * i + 1] = {v.r.x + L.x, v.r.y + 0.0}; points[9 * i + 2] = {v.r.x - L.x, v.r.y + 0.0}; @@ -134,27 +143,33 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) points[9 * i + 8] = {v.r.x - L.x, v.r.y - L.y}; } - jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram); + // run voronoi + jcv_diagram_generate(9 * nv, points.data(), &bounds, &diagram); + + // relax the sites by moving the vertices towards their centroids const jcv_site* sites = jcv_diagram_get_sites(&diagram); - for (int i = 0; i < diagram.numsites; i++) { const jcv_site* site = &sites[i]; unsigned int ind = site->index; if (ind % 9 == 0) { double Cx = 0.0; double Cy = 0.0; - unsigned int ne = 0; + double A = 0.0; const jcv_graphedge* e = site->edges; while (e) { - ne++; - Cx += e->pos[0].x; - Cy += e->pos[0].y; + double dA = e->pos[0].x * e->pos[1].y - e->pos[1].x * e->pos[0].y; + + A += dA; + Cx += (e->pos[0].x + e->pos[1].x) * dA; + Cy += (e->pos[0].y + e->pos[1].y) * dA; + e = e->next; } - Cx /= ne; - Cy /= ne; + A /= 2; + Cx /= 6 * A; + Cy /= 6 * A; double dx = Cx - vertices[ind / 9].r.x; double dy = Cy - vertices[ind / 9].r.y; @@ -183,7 +198,7 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) points[9 * i + 7] = {v.r.x + L.x, v.r.y - L.y}; points[9 * i + 8] = {v.r.x - L.x, v.r.y - L.y}; } - jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram); + jcv_diagram_generate(9 * nv, points.data(), &bounds, &diagram); const jcv_site* sites = jcv_diagram_get_sites(&diagram); diff --git a/src/fracture.cpp b/src/fracture.cpp index 742a3f9..47b4856 100644 --- a/src/fracture.cpp +++ b/src/fracture.cpp @@ -34,8 +34,9 @@ int main(int argc, char* argv[]) { double Lx = 16; double Ly = 16; double beta = 0.5; + double w = 0.01; - while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) { + while ((opt = getopt(argc, argv, "N:X:Y:b:w:")) != -1) { switch (opt) { case 'N': N = (unsigned int)atof(optarg); @@ -49,6 +50,9 @@ int main(int argc, char* argv[]) { case 'b': beta = atof(optarg); break; + case 'w': + w = atof(optarg); + break; default: exit(1); } @@ -65,12 +69,13 @@ int main(int argc, char* argv[]) { for (unsigned int trial = 0; trial < N; trial++) { while (true) { try { - graph G(Lx, Ly, rng); + graph G(Lx, Ly, rng, w); network network(G, &c); network.set_thresholds(beta, rng); network.fracture(meas); break; } catch (std::exception &e) { + std::cout << e.what() << '\n'; } } diff --git a/src/measurements.cpp b/src/measurements.cpp index 27f239d..15d37d1 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -112,6 +112,8 @@ 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), + 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), @@ -125,6 +127,8 @@ ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta) : N = 0; Nc = 0; Na = 0; + NC = 0; + NA = 0; // FFTW setup for correlation functions fftw_set_timelimit(1); @@ -150,6 +154,8 @@ ma::~ma() { update_distribution_file("sa", sa, Na, Lx, Ly, beta); update_distribution_file("sc", sc, Nc, Lx, Ly, 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); @@ -172,7 +178,9 @@ void ma::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]; if (c > lv && avalanches.back().size() > 0) { sa[avalanches.back().size() - 1]++; + sA[avalanches.back().size() - 1]++; Na++; + NA++; std::fill_n(fftw_forward_in, Mx * My, 0.0); @@ -224,8 +232,10 @@ void ma::post_fracture(network &n) { } sc[components[i].size() - 1]++; + sC[components[i].size() - 1]++; autocorrelation(Mx, My, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); Nc++; + NC++; 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; @@ -252,10 +262,10 @@ void ma::post_fracture(network &n) { // spanning cluster std::fill_n(fftw_forward_in, Mx * My, 0.0); - for (unsigned int i = 0; i < n.G.dual_vertices.size(); i++) { - if (component[i] == crack_component) { - fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[i].r, Lx, Ly, Mx, My)] += 1.0; - } + sC[components[crack_component].size() - 1]++; + NC++; + 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); @@ -277,6 +287,8 @@ void ma::post_fracture(network &n) { std::fill_n(fftw_forward_in, Mx * My, 0.0); // rewind the last avalanche + sA[avalanches.back().size() - 1]++; + NA++; 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); diff --git a/src/measurements.hpp b/src/measurements.hpp index bd93cfc..3d4ff3e 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -36,6 +36,8 @@ class ma : public hooks { // measurement storage std::vector sc; // cluster size distribution std::vector sa; // avalanche size distribution + std::vector sC; // cluster size distribution + std::vector sA; // avalanche size distribution std::vector sd; // avalanche size distribution std::vector sb; // bin size distribution std::vector Ccc; // cluster-cluster correlations @@ -46,7 +48,9 @@ class ma : public hooks { std::vector Cll; // damage-damage distribution std::vector Cee; // damage-damage distribution uint64_t Nc; + uint64_t NC; uint64_t Na; + uint64_t NA; public: long double lv; -- cgit v1.2.3-70-g09d2