summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/include/network.hpp10
-rw-r--r--lib/src/network.cpp130
-rw-r--r--src/analysis_tools.cpp20
-rw-r--r--src/animate_fracture.cpp28
-rw-r--r--src/fracture.cpp44
5 files changed, 157 insertions, 75 deletions
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index d562530..37f1912 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -31,9 +31,10 @@ class network {
private:
problem px;
problem py;
- std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs;
+ std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs_x;
+ std::unordered_map<std::array<unsigned, 2>, bool> seen_pairs_y;
- void update_backbone(const std::vector<double>& c);
+ void update_backbone(const std::vector<double>& c, const std::array<double, 2>& cond);
void break_edge(unsigned, bool unbreak = false);
void get_cycle_edges_helper(std::set<unsigned>&, std::set<unsigned>&, unsigned, unsigned) const;
bool find_cycle_helper(std::array<unsigned, 2>&, const std::set<unsigned>&, unsigned, unsigned,
@@ -47,6 +48,7 @@ private:
std::list<std::set<unsigned>> get_tie_flaps(unsigned) const;
public:
+ bool two_sides;
const graph& G;
ClusterTree C;
@@ -55,7 +57,7 @@ public:
std::vector<bool> backbone;
std::vector<long double> thresholds;
- network(const graph&, cholmod_common*);
+ network(const graph&, bool two_sides, cholmod_common*);
void set_thresholds(double, std::mt19937&);
void fracture(hooks&);
@@ -76,7 +78,7 @@ private:
double weight;
public:
- fuse_network(const graph&, cholmod_common*, double weight = 1.0);
+ fuse_network(const graph&, cholmod_common*);
fuse_network(const fuse_network&);
current_info get_current_info() override;
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 561796b..8944c15 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -1,14 +1,16 @@
#include "network.hpp"
#include <sstream>
+#include <iostream>
class nofuseException : public std::exception {
virtual const char* what() const throw() { return "No valid fuse was available to break."; }
} nofuseex;
-network::network(const graph& G, cholmod_common* 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()) {}
+network::network(const graph& G, bool two_sides, cholmod_common* c)
+ : px(G, 0, c), py(G, 1, c), two_sides(two_sides), 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) {
@@ -38,14 +40,16 @@ void network::fracture(hooks& m) {
current_info c = this->get_current_info();
- if (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) {
+ if ((c.conductivity[0] < min_cond || c.conductivity[1] < min_cond) && two_sides) {
+ break;
+ } else if (c.conductivity[0] < min_cond && !two_sides) {
break;
}
if (px.sol.conductivity[0] > min_cond) {
- this->update_backbone(px.sol.currents);
+ this->update_backbone(px.sol.currents, c.conductivity);
} else {
- this->update_backbone(py.sol.currents);
+ this->update_backbone(py.sol.currents, c.conductivity);
}
unsigned max_pos = UINT_MAX;
@@ -74,8 +78,8 @@ void network::fracture(hooks& m) {
}
void network::get_cycle_edges_helper(std::set<unsigned>& cycle_edges,
- std::set<unsigned>& seen_vertices, unsigned v_prev,
- unsigned v_cur) const {
+ 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]) {
@@ -103,7 +107,7 @@ 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 {
+ unsigned vPrev, unsigned vCur, unsigned vEnd) const {
for (unsigned ei : G.dual_vertices[vCur].ne) {
if (fuses[ei]) {
if (!cycle_edges.contains(ei)) {
@@ -134,14 +138,15 @@ 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 {
+ unsigned v1) const {
std::array<unsigned, 2> sig = {0, 0};
this->find_cycle_helper(sig, cycle_edges, v0, v0, v1);
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 {
+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]) {
if (!cycle_edges.contains(ei)) {
@@ -173,8 +178,8 @@ bool network::get_cycle_helper(std::array<unsigned, 2>& sig, std::set<unsigned>&
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::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);
@@ -202,7 +207,7 @@ std::set<unsigned> network::get_cluster_edges(unsigned v0) const {
}
void network::get_tie_flaps_helper(std::set<unsigned>& added_edges, unsigned v0,
- unsigned vCur) const {
+ unsigned vCur) const {
for (unsigned ei : G.vertices[vCur].ne) {
if (!backbone[ei]) {
if (!added_edges.contains(ei)) {
@@ -245,14 +250,25 @@ std::list<std::set<unsigned>> network::get_tie_flaps(unsigned v0) const {
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();
+void network::update_backbone(const std::vector<double>& c, const std::array<double, 2>& cc) {
+ bool done_x = cc[0] < 1.0 / G.edges.size();
+ bool done_y = cc[1] < 1.0 / G.edges.size();
+
+ auto cycle_check_current = [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));
+ };
+ auto cycle_check_wrap_x = [](std::array<unsigned, 2> c) {
+ return (c[0] % 2 == 0 && c[1] % 2 == 1);
+ };
+ auto cycle_check_wrap_y = [](std::array<unsigned, 2> c) {
+ return (c[0] % 2 == 1 && c[1] % 2 == 0);
+ };
// 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])) {
- if (!seen_pairs.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]})) {
+ if ((!seen_pairs_x.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]}) && !done_x) || (!seen_pairs_y.contains({G.dual_edges[i].v[0], G.dual_edges[i].v[1]}) && !done_y)) {
// 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.
@@ -261,7 +277,7 @@ void network::update_backbone(const std::vector<double>& c) {
// 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]);
+ this->find_cycle(cedges, G.dual_edges[i].v[0], G.dual_edges[i].v[1]);
if (G.dual_edges[i].crossings[0])
base_sig[0]++;
if (G.dual_edges[i].crossings[1])
@@ -272,11 +288,11 @@ void network::update_backbone(const std::vector<double>& c) {
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]);
+ 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]);
+ 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]};
+ new_sig_1[1] + new_sig_2[1]};
if (G.dual_edges[i].crossings[0])
new_sig[0]++;
@@ -292,13 +308,10 @@ void network::update_backbone(const std::vector<double>& c) {
// 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 (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_current)) {
// If none does, we remove it from the backbone!
- seen_pairs[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
+ seen_pairs_x[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
+ seen_pairs_y[{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
@@ -334,14 +347,13 @@ 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));
- })) {
+ } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_x)) {
// 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;
+ seen_pairs_x[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
+ } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_y)) {
+ seen_pairs_y[{G.dual_edges[i].v[0], G.dual_edges[i].v[1]}] = true;
}
}
}
@@ -371,7 +383,7 @@ void network::update_backbone(const std::vector<double>& c) {
if (C.same_component(G.vertices[v].nd[i], G.vertices[v].nd[l])) {
unsigned il = i < l ? i : l;
unsigned ig = i < l ? l : i;
- if (!seen_pairs.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]})) {
+ if ((!seen_pairs_x.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]}) && !done_x) || (!seen_pairs_y.contains({G.vertices[v].nd[il], G.vertices[v].nd[ig]}) && !done_y)) {
bool any_intervening1 = false;
bool any_intervening2 = false;
unsigned ie1 = 0;
@@ -449,14 +461,10 @@ void network::update_backbone(const std::vector<double>& c) {
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);
- })) {
+ if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_current)) {
// one of our pairs qualifies for trimming!
- seen_pairs[{d1, d2}] = true;
+ seen_pairs_x[{d1, d2}] = true;
+ seen_pairs_y[{d1, d2}] = true;
// We separate the flaps of the bowtie (there may be more than
// two!) into distinct sets.
@@ -490,11 +498,10 @@ 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;
+ } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_x)) {
+ seen_pairs_x[{d1, d2}] = true;
+ } else if (std::any_of(all_sigs.begin(), all_sigs.end(), cycle_check_wrap_y)) {
+ seen_pairs_y[{d1, d2}] = true;
}
}
}
@@ -520,7 +527,7 @@ std::string network::write() {
output += "\"fuses\"->{";
for (unsigned i = 0; i < G.edges.size(); i++) {
if (!fuses[i]) {
- output += std::to_string(i) +",";
+ output += std::to_string(i) + ",";
}
}
output.pop_back();
@@ -536,7 +543,8 @@ std::string network::write() {
output += std::to_string(t) + ",";
}
output.pop_back();
- output += "},\"conductivity\"->{" + std::to_string(c.conductivity[0]) + "," + std::to_string(c.conductivity[1]);
+ output += "},\"conductivity\"->{" + std::to_string(c.conductivity[0]) + "," +
+ std::to_string(c.conductivity[1]);
output += "},\"currents\"->{";
for (const double& t : c.currents) {
output += std::to_string(t) + ",";
@@ -547,35 +555,23 @@ std::string network::write() {
return output;
};
+fuse_network::fuse_network(const graph& G, cholmod_common* c) : network(G, false, c), weight(1.0) {}
-fuse_network::fuse_network(const graph& G, cholmod_common* c, double weight)
- : network(G, c), weight(weight) {}
-
-fuse_network::fuse_network(const fuse_network& n) : network(n), weight(n.weight) {}
+fuse_network::fuse_network(const fuse_network& n) : network(n), weight(1.0) {}
current_info fuse_network::get_current_info() {
px.solve(fuses);
py.solve(fuses);
bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size();
- bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();
current_info ctot;
ctot.currents.resize(G.edges.size());
ctot.conductivity = {px.sol.conductivity[0], py.sol.conductivity[1]};
- if (done_x && !done_y) {
- for (unsigned i = 0; i < G.edges.size(); i++) {
- ctot.currents[i] = fabs(weight * py.sol.currents[i] / py.sol.conductivity[1]);
- }
- } else if (done_y && !done_x) {
- for (unsigned i = 0; i < G.edges.size(); i++) {
- ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0]);
- }
- } else if (!done_x && !done_y) {
+ if (!done_x) {
for (unsigned i = 0; i < G.edges.size(); i++) {
- ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0] +
- weight * py.sol.currents[i] / py.sol.conductivity[1]);
+ ctot.currents[i] = fabs(px.sol.currents[i] / px.sol.conductivity[0]);
}
}
@@ -583,7 +579,7 @@ current_info fuse_network::get_current_info() {
}
elastic_network::elastic_network(const graph& G, cholmod_common* c, double weight)
- : network(G, c), weight(weight) {}
+ : network(G, true, c), weight(weight) {}
elastic_network::elastic_network(const elastic_network& n) : network(n), weight(n.weight) {}
@@ -617,7 +613,7 @@ current_info elastic_network::get_current_info() {
return ctot;
}
-percolation_network::percolation_network(const graph& G, cholmod_common* c) : network(G, c) {}
+percolation_network::percolation_network(const graph& G, cholmod_common* c) : network(G, true, c) {}
percolation_network::percolation_network(const percolation_network& n) : network(n) {}
diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp
index 2965915..8ad8235 100644
--- a/src/analysis_tools.cpp
+++ b/src/analysis_tools.cpp
@@ -30,10 +30,22 @@ std::pair<std::array<unsigned, 2>, std::set<unsigned>> find_minimal_crack(const
}
if (all_good) {
- if (cycles[0].second.size() > cycles[1].second.size()) {
- return cycles[1];
+ if (n.two_sides) {
+ if (cycles[0].second.size() > cycles[1].second.size()) {
+ return cycles[1];
+ } else {
+ return cycles[0];
+ }
} else {
- return cycles[0];
+ if (cycles[0].first[0] % 2 == 0) {
+ return cycles[0];
+ } else {
+ return cycles[1];
+ }
+ }
+ } else if (!n.two_sides) {
+ if (cycles[!not_good].first[0] % 2 == 0) {
+ return cycles[!not_good];
}
}
@@ -63,7 +75,7 @@ std::pair<std::array<unsigned, 2>, std::set<unsigned>> find_minimal_crack(const
pos++;
}
- if (cycles[!not_good].second.size() > new_cycle.size()) {
+ if (cycles[!not_good].second.size() > new_cycle.size() || !n.two_sides) {
return {{sum_sig_0, sum_sig_1}, new_cycle};
} else {
return cycles[!not_good];
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
index e8de8c1..fd1ffe2 100644
--- a/src/animate_fracture.cpp
+++ b/src/animate_fracture.cpp
@@ -69,6 +69,33 @@ int main(int argc, char* argv[]) {
randutils::auto_seed_128 seeds;
std::mt19937 rng{seeds};
+ if (w == 0.0 || w == 1.0) {
+ if (use_aN) {
+ animate meas(sqrt(2*n *a), sqrt( 2*n / a), width, argc, argv);
+
+ for (unsigned trial = 0; trial < N; trial++) {
+ graph G(n, a, rng);
+ fuse_network net(G, &c);
+ net.set_thresholds(beta, rng);
+ net.fracture(meas);
+ }
+ } else {
+ animate meas(Lx, Ly, width, argc, argv);
+ const graph G(Lx, Ly);
+ const fuse_network plain_net(G, &c);
+
+ for (unsigned trial = 0; trial < N; trial++) {
+ fuse_network net = plain_net;
+ net.set_thresholds(beta, rng);
+ try {
+ net.fracture(meas);
+ } catch (std::exception &e) {
+ std::cout << e.what() << std::endl;
+ getchar();
+ }
+ }
+ }
+ } else {
if (use_aN) {
animate meas(sqrt(2*n *a), sqrt( 2*n / a), width, argc, argv);
@@ -89,6 +116,7 @@ int main(int argc, char* argv[]) {
net.fracture(meas);
}
}
+ }
CHOL_F(finish)(&c);
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 96cdac3..22101c8 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -76,6 +76,49 @@ int main(int argc, char* argv[]) {
randutils::auto_seed_128 seeds;
std::mt19937 rng{seeds};
+ if (w == 0 || w == 1) {
+ if (use_aN) {
+ ma meas(n, a, beta, w);
+
+ for (unsigned trial = 0; trial < N; trial++) {
+ while (true) {
+ try {
+ graph G(n, a, rng);
+ fuse_network n(G, &c);
+ n.set_thresholds(beta, rng);
+ n.fracture(meas);
+ break;
+ } catch (std::exception &e) {
+ std::cout << e.what() << '\n';
+ }
+ }
+
+ if (quit.load())
+ break;
+ }
+ } else {
+ ma meas(Lx, Ly, beta, w);
+
+ const graph G(Lx, Ly);
+ const fuse_network n(G, &c);
+
+ for (unsigned trial = 0; trial < N; trial++) {
+ while (true) {
+ try {
+ fuse_network net = n;
+ net.set_thresholds(beta, rng);
+ net.fracture(meas);
+ break;
+ } catch (std::exception &e) {
+ std::cout << e.what() << '\n';
+ }
+ }
+
+ if (quit.load())
+ break;
+ }
+ }
+ } else {
if (use_aN) {
ma meas(n, a, beta, w);
@@ -117,6 +160,7 @@ int main(int argc, char* argv[]) {
break;
}
}
+ }
CHOL_F(finish)(&c);