summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2017-06-07 16:15:03 -0400
committerJaron Kent-Dobias <jaron@kent-dobias.com>2017-06-07 16:15:03 -0400
commitfa590b7542c11b417ee87eebdd162b7aa08b57d4 (patch)
treed4b884171ccd1072a78cfe0229b196523760b2c8
parent23682e445ba64452c453346f516f9101761438ba (diff)
downloadfuse_networks-fa590b7542c11b417ee87eebdd162b7aa08b57d4.tar.gz
fuse_networks-fa590b7542c11b417ee87eebdd162b7aa08b57d4.tar.bz2
fuse_networks-fa590b7542c11b417ee87eebdd162b7aa08b57d4.zip
ran clang-format on all files
-rw-r--r--lib/bound_set.c178
-rw-r--r--lib/break_edge.c91
-rw-r--r--lib/correlations.c75
-rw-r--r--lib/data.c40
-rw-r--r--lib/factor_update.c89
-rw-r--r--lib/fracture.h73
-rw-r--r--lib/gen_laplacian.c265
-rw-r--r--lib/gen_voltcurmat.c46
-rw-r--r--lib/geometry.c88
-rw-r--r--lib/get_dual_clusters.c35
-rw-r--r--lib/net.c251
-rw-r--r--lib/net_conductivity.c60
-rw-r--r--lib/net_currents.c95
-rw-r--r--lib/net_fracture.c89
-rw-r--r--lib/net_voltages.c69
-rw-r--r--lib/rand.c16
-rw-r--r--src/anal_process.c239
-rw-r--r--src/big_anal_process.c270
-rw-r--r--src/cen_anal_process.c269
-rw-r--r--src/corr_test.c30
-rw-r--r--src/fracture.c1017
-rw-r--r--src/fracture.h73
-rw-r--r--src/long_anal_process.c270
23 files changed, 1879 insertions, 1849 deletions
diff --git a/lib/bound_set.c b/lib/bound_set.c
index 65f1e6f..32d0fb7 100644
--- a/lib/bound_set.c
+++ b/lib/bound_set.c
@@ -2,101 +2,109 @@
#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;
-
+ 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)));
+ 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);
+ 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);
- }
+ 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;
+ 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;
+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
index 2f112c2..e53f7c1 100644
--- a/lib/break_edge.c
+++ b/lib/break_edge.c
@@ -1,51 +1,50 @@
#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;
+ 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
index 63532ec..98c5767 100644
--- a/lib/correlations.c
+++ b/lib/correlations.c
@@ -2,46 +2,45 @@
#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));
+ 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 < 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];
- }
- }
+ 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);
- }
+ if (nulldists) {
+ for (int i = 0; i < nv; i++) {
+ free(dists[i]);
+ }
+ free(dists);
+ }
- return corr;
+ return corr;
}
-
diff --git a/lib/data.c b/lib/data.c
index 2047c44..e39c450 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -2,34 +2,34 @@
#include "fracture.h"
data_t *data_create(uint_t ne) {
- data_t *data = malloc(1 * sizeof(data_t));
- assert(data != NULL);
+ data_t *data = malloc(1 * sizeof(data_t));
+ assert(data != NULL);
- data->num_broken = 0;
+ data->num_broken = 0;
- data->break_list = (uint_t *)malloc(ne * sizeof(uint_t));
- assert(data->break_list != NULL);
+ 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->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);
+ data->conductivity = (double *)malloc(ne * sizeof(double));
+ assert(data->conductivity != NULL);
- return data;
+ return data;
}
void data_free(data_t *data) {
- free(data->break_list);
- free(data->extern_field);
- free(data->conductivity);
- free(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++;
+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
index a51705a..ceaa01a 100644
--- a/lib/factor_update.c
+++ b/lib/factor_update.c
@@ -1,68 +1,69 @@
#include "fracture.h"
-void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c) {
- uint_t n = factor->n;
+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);
+ 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;
+ 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;
+ 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 = 0; i <= s1; i++) {
+ pp[i] = 0;
+ }
- for (uint_t i = s1 + 1; i <= n; i++) {
- pp[i] = 2;
- }
+ for (uint_t i = s1 + 1; i <= n; i++) {
+ pp[i] = 2;
+ }
- ii[0] = s1;
- ii[1] = s2;
- xx[0] = 1;
- xx[1] = -1;
+ 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);
+ 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(updown)(false, perm_update_mat, factor, c);
- CHOL_F(free_sparse)(&perm_update_mat, c);
- CHOL_F(free_sparse)(&update_mat, 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;
+ uint_t n = factor->n;
- cholmod_sparse *update_mat =
- CHOL_F(allocate_sparse)(n, n, 1, true, true, 0, CHOLMOD_REAL, c);
+ 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;
+ 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 = 0; i <= v; i++) {
+ pp[i] = 0;
+ }
- for (uint_t i = v + 1; i <= n; i++) {
- pp[i] = 1;
- }
+ for (uint_t i = v + 1; i <= n; i++) {
+ pp[i] = 1;
+ }
- ii[0] = v;
- xx[0] = 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);
+ 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(updown)(false, perm_update_mat, factor, c);
- CHOL_F(free_sparse)(&perm_update_mat, c);
- CHOL_F(free_sparse)(&update_mat, 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
index b1114fb..7b4f146 100644
--- a/lib/fracture.h
+++ b/lib/fracture.h
@@ -17,7 +17,6 @@
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
-#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
@@ -26,7 +25,6 @@
#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
@@ -36,68 +34,69 @@
#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;
+ 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;
+ 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);
+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_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 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);
+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);
+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_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_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);
+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);
+ 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);
+ 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);
+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);
@@ -112,11 +111,13 @@ 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);
+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);
+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);
diff --git a/lib/gen_laplacian.c b/lib/gen_laplacian.c
index 1cc93a4..75dccce 100644
--- a/lib/gen_laplacian.c
+++ b/lib/gen_laplacian.c
@@ -1,137 +1,142 @@
#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_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;
+ 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
index e870140..fede836 100644
--- a/lib/gen_voltcurmat.c
+++ b/lib/gen_voltcurmat.c
@@ -2,35 +2,35 @@
#include "fracture.h"
cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts,
- unsigned int *edges_to_verts,
- cholmod_common *c) {
+ 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);
+ 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;
+ 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;
- }
+ 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;
+ t_m->nnz = num_edges * 2;
- assert(CHOL_F(check_triplet)(t_m, c));
+ assert(CHOL_F(check_triplet)(t_m, c));
- cholmod_sparse *m = CHOL_F(triplet_to_sparse)(t_m, num_edges * 2, c);
+ cholmod_sparse *m = CHOL_F(triplet_to_sparse)(t_m, num_edges * 2, c);
- CHOL_F(free_triplet)(&t_m, c);
+ CHOL_F(free_triplet)(&t_m, c);
- return m;
+ return m;
}
diff --git a/lib/geometry.c b/lib/geometry.c
index ec788f1..2c1987b 100644
--- a/lib/geometry.c
+++ b/lib/geometry.c
@@ -2,54 +2,54 @@
#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;
- }
+ 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;
- }
+ 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);
- }
+ 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
index eaf7562..9336106 100644
--- a/lib/get_dual_clusters.c
+++ b/lib/get_dual_clusters.c
@@ -2,29 +2,30 @@
#include "fracture.h"
components_t *get_clusters(net_t *instance) {
- components_t *c = graph_components_get(instance->graph, instance->fuses, true);
- return c;
+ 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));
+ 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++;
- }
+ 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;
+ if (num_in_cluster == 0)
+ break;
- cluster_dist[num_in_cluster - 1]++;
- }
+ cluster_dist[num_in_cluster - 1]++;
+ }
- graph_components_free(c);
+ graph_components_free(c);
- return cluster_dist;
+ return cluster_dist;
}
diff --git a/lib/net.c b/lib/net.c
index d26b22c..ac89415 100644
--- a/lib/net.c
+++ b/lib/net.c
@@ -2,144 +2,145 @@
#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);
+ 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());
+ 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);
- }
+ for (uint_t i = 0; i < ne; i++) {
+ thres[i] = rand_dist_pow(r, beta);
+ }
- gsl_rng_free(r);
+ gsl_rng_free(r);
- return thres;
+ 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);
- }
- }
+ 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_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;
+ 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);
+ 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
index e9325bb..61148da 100644
--- a/lib/net_conductivity.c
+++ b/lib/net_conductivity.c
@@ -2,34 +2,34 @@
#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]);
- }
+ 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
index a078336..32dba04 100644
--- a/lib/net_currents.c
+++ b/lib/net_currents.c
@@ -1,51 +1,52 @@
#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;
+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
index e7f18fc..65ede9b 100644
--- a/lib/net_fracture.c
+++ b/lib/net_fracture.c
@@ -2,66 +2,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;
+ 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];
+ 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 (!broken && current > cutoff) {
+ long double val = current / net->thres[i];
- if (val > max_val) {
- max_val = val;
- max_pos = 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);
- }
+ 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;
+ return max_pos;
}
data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
- data_t *data = data_create(net->graph->ne);
+ 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);
+ 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);
+ double conductivity = net_conductivity(net, voltages);
- if (conductivity < cutoff) {
- free(voltages);
- free(currents);
- break;
- }
+ if (conductivity < cutoff) {
+ free(voltages);
+ free(currents);
+ break;
+ }
- uint_t last_broke = get_next_broken(net, currents, cutoff);
+ uint_t last_broke = get_next_broken(net, currents, cutoff);
- long double sim_current;
+ long double sim_current;
- if (net->voltage_bound) {
- sim_current = conductivity;
- } else {
- sim_current = 1;
- }
+ 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);
+ data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] /
+ ((long double)currents[last_broke])),
+ conductivity);
- free(voltages);
- free(currents);
+ free(voltages);
+ free(currents);
- break_edge(net, last_broke, c);
- }
+ break_edge(net, last_broke, c);
+ }
- return data;
+ return data;
}
-
diff --git a/lib/net_voltages.c b/lib/net_voltages.c
index 7b07201..d456a65 100644
--- a/lib/net_voltages.c
+++ b/lib/net_voltages.c
@@ -2,39 +2,38 @@
#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;
- }
+ 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
index 16aeba4..8ba2b2b 100644
--- a/lib/rand.c
+++ b/lib/rand.c
@@ -2,14 +2,14 @@
#include "fracture.h"
long double rand_dist_pow(const gsl_rng *r, double beta) {
- long double x = 0;
+ 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);
- }
+ // 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;
+ return x;
}
diff --git a/src/anal_process.c b/src/anal_process.c
index 0ebec84..de27571 100644
--- a/src/anal_process.c
+++ b/src/anal_process.c
@@ -1,134 +1,135 @@
#include "fracture.h"
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_laguerre.h>
-#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
-#include <gsl/gsl_blas.h>
void mom(uint_t len, uint_t n, uint32_t *data, double result[2]) {
- uint_t total = 0;
- double mom = 0;
- double mom_err = 0;
+ uint_t total = 0;
+ double mom = 0;
+ double mom_err = 0;
- for (uint_t i = 0; i < len; i++) {
- uint32_t datai = data[i];
- double in = pow(i, n);
+ for (uint_t i = 0; i < len; i++) {
+ uint32_t datai = data[i];
+ double in = pow(i, n);
- total += datai;
- mom += in * datai;
- mom_err += gsl_pow_2(in) * datai;
- }
+ total += datai;
+ mom += in * datai;
+ mom_err += gsl_pow_2(in) * datai;
+ }
- double momf = mom / total;
- double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);
+ double momf = mom / total;
+ double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);
- result[0] = momf;
- result[1] = momf_err;
+ result[0] = momf;
+ result[1] = momf_err;
}
int main(int argc, char *argv[]) {
- uint_t nc = argc - 1;
- uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
- double *betas_c = (double *)malloc(nc * sizeof(double));
- double *vals_c1 = (double *)malloc(nc * sizeof(double));
- double *errors_c1 = (double *)malloc(nc * sizeof(double));
- double *vals_c2 = (double *)malloc(nc * sizeof(double));
- double *errors_c2 = (double *)malloc(nc * sizeof(double));
- double *vals_c3 = (double *)malloc(nc * sizeof(double));
- double *errors_c3 = (double *)malloc(nc * sizeof(double));
- double *vals_c4 = (double *)malloc(nc * sizeof(double));
- double *errors_c4 = (double *)malloc(nc * sizeof(double));
-
- char *out_filename = (char *)malloc(100 * sizeof(char));
- uint_t out_filename_len = 0;
-
- for (uint_t i = 0; i < nc; i++) {
- char *fn = argv[1 + i];
- char *buff = (char *)malloc(20 * sizeof(char));
- uint_t pos = 0; uint_t c = 0;
- while (c < 5) {
- if (fn[pos] == '_') {
- c++;
- }
- if (i == 0) {
- out_filename[pos] = fn[pos];
- out_filename_len++;
- }
- pos++;
- }
- c = 0;
- while (fn[pos] != '_') {
- buff[c] = fn[pos];
- pos++;
- c++;
- }
- buff[c] = '\0';
- Ls_c[i] = atoi(buff);
- c = 0;
- pos++;
- while (fn[pos] != '_') {
- buff[c] = fn[pos];
- pos++;
- c++;
- }
- buff[c] = '\0';
- betas_c[i] = atof(buff);
- free(buff);
-
- uint_t dist_len = 4 * pow(Ls_c[i], 2);
- uint32_t *dist = malloc(dist_len * sizeof(uint32_t));
- FILE *dist_file = fopen(fn, "rb");
- fread(dist, sizeof(uint32_t), dist_len, dist_file);
- fclose(dist_file);
- {
- double mom1[2];
- mom(dist_len, 1, dist, mom1);
- vals_c1[i] = mom1[0];
- errors_c1[i] = mom1[1];
-
- double mom2[2];
- mom(dist_len, 2, dist, mom2);
- vals_c2[i] = mom2[0];
- errors_c2[i] = mom2[1];
-
- double mom3[2];
- mom(dist_len, 3, dist, mom3);
- vals_c3[i] = mom3[0];
- errors_c3[i] = mom3[1];
-
- double mom4[2];
- mom(dist_len, 4, dist, mom4);
- vals_c4[i] = mom4[0];
- errors_c4[i] = mom4[1];
- }
- free(dist);
- }
-
- out_filename[out_filename_len-1] = '.';
- out_filename[out_filename_len] = 't';
- out_filename[out_filename_len+1] = 'x';
- out_filename[out_filename_len+2] = 't';
- out_filename[out_filename_len+3] = '\0';
-
- FILE *out_file = fopen(out_filename, "w");
-
- for (uint_t i = 0; i < nc; i++) {
- fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
- }
-
- fclose(out_file);
-
-
- free(Ls_c);
- free(betas_c);
- free(vals_c1);
- free(errors_c1);
- free(vals_c2);
- free(errors_c2);
- free(vals_c3);
- free(errors_c3);
-
- return 0;
+ uint_t nc = argc - 1;
+ uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
+ double *betas_c = (double *)malloc(nc * sizeof(double));
+ double *vals_c1 = (double *)malloc(nc * sizeof(double));
+ double *errors_c1 = (double *)malloc(nc * sizeof(double));
+ double *vals_c2 = (double *)malloc(nc * sizeof(double));
+ double *errors_c2 = (double *)malloc(nc * sizeof(double));
+ double *vals_c3 = (double *)malloc(nc * sizeof(double));
+ double *errors_c3 = (double *)malloc(nc * sizeof(double));
+ double *vals_c4 = (double *)malloc(nc * sizeof(double));
+ double *errors_c4 = (double *)malloc(nc * sizeof(double));
+
+ char *out_filename = (char *)malloc(100 * sizeof(char));
+ uint_t out_filename_len = 0;
+
+ for (uint_t i = 0; i < nc; i++) {
+ char *fn = argv[1 + i];
+ char *buff = (char *)malloc(20 * sizeof(char));
+ uint_t pos = 0;
+ uint_t c = 0;
+ while (c < 5) {
+ if (fn[pos] == '_') {
+ c++;
+ }
+ if (i == 0) {
+ out_filename[pos] = fn[pos];
+ out_filename_len++;
+ }
+ pos++;
+ }
+ c = 0;
+ while (fn[pos] != '_') {
+ buff[c] = fn[pos];
+ pos++;
+ c++;
+ }
+ buff[c] = '\0';
+ Ls_c[i] = atoi(buff);
+ c = 0;
+ pos++;
+ while (fn[pos] != '_') {
+ buff[c] = fn[pos];
+ pos++;
+ c++;
+ }
+ buff[c] = '\0';
+ betas_c[i] = atof(buff);
+ free(buff);
+
+ uint_t dist_len = 4 * pow(Ls_c[i], 2);
+ uint32_t *dist = malloc(dist_len * sizeof(uint32_t));
+ FILE *dist_file = fopen(fn, "rb");
+ fread(dist, sizeof(uint32_t), dist_len, dist_file);
+ fclose(dist_file);
+ {
+ double mom1[2];
+ mom(dist_len, 1, dist, mom1);
+ vals_c1[i] = mom1[0];
+ errors_c1[i] = mom1[1];
+
+ double mom2[2];
+ mom(dist_len, 2, dist, mom2);
+ vals_c2[i] = mom2[0];
+ errors_c2[i] = mom2[1];
+
+ double mom3[2];
+ mom(dist_len, 3, dist, mom3);
+ vals_c3[i] = mom3[0];
+ errors_c3[i] = mom3[1];
+
+ double mom4[2];
+ mom(dist_len, 4, dist, mom4);
+ vals_c4[i] = mom4[0];
+ errors_c4[i] = mom4[1];
+ }
+ free(dist);
+ }
+
+ out_filename[out_filename_len - 1] = '.';
+ out_filename[out_filename_len] = 't';
+ out_filename[out_filename_len + 1] = 'x';
+ out_filename[out_filename_len + 2] = 't';
+ out_filename[out_filename_len + 3] = '\0';
+
+ FILE *out_file = fopen(out_filename, "w");
+
+ for (uint_t i = 0; i < nc; i++) {
+ fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i],
+ vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i],
+ errors_c3[i], vals_c4[i], errors_c4[i]);
+ }
+
+ fclose(out_file);
+
+ free(Ls_c);
+ free(betas_c);
+ free(vals_c1);
+ free(errors_c1);
+ free(vals_c2);
+ free(errors_c2);
+ free(vals_c3);
+ free(errors_c3);
+
+ return 0;
}
-
diff --git a/src/big_anal_process.c b/src/big_anal_process.c
index 0f7106f..8c1d8ba 100644
--- a/src/big_anal_process.c
+++ b/src/big_anal_process.c
@@ -1,156 +1,158 @@
#include "fracture.h"
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_laguerre.h>
-#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
-#include <gsl/gsl_blas.h>
#include <sys/stat.h>
void get_mean(uint_t len, double *data, double result[2]) {
- double total = 0;
+ double total = 0;
- for (uint_t i = 0; i < len; i++) {
- total += data[i];
- }
+ for (uint_t i = 0; i < len; i++) {
+ total += data[i];
+ }
- double mean = total / len;
- double total_sq_diff = 0;
+ double mean = total / len;
+ double total_sq_diff = 0;
- for (uint_t i = 0; i < len; i++) {
- total_sq_diff += pow(mean - data[i], 2);
- }
+ for (uint_t i = 0; i < len; i++) {
+ total_sq_diff += pow(mean - data[i], 2);
+ }
- double mean_err = sqrt(total_sq_diff) / len;
+ double mean_err = sqrt(total_sq_diff) / len;
- result[0] = mean;
- result[1] = mean_err;
+ result[0] = mean;
+ result[1] = mean_err;
}
-void get_mom(uint_t len, uint_t n, double *data, double mean[2], double result[2]) {
- double total_n_diff = 0;
- double total_n2_diff = 0;
+void get_mom(uint_t len, uint_t n, double *data, double mean[2],
+ double result[2]) {
+ double total_n_diff = 0;
+ double total_n2_diff = 0;
- for (uint_t i = 0; i < len; i++) {
- total_n_diff += pow(fabs(mean[0] - data[i]), n);
- total_n2_diff += pow(fabs(mean[0] - data[i]), 2 * n);
- }
+ for (uint_t i = 0; i < len; i++) {
+ total_n_diff += pow(fabs(mean[0] - data[i]), n);
+ total_n2_diff += pow(fabs(mean[0] - data[i]), 2 * n);
+ }
- double mom = total_n_diff / len;
- double mom_err = sqrt(total_n2_diff) / len;
+ double mom = total_n_diff / len;
+ double mom_err = sqrt(total_n2_diff) / len;
- result[0] = mom;
- result[1] = mom_err;
+ result[0] = mom;
+ result[1] = mom_err;
}
int main(int argc, char *argv[]) {
- uint_t nc = argc - 1;
- uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
- double *betas_c = (double *)malloc(nc * sizeof(double));
- double *vals_c1 = (double *)malloc(nc * sizeof(double));
- double *errors_c1 = (double *)malloc(nc * sizeof(double));
- double *vals_c2 = (double *)malloc(nc * sizeof(double));
- double *errors_c2 = (double *)malloc(nc * sizeof(double));
- double *vals_c3 = (double *)malloc(nc * sizeof(double));
- double *errors_c3 = (double *)malloc(nc * sizeof(double));
- double *vals_c4 = (double *)malloc(nc * sizeof(double));
- double *errors_c4 = (double *)malloc(nc * sizeof(double));
-
- char *out_filename = (char *)malloc(100 * sizeof(char));
- uint_t out_filename_len = 0;
-
- for (uint_t i = 0; i < nc; i++) {
- char *fn = argv[1 + i];
- char *buff = (char *)malloc(20 * sizeof(char));
- uint_t pos = 0; uint_t c = 0;
- while (c < 5) {
- if (fn[pos] == '_') {
- c++;
- }
- if (i == 0) {
- out_filename[pos] = fn[pos];
- out_filename_len++;
- }
- pos++;
- }
- c = 0;
- while (fn[pos] != '_') {
- buff[c] = fn[pos];
- pos++;
- c++;
- }
- buff[c] = '\0';
- Ls_c[i] = atoi(buff);
- c = 0;
- pos++;
- while (fn[pos] != '_') {
- buff[c] = fn[pos];
- pos++;
- c++;
- }
- buff[c] = '\0';
- betas_c[i] = atof(buff);
- free(buff);
-
- struct stat info;
- stat(fn, &info);
- uint_t num_bytes = info.st_size;
- uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(double);
-
- double *dist = malloc(dist_len * sizeof(double));
- FILE *dist_file = fopen(fn, "rb");
- fread(dist, sizeof(double), dist_len, dist_file);
- fclose(dist_file);
- {
- double mom1[2];
- get_mean(dist_len, dist, mom1);
- vals_c1[i] = mom1[0];
- errors_c1[i] = mom1[1];
-
- double mom2[2];
- get_mom(dist_len, 2, dist, mom1, mom2);
- vals_c2[i] = mom2[0];
- errors_c2[i] = mom2[1];
-
- double mom3[2];
- get_mom(dist_len, 3, dist, mom1, mom3);
- vals_c3[i] = mom3[0];
- errors_c3[i] = mom3[1];
-
- double mom4[2];
- get_mom(dist_len, 4, dist, mom1, mom4);
- vals_c4[i] = mom4[0];
- errors_c4[i] = mom4[1];
- }
- free(dist);
- }
-
- out_filename[out_filename_len-1] = '.';
- out_filename[out_filename_len] = 't';
- out_filename[out_filename_len+1] = 'x';
- out_filename[out_filename_len+2] = 't';
- out_filename[out_filename_len+3] = '\0';
-
- FILE *out_file = fopen(out_filename, "w");
-
- for (uint_t i = 0; i < nc; i++) {
- fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
- }
-
- fclose(out_file);
-
-
- free(Ls_c);
- free(betas_c);
- free(vals_c1);
- free(errors_c1);
- free(vals_c2);
- free(errors_c2);
- free(vals_c3);
- free(errors_c3);
- free(vals_c4);
- free(errors_c4);
-
- return 0;
+ uint_t nc = argc - 1;
+ uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
+ double *betas_c = (double *)malloc(nc * sizeof(double));
+ double *vals_c1 = (double *)malloc(nc * sizeof(double));
+ double *errors_c1 = (double *)malloc(nc * sizeof(double));
+ double *vals_c2 = (double *)malloc(nc * sizeof(double));
+ double *errors_c2 = (double *)malloc(nc * sizeof(double));
+ double *vals_c3 = (double *)malloc(nc * sizeof(double));
+ double *errors_c3 = (double *)malloc(nc * sizeof(double));
+ double *vals_c4 = (double *)malloc(nc * sizeof(double));
+ double *errors_c4 = (double *)malloc(nc * sizeof(double));
+
+ char *out_filename = (char *)malloc(100 * sizeof(char));
+ uint_t out_filename_len = 0;
+
+ for (uint_t i = 0; i < nc; i++) {
+ char *fn = argv[1 + i];
+ char *buff = (char *)malloc(20 * sizeof(char));
+ uint_t pos = 0;
+ uint_t c = 0;
+ while (c < 5) {
+ if (fn[pos] == '_') {
+ c++;
+ }
+ if (i == 0) {
+ out_filename[pos] = fn[pos];
+ out_filename_len++;
+ }
+ pos++;
+ }
+ c = 0;
+ while (fn[pos] != '_') {
+ buff[c] = fn[pos];
+ pos++;
+ c++;
+ }
+ buff[c] = '\0';
+ Ls_c[i] = atoi(buff);
+ c = 0;
+ pos++;
+ while (fn[pos] != '_') {
+ buff[c] = fn[pos];
+ pos++;
+ c++;
+ }
+ buff[c] = '\0';
+ betas_c[i] = atof(buff);
+ free(buff);
+
+ struct stat info;
+ stat(fn, &info);
+ uint_t num_bytes = info.st_size;
+ uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(double);
+
+ double *dist = malloc(dist_len * sizeof(double));
+ FILE *dist_file = fopen(fn, "rb");
+ fread(dist, sizeof(double), dist_len, dist_file);
+ fclose(dist_file);
+ {
+ double mom1[2];
+ get_mean(dist_len, dist, mom1);
+ vals_c1[i] = mom1[0];
+ errors_c1[i] = mom1[1];
+
+ double mom2[2];
+ get_mom(dist_len, 2, dist, mom1, mom2);
+ vals_c2[i] = mom2[0];
+ errors_c2[i] = mom2[1];
+
+ double mom3[2];
+ get_mom(dist_len, 3, dist, mom1, mom3);
+ vals_c3[i] = mom3[0];
+ errors_c3[i] = mom3[1];
+
+ double mom4[2];
+ get_mom(dist_len, 4, dist, mom1, mom4);
+ vals_c4[i] = mom4[0];
+ errors_c4[i] = mom4[1];
+ }
+ free(dist);
+ }
+
+ out_filename[out_filename_len - 1] = '.';
+ out_filename[out_filename_len] = 't';
+ out_filename[out_filename_len + 1] = 'x';
+ out_filename[out_filename_len + 2] = 't';
+ out_filename[out_filename_len + 3] = '\0';
+
+ FILE *out_file = fopen(out_filename, "w");
+
+ for (uint_t i = 0; i < nc; i++) {
+ fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i],
+ vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i],
+ errors_c3[i], vals_c4[i], errors_c4[i]);
+ }
+
+ fclose(out_file);
+
+ free(Ls_c);
+ free(betas_c);
+ free(vals_c1);
+ free(errors_c1);
+ free(vals_c2);
+ free(errors_c2);
+ free(vals_c3);
+ free(errors_c3);
+ free(vals_c4);
+ free(errors_c4);
+
+ return 0;
}
-
diff --git a/src/cen_anal_process.c b/src/cen_anal_process.c
index 3bf388c..ee2b029 100644
--- a/src/cen_anal_process.c
+++ b/src/cen_anal_process.c
@@ -1,154 +1,157 @@
#include "fracture.h"
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_laguerre.h>
-#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
-#include <gsl/gsl_blas.h>
void get_mean(uint_t len, uint32_t *data, double result[2]) {
- uint_t total = 0;
- double mean = 0;
- double mean_err = 0;
+ uint_t total = 0;
+ double mean = 0;
+ double mean_err = 0;
- for (uint_t i = 0; i < len; i++) {
- uint32_t datai = data[i];
+ for (uint_t i = 0; i < len; i++) {
+ uint32_t datai = data[i];
- total += datai;
- mean += i * datai;
- mean_err += gsl_pow_2(i) * datai;
- }
+ total += datai;
+ mean += i * datai;
+ mean_err += gsl_pow_2(i) * datai;
+ }
- double meanf = mean / total;
- double meanf_err = meanf * sqrt(mean_err / gsl_pow_2(mean) + 1 / total);
+ double meanf = mean / total;
+ double meanf_err = meanf * sqrt(mean_err / gsl_pow_2(mean) + 1 / total);
- result[0] = meanf;
- result[1] = meanf_err;
+ result[0] = meanf;
+ result[1] = meanf_err;
}
-void get_mom(uint_t len, uint_t n, uint32_t *data, double mean[2], double result[2]) {
- uint_t total = 0;
- double mom = 0;
- double mom_err = 0;
+void get_mom(uint_t len, uint_t n, uint32_t *data, double mean[2],
+ double result[2]) {
+ uint_t total = 0;
+ double mom = 0;
+ double mom_err = 0;
- for (uint_t i = 0; i < len; i++) {
- uint32_t datai = data[i];
- double in = pow(fabs(((double)i) - mean[0]), n);
+ for (uint_t i = 0; i < len; i++) {
+ uint32_t datai = data[i];
+ double in = pow(fabs(((double)i) - mean[0]), n);
- total += datai;
- mom += in * datai;
- mom_err += gsl_pow_2(in) * datai; // + gsl_pow_2(n * mean[1] / (((double)i) - mean[0])));
- }
+ total += datai;
+ mom += in * datai;
+ mom_err += gsl_pow_2(in) *
+ datai; // + gsl_pow_2(n * mean[1] / (((double)i) - mean[0])));
+ }
- double momf = mom / total;
- double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);
+ double momf = mom / total;
+ double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);
- result[0] = momf;
- result[1] = momf_err;
+ result[0] = momf;
+ result[1] = momf_err;
}
int main(int argc, char *argv[]) {
- uint_t nc = argc - 1;
- uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
- double *betas_c = (double *)malloc(nc * sizeof(double));
- double *vals_c1 = (double *)malloc(nc * sizeof(double));
- double *errors_c1 = (double *)malloc(nc * sizeof(double));
- double *vals_c2 = (double *)malloc(nc * sizeof(double));
- double *errors_c2 = (double *)malloc(nc * sizeof(double));
- double *vals_c3 = (double *)malloc(nc * sizeof(double));
- double *errors_c3 = (double *)malloc(nc * sizeof(double));
- double *vals_c4 = (double *)malloc(nc * sizeof(double));
- double *errors_c4 = (double *)malloc(nc * sizeof(double));
-
- char *out_filename = (char *)malloc(100 * sizeof(char));
- uint_t out_filename_len = 0;
-
- for (uint_t i = 0; i < nc; i++) {
- char *fn = argv[1 + i];
- char *buff = (char *)malloc(20 * sizeof(char));
- uint_t pos = 0; uint_t c = 0;
- while (c < 5) {
- if (fn[pos] == '_') {
- c++;
- }
- if (i == 0) {
- out_filename[pos] = fn[pos];
- out_filename_len++;
- }
- pos++;
- }
- c = 0;
- while (fn[pos] != '_') {
- buff[c] = fn[pos];
- pos++;
- c++;
- }
- buff[c] = '\0';
- Ls_c[i] = atoi(buff);
- c = 0;
- pos++;
- while (fn[pos] != '_') {
- buff[c] = fn[pos];
- pos++;
- c++;
- }
- buff[c] = '\0';
- betas_c[i] = atof(buff);
- free(buff);
-
- uint_t dist_len = 4 * pow(Ls_c[i], 2);
- uint32_t *dist = malloc(dist_len * sizeof(uint32_t));
- FILE *dist_file = fopen(fn, "rb");
- fread(dist, sizeof(uint32_t), dist_len, dist_file);
- fclose(dist_file);
- {
- double mom1[2];
- get_mean(dist_len, dist, mom1);
- vals_c1[i] = mom1[0];
- errors_c1[i] = mom1[1];
-
- double mom2[2];
- get_mom(dist_len, 2, dist, mom1, mom2);
- vals_c2[i] = mom2[0];
- errors_c2[i] = mom2[1];
-
- double mom3[2];
- get_mom(dist_len, 3, dist, mom1, mom3);
- vals_c3[i] = mom3[0];
- errors_c3[i] = mom3[1];
-
- double mom4[2];
- get_mom(dist_len, 4, dist, mom1, mom4);
- vals_c4[i] = mom4[0];
- errors_c4[i] = mom4[1];
- }
- free(dist);
- }
-
- out_filename[out_filename_len-1] = '.';
- out_filename[out_filename_len] = 't';
- out_filename[out_filename_len+1] = 'x';
- out_filename[out_filename_len+2] = 't';
- out_filename[out_filename_len+3] = '\0';
-
- FILE *out_file = fopen(out_filename, "w");
-
- for (uint_t i = 0; i < nc; i++) {
- fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
- }
-
- fclose(out_file);
-
-
- free(Ls_c);
- free(betas_c);
- free(vals_c1);
- free(errors_c1);
- free(vals_c2);
- free(errors_c2);
- free(vals_c3);
- free(errors_c3);
-
- return 0;
+ uint_t nc = argc - 1;
+ uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
+ double *betas_c = (double *)malloc(nc * sizeof(double));
+ double *vals_c1 = (double *)malloc(nc * sizeof(double));
+ double *errors_c1 = (double *)malloc(nc * sizeof(double));
+ double *vals_c2 = (double *)malloc(nc * sizeof(double));
+ double *errors_c2 = (double *)malloc(nc * sizeof(double));
+ double *vals_c3 = (double *)malloc(nc * sizeof(double));
+ double *errors_c3 = (double *)malloc(nc * sizeof(double));
+ double *vals_c4 = (double *)malloc(nc * sizeof(double));
+ double *errors_c4 = (double *)malloc(nc * sizeof(double));
+
+ char *out_filename = (char *)malloc(100 * sizeof(char));
+ uint_t out_filename_len = 0;
+
+ for (uint_t i = 0; i < nc; i++) {
+ char *fn = argv[1 + i];
+ char *buff = (char *)malloc(20 * sizeof(char));
+ uint_t pos = 0;
+ uint_t c = 0;
+ while (c < 5) {
+ if (fn[pos] == '_') {
+ c++;
+ }
+ if (i == 0) {
+ out_filename[pos] = fn[pos];
+ out_filename_len++;
+ }
+ pos++;
+ }
+ c = 0;
+ while (fn[pos] != '_') {
+ buff[c] = fn[pos];
+ pos++;
+ c++;
+ }
+ buff[c] = '\0';
+ Ls_c[i] = atoi(buff);
+ c = 0;
+ pos++;
+ while (fn[pos] != '_') {
+ buff[c] = fn[pos];
+ pos++;
+ c++;
+ }
+ buff[c] = '\0';
+ betas_c[i] = atof(buff);
+ free(buff);
+
+ uint_t dist_len = 4 * pow(Ls_c[i], 2);
+ uint32_t *dist = malloc(dist_len * sizeof(uint32_t));
+ FILE *dist_file = fopen(fn, "rb");
+ fread(dist, sizeof(uint32_t), dist_len, dist_file);
+ fclose(dist_file);
+ {
+ double mom1[2];
+ get_mean(dist_len, dist, mom1);
+ vals_c1[i] = mom1[0];
+ errors_c1[i] = mom1[1];
+
+ double mom2[2];
+ get_mom(dist_len, 2, dist, mom1, mom2);
+ vals_c2[i] = mom2[0];
+ errors_c2[i] = mom2[1];
+
+ double mom3[2];
+ get_mom(dist_len, 3, dist, mom1, mom3);
+ vals_c3[i] = mom3[0];
+ errors_c3[i] = mom3[1];
+
+ double mom4[2];
+ get_mom(dist_len, 4, dist, mom1, mom4);
+ vals_c4[i] = mom4[0];
+ errors_c4[i] = mom4[1];
+ }
+ free(dist);
+ }
+
+ out_filename[out_filename_len - 1] = '.';
+ out_filename[out_filename_len] = 't';
+ out_filename[out_filename_len + 1] = 'x';
+ out_filename[out_filename_len + 2] = 't';
+ out_filename[out_filename_len + 3] = '\0';
+
+ FILE *out_file = fopen(out_filename, "w");
+
+ for (uint_t i = 0; i < nc; i++) {
+ fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i],
+ vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i],
+ errors_c3[i], vals_c4[i], errors_c4[i]);
+ }
+
+ fclose(out_file);
+
+ free(Ls_c);
+ free(betas_c);
+ free(vals_c1);
+ free(errors_c1);
+ free(vals_c2);
+ free(errors_c2);
+ free(vals_c3);
+ free(errors_c3);
+
+ return 0;
}
-
diff --git a/src/corr_test.c b/src/corr_test.c
index 01b3e11..cb00212 100644
--- a/src/corr_test.c
+++ b/src/corr_test.c
@@ -2,26 +2,26 @@
#include "fracture.h"
int main() {
- cholmod_common c;
- CHOL_F(start)(&c);
+ cholmod_common c;
+ CHOL_F(start)(&c);
- unsigned int width = 64;
+ unsigned int width = 64;
- graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 128, false);
- net_t *instance = net_create(network, 1e-14, 3, 0, true, &c);
- net_fracture(instance, &c, 1e-10);
+ graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 128, false);
+ net_t *instance = net_create(network, 1e-14, 3, 0, true, &c);
+ net_fracture(instance, &c, 1e-10);
- double *corr = get_corr(instance, NULL, &c);
+ double *corr = get_corr(instance, NULL, &c);
- for (int i = 0; i < 2 * width; i++) {
- printf("%g ", corr[i]);
- }
- printf("\n");
+ for (int i = 0; i < 2 * width; i++) {
+ printf("%g ", corr[i]);
+ }
+ printf("\n");
- net_free(instance, &c);
- graph_free(network);
+ net_free(instance, &c);
+ graph_free(network);
- CHOL_F(finish)(&c);
+ CHOL_F(finish)(&c);
- return 0;
+ return 0;
}
diff --git a/src/fracture.c b/src/fracture.c
index 73059ff..ede7a24 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -3,510 +3,515 @@
#include "fracture.h"
int main(int argc, char *argv[]) {
- int opt;
-
- // defining variables to be (potentially) set by command line flags
- uint8_t filename_len;
- uint32_t N;
- uint_t L;
- double beta, inf, cutoff, crack_len;
- bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
- save_crit_stress, save_energy, save_conductivity,
- save_damage, save_threshold, save_current_load;
- bound_t boundary;
- lattice_t lattice;
-
-
- // assume filenames will be less than 100 characters
-
- filename_len = 100;
-
-
- // set default values
-
- N = 100;
- L = 16;
- crack_len = 0.;
- beta = .3;
- inf = 1e10;
- cutoff = 1e-9;
- boundary = FREE_BOUND;
- lattice = VORONOI_LATTICE;
- save_data = false;
- save_cluster_dist = false;
- use_voltage_boundaries = false;
- use_dual = false;
- save_network = false;
- save_crit_stress = false;
- save_damage = false;
- save_conductivity = false;
- save_energy = false;
- save_threshold = false;
- save_current_load = false;
-
-
- uint8_t bound_i;
- char boundc2 = 'f';
- uint8_t lattice_i;
- char lattice_c = 'v';
- char dual_c = 'o';
-
-
- // get commandline options
-
- while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TE")) != -1) {
- switch (opt) {
- case 'n':
- N = atoi(optarg);
- break;
- case 'L':
- L = atoi(optarg);
- break;
- case 'b':
- beta = atof(optarg);
- break;
- case 'l':
- crack_len = atof(optarg);
- break;
- case 'B':
- bound_i = atoi(optarg);
- switch (bound_i) {
- case 0:
- boundary = FREE_BOUND;
- boundc2 = 'f';
- break;
- case 1:
- boundary = CYLINDER_BOUND;
- boundc2 = 'c';
- break;
- case 2:
- boundary = TORUS_BOUND;
- use_voltage_boundaries = true;
- boundc2 = 't';
- break;
- case 3:
- boundary = EMBEDDED_BOUND;
- boundc2 = 'e';
- use_dual = true;
- use_voltage_boundaries = true;
- break;
- default:
- printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n");
- exit(EXIT_FAILURE);
- }
- break;
- case 'q':
- lattice_i = atoi(optarg);
- switch (lattice_i) {
- case 0:
- lattice = VORONOI_LATTICE;
- lattice_c = 'v';
- break;
- case 1:
- lattice = DIAGONAL_LATTICE;
- lattice_c = 'd';
- break;
- case 2:
- lattice = VORONOI_HYPERUNIFORM_LATTICE;
- lattice_c = 'h';
- break;
- case 3:
- lattice = TRIANGULAR_LATTICE;
- lattice_c = 't';
- break;
- case 4:
- lattice = SQUARE_LATTICE;
- lattice_c = 's';
- break;
- default:
- printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 (DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n");
- exit(EXIT_FAILURE);
- }
- break;
- case 'd':
- save_damage = true;
- break;
- case 'V':
- use_voltage_boundaries = true;
- break;
- case 'D':
- use_dual = true;
- dual_c = 'd';
- break;
- case 'c':
- save_cluster_dist = true;
- break;
- case 'o':
- save_data = true;
- break;
- case 'N':
- save_network = true;
- break;
- case 's':
- save_crit_stress = true;
- break;
- case 'r':
- save_conductivity = true;
- break;
- case 'E':
- save_energy = true;
- break;
- case 'T':
- save_threshold = true;
- break;
- case 'C':
- save_current_load = true;
- break;
- default: /* '?' */
- exit(EXIT_FAILURE);
- }
- }
-
-
- char boundc;
- if (use_voltage_boundaries) boundc = 'v';
- else boundc = 'c';
-
- FILE *data_out;
- if (save_data) {
- char *data_filename = (char *)malloc(filename_len * sizeof(char));
- snprintf(data_filename, filename_len, "data_%c_%c_%c_%c_%u_%g_%g.txt", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
- data_out = fopen(data_filename, "a");
- free(data_filename);
- }
-
- uint_t max_verts, max_edges;
-
- // these are very liberal estimates
- max_verts = 4 * pow(L, 2);
- max_edges = 4 * pow(L, 2);
-
- if (max_verts > CINT_MAX) {
- exit(EXIT_FAILURE);
- }
-
- // define arrays for saving cluster and avalanche distributions
- uint32_t *cluster_size_dist;
- uint32_t *avalanche_size_dist;
- char *c_filename;
- char *a_filename;
- if (save_cluster_dist) {
- cluster_size_dist =
- (uint32_t *)calloc(max_verts, sizeof(uint32_t));
- avalanche_size_dist =
- (uint32_t *)calloc(max_edges, sizeof(uint32_t));
-
- c_filename = (char *)malloc(filename_len * sizeof(char));
- a_filename = (char *)malloc(filename_len * sizeof(char));
- snprintf(c_filename, filename_len, "cstr_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
- snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-
- FILE *cluster_out = fopen(c_filename, "rb");
- FILE *avalanche_out = fopen(a_filename, "rb");
-
- if (cluster_out != NULL) {
- fread(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
- fclose(cluster_out);
- }
- if (avalanche_out != NULL) {
- fread(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);
- fclose(avalanche_out);
- }
- }
-
- long double *crit_stress;
- if (save_crit_stress) {
- crit_stress = (long double *)malloc(N * sizeof(long double));
- }
-
- double *conductivity;
- if (save_conductivity) {
- conductivity = (double *)malloc(N * sizeof(double));
- }
-
- // define arrays for saving damage distributions
- uint32_t *damage;
- char *d_filename;
- if (save_damage) {
- damage =
- (uint32_t *)calloc(max_edges, sizeof(uint32_t));
-
- d_filename = (char *)malloc(filename_len * sizeof(char));
- snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-
- FILE *damage_out = fopen(d_filename, "rb");
-
- if (damage_out != NULL) {
- fread(damage, sizeof(uint32_t), max_edges, damage_out);
- fclose(damage_out);
- }
- }
-
- long double *energy;
- if (save_energy) {
- energy = (long double *)malloc(N * sizeof(long double));
- }
-
- long double *thresholds;
- if (save_threshold) {
- thresholds = (long double *)malloc(N * sizeof(long double));
- }
-
- long double *loads;
- if (save_current_load) {
- loads = (long double *)malloc(N * sizeof(long double));
- }
-
-
- // start cholmod
- cholmod_common c;
- CHOL_F(start)(&c);
-
- /* if we use voltage boundary conditions, the laplacian matrix is positive
- * definite and we can use a supernodal LL decomposition. otherwise we need
- * to use the simplicial LDL decomposition
- */
- if (use_voltage_boundaries) {
- //(&c)->supernodal = CHOLMOD_SUPERNODAL;
- (&c)->supernodal = CHOLMOD_SIMPLICIAL;
- } else {
- (&c)->supernodal = CHOLMOD_SIMPLICIAL;
- }
-
-
- printf("\n");
- for (uint32_t i = 0; i < N; i++) {
- printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1, N);
-
- graph_t *g = graph_create(lattice, boundary, L, use_dual);
- net_t *net = net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c);
- net_t *tmp_net = net_copy(net, &c);
- data_t *data = net_fracture(tmp_net, &c, cutoff);
- net_free(tmp_net, &c);
-
- uint_t max_pos = 0;
- long double max_val = 0;
-
- double cond0;
- {
- double *tmp_voltages = net_voltages(net, &c);
- cond0 = net_conductivity(net, tmp_voltages);
- free(tmp_voltages);
- }
-
- for (uint_t j = 0; j < data->num_broken; j++) {
- long double val = data->extern_field[j];
-
- if (val > max_val) {
- max_pos = j;
- max_val = val;
- }
- }
-
- uint_t av_size = 0;
- long double cur_val = 0;
-
- for (uint_t j = 0; j < max_pos; j++) {
- uint_t next_broken = data->break_list[j];
-
- break_edge(net, next_broken, &c);
-
- long double val = data->extern_field[j];
- if (save_cluster_dist) {
- if (val < cur_val) {
- av_size++;
- }
-
- if (val > cur_val) {
- avalanche_size_dist[av_size]++;
- av_size = 0;
- cur_val = val;
- }
- }
- }
-
- if (save_crit_stress) crit_stress[i] = data->extern_field[max_pos];
-
- if (save_conductivity) {
- if (max_pos > 0) {
- conductivity[i] = data->conductivity[max_pos - 1];
- } else {
- conductivity[i] = cond0;
- }
- }
-
- if (save_damage) {
- uint_t would_break = 0;
- double *tmp_voltage = net_voltages(net, &c);
- double *tmp_current = net_currents(net, tmp_voltage, &c);
- free(tmp_voltage);
- for (uint_t j = 0; j < g->ne; j++) {
- bool broken = net->fuses[j];
- bool under_thres = net->thres[j] < net->thres[data->break_list[max_pos]];
- bool zero_field = fabs(tmp_current[j]) < cutoff;
- if (!broken && under_thres && zero_field) {
- break_edge(net, j, &c);
- }
- }
- damage[net->num_broken]++;
- free(tmp_current);
- }
-
- if (save_energy) {
- long double tmp_energy = 0;
- if (max_pos > 0) {
- long double sigma1 = data->extern_field[0];
- double cond1 = cond0;
- for (uint_t j = 0; j < max_pos - 1; j++) {
- long double sigma2 = data->extern_field[j+1];
- double cond2 = data->conductivity[j];
- if (sigma2 > sigma1) {
- tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1;
- sigma1 = sigma2;
- cond1 = cond2;
- }
- }
- }
- energy[i] = tmp_energy;
- }
-
- if (save_threshold) {
- thresholds[i] = net->thres[data->break_list[max_pos]];
- }
-
- if (save_current_load) {
- loads[i] = data->extern_field[max_pos] / net->thres[data->break_list[max_pos]];
- }
-
- if (save_data) {
- for (uint_t j = 0; j < data->num_broken; j++) {
- fprintf(data_out, "%u %Lg %g ", data->break_list[j],
- data->extern_field[j], data->conductivity[j]);
- }
- fprintf(data_out, "\n");
- }
-
- data_free(data);
- if (save_network) {
- FILE *net_out = fopen("network.txt", "w");
- for (uint_t j = 0; j < g->nv; j++) {
- fprintf(net_out, "%f %f ", g->vx[2 * j],
- g->vx[2 * j + 1]);
- }
- fprintf(net_out, "\n");
- for (uint_t j = 0; j < g->ne; j++) {
- fprintf(net_out, "%u %u ", g->ev[2 * j], g->ev[2 * j + 1]);
- }
- fprintf(net_out, "\n");
- for (uint_t j = 0; j < g->dnv; j++) {
- fprintf(net_out, "%f %f ", g->dvx[2 * j],
- g->dvx[2 * j + 1]);
- }
- fprintf(net_out, "\n");
- for (uint_t j = 0; j < g->ne; j++) {
- fprintf(net_out, "%u %u ", g->dev[2 * j], g->dev[2 * j + 1]);
- }
- fprintf(net_out, "\n");
- for (uint_t j = 0; j < g->ne; j++) {
- fprintf(net_out, "%d ", net->fuses[j]);
- }
- fclose(net_out);
- }
-
- if (save_cluster_dist) {
- uint_t *tmp_cluster_dist = get_cluster_dist(net);
- for (uint_t j = 0; j < g->dnv; j++) {
- cluster_size_dist[j] += tmp_cluster_dist[j];
- }
- free(tmp_cluster_dist);
- }
-
-
- net_free(net, &c);
- graph_free(g);
- }
-
- printf("\033[F\033[JFRACTURE: COMPLETE\n");
-
- if (save_cluster_dist) {
- FILE *cluster_out = fopen(c_filename, "wb");
- FILE *avalanche_out = fopen(a_filename, "wb");
-
- fwrite(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
- fwrite(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);
-
- fclose(cluster_out);
- fclose(avalanche_out);
-
- free(c_filename);
- free(a_filename);
- free(cluster_size_dist);
- free(avalanche_size_dist);
- }
-
- if (save_conductivity) {
- char *cond_filename = (char *)malloc(filename_len * sizeof(char));
- snprintf(cond_filename, filename_len, "cond_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
- FILE *cond_file = fopen(cond_filename, "ab");
- fwrite(conductivity, sizeof(double), N, cond_file);
- fclose(cond_file);
- free(cond_filename);
- free(conductivity);
- }
-
- if (save_energy) {
- char *tough_filename = (char *)malloc(filename_len * sizeof(char));
- snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
- FILE *tough_file = fopen(tough_filename, "ab");
- fwrite(energy, sizeof(long double), N, tough_file);
- fclose(tough_file);
- free(tough_filename);
- free(energy);
- }
-
- if (save_threshold) {
- char *thres_filename = (char *)malloc(filename_len * sizeof(char));
- snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
- FILE *thres_file = fopen(thres_filename, "ab");
- fwrite(thresholds, sizeof(long double), N, thres_file);
- fclose(thres_file);
- free(thres_filename);
- free(thresholds);
- }
-
- if (save_current_load) {
- char *load_filename = (char *)malloc(filename_len * sizeof(char));
- snprintf(load_filename, filename_len, "load_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
- FILE *load_file = fopen(load_filename, "ab");
- fwrite(loads, sizeof(long double), N, load_file);
- fclose(load_file);
- free(load_filename);
- free(loads);
- }
-
- if (save_damage) {
- FILE *hdam_file = fopen(d_filename, "wb");
- fwrite(damage, sizeof(uint32_t), max_edges, hdam_file);
- fclose(hdam_file);
- free(d_filename);
- free(damage);
- }
-
- if (save_data) {
- fclose(data_out);
- }
-
- if (save_crit_stress) {
- char *str_filename = (char *)malloc(filename_len * sizeof(char));
- snprintf(str_filename, filename_len, "strs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
- FILE *str_file = fopen(str_filename, "ab");
- fwrite(crit_stress, sizeof(long double), N, str_file);
- fclose(str_file);
- free(str_filename);
- free(crit_stress);
- }
-
- CHOL_F(finish)(&c);
-
- return 0;
+ int opt;
+
+ // defining variables to be (potentially) set by command line flags
+ uint8_t filename_len;
+ uint32_t N;
+ uint_t L;
+ double beta, inf, cutoff, crack_len;
+ bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual,
+ save_network, save_crit_stress, save_energy, save_conductivity,
+ save_damage, save_threshold, save_current_load;
+ bound_t boundary;
+ lattice_t lattice;
+
+ // assume filenames will be less than 100 characters
+
+ filename_len = 100;
+
+ // set default values
+
+ N = 100;
+ L = 16;
+ crack_len = 0.;
+ beta = .3;
+ inf = 1e10;
+ cutoff = 1e-9;
+ boundary = FREE_BOUND;
+ lattice = VORONOI_LATTICE;
+ save_data = false;
+ save_cluster_dist = false;
+ use_voltage_boundaries = false;
+ use_dual = false;
+ save_network = false;
+ save_crit_stress = false;
+ save_damage = false;
+ save_conductivity = false;
+ save_energy = false;
+ save_threshold = false;
+ save_current_load = false;
+
+ uint8_t bound_i;
+ char boundc2 = 'f';
+ uint8_t lattice_i;
+ char lattice_c = 'v';
+ char dual_c = 'o';
+
+ // get commandline options
+
+ while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TE")) != -1) {
+ switch (opt) {
+ case 'n':
+ N = atoi(optarg);
+ break;
+ case 'L':
+ L = atoi(optarg);
+ break;
+ case 'b':
+ beta = atof(optarg);
+ break;
+ case 'l':
+ crack_len = atof(optarg);
+ break;
+ case 'B':
+ bound_i = atoi(optarg);
+ switch (bound_i) {
+ case 0:
+ boundary = FREE_BOUND;
+ boundc2 = 'f';
+ break;
+ case 1:
+ boundary = CYLINDER_BOUND;
+ boundc2 = 'c';
+ break;
+ case 2:
+ boundary = TORUS_BOUND;
+ use_voltage_boundaries = true;
+ boundc2 = 't';
+ break;
+ case 3:
+ boundary = EMBEDDED_BOUND;
+ boundc2 = 'e';
+ use_dual = true;
+ use_voltage_boundaries = true;
+ break;
+ default:
+ printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), "
+ "or 2 (TORUS_BOUND).\n");
+ exit(EXIT_FAILURE);
+ }
+ break;
+ case 'q':
+ lattice_i = atoi(optarg);
+ switch (lattice_i) {
+ case 0:
+ lattice = VORONOI_LATTICE;
+ lattice_c = 'v';
+ break;
+ case 1:
+ lattice = DIAGONAL_LATTICE;
+ lattice_c = 'd';
+ break;
+ case 2:
+ lattice = VORONOI_HYPERUNIFORM_LATTICE;
+ lattice_c = 'h';
+ break;
+ case 3:
+ lattice = TRIANGULAR_LATTICE;
+ lattice_c = 't';
+ break;
+ case 4:
+ lattice = SQUARE_LATTICE;
+ lattice_c = 's';
+ break;
+ default:
+ printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 "
+ "(DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n");
+ exit(EXIT_FAILURE);
+ }
+ break;
+ case 'd':
+ save_damage = true;
+ break;
+ case 'V':
+ use_voltage_boundaries = true;
+ break;
+ case 'D':
+ use_dual = true;
+ dual_c = 'd';
+ break;
+ case 'c':
+ save_cluster_dist = true;
+ break;
+ case 'o':
+ save_data = true;
+ break;
+ case 'N':
+ save_network = true;
+ break;
+ case 's':
+ save_crit_stress = true;
+ break;
+ case 'r':
+ save_conductivity = true;
+ break;
+ case 'E':
+ save_energy = true;
+ break;
+ case 'T':
+ save_threshold = true;
+ break;
+ case 'C':
+ save_current_load = true;
+ break;
+ default: /* '?' */
+ exit(EXIT_FAILURE);
+ }
+ }
+
+ char boundc;
+ if (use_voltage_boundaries)
+ boundc = 'v';
+ else
+ boundc = 'c';
+
+ FILE *data_out;
+ if (save_data) {
+ char *data_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(data_filename, filename_len, "data_%c_%c_%c_%c_%u_%g_%g.txt",
+ lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+ data_out = fopen(data_filename, "a");
+ free(data_filename);
+ }
+
+ uint_t max_verts, max_edges;
+
+ // these are very liberal estimates
+ max_verts = 4 * pow(L, 2);
+ max_edges = 4 * pow(L, 2);
+
+ if (max_verts > CINT_MAX) {
+ exit(EXIT_FAILURE);
+ }
+
+ // define arrays for saving cluster and avalanche distributions
+ uint32_t *cluster_size_dist;
+ uint32_t *avalanche_size_dist;
+ char *c_filename;
+ char *a_filename;
+ if (save_cluster_dist) {
+ cluster_size_dist = (uint32_t *)calloc(max_verts, sizeof(uint32_t));
+ avalanche_size_dist = (uint32_t *)calloc(max_edges, sizeof(uint32_t));
+
+ c_filename = (char *)malloc(filename_len * sizeof(char));
+ a_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(c_filename, filename_len, "cstr_%c_%c_%c_%c_%d_%g_%g.dat",
+ lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+ snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%d_%g_%g.dat",
+ lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+
+ FILE *cluster_out = fopen(c_filename, "rb");
+ FILE *avalanche_out = fopen(a_filename, "rb");
+
+ if (cluster_out != NULL) {
+ fread(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
+ fclose(cluster_out);
+ }
+ if (avalanche_out != NULL) {
+ fread(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);
+ fclose(avalanche_out);
+ }
+ }
+
+ long double *crit_stress;
+ if (save_crit_stress) {
+ crit_stress = (long double *)malloc(N * sizeof(long double));
+ }
+
+ double *conductivity;
+ if (save_conductivity) {
+ conductivity = (double *)malloc(N * sizeof(double));
+ }
+
+ // define arrays for saving damage distributions
+ uint32_t *damage;
+ char *d_filename;
+ if (save_damage) {
+ damage = (uint32_t *)calloc(max_edges, sizeof(uint32_t));
+
+ d_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%d_%g_%g.dat",
+ lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+
+ FILE *damage_out = fopen(d_filename, "rb");
+
+ if (damage_out != NULL) {
+ fread(damage, sizeof(uint32_t), max_edges, damage_out);
+ fclose(damage_out);
+ }
+ }
+
+ long double *energy;
+ if (save_energy) {
+ energy = (long double *)malloc(N * sizeof(long double));
+ }
+
+ long double *thresholds;
+ if (save_threshold) {
+ thresholds = (long double *)malloc(N * sizeof(long double));
+ }
+
+ long double *loads;
+ if (save_current_load) {
+ loads = (long double *)malloc(N * sizeof(long double));
+ }
+
+ // start cholmod
+ cholmod_common c;
+ CHOL_F(start)(&c);
+
+ /* if we use voltage boundary conditions, the laplacian matrix is positive
+ * definite and we can use a supernodal LL decomposition. otherwise we need
+ * to use the simplicial LDL decomposition
+ */
+ if (use_voltage_boundaries) {
+ //(&c)->supernodal = CHOLMOD_SUPERNODAL;
+ (&c)->supernodal = CHOLMOD_SIMPLICIAL;
+ } else {
+ (&c)->supernodal = CHOLMOD_SIMPLICIAL;
+ }
+
+ printf("\n");
+ for (uint32_t i = 0; i < N; i++) {
+ printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1,
+ N);
+
+ graph_t *g = graph_create(lattice, boundary, L, use_dual);
+ net_t *net =
+ net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c);
+ net_t *tmp_net = net_copy(net, &c);
+ data_t *data = net_fracture(tmp_net, &c, cutoff);
+ net_free(tmp_net, &c);
+
+ uint_t max_pos = 0;
+ long double max_val = 0;
+
+ double cond0;
+ {
+ double *tmp_voltages = net_voltages(net, &c);
+ cond0 = net_conductivity(net, tmp_voltages);
+ free(tmp_voltages);
+ }
+
+ for (uint_t j = 0; j < data->num_broken; j++) {
+ long double val = data->extern_field[j];
+
+ if (val > max_val) {
+ max_pos = j;
+ max_val = val;
+ }
+ }
+
+ uint_t av_size = 0;
+ long double cur_val = 0;
+
+ for (uint_t j = 0; j < max_pos; j++) {
+ uint_t next_broken = data->break_list[j];
+
+ break_edge(net, next_broken, &c);
+
+ long double val = data->extern_field[j];
+ if (save_cluster_dist) {
+ if (val < cur_val) {
+ av_size++;
+ }
+
+ if (val > cur_val) {
+ avalanche_size_dist[av_size]++;
+ av_size = 0;
+ cur_val = val;
+ }
+ }
+ }
+
+ if (save_crit_stress)
+ crit_stress[i] = data->extern_field[max_pos];
+
+ if (save_conductivity) {
+ if (max_pos > 0) {
+ conductivity[i] = data->conductivity[max_pos - 1];
+ } else {
+ conductivity[i] = cond0;
+ }
+ }
+
+ if (save_damage) {
+ uint_t would_break = 0;
+ double *tmp_voltage = net_voltages(net, &c);
+ double *tmp_current = net_currents(net, tmp_voltage, &c);
+ free(tmp_voltage);
+ for (uint_t j = 0; j < g->ne; j++) {
+ bool broken = net->fuses[j];
+ bool under_thres =
+ net->thres[j] < net->thres[data->break_list[max_pos]];
+ bool zero_field = fabs(tmp_current[j]) < cutoff;
+ if (!broken && under_thres && zero_field) {
+ break_edge(net, j, &c);
+ }
+ }
+ damage[net->num_broken]++;
+ free(tmp_current);
+ }
+
+ if (save_energy) {
+ long double tmp_energy = 0;
+ if (max_pos > 0) {
+ long double sigma1 = data->extern_field[0];
+ double cond1 = cond0;
+ for (uint_t j = 0; j < max_pos - 1; j++) {
+ long double sigma2 = data->extern_field[j + 1];
+ double cond2 = data->conductivity[j];
+ if (sigma2 > sigma1) {
+ tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1;
+ sigma1 = sigma2;
+ cond1 = cond2;
+ }
+ }
+ }
+ energy[i] = tmp_energy;
+ }
+
+ if (save_threshold) {
+ thresholds[i] = net->thres[data->break_list[max_pos]];
+ }
+
+ if (save_current_load) {
+ loads[i] =
+ data->extern_field[max_pos] / net->thres[data->break_list[max_pos]];
+ }
+
+ if (save_data) {
+ for (uint_t j = 0; j < data->num_broken; j++) {
+ fprintf(data_out, "%u %Lg %g ", data->break_list[j],
+ data->extern_field[j], data->conductivity[j]);
+ }
+ fprintf(data_out, "\n");
+ }
+
+ data_free(data);
+ if (save_network) {
+ FILE *net_out = fopen("network.txt", "w");
+ for (uint_t j = 0; j < g->nv; j++) {
+ fprintf(net_out, "%f %f ", g->vx[2 * j], g->vx[2 * j + 1]);
+ }
+ fprintf(net_out, "\n");
+ for (uint_t j = 0; j < g->ne; j++) {
+ fprintf(net_out, "%u %u ", g->ev[2 * j], g->ev[2 * j + 1]);
+ }
+ fprintf(net_out, "\n");
+ for (uint_t j = 0; j < g->dnv; j++) {
+ fprintf(net_out, "%f %f ", g->dvx[2 * j], g->dvx[2 * j + 1]);
+ }
+ fprintf(net_out, "\n");
+ for (uint_t j = 0; j < g->ne; j++) {
+ fprintf(net_out, "%u %u ", g->dev[2 * j], g->dev[2 * j + 1]);
+ }
+ fprintf(net_out, "\n");
+ for (uint_t j = 0; j < g->ne; j++) {
+ fprintf(net_out, "%d ", net->fuses[j]);
+ }
+ fclose(net_out);
+ }
+
+ if (save_cluster_dist) {
+ uint_t *tmp_cluster_dist = get_cluster_dist(net);
+ for (uint_t j = 0; j < g->dnv; j++) {
+ cluster_size_dist[j] += tmp_cluster_dist[j];
+ }
+ free(tmp_cluster_dist);
+ }
+
+ net_free(net, &c);
+ graph_free(g);
+ }
+
+ printf("\033[F\033[JFRACTURE: COMPLETE\n");
+
+ if (save_cluster_dist) {
+ FILE *cluster_out = fopen(c_filename, "wb");
+ FILE *avalanche_out = fopen(a_filename, "wb");
+
+ fwrite(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
+ fwrite(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);
+
+ fclose(cluster_out);
+ fclose(avalanche_out);
+
+ free(c_filename);
+ free(a_filename);
+ free(cluster_size_dist);
+ free(avalanche_size_dist);
+ }
+
+ if (save_conductivity) {
+ char *cond_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(cond_filename, filename_len, "cond_%c_%c_%c_%c_%d_%g_%g.dat",
+ lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+ FILE *cond_file = fopen(cond_filename, "ab");
+ fwrite(conductivity, sizeof(double), N, cond_file);
+ fclose(cond_file);
+ free(cond_filename);
+ free(conductivity);
+ }
+
+ if (save_energy) {
+ char *tough_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat",
+ lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+ FILE *tough_file = fopen(tough_filename, "ab");
+ fwrite(energy, sizeof(long double), N, tough_file);
+ fclose(tough_file);
+ free(tough_filename);
+ free(energy);
+ }
+
+ if (save_threshold) {
+ char *thres_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat",
+ lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+ FILE *thres_file = fopen(thres_filename, "ab");
+ fwrite(thresholds, sizeof(long double), N, thres_file);
+ fclose(thres_file);
+ free(thres_filename);
+ free(thresholds);
+ }
+
+ if (save_current_load) {
+ char *load_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(load_filename, filename_len, "load_%c_%c_%c_%c_%d_%g_%g.dat",
+ lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+ FILE *load_file = fopen(load_filename, "ab");
+ fwrite(loads, sizeof(long double), N, load_file);
+ fclose(load_file);
+ free(load_filename);
+ free(loads);
+ }
+
+ if (save_damage) {
+ FILE *hdam_file = fopen(d_filename, "wb");
+ fwrite(damage, sizeof(uint32_t), max_edges, hdam_file);
+ fclose(hdam_file);
+ free(d_filename);
+ free(damage);
+ }
+
+ if (save_data) {
+ fclose(data_out);
+ }
+
+ if (save_crit_stress) {
+ char *str_filename = (char *)malloc(filename_len * sizeof(char));
+ snprintf(str_filename, filename_len, "strs_%c_%c_%c_%c_%d_%g_%g.dat",
+ lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+ FILE *str_file = fopen(str_filename, "ab");
+ fwrite(crit_stress, sizeof(long double), N, str_file);
+ fclose(str_file);
+ free(str_filename);
+ free(crit_stress);
+ }
+
+ CHOL_F(finish)(&c);
+
+ return 0;
}
diff --git a/src/fracture.h b/src/fracture.h
index b1114fb..7b4f146 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -17,7 +17,6 @@
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
-#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
@@ -26,7 +25,6 @@
#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
@@ -36,68 +34,69 @@
#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;
+ 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;
+ 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);
+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_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 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);
+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);
+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_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_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);
+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);
+ 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);
+ 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);
+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);
@@ -112,11 +111,13 @@ 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);
+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);
+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);
diff --git a/src/long_anal_process.c b/src/long_anal_process.c
index ba29152..d4ec4e0 100644
--- a/src/long_anal_process.c
+++ b/src/long_anal_process.c
@@ -1,156 +1,158 @@
#include "fracture.h"
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_laguerre.h>
-#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
-#include <gsl/gsl_blas.h>
#include <sys/stat.h>
void get_mean(uint_t len, long double *data, long double result[2]) {
- long double total = 0;
+ long double total = 0;
- for (uint_t i = 0; i < len; i++) {
- total += data[i];
- }
+ for (uint_t i = 0; i < len; i++) {
+ total += data[i];
+ }
- long double mean = total / len;
- long double total_sq_diff = 0;
+ long double mean = total / len;
+ long double total_sq_diff = 0;
- for (uint_t i = 0; i < len; i++) {
- total_sq_diff += pow(mean - data[i], 2);
- }
+ for (uint_t i = 0; i < len; i++) {
+ total_sq_diff += pow(mean - data[i], 2);
+ }
- long double mean_err = sqrt(total_sq_diff) / len;
+ long double mean_err = sqrt(total_sq_diff) / len;
- result[0] = mean;
- result[1] = mean_err;
+ result[0] = mean;
+ result[1] = mean_err;
}
-void get_mom(uint_t len, uint_t n, long double *data, long double mean[2], long double result[2]) {
- long double total_n_diff = 0;
- long double total_n2_diff = 0;
+void get_mom(uint_t len, uint_t n, long double *data, long double mean[2],
+ long double result[2]) {
+ long double total_n_diff = 0;
+ long double total_n2_diff = 0;
- for (uint_t i = 0; i < len; i++) {
- total_n_diff += pow(fabsl(mean[0] - data[i]), n);
- total_n2_diff += pow(fabsl(mean[0] - data[i]), 2 * n);
- }
+ for (uint_t i = 0; i < len; i++) {
+ total_n_diff += pow(fabsl(mean[0] - data[i]), n);
+ total_n2_diff += pow(fabsl(mean[0] - data[i]), 2 * n);
+ }
- long double mom = total_n_diff / len;
- long double mom_err = sqrt(total_n2_diff) / len;
+ long double mom = total_n_diff / len;
+ long double mom_err = sqrt(total_n2_diff) / len;
- result[0] = mom;
- result[1] = mom_err;
+ result[0] = mom;
+ result[1] = mom_err;
}
int main(int argc, char *argv[]) {
- uint_t nc = argc - 1;
- uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
- double *betas_c = (double *)malloc(nc * sizeof(double));
- long double *vals_c1 = (long double *)malloc(nc * sizeof(long double));
- long double *errors_c1 = (long double *)malloc(nc * sizeof(long double));
- long double *vals_c2 = (long double *)malloc(nc * sizeof(long double));
- long double *errors_c2 = (long double *)malloc(nc * sizeof(long double));
- long double *vals_c3 = (long double *)malloc(nc * sizeof(long double));
- long double *errors_c3 = (long double *)malloc(nc * sizeof(long double));
- long double *vals_c4 = (long double *)malloc(nc * sizeof(long double));
- long double *errors_c4 = (long double *)malloc(nc * sizeof(long double));
-
- char *out_filename = (char *)malloc(100 * sizeof(char));
- uint_t out_filename_len = 0;
-
- for (uint_t i = 0; i < nc; i++) {
- char *fn = argv[1 + i];
- char *buff = (char *)malloc(20 * sizeof(char));
- uint_t pos = 0; uint_t c = 0;
- while (c < 5) {
- if (fn[pos] == '_') {
- c++;
- }
- if (i == 0) {
- out_filename[pos] = fn[pos];
- out_filename_len++;
- }
- pos++;
- }
- c = 0;
- while (fn[pos] != '_') {
- buff[c] = fn[pos];
- pos++;
- c++;
- }
- buff[c] = '\0';
- Ls_c[i] = atoi(buff);
- c = 0;
- pos++;
- while (fn[pos] != '_') {
- buff[c] = fn[pos];
- pos++;
- c++;
- }
- buff[c] = '\0';
- betas_c[i] = atof(buff);
- free(buff);
-
- struct stat info;
- stat(fn, &info);
- uint_t num_bytes = info.st_size;
- uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(long double);
-
- long double *dist = malloc(dist_len * sizeof(long double));
- FILE *dist_file = fopen(fn, "rb");
- fread(dist, sizeof(long double), dist_len, dist_file);
- fclose(dist_file);
- {
- long double mom1[2];
- get_mean(dist_len, dist, mom1);
- vals_c1[i] = mom1[0];
- errors_c1[i] = mom1[1];
-
- long double mom2[2];
- get_mom(dist_len, 2, dist, mom1, mom2);
- vals_c2[i] = mom2[0];
- errors_c2[i] = mom2[1];
-
- long double mom3[2];
- get_mom(dist_len, 3, dist, mom1, mom3);
- vals_c3[i] = mom3[0];
- errors_c3[i] = mom3[1];
-
- long double mom4[2];
- get_mom(dist_len, 4, dist, mom1, mom4);
- vals_c4[i] = mom4[0];
- errors_c4[i] = mom4[1];
- }
- free(dist);
- }
-
- out_filename[out_filename_len-1] = '.';
- out_filename[out_filename_len] = 't';
- out_filename[out_filename_len+1] = 'x';
- out_filename[out_filename_len+2] = 't';
- out_filename[out_filename_len+3] = '\0';
-
- FILE *out_file = fopen(out_filename, "w");
-
- for (uint_t i = 0; i < nc; i++) {
- fprintf(out_file, "%u %g %Lg %Lg %Lg %Lg %Lg %Lg %Lg %Lg\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
- }
-
- fclose(out_file);
-
-
- free(Ls_c);
- free(betas_c);
- free(vals_c1);
- free(errors_c1);
- free(vals_c2);
- free(errors_c2);
- free(vals_c3);
- free(errors_c3);
- free(vals_c4);
- free(errors_c4);
-
- return 0;
+ uint_t nc = argc - 1;
+ uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
+ double *betas_c = (double *)malloc(nc * sizeof(double));
+ long double *vals_c1 = (long double *)malloc(nc * sizeof(long double));
+ long double *errors_c1 = (long double *)malloc(nc * sizeof(long double));
+ long double *vals_c2 = (long double *)malloc(nc * sizeof(long double));
+ long double *errors_c2 = (long double *)malloc(nc * sizeof(long double));
+ long double *vals_c3 = (long double *)malloc(nc * sizeof(long double));
+ long double *errors_c3 = (long double *)malloc(nc * sizeof(long double));
+ long double *vals_c4 = (long double *)malloc(nc * sizeof(long double));
+ long double *errors_c4 = (long double *)malloc(nc * sizeof(long double));
+
+ char *out_filename = (char *)malloc(100 * sizeof(char));
+ uint_t out_filename_len = 0;
+
+ for (uint_t i = 0; i < nc; i++) {
+ char *fn = argv[1 + i];
+ char *buff = (char *)malloc(20 * sizeof(char));
+ uint_t pos = 0;
+ uint_t c = 0;
+ while (c < 5) {
+ if (fn[pos] == '_') {
+ c++;
+ }
+ if (i == 0) {
+ out_filename[pos] = fn[pos];
+ out_filename_len++;
+ }
+ pos++;
+ }
+ c = 0;
+ while (fn[pos] != '_') {
+ buff[c] = fn[pos];
+ pos++;
+ c++;
+ }
+ buff[c] = '\0';
+ Ls_c[i] = atoi(buff);
+ c = 0;
+ pos++;
+ while (fn[pos] != '_') {
+ buff[c] = fn[pos];
+ pos++;
+ c++;
+ }
+ buff[c] = '\0';
+ betas_c[i] = atof(buff);
+ free(buff);
+
+ struct stat info;
+ stat(fn, &info);
+ uint_t num_bytes = info.st_size;
+ uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(long double);
+
+ long double *dist = malloc(dist_len * sizeof(long double));
+ FILE *dist_file = fopen(fn, "rb");
+ fread(dist, sizeof(long double), dist_len, dist_file);
+ fclose(dist_file);
+ {
+ long double mom1[2];
+ get_mean(dist_len, dist, mom1);
+ vals_c1[i] = mom1[0];
+ errors_c1[i] = mom1[1];
+
+ long double mom2[2];
+ get_mom(dist_len, 2, dist, mom1, mom2);
+ vals_c2[i] = mom2[0];
+ errors_c2[i] = mom2[1];
+
+ long double mom3[2];
+ get_mom(dist_len, 3, dist, mom1, mom3);
+ vals_c3[i] = mom3[0];
+ errors_c3[i] = mom3[1];
+
+ long double mom4[2];
+ get_mom(dist_len, 4, dist, mom1, mom4);
+ vals_c4[i] = mom4[0];
+ errors_c4[i] = mom4[1];
+ }
+ free(dist);
+ }
+
+ out_filename[out_filename_len - 1] = '.';
+ out_filename[out_filename_len] = 't';
+ out_filename[out_filename_len + 1] = 'x';
+ out_filename[out_filename_len + 2] = 't';
+ out_filename[out_filename_len + 3] = '\0';
+
+ FILE *out_file = fopen(out_filename, "w");
+
+ for (uint_t i = 0; i < nc; i++) {
+ fprintf(out_file, "%u %g %Lg %Lg %Lg %Lg %Lg %Lg %Lg %Lg\n", Ls_c[i],
+ betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i],
+ vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
+ }
+
+ fclose(out_file);
+
+ free(Ls_c);
+ free(betas_c);
+ free(vals_c1);
+ free(errors_c1);
+ free(vals_c2);
+ free(errors_c2);
+ free(vals_c3);
+ free(errors_c3);
+ free(vals_c4);
+ free(errors_c4);
+
+ return 0;
}
-