summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/CMakeLists.txt4
-rw-r--r--lib/bound_set.c110
-rw-r--r--lib/break_edge.c50
-rw-r--r--lib/correlations.c46
-rw-r--r--lib/data.c35
-rw-r--r--lib/factor_update.c69
-rw-r--r--lib/fracture.h123
-rw-r--r--lib/gen_laplacian.c142
-rw-r--r--lib/gen_voltcurmat.c36
-rw-r--r--lib/geometry.c55
-rw-r--r--lib/get_dual_clusters.c31
-rw-r--r--lib/include/graph.hpp29
-rw-r--r--lib/include/hooks.hpp12
-rw-r--r--lib/include/network.hpp48
-rw-r--r--lib/net.c146
-rw-r--r--lib/net_conductivity.c35
-rw-r--r--lib/net_currents.c52
-rw-r--r--lib/net_fracture.c68
-rw-r--r--lib/net_voltages.c39
-rw-r--r--lib/rand.c15
-rw-r--r--lib/src/graph.cpp38
-rw-r--r--lib/src/network.cpp280
22 files changed, 411 insertions, 1052 deletions
diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
new file mode 100644
index 0000000..bf3c3b5
--- /dev/null
+++ b/lib/CMakeLists.txt
@@ -0,0 +1,4 @@
+
+add_library(frac SHARED src/graph.cpp src/network.cpp)
+target_include_directories(frac PUBLIC include)
+
diff --git a/lib/bound_set.c b/lib/bound_set.c
deleted file mode 100644
index 32d0fb7..0000000
--- a/lib/bound_set.c
+++ /dev/null
@@ -1,110 +0,0 @@
-
-#include "fracture.h"
-
-double th_p(double x, double y, double th) {
- if (x >= 0 && y >= 0)
- return th;
- else if (x < 0 && y >= 0)
- return M_PI - th;
- else if (x < 0 && y < 0)
- return th - M_PI;
- else
- return -th;
-}
-
-double u_y(double x, double y) {
- double r = sqrt(pow(x, 2) + pow(y, 2));
- double th = th_p(x, y, atan(fabs(y / x)));
-
- return sqrt(r) * sin(th / 2);
-}
-
-void bound_set_embedded(double *bound, const graph_t *g, double notch_len) {
- uint_t L = g->L;
-
- for (uint_t i = 0; i < L / 2; i++) {
- double x1, y1, x2, y2, x3, y3, x4, y4;
- x1 = (2. * i + 1.) / L - notch_len;
- y1 = 0.5 - 1.;
- x2 = (2. * i + 1.) / L - notch_len;
- y2 = 0.5 - 0.;
- y3 = (2. * i + 1.) / L - 0.5;
- x3 = 0.5 - 1.;
- y4 = (2. * i + 1.) / L - 0.5;
- x4 = 0.5 - 0.;
-
- bound[g->b[g->bi[0] + i]] = u_y(x1, y1);
- bound[g->b[g->bi[1] + i]] = u_y(x2, y2);
- bound[g->b[g->bi[2] + i]] = u_y(x3, y3);
- bound[g->b[g->bi[3] + i]] = u_y(x4, y4);
- }
-}
-
-bool is_in(uint_t len, uint_t *list, uint_t element) {
- for (uint_t i = 0; i < len; i++) {
- if (list[i] == element) {
- return true;
- }
- }
- return false;
-}
-
-cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len,
- cholmod_common *c) {
-
- uint_t dim = g->nv;
-
- if (vb && g->boundary != TORUS_BOUND) {
- dim -= g->bi[g->nb];
- } else if (!vb) {
- dim += 2;
- }
-
- cholmod_dense *boundary = CHOL_F(zeros)(dim, 1, CHOLMOD_REAL, c);
- double *bound = (double *)boundary->x;
-
- switch (g->boundary) {
- case TORUS_BOUND:
- for (uint_t i = 0; i < g->bi[1]; i++) {
- uint_t be = g->b[i];
- uint_t v1 = g->ev[2 * be];
- uint_t v2 = g->ev[2 * be + 1];
- double v1y = g->vx[2 * v1 + 1];
- double v2y = g->vx[2 * v2 + 1];
-
- uint_t ind = v1y < v2y ? 0 : 1;
-
- bound[g->ev[2 * be + ind]] += 1;
- bound[g->ev[2 * be + !ind]] -= 1;
- }
- break;
- /*
-case EMBEDDED_BOUND:
- bound_set_embedded(bound, g, notch_len);
- break;
- */
- default:
- if (vb) {
- for (uint_t i = 0; i < dim; i++) {
- uint_t v = g->nbi[i];
- for (uint_t j = 0; j < g->vei[v + 1] - g->vei[v]; j++) {
- uint_t e = g->ve[g->vei[v] + j];
- uint_t v0 = g->ev[2 * e];
- uint_t v1 = g->ev[2 * e + 1];
-
- if (g->bq[v0] || g->bq[v1]) {
- uint_t vv = v0 == v ? v1 : v0;
- if (is_in(g->bi[1], g->b, vv)) {
- bound[i]++;
- }
- }
- }
- }
- } else {
- bound[g->nv] = 1;
- bound[g->nv + 1] = -1;
- }
- }
-
- return boundary;
-}
diff --git a/lib/break_edge.c b/lib/break_edge.c
deleted file mode 100644
index e53f7c1..0000000
--- a/lib/break_edge.c
+++ /dev/null
@@ -1,50 +0,0 @@
-
-#include "fracture.h"
-
-bool break_edge(net_t *net, uint_t edge, cholmod_common *c) {
- net->fuses[edge] = true;
- net->num_broken++;
-
- double *x = (double *)net->boundary_cond->x;
- const graph_t *g = net->graph;
-
- uint_t v1 = net->graph->ev[2 * edge];
- uint_t v2 = net->graph->ev[2 * edge + 1];
-
- uint_t s1 = v1 < v2 ? v1 : v2;
- uint_t s2 = v1 < v2 ? v2 : v1;
-
- if (net->graph->boundary == TORUS_BOUND) {
- if (net->graph->bq[edge]) {
- double v1y = g->vx[2 * v1 + 1];
- double v2y = g->vx[2 * v2 + 1];
- uint_t ind = v1y < v2y ? 0 : 1;
- x[g->ev[2 * edge + ind]] -= 1;
- x[g->ev[2 * edge + !ind]] += 1;
- }
- if (net->factor != NULL) {
- factor_update(net->factor, s1, s2, c);
- }
- } else if (net->voltage_bound) {
- if (g->bq[v1] || g->bq[v2]) {
- uint_t vv = g->bq[v1] ? v2 : v1;
- uint_t uu = v1 == vv ? v2 : v1;
- if (is_in(g->bi[1], g->b, uu)) {
- x[g->bni[vv]] -= 1;
- }
- if (net->factor != NULL) {
- factor_update2(net->factor, g->bni[vv], c);
- }
- } else {
- if (net->factor != NULL) {
- factor_update(net->factor, g->bni[s1], g->bni[s2], c);
- }
- }
- } else {
- if (net->factor != NULL) {
- factor_update(net->factor, s1, s2, c);
- }
- }
-
- return true;
-}
diff --git a/lib/correlations.c b/lib/correlations.c
deleted file mode 100644
index 98c5767..0000000
--- a/lib/correlations.c
+++ /dev/null
@@ -1,46 +0,0 @@
-
-#include "fracture.h"
-
-double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c) {
- uint_t nv = instance->graph->dnv;
- uint_t ne = instance->graph->ne;
- uint_t *ev = instance->graph->dev;
- bool nulldists = false;
- if (dists == NULL) {
- dists = graph_dists(instance->graph, false);
- nulldists = true;
- }
- double *corr = calloc(nv, sizeof(double));
- uint_t *numat = calloc(nv, sizeof(uint_t));
-
- for (uint_t i = 0; i < ne; i++) {
- uint_t v1 = ev[2 * i];
- uint_t v2 = ev[2 * i + 1];
- for (uint_t j = 0; j < ne; j++) {
- uint_t v3 = ev[2 * j];
- uint_t v4 = ev[2 * j + 1];
- uint_t dist1 = dists[v1][v3];
- uint_t dist2 = dists[v1][v4];
- uint_t dist3 = dists[v2][v3];
- uint_t dist4 = dists[v2][v4];
- uint_t dist = (dist1 + dist2 + dist3 + dist4) / 4;
- corr[dist] += instance->fuses[i] && instance->fuses[j];
- numat[dist]++;
- }
- }
-
- for (uint_t i = 0; i < nv; i++) {
- if (numat[i] > 0) {
- corr[i] /= numat[i];
- }
- }
-
- if (nulldists) {
- for (int i = 0; i < nv; i++) {
- free(dists[i]);
- }
- free(dists);
- }
-
- return corr;
-}
diff --git a/lib/data.c b/lib/data.c
deleted file mode 100644
index e39c450..0000000
--- a/lib/data.c
+++ /dev/null
@@ -1,35 +0,0 @@
-
-#include "fracture.h"
-
-data_t *data_create(uint_t ne) {
- data_t *data = malloc(1 * sizeof(data_t));
- assert(data != NULL);
-
- data->num_broken = 0;
-
- data->break_list = (uint_t *)malloc(ne * sizeof(uint_t));
- assert(data->break_list != NULL);
-
- data->extern_field = (long double *)malloc(ne * sizeof(long double));
- assert(data->extern_field != NULL);
-
- data->conductivity = (double *)malloc(ne * sizeof(double));
- assert(data->conductivity != NULL);
-
- return data;
-}
-
-void data_free(data_t *data) {
- free(data->break_list);
- free(data->extern_field);
- free(data->conductivity);
- free(data);
-}
-
-void data_update(data_t *data, uint_t last_broke, long double strength,
- double conductivity) {
- data->break_list[data->num_broken] = last_broke;
- data->extern_field[data->num_broken] = strength;
- data->conductivity[data->num_broken] = conductivity;
- data->num_broken++;
-}
diff --git a/lib/factor_update.c b/lib/factor_update.c
deleted file mode 100644
index ceaa01a..0000000
--- a/lib/factor_update.c
+++ /dev/null
@@ -1,69 +0,0 @@
-
-#include "fracture.h"
-
-void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2,
- cholmod_common *c) {
- uint_t n = factor->n;
-
- cholmod_sparse *update_mat =
- CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
-
- uint_t s1, s2;
- s1 = v1 < v2 ? v1 : v2;
- s2 = v1 > v2 ? v1 : v2;
-
- int_t *pp = (int_t *)update_mat->p;
- int_t *ii = (int_t *)update_mat->i;
- double *xx = (double *)update_mat->x;
-
- for (uint_t i = 0; i <= s1; i++) {
- pp[i] = 0;
- }
-
- for (uint_t 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, 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);
-}
-
-void factor_update2(cholmod_factor *factor, uint_t v, cholmod_common *c) {
- uint_t n = factor->n;
-
- cholmod_sparse *update_mat =
- CHOL_F(allocate_sparse)(n, n, 1, true, true, 0, CHOLMOD_REAL, c);
-
- int_t *pp = (int_t *)update_mat->p;
- int_t *ii = (int_t *)update_mat->i;
- double *xx = (double *)update_mat->x;
-
- for (uint_t i = 0; i <= v; i++) {
- pp[i] = 0;
- }
-
- for (uint_t i = v + 1; i <= n; i++) {
- pp[i] = 1;
- }
-
- ii[0] = v;
- xx[0] = 1;
-
- cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
- update_mat, 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);
-}
diff --git a/lib/fracture.h b/lib/fracture.h
deleted file mode 100644
index 5eb0a1d..0000000
--- a/lib/fracture.h
+++ /dev/null
@@ -1,123 +0,0 @@
-
-#pragma once
-
-#include <assert.h>
-#include <cholmod.h>
-#include <float.h>
-#include <getopt.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_randist.h>
-#include <gsl/gsl_rng.h>
-#include <gsl/gsl_sf_erf.h>
-#include <gsl/gsl_sf_exp.h>
-#include <gsl/gsl_sf_log.h>
-#include <inttypes.h>
-#include <math.h>
-#include <stdbool.h>
-#include <stdint.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <sys/types.h>
-#include <unistd.h>
-
-#include <jst/graph.h>
-#include <jst/rand.h>
-
-// these defs allow me to switch to long int cholmod in a sitch
-#define int_t int
-#define uint_t unsigned int
-#define CINT_MAX INT_MAX
-#define CHOL_F(x) cholmod_##x
-
-#define GSL_RAND_GEN gsl_rng_mt19937
-
-typedef struct {
- const graph_t *graph;
- bool *fuses;
- long double *thres;
- double inf;
- cholmod_dense *boundary_cond;
- cholmod_factor *factor;
- bool voltage_bound;
- uint_t num_broken;
- uint_t dim;
- uint_t nep;
- uint_t *evp;
- cholmod_sparse *voltcurmat;
-} net_t;
-
-typedef struct {
- uint_t num_broken;
- uint_t *break_list;
- double *conductivity;
- long double *extern_field;
-} data_t;
-
-intptr_t *run_voronoi(uint_t num_coords, double *coords, bool periodic,
- double xmin, double xmax, double ymin, double ymax);
-
-cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp,
- bool symmetric, cholmod_common *c);
-
-cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c);
-
-int edge_to_verts(uint_t width, bool periodic, uint_t edge, bool index);
-
-int dual_edge_to_verts(uint_t width, bool periodic, uint_t edge, bool index);
-
-double dual_vert_to_coord(uint_t width, bool periodic, uint_t vert, bool index);
-
-void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2,
- cholmod_common *c);
-void factor_update2(cholmod_factor *factor, uint_t v1, cholmod_common *c);
-
-void net_notch(net_t *net, double notch_len, cholmod_common *c);
-data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff);
-double *net_voltages(const net_t *net, cholmod_common *c);
-double *net_currents(const net_t *net, const double *voltages,
- cholmod_common *c);
-double net_conductivity(const net_t *net, const double *voltages);
-
-void update_boundary(net_t *instance, const double *avg_field);
-
-FILE *get_file(const char *prefix, uint_t width, uint_t crack, double beta,
- uint_t iter, uint_t num_iter, uint_t num, bool read);
-
-double update_beta(double beta, uint_t width, const double *stress,
- const double *damage, double bound_total);
-
-cholmod_sparse *gen_voltcurmat(uint_t num_edges, uint_t num_verts,
- uint_t *edges_to_verts, cholmod_common *c);
-
-net_t *net_copy(const net_t *net, cholmod_common *c);
-
-void net_free(net_t *instance, cholmod_common *c);
-
-net_t *net_create(const graph_t *g, double inf, double beta, double notch_len,
- bool vb, cholmod_common *c);
-
-bool break_edge(net_t *instance, uint_t edge, cholmod_common *c);
-
-components_t *get_clusters(net_t *instance);
-
-uint_t *get_cluster_dist(net_t *instance);
-
-void randfunc_flat(gsl_rng *r, double *x, double *y);
-void randfunc_gaus(gsl_rng *r, double *x, double *y);
-
-double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c);
-
-double *bin_values(graph_t *network, uint_t width, double *values);
-
-cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len,
- cholmod_common *c);
-
-data_t *data_create(uint_t num_edges);
-void data_free(data_t *data);
-void data_update(data_t *data, uint_t last_broke, long double strength,
- double conductivity);
-
-long double rand_dist_pow(const gsl_rng *r, double beta);
-
-bool is_in(uint_t len, uint_t *list, uint_t element);
diff --git a/lib/gen_laplacian.c b/lib/gen_laplacian.c
deleted file mode 100644
index 75dccce..0000000
--- a/lib/gen_laplacian.c
+++ /dev/null
@@ -1,142 +0,0 @@
-
-#include "fracture.h"
-
-cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp,
- bool symmetric, cholmod_common *c) {
- const graph_t *g = net->graph;
- uint_t nv;
- uint_t ne;
- uint_t nre;
- uint_t *ev;
-
- if (use_gp) {
- nv = net->dim;
- ne = net->nep;
- nre = (int_t)net->nep - (int_t)net->num_broken;
- ev = net->evp;
- } else if (dual) {
- nv = g->dnv;
- ne = g->ne;
- nre = net->num_broken;
- ev = g->dev;
- } else {
- nv = g->nv;
- ne = g->ne;
- nre = (int_t)g->ne - (int_t)net->num_broken;
- ev = g->ev;
- }
-
- uint_t nnz = nre;
-
- cholmod_triplet *t =
- CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
-
- int_t *ri = (int_t *)t->i;
- int_t *ci = (int_t *)t->j;
- double *ai = (double *)t->x;
-
- t->nnz = nnz;
-
- uint_t a = 0;
-
- for (uint_t i = 0; i < ne; i++) {
- if ((net->fuses[i] && dual) || (!net->fuses[i] && !dual)) {
- uint_t v1 = ev[2 * i];
- uint_t v2 = ev[2 * i + 1];
-
- uint_t s1 = v1 < v2 ? v1 : v2;
- uint_t s2 = v1 < v2 ? v2 : v1;
-
- ri[a] = s2;
- ci[a] = s1;
- ai[a] = 1;
-
- a++;
- }
- }
-
- cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c);
- CHOL_F(free_triplet)(&t, c);
-
- if (!symmetric) {
- cholmod_sparse *tmp_s = CHOL_F(copy)(s, 0, 1, c);
- CHOL_F(free_sparse)(&s, c);
- s = tmp_s;
- }
-
- return s;
-}
-
-cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c) {
- const graph_t *g = net->graph;
- uint_t nv = net->dim;
- uint_t ne = net->nep;
- uint_t *ev = net->evp;
-
- uint_t nnz = nv;
-
- cholmod_triplet *temp_m =
- CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
-
- int_t *rowind = (int_t *)temp_m->i;
- int_t *colind = (int_t *)temp_m->j;
- double *acoo = (double *)temp_m->x;
-
- temp_m->nnz = nnz;
-
- for (uint_t i = 0; i < nv; i++) {
- rowind[i] = i;
- colind[i] = i;
- acoo[i] = 0;
- }
-
- cholmod_sparse *adjacency = gen_adjacency(net, false, true, true, c);
-
- for (uint_t i = 0; i < ne; i++) {
- if (!net->fuses[i]) {
- uint_t v1 = ev[2 * i];
- uint_t v2 = ev[2 * i + 1];
-
- acoo[v1]++;
- acoo[v2]++;
- }
- }
-
- if (net->voltage_bound && g->boundary != TORUS_BOUND) {
- for (uint_t i = 0; i < net->dim; i++) {
- uint_t v = g->nbi[i];
- for (uint_t j = 0; j < g->vei[v + 1] - g->vei[v]; j++) {
- uint_t e = g->ve[g->vei[v] + j];
- uint_t v0 = g->ev[2 * e];
- uint_t v1 = g->ev[2 * e + 1];
-
- if (g->bq[v0] || g->bq[v1]) {
- acoo[i]++;
- }
- }
- }
- } else {
- acoo[0]++;
- }
-
- for (uint_t i = 0; i < nv; i++) {
- if (acoo[i] == 0)
- acoo[i]++;
- }
-
- // assert(CHOL_F(check_triplet)(temp_m, c));
-
- cholmod_sparse *t_out = CHOL_F(triplet_to_sparse)(temp_m, nnz, c);
- // assert(CHOL_F(check_sparse)(t_out, c));
-
- double alpha[2] = {1, 0};
- double beta[2] = {-1, 0};
- cholmod_sparse *laplacian =
- CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c);
-
- CHOL_F(free_sparse)(&t_out, c);
- CHOL_F(free_sparse)(&adjacency, c);
- CHOL_F(free_triplet)(&temp_m, c);
-
- return laplacian;
-}
diff --git a/lib/gen_voltcurmat.c b/lib/gen_voltcurmat.c
deleted file mode 100644
index fede836..0000000
--- a/lib/gen_voltcurmat.c
+++ /dev/null
@@ -1,36 +0,0 @@
-
-#include "fracture.h"
-
-cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts,
- unsigned int *edges_to_verts,
- cholmod_common *c) {
-
- cholmod_triplet *t_m = CHOL_F(allocate_triplet)(
- num_edges, num_verts, num_edges * 2, 0, CHOLMOD_REAL, c);
- assert(t_m != NULL);
-
- int_t *rowind = (int_t *)t_m->i;
- int_t *colind = (int_t *)t_m->j;
- double *acoo = (double *)t_m->x;
-
- for (unsigned int i = 0; i < num_edges; i++) {
- unsigned int v1 = edges_to_verts[2 * i];
- unsigned int v2 = edges_to_verts[2 * i + 1];
- rowind[2 * i] = i;
- rowind[2 * i + 1] = i;
- colind[2 * i] = v1;
- colind[2 * i + 1] = v2;
- acoo[2 * i] = 1;
- acoo[2 * i + 1] = -1;
- }
-
- t_m->nnz = num_edges * 2;
-
- assert(CHOL_F(check_triplet)(t_m, c));
-
- cholmod_sparse *m = CHOL_F(triplet_to_sparse)(t_m, num_edges * 2, c);
-
- CHOL_F(free_triplet)(&t_m, c);
-
- return m;
-}
diff --git a/lib/geometry.c b/lib/geometry.c
deleted file mode 100644
index 2c1987b..0000000
--- a/lib/geometry.c
+++ /dev/null
@@ -1,55 +0,0 @@
-
-#include "fracture.h"
-
-int edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
- bool index) {
- assert(edge < pow(width, 2));
-
- int x = edge / width + 1;
- int y = edge % width + 1;
-
- if (periodic) {
- return (((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) - 1) %
- width +
- (x - index) * width) /
- 2;
- } else {
- return ((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) +
- (x - index) * (width + 1) - 1) /
- 2;
- }
-}
-
-int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
- bool index) {
- assert(edge < pow(width, 2));
-
- int x = edge / width + 1;
- int y = edge % width + 1;
-
- if (periodic) {
- return (((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) - 1) %
- width +
- (x - index) * width) /
- 2;
- } else {
- return ((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) +
- (x - index) * (width + 1) - 1) /
- 2;
- }
-}
-
-double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert,
- bool index) {
- if (periodic) {
- if (index)
- return (2 * vert) % width + (2 * vert / width) % 2;
- else
- return 2 * vert / width;
- } else {
- if (index)
- return (2 * vert) % (width + 1);
- else
- return (2 * vert) / (width + 1);
- }
-}
diff --git a/lib/get_dual_clusters.c b/lib/get_dual_clusters.c
deleted file mode 100644
index 9336106..0000000
--- a/lib/get_dual_clusters.c
+++ /dev/null
@@ -1,31 +0,0 @@
-
-#include "fracture.h"
-
-components_t *get_clusters(net_t *instance) {
- components_t *c =
- graph_components_get(instance->graph, instance->fuses, true);
- return c;
-}
-
-unsigned int *get_cluster_dist(net_t *instance) {
- components_t *c = get_clusters(instance);
- unsigned int *cluster_dist =
- (unsigned int *)calloc(instance->graph->dnv, sizeof(unsigned int));
-
- for (uint32_t i = 1; i <= c->n; i++) {
- unsigned int num_in_cluster = 0;
- for (unsigned int j = 0; j < instance->graph->dnv; j++) {
- if (c->labels[j] == i)
- num_in_cluster++;
- }
-
- if (num_in_cluster == 0)
- break;
-
- cluster_dist[num_in_cluster - 1]++;
- }
-
- graph_components_free(c);
-
- return cluster_dist;
-}
diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
new file mode 100644
index 0000000..80cdd66
--- /dev/null
+++ b/lib/include/graph.hpp
@@ -0,0 +1,29 @@
+
+#pragma once
+
+#include <cmath>
+#include <vector>
+#include <array>
+
+class graph {
+ public:
+ typedef struct coordinate {
+ double x;
+ double y;
+ } coordinate;
+
+ typedef struct vertex {
+ coordinate r;
+ } vertex;
+
+ typedef std::array<unsigned int, 2> edge;
+
+ std::vector<vertex> vertices;
+ std::vector<edge> edges;
+
+ std::vector<vertex> dual_vertices;
+ std::vector<edge> dual_edges;
+
+ graph(unsigned int L);
+};
+
diff --git a/lib/include/hooks.hpp b/lib/include/hooks.hpp
new file mode 100644
index 0000000..350f318
--- /dev/null
+++ b/lib/include/hooks.hpp
@@ -0,0 +1,12 @@
+
+#pragma once
+
+class network;
+
+class hooks {
+ public:
+ virtual void pre_fracture(const network&) {};
+ virtual void bond_broken(const network&, const std::pair<double, std::vector<double>>&, unsigned int) {};
+ virtual void post_fracture(network&) {}; // post fracture hook can be destructive
+};
+
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
new file mode 100644
index 0000000..abf88cd
--- /dev/null
+++ b/lib/include/network.hpp
@@ -0,0 +1,48 @@
+
+#pragma once
+
+#include <vector>
+#include <functional>
+#include <utility>
+#include <random>
+
+#include "cholmod.h"
+
+#include "graph.hpp"
+#include "hooks.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 network {
+ private:
+ cholmod_dense *b;
+ cholmod_factor *factor;
+ cholmod_sparse *voltcurmat;
+ cholmod_common *c;
+
+ public:
+ const graph& G;
+ std::vector<bool> fuses;
+ std::vector<long double> thresholds;
+
+ network(const graph&, cholmod_common*);
+ network(const network &other);
+ ~network();
+
+ void set_thresholds(double, std::mt19937&);
+ void break_edge(unsigned int);
+ void add_edge(unsigned int);
+ std::pair<double, std::vector<double>> currents();
+ void fracture(hooks&, double cutoff = 1e-13);
+};
+
diff --git a/lib/net.c b/lib/net.c
deleted file mode 100644
index ac89415..0000000
--- a/lib/net.c
+++ /dev/null
@@ -1,146 +0,0 @@
-
-#include "fracture.h"
-
-long double *get_thres(uint_t ne, double beta) {
- long double *thres = (long double *)malloc(ne * sizeof(long double));
- assert(thres != NULL);
-
- gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN);
- gsl_rng_set(r, jst_rand_seed());
-
- for (uint_t i = 0; i < ne; i++) {
- thres[i] = rand_dist_pow(r, beta);
- }
-
- gsl_rng_free(r);
-
- return thres;
-}
-
-void net_notch(net_t *net, double notch_len, cholmod_common *c) {
- for (uint_t i = 0; i < net->graph->ne; i++) {
- uint_t v1, v2;
- double v1x, v1y, v2x, v2y, dy;
- bool crosses_center, not_wrapping, correct_length;
-
- v1 = net->graph->ev[2 * i];
- v2 = net->graph->ev[2 * i + 1];
-
- v1x = net->graph->vx[2 * v1];
- v1y = net->graph->vx[2 * v1 + 1];
- v2x = net->graph->vx[2 * v2];
- v2y = net->graph->vx[2 * v2 + 1];
-
- dy = v1y - v2y;
-
- crosses_center = (v1y >= 0.5 && v2y <= 0.5) || (v1y <= 0.5 && v2y >= 0.5);
- not_wrapping = fabs(dy) < 0.5;
- // correct_length = v1x + dx / dy * (v1y - 0.5) <= notch_len;
- correct_length = v1x < notch_len && v2x < notch_len;
-
- if (crosses_center && not_wrapping && correct_length) {
- break_edge(net, i, c);
- }
- }
-}
-
-net_t *net_create(const graph_t *g, double inf, double beta, double notch_len,
- bool vb, cholmod_common *c) {
- net_t *net = (net_t *)calloc(1, sizeof(net_t));
- assert(net != NULL);
-
- net->graph = g;
- net->num_broken = 0;
-
- net->fuses = (bool *)calloc(g->ne, sizeof(bool));
- assert(net->fuses != NULL);
-
- net->thres = get_thres(g->ne, beta);
- net->inf = inf;
-
- net->dim = g->nv;
-
- if (g->boundary == TORUS_BOUND) {
- net->nep = g->ne;
- net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
- memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
- } else {
- if (vb) {
- net->dim -= g->bi[g->nb];
- net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
- net->nep = 0;
- for (uint_t i = 0; i < g->ne; i++) {
- if (!(g->bq[g->ev[2 * i]] || g->bq[g->ev[2 * i + 1]])) {
- net->evp[2 * net->nep] = g->bni[g->ev[2 * i]];
- net->evp[2 * net->nep + 1] = g->bni[g->ev[2 * i + 1]];
- net->nep++;
- }
- }
- } else {
- net->dim += 2;
- net->evp = (uint_t *)malloc(2 * (g->ne + g->bi[2]) * sizeof(uint_t));
- memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
- net->nep = g->ne + g->bi[2];
-
- for (uint_t i = 0; i < 2; i++) {
- for (uint_t j = 0; j < g->bi[i + 1] - g->bi[i]; j++) {
- net->evp[2 * (g->ne + g->bi[i] + j)] = g->b[g->bi[i] + j];
- net->evp[2 * (g->ne + g->bi[i] + j) + 1] = g->nv + i;
- }
- }
- }
- }
-
- net->voltage_bound = vb;
- net->boundary_cond = bound_set(g, vb, notch_len, c);
-
- net_notch(net, notch_len, c);
-
- {
- cholmod_sparse *laplacian = gen_laplacian(net, c);
- net->factor = CHOL_F(analyze)(laplacian, c);
- CHOL_F(factorize)(laplacian, net->factor, c);
- CHOL_F(free_sparse)(&laplacian, c);
- }
-
- net->voltcurmat = gen_voltcurmat(g->ne, g->nv, g->ev, c);
-
- return net;
-}
-
-net_t *net_copy(const net_t *net, cholmod_common *c) {
- net_t *net_copy = (net_t *)calloc(1, sizeof(net_t));
- assert(net_copy != NULL);
- memcpy(net_copy, net, sizeof(net_t));
-
- size_t fuses_size = (net->graph)->ne * sizeof(bool);
- net_copy->fuses = (bool *)malloc(fuses_size);
- assert(net_copy->fuses != NULL);
- memcpy(net_copy->fuses, net->fuses, fuses_size);
-
- size_t thres_size = (net->graph)->ne * sizeof(long double);
- net_copy->thres = (long double *)malloc(thres_size);
- assert(net_copy->thres != NULL);
- memcpy(net_copy->thres, net->thres, thres_size);
-
- size_t evp_size = 2 * net->nep * sizeof(uint_t);
- net_copy->evp = (uint_t *)malloc(thres_size);
- assert(net_copy->evp != NULL);
- memcpy(net_copy->evp, net->evp, evp_size);
-
- net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c);
- net_copy->factor = CHOL_F(copy_factor)(net->factor, c);
- net_copy->voltcurmat = CHOL_F(copy_sparse)(net->voltcurmat, c);
-
- return net_copy;
-}
-
-void net_free(net_t *net, cholmod_common *c) {
- free(net->fuses);
- free(net->thres);
- CHOL_F(free_dense)(&(net->boundary_cond), c);
- CHOL_F(free_factor)(&(net->factor), c);
- CHOL_F(free_sparse)(&(net->voltcurmat), c);
- free(net->evp);
- free(net);
-}
diff --git a/lib/net_conductivity.c b/lib/net_conductivity.c
deleted file mode 100644
index 61148da..0000000
--- a/lib/net_conductivity.c
+++ /dev/null
@@ -1,35 +0,0 @@
-
-#include "fracture.h"
-
-double net_conductivity(const net_t *net, const double *voltages) {
- if (net->voltage_bound) {
- // the voltage drop across the network is fixed to one with voltage
- // boundary conditions, so the conductivity is the total current flowing
- double tot_cur = 0;
- for (uint_t i = 0; i < net->graph->num_spanning_edges; i++) {
- uint_t e = net->graph->spanning_edges[i];
-
- if (!net->fuses[e]) {
- uint_t v1, v2, s1, s2;
- double v1y, v2y;
-
- v1 = net->graph->ev[2 * e];
- v2 = net->graph->ev[2 * e + 1];
-
- v1y = net->graph->vx[2 * v1 + 1];
- v2y = net->graph->vx[2 * v2 + 1];
-
- s1 = v1y < v2y ? v1 : v2;
- s2 = v1y < v2y ? v2 : v1;
-
- tot_cur += voltages[s1] - voltages[s2];
- }
- }
-
- return fabs(tot_cur);
- } else {
- // the current across the network is fixed to one with current boundary
- // conditions, so the conductivity is the inverse of the total voltage drop
- return 1 / fabs(voltages[net->graph->nv] - voltages[net->graph->nv + 1]);
- }
-}
diff --git a/lib/net_currents.c b/lib/net_currents.c
deleted file mode 100644
index 32dba04..0000000
--- a/lib/net_currents.c
+++ /dev/null
@@ -1,52 +0,0 @@
-
-#include "fracture.h"
-
-double *net_currents(const net_t *net, const double *voltages,
- cholmod_common *c) {
- uint_t ne = net->graph->ne;
- uint_t dim = net->graph->nv;
- cholmod_sparse *voltcurmat = net->voltcurmat;
-
- cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c);
-
- double *tmp_x = x->x;
- x->x = (void *)voltages;
-
- cholmod_dense *y = CHOL_F(allocate_dense)(ne, 1, ne, CHOLMOD_REAL, c);
-
- double alpha[2] = {1, 0};
- double beta[2] = {0, 0};
- CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
-
- double *currents = (double *)malloc(ne * sizeof(double));
-
- for (int i = 0; i < ne; i++) {
- if (net->fuses[i]) {
- currents[i] = 0;
- } else {
- currents[i] = ((double *)y->x)[i];
- }
- }
-
- if (net->graph->boundary == TORUS_BOUND) {
- for (uint_t i = 0; i < net->graph->bi[1]; i++) {
- uint_t e = net->graph->b[i];
- uint_t v1 = net->graph->ev[2 * e];
- uint_t v2 = net->graph->ev[2 * e + 1];
- double v1y = net->graph->vx[2 * v1 + 1];
- double v2y = net->graph->vx[2 * v2 + 1];
-
- if (v1y > v2y) {
- currents[e] += 1;
- } else {
- currents[e] -= 1;
- }
- }
- }
-
- x->x = tmp_x;
- CHOL_F(free_dense)(&x, c);
- CHOL_F(free_dense)(&y, c);
-
- return currents;
-}
diff --git a/lib/net_fracture.c b/lib/net_fracture.c
deleted file mode 100644
index 65ede9b..0000000
--- a/lib/net_fracture.c
+++ /dev/null
@@ -1,68 +0,0 @@
-
-#include "fracture.h"
-
-uint_t get_next_broken(net_t *net, double *currents, double cutoff) {
- uint_t max_pos = UINT_MAX;
- long double max_val = 0;
-
- for (uint_t i = 0; i < net->graph->ne; i++) {
- long double current = fabs(currents[i]);
- bool broken = net->fuses[i];
-
- if (!broken && current > cutoff) {
- long double val = current / net->thres[i];
-
- if (val > max_val) {
- max_val = val;
- max_pos = i;
- }
- }
- }
-
- if (max_pos == UINT_MAX) {
- printf("GET_NEXT_BROKEN: currents is zero or NaN, no max_val found\n");
- exit(EXIT_FAILURE);
- }
-
- return max_pos;
-}
-
-data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
- data_t *data = data_create(net->graph->ne);
-
- uint_t n = 0;
- while (true) {
- n++;
- double *voltages = net_voltages(net, c);
- double *currents = net_currents(net, voltages, c);
-
- double conductivity = net_conductivity(net, voltages);
-
- if (conductivity < cutoff) {
- free(voltages);
- free(currents);
- break;
- }
-
- uint_t last_broke = get_next_broken(net, currents, cutoff);
-
- long double sim_current;
-
- if (net->voltage_bound) {
- sim_current = conductivity;
- } else {
- sim_current = 1;
- }
-
- data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] /
- ((long double)currents[last_broke])),
- conductivity);
-
- free(voltages);
- free(currents);
-
- break_edge(net, last_broke, c);
- }
-
- return data;
-}
diff --git a/lib/net_voltages.c b/lib/net_voltages.c
deleted file mode 100644
index d456a65..0000000
--- a/lib/net_voltages.c
+++ /dev/null
@@ -1,39 +0,0 @@
-
-#include "fracture.h"
-
-double *net_voltages(const net_t *net, cholmod_common *c) {
- cholmod_dense *b = net->boundary_cond;
- cholmod_factor *factor = net->factor;
-
- cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
-
- if (((double *)x->x)[0] != ((double *)x->x)[0]) {
- printf("GET_VOLTAGE: value is NaN\n");
- exit(EXIT_FAILURE);
- }
-
- double *t_voltages = (double *)x->x;
- x->x = NULL;
- CHOL_F(free_dense)(&x, c);
-
- const graph_t *g = net->graph;
-
- if (g->boundary == TORUS_BOUND) {
- return t_voltages;
- } else if (net->voltage_bound) {
- double *voltages = (double *)malloc(g->nv * sizeof(double));
- for (uint_t i = 0; i < g->nv - g->bi[g->nb]; i++) {
- voltages[g->nbi[i]] = t_voltages[i];
- }
- for (uint_t i = 0; i < 2; i++) {
- for (uint_t j = 0; j < g->bi[i + 1] - g->bi[i]; j++) {
- voltages[g->b[g->bi[i] + j]] = 1 - i;
- }
- }
-
- free(t_voltages);
- return voltages;
- } else {
- return t_voltages;
- }
-}
diff --git a/lib/rand.c b/lib/rand.c
deleted file mode 100644
index 8ba2b2b..0000000
--- a/lib/rand.c
+++ /dev/null
@@ -1,15 +0,0 @@
-
-#include "fracture.h"
-
-long double rand_dist_pow(const gsl_rng *r, double beta) {
- long double x = 0;
-
- // underflow means that for very small beta x is sometimes identically zero,
- // which causes problems
- while (x == 0.0) {
- long double y = logl(gsl_rng_uniform_pos(r)) / beta;
- x = expl(y);
- }
-
- return x;
-}
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);
+}
+