From 07906baa42470bad14d2c40f57967691f6118969 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 1 Nov 2018 12:33:37 -0400 Subject: revamped and simplied fracture code with c++ --- lib/src/graph.cpp | 38 +++++++ lib/src/network.cpp | 280 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 318 insertions(+) create mode 100644 lib/src/graph.cpp create mode 100644 lib/src/network.cpp (limited to 'lib/src') diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp new file mode 100644 index 0000000..e76947b --- /dev/null +++ b/lib/src/graph.cpp @@ -0,0 +1,38 @@ + +#include "graph.hpp" + +graph::graph(unsigned int L) { + double dx = 1.0 / L; + unsigned int nv = 2 * pow(L / 2, 2); + unsigned int ne = pow(L, 2); + + vertices.resize(nv); + edges.resize(ne); + + dual_vertices.resize(nv); + dual_edges.resize(ne); + + for (unsigned int i = 0; i < nv; i++) { + vertices[i].r.x = dx * ((1 + i / (L / 2)) % 2 + 2 * (i % (L / 2))); + vertices[i].r.y = dx * (i / (L / 2)); + + dual_vertices[i].r.x = dx * ((i / (L / 2)) % 2 + 2 * (i % (L / 2))); + dual_vertices[i].r.y = dx * (i / (L / 2)); + } + + for (unsigned int i = 0; i < ne; i++) { + unsigned int x = i / L; + unsigned int y = i % L; + + unsigned int v1 = (L * x) / 2 + ((y + x % 2) / 2) % (L / 2); + unsigned int v2 = ((L * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2)) % nv; + + edges[i] = {v1, v2}; + + unsigned int dv1 = (L * x) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2); + unsigned int dv2 = ((L * (x + 1)) / 2 + ((y + x % 2) / 2) % (L / 2)) % nv; + + dual_edges[i] = {dv1, dv2}; + } +} + diff --git a/lib/src/network.cpp b/lib/src/network.cpp new file mode 100644 index 0000000..8af140b --- /dev/null +++ b/lib/src/network.cpp @@ -0,0 +1,280 @@ + +#include "network.hpp" + +network::network(const graph& G, cholmod_common *c) : G(G), c(c), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) { + b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c); + for (unsigned int i = 0; i < G.edges.size(); i++) { + double v0y = G.vertices[G.edges[i][0]].r.y; + double v1y = G.vertices[G.edges[i][1]].r.y; + + if (fabs(v0y - v1y) > 0.5) { + bool ind = v0y < v1y ? 0 : 1; + + ((double *)b->x)[G.edges[i][ind]] += 1.0; + ((double *)b->x)[G.edges[i][!ind]] -= 1.0; + } + } + + unsigned int 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); + + t->nnz = nnz; + + for (unsigned int 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 int terms = G.vertices.size(); + + for (unsigned int i = 0; i < G.edges.size(); i++) { + unsigned int v0 = G.edges[i][0]; + unsigned int v1 = G.edges[i][1]; + + unsigned int s0 = v0 < v1 ? v0 : v1; + unsigned int s1 = v0 < v1 ? v1 : v0; + + ((CHOL_INT *)t->i)[terms] = v1; + ((CHOL_INT *)t->j)[terms] = v0; + ((double *)t->x)[terms] = -1.0; + + ((double *)t->x)[v0]++; + ((double *)t->x)[v1]++; + + terms++; + } + + ((double *)t->x)[0]++; + + cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, nnz, 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); + + 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 int i = 0; i < G.edges.size(); i++) { + ((CHOL_INT *)t->i)[2 * i] = i; + ((CHOL_INT *)t->j)[2 * i] = G.edges[i][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][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); +} + +network::network(const network &other) : G(other.G), thresholds(other.thresholds), fuses(other.fuses), 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); +} + +network::~network() { + CHOL_F(free_sparse)(&voltcurmat, c); + CHOL_F(free_dense)(&b, c); + CHOL_F(free_factor)(&factor, c); +} + +void network::set_thresholds(double beta, std::mt19937& rng) { + std::uniform_real_distribution d(0.0, 1.0); + + for (unsigned int i = 0; i < G.edges.size(); i++) { + long double x = 0.0; + + while (x == 0.0) { + x = expl(logl(d(rng)) / (long double)beta); + } + + thresholds[i] = x; + } +} + +void network::add_edge(unsigned int e) { + fuses[e] = false; + unsigned int v0 = G.edges[e][0]; + unsigned int v1 = G.edges[e][1]; + + unsigned int n = factor->n; + + cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); + + unsigned int 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 int i = 0; i <= s1; i++) { + pp[i] = 0; + } + + for (unsigned int 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)(false, perm_update_mat, factor, c); + + 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; + + if (fabs(v0y - v1y) > 0.5) { + bool ind = v0y < v1y ? 0 : 1; + + ((double *)b->x)[G.edges[e][ind]] += 1.0; + ((double *)b->x)[G.edges[e][!ind]] -= 1.0; + } +} + +void network::break_edge(unsigned int e) { + fuses[e] = true; + unsigned int v0 = G.edges[e][0]; + unsigned int v1 = G.edges[e][1]; + + unsigned int n = factor->n; + + cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); + + unsigned int 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 int i = 0; i <= s1; i++) { + pp[i] = 0; + } + + for (unsigned int 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)(false, perm_update_mat, factor, c); + + 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; + + if (fabs(v0y - v1y) > 0.5) { + bool ind = v0y < v1y ? 0 : 1; + + ((double *)b->x)[G.edges[e][ind]] -= 1.0; + ((double *)b->x)[G.edges[e][!ind]] += 1.0; + } +} + +std::pair> network::currents() { + cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); + + if (((double *)x->x)[0] != ((double *)x->x)[0]) { + exit(2); + } + + 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]; + + double v0y = G.vertices[G.edges[i][0]].r.y; + double v1y = G.vertices[G.edges[i][1]].r.y; + + if (fabs(v0y - v1y) > 0.5) { + if (v0y > v1y) { + 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 std::make_pair(total_current, currents); +} + +void network::fracture(hooks& m, double cutoff) { + m.pre_fracture(*this); + + while (true) { + auto currents = this->currents(); + + if (currents.first < cutoff) { + break; + } + + unsigned int max_pos = UINT_MAX; + long double max_val = 0; + + for (unsigned int i = 0; i < G.edges.size(); i++) { + if (!fuses[i] && fabs(currents.second[i]) > cutoff) { + long double val = (long double)fabs(currents.second[i]) / thresholds[i]; + if (val > max_val) { + max_val = val; + max_pos = i; + } + } + } + + if (max_pos == UINT_MAX) { + exit(3); + } + + this->break_edge(max_pos); + + m.bond_broken(*this, currents, max_pos); + } + + m.post_fracture(*this); +} + -- cgit v1.2.3-54-g00ecf