summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/include/network.hpp8
-rw-r--r--lib/src/network.cpp78
-rw-r--r--src/analysis_tools.cpp157
-rw-r--r--src/analysis_tools.hpp8
-rw-r--r--src/animate.cpp139
-rw-r--r--src/animate_fracture.cpp4
-rw-r--r--src/animate_fracture_square.cpp4
7 files changed, 144 insertions, 254 deletions
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 29bf55a..5b474bb 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -35,9 +35,10 @@ private:
void update_backbone(const std::vector<double>& c);
void break_edge(unsigned, bool unbreak = false);
void get_cycle_edges_helper(std::set<unsigned>&, std::set<unsigned>&, unsigned, unsigned) const;
- std::set<unsigned> get_cycle_edges(unsigned) const;
bool find_cycle_helper(std::array<unsigned, 2>&, const std::set<unsigned>&, unsigned, unsigned,
unsigned) const;
+ bool get_cycle_helper(std::array<unsigned, 2>&, std::set<unsigned>&, const std::set<unsigned>&, unsigned, unsigned,
+ unsigned) const;
std::array<unsigned, 2> find_cycle(const std::set<unsigned>&, unsigned, unsigned) const;
void get_cluster_edges_helper(std::set<unsigned>&, unsigned) const;
std::set<unsigned> get_cluster_edges(unsigned) const;
@@ -56,7 +57,10 @@ public:
network(const graph&, cholmod_common*);
void set_thresholds(double, std::mt19937&);
- void fracture(hooks&, bool one_axis = true);
+ void fracture(hooks&);
+
+ std::set<unsigned> get_cycle_edges(unsigned) const;
+ std::pair<std::array<unsigned, 2>, std::set<unsigned>> get_cycle(const std::set<unsigned>&, unsigned, unsigned) const;
virtual current_info get_current_info() {
current_info empty;
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index eef7e9a..3e460a9 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -6,8 +6,8 @@ class nofuseException : public std::exception {
} nofuseex;
network::network(const graph& G, cholmod_common* c)
- : G(G), C(G.dual_vertices.size()), fuses(G.edges.size()), backbone(G.edges.size()),
- thresholds(G.edges.size()), px(G, 0, c), py(G, 1, c) {}
+ : px(G, 0, c), py(G, 1, c), G(G), C(G.dual_vertices.size()), fuses(G.edges.size()),
+ backbone(G.edges.size()), thresholds(G.edges.size()) {}
void network::set_thresholds(double beta, std::mt19937& rng) {
if (beta == 0.0) {
@@ -29,25 +29,22 @@ void network::set_thresholds(double beta, std::mt19937& rng) {
}
}
-void network::fracture(hooks& m, bool one_axis) {
+void network::fracture(hooks& m) {
m.pre_fracture(*this);
while (true) {
- double min_cond = 1.0 / G.edges.size();
+ const double min_cond = 1.0 / G.edges.size();
current_info c = this->get_current_info();
- if (px.sol.conductivity[0] > min_cond) {
- this->update_backbone(px.sol.currents);
- } else {
- this->update_backbone(py.sol.currents);
- }
- if (c.conductivity[0] < min_cond && c.conductivity[1] < min_cond) {
+ if (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) {
break;
}
- if (one_axis && (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond)) {
- break;
+ if (px.sol.conductivity[0] > min_cond) {
+ this->update_backbone(px.sol.currents);
+ } else {
+ this->update_backbone(py.sol.currents);
}
unsigned max_pos = UINT_MAX;
@@ -144,6 +141,48 @@ std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edge
return sig;
}
+bool network::get_cycle_helper(std::array<unsigned, 2>& sig, std::set<unsigned>& edges, const std::set<unsigned>& cycle_edges,
+ unsigned vPrev, unsigned vCur, unsigned vEnd) const {
+ for (unsigned ei : G.dual_vertices[vCur].ne) {
+ if (fuses[ei]) {
+ auto it = cycle_edges.find(ei);
+ if (it == cycle_edges.end()) {
+ const std::array<unsigned, 2>& e = G.dual_edges[ei].v;
+ unsigned vn = e[0] == vCur ? e[1] : e[0];
+ if (vn != vPrev) {
+ if (vn == vEnd) {
+ edges.insert(ei);
+ if (G.dual_edges[ei].crossings[0])
+ sig[0]++;
+ if (G.dual_edges[ei].crossings[1])
+ sig[1]++;
+ return true;
+ } else {
+ if (this->get_cycle_helper(sig, edges, cycle_edges, vCur, vn, vEnd)) {
+ edges.insert(ei);
+ if (G.dual_edges[ei].crossings[0])
+ sig[0]++;
+ if (G.dual_edges[ei].crossings[1])
+ sig[1]++;
+ return true;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ return false;
+}
+
+std::pair<std::array<unsigned, 2>, std::set<unsigned>> network::get_cycle(const std::set<unsigned>& cycle_edges, unsigned v0,
+ unsigned v1) const {
+ std::set<unsigned> edges;
+ std::array<unsigned, 2> sig = {0, 0};
+ this->get_cycle_helper(sig, edges, cycle_edges, v0, v0, v1);
+ return {sig, edges};
+}
+
void network::get_cluster_edges_helper(std::set<unsigned>& seen_edges, unsigned v) const {
for (unsigned ei : G.vertices[v].ne) {
if (!backbone[ei]) {
@@ -264,7 +303,7 @@ void network::update_backbone(const std::vector<double>& c) {
return (c[0] % 2 == 0 && c[1] % 2 == 0) ||
((c[0] % 2 == 0 && done_x) || (c[1] % 2 == 0 && done_y));
})) {
- // If it does, we remove it from the backbone!
+ // If none does, we remove it from the backbone!
seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
backbone[i] = true;
@@ -301,6 +340,14 @@ void network::update_backbone(const std::vector<double>& c) {
}
}
}
+ } else if (std::any_of(all_sigs.begin(), all_sigs.end(),
+ [](std::array<unsigned, 2> c) {
+ return ((c[0] % 2 == 0 && c[1] % 2 == 1) || (c[0] % 2 == 1 && c[1] % 2 == 0));
+ })) {
+ // If the bond separates a wrapping path, breaking it would end the
+ // fracture and therefore it will never be removed from the backbone
+ // in this manner. We can skip it in the future.
+ seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
}
}
}
@@ -450,6 +497,11 @@ void network::update_backbone(const std::vector<double>& c) {
}
}
}
+ } else if (std::any_of(all_sigs.begin(), all_sigs.end(),
+ [](std::array<unsigned, 2> c) {
+ return ((c[0] % 2 == 0 && c[1] % 2 == 1) || (c[0] % 2 == 1 && c[1] % 2 == 0));
+ })) {
+ seen_pairs[{d1, d2}] = true;
}
}
}
diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp
index 5f180ee..2965915 100644
--- a/src/analysis_tools.cpp
+++ b/src/analysis_tools.cpp
@@ -1,125 +1,72 @@
#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<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> find_minimal_crack(const Graph& G, const network& n) {
- Graph Gtmp(n.G.vertices.size());
- std::list<unsigned> removed_edges;
-
- class add_tree_edges : public boost::default_dfs_visitor {
- public:
- Graph& G;
- std::list<unsigned>& E;
-
- add_tree_edges(Graph& G, std::list<unsigned>& 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>& E;
- unsigned end;
- struct done{};
-
- find_cycle(std::list<unsigned>& E, unsigned 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>> cycles;
+std::pair<std::array<unsigned, 2>, std::set<unsigned>> find_minimal_crack(const network& n, unsigned i0) {
+ std::set<unsigned> removed_edges = n.get_cycle_edges(n.G.dual_edges[i0].v[0]);
+
+ std::vector<std::pair<std::array<unsigned, 2>, std::set<unsigned>>> cycles;
+ cycles.reserve(removed_edges.size());
+
+ for (unsigned edge : removed_edges) {
+ cycles.push_back(n.get_cycle(removed_edges, n.G.dual_edges[edge].v[0], n.G.dual_edges[edge].v[1]));
+ cycles.back().second.insert(edge);
+ if (n.G.dual_edges[edge].crossings[0])
+ cycles.back().first[0]++;
+ if (n.G.dual_edges[edge].crossings[1])
+ cycles.back().first[1]++;
+ }
- for (auto edge : removed_edges) {
- std::list<unsigned> 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) {
+ return cycles.front();
}
- 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 all_good = true;
+ unsigned not_good;
+ for (unsigned i = 0; i < 2; i++) {
+ if (!((cycles[i].first[0] % 2 == 0 && cycles[i].first[1] % 2 == 1) || (cycles[i].first[0] % 2 == 1 && cycles[i].first[1] % 2 == 0))) {
+ all_good = false;
+ not_good = i;
}
- 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> new_cycle;
- unsigned pos = 0;
- for (uint8_t included : new_bool_cycle) {
- if (included) {
- new_cycle.push_back(pos);
- }
- pos++;
- }
- cycles.push_back(new_cycle);
+ if (all_good) {
+ if (cycles[0].second.size() > cycles[1].second.size()) {
+ return cycles[1];
+ } else {
+ return cycles[0];
}
}
- std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> output;
+ const auto& cycle1 = cycles[0];
+ const auto& cycle2 = cycles[1];
- // find the cycle representing the crack by counting boundary crossings
- for (auto cycle : cycles) {
- std::array<unsigned, 2> crossing_count{0,0};
+ unsigned sum_sig_0 = cycle1.first[0] + cycle2.first[0];
+ unsigned sum_sig_1 = cycle1.first[1] + cycle2.first[1];
- for (auto edge : cycle) {
- if (n.G.dual_edges[edge].crossings[0]) {
- crossing_count[0]++;
- }
- if (n.G.dual_edges[edge].crossings[1]) {
- crossing_count[1]++;
- }
- }
+ std::valarray<uint8_t> bool_cycle1(n.G.edges.size());
+ std::valarray<uint8_t> bool_cycle2(n.G.edges.size());
+ for (auto v : cycle1.second) {
+ bool_cycle1[v] = 1;
+ }
+ for (auto v : cycle2.second) {
+ bool_cycle2[v] = 1;
+ }
- if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
- output.push_back({{1, 0}, cycle});
- } else if (crossing_count[0] % 2 == 0 && crossing_count[1] % 2 == 1) {
- output.push_back({{0, 1}, cycle});
+ std::valarray<uint8_t> new_bool_cycle = bool_cycle1 ^ bool_cycle2;
+
+ std::set<unsigned> new_cycle;
+ unsigned pos = 0;
+ for (uint8_t included : new_bool_cycle) {
+ if (included) {
+ new_cycle.insert(pos);
}
+ pos++;
}
- return output;
+ if (cycles[!not_good].second.size() > new_cycle.size()) {
+ return {{sum_sig_0, sum_sig_1}, new_cycle};
+ } else {
+ return cycles[!not_good];
+ }
}
-*/
diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp
index 2a4c5d2..72d1e82 100644
--- a/src/analysis_tools.hpp
+++ b/src/analysis_tools.hpp
@@ -7,11 +7,5 @@
#include <network.hpp>
-
-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<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> find_minimal_crack(const Graph &, const network &);
+std::pair<std::array<unsigned, 2>, std::set<unsigned>> find_minimal_crack(const network &, unsigned i0);
diff --git a/src/animate.cpp b/src/animate.cpp
index 6c844d2..81bf4ac 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -1,5 +1,6 @@
#include "animate.hpp"
+#include <unistd.h>
#include <iostream>
animate::animate(double Lx, double Ly, unsigned window_size, int argc, char *argv[]) {
@@ -139,65 +140,27 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
}
glEnd();
glFlush();
- if (nw > 2 && !seen_guy) {seen_guy = true; getchar();}
+// if (nw > 2 && !seen_guy) {seen_guy = true; getchar();}
}
void animate::post_fracture(network &n) {
- /*
- std::list<unsigned> crack;
-// unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
- unsigned crack_component = 10000;
-
- std::default_random_engine gen;
- std::uniform_real_distribution<double> dis(0.0,1.0);
-
- auto av_it = avalanches.rbegin();
-
- while (true) {
- for (unsigned e : *av_it) {
- boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
- n.fuses[e] = false;
- }
-
- auto cracks = find_minimal_crack(G, n);
-
- if (cracks.size() == 0) {
- break;
- }
+ auto crack = find_minimal_crack(n, avalanches.back().back());
- av_it++;
- }
-
- std::vector<unsigned> component(boost::num_vertices(G));
- unsigned num = boost::connected_components(G, &component[0]);
-
- std::vector<std::list<unsigned>> components(num);
-
- for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
- components[component[i]].push_back(i);
- }
-
-
- 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 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;
+ for (unsigned j : crack.second) {
+ glBegin(GL_LINES);
+ glColor3d(1.0, 0.0, 0.0);
+ graph::coordinate r1 = n.G.dual_vertices[n.G.dual_edges[j].v[0]].r;
+ graph::coordinate r2 = n.G.dual_vertices[n.G.dual_edges[j].v[1]].r;
- if (n.G.edges[i].crossings[0]) {
+ if (n.G.dual_edges[j].crossings[0]) {
if (r1.x < r2.x) {
r1.x += n.G.L.x;
} else {
r2.x += n.G.L.x;
}
}
- if (n.G.edges[i].crossings[1]) {
+ if (n.G.dual_edges[j].crossings[1]) {
if (r1.y < r2.y) {
r1.y += n.G.L.y;
} else {
@@ -205,82 +168,12 @@ void animate::post_fracture(network &n) {
}
}
- glVertex2d(r1.x, r1.y);
- glVertex2d(r2.x, r2.y);
- }
- }
- glEnd();
-
- switch (key) {
- case 's' :
- for (auto edge : crack) {
- glColor3f(1.0f, 0.0f, 0.0f);
- glBegin(GL_POLYGON);
- for (const graph::coordinate &r : n.G.dual_vertices[n.G.dual_edges[edge].v[0]].polygon) {
- glVertex2d(r.x, r.y);
- }
- glEnd();
- glBegin(GL_POLYGON);
- for (const graph::coordinate &r : n.G.dual_vertices[n.G.dual_edges[edge].v[1]].polygon) {
- glVertex2d(r.x, r.y);
- }
- glEnd();
- }
- glFlush();
- break;
-
- case 'c' :
- for (unsigned 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++) {
- glBegin(GL_POLYGON); // Each set of 3 vertices form a triangle
- for (const graph::coordinate &r: n.G.dual_vertices[*it].polygon) {
- glVertex2d(r.x, r.y);
- }
- glEnd();
- }
- }
- glFlush();
- break;
-
- case 'C' :
- for (unsigned 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++) {
- glBegin(GL_POLYGON); // Each set of 3 vertices form a triangle
- for (const graph::coordinate &r :n.G.dual_vertices[*it].polygon) {
- glVertex2d(r.x, r.y);
- }
- glEnd();
- }
- }
- }
- glFlush();
- break;
-
- case 'm' :
- for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
- glBegin(GL_POLYGON);
- glColor3d(1.0, 0.0, 0.0);
- for (const graph::coordinate &r :n.G.dual_vertices[*it].polygon) {
- glVertex2d(r.x, r.y);
- }
- glEnd();
- }
- glFlush();
- }
+ glVertex2d(r1.x, r1.y);
+ glVertex2d(r2.x, r2.y);
+ glEnd();
}
-*/
+ glFlush();
+
+ getchar();
}
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
index 9d875e6..2b5e55c 100644
--- a/src/animate_fracture.cpp
+++ b/src/animate_fracture.cpp
@@ -57,7 +57,7 @@ int main(int argc, char* argv[]) {
cholmod_common c;
CHOL_F(start)(&c);
- animate meas(sqrt(2*n *a), sqrt( 2*n / a), 700, argc, argv);
+ animate meas(sqrt(2*n *a), sqrt( 2*n / a), 1000, argc, argv);
randutils::auto_seed_128 seeds;
std::mt19937 rng{seeds};
@@ -66,7 +66,7 @@ int main(int argc, char* argv[]) {
graph G(n, a, rng);
elastic_network elastic_network(G, &c);
elastic_network.set_thresholds(beta, rng);
- elastic_network.fracture(meas, true);
+ elastic_network.fracture(meas);
if (quit.load())
break;
diff --git a/src/animate_fracture_square.cpp b/src/animate_fracture_square.cpp
index 12de7da..ad43254 100644
--- a/src/animate_fracture_square.cpp
+++ b/src/animate_fracture_square.cpp
@@ -57,7 +57,7 @@ int main(int argc, char* argv[]) {
cholmod_common c;
CHOL_F(start)(&c);
- animate meas(Lx, Ly, 512, argc, argv);
+ animate meas(Lx, Ly, 1000, argc, argv);
randutils::auto_seed_128 seeds;
std::mt19937 rng{seeds};
@@ -68,7 +68,7 @@ int main(int argc, char* argv[]) {
for (unsigned trial = 0; trial < N; trial++) {
elastic_network tmp_network(perm_network);
tmp_network.set_thresholds(beta, rng);
- tmp_network.fracture(meas, true);
+ tmp_network.fracture(meas);
if (quit.load())
break;