summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/CMakeLists.txt15
-rw-r--r--src/analysis_tools.cpp125
-rw-r--r--src/analysis_tools.hpp29
-rw-r--r--src/animate.cpp205
-rw-r--r--src/animate.hpp24
-rw-r--r--src/animate_fracture.cpp79
-rw-r--r--src/fracture.cpp17
-rw-r--r--src/measurements.cpp303
-rw-r--r--src/measurements.hpp24
-rw-r--r--src/test.cpp83
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;
+ */
+}
+