From 873a9f9bedbbfb07d475e271923a7b86464e515f Mon Sep 17 00:00:00 2001 From: pants Date: Wed, 7 Sep 2016 14:55:30 -0400 Subject: more major refactoring --- makefile | 4 +- makefile_hal | 4 +- src/bin_values.c | 4 +- src/bound_set.c | 62 +++++ src/break_edge.c | 9 +- src/compare_voronoi_fracture.c | 12 +- src/corr_test.c | 10 +- src/correlations.c | 2 +- src/course_grain_square.c | 2 +- src/coursegrain.c | 2 +- src/cracked_voronoi_fracture.c | 504 ----------------------------------------- src/cracking_ini.c | 77 ------- src/current_scaling.c | 6 +- src/fracture.c | 16 +- src/fracture.h | 44 ++-- src/fracture_network.c | 67 ------ src/free_network.c | 8 +- src/gen_laplacian.c | 15 +- src/get_conductivity.c | 4 +- src/homo_square_fracture.c | 10 +- src/homo_voronoi_fracture.c | 12 +- src/ini_network.c | 121 +++++----- src/instance.c | 124 ---------- src/net.c | 121 ++++++++++ src/net_fracture.c | 66 ++++++ src/net_notch.c | 29 +++ src/voro_fracture.c | 492 ++++++++++++++++++++++++++++++++++++++++ src/voronoi_bound_ini.c | 33 --- 28 files changed, 903 insertions(+), 957 deletions(-) create mode 100644 src/bound_set.c delete mode 100644 src/cracked_voronoi_fracture.c delete mode 100644 src/cracking_ini.c delete mode 100644 src/fracture_network.c delete mode 100644 src/instance.c create mode 100644 src/net.c create mode 100644 src/net_fracture.c create mode 100644 src/net_notch.c create mode 100644 src/voro_fracture.c delete mode 100644 src/voronoi_bound_ini.c diff --git a/makefile b/makefile index 0304fe6..14954e7 100644 --- a/makefile +++ b/makefile @@ -3,8 +3,8 @@ CC = clang CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -march=native #-fopenmp LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler -ltcmalloc -OBJ = break_data voronoi_bound_ini bin_values correlations beta_scales randfuncs instance get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry fracture_network get_current cracking_ini update_factor update_boundary get_file update_beta gen_voltcurmat ini_network free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity -BIN = current_scaling course_grain_square corr_test homo_voronoi_fracture homo_square_fracture compare_voronoi_fracture cracked_voronoi_fracture +OBJ = break_data bound_set bin_values correlations beta_scales randfuncs net get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry net_fracture get_current update_factor update_boundary get_file update_beta gen_voltcurmat ini_network free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity net_notch +BIN = corr_test voro_fracture all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} diff --git a/makefile_hal b/makefile_hal index 4559e77..f2bc954 100644 --- a/makefile_hal +++ b/makefile_hal @@ -3,8 +3,8 @@ CC = clang CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto #-fopenmp LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler -ltcmalloc -OBJ = break_data voronoi_bound_ini bin_values correlations beta_scales randfuncs instance get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry fracture_network get_current cracking_ini update_factor update_boundary get_file update_beta gen_voltcurmat ini_network free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity -BIN = current_scaling course_grain_square corr_test homo_voronoi_fracture homo_square_fracture compare_voronoi_fracture cracked_voronoi_fracture +OBJ = break_data bound_set bin_values correlations beta_scales randfuncs net get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry net_fracture get_current update_factor update_boundary get_file update_beta gen_voltcurmat ini_network free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity net_notch +BIN = corr_test voro_fracture all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} diff --git a/src/bin_values.c b/src/bin_values.c index 8cde548..70009c1 100644 --- a/src/bin_values.c +++ b/src/bin_values.c @@ -6,8 +6,8 @@ double *bin_values(graph_t *network, unsigned int width, double *values) { unsigned int *num_binned = calloc(pow(width, 2), sizeof(unsigned int)); for (unsigned int i = 0; i < network->ne; i++) { if (values[i] != 0) { - unsigned int x = ((int)(network->edge_coords[2 * i] * width)) % width; - unsigned int y = ((int)(network->edge_coords[2 * i + 1] * width)) % width; + unsigned int x = ((int)(network->ex[2 * i] * width)) % width; + unsigned int y = ((int)(network->ex[2 * i + 1] * width)) % width; binned[width * x + y] += fabs(values[i]); num_binned[width * x + y]++; } diff --git a/src/bound_set.c b/src/bound_set.c new file mode 100644 index 0000000..9c67f82 --- /dev/null +++ b/src/bound_set.c @@ -0,0 +1,62 @@ + +#include "fracture.h" + +double th_p(double x, double y, double th) { + if (x >= 0 && y >= 0) return th; + else if (x < 0 && y >= 0) return M_PI - th; + else if (x < 0 && y < 0) return th - M_PI; + else return -th; + +} + +double u_y(double x, double y) { + double r = sqrt(pow(x, 2) + pow(y, 2)); + double th = th_p(x, y, atan(fabs(y / x))); + + return sqrt(r) * sin(th / 2); +} + +void bound_set_embedded(double *bound, uint_t L, double notch_len) { + 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[i] = u_y(x1, y1); + bound[L / 2 + i] = u_y(x2, y2); + bound[L + i] = u_y(x3, y3); + bound[3 * L / 2 + i] = u_y(x4, y4); + } +} + +cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_common *c) { + cholmod_dense *boundary = CHOL_F(zeros)(g->break_dim, 1, CHOLMOD_REAL, c); + double *bound = (double *)boundary->x; + + switch (g->boundary) { + case TORUS_BOUND: + for (uint_t i = 0; i < g->bound_inds[1]; i++) { + bound[g->bound_verts[i]] = 1; + bound[g->nv + i] = -1; + } + bound[g->bound_verts[0]] = 1; + break; + case EMBEDDED_BOUND: + bound_set_embedded(bound, g->L, notch_len); + break; + default: + if (vb) { + for (uint_t i = 0; i < g->bound_inds[1]; i++) { + bound[g->bound_verts[i]] = 1; + } + } else { + bound[0] = 1; + bound[g->nv] = 1; + bound[g->nv + 1] = -1; + } + } + + return boundary; +} diff --git a/src/break_edge.c b/src/break_edge.c index 7b94160..4e1559b 100644 --- a/src/break_edge.c +++ b/src/break_edge.c @@ -3,7 +3,6 @@ bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) { instance->fuses[edge] = true; - instance->num_remaining_edges--; unsigned int v1 = instance->graph->ev_break[2 * edge]; unsigned int v2 = instance->graph->ev_break[2 * edge + 1]; @@ -56,10 +55,10 @@ bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) { double v1x, v1y, v2x, v2y; v1 = instance->graph->dev[2 * ee + !side]; v2 = instance->graph->dev[2 * ee + side]; - v1x = instance->graph->dual_vert_coords[2 * v1]; - v1y = instance->graph->dual_vert_coords[2 * v1 + 1]; - v2x = instance->graph->dual_vert_coords[2 * v2]; - v2y = instance->graph->dual_vert_coords[2 * v2 + 1]; + v1x = instance->graph->dvx[2 * v1]; + v1y = instance->graph->dvx[2 * v1 + 1]; + v2x = instance->graph->dvx[2 * v2]; + v2y = instance->graph->dvx[2 * v2 + 1]; double dx = v1x - v2x; double dy = v1y - v2y; if (((v1x > 0.5 && v2x < 0.5) || (v1x < 0.5 && v2x > 0.5)) && fabs(dx) < 0.5) { diff --git a/src/compare_voronoi_fracture.c b/src/compare_voronoi_fracture.c index 2edcae2..91fdcea 100644 --- a/src/compare_voronoi_fracture.c +++ b/src/compare_voronoi_fracture.c @@ -66,10 +66,10 @@ int main(int argc, char *argv[]) { (&c)->supernodal = CHOLMOD_SIMPLICIAL; - graph_t *network = ini_voronoi_network(L, false, boundary, genfunc_hyperuniform, &c); + graph_t *network = ini_voro_graph(L, false, boundary, genfunc_hyperuniform, &c); net_t *perm_voltage_instance = create_instance(network, inf, true, true, &c); net_t *perm_current_instance = create_instance(network, inf, false, true, &c); - double *fuse_thres = gen_fuse_thres(network->ne, network->edge_coords, beta, beta_scaling_flat); + double *fuse_thres = gen_fuse_thres(network->ne, network->ex, beta, beta_scaling_flat); net_t *voltage_instance = copy_instance(perm_voltage_instance, &c); net_t *current_instance = copy_instance(perm_current_instance, &c); data_t *breaking_data_voltage = fracture_network(voltage_instance, fuse_thres, &c, cutoff); @@ -82,8 +82,8 @@ int main(int argc, char *argv[]) { FILE *net_out = fopen("network.txt", "w"); for (unsigned int j = 0; j < network->nv; j++) { - fprintf(net_out, "%f %f ", network->vert_coords[2 * j], - network->vert_coords[2 * j + 1]); + fprintf(net_out, "%f %f ", network->vx[2 * j], + network->vx[2 * j + 1]); } fprintf(net_out, "\n"); for (unsigned int j = 0; j < network->ne; j++) { @@ -92,8 +92,8 @@ int main(int argc, char *argv[]) { } fprintf(net_out, "\n"); for (unsigned int j = 0; j < network->dnv; j++) { - fprintf(net_out, "%f %f ", network->dual_vert_coords[2 * j], - network->dual_vert_coords[2 * j + 1]); + fprintf(net_out, "%f %f ", network->dvx[2 * j], + network->dvx[2 * j + 1]); } fprintf(net_out, "\n"); for (unsigned int j = 0; j < network->ne; j++) { diff --git a/src/corr_test.c b/src/corr_test.c index e0de2d8..a8eaff5 100644 --- a/src/corr_test.c +++ b/src/corr_test.c @@ -9,10 +9,8 @@ int main() { unsigned int n = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2); graph_t *network = ini_square_network(width, true, false, &c); - net_t *instance = create_instance(network, 1e-14, true, true, &c); - double *fuse_thres = gen_fuse_thres(network->ne, network->edge_coords, - 0.001, beta_scaling_flat); - fracture_network(instance, fuse_thres, &c, 1e-10); + net_t *instance = net_create(network, 1e-14, 0.001, 0, true, &c); + net_fracture(instance, &c, 1e-10); double *corr = get_corr(instance, NULL, &c); @@ -21,8 +19,8 @@ int main() { } printf("\n"); - free_instance(instance, &c); - free_net(network, &c); + net_free(instance, &c); + graph_free(network, &c); CHOL_F(finish)(&c); diff --git a/src/correlations.c b/src/correlations.c index 4fc44ac..7fedfbf 100644 --- a/src/correlations.c +++ b/src/correlations.c @@ -298,7 +298,7 @@ double *get_corr(net_t *instance, unsigned int **dists, cholmod_common *c) { /* double *get_space_corr(net_t *instance, cholmod_common *c) { unsigned int nv = instance->graph->dnv; - double *vc = instance->graph->dual_vert_coords; + double *vc = instance->graph->dvx; unsigned int ne = instance->graph->ne; unsigned int *ev = instance->graph->dev; double *corr = calloc(nv, sizeof(unsigned int)); diff --git a/src/course_grain_square.c b/src/course_grain_square.c index 9c9ca6e..6587e83 100644 --- a/src/course_grain_square.c +++ b/src/course_grain_square.c @@ -84,7 +84,7 @@ int main(int argc, char *argv[]) { net_t *instance = NULL; while (breaking_data == NULL) { double *fuse_thres = gen_fuse_thres( - network->ne, network->edge_coords, beta, beta_scaling_flat); + network->ne, network->ex, beta, beta_scaling_flat); instance = copy_instance(perm_instance, &c); breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); free_instance(instance, &c); diff --git a/src/coursegrain.c b/src/coursegrain.c index 1e0d5a7..3e9db1a 100644 --- a/src/coursegrain.c +++ b/src/coursegrain.c @@ -5,7 +5,7 @@ net_t *coursegrain_square(net_t *instance, graph_t *network_p, cholmod_common *c unsigned int width = sqrt(instance->graph->ne); assert(width % 4 == 0); - net_t *instance_p = create_instance(network_p, instance->inf, + net_t *instance_p = create_instance(network_p, instance->inf, 1, instance->voltage_bound, true, c); unsigned int width_p = width / 2; diff --git a/src/cracked_voronoi_fracture.c b/src/cracked_voronoi_fracture.c deleted file mode 100644 index 93eb112..0000000 --- a/src/cracked_voronoi_fracture.c +++ /dev/null @@ -1,504 +0,0 @@ - - -#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 include_breaking, save_cluster_dist, use_voltage_boundaries, use_dual, save_network, - save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_conductivity, - save_damage, save_damage_field; - bound_t boundary; - - - // 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; - include_breaking = false; - save_cluster_dist = false; - use_voltage_boundaries = false; - use_dual = false; - save_network = false; - save_crit_stress = false; - save_stress_field = false; - save_voltage_field = false; - save_damage = false; - save_damage_field = false; - save_conductivity = false; - save_toughness = false; - - - uint8_t bound_i; - char boundc2 = 'f'; - - - // get commandline options - - while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrtDSvel:")) != -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 'd': - save_damage = true; - break; - case 'e': - save_damage_field = true; - break; - case 'V': - use_voltage_boundaries = true; - break; - case 'D': - use_dual = true; - break; - case 'c': - save_cluster_dist = true; - break; - case 'o': - include_breaking = true; - break; - case 'N': - save_network = true; - break; - case 's': - save_crit_stress = true; - break; - case 'S': - save_stress_field = true; - break; - case 'v': - save_voltage_field = true; - break; - case 'r': - save_conductivity = true; - //inf = 1; - break; - case 't': - save_toughness = true; - break; - default: /* '?' */ - exit(EXIT_FAILURE); - } - } - - - char boundc; - if (use_voltage_boundaries) boundc = 'v'; - else boundc = 'c'; - - FILE *break_out; - if (include_breaking) { - char *break_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(break_filename, filename_len, "breaks_v_%c_%c_%u_%g_%g.txt", boundc, boundc2, L, beta, crack_len); - break_out = fopen(break_filename, "a"); - free(break_filename); - } - - uint_t voronoi_max_verts, c_dist_size, a_dist_size; - - voronoi_max_verts = 4 * pow(L, 2); - c_dist_size = voronoi_max_verts; - a_dist_size = voronoi_max_verts; - - if (voronoi_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 *)malloc(c_dist_size * sizeof(uint32_t)); - avalanche_size_dist = - (uint32_t *)malloc(a_dist_size * 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_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); - snprintf(a_filename, filename_len, "avln_v_%c_%c_%d_%g_%g.dat", 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), c_dist_size, cluster_out); - fclose(cluster_out); - } - if (avalanche_out != NULL) { - fread(avalanche_size_dist, sizeof(uint32_t), a_dist_size, avalanche_out); - fclose(avalanche_out); - } - } - - double *crit_stress; - if (save_crit_stress) { - crit_stress = (double *)malloc(N * sizeof(double)); - } - - double *stress_field; - unsigned int stress_pos = 0; - if (save_stress_field) { - stress_field = (double *)malloc(3 * N * voronoi_max_verts * sizeof(double)); - } - - double *voltage_field; - unsigned int voltage_pos = 0; - if (save_voltage_field) { - voltage_field = (double *)malloc(3 * N * voronoi_max_verts * sizeof(double)); - } - - double *damage_field; - unsigned int damage_pos = 0; - if (save_damage_field) { - damage_field = (double *)malloc(2 * N * voronoi_max_verts * sizeof(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 *)malloc(a_dist_size * sizeof(uint32_t)); - - d_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(d_filename, filename_len, "damg_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); - - FILE *damage_out = fopen(d_filename, "rb"); - - if (damage_out != NULL) { - fread(damage, sizeof(uint32_t), a_dist_size, damage_out); - fclose(damage_out); - } - } - - double *toughness; - if (save_toughness) { - toughness = (double *)malloc(N * sizeof(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 *graph = ini_voronoi_network(L, boundary, use_dual, genfunc_hyperuniform, &c); - net_t *perm_instance = create_instance(graph, inf, use_voltage_boundaries, false, &c); - if (boundary == EMBEDDED_BOUND) { - voronoi_bound_ini(perm_instance, L, crack_len); - } - gen_voro_crack(perm_instance, crack_len, &c); - finish_instance(perm_instance, &c); - double *fuse_thres = gen_fuse_thres(graph->ne, graph->edge_coords, beta, beta_scaling_flat); - net_t *instance = copy_instance(perm_instance, &c); - data_t *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); - free_instance(instance, &c); - free(fuse_thres); - - uint_t max_pos = 0; - double max_val = 0; - - for (uint_t j = 0; j < breaking_data->num_broken; j++) { - double val = breaking_data->extern_field[j]; - - if (val > max_val) { - max_pos = j; - max_val = val; - } - } - - if (save_crit_stress) crit_stress[i] = breaking_data->extern_field[max_pos]; - - if (save_conductivity) conductivity[i] = breaking_data->conductivity[max_pos]; - - if (save_damage) damage[max_pos]++; - - net_t *tmp_instance = copy_instance(perm_instance, &c); - - uint_t av_size = 0; - double cur_val = 0; - for (uint_t j = 0; j < max_pos; j++) { - break_edge(tmp_instance, breaking_data->break_list[j], &c); - - double val = breaking_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_stress_field || save_voltage_field) { - double *tmp_voltages = get_voltage(tmp_instance, &c); - if (save_voltage_field) { - for (uint_t j = 0; j < tmp_instance->graph->nv; j++) { - voltage_field[3 * voltage_pos] = tmp_instance->graph->vert_coords[2 * j]; - voltage_field[3 * voltage_pos + 1] = tmp_instance->graph->vert_coords[2 * j + 1]; - voltage_field[3 * voltage_pos + 2] = tmp_voltages[j]; - voltage_pos++; - } - } - if (save_stress_field) { - double *tmp_currents = get_current_v(tmp_instance, tmp_voltages, &c); - for (uint_t j = 0; j < tmp_instance->graph->ne; j++) { - stress_field[3 * stress_pos] = tmp_instance->graph->edge_coords[2 * j]; - stress_field[3 * stress_pos + 1] = tmp_instance->graph->edge_coords[2 * j + 1]; - stress_field[3 * stress_pos + 2] = tmp_currents[j]; - stress_pos++; - } - free(tmp_currents); - } - free(tmp_voltages); - } - - if (save_damage_field) { - for (uint_t j = 0; j < max_pos; j++) { - damage_field[2 * damage_pos] = tmp_instance->graph->edge_coords[2 * breaking_data->break_list[j]]; - damage_field[2 * damage_pos + 1] = tmp_instance->graph->edge_coords[2 * breaking_data->break_list[j] + 1]; - damage_pos++; - } - } - - if (save_toughness) { - double tmp_toughness = 0; - if (max_pos > 0) { - double sigma1 = breaking_data->extern_field[0]; - double epsilon1 = sigma1 / breaking_data->conductivity[0]; - for (uint_t j = 0; j < max_pos - 1; j++) { - double sigma2 = breaking_data->extern_field[j+1]; - double epsilon2 = sigma2 / breaking_data->conductivity[j+1]; - if (epsilon2 > epsilon1) { - tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2; - sigma1 = sigma2; epsilon1 = epsilon2; - } - } - } - toughness[i] = tmp_toughness; - } - - if (save_cluster_dist) { - uint_t *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c); - for (uint_t j = 0; j < tmp_instance->graph->dnv; j++) { - cluster_size_dist[j] += tmp_cluster_dist[j]; - } - free(tmp_cluster_dist); - } - - if (save_network) { - FILE *net_out = fopen("network.txt", "w"); - for (uint_t j = 0; j < graph->nv; j++) { - fprintf(net_out, "%f %f ", graph->vert_coords[2 * j], - tmp_instance->graph->vert_coords[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < tmp_instance->graph->ne; j++) { - fprintf(net_out, "%u %u ", tmp_instance->graph->ev[2 * j], - tmp_instance->graph->ev[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < tmp_instance->graph->dnv; j++) { - fprintf(net_out, "%f %f ", tmp_instance->graph->dual_vert_coords[2 * j], - tmp_instance->graph->dual_vert_coords[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < tmp_instance->graph->ne; j++) { - fprintf(net_out, "%u %u ", tmp_instance->graph->dev[2 * j], - tmp_instance->graph->dev[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < tmp_instance->graph->ne; j++) { - fprintf(net_out, "%d ", tmp_instance->fuses[j]); - } - fclose(net_out); - } - - free_instance(tmp_instance, &c); - free_instance(perm_instance, &c); - free_net(graph, &c); - - if (include_breaking) { - for (uint_t j = 0; j < breaking_data->num_broken; j++) { - fprintf(break_out, "%u %g %g ", breaking_data->break_list[j], - breaking_data->extern_field[j], breaking_data->conductivity[j]); - } - fprintf(break_out, "\n"); - } - - free_break_data(breaking_data); - } - - 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), c_dist_size, cluster_out); - fwrite(avalanche_size_dist, sizeof(uint32_t), a_dist_size, avalanche_out); - - fclose(cluster_out); - fclose(avalanche_out); - - free(c_filename); - free(a_filename); - free(cluster_size_dist); - free(avalanche_size_dist); - } - - if (save_voltage_field) { - char *vfld_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(vfld_filename, filename_len, "vfld_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); - FILE *vfld_file = fopen(vfld_filename, "ab"); - fwrite(voltage_field, sizeof(double), 3 * voltage_pos, vfld_file); - fclose(vfld_file); - free(vfld_filename); - free(voltage_field); - } - - if (save_stress_field) { - char *cfld_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(cfld_filename, filename_len, "cfld_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); - FILE *cfld_file = fopen(cfld_filename, "ab"); - fwrite(stress_field, sizeof(double), 3 * stress_pos, cfld_file); - fclose(cfld_file); - free(cfld_filename); - free(stress_field); - } - - if (save_damage_field) { - char *dfld_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(dfld_filename, filename_len, "dfld_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); - FILE *dfld_file = fopen(dfld_filename, "ab"); - fwrite(damage_field, sizeof(double), 2 * damage_pos, dfld_file); - fclose(dfld_file); - free(dfld_filename); - free(damage_field); - } - - if (save_conductivity) { - char *cond_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(cond_filename, filename_len, "cond_v_%c_%c_%d_%g_%g.dat", 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_toughness) { - char *tough_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(tough_filename, filename_len, "tuff_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); - FILE *tough_file = fopen(tough_filename, "ab"); - fwrite(toughness, sizeof(double), N, tough_file); - fclose(tough_file); - free(tough_filename); - free(toughness); - } - - if (save_damage) { - FILE *hdam_file = fopen(d_filename, "wb"); - fwrite(damage, sizeof(uint32_t), a_dist_size, hdam_file); - fclose(hdam_file); - free(d_filename); - free(damage); - } - - if (include_breaking) { - fclose(break_out); - } - - if (save_crit_stress) { - char *str_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(str_filename, filename_len, "strs_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); - FILE *str_file = fopen(str_filename, "ab"); - fwrite(crit_stress, sizeof(double), N, str_file); - fclose(str_file); - free(str_filename); - free(crit_stress); - } - - CHOL_F(finish)(&c); - - return 0; -} diff --git a/src/cracking_ini.c b/src/cracking_ini.c deleted file mode 100644 index 93e5765..0000000 --- a/src/cracking_ini.c +++ /dev/null @@ -1,77 +0,0 @@ - -#include "fracture.h" - -double *gen_fuse_thres(unsigned int num_edges, double *edge_coords, double beta, - double (*beta_scale)(double, double, double)) { - unsigned int size = num_edges; - assert(beta > 0); - - double *fuse_thres = (double *)malloc(size * sizeof(double)); - assert(fuse_thres != NULL); - - gsl_set_error_handler_off(); - - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - FILE *rf = fopen("/dev/urandom", "r"); - unsigned long int seed; - fread(&seed, sizeof(unsigned long int), 1, rf); - fclose(rf); - gsl_rng_set(r, seed); - - for (unsigned int i = 0; i < size; i++) { - double x = edge_coords[2 * i]; - double y = edge_coords[2 * i + 1]; - while ((fuse_thres[i] = gsl_sf_exp(gsl_sf_log(gsl_ran_flat(r, 0, 1)) / - beta_scale(beta, x, y))) == 0.0) - ; - } - - gsl_rng_free(r); - - return fuse_thres; -} - -void gen_voro_crack(net_t *instance, double crack_len, cholmod_common *c) { - for (uint_t i = 0; i < instance->graph->ne; i++) { - uint_t v1, v2; - double v1x, v1y, v2x, v2y, dx, dy; - - v1 = instance->graph->ev[2 * i]; - v2 = instance->graph->ev[2 * i + 1]; - - v1x = instance->graph->vert_coords[2 * v1]; - v1y = instance->graph->vert_coords[2 * v1 + 1]; - v2x = instance->graph->vert_coords[2 * v2]; - v2y = instance->graph->vert_coords[2 * v2 + 1]; - - dx = v1x - v2x; - dy = v1y - v2y; - - if (((v1y > 0.5 && v2y < 0.5) || (v1y < 0.5 && v2y > 0.5)) && fabs(dy) < 0.5) { - if (v1x + dx / dy * (v1y - 0.5) <= crack_len) { - break_edge(instance, i, c); - } - } - } -} - -bool gen_crack(net_t *instance, double crack_len, double crack_width, - cholmod_common *c) { - assert(instance != NULL); - bool *fuses = instance->fuses; - assert(fuses != NULL); - const graph_t *network = instance->graph; - unsigned int num_edges = network->ne; - double *edge_coords = network->edge_coords; - - for (unsigned int j = 0; j < num_edges; j++) { - if (edge_coords[2 * j + 1] < crack_len && - fabs(edge_coords[2 * j] - network->L / 2) < crack_width) { - instance->fuses[j] = true; - instance->num_remaining_edges--; - } else - fuses[j] = false; - } - - return true; -} diff --git a/src/current_scaling.c b/src/current_scaling.c index f750cb1..ef36b0e 100644 --- a/src/current_scaling.c +++ b/src/current_scaling.c @@ -90,13 +90,13 @@ int main(int argc, char *argv[]) { graph_t *network = ini_square_network(width, false, true, &c); net_t *perm_instance = create_instance(network, inf, voltage_bound, false, &c); - gen_crack(perm_instance, crack_len, 1, &c); + gen_crack(perm_instance, crack_len, &c); finish_instance(perm_instance, &c); if (voltage_bound) { (&c)->supernodal = CHOLMOD_SIMPLICIAL; net_t *tmp_instance = create_instance(network, inf, false, false, &c); - gen_crack(tmp_instance, crack_len, 1, &c); + gen_crack(tmp_instance, crack_len, &c); finish_instance(tmp_instance, &c); double *voltage = get_voltage(tmp_instance, &c); @@ -151,7 +151,7 @@ int main(int argc, char *argv[]) { data_t *breaking_data = NULL; while (breaking_data == NULL) { double *fuse_thres = gen_fuse_thres( - network->ne, network->edge_coords, beta, beta_scaling_flat); + network->ne, network->ex, beta, beta_scaling_flat); net_t *instance = copy_instance(perm_instance, &c); breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); free_instance(instance, &c); diff --git a/src/fracture.c b/src/fracture.c index 23bce03..6053146 100644 --- a/src/fracture.c +++ b/src/fracture.c @@ -166,14 +166,14 @@ int main(int argc, char *argv[]) { genfunc_uniform, &c)) == NULL) ; perm_instance = create_instance(network, inf, voltage_bound, false, &c); - gen_crack(perm_instance, crack_len, crack_width, &c); + gen_crack(perm_instance, crack_len, &c); finish_instance(perm_instance, &c); } } else { network = ini_square_network(width, periodic, side_bounds, &c); crack_width = 1; perm_instance = create_instance(network, inf, voltage_bound, false, &c); - gen_crack(perm_instance, crack_len, crack_width, &c); + gen_crack(perm_instance, crack_len, &c); finish_instance(perm_instance, &c); c_dist_size = network->dnv; a_dist_size = network->nv; @@ -286,14 +286,14 @@ int main(int argc, char *argv[]) { genfunc_uniform, &c)) == NULL) ; perm_instance = create_instance(network, inf, voltage_bound, false, &c); - gen_crack(perm_instance, crack_len, crack_width, &c); + gen_crack(perm_instance, crack_len, &c); finish_instance(perm_instance, &c); if (supplied_bound) { voronoi_bound_ini(perm_instance, square_bound, width); } } double *fuse_thres = gen_fuse_thres( - network->ne, network->edge_coords, beta, beta_scaling_flat); + network->ne, network->ex, beta, beta_scaling_flat); finst *instance = copy_instance(perm_instance, &c); breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); free_instance(instance, &c); @@ -414,8 +414,8 @@ int main(int argc, char *argv[]) { if (network_out) { FILE *net_out = fopen("network.txt", "w"); for (unsigned int i = 0; i < network->nv; i++) { - fprintf(net_out, "%f %f ", network->vert_coords[2 * i], - network->vert_coords[2 * i + 1]); + fprintf(net_out, "%f %f ", network->vx[2 * i], + network->vx[2 * i + 1]); } fprintf(net_out, "\n"); for (unsigned int i = 0; i < network->ne; i++) { @@ -424,8 +424,8 @@ int main(int argc, char *argv[]) { } fprintf(net_out, "\n"); for (unsigned int i = 0; i < network->dnv; i++) { - fprintf(net_out, "%f %f ", network->dual_vert_coords[2 * i], - network->dual_vert_coords[2 * i + 1]); + fprintf(net_out, "%f %f ", network->dvx[2 * i], + network->dvx[2 * i + 1]); } fprintf(net_out, "\n"); for (unsigned int i = 0; i < network->ne; i++) { diff --git a/src/fracture.h b/src/fracture.h index a9973ea..17ab8fd 100644 --- a/src/fracture.h +++ b/src/fracture.h @@ -37,38 +37,35 @@ typedef enum bound_t { } bound_t; typedef struct { - uint_t ne; // number of edges - uint_t nv; // number of vertices - uint_t dnv; // number of dual vertices + uint_t ne; + uint_t nv; + uint_t dnv; uint_t nv_break; uint_t num_bounds; - // 2*ne length array. edge i connects vertices ev[2*i], ev[2*i+1] < nv uint_t *ev; uint_t *ev_break; - // nv+1 length array. vertex i has degree vei[i+1]-vei[i], vei[0] = 0 uint_t *vei; - // vei[nv] length array. vertex i is connected to edges ve[vei[i]+j], 0 <= j < vei[i+1]-vei[i] uint_t *ve; uint_t *bound_inds; uint_t *bound_verts; - double *vert_coords; - double *edge_coords; + double *vx; + double *ex; uint_t *dev; uint_t *dvei; uint_t *dve; - double *dual_vert_coords; + double *dvx; uint_t num_spanning_edges; uint_t *spanning_edges; - double L; + uint_t L; uint_t break_dim; cholmod_sparse *voltcurmat; bound_t boundary; } graph_t; typedef struct { - graph_t *graph; - unsigned int num_remaining_edges; + const graph_t *graph; bool *fuses; + double *thres; double inf; cholmod_dense *boundary_cond; cholmod_factor *factor; @@ -112,18 +109,13 @@ double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert, bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2, cholmod_common *c); -data_t *fracture_network(net_t *instance, double *fuse_thres, - cholmod_common *c, double cutoff); +data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff); double *get_current(const net_t *instance, cholmod_common *c); double *get_current_v(const net_t *instance, double *voltages, cholmod_common *c); double *get_voltage(const net_t *instance, cholmod_common *c); -double *gen_fuse_thres(unsigned int num_edges, double *edge_coords, double beta, - double (*beta_scale)(double, double, double)); - -bool gen_crack(net_t *instance, double crack_len, double crack_width, - cholmod_common *c); +void net_notch(net_t *net, double notch_len, cholmod_common *c); void update_boundary(net_t *instance, const double *avg_field); @@ -137,18 +129,17 @@ double update_beta(double beta, unsigned int width, const double *stress, cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts, unsigned int *edges_to_verts, cholmod_common *c); -net_t *copy_instance(const net_t *instance, cholmod_common *c); +net_t *net_copy(const net_t *net, cholmod_common *c); graph_t *ini_square_network(unsigned int width, bound_t boundary, bool side_bounds, cholmod_common *c); -void free_net(graph_t *network, cholmod_common *c); -void free_instance(net_t *instance, cholmod_common *c); +void graph_free(graph_t *network, cholmod_common *c); +void net_free(net_t *instance, cholmod_common *c); -net_t *create_instance(graph_t *network, double inf, bool voltage_bound, - bool startnow, cholmod_common *c); +net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c); -graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual, +graph_t *ini_voro_graph(unsigned int L, bound_t boundary, bool use_dual, double *(*genfunc)(unsigned int, bound_t, gsl_rng *, unsigned int *), cholmod_common *c); @@ -178,7 +169,7 @@ double *get_corr(net_t *instance, unsigned int **dists, cholmod_common *c); double *bin_values(graph_t *network, unsigned int width, double *values); -void voronoi_bound_ini(net_t *instance, uint_t L, double crack_len); +cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_common *c); data_t *alloc_break_data(unsigned int num_edges); void free_break_data(data_t *data); @@ -186,4 +177,3 @@ void update_break_data(data_t *data, unsigned int last_broke, double strength, d double get_conductivity(net_t *inst, double *current, cholmod_common *c); -void gen_voro_crack(net_t *instance, double crack_len, cholmod_common *c); diff --git a/src/fracture_network.c b/src/fracture_network.c deleted file mode 100644 index 428b092..0000000 --- a/src/fracture_network.c +++ /dev/null @@ -1,67 +0,0 @@ - -#include "fracture.h" - -int inc_break_fuses(net_t *instance, double *thres, double *field, - double cutoff) { - unsigned int size = (instance->graph)->ne; - - int min_pos = -1; - long double min_val = -1; - - for (unsigned int i = 0; i < size; i++) { - if (!instance->fuses[i] && fabs(field[i]) > cutoff) { - double val = fabs(field[i]) / thres[i]; - if (min_val < val) { - min_val = val; min_pos = i; - } - } - } - - return min_pos; -} - -data_t *fracture_network(net_t *instance, double *fuse_thres, - cholmod_common *c, double cutoff) { - unsigned int num_edges = instance->graph->ne; - unsigned int num_verts = instance->graph->nv; - - data_t *breaking_data = alloc_break_data(num_edges); - - while (true) { - double *voltages = get_voltage(instance, c); - double *field = get_current_v(instance, voltages, c); - - double conductivity = get_conductivity(instance, voltages, c); - if (conductivity < 1e-12 && instance->voltage_bound) { - free(voltages); - free(field); - break; - } - - int last_broke = inc_break_fuses(instance, fuse_thres, field, cutoff); - if (last_broke > num_edges || last_broke < -1 || conductivity < 1e-8) { - printf("whoops %u\n\n", breaking_data->num_broken); - free(voltages); - free(field); - break; - } - - update_break_data(breaking_data, last_broke, fabs(conductivity * fuse_thres[last_broke] / field[last_broke]), conductivity); - - free(voltages); - free(field); - - break_edge(instance, last_broke, c); - - if (instance->num_components > 1 && instance->graph->boundary == TORUS_BOUND) { - break; - } - - if (instance->marks[num_verts] != instance->marks[num_verts + 1] && instance->graph->boundary != TORUS_BOUND) { - break; - } - } - - return breaking_data; -} - diff --git a/src/free_network.c b/src/free_network.c index 5fd3290..2001479 100644 --- a/src/free_network.c +++ b/src/free_network.c @@ -1,7 +1,7 @@ #include "fracture.h" -void free_net(graph_t *network, cholmod_common *c) { +void graph_free(graph_t *network, cholmod_common *c) { free(network->ev); if (network->ev_break != network->ev) { free(network->ev_break); @@ -10,10 +10,10 @@ void free_net(graph_t *network, cholmod_common *c) { free(network->ve); free(network->bound_inds); free(network->bound_verts); - free(network->vert_coords); - free(network->edge_coords); + free(network->vx); + free(network->ex); free(network->dev); - free(network->dual_vert_coords); + free(network->dvx); free(network->dvei); free(network->dve); free(network->spanning_edges); diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c index 6034a82..4196937 100644 --- a/src/gen_laplacian.c +++ b/src/gen_laplacian.c @@ -3,14 +3,9 @@ cholmod_sparse *gen_adjacency(const net_t *instance, bool dual, bool breakv, unsigned int pad, cholmod_common *c) { - unsigned int ne; - if (dual) - ne = ((int)instance->graph->ne); - else { - ne = instance->num_remaining_edges; - if (!breakv && instance->graph->boundary != TORUS_BOUND) { - ne += instance->graph->bound_inds[2]; - } + unsigned int ne = instance->graph->ne; + if (!breakv && instance->graph->boundary != TORUS_BOUND && !dual) { + ne += instance->graph->bound_inds[2]; } unsigned int nnz = 2 * ne; @@ -58,7 +53,7 @@ cholmod_sparse *gen_adjacency(const net_t *instance, bool dual, bool breakv, ai[2 * count + 1] = 1; count++; - } else if (dual) { + } else { unsigned int v1 = etv[2 * i]; unsigned int v2 = etv[2 * i + 1]; @@ -99,7 +94,7 @@ cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c, bool symmetric) { const graph_t *network = instance->graph; unsigned int num_verts = network->nv_break; - double *vert_coords = network->vert_coords; + double *vert_coords = network->vx; unsigned int num_bounds = network->num_bounds; double inf = instance->inf; bool voltage_bound = instance->voltage_bound; diff --git a/src/get_conductivity.c b/src/get_conductivity.c index 793987e..8c4d228 100644 --- a/src/get_conductivity.c +++ b/src/get_conductivity.c @@ -9,8 +9,8 @@ double get_conductivity(net_t *inst, double *voltage, cholmod_common *c) { if (!inst->fuses[e]) { unsigned int v1 = inst->graph->ev[2*e]; unsigned int v2 = inst->graph->ev[2*e+1]; - double v1y = inst->graph->vert_coords[2 * v1 + 1]; - double v2y = inst->graph->vert_coords[2 * v2 + 1]; + double v1y = inst->graph->vx[2 * v1 + 1]; + double v2y = inst->graph->vx[2 * v2 + 1]; unsigned int s1 = v1y < v2y ? v1 : v2; unsigned int s2 = v1y < v2y ? v2 : v1; tot_cur += voltage[s1] - voltage[s2]; diff --git a/src/homo_square_fracture.c b/src/homo_square_fracture.c index e3e8ad3..b301136 100644 --- a/src/homo_square_fracture.c +++ b/src/homo_square_fracture.c @@ -193,7 +193,7 @@ int main(int argc, char *argv[]) { for (unsigned int i = 0; i < N; i++) { printf("\033[F\033[JFRACTURE: %0*d / %d\n", (int)log10(N) + 1, i + 1, N); - double *fuse_thres = gen_fuse_thres(network->ne, network->edge_coords, beta, beta_scaling_flat); + double *fuse_thres = gen_fuse_thres(network->ne, network->ex, beta, beta_scaling_flat); net_t *instance = copy_instance(perm_instance, &c); data_t *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); free_instance(instance, &c); @@ -287,8 +287,8 @@ int main(int argc, char *argv[]) { if (save_network) { FILE *net_out = fopen("network.txt", "w"); for (unsigned int j = 0; j < network->nv; j++) { - fprintf(net_out, "%f %f ", network->vert_coords[2 * j], - tmp_instance->graph->vert_coords[2 * j + 1]); + fprintf(net_out, "%f %f ", network->vx[2 * j], + tmp_instance->graph->vx[2 * j + 1]); } fprintf(net_out, "\n"); for (unsigned int j = 0; j < tmp_instance->graph->ne; j++) { @@ -297,8 +297,8 @@ int main(int argc, char *argv[]) { } fprintf(net_out, "\n"); for (unsigned int j = 0; j < tmp_instance->graph->dnv; j++) { - fprintf(net_out, "%f %f ", tmp_instance->graph->dual_vert_coords[2 * j], - tmp_instance->graph->dual_vert_coords[2 * j + 1]); + fprintf(net_out, "%f %f ", tmp_instance->graph->dvx[2 * j], + tmp_instance->graph->dvx[2 * j + 1]); } fprintf(net_out, "\n"); for (unsigned int j = 0; j < tmp_instance->graph->ne; j++) { diff --git a/src/homo_voronoi_fracture.c b/src/homo_voronoi_fracture.c index 6171480..85a9c2b 100644 --- a/src/homo_voronoi_fracture.c +++ b/src/homo_voronoi_fracture.c @@ -211,9 +211,9 @@ int main(int argc, char *argv[]) { 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 *network = ini_voronoi_network(L, boundary, use_dual, genfunc_hyperuniform, &c); + graph_t *network = ini_voro_graph(L, boundary, use_dual, genfunc_hyperuniform, &c); net_t *perm_instance = create_instance(network, inf, use_voltage_boundaries, true, &c); - double *fuse_thres = gen_fuse_thres(network->ne, network->edge_coords, beta, beta_scaling_flat); + double *fuse_thres = gen_fuse_thres(network->ne, network->ex, beta, beta_scaling_flat); net_t *instance = copy_instance(perm_instance, &c); data_t *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff); free_instance(instance, &c); @@ -299,8 +299,8 @@ int main(int argc, char *argv[]) { if (save_network) { FILE *net_out = fopen("network.txt", "w"); for (uint_t j = 0; j < network->nv; j++) { - fprintf(net_out, "%f %f ", network->vert_coords[2 * j], - tmp_instance->graph->vert_coords[2 * j + 1]); + fprintf(net_out, "%f %f ", network->vx[2 * j], + tmp_instance->graph->vx[2 * j + 1]); } fprintf(net_out, "\n"); for (uint_t j = 0; j < tmp_instance->graph->ne; j++) { @@ -309,8 +309,8 @@ int main(int argc, char *argv[]) { } fprintf(net_out, "\n"); for (uint_t j = 0; j < tmp_instance->graph->dnv; j++) { - fprintf(net_out, "%f %f ", tmp_instance->graph->dual_vert_coords[2 * j], - tmp_instance->graph->dual_vert_coords[2 * j + 1]); + fprintf(net_out, "%f %f ", tmp_instance->graph->dvx[2 * j], + tmp_instance->graph->dvx[2 * j + 1]); } fprintf(net_out, "\n"); for (uint_t j = 0; j < tmp_instance->graph->ne; j++) { diff --git a/src/ini_network.c b/src/ini_network.c index ebfea76..cfb8d53 100644 --- a/src/ini_network.c +++ b/src/ini_network.c @@ -219,21 +219,21 @@ graph_t *ini_square_network(unsigned int width, bound_t boundary, bool side_boun } free(vert_counts); - network->vert_coords = + network->vx = (double *)malloc(2 * network->nv * sizeof(double)); for (unsigned int i = 0; i < network->nv; i++) { if (!periodic) { - network->vert_coords[2 * i] = ((double)((2 * i + 1) / (width + 1)))/width; - network->vert_coords[2 * i + 1] = ((double)((2 * i + 1) % (width + 1)))/width; + network->vx[2 * i] = ((double)((2 * i + 1) / (width + 1)))/width; + network->vx[2 * i + 1] = ((double)((2 * i + 1) % (width + 1)))/width; } else { - network->vert_coords[2 * i] = ((double)((2 * i + 1) / (width)))/width; - network->vert_coords[2 * i + 1] = ((double)((((2 * i + 1) / (width)+1) % 2) + ((2 * i) % (width))))/width; + network->vx[2 * i] = ((double)((2 * i + 1) / (width)))/width; + network->vx[2 * i + 1] = ((double)((((2 * i + 1) / (width)+1) % 2) + ((2 * i) % (width))))/width; } } network->L = width; - network->edge_coords = get_edge_coords( - network->ne, network->vert_coords, network->ev); + network->ex = get_edge_coords( + network->ne, network->vx, network->ev); network->dev = (unsigned int *)malloc(2 * network->ne * sizeof(unsigned int)); @@ -243,12 +243,12 @@ graph_t *ini_square_network(unsigned int width, bound_t boundary, bool side_boun network->dev[2 * i + 1] = dual_edge_to_verts(width, periodic, i, 1) % network->dnv; } - network->dual_vert_coords = + network->dvx = (double *)malloc(2 * network->dnv * sizeof(double)); for (unsigned int i = 0; i < network->dnv; i++) { - network->dual_vert_coords[2 * i] = + network->dvx[2 * i] = 2*dual_vert_to_coord(width, periodic, i, 0); - network->dual_vert_coords[2 * i + 1] = + network->dvx[2 * i + 1] = 2*dual_vert_to_coord(width, periodic, i, 1); } @@ -262,7 +262,7 @@ graph_t *ini_square_network(unsigned int width, bound_t boundary, bool side_boun network->dve = get_verts_to_edges( network->dnv, network->ne, network->dev, network->dvei); - network->spanning_edges = get_spanning_edges(network->ne, network->ev_break, network->vert_coords, 0.51, &(network->num_spanning_edges)); + network->spanning_edges = get_spanning_edges(network->ne, network->ev_break, network->vx, 0.51, &(network->num_spanning_edges)); return network; } @@ -304,10 +304,10 @@ unsigned int *get_voro_dual_edges(unsigned int num_edges, return dual_edges; } -graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual, +graph_t *ini_voro_graph(unsigned int L, bound_t boundary, bool use_dual, double *(*genfunc)(unsigned int, bound_t, gsl_rng *, unsigned int *), cholmod_common *c) { - graph_t *network = (graph_t *)calloc(1, sizeof(graph_t)); + graph_t *g = (graph_t *)calloc(1, sizeof(graph_t)); // generate the dual lattice double *lattice; @@ -459,11 +459,11 @@ graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual, free(tmp_edges); free(tmp_dual_edges); num_bounds = 2; - network->ev_break = edges; - network->ev = edges; - network->nv_break = num_verts; - network->nv = num_verts; - network->break_dim = num_verts + num_bounds; + g->ev_break = edges; + g->ev = edges; + g->nv_break = num_verts; + g->nv = num_verts; + g->break_dim = num_verts + num_bounds; break; } case CYLINDER_BOUND: { @@ -526,11 +526,11 @@ graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual, free(bound_top); free(bound_b); free(tmp_edges); free(tmp_dual_edges); - network->ev_break = edges; - network->ev = edges; - network->nv_break = num_verts; - network->nv = num_verts; - network->break_dim = num_verts + num_bounds; + g->ev_break = edges; + g->ev = edges; + g->nv_break = num_verts; + g->nv = num_verts; + g->break_dim = num_verts + num_bounds; break; } case TORUS_BOUND: { @@ -596,11 +596,11 @@ graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual, free(tmp_vert_coords); free(bound_top); free(edge_change); - network->nv_break = num_verts; - network->nv = tmp_num_verts; - network->ev_break = edges; - network->ev = tmp_edges; - network->break_dim = num_verts; + g->nv_break = num_verts; + g->nv = tmp_num_verts; + g->ev_break = edges; + g->ev = tmp_edges; + g->break_dim = num_verts; break; } case EMBEDDED_BOUND: { @@ -618,50 +618,49 @@ graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual, edges = tmp_edges; dual_edges = tmp_dual_edges; vert_coords = tmp_vert_coords; - network->nv_break = num_verts; - network->nv = num_verts; - network->ev_break = edges; - network->ev = edges; - network->break_dim = num_verts + 2; + g->nv_break = num_verts; + g->nv = num_verts; + g->ev_break = edges; + g->ev = edges; + g->break_dim = num_verts + 2; } } - network->boundary = boundary; - network->num_bounds = num_bounds; - network->bound_inds = bound_inds; - network->bound_verts = bound_verts; - network->ne = num_edges; - network->dev = dual_edges; - network->vert_coords = vert_coords; + g->boundary = boundary; + g->num_bounds = num_bounds; + g->bound_inds = bound_inds; + g->bound_verts = bound_verts; + g->ne = num_edges; + g->dev = dual_edges; + g->vx = vert_coords; - network->vei = get_verts_to_edges_ind( - network->nv, network->ne, network->ev); - network->ve = - get_verts_to_edges(network->nv, network->ne, - network->ev, network->vei); + g->vei = get_verts_to_edges_ind(g->nv, g->ne, g->ev); + g->ve = + get_verts_to_edges(g->nv, g->ne, + g->ev, g->vei); - network->L = 1; + g->L = L; - network->edge_coords = get_edge_coords( - network->ne, network->vert_coords, network->ev); + g->ex = get_edge_coords( + g->ne, g->vx, g->ev); free(tmp_tris); - network->dual_vert_coords = lattice; - network->dnv = num; + g->dvx = lattice; + g->dnv = num; - network->voltcurmat = gen_voltcurmat(network->ne, - network->break_dim, - network->ev_break, c); + g->voltcurmat = gen_voltcurmat(g->ne, + g->break_dim, + g->ev_break, c); - network->dvei = - get_verts_to_edges_ind(network->dnv, network->ne, - network->dev); - network->dve = get_verts_to_edges( - network->dnv, network->ne, network->dev, - network->dvei); + g->dvei = + get_verts_to_edges_ind(g->dnv, g->ne, + g->dev); + g->dve = get_verts_to_edges( + g->dnv, g->ne, g->dev, + g->dvei); - network->spanning_edges = get_spanning_edges(network->ne, network->ev_break, network->vert_coords, 0.5, &(network->num_spanning_edges)); + g->spanning_edges = get_spanning_edges(g->ne, g->ev_break, g->vx, 0.5, &(g->num_spanning_edges)); - return network; + return g; } diff --git a/src/instance.c b/src/instance.c deleted file mode 100644 index bb1ac8c..0000000 --- a/src/instance.c +++ /dev/null @@ -1,124 +0,0 @@ - -#include "fracture.h" - -net_t *create_instance(graph_t *network, double inf, bool voltage_bound, - bool startnow, cholmod_common *c) { - net_t *instance = (net_t *)calloc(1, sizeof(net_t)); - assert(instance != NULL); - - instance->graph = network; - instance->num_remaining_edges = network->ne; - instance->fuses = (bool *)calloc(network->ne, sizeof(bool)); - assert(instance->fuses != NULL); - instance->inf = inf; - instance->voltage_bound = voltage_bound; - instance->boundary_cond = CHOL_F(zeros)( - network->break_dim, 1, CHOLMOD_REAL, c); - if (network->boundary == TORUS_BOUND) { - for (unsigned int i = 0; i < network->bound_inds[1]; i++) { - ((double *)instance->boundary_cond->x)[network->bound_verts[i]] = 1; - ((double *)instance->boundary_cond->x)[network->nv + i] = -1; - } - ((double *)instance->boundary_cond->x)[network->bound_verts[0]] = 1; - } else if (network->boundary == EMBEDDED_BOUND) { - // do nothing - for (unsigned int i = 0; i < network->bound_inds[1]; i++) { - ((double *)instance->boundary_cond->x)[network->bound_verts[i]] =1; - } - } else { - if (voltage_bound) { - for (unsigned int i = 0; i < network->bound_inds[1]; i++) { - ((double *)instance->boundary_cond->x)[network->bound_verts[i]] =1; - } - } else { - ((double *)instance->boundary_cond->x)[0] = 1; - ((double *)instance->boundary_cond->x)[network->nv] = 1; - ((double *)instance->boundary_cond->x)[network->nv + 1] = -1; - } - } - - if (network->boundary != TORUS_BOUND) instance->adjacency = gen_adjacency(instance, false, false, 0, c); - else instance->adjacency = gen_adjacency(instance, true, false, 0, c); - - if (startnow) { - cholmod_sparse *laplacian = gen_laplacian(instance, c, true); - instance->factor = CHOL_F(analyze)(laplacian, c); - CHOL_F(factorize)(laplacian, instance->factor, c); - CHOL_F(free_sparse)(&laplacian, c); - } - - instance->marks = (unsigned int *)malloc( - (instance->graph->break_dim) * - sizeof(unsigned int)); - instance->dual_marks = (unsigned int *)malloc( - (instance->graph->dnv) * - sizeof(unsigned int)); - assert(instance->marks != NULL); - - for (unsigned int i = 0; - i < (instance->graph->break_dim); - i++) { - instance->marks[i] = 1; - } - for (unsigned int i = 0; - i < (instance->graph->dnv); - i++) { - instance->dual_marks[i] = i+1; - } - instance->num_components = 1; - - return instance; -} - -void finish_instance(net_t *instance, cholmod_common *c) { - cholmod_sparse *laplacian = gen_laplacian(instance, c, true); - instance->factor = CHOL_F(analyze)(laplacian, c); - CHOL_F(factorize)(laplacian, instance->factor, c); - CHOL_F(free_sparse)(&laplacian, c); -} - -net_t *copy_instance(const net_t *instance, cholmod_common *c) { - net_t *instance_copy = (net_t *)calloc(1, sizeof(net_t)); - memcpy(instance_copy, instance, sizeof(net_t)); - - size_t fuses_size = (instance->graph)->ne * sizeof(bool); - instance_copy->fuses = (bool *)malloc(fuses_size); - memcpy(instance_copy->fuses, instance->fuses, fuses_size); - - size_t marks_size = - (instance->graph->break_dim) * - sizeof(unsigned int); - instance_copy->marks = (unsigned int *)malloc(marks_size); - memcpy(instance_copy->marks, instance->marks, marks_size); - - size_t dual_marks_size = - (instance->graph->dnv) * - sizeof(unsigned int); - instance_copy->dual_marks = (unsigned int *)malloc(dual_marks_size); - memcpy(instance_copy->dual_marks, instance->dual_marks, dual_marks_size); - - instance_copy->adjacency = CHOL_F(copy_sparse)(instance->adjacency, c); - instance_copy->boundary_cond = CHOL_F(copy_dense)(instance->boundary_cond, c); - instance_copy->factor = CHOL_F(copy_factor)(instance->factor, c); - - return instance_copy; -} - -void free_instance(net_t *instance, cholmod_common *c) { - free(instance->fuses); - CHOL_F(free_dense)(&(instance->boundary_cond), c); - CHOL_F(free_sparse)(&(instance->adjacency), c); - CHOL_F(free_factor)(&(instance->factor), c); - free(instance->marks); - free(instance->dual_marks); - free(instance); -} - -bool check_instance(const net_t *instance, cholmod_common *c) { - assert(instance != NULL); - assert(instance->fuses != NULL); - assert(CHOL_F(check_dense)(instance->boundary_cond, c)); - assert(CHOL_F(check_factor)(instance->factor, c)); - - return true; -} diff --git a/src/net.c b/src/net.c new file mode 100644 index 0000000..d090a9a --- /dev/null +++ b/src/net.c @@ -0,0 +1,121 @@ + +#include "fracture.h" + +double *get_thres(uint_t ne, double beta) { + assert(beta > 0); + + double *thres = (double *)malloc(ne * sizeof(double)); + assert(thres != NULL); + + gsl_set_error_handler_off(); + + gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); + { + FILE *rf = fopen("/dev/urandom", "r"); + unsigned long int seed; + fread(&seed, sizeof(unsigned long int), 1, rf); + fclose(rf); + gsl_rng_set(r, seed); + } + + for (uint_t i = 0; i < ne; i++) { + while ((thres[i] = exp(log(gsl_ran_flat(r, 0, 1)) / beta)) == 0.0); + } + + gsl_rng_free(r); + + return thres; +} + + +net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c) { + net_t *net = (net_t *)calloc(1, sizeof(net_t)); + assert(net != NULL); + + net->graph = g; + net->fuses = (bool *)calloc(g->ne, sizeof(bool)); + assert(net->fuses != NULL); + net->thres = get_thres(g->ne, beta); + net->inf = inf; + + net->voltage_bound = vb; + net->boundary_cond = bound_set(g, vb, notch_len, c); + + if (g->boundary != TORUS_BOUND) net->adjacency = gen_adjacency(net, false, false, 0, c); + else net->adjacency = gen_adjacency(net, true, false, 0, c); + + net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t)); + net->dual_marks = (uint_t *)malloc((net->graph->dnv) * sizeof(uint_t)); + assert(net->marks != NULL); + + for (uint_t i = 0; i < (net->graph->break_dim); i++) { + net->marks[i] = 1; + } + for (uint_t i = 0; i < (net->graph->dnv); i++) { + net->dual_marks[i] = i+1; + } + net->num_components = 1; + + net_notch(net, notch_len, c); + + { + cholmod_sparse *laplacian = gen_laplacian(net, c, true); + net->factor = CHOL_F(analyze)(laplacian, c); + CHOL_F(factorize)(laplacian, net->factor, c); + CHOL_F(free_sparse)(&laplacian, c); + } + + return net; +} + +net_t *net_copy(const net_t *net, cholmod_common *c) { + net_t *net_copy = (net_t *)calloc(1, sizeof(net_t)); + assert(net_copy != NULL); + memcpy(net_copy, net, sizeof(net_t)); + + size_t fuses_size = (net->graph)->ne * sizeof(bool); + net_copy->fuses = (bool *)malloc(fuses_size); + assert(net_copy->fuses != NULL); + memcpy(net_copy->fuses, net->fuses, fuses_size); + + size_t thres_size = (net->graph)->ne * sizeof(double); + net_copy->thres = (double *)malloc(thres_size); + assert(net_copy->thres != NULL); + memcpy(net_copy->thres, net->thres, thres_size); + + size_t marks_size = (net->graph->break_dim) * sizeof(uint_t); + net_copy->marks = (uint_t *)malloc(marks_size); + assert(net_copy->marks != NULL); + memcpy(net_copy->marks, net->marks, marks_size); + + size_t dual_marks_size = (net->graph->dnv) * sizeof(uint_t); + net_copy->dual_marks = (uint_t *)malloc(dual_marks_size); + assert(net_copy->dual_marks != NULL); + memcpy(net_copy->dual_marks, net->dual_marks, dual_marks_size); + + net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c); + net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c); + net_copy->factor = CHOL_F(copy_factor)(net->factor, c); + + return net_copy; +} + +void net_free(net_t *net, cholmod_common *c) { + free(net->fuses); + free(net->thres); + CHOL_F(free_dense)(&(net->boundary_cond), c); + CHOL_F(free_sparse)(&(net->adjacency), c); + CHOL_F(free_factor)(&(net->factor), c); + free(net->marks); + free(net->dual_marks); + free(net); +} + +bool check_instance(const net_t *instance, cholmod_common *c) { + assert(instance != NULL); + assert(instance->fuses != NULL); + assert(CHOL_F(check_dense)(instance->boundary_cond, c)); + assert(CHOL_F(check_factor)(instance->factor, c)); + + return true; +} diff --git a/src/net_fracture.c b/src/net_fracture.c new file mode 100644 index 0000000..b7fa61d --- /dev/null +++ b/src/net_fracture.c @@ -0,0 +1,66 @@ + +#include "fracture.h" + +uint_t get_next_broken(net_t *net, double *currents, double cutoff) { + uint_t max_pos = UINT_MAX; + double max_val = 0; + + for (uint_t i = 0; i < net->graph->ne; i++) { + double current = fabs(currents[i]); + bool broken = net->fuses[i]; + + if (!broken && current > cutoff) { + double val = current / net->thres[i]; + + if (val > max_val) { + max_val = val; + max_pos = i; + } + } + } + + if (max_pos == UINT_MAX) { + printf("GET_NEXT_BROKEN: currents is zero or NaN, no max_val found\n"); + exit(EXIT_FAILURE); + } + + return max_pos; +} + + +data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) { + data_t *data = alloc_break_data(net->graph->ne); + + while (true) { + double *voltages = get_voltage(net, c); + double *currents = get_current_v(net, voltages, c); + + double conductivity = get_conductivity(net, voltages, c); + + if (conductivity < 1e-12 && net->graph->boundary == TORUS_BOUND) { + free(voltages); + free(currents); + break; + } + + uint_t last_broke = get_next_broken(net, currents, cutoff); + + update_break_data(data, last_broke, fabs(conductivity * (net->thres)[last_broke] / currents[last_broke]), conductivity); + + free(voltages); + free(currents); + + break_edge(net, last_broke, c); + + if (net->num_components > 1 && net->graph->boundary == TORUS_BOUND) { + break; + } + + if (net->marks[net->graph->nv] != net->marks[net->graph->nv + 1] && net->graph->boundary != TORUS_BOUND) { + break; + } + } + + return data; +} + diff --git a/src/net_notch.c b/src/net_notch.c new file mode 100644 index 0000000..608d0b5 --- /dev/null +++ b/src/net_notch.c @@ -0,0 +1,29 @@ + +#include "fracture.h" + +void net_notch(net_t *net, double notch_len, cholmod_common *c) { + for (uint_t i = 0; i < net->graph->ne; i++) { + uint_t v1, v2; + double v1x, v1y, v2x, v2y, dx, 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]; + + dx = v1x - v2x; + 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; + + if (crosses_center && not_wrapping && correct_length) { + break_edge(net, i, c); + } + } +} diff --git a/src/voro_fracture.c b/src/voro_fracture.c new file mode 100644 index 0000000..237192a --- /dev/null +++ b/src/voro_fracture.c @@ -0,0 +1,492 @@ + + +#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_stress_field, save_voltage_field, save_toughness, save_conductivity, + save_damage, save_damage_field; + bound_t boundary; + + + // 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; + save_data = false; + save_cluster_dist = false; + use_voltage_boundaries = false; + use_dual = false; + save_network = false; + save_crit_stress = false; + save_stress_field = false; + save_voltage_field = false; + save_damage = false; + save_damage_field = false; + save_conductivity = false; + save_toughness = false; + + + uint8_t bound_i; + char boundc2 = 'f'; + + + // get commandline options + + while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrtDSvel:")) != -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 'd': + save_damage = true; + break; + case 'e': + save_damage_field = true; + break; + case 'V': + use_voltage_boundaries = true; + break; + case 'D': + use_dual = true; + 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 'S': + save_stress_field = true; + break; + case 'v': + save_voltage_field = true; + break; + case 'r': + save_conductivity = true; + //inf = 1; + break; + case 't': + save_toughness = 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_v_%c_%c_%u_%g_%g.txt", boundc, boundc2, L, beta, crack_len); + data_out = fopen(data_filename, "a"); + free(data_filename); + } + + uint_t voronoi_max_verts, c_dist_size, a_dist_size; + + voronoi_max_verts = 4 * pow(L, 2); + c_dist_size = voronoi_max_verts; + a_dist_size = voronoi_max_verts; + + if (voronoi_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 *)malloc(c_dist_size * sizeof(uint32_t)); + avalanche_size_dist = + (uint32_t *)malloc(a_dist_size * 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_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + snprintf(a_filename, filename_len, "avln_v_%c_%c_%d_%g_%g.dat", 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), c_dist_size, cluster_out); + fclose(cluster_out); + } + if (avalanche_out != NULL) { + fread(avalanche_size_dist, sizeof(uint32_t), a_dist_size, avalanche_out); + fclose(avalanche_out); + } + } + + double *crit_stress; + if (save_crit_stress) { + crit_stress = (double *)malloc(N * sizeof(double)); + } + + double *stress_field; + unsigned int stress_pos = 0; + if (save_stress_field) { + stress_field = (double *)malloc(3 * N * voronoi_max_verts * sizeof(double)); + } + + double *voltage_field; + unsigned int voltage_pos = 0; + if (save_voltage_field) { + voltage_field = (double *)malloc(3 * N * voronoi_max_verts * sizeof(double)); + } + + double *damage_field; + unsigned int damage_pos = 0; + if (save_damage_field) { + damage_field = (double *)malloc(2 * N * voronoi_max_verts * sizeof(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 *)malloc(a_dist_size * sizeof(uint32_t)); + + d_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(d_filename, filename_len, "damg_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + + FILE *damage_out = fopen(d_filename, "rb"); + + if (damage_out != NULL) { + fread(damage, sizeof(uint32_t), a_dist_size, damage_out); + fclose(damage_out); + } + } + + double *toughness; + if (save_toughness) { + toughness = (double *)malloc(N * sizeof(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 = ini_voro_graph(L, boundary, use_dual, genfunc_hyperuniform, &c); + 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; + double max_val = 0; + + for (uint_t j = 0; j < data->num_broken; j++) { + double val = data->extern_field[j]; + + if (val > max_val) { + max_pos = j; + max_val = val; + } + } + + if (save_crit_stress) crit_stress[i] = data->extern_field[max_pos]; + + if (save_conductivity) conductivity[i] = data->conductivity[max_pos]; + + if (save_damage) damage[max_pos]++; + + uint_t av_size = 0; + double cur_val = 0; + for (uint_t j = 0; j < max_pos; j++) { + break_edge(net, data->break_list[j], &c); + + 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_stress_field || save_voltage_field) { + double *tmp_voltages = get_voltage(net, &c); + if (save_voltage_field) { + for (uint_t j = 0; j < g->nv; j++) { + voltage_field[3 * voltage_pos] = g->vx[2 * j]; + voltage_field[3 * voltage_pos + 1] = g->vx[2 * j + 1]; + voltage_field[3 * voltage_pos + 2] = tmp_voltages[j]; + voltage_pos++; + } + } + if (save_stress_field) { + double *tmp_currents = get_current_v(net, tmp_voltages, &c); + for (uint_t j = 0; j < g->ne; j++) { + stress_field[3 * stress_pos] = g->ex[2 * j]; + stress_field[3 * stress_pos + 1] = g->ex[2 * j + 1]; + stress_field[3 * stress_pos + 2] = tmp_currents[j]; + stress_pos++; + } + free(tmp_currents); + } + free(tmp_voltages); + } + + if (save_damage_field) { + for (uint_t j = 0; j < max_pos; j++) { + damage_field[2 * damage_pos] = g->ex[2 * data->break_list[j]]; + damage_field[2 * damage_pos + 1] = g->ex[2 * data->break_list[j] + 1]; + damage_pos++; + } + } + + if (save_toughness) { + double tmp_toughness = 0; + if (max_pos > 0) { + double sigma1 = data->extern_field[0]; + double epsilon1 = sigma1 / data->conductivity[0]; + for (uint_t j = 0; j < max_pos - 1; j++) { + double sigma2 = data->extern_field[j+1]; + double epsilon2 = sigma2 / data->conductivity[j+1]; + if (epsilon2 > epsilon1) { + tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2; + sigma1 = sigma2; epsilon1 = epsilon2; + } + } + } + toughness[i] = tmp_toughness; + } + + if (save_cluster_dist) { + uint_t *tmp_cluster_dist = get_cluster_dist(net, &c); + for (uint_t j = 0; j < g->dnv; j++) { + cluster_size_dist[j] += tmp_cluster_dist[j]; + } + free(tmp_cluster_dist); + } + + 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); + } + + net_free(net, &c); + graph_free(g, &c); + + if (save_data) { + for (uint_t j = 0; j < data->num_broken; j++) { + fprintf(data_out, "%u %g %g ", data->break_list[j], + data->extern_field[j], data->conductivity[j]); + } + fprintf(data_out, "\n"); + } + + free_break_data(data); + } + + 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), c_dist_size, cluster_out); + fwrite(avalanche_size_dist, sizeof(uint32_t), a_dist_size, avalanche_out); + + fclose(cluster_out); + fclose(avalanche_out); + + free(c_filename); + free(a_filename); + free(cluster_size_dist); + free(avalanche_size_dist); + } + + if (save_voltage_field) { + char *vfld_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(vfld_filename, filename_len, "vfld_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + FILE *vfld_file = fopen(vfld_filename, "ab"); + fwrite(voltage_field, sizeof(double), 3 * voltage_pos, vfld_file); + fclose(vfld_file); + free(vfld_filename); + free(voltage_field); + } + + if (save_stress_field) { + char *cfld_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(cfld_filename, filename_len, "cfld_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + FILE *cfld_file = fopen(cfld_filename, "ab"); + fwrite(stress_field, sizeof(double), 3 * stress_pos, cfld_file); + fclose(cfld_file); + free(cfld_filename); + free(stress_field); + } + + if (save_damage_field) { + char *dfld_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(dfld_filename, filename_len, "dfld_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + FILE *dfld_file = fopen(dfld_filename, "ab"); + fwrite(damage_field, sizeof(double), 2 * damage_pos, dfld_file); + fclose(dfld_file); + free(dfld_filename); + free(damage_field); + } + + if (save_conductivity) { + char *cond_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(cond_filename, filename_len, "cond_v_%c_%c_%d_%g_%g.dat", 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_toughness) { + char *tough_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(tough_filename, filename_len, "tuff_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + FILE *tough_file = fopen(tough_filename, "ab"); + fwrite(toughness, sizeof(double), N, tough_file); + fclose(tough_file); + free(tough_filename); + free(toughness); + } + + if (save_damage) { + FILE *hdam_file = fopen(d_filename, "wb"); + fwrite(damage, sizeof(uint32_t), a_dist_size, 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_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + FILE *str_file = fopen(str_filename, "ab"); + fwrite(crit_stress, sizeof(double), N, str_file); + fclose(str_file); + free(str_filename); + free(crit_stress); + } + + CHOL_F(finish)(&c); + + return 0; +} diff --git a/src/voronoi_bound_ini.c b/src/voronoi_bound_ini.c deleted file mode 100644 index 0b65ef5..0000000 --- a/src/voronoi_bound_ini.c +++ /dev/null @@ -1,33 +0,0 @@ - -#include "fracture.h" - -double th_p(double x, double y, double th) { - if (x >= 0 && y >= 0) return th; - if (x < 0 && y >= 0) return M_PI - th; - if (x < 0 && y < 0) return th - M_PI; - if (x >= 0 && y < 0) return -th; -} - -double u_y(double x, double y) { - double r = sqrt(pow(x, 2) + pow(y, 2)); - double th = th_p(x, y, atan(fabs(y / x))); - - return sqrt(r) * sin(th / 2); -} - -void voronoi_bound_ini(net_t *instance, uint_t L, double crack_len) { - double *bound = (double *)instance->boundary_cond->x; - for (uint_t i = 0; i < L / 2; i++) { - double x1, y1, x2, y2, x3, y3, x4, y4; - x1 = (2. * i + 1.) / L - crack_len; y1 = 0.5 - 1.; - x2 = (2. * i + 1.) / L - crack_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[i] = u_y(x1, y1); - bound[L / 2 + i] = u_y(x2, y2); - bound[L + i] = u_y(x3, y3); - bound[3 * L / 2 + i] = u_y(x4, y4); - } -} - -- cgit v1.2.3-70-g09d2