summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/include/network.hpp56
-rw-r--r--lib/src/graph.cpp22
-rw-r--r--lib/src/network.cpp472
-rw-r--r--src/analysis_tools.cpp2
-rw-r--r--src/analysis_tools.hpp4
-rw-r--r--src/animate.cpp9
-rw-r--r--src/animate.hpp1
-rw-r--r--src/animate_fracture_square.cpp2
-rw-r--r--src/measurements.cpp13
-rw-r--r--src/measurements.hpp2
10 files changed, 258 insertions, 325 deletions
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index bf9015c..8748ea9 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -18,34 +18,46 @@
#define CURRENT_CUTOFF 1e-10
+class fuse_network;
+class elastic_network;
+class percolation_network;
+
class network {
- private:
- problem px;
- problem py;
- std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs;
+ friend class fuse_network;
+ friend class elastic_network;
+ friend class percolation_network;
- 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;
- std::array<unsigned, 2> find_cycle(const std::set<unsigned>&, unsigned, unsigned) const;
+ private:
+ problem px;
+ problem py;
+ std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs;
+
+ 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;
+ 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;
+ void get_tie_flaps_helper(std::set<unsigned>&, unsigned, unsigned) const;
+ std::list<std::set<unsigned>> get_tie_flaps(unsigned) const;
public:
- const graph& G;
+ const graph& G;
- ClusterTree C;
+ ClusterTree C;
- std::vector<bool> fuses;
- std::vector<bool> backbone;
- std::vector<long double> thresholds;
+ std::vector<bool> fuses;
+ std::vector<bool> backbone;
+ std::vector<long double> thresholds;
- network(const graph&, cholmod_common*);
- network(const network&);
+ network(const graph&, cholmod_common*);
- void set_thresholds(double, std::mt19937&);
- void fracture(hooks&, bool one_axis = true);
+ void set_thresholds(double, std::mt19937&);
+ void fracture(hooks&, bool one_axis = true);
- virtual current_info get_current_info() const {current_info empty; return empty;};
+ virtual current_info get_current_info() {current_info empty; return empty;};
};
class fuse_network : public network {
@@ -56,7 +68,7 @@ class fuse_network : public network {
fuse_network(const graph&, cholmod_common*, double weight = 1.0);
fuse_network(const fuse_network&);
- current_info get_current_info() const override;
+ current_info get_current_info() override;
};
class elastic_network : public network {
@@ -67,7 +79,7 @@ class elastic_network : public network {
elastic_network(const graph&, cholmod_common*, double weight = 0.5);
elastic_network(const elastic_network&);
- current_info get_current_info() const override;
+ current_info get_current_info() override;
};
class percolation_network : public network {
@@ -75,6 +87,6 @@ class percolation_network : public network {
percolation_network(const graph&, cholmod_common*);
percolation_network(const percolation_network&);
- current_info get_current_info() const override;
+ current_info get_current_info() override;
};
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index bfdb952..a5063de 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -44,7 +44,7 @@ graph::graph(unsigned Nx, unsigned Ny) {
dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
dual_vertices[i].r.y = (double)(i / (Nx / 2));
- dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx * (i / (Nx / 2)) + (i % (Nx / 2)), (nv + i - Nx / 2) % nv};
+ dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv};
dual_vertices[i].polygon = {
{dual_vertices[i].r.x - 1.0, vertices[i].r.y},
{dual_vertices[i].r.x, vertices[i].r.y - 1.0},
@@ -93,6 +93,16 @@ graph::graph(unsigned Nx, unsigned Ny) {
}
}
+ for (vertex& v : dual_vertices) {
+ v.ne.reserve(v.nd.size());
+ }
+
+ for (unsigned i = 0; i < dual_edges.size(); i++) {
+ for (unsigned vi : dual_edges[i].v) {
+ dual_vertices[vi].ne.push_back(i);
+ }
+ }
+
}
class eulerException: public std::exception
@@ -411,6 +421,16 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
}
}
+ for (vertex& v : dual_vertices) {
+ v.ne.reserve(v.nd.size());
+ }
+
+ for (unsigned i = 0; i < dual_edges.size(); i++) {
+ for (unsigned vi : dual_edges[i].v) {
+ dual_vertices[vi].ne.push_back(i);
+ }
+ }
+
jcv_diagram_free(&diagram);
}
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 2a187c5..0cf50c7 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -77,6 +77,7 @@ void network::fracture(hooks& m, bool one_axis) {
void network::get_cycle_edges_helper(std::set<unsigned>& cycle_edges,
std::set<unsigned>& seen_vertices, unsigned v_prev,
unsigned v_cur) const {
+ seen_vertices.insert(v_cur);
for (unsigned ei : G.dual_vertices[v_cur].ne) {
if (fuses[ei]) {
const std::array<unsigned, 2>& e = G.dual_edges[ei].v;
@@ -84,9 +85,8 @@ void network::get_cycle_edges_helper(std::set<unsigned>& cycle_edges,
if (vn != v_prev) {
auto it = seen_vertices.find(vn);
- if (it == seen_vertices.end()) {
- // if (seen_vertices.contains()) { // need to wait for better c++20 support...
- cycle_edges.insert(vn);
+ if (it != seen_vertices.end()) {
+ cycle_edges.insert(ei);
} else {
this->get_cycle_edges_helper(cycle_edges, seen_vertices, v_cur, vn);
}
@@ -107,16 +107,21 @@ std::set<unsigned> network::get_cycle_edges(unsigned v0) const {
bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<unsigned>& cycle_edges, unsigned vPrev, unsigned vCur, unsigned vEnd) const {
for (unsigned ei : G.dual_vertices[vCur].ne) {
if (fuses[ei]) {
- 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) {
- return true;
- } else {
- if (find_cycle_helper(sig, cycle_edges, vCur, vn, vEnd)) {
+ 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) {
if (G.dual_edges[ei].crossings[0]) sig[0]++;
if (G.dual_edges[ei].crossings[1]) sig[1]++;
return true;
+ } else {
+ if (this->find_cycle_helper(sig, cycle_edges, vCur, vn, vEnd)) {
+ if (G.dual_edges[ei].crossings[0]) sig[0]++;
+ if (G.dual_edges[ei].crossings[1]) sig[1]++;
+ return true;
+ }
}
}
}
@@ -127,202 +132,169 @@ bool network::find_cycle_helper(std::array<unsigned, 2>& sig, const std::set<uns
}
std::array<unsigned, 2> network::find_cycle(const std::set<unsigned>& cycle_edges, unsigned v0, unsigned v1) const {
- std::array<unsigned, 2> sig;
- this->find_cycle_helper(sig, cycle_edges, v0, v0,
+ std::array<unsigned, 2> sig = {0, 0};
+ this->find_cycle_helper(sig, cycle_edges, v0, v0, v1);
+ return sig;
}
-void network::update_backbone(const std::vector<double>& c) {
- bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size();
- bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();
-
- class find_cycle : public boost::default_dfs_visitor {
- public:
- const graph& G;
- std::array<unsigned, 2>& E;
- unsigned end;
- struct done {};
-
- find_cycle(const graph& G, std::array<unsigned, 2>& E, unsigned end) : G(G), E(E), end(end) {}
+void network::get_cluster_edges_helper(std::set<unsigned>& seen_edges, unsigned v) const {
+ for (unsigned ei : G.vertices[v].ne) {
+ if (!backbone[ei]) {
+ auto it = seen_edges.find(ei);
+ if (it == seen_edges.end()) {
+ const std::array<unsigned, 2>& e = G.edges[ei].v;
+ unsigned vn = e[0] == v ? e[1] : e[0];
- void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
- if (v == end) {
- throw done{};
+ seen_edges.insert(ei);
+ this->get_cluster_edges_helper(seen_edges, vn);
}
}
+ }
+}
- void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
- unsigned e = g[ed].index;
- if (G.dual_edges[e].crossings[0])
- E[0]++;
- if (G.dual_edges[e].crossings[1])
- E[1]++;
- }
+std::set<unsigned> network::get_cluster_edges(unsigned v0) const {
+ std::set<unsigned> cluster_edges;
+ this->get_cluster_edges_helper(cluster_edges, v0);
+ return cluster_edges;
+}
- void finish_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
- unsigned e = g[ed].index;
- if (G.dual_edges[e].crossings[0])
- E[0]--;
- if (G.dual_edges[e].crossings[1])
- E[1]--;
+void network::get_tie_flaps_helper(std::set<unsigned>& added_edges, unsigned v0, unsigned vCur) const {
+ for (unsigned ei : G.vertices[vCur].ne) {
+ if (!backbone[ei]) {
+ auto it = added_edges.find(ei);
+ if (it == added_edges.end()) {
+ const std::array<unsigned, 2>& e = G.edges[ei].v;
+ unsigned vn = e[0] == vCur ? e[1] : e[0];
+
+ if (vn != v0) {
+ added_edges.insert(ei);
+ this->get_tie_flaps_helper(added_edges, v0, vn);
+ }
+ }
}
- };
+ }
+}
- class get_cluster_edges : public boost::default_dfs_visitor {
- public:
- std::set<unsigned>& E;
+std::list<std::set<unsigned>> network::get_tie_flaps(unsigned v0) const {
+ std::list<std::set<unsigned>> tie_flaps;
+
+ for (unsigned ei : G.vertices[v0].ne) {
+ if (!backbone[ei]) {
+ bool seen_edge = false;
+ for (const std::set<unsigned>& flap : tie_flaps) {
+ auto it = flap.find(ei);
+ if (it != flap.end()) {
+ seen_edge = true;
+ break;
+ }
+ }
- get_cluster_edges(std::set<unsigned>& E) : E(E) {}
+ if (!seen_edge) {
+ tie_flaps.push_back({ei});
- void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
- unsigned e = g[ed].index;
- E.insert(e);
+ const std::array<unsigned, 2>& e = G.edges[ei].v;
+ unsigned vn = e[0] == v0 ? e[1] : e[0];
+
+ this->get_tie_flaps_helper(tie_flaps.back(), v0, vn);
+ }
}
- };
+ }
+
+ return tie_flaps;
+}
+
+void network::update_backbone(const std::vector<double>& c) {
+ bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size();
+ bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();
+ // First, we will check for lollipops!
for (unsigned i = 0; i < G.edges.size(); i++) {
if ((!backbone[i]) && C.same_component(G.dual_edges[i].v[0], G.dual_edges[i].v[1])) {
auto it_search = seen_pairs.find({G.dual_edges[i].v[0], G.dual_edges[i].v[1]});
if (it_search == seen_pairs.end()) {
- std::list<unsigned> cedges = this->get_cycle_edges(G.dual_edges[i].v[0], G.dual_edges[i].v[1]);
-
- std::array<unsigned, 2> cycle = {0, 0};
+ // This is a candidate lollipop stem. First, we will identify any
+ // cycles in the dual cluster that impinges on the stem and mark them
+ // by an edge that uniquely severs each.
+ std::set<unsigned> cedges = this->get_cycle_edges(G.dual_edges[i].v[0]);
+
+ // Now, we create a crossing signature for each cycle. First, we do it
+ // for the new cycle that would form by removing the stem.
+ std::array<unsigned, 2> base_sig = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[i].v[1]);
if (G.dual_edges[i].crossings[0])
- cycle[0]++;
+ base_sig[0]++;
if (G.dual_edges[i].crossings[1])
- cycle[1]++;
- find_cycle vis2(G, cycle, G.dual_edges[i].v[1]);
- std::vector<boost::default_color_type> cm2(G.dual_vertices.size());
- try {
- boost::depth_first_visit(dG, G.dual_edges[i].v[0], vis2,
- boost::make_iterator_property_map(
- cm2.begin(), boost::get(boost::vertex_index, dG), cm2[0]));
- } catch (find_cycle::done const&) {
- std::list<std::array<unsigned, 2>> other_cycles = {cycle};
- for (unsigned e : cedges) {
- std::array<unsigned, 2> tcycle = {0, 0};
- if (G.dual_edges[i].crossings[0])
- tcycle[0]++;
- if (G.dual_edges[i].crossings[1])
- tcycle[1]++;
- if (G.dual_edges[e].crossings[0])
- tcycle[0]++;
- if (G.dual_edges[e].crossings[1])
- tcycle[1]++;
- find_cycle vis3(G, tcycle, G.dual_edges[e].v[0]);
- std::vector<boost::default_color_type> cm3(G.dual_vertices.size());
- try {
- boost::depth_first_visit(
- dG, G.dual_edges[i].v[0], vis3,
- boost::make_iterator_property_map(cm3.begin(),
- boost::get(boost::vertex_index, dG), cm3[0]));
- } catch (find_cycle::done const&) {
- }
- find_cycle vis4(G, tcycle, G.dual_edges[e].v[1]);
- std::vector<boost::default_color_type> cm4(G.dual_vertices.size());
- try {
- boost::depth_first_visit(
- dG, G.dual_edges[i].v[1], vis4,
- boost::make_iterator_property_map(cm4.begin(),
- boost::get(boost::vertex_index, dG), cm4[0]));
- } catch (find_cycle::done const&) {
- }
- other_cycles.push_back(tcycle);
- }
+ base_sig[1]++;
- if (std::any_of(other_cycles.begin(), other_cycles.end(),
- [done_x, done_y](std::array<unsigned, 2> c) {
- return (c[0] % 2 == 0 && c[1] % 2 == 0) ||
- ((c[0] % 2 == 0 && done_x) || (c[1] % 2 == 0 && done_y));
- })) {
- seen_pairs[{G.edges[i].v[0], G.edges[i].v[1]}] = true;
- backbone[i] = true;
- boost::remove_edge(G.edges[i].v[0], G.edges[i].v[1], bG);
-
- for (unsigned j = 0; j < 2; j++) {
- std::set<unsigned> component_edges = {};
- get_cluster_edges visc(component_edges);
- std::vector<boost::default_color_type> cm2(G.vertices.size());
- boost::depth_first_visit(
- bG, G.edges[i].v[j], visc,
- boost::make_iterator_property_map(cm2.begin(),
- boost::get(boost::vertex_index, bG), cm2[0]));
-
- std::array<double, 2> total_cur = {0.0, 0.0};
- for (unsigned e : component_edges) {
- if (G.edges[e].crossings[0]) {
- if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) {
- total_cur[0] += c[e];
- } else {
- total_cur[0] -= c[e];
- }
+ // Then, we do it for each of the existing cycles we found in the
+ // previous step.
+ std::list<std::array<unsigned, 2>> all_sigs = {base_sig};
+ for (unsigned e : cedges) {
+ std::array<unsigned, 2> new_sig_1 = this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[e].v[0]);
+ std::array<unsigned, 2> new_sig_2 = this->find_cycle(cedges, G.dual_edges[i].v[1], G.dual_edges[e].v[1]);
+ std::array<unsigned, 2> new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]};
+
+ if (G.dual_edges[i].crossings[0])
+ new_sig[0]++;
+ if (G.dual_edges[i].crossings[1])
+ new_sig[1]++;
+ if (G.dual_edges[e].crossings[0])
+ new_sig[0]++;
+ if (G.dual_edges[e].crossings[1])
+ new_sig[1]++;
+
+ all_sigs.push_back(new_sig);
+ }
+
+ // Now, having found all cycles involving the candidate stem, we check
+ // if any of them spans the torus and therefore can carry current.
+ if (std::any_of(all_sigs.begin(), all_sigs.end(),
+ [done_x, done_y](std::array<unsigned, 2> 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!
+ seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
+ backbone[i] = true;
+
+ // We're not done yet: we've severed the stem but the pop remains! We
+ // check each side of the lollipop and sum up all of the currents
+ // that span an edge. Any component without a sufficiently large net
+ // current must be disconnected from the current-carrying cluster and
+ // can also be removed from the backbone.
+ for (unsigned j = 0; j < 2; j++) {
+ std::set<unsigned> cluster_edges = this->get_cluster_edges(G.edges[i].v[j]);
+ std::array<double, 2> total_current = {0.0, 0.0};
+
+ for (unsigned e : cluster_edges) {
+ if (G.edges[e].crossings[0]) {
+ if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) {
+ total_current[0] += c[e];
+ } else {
+ total_current[0] -= c[e];
}
- if (G.edges[e].crossings[1]) {
- if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) {
- total_cur[1] += c[e];
- } else {
- total_cur[1] -= c[e];
- }
+ }
+ if (G.edges[e].crossings[1]) {
+ if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) {
+ total_current[1] += c[e];
+ } else {
+ total_current[1] -= c[e];
}
}
+ }
- if (fabs(total_cur[0]) < 1.0 / G.edges.size() &&
- fabs(total_cur[1]) < 1.0 / G.edges.size()) {
- for (unsigned e : component_edges) {
- backbone[e] = true;
- boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
- }
+ if (fabs(total_current[0]) < 1.0 / G.edges.size() &&
+ fabs(total_current[1]) < 1.0 / G.edges.size()) {
+ for (unsigned e : cluster_edges) {
+ backbone[e] = true;
}
}
}
}
- for (unsigned e : cedges) {
- boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
- }
}
}
}
- class get_tie_edges : public boost::default_dfs_visitor {
- public:
- std::set<unsigned>& LE;
- std::set<unsigned>& RE;
- unsigned start_edge;
- bool initialized;
- bool done_left;
- unsigned disc;
- unsigned dr;
-
- get_tie_edges(std::set<unsigned>& LE, std::set<unsigned>& RE) : LE(LE), RE(RE) {
- initialized = false;
- done_left = false;
- disc = 0;
- dr = 0;
- }
-
- void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor vd, const Graph& g) {
- if (!initialized && disc > 0) {
- start_edge = vd;
- initialized = true;
- }
- disc++;
- }
-
- void finish_vertex(boost::graph_traits<Graph>::vertex_descriptor vd, const Graph& g) {
- if (vd == start_edge) {
- done_left = true;
- }
- }
-
- void examine_edge(boost::graph_traits<Graph>::edge_descriptor ed, const Graph& g) {
- unsigned e = g[ed].index;
- if (!done_left) {
- LE.insert(e);
- } else if (LE.find(e) == LE.end()) {
- RE.insert(e);
- }
- }
- };
-
+ // Now, we will check for bow ties!
std::set<unsigned> bb_verts;
for (unsigned i = 0; i < G.edges.size(); i++) {
@@ -385,139 +357,75 @@ void network::update_backbone(const std::vector<double>& c) {
d1 = G.vertices[v].nd[il];
d2 = G.vertices[v].nd[ig];
- std::list<unsigned> cedges = {};
- get_cycle_edges vis(cedges);
- std::vector<boost::default_color_type> cm1(G.dual_vertices.size());
- boost::depth_first_visit(
- dG, d1, vis,
- boost::make_iterator_property_map(cm1.begin(),
- boost::get(boost::vertex_index, dG), cm1[0]));
-
- for (unsigned e : cedges) {
- boost::remove_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], dG);
- }
+ std::set<unsigned> cedges = this->get_cycle_edges(d1);
- std::array<unsigned, 2> cycle = {0, 0};
+ std::array<unsigned, 2> base_sig = this->find_cycle(cedges, d1, d2);
for (unsigned k = p1; k < p2; k++) {
if (G.dual_edges[G.vertices[v].ne[k]].crossings[0])
- cycle[0]++;
+ base_sig[0]++;
if (G.dual_edges[G.vertices[v].ne[k]].crossings[1])
- cycle[1]++;
+ base_sig[1]++;
}
- find_cycle vis2(G, cycle, d2);
- std::vector<boost::default_color_type> cm2(G.dual_vertices.size());
- try {
- boost::depth_first_visit(
- dG, d1, vis2,
- boost::make_iterator_property_map(cm2.begin(),
- boost::get(boost::vertex_index, dG), cm2[0]));
- } catch (find_cycle::done const&) {
- std::list<std::array<unsigned, 2>> other_cycles = {cycle};
- for (unsigned e : cedges) {
- std::array<unsigned, 2> tcycle = {0, 0};
- for (unsigned k = p1; k < p2; k++) {
- if (G.dual_edges[G.vertices[v].ne[k]].crossings[0])
- tcycle[0]++;
- if (G.dual_edges[G.vertices[v].ne[k]].crossings[1])
- tcycle[1]++;
- }
- if (G.dual_edges[e].crossings[0])
- tcycle[0]++;
- if (G.dual_edges[e].crossings[1])
- tcycle[1]++;
- find_cycle vis3(G, tcycle, G.dual_edges[e].v[0]);
- std::vector<boost::default_color_type> cm3(G.dual_vertices.size());
- try {
- boost::depth_first_visit(
- dG, d1, vis3,
- boost::make_iterator_property_map(
- cm3.begin(), boost::get(boost::vertex_index, dG), cm3[0]));
- } catch (find_cycle::done const&) {
- }
- find_cycle vis4(G, tcycle, G.dual_edges[e].v[1]);
- std::vector<boost::default_color_type> cm4(G.dual_vertices.size());
- try {
- boost::depth_first_visit(
- dG, d2, vis4,
- boost::make_iterator_property_map(
- cm4.begin(), boost::get(boost::vertex_index, dG), cm4[0]));
- } catch (find_cycle::done const&) {
- }
- other_cycles.push_back(tcycle);
+
+ // Then, we do it for each of the existing cycles we found in the
+ // previous step.
+ std::list<std::array<unsigned, 2>> all_sigs = {base_sig};
+ for (unsigned e : cedges) {
+ std::array<unsigned, 2> new_sig_1 = this->find_cycle(cedges, d1, G.dual_edges[e].v[0]);
+ std::array<unsigned, 2> new_sig_2 = this->find_cycle(cedges, d2, G.dual_edges[e].v[1]);
+ std::array<unsigned, 2> new_sig = {new_sig_1[0] + new_sig_2[0], new_sig_1[1] + new_sig_2[1]};
+
+ for (unsigned k = p1; k < p2; k++) {
+ if (G.dual_edges[G.vertices[v].ne[k]].crossings[0])
+ new_sig[0]++;
+ if (G.dual_edges[G.vertices[v].ne[k]].crossings[1])
+ new_sig[1]++;
}
+ if (G.dual_edges[e].crossings[0])
+ new_sig[0]++;
+ if (G.dual_edges[e].crossings[1])
+ new_sig[1]++;
- if (std::any_of(other_cycles.begin(), other_cycles.end(),
- [done_x, done_y](std::array<unsigned, 2> c) {
- return ((c[0] % 2 == 0 && c[1] % 2 == 0) ||
- (c[0] % 2 == 0 && done_x)) ||
- (c[1] % 2 == 0 && done_y);
- })) {
- seen_pairs[{d1, d2}] = true;
-
- std::set<unsigned> left_edges = {};
- std::set<unsigned> right_edges = {};
- get_tie_edges visc(left_edges, right_edges);
- std::vector<boost::default_color_type> cm2(G.vertices.size());
- boost::depth_first_visit(
- bG, v, visc,
- boost::make_iterator_property_map(
- cm2.begin(), boost::get(boost::vertex_index, bG), cm2[0]));
-
- std::array<double, 2> total_cur_left = {0.0, 0.0};
- for (unsigned e : left_edges) {
- if (G.edges[e].crossings[0]) {
- if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) {
- total_cur_left[0] += c[e];
- } else {
- total_cur_left[0] -= c[e];
- }
- }
- if (G.edges[e].crossings[1]) {
- if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) {
- total_cur_left[1] += c[e];
- } else {
- total_cur_left[1] -= c[e];
- }
- }
- }
+ all_sigs.push_back(new_sig);
+ }
+
+ if (std::any_of(all_sigs.begin(), all_sigs.end(),
+ [done_x, done_y](std::array<unsigned, 2> c) {
+ return ((c[0] % 2 == 0 && c[1] % 2 == 0) ||
+ (c[0] % 2 == 0 && done_x)) ||
+ (c[1] % 2 == 0 && done_y);
+ })) {
+ seen_pairs[{d1, d2}] = true;
- std::array<double, 2> total_cur_right = {0.0, 0.0};
- for (unsigned e : right_edges) {
+ std::list<std::set<unsigned>> flaps = this->get_tie_flaps(v);
+
+ for (const std::set<unsigned>& flap : flaps) {
+ std::array<double, 2> total_current = {0.0, 0.0};
+ for (unsigned e : flap) {
if (G.edges[e].crossings[0]) {
if (G.vertices[G.edges[e].v[0]].r.x < G.vertices[G.edges[e].v[1]].r.x) {
- total_cur_right[0] += c[e];
+ total_current[0] += c[e];
} else {
- total_cur_right[0] -= c[e];
+ total_current[0] -= c[e];
}
}
if (G.edges[e].crossings[1]) {
if (G.vertices[G.edges[e].v[0]].r.y < G.vertices[G.edges[e].v[1]].r.y) {
- total_cur_right[1] += c[e];
+ total_current[1] += c[e];
} else {
- total_cur_right[1] -= c[e];
+ total_current[1] -= c[e];
}
}
}
- if (fabs(total_cur_left[0]) < 1.0 / G.edges.size() &&
- fabs(total_cur_left[1]) < 1.0 / G.edges.size()) {
- for (unsigned e : left_edges) {
+ if (fabs(total_current[0]) < 1.0 / G.edges.size() &&
+ fabs(total_current[1]) < 1.0 / G.edges.size()) {
+ for (unsigned e : flap) {
backbone[e] = true;
- boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
- }
- }
- if (fabs(total_cur_right[0]) < 1.0 / G.edges.size() &&
- fabs(total_cur_right[1]) < 1.0 / G.edges.size()) {
- for (unsigned e : right_edges) {
- backbone[e] = true;
- boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG);
}
}
}
}
- for (unsigned e : cedges) {
- boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG);
- }
}
}
}
diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp
index 2ef74c6..5f180ee 100644
--- a/src/analysis_tools.cpp
+++ b/src/analysis_tools.cpp
@@ -6,6 +6,7 @@ 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;
}
@@ -120,4 +121,5 @@ std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> find_minimal_
return output;
}
+*/
diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp
index c3bffa2..2a4c5d2 100644
--- a/src/analysis_tools.hpp
+++ b/src/analysis_tools.hpp
@@ -11,7 +11,7 @@
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>);
+//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::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> find_minimal_crack(const Graph &, const network &);
diff --git a/src/animate.cpp b/src/animate.cpp
index a7b6173..6c844d2 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -2,7 +2,7 @@
#include "animate.hpp"
#include <iostream>
-animate::animate(double Lx, double Ly, unsigned window_size, int argc, char *argv[]) : G(2 * (unsigned)ceil(Lx * Ly / 2)) {
+animate::animate(double Lx, double Ly, unsigned window_size, int argc, char *argv[]) {
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
glutInitWindowSize((unsigned)(Lx / Ly * window_size), window_size);
@@ -16,7 +16,6 @@ animate::animate(double Lx, double Ly, unsigned window_size, int argc, char *arg
void animate::pre_fracture(const network &n) {
lv = std::numeric_limits<long double>::lowest();
avalanches = {};
- boost::remove_edge_if(trivial, G);
seen_guy = false;
glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
@@ -59,7 +58,6 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
avalanches.back().push_back(i);
}
- boost::add_edge(n.G.dual_edges[i].v[0], n.G.dual_edges[i].v[1], {i}, G);
glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
glClear(GL_COLOR_BUFFER_BIT);
@@ -110,6 +108,11 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
draw = true;
}
+ if (cur.currents[j] > 1e-9 && n.backbone[j]) {
+ weird = true;
+ }
+
+
if (draw) {
graph::coordinate r1 = n.G.vertices[n.G.edges[j].v[0]].r;
graph::coordinate r2 = n.G.vertices[n.G.edges[j].v[1]].r;
diff --git a/src/animate.hpp b/src/animate.hpp
index 865b3bd..d14d755 100644
--- a/src/animate.hpp
+++ b/src/animate.hpp
@@ -9,7 +9,6 @@
class animate : public hooks {
private:
- Graph G;
bool seen_guy;
public:
long double lv;
diff --git a/src/animate_fracture_square.cpp b/src/animate_fracture_square.cpp
index bc4c3a2..12de7da 100644
--- a/src/animate_fracture_square.cpp
+++ b/src/animate_fracture_square.cpp
@@ -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, false);
+ tmp_network.fracture(meas, true);
if (quit.load())
break;
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 984b80c..3004d7e 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -144,7 +144,6 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u
}
ma::ma(unsigned n, double a, double beta, double weight, bool one) :
- G(2 * n),
sn(2 * n),
ss(2 * n),
sm(2 * n),
@@ -179,7 +178,6 @@ ma::ma(unsigned n, double a, double beta, double weight, bool one) :
}
ma::ma(unsigned Lx, unsigned Ly, double beta, double weight, bool one) :
- G(Lx * Ly / 2),
sn(Lx * Ly / 2),
ss(Lx * Ly / 2),
sm(Lx * Ly / 2),
@@ -228,7 +226,6 @@ ma::~ma() {
void ma::pre_fracture(const network&) {
lv = std::numeric_limits<long double>::lowest();
- boost::remove_edge_if(trivial, G);
avalanches = {};
num = 0;
}
@@ -250,14 +247,11 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
}
}
- boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
num++;
}
void ma::post_fracture(network &n) {
- auto post_cracks = find_minimal_crack(G, n);
- std::vector<unsigned> component(boost::num_vertices(G));
- unsigned num = boost::connected_components(G, &component[0]);
+/* auto post_cracks = find_minimal_crack(G, n);
if (post_cracks.size() > 2 || post_cracks.size() == 0) {
throw badcycleex;
}
@@ -271,10 +265,6 @@ void ma::post_fracture(network &n) {
autocorrelation2(n.G.L.x, n.G.L.y, 2 * sqrt(cl.size()), 2 * sqrt(cl.size()), cl, cl_cs, c.first);
}
- unsigned crack_component = component[n.G.dual_edges[post_cracks.front().second.front()].v[0]];
-
- std::vector<std::list<graph::coordinate>> components(num);
-
for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) {
components[component[i]].push_back(n.G.dual_vertices[i].r);
}
@@ -336,6 +326,7 @@ void ma::post_fracture(network &n) {
}
autocorrelation2(n.G.L.x, n.G.L.y, 2 * sqrt(cA.size()), 2 * sqrt(cA.size()), cA, cA_co, post_cracks.front().first);
+*/
sd[num - 1]++;
}
diff --git a/src/measurements.hpp b/src/measurements.hpp
index b22c327..879f511 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -19,8 +19,6 @@ class ma : public hooks {
// - interface for turning on and off specific measurements
//
private:
- Graph G;
-
unsigned num;
// measurement storage