diff options
| -rw-r--r-- | lib/bound_set.c | 178 | ||||
| -rw-r--r-- | lib/break_edge.c | 91 | ||||
| -rw-r--r-- | lib/correlations.c | 75 | ||||
| -rw-r--r-- | lib/data.c | 40 | ||||
| -rw-r--r-- | lib/factor_update.c | 89 | ||||
| -rw-r--r-- | lib/fracture.h | 73 | ||||
| -rw-r--r-- | lib/gen_laplacian.c | 265 | ||||
| -rw-r--r-- | lib/gen_voltcurmat.c | 46 | ||||
| -rw-r--r-- | lib/geometry.c | 88 | ||||
| -rw-r--r-- | lib/get_dual_clusters.c | 35 | ||||
| -rw-r--r-- | lib/net.c | 251 | ||||
| -rw-r--r-- | lib/net_conductivity.c | 60 | ||||
| -rw-r--r-- | lib/net_currents.c | 95 | ||||
| -rw-r--r-- | lib/net_fracture.c | 89 | ||||
| -rw-r--r-- | lib/net_voltages.c | 69 | ||||
| -rw-r--r-- | lib/rand.c | 16 | ||||
| -rw-r--r-- | src/anal_process.c | 239 | ||||
| -rw-r--r-- | src/big_anal_process.c | 270 | ||||
| -rw-r--r-- | src/cen_anal_process.c | 269 | ||||
| -rw-r--r-- | src/corr_test.c | 30 | ||||
| -rw-r--r-- | src/fracture.c | 1017 | ||||
| -rw-r--r-- | src/fracture.h | 73 | ||||
| -rw-r--r-- | src/long_anal_process.c | 270 |
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; } - @@ -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; } @@ -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; + } } - @@ -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; } - |