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/CMakeLists.txt | 4 + lib/bound_set.c | 110 ------------------- lib/break_edge.c | 50 --------- lib/correlations.c | 46 -------- lib/data.c | 35 ------ lib/factor_update.c | 69 ------------ lib/fracture.h | 123 --------------------- lib/gen_laplacian.c | 142 ------------------------ lib/gen_voltcurmat.c | 36 ------- lib/geometry.c | 55 ---------- lib/get_dual_clusters.c | 31 ------ lib/include/graph.hpp | 29 +++++ lib/include/hooks.hpp | 12 +++ lib/include/network.hpp | 48 +++++++++ lib/net.c | 146 ------------------------- lib/net_conductivity.c | 35 ------ lib/net_currents.c | 52 --------- lib/net_fracture.c | 68 ------------ lib/net_voltages.c | 39 ------- lib/rand.c | 15 --- lib/src/graph.cpp | 38 +++++++ lib/src/network.cpp | 280 ++++++++++++++++++++++++++++++++++++++++++++++++ 22 files changed, 411 insertions(+), 1052 deletions(-) create mode 100644 lib/CMakeLists.txt delete mode 100644 lib/bound_set.c delete mode 100644 lib/break_edge.c delete mode 100644 lib/correlations.c delete mode 100644 lib/data.c delete mode 100644 lib/factor_update.c delete mode 100644 lib/fracture.h delete mode 100644 lib/gen_laplacian.c delete mode 100644 lib/gen_voltcurmat.c delete mode 100644 lib/geometry.c delete mode 100644 lib/get_dual_clusters.c create mode 100644 lib/include/graph.hpp create mode 100644 lib/include/hooks.hpp create mode 100644 lib/include/network.hpp delete mode 100644 lib/net.c delete mode 100644 lib/net_conductivity.c delete mode 100644 lib/net_currents.c delete mode 100644 lib/net_fracture.c delete mode 100644 lib/net_voltages.c delete mode 100644 lib/rand.c create mode 100644 lib/src/graph.cpp create mode 100644 lib/src/network.cpp (limited to 'lib') 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 -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -// 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 +#include +#include + +class graph { + public: + typedef struct coordinate { + double x; + double y; + } coordinate; + + typedef struct vertex { + coordinate r; + } vertex; + + typedef std::array edge; + + std::vector vertices; + std::vector edges; + + std::vector dual_vertices; + std::vector 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>&, 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 +#include +#include +#include + +#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 fuses; + std::vector 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> 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 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-70-g09d2