From 09200a607661f739782a966807d31345485e2c41 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 20 Dec 2018 12:20:19 -0500 Subject: added animation example, and did many fixes to the voronoi system --- CMakeLists.txt | 4 +- lib/include/graph.hpp | 11 +- lib/src/graph.cpp | 217 ++++++++++++++++++++++++--------- lib/src/network.cpp | 28 ++--- src/CMakeLists.txt | 15 ++- src/analysis_tools.cpp | 125 +++++++++++++++++++ src/analysis_tools.hpp | 29 +++++ src/animate.cpp | 205 ++++++++++++++++++++++++++++++++ src/animate.hpp | 24 ++++ src/animate_fracture.cpp | 79 ++++++++++++ src/fracture.cpp | 17 ++- src/measurements.cpp | 303 ++++++++++++++--------------------------------- src/measurements.hpp | 24 +--- src/test.cpp | 83 +++++++++++++ 14 files changed, 849 insertions(+), 315 deletions(-) create mode 100644 src/analysis_tools.cpp create mode 100644 src/analysis_tools.hpp create mode 100644 src/animate.cpp create mode 100644 src/animate.hpp create mode 100644 src/animate_fracture.cpp create mode 100644 src/test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index a804ac5..35e8436 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,8 +1,8 @@ cmake_minimum_required(VERSION 3.8) -set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall -stdlib=libc++") -set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -stdlib=libc++") +set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall -stdlib=libc++ -fuse-ld=lld") +set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto -Wall -stdlib=libc++ -march=native -fuse-ld=lld") set (CMAKE_CXX_STANDARD 17) set (CMAKE_C_STANDARD 11) diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp index 75d54f4..d746ddd 100644 --- a/lib/include/graph.hpp +++ b/lib/include/graph.hpp @@ -4,8 +4,11 @@ #include #include #include -#include +#include #include +#include +#include + class graph { public: @@ -16,9 +19,13 @@ class graph { typedef struct vertex { coordinate r; + std::vector polygon; } vertex; - typedef std::array edge; + typedef struct edge { + std::array v; + coordinate r; + } edge; coordinate L; diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp index f25cf87..10f62b7 100644 --- a/lib/src/graph.cpp +++ b/lib/src/graph.cpp @@ -1,7 +1,6 @@ #include "graph.hpp" #include -#include #define JC_VORONOI_IMPLEMENTATION #define JCV_REAL_TYPE double @@ -24,9 +23,21 @@ graph::graph(unsigned int Nx, unsigned int Ny) { for (unsigned int i = 0; i < nv; i++) { vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))); vertices[i].r.y = (double)(i / (Nx / 2)); + vertices[i].polygon = { + {vertices[i].r.x - 0.5, vertices[i].r.y}, + {vertices[i].r.x, vertices[i].r.y - 0.5}, + {vertices[i].r.x + 0.5, vertices[i].r.y}, + {vertices[i].r.x, vertices[i].r.y + 0.5} + }; dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))); dual_vertices[i].r.y = (double)(i / (Nx / 2)); + dual_vertices[i].polygon = { + {dual_vertices[i].r.x - 0.5, vertices[i].r.y}, + {dual_vertices[i].r.x, vertices[i].r.y - 0.5}, + {dual_vertices[i].r.x + 0.5, vertices[i].r.y}, + {dual_vertices[i].r.x, vertices[i].r.y + 0.5} + }; } for (unsigned int x = 0; x < Ny; x++) { @@ -34,40 +45,70 @@ graph::graph(unsigned int Nx, unsigned int Ny) { unsigned int v1 = (Nx * x) / 2 + ((y + x % 2) / 2) % (Nx / 2); unsigned int v2 = ((Nx * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2)) % nv; - edges.push_back({v1, v2}); + edges.push_back({{v1, v2}, {0.5 + (double)y, 0.5 + (double)x}}); unsigned int dv1 = (Nx * x) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2); unsigned int dv2 = ((Nx * (x + 1)) / 2 + ((y + x % 2) / 2) % (Nx / 2)) % nv; - dual_edges.push_back({dv1, dv2}); + dual_edges.push_back({{dv1, dv2}, {0.5 + (double)y, 0.5 + (double)x}}); } } } +namespace std +{ + template + struct hash > + { + typedef array argument_type; + typedef size_t result_type; + + result_type operator()(const argument_type& a) const + { + hash hasher; + result_type h = 0; + for (result_type i = 0; i < N; ++i) + { + h = h * 31 + hasher(a[i]); + } + return h; + } + }; +} + +class eulerException: public std::exception +{ + virtual const char* what() const throw() + { + return "The voronoi tessellation generated has the incorrect Euler characteristic for a torus and is malformed."; + } +} eulerex; + + graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) { L = {(double)Nx, (double)Ny}; - unsigned int dnv = Nx * Ny / 2; + unsigned int nv = Nx * Ny / 2; - dual_vertices.resize(dnv); + vertices.resize(nv); std::normal_distribution d(0.0, spread); - // the coordinates of the dual lattice, from which a delaunay triangulation - // and corresponding voronoi tesselation will be built - for (unsigned int i = 0; i < dnv; i++) { + // the coordinates of the lattice, from which a delaunay triangulation + // and corresponding voronoi tessellation will be built + for (unsigned int i = 0; i < nv; i++) { double rx = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))) + d(rng); double ry = (double)(i / (Nx / 2)) + d(rng); - dual_vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}}; + vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}}; } - // to make the resulting tesselation periodic, we tile eight (!) copies of - // the original dual points for a total of nine. note that the index of each + // to make the resulting tessellation periodic, we tile eight (!) copies of + // the original points for a total of nine. note that the index of each // point quotient 9 is equal to the original index (we'll use this to // translate back) std::vector points; - points.reserve(9 * dnv); - for (const vertex &v : dual_vertices) { + points.reserve(9 * nv); + for (const vertex &v : vertices) { points.push_back({v.r.x, v.r.y}); points.push_back({v.r.x + L.x, v.r.y + 0.0}); points.push_back({v.r.x - L.x, v.r.y + 0.0}); @@ -82,81 +123,151 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) jcv_diagram diagram; memset(&diagram, 0, sizeof(jcv_diagram)); - jcv_diagram_generate(9 * dnv, points.data(), NULL, &diagram); + jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram); const jcv_site* sites = jcv_diagram_get_sites(&diagram); - struct coorcomp { - bool operator() (const coordinate& lhs, const coordinate& rhs) const - {return lhs.x& lhs, const std::array& rhs) const + {if (lhs[0] != lhs[1]) return lhs[0] < lhs[1]; + else return rhs[0] < rhs[1]; + } }; - std::map known_vertices; + std::unordered_map, unsigned int> known_vertices; for (int i = 0; i < diagram.numsites; i++) { const jcv_site* site = &sites[i]; // we only care about processing the cells of our original, central sites if (site->index % 9 == 0) { + unsigned int i1 = (unsigned int)(site->index / 9); const jcv_graphedge* e = site->edges; + const jcv_graphedge* ep = site->edges; + while (ep->next) { + ep = ep->next; + } // for each edge on the site's cell while(e) { - // assess whether the dual edge needs to be added. only neighboring - // dual sites whose index is larger than ours are given edges, so we + // assess whether the edge needs to be added. only neighboring + // sites whose index is larger than ours are given edges, so we // don't double up later const jcv_site* neighbor = e->neighbor; - unsigned int index = (unsigned int)(site->index / 9); - unsigned int real_index = (unsigned int)(neighbor->index / 9); - - if (index < real_index) { - dual_edges.push_back({index, real_index}); - coordinate r1 = {1e-10 * round(1e10 * fmod(L.x + e->pos[0].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[0].y, L.y))}; - coordinate r2 = {1e-10 * round(1e10 * fmod(L.x + e->pos[1].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[1].y, L.y))}; + unsigned int i2 = (unsigned int)(neighbor->index / 9); + + vertices[i1].polygon.push_back({e->pos[0].x, e->pos[0].y}); + + unsigned int i3p = (unsigned int)(ep->neighbor->index / 9); + std::array t1 = {i1, i2, i3p}; + std::sort(t1.begin(), t1.end()); - std::map::iterator it1 = known_vertices.find(r1); - std::map::iterator it2 = known_vertices.find(r2); + auto it1 = known_vertices.find(t1); - unsigned int vi1, vi2; + unsigned int vi1; - if (it1 == known_vertices.end()) { - vi1 = vertices.size(); - vertices.push_back({r1}); - known_vertices[r1] = vi1; + if (it1 == known_vertices.end()) { + vi1 = dual_vertices.size(); + dual_vertices.push_back({{fmod(L.x + e->pos[0].x, L.x), fmod(L.y + e->pos[0].y, L.y)},{{site->p.x, site->p.y}}}); + known_vertices[t1] = vi1; + } else { + vi1 = it1->second; + dual_vertices[vi1].polygon.push_back({site->p.x, site->p.y}); + } + + if (i1 < i2) { + edges.push_back({{i1, i2}, + {fmod(L.x + (site->p.x + neighbor->p.x) / 2, L.x), + fmod(L.y + (site->p.y + neighbor->p.y) / 2, L.y)} + }); + + jcv_graphedge *en; + if (e->next == NULL) { + en = site->edges; } else { - vi1 = it1->second; + en = e->next; } + + unsigned int i3n = (unsigned int)(en->neighbor->index / 9); + std::array t2 = {i1, i2, i3n}; + std::sort(t2.begin(), t2.end()); + + auto it2 = known_vertices.find(t2); + + unsigned int vi2; + if (it2 == known_vertices.end()) { - vi2 = vertices.size(); - vertices.push_back({r2}); - known_vertices[r2] = vi2; + vi2 = dual_vertices.size(); + dual_vertices.push_back({{fmod(L.x + e->pos[1].x, L.x), fmod(L.y + e->pos[1].y, L.y)},{}}); + known_vertices[t2] = vi2; } else { vi2 = it2->second; } - edges.push_back({vi1, vi2}); + dual_edges.push_back({{vi1, vi2}, + {fmod(L.x + (e->pos[0].x + e->pos[1].x) / 2, L.x), + fmod(L.y + (e->pos[0].y + e->pos[1].y) / 2, L.y)} + }); } + ep = e; e = e->next; } } } - for (edge &e : edges) { - std::cout << e[0] << " " << e[1] << " "; - } - std::cout << "\n"; - for (vertex &v : vertices) { - std::cout << v.r.x << " " << v.r.y << " "; + if (edges.size() != vertices.size() + dual_vertices.size()) { + throw eulerex; } - std::cout << "\n"; - for (edge &e : dual_edges) { - std::cout << e[0] << " " << e[1] << " "; - } - std::cout << "\n"; - for (vertex &v : dual_vertices) { - std::cout << v.r.x << " " << v.r.y << " "; + + for (vertex &v : dual_vertices) { + if (fabs(v.polygon[0].x - v.polygon[1].x) > L.x / 2) { + if (v.polygon[0].x < L.x / 2) { + v.polygon[0].x += L.x; + } else { + v.polygon[1].x += L.x; + } + } + + if (fabs(v.polygon[0].y - v.polygon[1].y) > L.y / 2) { + if (v.polygon[0].y < L.y / 2) { + v.polygon[0].y += L.y; + } else { + v.polygon[1].y += L.y; + } + } + + if (fabs(v.polygon[2].x - v.polygon[1].x) > L.x / 2) { + if (v.polygon[2].x < L.x / 2) { + v.polygon[2].x += L.x; + } else { + v.polygon[1].x += L.x; + } + } + + if (fabs(v.polygon[2].y - v.polygon[1].y) > L.y / 2) { + if (v.polygon[2].y < L.y / 2) { + v.polygon[2].y += L.y; + } else { + v.polygon[1].y += L.y; + } + } + + if (fabs(v.polygon[2].x - v.polygon[0].x) > L.x / 2) { + if (v.polygon[2].x < L.x / 2) { + v.polygon[2].x += L.x; + } else { + v.polygon[0].x += L.x; + } + } + + if (fabs(v.polygon[2].y - v.polygon[0].y) > L.y / 2) { + if (v.polygon[2].y < L.y / 2) { + v.polygon[2].y += L.y; + } else { + v.polygon[0].y += L.y; + } + } } - std::cout << "\n"; jcv_diagram_free(&diagram); } diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 9cb1007..4d4ed2d 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -4,14 +4,14 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) { b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c); for (unsigned int i = 0; i < G.edges.size(); i++) { - double v0y = G.vertices[G.edges[i][0]].r.y; - double v1y = G.vertices[G.edges[i][1]].r.y; + double v0y = G.vertices[G.edges[i].v[0]].r.y; + double v1y = G.vertices[G.edges[i].v[1]].r.y; if (fabs(v0y - v1y) > G.L.y / 2) { bool ind = v0y < v1y ? 0 : 1; - ((double *)b->x)[G.edges[i][ind]] += 1.0; - ((double *)b->x)[G.edges[i][!ind]] -= 1.0; + ((double *)b->x)[G.edges[i].v[ind]] += 1.0; + ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0; } } @@ -30,8 +30,8 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges. unsigned int terms = G.vertices.size(); for (unsigned int i = 0; i < G.edges.size(); i++) { - unsigned int v0 = G.edges[i][0]; - unsigned int v1 = G.edges[i][1]; + unsigned int v0 = G.edges[i].v[0]; + unsigned int v1 = G.edges[i].v[1]; unsigned int s0 = v0 < v1 ? v0 : v1; unsigned int s1 = v0 < v1 ? v1 : v0; @@ -60,11 +60,11 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges. for (unsigned int i = 0; i < G.edges.size(); i++) { ((CHOL_INT *)t->i)[2 * i] = i; - ((CHOL_INT *)t->j)[2 * i] = G.edges[i][0]; + ((CHOL_INT *)t->j)[2 * i] = G.edges[i].v[0]; ((double *)t->x)[2 * i] = 1.0; ((CHOL_INT *)t->i)[2 * i + 1] = i; - ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i][1]; + ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1]; ((double *)t->x)[2 * i + 1] = -1.0; } @@ -99,8 +99,8 @@ void network::set_thresholds(double beta, std::mt19937& rng) { void network::break_edge(unsigned int e, bool unbreak) { fuses[e] = !unbreak; - unsigned int v0 = G.edges[e][0]; - unsigned int v1 = G.edges[e][1]; + unsigned int v0 = G.edges[e].v[0]; + unsigned int v1 = G.edges[e].v[1]; unsigned int n = factor->n; @@ -141,8 +141,8 @@ void network::break_edge(unsigned int e, bool unbreak) { if (fabs(v0y - v1y) > G.L.y / 2) { bool ind = v0y < v1y ? unbreak : !unbreak; - ((double *)b->x)[G.edges[e][ind]] -= 1.0; - ((double *)b->x)[G.edges[e][!ind]] += 1.0; + ((double *)b->x)[G.edges[e].v[ind]] -= 1.0; + ((double *)b->x)[G.edges[e].v[!ind]] += 1.0; } } @@ -169,8 +169,8 @@ current_info network::get_current_info() { } else { currents[i] = ((double *)y->x)[i]; - double v0y = G.vertices[G.edges[i][0]].r.y; - double v1y = G.vertices[G.edges[i][1]].r.y; + double v0y = G.vertices[G.edges[i].v[0]].r.y; + double v1y = G.vertices[G.edges[i].v[1]].r.y; if (fabs(v0y - v1y) > G.L.y / 2) { if (v0y > v1y) { diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d93783c..283eb10 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,7 +1,18 @@ +add_library(analysis_tools analysis_tools.cpp) +target_link_libraries(analysis_tools frac) + add_library(measurements measurements.cpp) -target_link_libraries(measurements frac) +target_link_libraries(measurements frac analysis_tools) + +add_library(animate animate.cpp) +target_link_libraries(animate frac analysis_tools GL GLU glut) add_executable(fracture fracture.cpp) -target_link_libraries(fracture frac measurements fftw3 cholmod profiler) +target_link_libraries(fracture frac measurements fftw3 cholmod profiler GL GLU glut) + +add_executable(animate_fracture animate_fracture.cpp) +target_link_libraries(animate_fracture frac animate fftw3 cholmod profiler GL GLU glut) +add_executable(test test.cpp) +target_link_libraries(test frac measurements fftw3 cholmod profiler) diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp new file mode 100644 index 0000000..1e799bb --- /dev/null +++ b/src/analysis_tools.cpp @@ -0,0 +1,125 @@ + +#include "analysis_tools.hpp" + +template +bool is_shorter(const std::list &l1, const std::list &l2) { + return l1.size() < l2.size(); +} + +bool trivial(boost::detail::edge_desc_impl) { + return true; +} + +std::list find_minimal_crack(const Graph& G, const network& n) { + Graph Gtmp(n.G.vertices.size()); + std::list removed_edges; + + class add_tree_edges : public boost::default_dfs_visitor { + public: + Graph& G; + std::list& E; + + add_tree_edges(Graph& G, std::list& E) : G(G), E(E) {} + + void tree_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { + boost::add_edge(boost::source(e, g), boost::target(e, g), g[e], G); + } + + void back_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { + if (!(boost::edge(boost::source(e, g), boost::target(e, g), G).second)) { + E.push_back(g[e].index); + } + } + }; + + add_tree_edges ate(Gtmp, removed_edges); + boost::depth_first_search(G, visitor(ate)); + + class find_cycle : public boost::default_dfs_visitor { + public: + std::list& E; + unsigned int end; + struct done{}; + + find_cycle(std::list& E, unsigned int end) : E(E), end(end) {} + + void discover_vertex(boost::graph_traits::vertex_descriptor v, const Graph& g) { + if (v == end) { + throw done{}; + } + } + + void examine_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { + E.push_back(g[e].index); + } + + void finish_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { + E.erase(std::find(E.begin(), E.end(), g[e].index)); + } + }; + + std::list> cycles; + + for (auto edge : removed_edges) { + std::list cycle = {edge}; + find_cycle vis(cycle, n.G.dual_edges[edge].v[1]); + std::vector new_color_map(boost::num_vertices(Gtmp)); + try { + boost::depth_first_visit(Gtmp, n.G.dual_edges[edge].v[0], vis, boost::make_iterator_property_map(new_color_map.begin(), boost::get(boost::vertex_index, Gtmp), new_color_map[0])); + } catch(find_cycle::done const&) { + cycles.push_back(cycle); + } + } + + if (cycles.size() > 1) { + std::list> bool_cycles; + for (auto cycle : cycles) { + std::valarray bool_cycle(n.G.edges.size()); + for (auto v : cycle) { + bool_cycle[v] = 1; + } + bool_cycles.push_back(bool_cycle); + } + + // generate all possible cycles by taking xor of the edge sets of the known cycles + for (auto it1 = bool_cycles.begin(); it1 != std::prev(bool_cycles.end()); it1++) { + for (auto it2 = std::next(it1); it2 != bool_cycles.end(); it2++) { + std::valarray new_bool_cycle = (*it1) ^ (*it2); + std::list new_cycle; + unsigned int pos = 0; + for (uint8_t included : new_bool_cycle) { + if (included) { + new_cycle.push_back(pos); + } + pos++; + } + cycles.push_back(new_cycle); + } + } + + // find the cycle representing the crack by counting boundary crossings + for (auto cycle : cycles) { + std::array crossing_count{0,0}; + + for (auto edge : cycle) { + double dx = fabs(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r.x - n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r.x); + if (dx > n.G.L.x / 2) { + crossing_count[0]++; + } + double dy = fabs(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r.y - n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r.y); + if (dy > n.G.L.y / 2) { + crossing_count[1]++; + } + } + + if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) { + return cycle; + } + } + } else { + return cycles.front(); + } + + exit(5); +} + diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp new file mode 100644 index 0000000..34ed967 --- /dev/null +++ b/src/analysis_tools.hpp @@ -0,0 +1,29 @@ + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +struct EdgeProperties { + unsigned int index; +}; + +typedef boost::adjacency_list Graph; + + +template +bool is_shorter(const std::list &, const std::list &); + +bool trivial(boost::detail::edge_desc_impl); + +std::list find_minimal_crack(const Graph &, const network &); + diff --git a/src/animate.cpp b/src/animate.cpp new file mode 100644 index 0000000..8c2d12f --- /dev/null +++ b/src/animate.cpp @@ -0,0 +1,205 @@ + +#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) { + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); + glutInitWindowSize(window_size, (unsigned int )(window_size * ((double)Ly / (double)Lx))); + 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); +} + +void animate::pre_fracture(const network &n) { + lv = std::numeric_limits::lowest(); + avalanches = {{}}; + boost::remove_edge_if(trivial, G); + + glClearColor(1.0f, 1.0f, 1.0f, 1.0f ); + glClear(GL_COLOR_BUFFER_BIT); + glBegin(GL_LINES); + glColor3f(0.0f, 0.0f, 0.0f); + for (unsigned int i = 0; i < n.G.edges.size(); i++) { + 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; + } else { + r2.x += Lx; + } + } + if (fabs(r1.y - r2.y) > Ly / 2) { + if (r1.y < Ly / 2) { + r1.y += Ly; + } else { + r2.y += Ly; + } + } + + glVertex2d(r1.x, r1.y); + glVertex2d(r2.x, r2.y); + } + glEnd(); + 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]; + if (c > lv && avalanches.back().size() > 0) { + lv = c; + avalanches.push_back({i}); + } else { + avalanches.back().push_back(i); + } + + boost::add_edge(net.G.dual_edges[i].v[0], net.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; + + if (fabs(r1.x - r2.x) > Lx / 2) { + if (r1.x < Lx / 2) { + r2.x -= Lx; + } else { + r2.x += Lx; + } + } + if (fabs(r1.y - r2.y) > Ly / 2) { + if (r1.y < Ly / 2) { + r2.y -= Ly; + } else { + r2.y += Ly; + } + } + + glVertex2d(r1.x, r1.y); + glVertex2d(r2.x, r2.y); + glEnd(); + glFlush(); +} + +void animate::post_fracture(network &n) { + std::vector component(boost::num_vertices(G)); + unsigned int num = boost::connected_components(G, &component[0]); + + std::list crack = find_minimal_crack(G, n); + unsigned int crack_component = component[n.G.dual_edges[crack.front()].v[0]]; + + std::vector> components(num); + + for (unsigned int i = 0; i < n.G.dual_vertices.size(); i++) { + components[component[i]].push_back(i); + } + + std::default_random_engine gen; + std::uniform_real_distribution dis(0.0,1.0); + + char key; + while ((key = getchar()) != 'n') { + glClearColor(1.0f, 1.0f, 1.0f, 1.0f ); + glClear(GL_COLOR_BUFFER_BIT); + glBegin(GL_LINES); + glColor3f(0.0f, 0.0f, 0.0f); + for (unsigned int i = 0; i < n.G.edges.size(); i++) { + if (!n.fuses[i]) { + 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; + } else { + r2.x += Lx; + } + } + if (fabs(r1.y - r2.y) > Ly / 2) { + if (r1.y < Ly / 2) { + r1.y += Ly; + } else { + r2.y += Ly; + } + } + + glVertex2d(r1.x, r1.y); + glVertex2d(r2.x, r2.y); + } + } + glEnd(); + + switch (key) { + case 's' : + for (auto edge : crack) { + glBegin(GL_TRIANGLES); + glColor3f(1.0f, 0.0f, 0.0f); + glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[0].x, n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[0].y); + glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[1].x, n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[1].y); + glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[2].x, n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon[2].y); + + glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[0].x, n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[0].y); + glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[1].x, n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[1].y); + glVertex2d(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[2].x, n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon[2].y); + glEnd(); + } + glFlush(); + break; + + case 'c' : + glBegin(GL_TRIANGLES); // Each set of 3 vertices form a triangle + for (unsigned int i = 0; i < num; i++) { + if (i == crack_component) { + glColor3d(1.0, 0.0, 0.0); + } else { + glColor3d(dis(gen), dis(gen), dis(gen)); + } + + for (auto it = components[i].begin(); it != components[i].end(); it++) { + glVertex2d(n.G.dual_vertices[*it].polygon[0].x, n.G.dual_vertices[*it].polygon[0].y); + glVertex2d(n.G.dual_vertices[*it].polygon[1].x, n.G.dual_vertices[*it].polygon[1].y); + glVertex2d(n.G.dual_vertices[*it].polygon[2].x, n.G.dual_vertices[*it].polygon[2].y); + } + } + glEnd(); + glFlush(); + break; + + case 'C' : + glBegin(GL_TRIANGLES); // Each set of 3 vertices form a triangle + for (unsigned int i = 0; i < num; i++) { + if (components[i].size() > 1) { + if (i == crack_component) { + glColor3d(1.0, 0.0, 0.0); + } else { + glColor3d(dis(gen), dis(gen), dis(gen)); + } + + for (auto it = components[i].begin(); it != components[i].end(); it++) { + glVertex2d(n.G.dual_vertices[*it].polygon[0].x, n.G.dual_vertices[*it].polygon[0].y); + glVertex2d(n.G.dual_vertices[*it].polygon[1].x, n.G.dual_vertices[*it].polygon[1].y); + glVertex2d(n.G.dual_vertices[*it].polygon[2].x, n.G.dual_vertices[*it].polygon[2].y); + } + } + } + glEnd(); + glFlush(); + break; + + case 'm' : + glBegin(GL_TRIANGLES); + glColor3d(1.0, 0.0, 0.0); + for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) { + glVertex2d(n.G.dual_vertices[*it].polygon[0].x, n.G.dual_vertices[*it].polygon[0].y); + glVertex2d(n.G.dual_vertices[*it].polygon[1].x, n.G.dual_vertices[*it].polygon[1].y); + glVertex2d(n.G.dual_vertices[*it].polygon[2].x, n.G.dual_vertices[*it].polygon[2].y); + } + glEnd(); + glFlush(); + } + } +} + diff --git a/src/animate.hpp b/src/animate.hpp new file mode 100644 index 0000000..793771d --- /dev/null +++ b/src/animate.hpp @@ -0,0 +1,24 @@ + + +#include +#include + +#include "analysis_tools.hpp" + +#include + +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[]); + + void pre_fracture(const network &) override; + void bond_broken(const network& net, const current_info& cur, unsigned int i) override; + void post_fracture(network &n) override; +}; diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp new file mode 100644 index 0000000..10b63d0 --- /dev/null +++ b/src/animate_fracture.cpp @@ -0,0 +1,79 @@ + +#include +#include + +#include + +#include "randutils/randutils.hpp" + +#include +#include +#include +#include "animate.hpp" + +#include +#include +#include + +std::atomic quit(false); // signal flag + +void got_signal(int) { + quit.store(true); +} + +int main(int argc, char* argv[]) { + struct sigaction sa; + memset( &sa, 0, sizeof(sa) ); + sa.sa_handler = got_signal; + sigfillset(&sa.sa_mask); + sigaction(SIGINT, &sa, NULL); + + int opt; + + unsigned int N = 1; + unsigned int Lx = 16; + unsigned int Ly = 16; + double beta = 0.5; + + while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) { + switch (opt) { + case 'N': + N = (unsigned int)atof(optarg); + break; + case 'X': + Lx = atoi(optarg); + break; + case 'Y': + Ly = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + default: + exit(1); + } + } + + cholmod_common c; + CHOL_F(start)(&c); + + animate meas(Lx, Ly, 512, argc, argv); + + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + for (unsigned int trial = 0; trial < N; trial++) { + graph G(Lx, Ly, rng, 0.4); + network network(G, &c); + network.set_thresholds(beta, rng); + network.fracture(meas); + + if (quit.load()) + break; + } + + CHOL_F(finish)(&c); + + return 0; +} + diff --git a/src/fracture.cpp b/src/fracture.cpp index 317d91b..4b6fb5b 100644 --- a/src/fracture.cpp +++ b/src/fracture.cpp @@ -57,17 +57,22 @@ int main(int argc, char* argv[]) { cholmod_common c; CHOL_F(start)(&c); - graph G(Lx, Ly); - network base_network(G, &c); - ma meas(Lx, Ly, beta); + ma meas(Lx, Ly, Lx, Ly, beta); randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; for (unsigned int trial = 0; trial < N; trial++) { - network tmp_network(base_network); - tmp_network.set_thresholds(beta, rng); - tmp_network.fracture(meas); + while (true) { + try { + graph G(Lx, Ly, rng, 0.4); + network network(G, &c); + network.set_thresholds(beta, rng); + network.fracture(meas); + break; + } catch (std::exception &e) { + } + } if (quit.load()) break; diff --git a/src/measurements.cpp b/src/measurements.cpp index c8cc73c..43483be 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -2,128 +2,6 @@ #include "measurements.hpp" #include -template -bool is_shorter(std::list l1, std::list l2) { - return l1.size() < l2.size(); -} - -bool trivial(boost::detail::edge_desc_impl) { - return true; -} - -std::list find_minimal_crack(const Graph& G, const network& n) { - Graph Gtmp(n.G.vertices.size()); - std::list removed_edges; - - class add_tree_edges : public boost::default_dfs_visitor { - public: - Graph& G; - std::list& E; - - add_tree_edges(Graph& G, std::list& E) : G(G), E(E) {} - - void tree_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { - boost::add_edge(boost::source(e, g), boost::target(e, g), g[e], G); - } - - void back_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { - if (!(boost::edge(boost::source(e, g), boost::target(e, g), G).second)) { - E.push_back(g[e].index); - } - } - }; - - add_tree_edges ate(Gtmp, removed_edges); - boost::depth_first_search(G, visitor(ate)); - - class find_cycle : public boost::default_dfs_visitor { - public: - std::list& E; - unsigned int end; - struct done{}; - - find_cycle(std::list& E, unsigned int end) : E(E), end(end) {} - - void discover_vertex(boost::graph_traits::vertex_descriptor v, const Graph& g) { - if (v == end) { - throw done{}; - } - } - - void examine_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { - E.push_back(g[e].index); - } - - void finish_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { - E.erase(std::find(E.begin(), E.end(), g[e].index)); - } - }; - - std::list> cycles; - - for (auto edge : removed_edges) { - std::list cycle = {edge}; - find_cycle vis(cycle, n.G.dual_edges[edge][1]); - std::vector new_color_map(boost::num_vertices(Gtmp)); - try { - boost::depth_first_visit(Gtmp, n.G.dual_edges[edge][0], vis, boost::make_iterator_property_map(new_color_map.begin(), boost::get(boost::vertex_index, Gtmp), new_color_map[0])); - } catch(find_cycle::done const&) { - cycles.push_back(cycle); - } - } - - if (cycles.size() > 1) { - std::list> bool_cycles; - for (auto cycle : cycles) { - std::valarray bool_cycle(n.G.edges.size()); - for (auto v : cycle) { - bool_cycle[v] = 1; - } - bool_cycles.push_back(bool_cycle); - } - - // generate all possible cycles by taking xor of the edge sets of the known cycles - for (auto it1 = bool_cycles.begin(); it1 != std::prev(bool_cycles.end()); it1++) { - for (auto it2 = std::next(it1); it2 != bool_cycles.end(); it2++) { - std::valarray new_bool_cycle = (*it1) ^ (*it2); - std::list new_cycle; - unsigned int pos = 0; - for (uint8_t included : new_bool_cycle) { - if (included) { - new_cycle.push_back(pos); - } - pos++; - } - cycles.push_back(new_cycle); - } - } - - // find the cycle representing the crack by counting boundary crossings - for (auto cycle : cycles) { - std::array crossing_count{0,0}; - - for (auto edge : cycle) { - double dx = fabs(n.G.dual_vertices[n.G.dual_edges[edge][0]].r.x - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.x); - if (dx > n.G.L.x / 2) { - crossing_count[0]++; - } - double dy = fabs(n.G.dual_vertices[n.G.dual_edges[edge][0]].r.y - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.y); - if (dy > n.G.L.y / 2) { - crossing_count[1]++; - } - } - - if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) { - return cycle; - } - } - } else { - return cycles.front(); - } - - exit(5); -} - 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); @@ -153,8 +31,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 Lx, unsigned int Ly, double beta) { - std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + 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, 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); uint64_t N_old = 0; @@ -226,18 +104,22 @@ void autocorrelation(unsigned int Lx, unsigned int Ly, std::vector& out_data, } } -ma::ma(unsigned int Lx, unsigned int Ly, double beta) : Lx(Lx), Ly(Ly), beta(beta), G(Lx * Ly / 2), +unsigned int edge_r_to_ind(graph::coordinate r, unsigned int Lx, unsigned int 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), - 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) + sb(log2(Mx < My ? Mx : My) + 1, 0), + Ccc((Mx / 2 + 1) * (My / 2 + 1), 0), + Css((Mx / 2 + 1) * (My / 2 + 1), 0), + Cmm((Mx / 2 + 1) * (My / 2 + 1), 0), + Caa((Mx / 2 + 1) * (My / 2 + 1), 0), + Cdd((Mx / 2 + 1) * (My / 2 + 1), 0), + Cll((Mx / 2 + 1) * (My / 2 + 1), 0), + Cee((Mx / 2 + 1) * (My / 2 + 1), 0) { N = 0; Nc = 0; @@ -246,17 +128,13 @@ ma::ma(unsigned int Lx, unsigned int Ly, double beta) : Lx(Lx), Ly(Ly), beta(bet // FFTW setup for correlation functions fftw_set_timelimit(1); - 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(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); - - std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_broken.dat"; + 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)); - bondfile.open(filename, std::ifstream::app); + 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() { @@ -274,18 +152,16 @@ ma::~ma() { 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); - - bondfile.close(); + 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("Cdd", Cdd, 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("Cee", Cee, N, Lx, Ly, beta, Mx, My); } -void ma::pre_fracture(const network &) { +void ma::pre_fracture(const network&) { lv = std::numeric_limits::lowest(); avalanches = {{}}; boost::remove_edge_if(trivial, G); @@ -297,13 +173,13 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned int i sa[avalanches.back().size() - 1]++; Na++; - std::fill_n(fftw_forward_in, net.G.edges.size(), 0.0); + std::fill_n(fftw_forward_in, Mx * My, 0.0); for (auto e : avalanches.back()) { - fftw_forward_in[e] = 1.0; + fftw_forward_in[edge_r_to_ind(net.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0; } - autocorrelation(Lx, Ly, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(Mx, My, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); lv = c; avalanches.push_back({i}); @@ -311,121 +187,116 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned int i avalanches.back().push_back(i); } - boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], {i}, G); - - bondfile << i << " " << c << std::scientific << " "; + boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G); } void ma::post_fracture(network &n) { - bondfile << "\n"; std::vector component(boost::num_vertices(G)); unsigned int num = boost::connected_components(G, &component[0]); std::list crack = find_minimal_crack(G, n); - unsigned int crack_component = component[n.G.dual_edges[crack.front()][0]]; + // crack surface correlations + std::fill_n(fftw_forward_in, Mx * My, 0.0); + + for (auto edge : crack) { + fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r, Lx, Ly, Mx, My)] = 0.5; + fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My)] = 0.5; + } + + autocorrelation(Mx, My, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + + unsigned int crack_component = component[n.G.dual_edges[crack.front()].v[0]]; + + std::vector> components(num); + + for (unsigned int i = 0; i < n.G.dual_vertices.size(); i++) { + components[component[i]].push_back(i); + } // 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) { - unsigned int size = 0; - - for (unsigned int j = 0; j < n.G.edges.size(); j++) { - if (component[n.G.dual_edges[j][0]] == i && n.fuses[j]) { - size++; - fftw_forward_in[j] = 1.0; - } else{ - fftw_forward_in[j] = 0.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; } - if (size > 0) { - sc[size - 1]++; - autocorrelation(Lx, Ly, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); - Nc++; - } - } - } + 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++; - // bin counting - for (unsigned int be = 0; be < log2(Ly); be++) { - unsigned int bin = pow(2, be); - - 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 = 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; - } - } - - if (in_bin) { - sb[be]++; + 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; } } } - sb[log2(Ly)]++; + unsigned int max_factor = log2(Mx < My ? Mx : My) + 1; + std::vector> bins(max_factor); + for (unsigned int i = 0; i < max_factor; i++) { + bins[i].resize(Mx * My / (unsigned int)pow(2, 2 * i)); + } - for (unsigned int i = 0; i < n.G.edges.size(); i++) { - if (n.fuses[i] && component[n.G.dual_edges[i][0]] == crack_component) { - fftw_forward_in[i] = 1.0; - } else { - fftw_forward_in[i] = 0.0; + for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) { + for (unsigned int i = 0; i < max_factor; i++) { + bins[i][edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx / pow(2, i), My / pow(2, i))] = 1; } } - autocorrelation(Lx, Ly, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + for (unsigned int i =0 ; i < max_factor; i++) { + sb[i] += bins[i].sum(); + } - // crack surface correlations - std::fill_n(fftw_forward_in, n.G.edges.size(), 0.0); + // spanning cluster + std::fill_n(fftw_forward_in, Mx * My, 0.0); - for (auto edge : crack) { - fftw_forward_in[edge] = 1.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; + } } - autocorrelation(Lx, Ly, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(Mx, My, Cmm, 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]; + return component[n.G.dual_edges[i].v[0]] == crack_component; }; for (auto avalanche : avalanches) { if (avalanche.end() != std::find_if(avalanche.begin(), avalanche.end(), inCrack)) { for (auto edge : avalanche) { - fftw_forward_in[edge] = 1.0; + fftw_forward_in[edge_r_to_ind(n.G.edges[edge].r, Lx, Ly, Mx, My)] += 1.0; } } } - autocorrelation(Lx, Ly, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(Mx, My, 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); + std::fill_n(fftw_forward_in, Mx * My, 0.0); // rewind the last avalanche for (auto e : avalanches.back()) { - boost::remove_edge(n.G.dual_edges[e][0], n.G.dual_edges[e][1], G); + boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G); n.break_edge(e, true); - fftw_forward_in[e] = 1; + fftw_forward_in[edge_r_to_ind(n.G.edges[e].r, Lx, Ly, Mx, My)] += 1.0; } - autocorrelation(Lx, Ly, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(Mx, My, 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; + std::fill_n(fftw_forward_in, Mx * My, 0.0); + for (unsigned int i = 0; i < n.G.edges.size(); i++) { if (n.fuses[i]) { total_broken++; - fftw_forward_in[i] = 1.0; - } else { - fftw_forward_in[i] = 0.0; + fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0; } } - autocorrelation(Lx, Ly, Cdd, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + autocorrelation(Mx, My, 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 2bc5bae..685806b 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -1,34 +1,18 @@ -#include -#include -#include #include -#include #include #include #include #include #include #include -#include #include -#include -#include -#include -#include -#include - #include -#include #include -struct EdgeProperties { - unsigned int index; -}; - -typedef boost::adjacency_list Graph; +#include "analysis_tools.hpp" class ma : public hooks { // need: @@ -44,6 +28,8 @@ class ma : public hooks { unsigned int N; unsigned int Lx; unsigned int Ly; + unsigned int Mx; + unsigned int My; double beta; Graph G; @@ -62,15 +48,13 @@ class ma : public hooks { uint64_t Nc; uint64_t Na; - std::ofstream bondfile; - public: long double lv; std::list> avalanches; - ma(unsigned int Lx, unsigned int Ly, double beta); + ma(unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My, double beta); ~ma(); void pre_fracture(const network &) override; diff --git a/src/test.cpp b/src/test.cpp new file mode 100644 index 0000000..2ebf3ac --- /dev/null +++ b/src/test.cpp @@ -0,0 +1,83 @@ + +#include +#include + +#include + +#include "randutils/randutils.hpp" + +#include +#include +#include +#include "measurements.hpp" + +#include +#include +#include + +std::atomic quit(false); // signal flag + +void got_signal(int) { + quit.store(true); +} + +int main(int argc, char* argv[]) { + struct sigaction sa; + memset( &sa, 0, sizeof(sa) ); + sa.sa_handler = got_signal; + sigfillset(&sa.sa_mask); + sigaction(SIGINT, &sa, NULL); + + int opt; + + unsigned int N = 1; + unsigned int Lx = 16; + unsigned int Ly = 16; + double beta = 0.5; + + while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) { + switch (opt) { + case 'N': + N = (unsigned int)atof(optarg); + break; + case 'X': + Lx = atoi(optarg); + break; + case 'Y': + Ly = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + default: + exit(1); + } + } + + cholmod_common c; + CHOL_F(start)(&c); + + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + graph G(Lx, Ly, rng); + return 0; + /* + network base_network(G, &c); + ma meas(Lx, Ly, beta); + + for (unsigned int trial = 0; trial < N; trial++) { + network tmp_network(base_network); + tmp_network.set_thresholds(beta, rng); + tmp_network.fracture(meas); + + if (quit.load()) + break; + } + + CHOL_F(finish)(&c); + + return 0; + */ +} + -- cgit v1.2.3-70-g09d2