From 948f90b6493da83d10e18f30b0fbb8e937e29c0b Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 24 Jun 2019 21:41:34 -0400 Subject: mostly implemented the ability to find dead bonds using topological properties --- lib/include/graph.hpp | 2 + lib/include/network.hpp | 46 ++-- lib/include/problem.hpp | 4 +- lib/src/graph.cpp | 54 ++++- lib/src/network.cpp | 479 +++++++++++++++++++++++++++++++++++----- lib/src/problem.cpp | 26 +-- src/analysis_tools.hpp | 14 -- src/animate.cpp | 62 +++++- src/animate_fracture.cpp | 20 +- src/animate_fracture_square.cpp | 2 +- src/fracture.cpp | 14 +- src/measurements.cpp | 21 +- src/measurements.hpp | 1 + 13 files changed, 622 insertions(+), 123 deletions(-) diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp index 45cc478..0452339 100644 --- a/lib/include/graph.hpp +++ b/lib/include/graph.hpp @@ -21,6 +21,8 @@ class graph { typedef struct vertex { coordinate r; + std::vector nd; + std::vector ne; std::vector polygon; } vertex; diff --git a/lib/include/network.hpp b/lib/include/network.hpp index 544f828..8211cfc 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include "problem.hpp" #include "graph.hpp" @@ -16,60 +17,79 @@ #define CURRENT_CUTOFF 1e-10 +#include +#include +#include +#include +#include +#include +#include +#include + +struct EdgeProperties { + unsigned index; +}; + +typedef boost::adjacency_list Graph; +typedef boost::graph_traits::vertex_descriptor Vertex; +typedef boost::graph_traits::vertices_size_type VertexIndex; + class network { public: const graph& G; + + Graph bG; + Graph dG; + std::vector rank; + std::vector parent; + boost::disjoint_sets ds; + std::unordered_map, bool> seen_pairs; + std::vector fuses; + std::vector backbone; std::vector thresholds; - network(const graph&); + problem px; + problem py; + + network(const graph&, cholmod_common*); network(const network&); void set_thresholds(double, std::mt19937&); void fracture(hooks&, bool one_axis = false); + void update_backbone(const std::vector& c); - virtual void break_edge(unsigned e, bool unbreak = false) {fuses[e] = !unbreak;}; + void break_edge(unsigned, bool unbreak = false); virtual current_info get_current_info() {current_info empty; return empty;}; }; class fuse_network : public network { private: double weight; - problem ohm_x; - problem ohm_y; public: fuse_network(const graph&, cholmod_common*, double weight = 1.0); fuse_network(const fuse_network&); - void break_edge(unsigned, bool unbreak = false) override; current_info get_current_info() override; }; class elastic_network : public network { private: double weight; - problem hook_x; - problem hook_y; public: elastic_network(const graph&, cholmod_common*, double weight = 0.5); elastic_network(const elastic_network&); - void break_edge(unsigned, bool unbreak = false) override; current_info get_current_info() override; }; class percolation_network : public network { - private: - problem px; - problem py; - public: percolation_network(const graph&, cholmod_common*); percolation_network(const percolation_network&); - void break_edge(unsigned, bool unbreak = false) override; current_info get_current_info() override; }; diff --git a/lib/include/problem.hpp b/lib/include/problem.hpp index c2c9e09..740751d 100644 --- a/lib/include/problem.hpp +++ b/lib/include/problem.hpp @@ -28,11 +28,13 @@ class problem { cholmod_common* c; public: + current_info sol; + problem(const graph&, unsigned, cholmod_common*); problem(const graph&, unsigned, cholmod_sparse*, cholmod_common*); problem(const problem&); ~problem(); - current_info solve(std::vector& fuses); + void solve(std::vector& fuses); void break_edge(unsigned, bool unbreak = false); }; diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp index 974d3b3..bfdb952 100644 --- a/lib/src/graph.cpp +++ b/lib/src/graph.cpp @@ -33,6 +33,8 @@ graph::graph(unsigned Nx, unsigned Ny) { for (unsigned i = 0; i < nv; i++) { vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))); vertices[i].r.y = (double)(i / (Nx / 2)); + signed f = (1 + i / (Nx / 2)) % 2 == 1 ? 1 : -1; + vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv}; vertices[i].polygon = { {vertices[i].r.x - 1.0, vertices[i].r.y}, {vertices[i].r.x, vertices[i].r.y - 1.0}, @@ -42,6 +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].polygon = { {dual_vertices[i].r.x - 1.0, vertices[i].r.y}, {dual_vertices[i].r.x, vertices[i].r.y - 1.0}, @@ -67,6 +70,29 @@ graph::graph(unsigned Nx, unsigned Ny) { } } + for (vertex& v : vertices) { + v.ne.resize(v.nd.size()); + } + + for (unsigned i = 0; i < edges.size(); i++) { + for (unsigned vi : edges[i].v) { + auto it1 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[0]); + auto it2 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[1]); + + unsigned d1 = std::distance(vertices[vi].nd.begin(), it1); + unsigned d2 = std::distance(vertices[vi].nd.begin(), it2); + + unsigned dv1 = d1 < d2 ? d1 : d2; + unsigned dv2 = d1 < d2 ? d2 : d1; + + if (dv2 - dv1 == 1) { + vertices[vi].ne[dv1] = i; + } else { + vertices[vi].ne[dv2] = i; + } + } + } + } class eulerException: public std::exception @@ -239,12 +265,15 @@ void graph::helper(unsigned nv, std::mt19937& rng) { if (it1 == known_vertices.end()) { vi1 = dual_vertices.size(); - dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{{site->p.x, site->p.y}}}); + dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{i1},{},{{site->p.x, site->p.y}}}); known_vertices[t1] = vi1; } else { vi1 = it1->second; + dual_vertices[vi1].nd.push_back(i1); dual_vertices[vi1].polygon.push_back({site->p.x, site->p.y}); } + + vertices[i1].nd.push_back(vi1); if (i1 < i2 || (i1 == i2 && !self_bonded)) { if (i1 == i2) { @@ -359,6 +388,29 @@ void graph::helper(unsigned nv, std::mt19937& rng) { } } + for (vertex& v : vertices) { + v.ne.resize(v.nd.size()); + } + + for (unsigned i = 0; i < edges.size(); i++) { + for (unsigned vi : edges[i].v) { + auto it1 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[0]); + auto it2 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[1]); + + unsigned d1 = std::distance(vertices[vi].nd.begin(), it1); + unsigned d2 = std::distance(vertices[vi].nd.begin(), it2); + + unsigned dv1 = d1 < d2 ? d1 : d2; + unsigned dv2 = d1 < d2 ? d2 : d1; + + if (dv2 - dv1 == 1) { + vertices[vi].ne[dv1] = i; + } else { + vertices[vi].ne[dv2] = i; + } + } + } + jcv_diagram_free(&diagram); } diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 3424b94..dcc5e0e 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -9,9 +9,20 @@ class nofuseException: public std::exception } } nofuseex; -network::network(const graph& G) : G(G), fuses(G.edges.size()), thresholds(G.edges.size()) {} +network::network(const graph& G, cholmod_common* c) : + G(G), bG(G.vertices.size()), dG(G.dual_vertices.size()), + rank(G.dual_vertices.size()), parent(G.dual_vertices.size()), ds(&rank[0], &parent[0]), + fuses(G.edges.size()), backbone(G.edges.size()), thresholds(G.edges.size()), + px(G, 0, c), py(G, 1, c) { + for (unsigned i = 0; i < G.edges.size(); i++) { + boost::add_edge(G.edges[i].v[0], G.edges[i].v[1], {i}, bG); + } + initialize_incremental_components(dG, ds); +} -network::network(const network& n) : G(n.G), fuses(n.fuses), thresholds(n.thresholds) {} +network::network(const network& n) : G(n.G), bG(n.bG), dG(n.dG), rank(n.rank), parent(n.parent), ds(&rank[0], &parent[0]), fuses(n.fuses), backbone(n.backbone), thresholds(n.thresholds), px(n.px), py(n.py) { + initialize_incremental_components(dG, ds); +} void network::set_thresholds(double beta, std::mt19937& rng) { if (beta == 0.0) { @@ -36,10 +47,15 @@ void network::fracture(hooks& m, bool one_axis) { m.pre_fracture(*this); while (true) { - current_info c = this->get_current_info(); - 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) { break; } @@ -52,7 +68,7 @@ void network::fracture(hooks& m, bool one_axis) { long double max_val = std::numeric_limits::lowest(); for (unsigned i = 0; i < G.edges.size(); i++) { - if (!fuses[i] && c.currents[i] > CURRENT_CUTOFF) { + if (!backbone[i]) { long double val = logl(c.currents[i]) - thresholds[i]; if (val > max_val) { @@ -73,43 +89,413 @@ void network::fracture(hooks& m, bool one_axis) { m.post_fracture(*this); } +void network::update_backbone(const std::vector& c) { + class get_cycle_edges : public boost::default_dfs_visitor { + public: + std::set tE; + std::list& bE; + + get_cycle_edges(std::list& bE) : bE(bE) {}; + + void tree_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { + tE.insert(g[e].index); + } + + void back_edge(boost::graph_traits::edge_descriptor e, const Graph& g) { + if (tE.find(g[e].index) == tE.end()) { + bE.push_back(g[e].index); + } + } + }; + + class find_cycle : public boost::default_dfs_visitor { + public: + const graph& G; + std::array& E; + unsigned end; + struct done{}; + + find_cycle(const graph& G, std::array& E, unsigned end) : G(G), E(E), end(end) {} + + void discover_vertex(boost::graph_traits::vertex_descriptor v, const Graph& g) { + if (v == end) { + throw done{}; + } + } + + void examine_edge(boost::graph_traits::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 finish_edge(boost::graph_traits::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]--; + } + }; + + class get_cluster_edges : public boost::default_dfs_visitor { + public: + std::set& E; + + get_cluster_edges(std::set& E) : E(E) {} + + void examine_edge(boost::graph_traits::edge_descriptor ed, const Graph& g) { + unsigned e = g[ed].index; + E.insert(e); + } + }; + + for (unsigned i = 0; i < G.edges.size(); i++) { + if ((!backbone[i]) && boost::same_component(G.dual_edges[i].v[0], G.dual_edges[i].v[1], ds)) { + 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 cedges = {}; + get_cycle_edges vis(cedges); + std::vector cm1(G.dual_vertices.size()); + boost::depth_first_visit(dG, G.dual_edges[i].v[0], 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::array cycle = {0, 0}; + if (G.dual_edges[i].crossings[0]) cycle[0]++; + if (G.dual_edges[i].crossings[1]) cycle[1]++; + find_cycle vis2(G, cycle, G.dual_edges[i].v[1]); + std::vector 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> other_cycles = {cycle}; + for (unsigned e : cedges) { + std::array 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 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 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); + } + + if (std::any_of(other_cycles.begin(), other_cycles.end(), [](std::array c){return c[0] % 2 == 0 && c[1] % 2 == 0;})) { + 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 component_edges = {}; + get_cluster_edges visc(component_edges); + std::vector 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 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]; + } + } + 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 (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); + } + } + } + } + } + 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& LE; + std::set& RE; + unsigned start_edge; + bool initialized; + bool done_left; + unsigned disc; + unsigned dr; + + get_tie_edges(std::set& LE, std::set& RE) : LE(LE), RE(RE) { + initialized = false; + done_left = false; + disc = 0; + dr = 0; + } -fuse_network::fuse_network(const graph& G, cholmod_common* c, double weight) : - network(G), ohm_x(G, 0, c), ohm_y(G, 1, c) {} + void discover_vertex(boost::graph_traits::vertex_descriptor vd, const Graph& g) { + if (!initialized && disc > 0) { + start_edge = vd; + initialized = true; + } + disc++; + } -fuse_network::fuse_network(const fuse_network& n) : - network(n), ohm_x(n.ohm_x), ohm_y(n.ohm_y) { + void finish_vertex(boost::graph_traits::vertex_descriptor vd, const Graph& g) { + if (vd == start_edge) { + done_left = true; + } + } + + void examine_edge(boost::graph_traits::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); + } + } + + + }; + + std::set bb_verts; + + for (unsigned i = 0; i < G.edges.size(); i++) { + if (!backbone[i]) { + bb_verts.insert(G.edges[i].v[0]); + bb_verts.insert(G.edges[i].v[1]); + } + } + + for (unsigned v : bb_verts) { + bool found_pair = false; + unsigned d1, d2, p1, p2; + for (unsigned i = 0; i < G.vertices[v].nd.size(); i++) { + for (unsigned j = 2; j < G.vertices[v].nd.size() - 1; j++) { + unsigned l = (i+j)%G.vertices[v].nd.size(); + if (boost::same_component(G.vertices[v].nd[i], G.vertices[v].nd[l], ds)) { + unsigned il = i < l ? i : l; + unsigned ig = i < l ? l : i; + auto it_search = seen_pairs.find({G.vertices[v].nd[il], G.vertices[v].nd[ig]}); + if (it_search == seen_pairs.end()) { + bool any_intervening1 = false; + bool any_intervening_e_1 = false; + bool any_intervening2 = false; + bool any_intervening_e_2 = false; + + for (unsigned k = il + 1; k < ig; k++) { + if (!boost::same_component(G.vertices[v].nd[i], G.vertices[v].nd[k], ds)) { + any_intervening1 = true; + break; + } + } + for (unsigned k = (ig + 1); k%G.vertices[v].nd.size() != il; k++) { + if (!boost::same_component(G.vertices[v].nd[i], G.vertices[v].nd[k%G.vertices[v].nd.size()], ds)) { + any_intervening2 = true; + break; + } + } + if (any_intervening2 && !any_intervening1) { + for (unsigned k = il + 1; k <= ig; k++) { + if (!backbone[G.vertices[v].ne[k]]) { + any_intervening_e_1 = true; + break; + } + } + } + if (any_intervening1 && !any_intervening2) { + for (unsigned k = ig + 1; k%G.vertices[v].nd.size() != il + 1; k++) { + if (!backbone[G.vertices[v].ne[k%G.vertices[v].nd.size()]]) { + any_intervening_e_2 = true; + break; + } + } + } + if ((any_intervening1 && any_intervening2) || (any_intervening1 && any_intervening_e_2) || (any_intervening2 && any_intervening_e_1)) { + found_pair = true; + p1 = il; + p2 = ig; + d1 = G.vertices[v].nd[il]; + d2 = G.vertices[v].nd[ig]; + break; + } + } + } + } + } + + if (found_pair) { + std::list cedges = {}; + get_cycle_edges vis(cedges); + std::vector 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::array cycle = {0, 0}; + for (unsigned k = p1; k < p2; k++) { + if (G.dual_edges[G.vertices[v].ne[k]].crossings[0]) cycle[0]++; + if (G.dual_edges[G.vertices[v].ne[k]].crossings[1]) cycle[1]++; + } + find_cycle vis2(G, cycle, d2); + std::vector 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> other_cycles = {cycle}; + for (unsigned e : cedges) { + std::array 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 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 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); + } + + if (std::any_of(other_cycles.begin(), other_cycles.end(), [](std::array c){return c[0] % 2 == 0 && c[1] % 2 == 0;})) { + seen_pairs[{d1, d2}] = true; + + std::set left_edges = {}; + std::set right_edges = {}; + get_tie_edges visc(left_edges, right_edges); + std::vector 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 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]; + } + } + } + + std::array total_cur_right = {0.0, 0.0}; + for (unsigned e : right_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_right[0] += c[e]; + } else { + total_cur_right[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]; + } else { + total_cur_right[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) { + 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); + } + } + } } -void fuse_network::break_edge(unsigned e, bool unbreak) { +void network::break_edge(unsigned e, bool unbreak) { fuses[e] = !unbreak; - ohm_x.break_edge(e, unbreak); - ohm_y.break_edge(e, unbreak); + backbone[e] = !unbreak; + boost::remove_edge(G.edges[e].v[0], G.edges[e].v[1], bG); + boost::add_edge(G.dual_edges[e].v[0], G.dual_edges[e].v[1], {e}, dG); + ds.union_set(G.dual_edges[e].v[0], G.dual_edges[e].v[1]); + px.break_edge(e, unbreak); + py.break_edge(e, unbreak); } + +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) {} + current_info fuse_network::get_current_info() { - current_info cx = ohm_x.solve(fuses); - current_info cy = ohm_y.solve(fuses); + px.solve(fuses); + py.solve(fuses); - bool done_x = cx.conductivity[0] < 1.0 / G.edges.size(); - bool done_y = cy.conductivity[1] < 1.0 / G.edges.size(); + 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 = {cx.conductivity[0], cy.conductivity[1]}; + 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 * cy.currents[i] / cy.conductivity[1]); + 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) * cx.currents[i] / cx.conductivity[0]); + ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0]); } } else if (!done_x && !done_y) { for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = fabs((1 - weight) * cx.currents[i] / cx.conductivity[0] + - weight * cy.currents[i] / cy.conductivity[1]); + ctot.currents[i] = fabs((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0] + + weight * py.sol.currents[i] / py.sol.conductivity[1]); } } @@ -118,42 +504,35 @@ current_info fuse_network::get_current_info() { elastic_network::elastic_network(const graph& G, cholmod_common* c, double weight) : - network(G), weight(weight), - hook_x(G, 0, c), hook_y(G, 1, c) {} + network(G, c), weight(weight) {} elastic_network::elastic_network(const elastic_network& n) : - network(n), weight(n.weight), hook_x(n.hook_x), hook_y(n.hook_y) {} - -void elastic_network::break_edge(unsigned e, bool unbreak) { - fuses[e] = !unbreak; - hook_x.break_edge(e, unbreak); - hook_y.break_edge(e, unbreak); -} + network(n), weight(n.weight) {} current_info elastic_network::get_current_info() { - current_info cx = hook_x.solve(fuses); - current_info cy = hook_y.solve(fuses); + px.solve(fuses); + py.solve(fuses); - bool done_x = cx.conductivity[0] < 1.0 / G.edges.size(); - bool done_y = cy.conductivity[1] < 1.0 / G.edges.size(); + 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[0] = cx.conductivity[0]; - ctot.conductivity[1] = cy.conductivity[1]; + ctot.conductivity[0] = px.sol.conductivity[0]; + ctot.conductivity[1] = py.sol.conductivity[1]; if (done_x && !done_y) { for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity[1]; + ctot.currents[i] = weight * fabs(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] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity[0]; + ctot.currents[i] = (1 - weight) * fabs(px.sol.currents[i]) / px.sol.conductivity[0]; } } else if (!done_x && !done_y) { for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity[0], 2) + - pow(weight * cy.currents[i] / cy.conductivity[1], 2)); + ctot.currents[i] = sqrt(pow((1 - weight) * px.sol.currents[i] / px.sol.conductivity[0], 2) + + pow(weight * py.sol.currents[i] / py.sol.conductivity[1], 2)); } } @@ -161,23 +540,21 @@ current_info elastic_network::get_current_info() { } -percolation_network::percolation_network(const graph& G, cholmod_common* c) : - network(G), px(G, 0, c), py(G, 1, c) {} +percolation_network::percolation_network(const graph& G, cholmod_common* c) : network(G, c) {} -percolation_network::percolation_network(const percolation_network& n) : - network(n), px(n.px), py(n.py) {} +percolation_network::percolation_network(const percolation_network& n) : network(n) {} current_info percolation_network::get_current_info() { current_info ctot; ctot.currents.resize(G.edges.size(), 0.0); - current_info cx = px.solve(fuses); - current_info cy = py.solve(fuses); + px.solve(fuses); + py.solve(fuses); - ctot.conductivity = {cx.conductivity[0], cy.conductivity[1]}; + ctot.conductivity = {px.sol.conductivity[0], py.sol.conductivity[1]}; for (unsigned i = 0; i < G.edges.size(); i++) { - if (fabs(cx.currents[i]) > CURRENT_CUTOFF || fabs(cy.currents[i]) > CURRENT_CUTOFF) { + if (fabs(px.sol.currents[i]) > CURRENT_CUTOFF || fabs(py.sol.currents[i]) > CURRENT_CUTOFF) { ctot.currents[i] = 1.0; } } @@ -185,9 +562,3 @@ current_info percolation_network::get_current_info() { return ctot; } -void percolation_network::break_edge(unsigned e, bool unbreak) { - fuses[e] = !unbreak; - px.break_edge(e, unbreak); - py.break_edge(e, unbreak); -} - diff --git a/lib/src/problem.cpp b/lib/src/problem.cpp index f66a35c..0645244 100644 --- a/lib/src/problem.cpp +++ b/lib/src/problem.cpp @@ -75,6 +75,8 @@ problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_c factor = CHOL_F(analyze)(laplacian, c); CHOL_F(factorize)(laplacian, factor, c); CHOL_F(free_sparse)(&laplacian, c); + + sol.currents.resize(G.edges.size()); } problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) { @@ -97,7 +99,7 @@ problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, CHOL_F(free_triplet)(&t, c); } -problem::problem(const problem& other) : G(other.G), axis(other.axis), c(other.c) { +problem::problem(const problem& other) : G(other.G), axis(other.axis), c(other.c), sol(other.sol) { b = CHOL_F(copy_dense)(other.b, c); factor = CHOL_F(copy_factor)(other.factor, c); voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c); @@ -109,7 +111,7 @@ problem::~problem() { CHOL_F(free_sparse)(&voltcurmat, c); } -current_info problem::solve(std::vector& fuses) { +void problem::solve(std::vector& fuses) { cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); if (((double *)x->x)[0] != ((double *)x->x)[0]) { @@ -122,15 +124,13 @@ current_info problem::solve(std::vector& fuses) { double beta[2] = {0, 0}; CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); - std::vector currents(G.edges.size()); - - std::array total_current = {0, 0}; + sol.conductivity = {0, 0}; for (int i = 0; i < G.edges.size(); i++) { if (fuses[i]) { - currents[i] = 0; + sol.currents[i] = 0; } else { - currents[i] = ((double *)y->x)[i]; + sol.currents[i] = ((double *)y->x)[i]; graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r; graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r; @@ -144,22 +144,20 @@ current_info problem::solve(std::vector& fuses) { } if (comp) { - currents[i] += 1.0; - total_current[axis] += currents[i]; + sol.currents[i] += 1.0; + sol.conductivity[axis] += sol.currents[i]; } else { - currents[i] -= 1.0; - total_current[axis] -= currents[i]; + sol.currents[i] -= 1.0; + sol.conductivity[axis] -= sol.currents[i]; } } } } - total_current[!axis] = 0.0; + sol.conductivity[!axis] = 0.0; CHOL_F(free_dense)(&x, c); CHOL_F(free_dense)(&y, c); - - return {total_current, currents}; } void problem::break_edge(unsigned e, bool unbreak) { diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp index 659cacf..c3bffa2 100644 --- a/src/analysis_tools.hpp +++ b/src/analysis_tools.hpp @@ -1,10 +1,4 @@ -#include -#include -#include -#include -#include - #include #include #include @@ -13,14 +7,6 @@ #include -struct EdgeProperties { - unsigned index; -}; - -typedef boost::adjacency_list Graph; -typedef boost::graph_traits::vertex_descriptor Vertex; -typedef boost::graph_traits::vertices_size_type VertexIndex; - template bool is_shorter(const std::list &, const std::list &); diff --git a/src/animate.cpp b/src/animate.cpp index 9ba335b..532bd70 100644 --- a/src/animate.cpp +++ b/src/animate.cpp @@ -19,9 +19,10 @@ void animate::pre_fracture(const network &n) { boost::remove_edge_if(trivial, G); glClearColor(1.0f, 1.0f, 1.0f, 1.0f ); + glLineWidth(5); glClear(GL_COLOR_BUFFER_BIT); glBegin(GL_LINES); - glColor3f(0.0f, 0.0f, 0.0f); + glColor3f(0.9f, 0.9f, 0.9f); for (unsigned 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; @@ -59,19 +60,62 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned 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); glBegin(GL_LINES); // Each set of 3 vertices form a triangle - glColor3f(1.0f, 1.0f, 1.0f); // Blue - 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 (n.G.edges[i].crossings[0]) { + /* + glColor3f(0.0f, 0.0f, 0.0f); // Blue + graph::coordinate r1 = n.G.dual_vertices[n.G.dual_edges[i].v[0]].r; + graph::coordinate r2 = n.G.dual_vertices[n.G.dual_edges[i].v[1]].r; + + if (n.G.dual_edges[i].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[i].crossings[1]) { + if (r1.y < r2.y) { + r1.y += n.G.L.y; + } else { + r2.y += n.G.L.y; + } + } + + glVertex2d(r1.x, r1.y); + glVertex2d(r2.x, r2.y); + */ + + bool weird = false; + + for (unsigned j = 0; j < n.G.edges.size(); j++) { + bool draw = false; + if (cur.currents[j] < 1e-9 && !n.backbone[j]) { + glColor3f(1.0f, 0.0f, 0.0f); // Blue + weird = true; + draw = true; + } else if (n.backbone[j] && !n.fuses[j] && j != i) { + glColor3f(0.8f, 0.8f, 0.8f); // Blue + draw = true; + } else if (!n.fuses[j]) { + glColor3f(0.0f, 0.0f, 0.0f); // Blue + draw = 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; + + if (n.G.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[j].crossings[1]) { if (r1.y < r2.y) { r1.y += n.G.L.y; } else { @@ -81,12 +125,17 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i) glVertex2d(r1.x, r1.y); glVertex2d(r2.x, r2.y); + } + + } glEnd(); glFlush(); + if (weird) {std::cout << "\n"; getchar();} } void animate::post_fracture(network &n) { + /* std::list crack; // unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]]; unsigned crack_component = 10000; @@ -223,5 +272,6 @@ void animate::post_fracture(network &n) { glFlush(); } } +*/ } diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp index c3b4a69..9d875e6 100644 --- a/src/animate_fracture.cpp +++ b/src/animate_fracture.cpp @@ -31,20 +31,20 @@ int main(int argc, char* argv[]) { int opt; unsigned N = 1; - double Lx = 16.0; - double Ly = 16.0; + unsigned n = 16; + double a = 16.0; double beta = 0.5; - while ((opt = getopt(argc, argv, "X:Y:N:b:")) != -1) { + while ((opt = getopt(argc, argv, "n:a:N:b:")) != -1) { switch (opt) { case 'N': N = (unsigned)atof(optarg); break; - case 'X': - Lx = atof(optarg); + case 'n': + n = atoi(optarg); break; - case 'Y': - Ly = atof(optarg); + case 'a': + a = atof(optarg); break; case 'b': beta = atof(optarg); @@ -57,16 +57,16 @@ int main(int argc, char* argv[]) { cholmod_common c; CHOL_F(start)(&c); - animate meas(Lx, Ly, 700, argc, argv); + animate meas(sqrt(2*n *a), sqrt( 2*n / a), 700, argc, argv); randutils::auto_seed_128 seeds; std::mt19937 rng{seeds}; for (unsigned trial = 0; trial < N; trial++) { - graph G(Lx, Ly, rng); + graph G(n, a, rng); elastic_network elastic_network(G, &c); elastic_network.set_thresholds(beta, rng); - elastic_network.fracture(meas); + elastic_network.fracture(meas, true); if (quit.load()) break; diff --git a/src/animate_fracture_square.cpp b/src/animate_fracture_square.cpp index 84cc8a5..bc4c3a2 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); + tmp_network.fracture(meas, false); if (quit.load()) break; diff --git a/src/fracture.cpp b/src/fracture.cpp index 7df3f67..483a3d2 100644 --- a/src/fracture.cpp +++ b/src/fracture.cpp @@ -38,8 +38,9 @@ int main(int argc, char* argv[]) { unsigned n = 128; double a = 1.0; bool use_aN = false; + double w = 0.5; - while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:")) != -1) { + while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:w:")) != -1) { switch (opt) { case 'N': N = (unsigned)atof(optarg); @@ -61,6 +62,9 @@ int main(int argc, char* argv[]) { a = atof(optarg); use_aN = true; break; + case 'w': + w = atof(optarg); + break; default: exit(1); } @@ -82,9 +86,9 @@ int main(int argc, char* argv[]) { while (true) { try { graph G(n, a, rng); - elastic_network fuse_network(G, &c); + elastic_network fuse_network(G, &c, w); fuse_network.set_thresholds(beta, rng); - fuse_network.fracture(meas); + fuse_network.fracture(meas, false); break; } catch (std::exception &e) { std::cout << e.what() << '\n'; @@ -105,9 +109,9 @@ int main(int argc, char* argv[]) { while (true) { try { graph G(Lx, Ly); - elastic_network fuse_network(G, &c); + elastic_network fuse_network(G, &c, w); fuse_network.set_thresholds(beta, rng); - fuse_network.fracture(meas); + fuse_network.fracture(meas, false); break; } catch (std::exception &e) { std::cout << e.what() << '\n'; diff --git a/src/measurements.cpp b/src/measurements.cpp index 0732d24..c58fd84 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -153,7 +153,8 @@ ma::ma(unsigned n, double a, double beta) : sc(2 * n), sa(3 * n), sA(3 * n), - si(10000) + si(10000), + sI(10000) { if (beta != 0.0) { model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_"; @@ -173,7 +174,8 @@ ma::ma(unsigned Lx, unsigned Ly, double beta) : sc(Lx * Ly / 2), sa(Lx * Ly), sA(Lx * Ly), - si(10000) + si(10000), + sI(10000) { if (beta != 0.0) { model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_"; @@ -193,6 +195,7 @@ ma::~ma() { update_distribution_file("sa", sa, model_string); update_distribution_file("sA", sA, model_string); update_distribution_file("si", si, model_string); + update_distribution_file("sI", sI, model_string); } void ma::pre_fracture(const network&) { @@ -212,8 +215,10 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) { } for (unsigned j = 0; j < cur.currents.size(); j++) { - if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0) { + if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0 && (!net.backbone[j] || j == i)) { si[(unsigned)(10000 * (logl(cur.currents[j]) + 100) / 100)]++; + } else if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0 && net.backbone[j] && j != i) { + sI[(unsigned)(10000 * (logl(cur.currents[j]) + 100) / 100)]++; } } @@ -226,7 +231,13 @@ void ma::post_fracture(network &n) { std::vector component(boost::num_vertices(G)); unsigned num = boost::connected_components(G, &component[0]); if (post_cracks.size() > 2 || post_cracks.size() == 0) { - std::cout << post_cracks.size() << "\n"; + for (auto c : post_cracks) { + for (unsigned e : c.second) { + std::cout << e << " "; + } + std::cout << "\n"; + } + getchar(); throw badcycleex; } for (auto c : post_cracks) { @@ -278,6 +289,7 @@ void ma::post_fracture(network &n) { sc[new_components[i].size() - 1]++; } + /* current_info ct = n.get_current_info(); @@ -297,6 +309,7 @@ void ma::post_fracture(network &n) { } sb[bb_size - 1]++; + */ av_it++; diff --git a/src/measurements.hpp b/src/measurements.hpp index 5b76e26..274a550 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -35,6 +35,7 @@ class ma : public hooks { std::vector sA; // non-final avalanche size distribution std::vector si; + std::vector sI; public: long double lv; -- cgit v1.2.3-54-g00ecf