summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/CMakeLists.txt2
-rw-r--r--src/measurements.cpp344
-rw-r--r--src/measurements.hpp6
3 files changed, 184 insertions, 168 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 57a817b..d93783c 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,5 +3,5 @@ add_library(measurements measurements.cpp)
target_link_libraries(measurements frac)
add_executable(fracture fracture.cpp)
-target_link_libraries(fracture frac measurements fftw3 cholmod)
+target_link_libraries(fracture frac measurements fftw3 cholmod profiler)
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 578c70a..0bd72c4 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -11,6 +11,119 @@ 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:
+ unsigned int end;
+ std::list<unsigned int>& E;
+ 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.vertices[n.G.dual_edges[edge][0]].r.x - n.G.vertices[n.G.dual_edges[edge][1]].r.x);
+ if (dx > 0.5) {
+ crossing_count[0]++;
+ }
+ double dy = fabs(n.G.vertices[n.G.dual_edges[edge][0]].r.y - n.G.vertices[n.G.dual_edges[edge][1]].r.y);
+ if (dy > 0.5) {
+ 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 L, double beta) {
std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat";
std::ifstream file(filename);
@@ -97,12 +210,8 @@ void correlation(unsigned int L, std::vector<T>& data, const std::vector<fftw_co
}
}
-template <class T, class U>
-void autocorrelation(unsigned int L, std::vector<T>& out_data, const std::vector<U>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
- for (unsigned int i = 0; i < pow(L, 2); i++) {
- fftw_forward_in[i] = (double)data[i];
- }
-
+template <class T>
+void autocorrelation(unsigned int L, std::vector<T>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
fftw_execute(forward_plan);
for (unsigned int i = 0; i < L * (L / 2 + 1); i++) {
@@ -124,8 +233,11 @@ ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2))
sb(log2(L) + 1, 0),
Ccc(pow(L / 2 + 1, 2), 0),
Css(pow(L / 2 + 1, 2), 0),
+ Cmm(pow(L / 2 + 1, 2), 0),
+ Caa(pow(L / 2 + 1, 2), 0),
Cdd(pow(L / 2 + 1, 2), 0),
- Caa(pow(L / 2 + 1, 2), 0)
+ Cll(pow(L / 2 + 1, 2), 0),
+ Cee(pow(L / 2 + 1, 2), 0)
{
Nc = 0;
Na = 0;
@@ -159,34 +271,36 @@ ma::~ma() {
update_field_file("Ccc", Ccc, Nc, L, beta);
update_field_file("Css", Css, N, L, beta);
+ update_field_file("Cmm", Cmm, N, L, beta);
update_field_file("Cdd", Cdd, N, L, beta);
update_field_file("Caa", Caa, Na, L, beta);
-
+ update_field_file("Cll", Cll, N, L, beta);
+ update_field_file("Cee", Cee, N, L, beta);
}
void ma::pre_fracture(const network &) {
lv = 0;
- avalanche.clear();
+ avalanches = {{}};
boost::remove_edge_if(trivial, G);
}
void ma::bond_broken(const network& net, const std::pair<double, std::vector<double>>& cur, unsigned int i) {
if (cur.first / fabs(cur.second[i]) * net.thresholds[i] > lv) {
- sa[avalanche.size()]++;
+ sa[avalanches.back().size()]++;
Na++;
- std::vector<bool>avalanche_damage(pow(L, 2), 0);
+ memset(fftw_forward_in, 0.0, net.G.edges.size());
- for (auto it = avalanche.begin(); it != avalanche.end(); it++) {
- avalanche_damage[*it] = 1;
+ for (auto e : avalanches.back()) {
+ fftw_forward_in[e] = 1.0;
}
- autocorrelation(L, Caa, avalanche_damage, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation(L, Caa, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
lv = cur.first / fabs(cur.second[i]) * net.thresholds[i];
- avalanche.clear();
+ avalanches.push_back({i});
} else {
- avalanche.push_back(i);
+ avalanches.back().push_back(i);
}
boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], {i}, G);
@@ -196,135 +310,9 @@ void ma::post_fracture(network &n) {
std::vector<unsigned int> component(boost::num_vertices(G));
unsigned int num = boost::connected_components(G, &component[0]);
- std::vector<unsigned int> comp_sizes(num, 0);
-
- for (unsigned int c : component) {
- comp_sizes[c]++;
- }
+ std::list<unsigned int> crack = find_minimal_crack(G, n);
- unsigned int max_i = 0;
- unsigned int max_size = 0;
-
- for (unsigned int i = 0; i < num; i++) {
- if (comp_sizes[i] > max_size) {
- max_i = i;
- max_size = comp_sizes[i];
- }
- }
-
- unsigned int vertex_in_largest = 0;
-
- while (component[vertex_in_largest] != max_i) {
- vertex_in_largest++;
- }
-
- Graph Gtmp(2 * pow(L / 2, 2));
- 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:
- unsigned int end;
- std::list<unsigned int>& E;
- 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);
- }
- }
-
- std::list<unsigned int> crack;
-
- 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);
- }
-
- 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);
- }
- }
-
- for (auto cycle : cycles) {
- std::array<unsigned int, 2> crossing_count{0,0};
-
- for (auto edge : cycle) {
- double dx = fabs(n.G.vertices[n.G.edges[edge][0]].r.x - n.G.vertices[n.G.edges[edge][1]].r.x);
- if (dx > 0.5) {
- crossing_count[0]++;
- }
- double dy = fabs(n.G.vertices[n.G.edges[edge][0]].r.y - n.G.vertices[n.G.edges[edge][1]].r.y);
- if (dy > 0.5) {
- crossing_count[1]++;
- }
- }
-
- if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
- crack = cycle;
- break;
- }
- }
- } else {
- crack = cycles.front();
- }
+ unsigned int crack_component = component[n.G.dual_edges[crack.front()][0]];
for (unsigned int be = 0; be < log2(L); be++) {
unsigned int bin = pow(2, be);
@@ -333,7 +321,7 @@ void ma::post_fracture(network &n) {
bool in_bin = false;
for (unsigned int j = 0; j < pow(bin, 2); j++) {
unsigned int edge = L * (bin * ((i * bin) / L) + j / bin) + (i * bin) % L + j % bin;
- if (!n.fuses[edge] && max_i == component[n.G.dual_edges[edge][0]]) {
+ if (!n.fuses[edge] && crack_component == component[n.G.dual_edges[edge][0]]) {
in_bin = true;
break;
}
@@ -347,62 +335,86 @@ void ma::post_fracture(network &n) {
sb[log2(L)]++;
+ 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;
+ }
+ }
+
+ autocorrelation(L, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+
// crack surface correlations
- std::vector<double> crack_damage(pow(L, 2), 0.0);
+ memset(fftw_forward_in, 0.0, n.G.edges.size());
for (auto edge : crack) {
- crack_damage[edge] = 1.0;
+ fftw_forward_in[edge] = 1.0;
+ }
+
+ autocorrelation(L, Css, 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];
+ };
+
+ 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;
+ }
+ }
}
- std::vector<fftw_complex> t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out);
- correlation(L, Css, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation(L, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
- std::vector<bool> last_crack(pow(L, 2), 0);
+ memset(fftw_forward_in, 0.0, n.G.edges.size());
// rewind the last avalanche
- for (auto e = avalanche.rbegin(); e != avalanche.rend(); e++) {
- boost::remove_edge(n.G.dual_edges[*e][0], n.G.dual_edges[*e][1], G);
- n.add_edge(*e);
- last_crack[*e] = 1;
+ for (auto e : avalanches.back()) {
+ boost::remove_edge(n.G.dual_edges[e][0], n.G.dual_edges[e][1], G);
+ n.add_edge(e);
+ fftw_forward_in[e] = 1;
}
+ autocorrelation(L, Cll, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
// cluster size distribution and cluster-cluster correlation
num = boost::connected_components(G, &component[0]);
for (unsigned int i = 0; i < num; i++) {
- double size = 0;
- std::fill(crack_damage.begin(), crack_damage.end(), 0.0);
+ unsigned int size = 0;
- for (unsigned int j = 0; j < pow(L, 2); j++) {
- if (component[n.G.edges[j][0]] == i && n.fuses[j]) {
+ for (unsigned int j = 0; j < n.G.edges.size(); j++) {
+ if (component[n.G.dual_edges[j][0]] == i && n.fuses[j]) {
size++;
- crack_damage[j] = 1.0;
+ fftw_forward_in[j] = 1.0;
+ } else{
+ fftw_forward_in[j] = 0.0;
}
}
if (size > 0) {
sc[size - 1]++;
- t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out);
- correlation(L, Ccc, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out);
+ autocorrelation(L, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
Nc++;
}
}
- // damage-damage correlations
- std::vector<fftw_complex> t_damage = data_transform(L, n.fuses, forward_plan, fftw_forward_in, fftw_forward_out);
- correlation(L, Cdd, t_damage, t_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out);
-
-
// damage size distribution
unsigned int total_broken = 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;
}
}
+ autocorrelation(L, 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 4c05d93..621c425 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -2,6 +2,7 @@
#include <vector>
#include <algorithm>
#include <cmath>
+#include <cstring>
#include <list>
#include <fstream>
#include <string>
@@ -49,8 +50,11 @@ class ma : public hooks {
std::vector<uint64_t> sb; // bin size distribution
std::vector<uint64_t> Ccc; // cluster-cluster correlations
std::vector<uint64_t> Css; // surface-surface correlations
+ std::vector<uint64_t> Cmm; // surface-surface correlations
std::vector<uint64_t> Caa; // avalanche-avalanche correlations
std::vector<uint64_t> Cdd; // damage-damage distribution
+ std::vector<uint64_t> Cll; // damage-damage distribution
+ std::vector<uint64_t> Cee; // damage-damage distribution
uint64_t Nc;
uint64_t Na;
@@ -61,7 +65,7 @@ class ma : public hooks {
double beta;
- std::list<unsigned int> avalanche;
+ std::list<std::list<unsigned int>> avalanches;
ma(unsigned int N, unsigned int L, double beta);
~ma();