From 901b9f16494f37890be17ef4bb66e6efc6873340 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 10 Feb 2017 12:18:11 -0500 Subject: changed code to rely on jst --- lib/bound_set.c | 102 ++++++++++++++++++++++++++++++++++ lib/break_edge.c | 51 +++++++++++++++++ lib/correlations.c | 47 ++++++++++++++++ lib/data.c | 35 ++++++++++++ lib/factor_update.c | 68 +++++++++++++++++++++++ lib/fracture.h | 123 +++++++++++++++++++++++++++++++++++++++++ lib/gen_laplacian.c | 137 +++++++++++++++++++++++++++++++++++++++++++++ lib/gen_voltcurmat.c | 36 ++++++++++++ lib/geometry.c | 55 ++++++++++++++++++ lib/get_dual_clusters.c | 30 ++++++++++ lib/net.c | 144 ++++++++++++++++++++++++++++++++++++++++++++++++ lib/net_conductivity.c | 35 ++++++++++++ lib/net_currents.c | 51 +++++++++++++++++ lib/net_fracture.c | 67 ++++++++++++++++++++++ lib/net_voltages.c | 40 ++++++++++++++ lib/rand.c | 15 +++++ 16 files changed, 1036 insertions(+) create mode 100644 lib/bound_set.c create mode 100644 lib/break_edge.c create mode 100644 lib/correlations.c create mode 100644 lib/data.c create mode 100644 lib/factor_update.c create mode 100644 lib/fracture.h create mode 100644 lib/gen_laplacian.c create mode 100644 lib/gen_voltcurmat.c create mode 100644 lib/geometry.c create mode 100644 lib/get_dual_clusters.c create mode 100644 lib/net.c create mode 100644 lib/net_conductivity.c create mode 100644 lib/net_currents.c create mode 100644 lib/net_fracture.c create mode 100644 lib/net_voltages.c create mode 100644 lib/rand.c (limited to 'lib') diff --git a/lib/bound_set.c b/lib/bound_set.c new file mode 100644 index 0000000..65f1e6f --- /dev/null +++ b/lib/bound_set.c @@ -0,0 +1,102 @@ + +#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 new file mode 100644 index 0000000..2f112c2 --- /dev/null +++ b/lib/break_edge.c @@ -0,0 +1,51 @@ + +#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 new file mode 100644 index 0000000..63532ec --- /dev/null +++ b/lib/correlations.c @@ -0,0 +1,47 @@ + +#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 new file mode 100644 index 0000000..2047c44 --- /dev/null +++ b/lib/data.c @@ -0,0 +1,35 @@ + +#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 new file mode 100644 index 0000000..a51705a --- /dev/null +++ b/lib/factor_update.c @@ -0,0 +1,68 @@ + +#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 new file mode 100644 index 0000000..b1114fb --- /dev/null +++ b/lib/fracture.h @@ -0,0 +1,123 @@ + +#pragma once + +#include +#include +#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 new file mode 100644 index 0000000..1cc93a4 --- /dev/null +++ b/lib/gen_laplacian.c @@ -0,0 +1,137 @@ + +#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 new file mode 100644 index 0000000..e870140 --- /dev/null +++ b/lib/gen_voltcurmat.c @@ -0,0 +1,36 @@ + +#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 new file mode 100644 index 0000000..ec788f1 --- /dev/null +++ b/lib/geometry.c @@ -0,0 +1,55 @@ + +#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 new file mode 100644 index 0000000..eaf7562 --- /dev/null +++ b/lib/get_dual_clusters.c @@ -0,0 +1,30 @@ + +#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/net.c b/lib/net.c new file mode 100644 index 0000000..eb44721 --- /dev/null +++ b/lib/net.c @@ -0,0 +1,144 @@ + +#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); + + 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 new file mode 100644 index 0000000..e9325bb --- /dev/null +++ b/lib/net_conductivity.c @@ -0,0 +1,35 @@ + +#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 new file mode 100644 index 0000000..a078336 --- /dev/null +++ b/lib/net_currents.c @@ -0,0 +1,51 @@ + +#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 new file mode 100644 index 0000000..e7f18fc --- /dev/null +++ b/lib/net_fracture.c @@ -0,0 +1,67 @@ + +#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 new file mode 100644 index 0000000..7b07201 --- /dev/null +++ b/lib/net_voltages.c @@ -0,0 +1,40 @@ + +#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 new file mode 100644 index 0000000..16aeba4 --- /dev/null +++ b/lib/rand.c @@ -0,0 +1,15 @@ + +#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; +} -- cgit v1.2.3-54-g00ecf