summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/break_edge.c4
-rw-r--r--src/corr_test.c4
-rw-r--r--src/factor_update.c38
-rw-r--r--src/fracture.c4
-rw-r--r--src/fracture.h21
-rw-r--r--src/get_current.c98
-rw-r--r--src/net.c130
-rw-r--r--src/net_conductivity.c (renamed from src/get_conductivity.c)2
-rw-r--r--src/net_copy.c35
-rw-r--r--src/net_create.c69
-rw-r--r--src/net_currents.c35
-rw-r--r--src/net_fracture.c6
-rw-r--r--src/net_free.c14
-rw-r--r--src/net_notch.c29
-rw-r--r--src/net_voltages.c22
-rw-r--r--src/rand.c20
-rw-r--r--src/update_factor.c44
17 files changed, 266 insertions, 309 deletions
diff --git a/src/break_edge.c b/src/break_edge.c
index daed1c5..54eaf34 100644
--- a/src/break_edge.c
+++ b/src/break_edge.c
@@ -7,7 +7,7 @@ bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) {
unsigned int v1 = instance->graph->ev_break[2 * edge];
unsigned int v2 = instance->graph->ev_break[2 * edge + 1];
- if (instance->factor != NULL) update_factor(instance->factor, v1, v2, c);
+ if (instance->factor != NULL) factor_update(instance->factor, v1, v2, c);
if (instance->graph->boundary != TORUS_BOUND) {
unsigned int w1 = instance->graph->ev[2 * edge];
@@ -93,7 +93,7 @@ bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) {
lap_x[lap_p[tw2] + i] = 1;
}
- set_connected(instance->adjacency, instance->dual_marks, dw1, instance->dual_marks[v1], -1, 0);
+ set_connected(instance->adjacency, instance->dual_marks, dw1, instance->dual_marks[dw2], -1, 0);
}
diff --git a/src/corr_test.c b/src/corr_test.c
index 4805fb0..3d56c07 100644
--- a/src/corr_test.c
+++ b/src/corr_test.c
@@ -7,8 +7,8 @@ int main() {
unsigned int width = 64;
- graph_t *network = ini_square_network(width, true, false, &c);
- net_t *instance = net_create(network, 1e-14, 0.001, 0, true, &c);
+ graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 64, false, &c);
+ net_t *instance = net_create(network, 1e-14, 0.5, 0, true, &c);
net_fracture(instance, &c, 1e-10);
double *corr = get_corr(instance, NULL, &c);
diff --git a/src/factor_update.c b/src/factor_update.c
new file mode 100644
index 0000000..f27971b
--- /dev/null
+++ b/src/factor_update.c
@@ -0,0 +1,38 @@
+
+#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);
+}
diff --git a/src/fracture.c b/src/fracture.c
index f36cfdb..f38568f 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -320,7 +320,7 @@ int main(int argc, char *argv[]) {
}
if (save_stress_field || save_voltage_field) {
- double *tmp_voltages = get_voltage(net, &c);
+ double *tmp_voltages = net_voltages(net, &c);
if (save_voltage_field) {
for (uint_t j = 0; j < g->nv; j++) {
voltage_field[3 * voltage_pos] = g->vx[2 * j];
@@ -330,7 +330,7 @@ int main(int argc, char *argv[]) {
}
}
if (save_stress_field) {
- double *tmp_currents = get_current_v(net, tmp_voltages, &c);
+ double *tmp_currents = net_currents(net, tmp_voltages, &c);
for (uint_t j = 0; j < g->ne; j++) {
stress_field[3 * stress_pos] = g->ex[2 * j];
stress_field[3 * stress_pos + 1] = g->ex[2 * j + 1];
diff --git a/src/fracture.h b/src/fracture.h
index 2401d25..5e1b3e4 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -29,6 +29,8 @@
#define CINT_MAX INT_MAX
#define CHOL_F(x) cholmod_##x
+#define GSL_RAND_GEN gsl_rng_mt19937
+
typedef enum lattice_t {
VORONOI_LATTICE,
SQUARE_LATTICE
@@ -111,16 +113,13 @@ int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert,
bool index);
-bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2,
- cholmod_common *c);
-
-data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff);
-
-double *get_current(const net_t *instance, cholmod_common *c);
-double *get_current_v(const net_t *instance, double *voltages, cholmod_common *c);
-double *get_voltage(const net_t *instance, cholmod_common *c);
+void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, 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);
@@ -174,10 +173,12 @@ data_t *data_create(uint_t num_edges);
void data_free(data_t *data);
void data_update(data_t *data, uint_t last_broke, double strength, double conductivity);
-double get_conductivity(net_t *inst, double *current, cholmod_common *c);
-
graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c);
uint_t find_cycles(uint_t num_edges, const bool *fuses, const uint_t *ev, const uint_t *vei, const uint_t *ve, int **cycles);
bool set_connected(const cholmod_sparse *laplacian, uint_t *marks, int vertex, int label, int stop_at, int exclude);
+
+unsigned long int rand_seed();
+
+double rand_dist_pow(const gsl_rng *r, double beta);
diff --git a/src/get_current.c b/src/get_current.c
deleted file mode 100644
index a83c399..0000000
--- a/src/get_current.c
+++ /dev/null
@@ -1,98 +0,0 @@
-
-#include "fracture.h"
-
-double *get_voltage(const net_t *instance, cholmod_common *c) {
- cholmod_dense *b = instance->boundary_cond;
- cholmod_factor *factor = instance->factor;
-
- cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
-
- if (((double *)x->x)[0] != ((double *)x->x)[0]) {
- for (uint_t i = 0; i < instance->graph->ne; i++) {
- printf("%d ", instance->fuses[i]);
- }
- printf("\n");
- printf("GET_VOLTAGE: value is NaN\n");
- exit(EXIT_FAILURE);
- }
-
- double *field = (double *)x->x;
- x->x = NULL;
-
- CHOL_F(free_dense)(&x, c);
-
-
- return field;
-}
-
-double *get_current(const net_t *instance, cholmod_common *c) {
- unsigned int num_edges = instance->graph->ne;
- unsigned int num_gverts = instance->graph->break_dim;
- cholmod_sparse *voltcurmat = instance->graph->voltcurmat;
-
- double *voltages = get_voltage(instance, c);
- if (voltages == NULL) {
- return NULL;
- }
- cholmod_dense *x = CHOL_F(allocate_dense)(
- num_gverts, 1, num_gverts, CHOLMOD_REAL, c);
- double *tmp_x = x->x;
- x->x = voltages;
-
- cholmod_dense *y =
- CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c);
-
- double alpha[2] = {1, 0};
- double beta[2] = {0, 0};
- CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
-
- double *field = (double *)malloc(num_edges * sizeof(double));
-
- for (int i = 0; i < num_edges; i++) {
- if (instance->fuses[i])
- field[i] = 0;
- else
- field[i] = ((double *)y->x)[i];
- }
-
- x->x = tmp_x;
- free(voltages);
- CHOL_F(free_dense)(&x, c);
- CHOL_F(free_dense)(&y, c);
-
- return field;
-}
-
-
-double *get_current_v(const net_t *instance, double *voltages, cholmod_common *c) {
- unsigned int num_edges = instance->graph->ne;
- unsigned int num_gverts = instance->graph->break_dim;
- cholmod_sparse *voltcurmat = instance->graph->voltcurmat;
-
- cholmod_dense *x = CHOL_F(allocate_dense)(
- num_gverts, 1, num_gverts, CHOLMOD_REAL, c);
- double *tmp_x = x->x;
- x->x = voltages;
-
- cholmod_dense *y =
- CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c);
-
- double alpha[2] = {1, 0};
- double beta[2] = {0, 0};
- CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
-
- double *field = (double *)malloc(num_edges * sizeof(double));
-
- for (int i = 0; i < num_edges; i++) {
- if (instance->fuses[i])
- field[i] = 0;
- else
- field[i] = ((double *)y->x)[i];
- }
-
- x->x = tmp_x;
- CHOL_F(free_dense)(&x, c);
- CHOL_F(free_dense)(&y, c);
-
- return field;
-}
diff --git a/src/net.c b/src/net.c
new file mode 100644
index 0000000..9f6965b
--- /dev/null
+++ b/src/net.c
@@ -0,0 +1,130 @@
+
+#include "fracture.h"
+
+double *get_thres(uint_t ne, double beta) {
+ double *thres = (double *)malloc(ne * sizeof(double));
+ assert(thres != NULL);
+
+ gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN);
+ gsl_rng_set(r, 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->fuses = (bool *)calloc(g->ne, sizeof(bool));
+ assert(net->fuses != NULL);
+ net->thres = get_thres(g->ne, beta);
+ net->inf = inf;
+
+ net->voltage_bound = vb;
+ net->boundary_cond = bound_set(g, vb, notch_len, c);
+
+ if (g->boundary != TORUS_BOUND) net->adjacency = gen_adjacency(net, false, false, 0, c);
+ else net->adjacency = gen_adjacency(net, true, false, 0, c);
+
+ net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t));
+ assert(net->marks != NULL);
+
+ net->dual_marks = (uint_t *)malloc((net->graph->dnv) * sizeof(uint_t));
+ assert(net->dual_marks != NULL);
+
+ for (uint_t i = 0; i < (net->graph->break_dim); i++) {
+ net->marks[i] = 1;
+ }
+ for (uint_t i = 0; i < (net->graph->dnv); i++) {
+ net->dual_marks[i] = i+1;
+ }
+ net->num_components = 1;
+
+ net_notch(net, notch_len, c);
+
+ {
+ cholmod_sparse *laplacian = gen_laplacian(net, c, true);
+ net->factor = CHOL_F(analyze)(laplacian, c);
+ CHOL_F(factorize)(laplacian, net->factor, c);
+ CHOL_F(free_sparse)(&laplacian, 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(double);
+ net_copy->thres = (double *)malloc(thres_size);
+ assert(net_copy->thres != NULL);
+ memcpy(net_copy->thres, net->thres, thres_size);
+
+ size_t marks_size = (net->graph->break_dim) * sizeof(uint_t);
+ net_copy->marks = (uint_t *)malloc(marks_size);
+ assert(net_copy->marks != NULL);
+ memcpy(net_copy->marks, net->marks, marks_size);
+
+ size_t dual_marks_size = (net->graph->dnv) * sizeof(uint_t);
+ net_copy->dual_marks = (uint_t *)malloc(dual_marks_size);
+ assert(net_copy->dual_marks != NULL);
+ memcpy(net_copy->dual_marks, net->dual_marks, dual_marks_size);
+
+ net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c);
+ 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_sparse)(&(net->adjacency), c);
+ CHOL_F(free_factor)(&(net->factor), c);
+ free(net->marks);
+ free(net->dual_marks);
+ free(net);
+}
diff --git a/src/get_conductivity.c b/src/net_conductivity.c
index 23b7056..e9325bb 100644
--- a/src/get_conductivity.c
+++ b/src/net_conductivity.c
@@ -1,7 +1,7 @@
#include "fracture.h"
-double get_conductivity(net_t *net, double *voltages, cholmod_common *c) {
+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
diff --git a/src/net_copy.c b/src/net_copy.c
deleted file mode 100644
index dcb4080..0000000
--- a/src/net_copy.c
+++ /dev/null
@@ -1,35 +0,0 @@
-
-#include "fracture.h"
-
-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(double);
- net_copy->thres = (double *)malloc(thres_size);
- assert(net_copy->thres != NULL);
- memcpy(net_copy->thres, net->thres, thres_size);
-
- size_t marks_size = (net->graph->break_dim) * sizeof(uint_t);
- net_copy->marks = (uint_t *)malloc(marks_size);
- assert(net_copy->marks != NULL);
- memcpy(net_copy->marks, net->marks, marks_size);
-
- size_t dual_marks_size = (net->graph->dnv) * sizeof(uint_t);
- net_copy->dual_marks = (uint_t *)malloc(dual_marks_size);
- assert(net_copy->dual_marks != NULL);
- memcpy(net_copy->dual_marks, net->dual_marks, dual_marks_size);
-
- net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c);
- 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;
-}
-
diff --git a/src/net_create.c b/src/net_create.c
deleted file mode 100644
index fedcb3a..0000000
--- a/src/net_create.c
+++ /dev/null
@@ -1,69 +0,0 @@
-
-#include "fracture.h"
-
-double *get_thres(uint_t ne, double beta) {
- assert(beta > 0);
-
- double *thres = (double *)malloc(ne * sizeof(double));
- assert(thres != NULL);
-
- gsl_set_error_handler_off();
-
- gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
- {
- FILE *rf = fopen("/dev/urandom", "r");
- unsigned long int seed;
- fread(&seed, sizeof(unsigned long int), 1, rf);
- fclose(rf);
- gsl_rng_set(r, seed);
- }
-
- for (uint_t i = 0; i < ne; i++) {
- while ((thres[i] = exp(log(gsl_ran_flat(r, 0, 1)) / beta)) == 0.0);
- }
-
- gsl_rng_free(r);
-
- return thres;
-}
-
-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->fuses = (bool *)calloc(g->ne, sizeof(bool));
- assert(net->fuses != NULL);
- net->thres = get_thres(g->ne, beta);
- net->inf = inf;
-
- net->voltage_bound = vb;
- net->boundary_cond = bound_set(g, vb, notch_len, c);
-
- if (g->boundary != TORUS_BOUND) net->adjacency = gen_adjacency(net, false, false, 0, c);
- else net->adjacency = gen_adjacency(net, true, false, 0, c);
-
- net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t));
- net->dual_marks = (uint_t *)malloc((net->graph->dnv) * sizeof(uint_t));
- assert(net->marks != NULL);
-
- for (uint_t i = 0; i < (net->graph->break_dim); i++) {
- net->marks[i] = 1;
- }
- for (uint_t i = 0; i < (net->graph->dnv); i++) {
- net->dual_marks[i] = i+1;
- }
- net->num_components = 1;
-
- net_notch(net, notch_len, c);
-
- {
- cholmod_sparse *laplacian = gen_laplacian(net, c, true);
- net->factor = CHOL_F(analyze)(laplacian, c);
- CHOL_F(factorize)(laplacian, net->factor, c);
- CHOL_F(free_sparse)(&laplacian, c);
- }
-
- return net;
-}
-
diff --git a/src/net_currents.c b/src/net_currents.c
new file mode 100644
index 0000000..431818f
--- /dev/null
+++ b/src/net_currents.c
@@ -0,0 +1,35 @@
+
+#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->break_dim;
+ cholmod_sparse *voltcurmat = net->graph->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];
+ }
+ }
+
+ x->x = tmp_x;
+ CHOL_F(free_dense)(&x, c);
+ CHOL_F(free_dense)(&y, c);
+
+ return currents;
+}
diff --git a/src/net_fracture.c b/src/net_fracture.c
index 48e81f9..54f37dd 100644
--- a/src/net_fracture.c
+++ b/src/net_fracture.c
@@ -31,10 +31,10 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
data_t *data = data_create(net->graph->ne);
while (true) {
- double *voltages = get_voltage(net, c);
- double *currents = get_current_v(net, voltages, c);
+ double *voltages = net_voltages(net, c);
+ double *currents = net_currents(net, voltages, c);
- double conductivity = get_conductivity(net, voltages, c);
+ double conductivity = net_conductivity(net, voltages);
if (conductivity < 1e-12 && net->graph->boundary == TORUS_BOUND) {
free(voltages);
diff --git a/src/net_free.c b/src/net_free.c
deleted file mode 100644
index 8b5af50..0000000
--- a/src/net_free.c
+++ /dev/null
@@ -1,14 +0,0 @@
-
-#include "fracture.h"
-
-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_sparse)(&(net->adjacency), c);
- CHOL_F(free_factor)(&(net->factor), c);
- free(net->marks);
- free(net->dual_marks);
- free(net);
-}
-
diff --git a/src/net_notch.c b/src/net_notch.c
deleted file mode 100644
index ccbb387..0000000
--- a/src/net_notch.c
+++ /dev/null
@@ -1,29 +0,0 @@
-
-#include "fracture.h"
-
-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);
- }
- }
-}
diff --git a/src/net_voltages.c b/src/net_voltages.c
new file mode 100644
index 0000000..dedf5b2
--- /dev/null
+++ b/src/net_voltages.c
@@ -0,0 +1,22 @@
+
+#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 *voltages = (double *)x->x;
+ x->x = NULL;
+
+ CHOL_F(free_dense)(&x, c);
+
+ return voltages;
+}
+
diff --git a/src/rand.c b/src/rand.c
new file mode 100644
index 0000000..1a6a3d4
--- /dev/null
+++ b/src/rand.c
@@ -0,0 +1,20 @@
+
+#include "fracture.h"
+
+unsigned long int rand_seed() {
+ FILE *f = fopen("/dev/urandom", "r");
+ unsigned long int seed;
+ fread(&seed, sizeof(unsigned long int), 1, f);
+ fclose(f);
+ return seed;
+}
+
+double rand_dist_pow(const gsl_rng *r, double beta) {
+ double x = 0;
+
+ // underflow means that for very small beta x is sometimes identically zero,
+ // which causes problems
+ while ((x = exp(log(gsl_rng_uniform_pos(r)) / beta)) == 0.0);
+
+ return x;
+}
diff --git a/src/update_factor.c b/src/update_factor.c
deleted file mode 100644
index debfe7b..0000000
--- a/src/update_factor.c
+++ /dev/null
@@ -1,44 +0,0 @@
-
-#include "fracture.h"
-
-bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2,
- cholmod_common *c) {
- int n = factor->n;
- assert(v1 < n);
- assert(v2 < n);
-
- cholmod_sparse *update_mat =
- CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
-
- unsigned int v3, v4;
- v3 = v1 < v2 ? v1 : v2;
- v4 = v1 > v2 ? v1 : v2;
-
- for (int i = 0; i < n; i++) {
- if (i <= v3)
- ((int_t *)update_mat->p)[i] = 0;
- else
- ((int_t *)update_mat->p)[i] = 2;
- }
- ((int_t *)update_mat->p)[n] = 2;
- ((int_t *)update_mat->i)[0] = v3;
- ((int_t *)update_mat->i)[1] = v4;
- ((double *)update_mat->x)[0] = 1;
- ((double *)update_mat->x)[1] = -1;
-
- // assert(CHOL_F(check_sparse)(update_mat, c));
-
- cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
- update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
-
- // assert(CHOL_F(check_sparse)(perm_update_mat, c));
-
- CHOL_F(updown)(false, perm_update_mat, factor, c);
-
- // assert(CHOL_F(check_factor)(factor, c));
-
- CHOL_F(free_sparse)(&perm_update_mat, c);
- CHOL_F(free_sparse)(&update_mat, c);
-
- return true;
-}