revamped and simplied fracture code with c++

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);
+}
+