diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/CMakeLists.txt | 15 | ||||
| -rw-r--r-- | src/analysis_tools.cpp | 125 | ||||
| -rw-r--r-- | src/analysis_tools.hpp | 29 | ||||
| -rw-r--r-- | src/animate.cpp | 205 | ||||
| -rw-r--r-- | src/animate.hpp | 24 | ||||
| -rw-r--r-- | src/animate_fracture.cpp | 79 | ||||
| -rw-r--r-- | src/fracture.cpp | 17 | ||||
| -rw-r--r-- | src/measurements.cpp | 303 | ||||
| -rw-r--r-- | src/measurements.hpp | 24 | ||||
| -rw-r--r-- | src/test.cpp | 83 |
10 files changed, 660 insertions, 244 deletions
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 <class T> +bool is_shorter(const std::list<T> &l1, const std::list<T> &l2) { + return l1.size() < l2.size(); +} + +bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) { + return true; +} + +std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) { + Graph Gtmp(n.G.vertices.size()); + std::list<unsigned int> removed_edges; + + class add_tree_edges : public boost::default_dfs_visitor { + public: + Graph& G; + std::list<unsigned int>& E; + + add_tree_edges(Graph& G, std::list<unsigned int>& E) : G(G), E(E) {} + + void tree_edge(boost::graph_traits<Graph>::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<Graph>::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<unsigned int>& E; + unsigned int end; + struct done{}; + + find_cycle(std::list<unsigned int>& E, unsigned int end) : E(E), end(end) {} + + void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) { + if (v == end) { + throw done{}; + } + } + + void examine_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { + E.push_back(g[e].index); + } + + void finish_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { + E.erase(std::find(E.begin(), E.end(), g[e].index)); + } + }; + + std::list<std::list<unsigned int>> cycles; + + for (auto edge : removed_edges) { + std::list<unsigned int> cycle = {edge}; + find_cycle vis(cycle, n.G.dual_edges[edge].v[1]); + std::vector<boost::default_color_type> 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<std::valarray<uint8_t>> bool_cycles; + for (auto cycle : cycles) { + std::valarray<uint8_t> 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<uint8_t> new_bool_cycle = (*it1) ^ (*it2); + std::list<unsigned int> 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<unsigned int, 2> 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 <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> +#include <boost/graph/depth_first_search.hpp> +#include <boost/range/combine.hpp> +#include <boost/foreach.hpp> + +#include <vector> +#include <algorithm> +#include <cmath> +#include <list> +#include <valarray> + +#include <network.hpp> + +struct EdgeProperties { + unsigned int index; +}; + +typedef boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, boost::no_property, EdgeProperties> Graph; + + +template<class T> +bool is_shorter(const std::list<T> &, const std::list<T> &); + +bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>); + +std::list<unsigned int> 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<long double>::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<unsigned int> component(boost::num_vertices(G)); + unsigned int num = boost::connected_components(G, &component[0]); + + std::list<unsigned int> crack = find_minimal_crack(G, n); + unsigned int crack_component = component[n.G.dual_edges[crack.front()].v[0]]; + + std::vector<std::list<unsigned int>> 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<double> 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 <functional> +#include <hooks.hpp> + +#include "analysis_tools.hpp" + +#include <GL/glut.h> + +class animate : public hooks { + private: + unsigned int Lx; + unsigned int Ly; + Graph G; + public: + long double lv; + std::list<std::list<unsigned int>> 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 <random> +#include <iostream> + +#include <cholmod.h> + +#include "randutils/randutils.hpp" + +#include <graph.hpp> +#include <network.hpp> +#include <hooks.hpp> +#include "animate.hpp" + +#include <csignal> +#include <cstring> +#include <atomic> + +std::atomic<bool> 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 <iostream> -template <class T> -bool is_shorter(std::list<T> l1, std::list<T> l2) { - return l1.size() < l2.size(); -} - -bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) { - return true; -} - -std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) { - Graph Gtmp(n.G.vertices.size()); - std::list<unsigned int> removed_edges; - - class add_tree_edges : public boost::default_dfs_visitor { - public: - Graph& G; - std::list<unsigned int>& E; - - add_tree_edges(Graph& G, std::list<unsigned int>& E) : G(G), E(E) {} - - void tree_edge(boost::graph_traits<Graph>::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<Graph>::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<unsigned int>& E; - unsigned int end; - struct done{}; - - find_cycle(std::list<unsigned int>& E, unsigned int end) : E(E), end(end) {} - - void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) { - if (v == end) { - throw done{}; - } - } - - void examine_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { - E.push_back(g[e].index); - } - - void finish_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) { - E.erase(std::find(E.begin(), E.end(), g[e].index)); - } - }; - - std::list<std::list<unsigned int>> cycles; - - for (auto edge : removed_edges) { - std::list<unsigned int> cycle = {edge}; - find_cycle vis(cycle, n.G.dual_edges[edge][1]); - std::vector<boost::default_color_type> 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<std::valarray<uint8_t>> bool_cycles; - for (auto cycle : cycles) { - std::valarray<uint8_t> 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<uint8_t> new_bool_cycle = (*it1) ^ (*it2); - std::list<unsigned int> 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<unsigned int, 2> 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<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); @@ -153,8 +31,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 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<T>& 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<T>& 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<long double>::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<unsigned int> component(boost::num_vertices(G)); unsigned int num = boost::connected_components(G, &component[0]); std::list<unsigned int> 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<std::list<unsigned int>> 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<std::valarray<unsigned int>> 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<bool(unsigned int)> 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 <vector> -#include <algorithm> -#include <cmath> #include <cstring> -#include <list> #include <fstream> #include <string> #include <cinttypes> #include <sstream> #include <functional> #include <iostream> -#include <valarray> #include <array> -#include <boost/graph/adjacency_list.hpp> -#include <boost/graph/connected_components.hpp> -#include <boost/graph/depth_first_search.hpp> -#include <boost/range/combine.hpp> -#include <boost/foreach.hpp> - #include <fftw3.h> -#include <network.hpp> #include <hooks.hpp> -struct EdgeProperties { - unsigned int index; -}; - -typedef boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, boost::no_property, EdgeProperties> 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<std::list<unsigned int>> 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 <random> +#include <iostream> + +#include <cholmod.h> + +#include "randutils/randutils.hpp" + +#include <graph.hpp> +#include <network.hpp> +#include <hooks.hpp> +#include "measurements.hpp" + +#include <csignal> +#include <cstring> +#include <atomic> + +std::atomic<bool> 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; + */ +} + |