From cb1b2e6822bdd1d1644ff2dad2d6157858e105b0 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 24 Apr 2019 23:31:40 -0400 Subject: many changes to introduce two-component, elastic-like fracture --- lib/src/network.cpp | 204 ++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 141 insertions(+), 63 deletions(-) (limited to 'lib/src') diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 1edc973..8b0e8df 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -17,14 +17,19 @@ class nofuseException: public std::exception } } nofuseex; -network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) { +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++) { - double v0y = G.vertices[G.edges[i].v[0]].r.y; - double v1y = G.vertices[G.edges[i].v[1]].r.y; - - if (G.edges[i].crossings[1]) { - bool ind = v0y < v1y ? 0 : 1; + 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; @@ -78,8 +83,10 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges. factor = CHOL_F(analyze)(laplacian, c); CHOL_F(factorize)(laplacian, factor, c); CHOL_F(free_sparse)(&laplacian, c); +} - t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, 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(); @@ -98,39 +105,70 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges. CHOL_F(free_triplet)(&t, c); } -network::network(const network &other) : c(other.c), G(other.G), fuses(other.fuses), thresholds(other.thresholds) { +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); } -network::~network() { - CHOL_F(free_sparse)(&voltcurmat, c); +problem::~problem() { CHOL_F(free_dense)(&b, c); CHOL_F(free_factor)(&factor, c); + CHOL_F(free_sparse)(&voltcurmat, c); } -void network::set_thresholds(double beta, std::mt19937& rng) { - if (beta == 0.0) { - /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */ - for (long double& threshold : thresholds) { - threshold = 1.0; - } - } else { - std::uniform_real_distribution d(0.0, 1.0); +current_info problem::solve(std::vector& fuses) { + cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); - for (long double& threshold : thresholds) { - threshold = std::numeric_limits::lowest(); + if (((double *)x->x)[0] != ((double *)x->x)[0]) { + throw nanex; + } - while (threshold == std::numeric_limits::lowest()) { - threshold = logl(d(rng)) / (long double)beta; + 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 network::break_edge(unsigned e, bool unbreak) { - fuses[e] = !unbreak; +void problem::break_edge(unsigned e, bool unbreak) { unsigned v0 = G.edges[e].v[0]; unsigned v1 = G.edges[e].v[1]; @@ -167,78 +205,116 @@ void network::break_edge(unsigned e, bool unbreak) { CHOL_F(free_sparse)(&perm_update_mat, c); CHOL_F(free_sparse)(&update_mat, c); - double v0y = G.vertices[v0].r.y; - double v1y = G.vertices[v1].r.y; + graph::coordinate r0 = G.vertices[v0].r; + graph::coordinate r1 = G.vertices[v1].r; - if (G.edges[e].crossings[1]) { - bool ind = v0y < v1y ? unbreak : !unbreak; + 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; } } -current_info network::get_current_info() { - cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); +network::network(const graph& G) : G(G), fuses(G.edges.size()), thresholds(G.edges.size()) {} - if (((double *)x->x)[0] != ((double *)x->x)[0]) { - throw nanex; - } +network::network(const network& n) : G(n.G), fuses(n.fuses), thresholds(n.thresholds) {} - cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c); +void network::set_thresholds(double beta, std::mt19937& rng) { + if (beta == 0.0) { + /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */ + for (long double& threshold : thresholds) { + threshold = 1.0; + } + } else { + std::uniform_real_distribution d(0.0, 1.0); - double alpha[2] = {1, 0}; - double beta[2] = {0, 0}; - CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); + for (long double& threshold : thresholds) { + threshold = std::numeric_limits::lowest(); - std::vector currents(G.edges.size()); + while (threshold == std::numeric_limits::lowest()) { + threshold = logl(d(rng)) / (long double)beta; + } + } + } +} - double total_current = 0; +fuse_network::fuse_network(const graph& G, cholmod_common *c) : network(G), ohm(G, 1, c) { +} - for (int i = 0; i < G.edges.size(); i++) { - if (fuses[i]) { - currents[i] = 0; - } else { - currents[i] = ((double *)y->x)[i]; +void fuse_network::fracture(hooks& m, double cutoff) { + m.pre_fracture(*this); - double v0y = G.vertices[G.edges[i].v[0]].r.y; - double v1y = G.vertices[G.edges[i].v[1]].r.y; + while (true) { + current_info ci = ohm.solve(fuses); - if (G.edges[i].crossings[1]) { - if (v0y > v1y) { - currents[i] += 1.0; - total_current += currents[i]; - } else { - currents[i] -= 1.0; - total_current -= currents[i]; + if (ci.conductivity < cutoff * G.vertices.size()) { + 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] && fabs(ci.currents[i]) > cutoff) { + long double val = logl(fabs(ci.currents[i])) - thresholds[i]; + if (val > max_val) { + max_val = val; + max_pos = i; } } } + if (max_pos == UINT_MAX) { + throw nofuseex; + } + + fuses[max_pos] = true; + ohm.break_edge(max_pos); + + m.bond_broken(*this, ci, max_pos); } - CHOL_F(free_dense)(&x, c); - CHOL_F(free_dense)(&y, c); + m.post_fracture(*this); +} - return {total_current, currents}; +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 network::fracture(hooks& m, double cutoff) { +void elastic_network::fracture(hooks& m, double weight, double cutoff) { m.pre_fracture(*this); while (true) { - current_info ci = this->get_current_info(); + current_info cx = hook_x.solve(fuses); + current_info cy = hook_y.solve(fuses); - if (ci.conductivity < cutoff * G.vertices.size()) { + current_info ctot; + + ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); + + if (ctot.conductivity < cutoff * G.vertices.size()) { break; } + ctot.currents.resize(G.edges.size()); + for (unsigned i = 0; i < G.edges.size(); i++) { + ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i], 2) + pow(weight * cy.currents[i], 2)); + } + 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] && fabs(ci.currents[i]) > cutoff) { - long double val = logl(fabs(ci.currents[i])) - thresholds[i]; + if (!fuses[i] && fabs(ctot.currents[i]) > cutoff) { + long double val = logl(fabs(ctot.currents[i])) - thresholds[i]; if (val > max_val) { max_val = val; max_pos = i; @@ -250,9 +326,11 @@ void network::fracture(hooks& m, double cutoff) { throw nofuseex; } - this->break_edge(max_pos); + fuses[max_pos] = true; + hook_x.break_edge(max_pos); + hook_y.break_edge(max_pos); - m.bond_broken(*this, ci, max_pos); + m.bond_broken(*this, ctot, max_pos); } m.post_fracture(*this); -- cgit v1.2.3-54-g00ecf