summaryrefslogtreecommitdiff
path: root/lib/src
diff options
context:
space:
mode:
Diffstat (limited to 'lib/src')
-rw-r--r--lib/src/graph.cpp38
-rw-r--r--lib/src/network.cpp280
2 files changed, 318 insertions, 0 deletions
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<long double> 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<double, std::vector<double>> 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<double> 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);
+}
+