From ae6ad8569615d81fd4b4a8f13318c8f90a768a37 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 27 May 2019 22:26:21 -0400 Subject: refactored much of the fracture library and fixed some percolation measurements --- lib/CMakeLists.txt | 2 +- lib/include/current_info.hpp | 2 +- lib/include/network.hpp | 62 ++----- lib/include/problem.hpp | 38 ++++ lib/src/network.cpp | 417 +++++++++---------------------------------- lib/src/problem.cpp | 217 ++++++++++++++++++++++ src/analysis_tools.hpp | 2 + src/perc_meas.cpp | 4 +- src/percolation.cpp | 5 +- 9 files changed, 368 insertions(+), 381 deletions(-) create mode 100644 lib/include/problem.hpp create mode 100644 lib/src/problem.cpp diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 49a9602..2592f3d 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -1,4 +1,4 @@ -add_library(frac SHARED src/graph.cpp src/network.cpp) +add_library(frac SHARED src/graph.cpp src/network.cpp src/problem.cpp) target_include_directories(frac PUBLIC include PRIVATE voronoi/src) diff --git a/lib/include/current_info.hpp b/lib/include/current_info.hpp index 90dd3c7..a84fb5b 100644 --- a/lib/include/current_info.hpp +++ b/lib/include/current_info.hpp @@ -4,7 +4,7 @@ #include typedef struct current_info { - double conductivity; + std::array conductivity; std::vector currents; } current_info; diff --git a/lib/include/network.hpp b/lib/include/network.hpp index 15cb8ea..1ffe742 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -6,44 +6,14 @@ #include #include #include +#include -#include - +#include "problem.hpp" #include "graph.hpp" #include "hooks.hpp" #include "current_info.hpp" #include "array_hash.hpp" -#ifdef FRACTURE_LONGINT - -#define CHOL_F(x) cholmod_l_##x -#define CHOL_INT long int - -#else - -#define CHOL_F(x) cholmod_##x -#define CHOL_INT int - -#endif - -class problem { - private: - const graph& G; - cholmod_dense* b; - cholmod_factor* factor; - cholmod_sparse* voltcurmat; - cholmod_common* c; - unsigned axis; - - public: - problem(const graph&, unsigned axis, cholmod_common*); - problem(const graph&, unsigned axis, cholmod_sparse*, cholmod_common*); - problem(const problem&); - ~problem(); - current_info solve(std::vector& fuses); - void break_edge(unsigned, bool unbreak = false); -}; - class network { public: const graph& G; @@ -54,45 +24,49 @@ class network { network(const network&); void set_thresholds(double, std::mt19937&); + void fracture(hooks&, bool one_axis = false); + virtual void break_edge(unsigned e, bool unbreak = false) {fuses[e] = !unbreak;}; virtual current_info get_current_info() {current_info empty; return empty;}; }; class fuse_network : public network { - public: - problem ohm; + private: + double weight; + problem ohm_x; + problem ohm_y; - fuse_network(const graph&, cholmod_common*); + public: + fuse_network(const graph&, cholmod_common*, double weight = 1.0); + fuse_network(const fuse_network&); - void fracture(hooks&, double cutoff = 1e-13); + void break_edge(unsigned, bool unbreak = false) override; + current_info get_current_info() override; }; class elastic_network : public network { private: double weight; - public: problem hook_x; problem hook_y; - elastic_network(const graph&, cholmod_common*); + public: + elastic_network(const graph&, cholmod_common*, double weight = 0.5); elastic_network(const elastic_network&); - void fracture(hooks&, double weight = 0.5, double cutoff = 1e-10); void break_edge(unsigned, bool unbreak = false) override; current_info get_current_info() override; }; class percolation_network : public network { private: - double weight; - public: - problem hook_x; - problem hook_y; + problem px; + problem py; + public: percolation_network(const graph&, cholmod_common*); percolation_network(const percolation_network&); - void fracture(hooks&, double weight = 0.5, double cutoff = 1e-10); 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 new file mode 100644 index 0000000..c2c9e09 --- /dev/null +++ b/lib/include/problem.hpp @@ -0,0 +1,38 @@ + +#include +#include + +#include +#include "graph.hpp" +#include "current_info.hpp" + +#ifdef FRACTURE_LONGINT + +#define CHOL_F(x) cholmod_l_##x +#define CHOL_INT long int + +#else + +#define CHOL_F(x) cholmod_##x +#define CHOL_INT int + +#endif + +class problem { + private: + const graph& G; + unsigned axis; + cholmod_dense* b; + cholmod_factor* factor; + cholmod_sparse* voltcurmat; + cholmod_common* c; + + public: + 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 break_edge(unsigned, bool unbreak = false); +}; + diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 1d9b3f9..b081a3c 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -1,13 +1,6 @@ #include "network.hpp" - -class nanException: public std::exception -{ - virtual const char* what() const throw() - { - return "The linear problem returned NaN."; - } -} nanex; +#include class nofuseException: public std::exception { @@ -17,210 +10,6 @@ class nofuseException: public std::exception } } nofuseex; -problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), voltcurmat(vcmat), c(c), axis(axis) { - b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c); - for (unsigned i = 0; i < G.edges.size(); i++) { - graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r; - graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r; - - if (G.edges[i].crossings[axis]) { - bool ind; - if (axis == 1) { - ind = v0.y < v1.y ? 0 : 1; - } else { - ind = v0.x < v1.x ? 0 : 1; - } - - ((double *)b->x)[G.edges[i].v[ind]] += 1.0; - ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0; - } - } - - unsigned nnz = G.vertices.size() + G.edges.size(); - - cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz, 1, CHOLMOD_REAL, c); - - for (unsigned i = 0; i < G.vertices.size(); i++) { - ((CHOL_INT *)t->i)[i] = i; - ((CHOL_INT *)t->j)[i] = i; - ((double *)t->x)[i] = 0.0; - } - - unsigned terms = G.vertices.size(); - - std::unordered_map, unsigned> known_edges; - - for (unsigned i = 0; i < G.edges.size(); i++) { - unsigned v0 = G.edges[i].v[0]; - unsigned v1 = G.edges[i].v[1]; - - ((double *)t->x)[v0]++; - ((double *)t->x)[v1]++; - - unsigned s0 = v0 < v1 ? v0 : v1; - unsigned s1 = v0 < v1 ? v1 : v0; - - auto it = known_edges.find({s0, s1}); - - if (it == known_edges.end()) { - ((CHOL_INT *)t->i)[terms] = s1; - ((CHOL_INT *)t->j)[terms] = s0; - ((double *)t->x)[terms] = -1.0; - - known_edges[{s0, s1}] = terms; - terms++; - } else { - ((double *)t->x)[it->second] -= 1.0; - } - } - - ((double *)t->x)[0]++; - - t->nnz = terms; - - cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, terms, c); - CHOL_F(free_triplet)(&t, c); - factor = CHOL_F(analyze)(laplacian, c); - CHOL_F(factorize)(laplacian, factor, c); - CHOL_F(free_sparse)(&laplacian, c); -} - -problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) { - cholmod_triplet *t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c); - - t->nnz = 2 * G.edges.size(); - - for (unsigned i = 0; i < G.edges.size(); i++) { - ((CHOL_INT *)t->i)[2 * i] = i; - ((CHOL_INT *)t->j)[2 * i] = G.edges[i].v[0]; - ((double *)t->x)[2 * i] = 1.0; - - ((CHOL_INT *)t->i)[2 * i + 1] = i; - ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1]; - ((double *)t->x)[2 * i + 1] = -1.0; - } - - voltcurmat = CHOL_F(triplet_to_sparse)(t, 2 * G.edges.size(), c); - - CHOL_F(free_triplet)(&t, c); -} - -problem::problem(const problem& other) : G(other.G), c(other.c), axis(other.axis) { - b = CHOL_F(copy_dense)(other.b, c); - factor = CHOL_F(copy_factor)(other.factor, c); - voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c); -} - -problem::~problem() { - CHOL_F(free_dense)(&b, c); - CHOL_F(free_factor)(&factor, c); - CHOL_F(free_sparse)(&voltcurmat, c); -} - -current_info 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]) { - throw nanex; - } - - cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c); - - double alpha[2] = {1, 0}; - double beta[2] = {0, 0}; - CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); - - std::vector currents(G.edges.size()); - - double total_current = 0; - - for (int i = 0; i < G.edges.size(); i++) { - if (fuses[i]) { - currents[i] = 0; - } else { - 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; - - if (G.edges[i].crossings[axis]) { - bool comp; - if (axis == 1) { - comp = v0.y > v1.y; - } else { - comp = v0.x > v1.x; - } - - if (comp) { - currents[i] += 1.0; - total_current += currents[i]; - } else { - currents[i] -= 1.0; - total_current -= currents[i]; - } - } - } - } - - CHOL_F(free_dense)(&x, c); - CHOL_F(free_dense)(&y, c); - - return {total_current, currents}; -} - -void problem::break_edge(unsigned e, bool unbreak) { - unsigned v0 = G.edges[e].v[0]; - unsigned v1 = G.edges[e].v[1]; - - unsigned n = factor->n; - - cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); - - unsigned s1, s2; - s1 = v0 < v1 ? v0 : v1; - s2 = v0 < v1 ? v1 : v0; - - CHOL_INT *pp = (CHOL_INT *)update_mat->p; - CHOL_INT *ii = (CHOL_INT *)update_mat->i; - double *xx = (double *)update_mat->x; - - for (unsigned i = 0; i <= s1; i++) { - pp[i] = 0; - } - - for (unsigned i = s1 + 1; i <= n; i++) { - pp[i] = 2; - } - - ii[0] = s1; - ii[1] = s2; - xx[0] = 1; - xx[1] = -1; - - cholmod_sparse *perm_update_mat = CHOL_F(submatrix)( - update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c); - - CHOL_F(updown)(unbreak, perm_update_mat, factor, c); - - CHOL_F(free_sparse)(&perm_update_mat, c); - CHOL_F(free_sparse)(&update_mat, c); - - graph::coordinate r0 = G.vertices[v0].r; - graph::coordinate r1 = G.vertices[v1].r; - - if (G.edges[e].crossings[axis]) { - bool ind; - if (axis == 1) { - ind = r0.y < r1.y ? unbreak : !unbreak; - } else { - ind = r0.x < r1.x ? unbreak : !unbreak; - } - - ((double *)b->x)[G.edges[e].v[ind]] -= 1.0; - ((double *)b->x)[G.edges[e].v[!ind]] += 1.0; - } -} - network::network(const graph& G) : G(G), fuses(G.edges.size()), thresholds(G.edges.size()) {} network::network(const network& n) : G(n.G), fuses(n.fuses), thresholds(n.thresholds) {} @@ -244,16 +33,19 @@ void network::set_thresholds(double beta, std::mt19937& rng) { } } -fuse_network::fuse_network(const graph& G, cholmod_common *c) : network(G), ohm(G, 1, c) { -} - -void fuse_network::fracture(hooks& m, double cutoff) { +void network::fracture(hooks& m, bool one_axis) { m.pre_fracture(*this); while (true) { - current_info ci = ohm.solve(fuses); + current_info c = this->get_current_info(); + + double min_cond = 1.0 / G.edges.size(); - if (ci.conductivity < cutoff * G.vertices.size()) { + 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; } @@ -261,8 +53,9 @@ void fuse_network::fracture(hooks& m, double cutoff) { long double max_val = std::numeric_limits::lowest(); for (unsigned i = 0; i < G.edges.size(); i++) { - if (!fuses[i] && fabs(ci.currents[i]) > cutoff) { - long double val = logl(fabs(ci.currents[i])) - thresholds[i]; + if (!fuses[i] && c.currents[i] > min_cond) { + long double val = logl(c.currents[i]) - thresholds[i]; + if (val > max_val) { max_val = val; max_pos = i; @@ -274,172 +67,130 @@ void fuse_network::fracture(hooks& m, double cutoff) { throw nofuseex; } - fuses[max_pos] = true; - ohm.break_edge(max_pos); - - m.bond_broken(*this, ci, max_pos); + this->break_edge(max_pos); + m.bond_broken(*this, c, max_pos); } m.post_fracture(*this); } -elastic_network::elastic_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) { -} - -elastic_network::elastic_network(const elastic_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) { -} - -void elastic_network::fracture(hooks& m, double weight, double cutoff) { - this->weight = weight; - m.pre_fracture(*this); - - while (true) { - current_info ctot = this->get_current_info(); - - if (ctot.conductivity == 0) { - break; - } - - unsigned max_pos = UINT_MAX; - long double max_val = std::numeric_limits::lowest(); - - for (unsigned i = 0; i < G.edges.size(); i++) { - if (!fuses[i] && ctot.currents[i] > 1.0 / G.edges.size()) { - long double val = logl(ctot.currents[i]) - thresholds[i]; - if (val > max_val) { - max_val = val; - max_pos = i; - } - } - } - - if (max_pos == UINT_MAX) { - throw nofuseex; - } - this->break_edge(max_pos); +fuse_network::fuse_network(const graph& G, cholmod_common* c, double weight) : + network(G), ohm_x(G, 0, c), ohm_y(G, 1, c) {} - m.bond_broken(*this, ctot, max_pos); - } +fuse_network::fuse_network(const fuse_network& n) : + network(n), ohm_x(n.ohm_x), ohm_y(n.ohm_y) { +} - m.post_fracture(*this); +void fuse_network::break_edge(unsigned e, bool unbreak) { + fuses[e] = !unbreak; + ohm_x.break_edge(e, unbreak); + ohm_y.break_edge(e, unbreak); } -current_info elastic_network::get_current_info() { - current_info cx = hook_x.solve(fuses); - current_info cy = hook_y.solve(fuses); +current_info fuse_network::get_current_info() { + current_info cx = ohm_x.solve(fuses); + current_info cy = ohm_y.solve(fuses); - bool done_x = cx.conductivity < 1.0 / G.edges.size(); - bool done_y = cy.conductivity < 1.0 / G.edges.size(); + bool done_x = cx.conductivity[0] < 1.0 / G.edges.size(); + bool done_y = cy.conductivity[1] < 1.0 / G.edges.size(); current_info ctot; ctot.currents.resize(G.edges.size()); + ctot.conductivity = {cx.conductivity[0], cy.conductivity[1]}; - if (done_x && done_y) { - ctot.conductivity = 0; - } else if (done_x) { - ctot.conductivity = cy.conductivity; + if (done_x && !done_y) { for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity; + ctot.currents[i] = fabs(weight * cy.currents[i] / cy.conductivity[1]); } - } else if (done_y) { - ctot.conductivity = cx.conductivity; + } 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; + ctot.currents[i] = fabs((1 - weight) * cx.currents[i] / cx.conductivity[0]); } - } else { - ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); + } 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, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2)); + ctot.currents[i] = fabs((1 - weight) * cx.currents[i] / cx.conductivity[0] + + weight * cy.currents[i] / cy.conductivity[1]); } } return ctot; } + +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) {} + +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); } -percolation_network::percolation_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) { -} +current_info elastic_network::get_current_info() { + current_info cx = hook_x.solve(fuses); + current_info cy = hook_y.solve(fuses); -percolation_network::percolation_network(const percolation_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) { -} + bool done_x = cx.conductivity[0] < 1.0 / G.edges.size(); + bool done_y = cy.conductivity[1] < 1.0 / G.edges.size(); -void percolation_network::fracture(hooks& m, double weight, double cutoff) { - this->weight = weight; - m.pre_fracture(*this); - - while (true) { - current_info ctot = this->get_current_info(); + current_info ctot; + ctot.currents.resize(G.edges.size()); + ctot.conductivity[0] = cx.conductivity[0]; + ctot.conductivity[1] = cy.conductivity[1]; - if (ctot.conductivity == 0) { - break; + 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]; } - - unsigned max_pos = UINT_MAX; - long double max_val = std::numeric_limits::lowest(); - + } else if (done_y && !done_x) { for (unsigned i = 0; i < G.edges.size(); i++) { - if (!fuses[i] && ctot.currents[i] > 1.0 / G.edges.size()) { - long double val = - thresholds[i]; - if (val > max_val) { - max_val = val; - max_pos = i; - } - } + ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity[0]; } - - if (max_pos == UINT_MAX) { - throw nofuseex; + } 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)); } - - this->break_edge(max_pos); - - m.bond_broken(*this, ctot, max_pos); } - m.post_fracture(*this); + return ctot; } -current_info percolation_network::get_current_info() { - current_info cx = hook_x.solve(fuses); - current_info cy = hook_y.solve(fuses); - bool done_x = cx.conductivity < 1.0 / G.edges.size(); - bool done_y = cy.conductivity < 1.0 / G.edges.size(); +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 percolation_network& n) : + network(n), px(n.px), py(n.py) {} +current_info percolation_network::get_current_info() { current_info ctot; - ctot.currents.resize(G.edges.size()); + ctot.currents.resize(G.edges.size(), 0); - if (done_x && done_y) { - ctot.conductivity = 0; - } else if (done_x) { - ctot.conductivity = cy.conductivity; - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity; - } - } else if (done_y) { - ctot.conductivity = cx.conductivity; - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity; - } - } else { - ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2)); + current_info cx = px.solve(fuses); + current_info cy = py.solve(fuses); + + ctot.conductivity = {cx.conductivity[0], cy.conductivity[1]}; + + double min_cur = 1.0 / G.edges.size(); + + for (unsigned i = 0; i < G.edges.size(); i++) { + if (fabs(cx.currents[i]) > min_cur || fabs(cy.currents[i]) > min_cur) { + ctot.currents[i] = 1.0; } } - + return ctot; } void percolation_network::break_edge(unsigned e, bool unbreak) { fuses[e] = !unbreak; - hook_x.break_edge(e, unbreak); - hook_y.break_edge(e, unbreak); + px.break_edge(e, unbreak); + py.break_edge(e, unbreak); } diff --git a/lib/src/problem.cpp b/lib/src/problem.cpp new file mode 100644 index 0000000..f66a35c --- /dev/null +++ b/lib/src/problem.cpp @@ -0,0 +1,217 @@ + +#include "problem.hpp" + +class nanException: public std::exception +{ + virtual const char* what() const throw() + { + return "The linear problem returned NaN."; + } +} nanex; + +problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), axis(axis), voltcurmat(vcmat), c(c) { + b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c); + for (unsigned i = 0; i < G.edges.size(); i++) { + graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r; + graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r; + + if (G.edges[i].crossings[axis]) { + bool ind; + if (axis == 1) { + ind = v0.y < v1.y ? 0 : 1; + } else { + ind = v0.x < v1.x ? 0 : 1; + } + + ((double *)b->x)[G.edges[i].v[ind]] += 1.0; + ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0; + } + } + + unsigned nnz = G.vertices.size() + G.edges.size(); + + cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz, 1, CHOLMOD_REAL, c); + + for (unsigned i = 0; i < G.vertices.size(); i++) { + ((CHOL_INT *)t->i)[i] = i; + ((CHOL_INT *)t->j)[i] = i; + ((double *)t->x)[i] = 0.0; + } + + unsigned terms = G.vertices.size(); + + std::unordered_map, unsigned> known_edges; + + for (unsigned i = 0; i < G.edges.size(); i++) { + unsigned v0 = G.edges[i].v[0]; + unsigned v1 = G.edges[i].v[1]; + + ((double *)t->x)[v0]++; + ((double *)t->x)[v1]++; + + unsigned s0 = v0 < v1 ? v0 : v1; + unsigned s1 = v0 < v1 ? v1 : v0; + + auto it = known_edges.find({s0, s1}); + + if (it == known_edges.end()) { + ((CHOL_INT *)t->i)[terms] = s1; + ((CHOL_INT *)t->j)[terms] = s0; + ((double *)t->x)[terms] = -1.0; + + known_edges[{s0, s1}] = terms; + terms++; + } else { + ((double *)t->x)[it->second] -= 1.0; + } + } + + ((double *)t->x)[0]++; + + t->nnz = terms; + + cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, terms, c); + CHOL_F(free_triplet)(&t, c); + factor = CHOL_F(analyze)(laplacian, c); + CHOL_F(factorize)(laplacian, factor, c); + CHOL_F(free_sparse)(&laplacian, c); +} + +problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) { + cholmod_triplet *t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c); + + t->nnz = 2 * G.edges.size(); + + for (unsigned i = 0; i < G.edges.size(); i++) { + ((CHOL_INT *)t->i)[2 * i] = i; + ((CHOL_INT *)t->j)[2 * i] = G.edges[i].v[0]; + ((double *)t->x)[2 * i] = 1.0; + + ((CHOL_INT *)t->i)[2 * i + 1] = i; + ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1]; + ((double *)t->x)[2 * i + 1] = -1.0; + } + + voltcurmat = CHOL_F(triplet_to_sparse)(t, 2 * G.edges.size(), c); + + CHOL_F(free_triplet)(&t, c); +} + +problem::problem(const problem& other) : G(other.G), axis(other.axis), c(other.c) { + b = CHOL_F(copy_dense)(other.b, c); + factor = CHOL_F(copy_factor)(other.factor, c); + voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c); +} + +problem::~problem() { + CHOL_F(free_dense)(&b, c); + CHOL_F(free_factor)(&factor, c); + CHOL_F(free_sparse)(&voltcurmat, c); +} + +current_info 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]) { + throw nanex; + } + + cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c); + + double alpha[2] = {1, 0}; + 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}; + + for (int i = 0; i < G.edges.size(); i++) { + if (fuses[i]) { + currents[i] = 0; + } else { + 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; + + if (G.edges[i].crossings[axis]) { + bool comp; + if (axis == 1) { + comp = v0.y > v1.y; + } else { + comp = v0.x > v1.x; + } + + if (comp) { + currents[i] += 1.0; + total_current[axis] += currents[i]; + } else { + currents[i] -= 1.0; + total_current[axis] -= currents[i]; + } + } + } + } + + total_current[!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) { + unsigned v0 = G.edges[e].v[0]; + unsigned v1 = G.edges[e].v[1]; + + unsigned n = factor->n; + + cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); + + unsigned s1, s2; + s1 = v0 < v1 ? v0 : v1; + s2 = v0 < v1 ? v1 : v0; + + CHOL_INT *pp = (CHOL_INT *)update_mat->p; + CHOL_INT *ii = (CHOL_INT *)update_mat->i; + double *xx = (double *)update_mat->x; + + for (unsigned i = 0; i <= s1; i++) { + pp[i] = 0; + } + + for (unsigned i = s1 + 1; i <= n; i++) { + pp[i] = 2; + } + + ii[0] = s1; + ii[1] = s2; + xx[0] = 1; + xx[1] = -1; + + cholmod_sparse *perm_update_mat = CHOL_F(submatrix)( + update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c); + + CHOL_F(updown)(unbreak, perm_update_mat, factor, c); + + CHOL_F(free_sparse)(&perm_update_mat, c); + CHOL_F(free_sparse)(&update_mat, c); + + graph::coordinate r0 = G.vertices[v0].r; + graph::coordinate r1 = G.vertices[v1].r; + + if (G.edges[e].crossings[axis]) { + bool ind; + if (axis == 1) { + ind = r0.y < r1.y ? unbreak : !unbreak; + } else { + ind = r0.x < r1.x ? unbreak : !unbreak; + } + + ((double *)b->x)[G.edges[e].v[ind]] -= 1.0; + ((double *)b->x)[G.edges[e].v[!ind]] += 1.0; + } +} + diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp index 3beb1a8..659cacf 100644 --- a/src/analysis_tools.hpp +++ b/src/analysis_tools.hpp @@ -18,6 +18,8 @@ struct EdgeProperties { }; typedef boost::adjacency_list Graph; +typedef boost::graph_traits::vertex_descriptor Vertex; +typedef boost::graph_traits::vertices_size_type VertexIndex; template diff --git a/src/perc_meas.cpp b/src/perc_meas.cpp index 4220a93..ec7746e 100644 --- a/src/perc_meas.cpp +++ b/src/perc_meas.cpp @@ -159,7 +159,9 @@ void pm::post_fracture(network &n) { } else { ss[components[i].size() - 1]++; } - sn[p][components[i].size() - 1]++; + for (unsigned j = p; j < sn.size(); j++) { + sn[j][components[i].size() - 1]++; + } } sp[p - 1]++; diff --git a/src/percolation.cpp b/src/percolation.cpp index 03da117..757ab44 100644 --- a/src/percolation.cpp +++ b/src/percolation.cpp @@ -76,20 +76,23 @@ int main(int argc, char* argv[]) { pm meas(n, a); for (unsigned trial = 0; trial < N; trial++) { + /* while (true) { try { + */ graph G(n, a, rng); percolation_network fuse_network(G, &c); fuse_network.set_thresholds(beta, rng); fuse_network.fracture(meas); + /* break; } catch (std::exception &e) { std::cout << e.what() << '\n'; } } - if (quit.load()) break; + */ } } -- cgit v1.2.3-70-g09d2