diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/beta_scales.c | 21 | ||||
| -rw-r--r-- | src/bin_values.c | 23 | ||||
| -rw-r--r-- | src/compare_voronoi_fracture.c | 124 | ||||
| -rw-r--r-- | src/course_grain_square.c | 168 | ||||
| -rw-r--r-- | src/coursegrain.c | 40 | ||||
| -rw-r--r-- | src/current_scaling.c | 300 | ||||
| -rw-r--r-- | src/get_file.c | 38 | ||||
| -rw-r--r-- | src/homo_square_fracture.c | 418 | ||||
| -rw-r--r-- | src/homo_voronoi_fracture.c | 405 | ||||
| -rw-r--r-- | src/toy_setup.c | 37 | ||||
| -rw-r--r-- | src/update_beta.c | 34 | ||||
| -rw-r--r-- | src/update_boundary.c | 86 | ||||
| -rw-r--r-- | src/voro_fracture.c | 492 |
13 files changed, 0 insertions, 2186 deletions
diff --git a/src/beta_scales.c b/src/beta_scales.c deleted file mode 100644 index 8cbc5d1..0000000 --- a/src/beta_scales.c +++ /dev/null @@ -1,21 +0,0 @@ - -#include "fracture.h" - -double beta_scaling_flat(double beta, double x, double y) { return beta; } - -double beta_mag(double beta) { - double aa = -7.52579; - double bb = 9.63706; - double cc = 0.692515; - double dd = -2.47638; - - return gsl_sf_exp(aa + bb * gsl_sf_erf(cc * (gsl_sf_log(beta) - dd))); -} - -double beta_scaling_gaus(double beta, double x, double y) { - double sigma = 0.25; - double nu_f = 1.56; - double betap = beta * gsl_sf_exp((gsl_pow_2(x - 0.5) + gsl_pow_2(y - 0.5)) / - (2 * sigma * 2 * nu_f)); - return betap; -} diff --git a/src/bin_values.c b/src/bin_values.c deleted file mode 100644 index 70009c1..0000000 --- a/src/bin_values.c +++ /dev/null @@ -1,23 +0,0 @@ - -#include "fracture.h" - -double *bin_values(graph_t *network, unsigned int width, double *values) { - double *binned = calloc(pow(width, 2), sizeof(double)); - 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->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]++; - } - } - - for (unsigned int i = 0; i < pow(width, 2); i++) { - if (num_binned[i] != 0) { - binned[i] /= num_binned[i]; - } - } - - return binned; -} diff --git a/src/compare_voronoi_fracture.c b/src/compare_voronoi_fracture.c deleted file mode 100644 index 91fdcea..0000000 --- a/src/compare_voronoi_fracture.c +++ /dev/null @@ -1,124 +0,0 @@ - - -#include "fracture.h" - -int main(int argc, char *argv[]) { - int opt; - - // defining variables to be (potentially) set by command line flags - unsigned int N, L, filename_len; - double beta, inf, cutoff; - bound_t boundary; - - filename_len = 100; - - N = 100; - L = 16; - beta = .3; - inf = 1e10; - cutoff = 1e-10; - boundary = FREE_BOUND; - - int boundary_int; - char boundc2 = 'f'; - - while ((opt = getopt(argc, argv, "n:L:b:B:dcoNsCrt")) != -1) { - switch (opt) { - case 'n': - N = atoi(optarg); - break; - case 'L': - L = atoi(optarg); - break; - case 'b': - beta = atof(optarg); - break; - case 'B': - boundary_int = atoi(optarg); - switch (boundary_int) { - case 0: - boundary = FREE_BOUND; - boundc2 = 'f'; - break; - case 1: - boundary = CYLINDER_BOUND; - boundc2 = 'c'; - break; - default: - printf("boundary specifier must be 0 (FREE_BOUND) or 1 (CYLINDER_BOUND).\n"); - } - break; - default: /* '?' */ - exit(EXIT_FAILURE); - } - } - - char *break_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(break_filename, filename_len, "breaks_v_vc_%c_%u_%g.txt", boundc2, L, beta); - FILE *break_out = fopen(break_filename, "a"); - free(break_filename); - - - // start cholmod - cholmod_common c; - CHOL_F(start)(&c); - - (&c)->supernodal = CHOLMOD_SIMPLICIAL; - - - 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->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); - data_t *breaking_data_current = fracture_network(current_instance, fuse_thres, &c, cutoff); - free_instance(voltage_instance, &c); - free_instance(current_instance, &c); - free_instance(perm_voltage_instance, &c); - free_instance(perm_current_instance, &c); - free(fuse_thres); - - FILE *net_out = fopen("network.txt", "w"); - for (unsigned int j = 0; j < network->nv; j++) { - 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++) { - fprintf(net_out, "%u %u ", network->ev[2 * j], - network->ev[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (unsigned int j = 0; j < network->dnv; j++) { - 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++) { - fprintf(net_out, "%u %u ", network->dev[2 * j], - network->dev[2 * j + 1]); - } - - free_net(network, &c); - - for (unsigned int j = 0; j < breaking_data_voltage->num_broken; j++) { - fprintf(break_out, "%u %g %g ", breaking_data_voltage->break_list[j], - breaking_data_voltage->extern_field[j], breaking_data_voltage->conductivity[j]); - } - fprintf(break_out, "\n"); - for (unsigned int j = 0; j < breaking_data_current->num_broken; j++) { - fprintf(break_out, "%u %g %g ", breaking_data_current->break_list[j], - breaking_data_current->extern_field[j], breaking_data_current->conductivity[j]); - } - fprintf(break_out, "\n"); - - free_break_data(breaking_data_voltage); - free_break_data(breaking_data_current); - fclose(break_out); - - CHOL_F(finish)(&c); - - return 0; -} diff --git a/src/course_grain_square.c b/src/course_grain_square.c deleted file mode 100644 index 6587e83..0000000 --- a/src/course_grain_square.c +++ /dev/null @@ -1,168 +0,0 @@ - -#include "fracture.h" - -int main(int argc, char *argv[]) { - int opt; - - // defining variables to be (potentially) set by command line flags - int num = 100; - int width = 16; - double beta = .3; - bool save_clusters = false; - bool voltage_bound = false; - bool output_break_data = false; - - while ((opt = getopt(argc, argv, "n:w:b:cVo")) != -1) { - switch (opt) { - case 'n': - num = atoi(optarg); - break; - case 'w': - width = atoi(optarg); - break; - case 'b': - beta = atof(optarg); - break; - case 'c': - save_clusters = true; - break; - case 'V': - voltage_bound = true; - break; - case 'o': - output_break_data = true; - break; - default: /* '?' */ - exit(EXIT_FAILURE); - } - } - - FILE *break_out; - if (output_break_data) { - char *break_filename = (char *)malloc(100 * sizeof(char)); - snprintf(break_filename, 100, "breaks_%d_%.3f_%d.txt", width, beta, num); - break_out = fopen(break_filename, "w"); - free(break_filename); - } - - bool periodic = true; - double inf = 1; - double cutoff = 1e-10; - - // 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 (voltage_bound) { - (&c)->supernodal = CHOLMOD_SUPERNODAL; - } else { - (&c)->supernodal = CHOLMOD_SIMPLICIAL; - } - - graph_t *network = ini_square_network(width, periodic, false, &c); - graph_t *network_p = ini_square_network(width / 2, periodic, false, &c); - net_t *perm_instance = create_instance(network, inf, voltage_bound, true, &c); - unsigned int c_dist_size = network->dnv; - unsigned int c_p_dist_size = network_p->dnv; - - // define arrays for saving cluster and avalanche distributions - unsigned int *cluster_size_dist = - (unsigned int *)calloc(c_dist_size, sizeof(unsigned int)); - unsigned int *cluster_p_size_dist = - (unsigned int *)calloc(c_p_dist_size, sizeof(unsigned int)); - - printf("\n"); - for (int DUMB = 0; DUMB < num; DUMB++) { - printf("\033[F\033[JCOURSEGRAIN_SQUARE: %0*d / %d\n", (int)log10(num) + 1, - DUMB + 1, num); - - data_t *breaking_data = NULL; - net_t *instance = NULL; - while (breaking_data == NULL) { - double *fuse_thres = gen_fuse_thres( - 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); - free(fuse_thres); - } - - unsigned int min_pos = 0; - double min_val = DBL_MAX; - - for (unsigned int j = 0; j < breaking_data->num_broken; j++) { - double val = fabs(breaking_data->extern_field[j]); - if (val < min_val) { - min_pos = j; - min_val = val; - } - if (val > 10 * min_val) - break; - } - - net_t *tmp_instance = copy_instance(perm_instance, &c); - - for (unsigned int i = 0; i < breaking_data->num_broken; i++) { - break_edge(tmp_instance, breaking_data->break_list[i], &c); - } - - unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c); - for (unsigned int i = 0; i < network->dnv; i++) { - cluster_size_dist[i] += tmp_cluster_dist[i]; - } - free(tmp_cluster_dist); - - net_t *instance_p = coursegrain_square(tmp_instance, network_p, &c); - - unsigned int *tmp_cluster_p_dist = get_cluster_dist(instance_p, &c); - for (unsigned int i = 0; i < network_p->dnv; i++) { - cluster_p_size_dist[i] += tmp_cluster_p_dist[i]; - } - free(tmp_cluster_p_dist); - - free_instance(tmp_instance, &c); - free_instance(instance_p, &c); - if (output_break_data) { - for (unsigned int i = 0; i < breaking_data->num_broken; i++) { - fprintf(break_out, "%u %f ", breaking_data->break_list[i], - breaking_data->extern_field[i]); - } - fprintf(break_out, "\n"); - } - free(breaking_data->break_list); - free(breaking_data->extern_field); - free(breaking_data); - } - - printf("\033[F\033[JCURRENT_SCALING: COMPLETE"); - - if (save_clusters) { - FILE *cluster_out = get_file("cluster", width, 0, beta, 1, 1, num, false); - for (int i = 0; i < c_dist_size; i++) { - fprintf(cluster_out, "%u ", cluster_size_dist[i]); - } - fprintf(cluster_out, "\n"); - for (int i = 0; i < c_p_dist_size; i++) { - fprintf(cluster_out, "%u ", cluster_p_size_dist[i]); - } - fclose(cluster_out); - } - - if (output_break_data) { - fclose(break_out); - } - - free(cluster_size_dist); - free(cluster_p_size_dist); - free_instance(perm_instance, &c); - free_net(network, &c); - free_net(network_p, &c); - - CHOL_F(finish)(&c); - - return 0; -} diff --git a/src/coursegrain.c b/src/coursegrain.c deleted file mode 100644 index 3e9db1a..0000000 --- a/src/coursegrain.c +++ /dev/null @@ -1,40 +0,0 @@ - -#include "fracture.h" - -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, 1, - instance->voltage_bound, true, c); - - unsigned int width_p = width / 2; - bool *fuses = instance->fuses; - - for (unsigned int i = 0; i < network_p->ne; i++) { - int xp = i / width_p; - int yp = i % width_p; - unsigned int x1 = 2 * xp; - unsigned int y1 = (2 * yp - 1) % width; - unsigned int x2 = 2 * xp; - unsigned int y2 = 2 * yp; - unsigned int x3 = 2 * xp + 1; - unsigned int y3 = (2 * yp - 1) % width; - unsigned int x4 = 2 * xp + 1; - unsigned int y4 = 2 * yp; - bool f1 = fuses[width * x1 + y1]; - bool f2 = fuses[width * x2 + y2]; - bool f3 = fuses[width * x3 + y3]; - bool f4 = fuses[width * x4 + y4]; - - if ((f1 && f2) || (f3 && f4)) { - // instance_p->fuses[i] = true; - // instance_p->num_remaining_edges--; - break_edge(instance_p, i, c); - } - } - - // fin_instance(instance_p, c); - - return instance_p; -} diff --git a/src/current_scaling.c b/src/current_scaling.c deleted file mode 100644 index ef36b0e..0000000 --- a/src/current_scaling.c +++ /dev/null @@ -1,300 +0,0 @@ - -#include "fracture.h" - -int main(int argc, char *argv[]) { - int opt; - - // defining variables to be (potentially) set by command line flags - int iter = 1; - int num = 100; - int width = 16; - double crack_len = 8; - double beta = .3; - double inf = 1e-8; - double cutoff = 1e-8; - bool beta_shift = false; - bool supplied_bound = false; - bool ash_beta = false; - char *bound_file; - bool voltage_bound = false; - bool use_first = false; - bool save_stress = false; - bool save_bound = false; - bool save_damage = false; - bool save_strength = false; - - while ((opt = getopt(argc, argv, "n:w:b:l:i:Bf:aVFsSed")) != -1) { - switch (opt) { - case 'n': - num = atoi(optarg); - break; - case 'w': - width = atoi(optarg); - break; - case 'b': - beta = atof(optarg); - break; - case 'l': - crack_len = atof(optarg); - break; - case 'i': - iter = atoi(optarg); - break; - case 'B': - beta_shift = true; - break; - case 'a': - ash_beta = true; - break; - case 'V': - voltage_bound = true; - break; - case 'F': - use_first = true; - break; - case 's': - save_stress = true; - break; - case 'd': - save_damage = true; - break; - case 'e': - save_bound = true; - break; - case 'S': - save_strength = true; - break; - case 'f': - supplied_bound = true; - bound_file = optarg; - break; - default: /* '?' */ - exit(EXIT_FAILURE); - } - } - - // 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 (voltage_bound) { - (&c)->supernodal = CHOLMOD_SUPERNODAL; - } else { - (&c)->supernodal = CHOLMOD_SIMPLICIAL; - } - - 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, &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, &c); - finish_instance(tmp_instance, &c); - double *voltage = get_voltage(tmp_instance, &c); - - for (int i = 0; i < network->num_bounds; i++) { - for (int j = 0; j < network->bound_inds[i + 1] - network->bound_inds[i]; - j++) { - ((double *)perm_instance->boundary_cond - ->x)[network->bound_verts[network->bound_inds[i] + j]] = - voltage[network->bound_verts[network->bound_inds[i] + j]]; - } - } - (&c)->supernodal = CHOLMOD_SUPERNODAL; - } - - if (supplied_bound) { - FILE *bound_f = fopen(bound_file, "r"); - for (int i = 0; i < network->nv; i++) { - double tmp; - fscanf(bound_f, "%lg ", &tmp); - ((double *)perm_instance->boundary_cond->x)[i] = tmp; - } - - ((double *)perm_instance->boundary_cond->x)[network->nv] = 0; - ((double *)perm_instance->boundary_cond->x)[network->nv + 1] = 0; - fclose(bound_f); - } - - printf("\n"); - for (int DUMB2 = 0; DUMB2 < iter; DUMB2++) { - - double *strength; - if (save_strength) { - strength = (double *)malloc(num * sizeof(double)); - } - - double *damage; - if (save_damage) { - damage = (double *)calloc(network->ne, sizeof(double)); - } - double *avg_current = (double *)calloc(network->ne, sizeof(double)); - unsigned int *num_current_skipped = - (unsigned int *)calloc(network->ne, sizeof(unsigned int)); - double *avg_voltage = (double *)calloc(network->nv, sizeof(double)); - unsigned int *num_voltage_skipped = - (unsigned int *)calloc(network->nv, sizeof(unsigned int)); - - for (int DUMB = 0; DUMB < num; DUMB++) { - printf("\033[F\033[JCURRENT_SCALING: ITERATION %0*d: %0*d / %d\n", - (int)log10(iter) + 1, DUMB2 + 1, (int)log10(num) + 1, DUMB + 1, - num); - - data_t *breaking_data = NULL; - while (breaking_data == NULL) { - double *fuse_thres = gen_fuse_thres( - 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); - free(fuse_thres); - } - - unsigned int min_pos = 0; - double min_val = DBL_MAX; - - for (unsigned int j = 0; j < breaking_data->num_broken; j++) { - double val = fabs(breaking_data->extern_field[j]); - if (val < min_val) { - min_pos = j; - min_val = val; - } - } - - if (save_strength) { - strength[DUMB] = fabs(breaking_data->extern_field[min_pos]); - } - - net_t *tmp_instance = copy_instance(perm_instance, &c); - - unsigned int until = min_pos; - if (use_first) - until = 1; - for (unsigned int i = 0; i < until; i++) { - break_edge(tmp_instance, breaking_data->break_list[i], &c); - if (save_damage) { - damage[breaking_data->break_list[i]] += 1. / num; - } - } - - double *voltage = get_voltage(tmp_instance, &c); - double *current = get_current(tmp_instance, &c); - - for (unsigned int i = 0; i < network->ne; i++) { - avg_current[i] += current[i]; - if (current[i] == 0) - num_current_skipped[i]++; - } - - for (unsigned int i = 0; i < network->nv; i++) { - if (tmp_instance->marks[i] == tmp_instance->marks[network->nv]) { - avg_voltage[i] += voltage[i]; - } else { - num_voltage_skipped[i]++; - } - } - - free(current); - free(voltage); - free_instance(tmp_instance, &c); - free(breaking_data->break_list); - free(breaking_data->extern_field); - free(breaking_data); - } - - for (int i = 0; i < network->ne; i++) { - if (num_current_skipped[i] < num) { - avg_current[i] /= num - num_current_skipped[i]; - } - } - - for (int i = 0; i < network->nv; i++) { - if (num_voltage_skipped[i] < num) { - avg_voltage[i] /= num - num_voltage_skipped[i]; - } - } - - double *avg_field; - if (voltage_bound) - avg_field = avg_voltage; - else - avg_field = avg_current; - - update_boundary(perm_instance, avg_field); - - if (save_stress) { - char *c_filename = (char *)malloc(100 * sizeof(char)); - snprintf(c_filename, 100, "current_%d_%g_%d_%g.txt", width, crack_len, - iter, beta); - FILE *outfile = fopen(c_filename, "w"); - for (int i = 0; i < network->ne; i++) { - fprintf(outfile, "%g ", avg_current[i]); - } - fclose(outfile); - free(c_filename); - } - - if (save_damage) { - char *c_filename = (char *)malloc(100 * sizeof(char)); - snprintf(c_filename, 100, "damage_%d_%g_%d_%g.txt", width, crack_len, - iter, beta); - FILE *outfile = fopen(c_filename, "w"); - for (int i = 0; i < network->ne; i++) { - fprintf(outfile, "%g ", damage[i]); - } - fclose(outfile); - free(c_filename); - } - - if (save_strength) { - char *s_filename = (char *)malloc(100 * sizeof(char)); - snprintf(s_filename, 100, "strength_%d_%g_%d_%g.txt", width, crack_len, iter, beta); - FILE *f = fopen(s_filename, "a"); - for (int i = 0; i < num; i++) { - fprintf(f, "%g ", strength[i]); - } - fclose(f); - free(s_filename); - } - - if (save_bound) { - char *b_filename = (char *)malloc(100 * sizeof(char)); - snprintf(b_filename, 100, "bounds_%d_%g_%d_%g.txt", width, crack_len, - iter, beta); - FILE *outfile = fopen(b_filename, "w"); - for (int i = 0; i < network->nv; i++) { - fprintf(outfile, "%g ", ((double *)perm_instance->boundary_cond->x)[i]); - } - fclose(outfile); - free(b_filename); - } - - free(avg_current); - free(avg_voltage); - if (save_damage) free(damage); - free(num_current_skipped); - free(num_voltage_skipped); - if (save_strength) { - free(strength); - } - - printf( - "\033[F\033[JCURRENT_SCALING: ITERATION %0*d COMPLETE, BETA = %.2f\n\n", - (int)log10(iter) + 1, DUMB2 + 1, beta); - } - - free_instance(perm_instance, &c); - free_net(network, &c); - - CHOL_F(finish)(&c); - - return 0; -} diff --git a/src/get_file.c b/src/get_file.c deleted file mode 100644 index 9141631..0000000 --- a/src/get_file.c +++ /dev/null @@ -1,38 +0,0 @@ - -#include "fracture.h" - -FILE *get_file(const char *prefix, unsigned int width, unsigned int crack, - double beta, unsigned int iter, unsigned int num_iter, - unsigned int num, bool read) { - int prefix_len = strlen(prefix); - int width_len = 1 + (int)log10(width); - int crack_len; - if (crack != 0) { - crack_len = 1 + (int)log10(crack); - } else { - crack_len = 1; - } - int beta_len; - if (beta > 1) { - beta_len = 1 + (int)log10(beta) + 3; - } else { - beta_len = 4; - } - int iter_len = 1 + (int)log10(num_iter); - int num_len = 1 + (int)log10(num); - int num_underscores = 5; - - int len = prefix_len + width_len + crack_len + beta_len + iter_len + num_len + - num_underscores + 4; - - char filename[len + 1]; - snprintf(filename, sizeof(filename), "%s_%u_%u_%.2f_%0*u_%u.txt", prefix, - width, crack, beta, iter_len, iter + 1, num); - - char *mode = "w"; - if (read) { - mode = "r"; - } - FILE *file = fopen(filename, mode); - return file; -} diff --git a/src/homo_square_fracture.c b/src/homo_square_fracture.c deleted file mode 100644 index b301136..0000000 --- a/src/homo_square_fracture.c +++ /dev/null @@ -1,418 +0,0 @@ - - -#include "fracture.h" - -int main(int argc, char *argv[]) { - int opt; - - // defining variables to be (potentially) set by command line flags - unsigned int N, L, filename_len; - double beta, inf, cutoff; - bound_t boundary; - bool include_breaking, save_cluster_dist, use_voltage_boundaries, save_network, - save_crit_stress, save_toughness, save_corr, save_conductivity, - save_damage; - - filename_len = 100; - - N = 100; - L = 16; - beta = .3; - inf = 1e-8; - cutoff = 1e-10; - boundary = FREE_BOUND; - include_breaking = false; - save_cluster_dist = false; - use_voltage_boundaries = false; - save_network = false; - save_crit_stress = false; - save_damage = false; - save_corr = false; - save_conductivity = false; - save_toughness = false; - - int boundary_int; - char boundc2 = 'f'; - - while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrt")) != -1) { - switch (opt) { - case 'n': - N = atoi(optarg); - break; - case 'L': - L = atoi(optarg); - break; - case 'b': - beta = atof(optarg); - break; - case 'B': - boundary_int = atoi(optarg); - switch (boundary_int) { - 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; - default: - printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n"); - } - break; - case 'd': - save_damage = true; - break; - case 'V': - use_voltage_boundaries = 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 'C': - save_corr = 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_s_%c_%c_%u_%g.txt", boundc, boundc2, L, beta); - break_out = fopen(break_filename, "a"); - free(break_filename); - } - - // 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; - } - - - graph_t *network = ini_square_network(L, boundary, false, &c); - net_t *perm_instance = create_instance(network, inf, use_voltage_boundaries, true, &c); - unsigned int c_dist_size = network->dnv; - unsigned int a_dist_size = network->nv; - - // define arrays for saving cluster and avalanche distributions - unsigned int *cluster_size_dist; - unsigned int *avalanche_size_dist; - char *c_filename; - if (save_cluster_dist) { - cluster_size_dist = - (unsigned int *)calloc(c_dist_size, sizeof(unsigned int)); - avalanche_size_dist = - (unsigned int *)calloc(a_dist_size, sizeof(unsigned int)); - - c_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(c_filename, filename_len, "cstr_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); - - FILE *cluster_out = fopen(c_filename, "r"); - - if (cluster_out != NULL) { - for (unsigned int i = 0; i < c_dist_size; i++) { - unsigned int tmp; - fscanf(cluster_out, "%u ", &tmp); - cluster_size_dist[i] = tmp; - } - fscanf(cluster_out, "\n"); - for (unsigned int i = 0; i < a_dist_size; i++) { - unsigned int tmp; - fscanf(cluster_out, "%u ", &tmp); - avalanche_size_dist[i] = tmp; - } - fclose(cluster_out); - } - } - - double *crit_stress; - if (save_crit_stress) { - crit_stress = (double *)malloc(N * sizeof(double)); - } - - double *avg_corr; - unsigned int **dists = NULL; - if (save_corr) { - avg_corr = (double *)calloc(c_dist_size, sizeof(double)); - } - - double *conductivity; - if (save_conductivity) { - conductivity = (double *)malloc(N * sizeof(double)); - } - - double *damage; - if (save_damage) { - damage = (double *)malloc(N * sizeof(double)); - } - - double *toughness; - if (save_toughness) { - toughness = (double *)malloc(N * sizeof(double)); - } - - - printf("\n"); - 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->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); - free(fuse_thres); - - unsigned int max_pos = 0; - double max_val = 0; - - for (unsigned int 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]; - } - - net_t *tmp_instance = copy_instance(perm_instance, &c); - - unsigned int av_size = 0; - double cur_val = DBL_MAX; - for (unsigned int j = 0; j < max_pos; j++) { - break_edge(tmp_instance, breaking_data->break_list[j], &c); - - double val = fabs(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_conductivity) { - if (!use_voltage_boundaries) { - double *tmp_voltage = get_voltage(tmp_instance, &c); - conductivity[i] = 1/fabs(tmp_voltage[tmp_instance->graph->nv + 1] - tmp_voltage[tmp_instance->graph->nv]); - free(tmp_voltage); - } else { - conductivity[i] = breaking_data->conductivity[max_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 (unsigned int 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_damage) { - damage[i] = ((double)max_pos) / tmp_instance->graph->ne; - } - - if (save_cluster_dist) { - unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c); - for (unsigned int j = 0; j < tmp_instance->graph->dnv; j++) { - cluster_size_dist[j] += tmp_cluster_dist[j]; - } - free(tmp_cluster_dist); - } - - if (save_corr) { - double *tmp_corr = get_corr(tmp_instance, dists, &c); - for (unsigned int j = 0; j < tmp_instance->graph->dnv; j++) { - avg_corr[i] += tmp_corr[j] / N; - } - free(tmp_corr); - } - - 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->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++) { - fprintf(net_out, "%u %u ", tmp_instance->graph->ev[2 * j], - tmp_instance->graph->ev[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (unsigned int j = 0; j < tmp_instance->graph->dnv; j++) { - 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++) { - fprintf(net_out, "%u %u ", tmp_instance->graph->dev[2 * j], - tmp_instance->graph->dev[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (unsigned int j = 0; j < tmp_instance->graph->ne; j++) { - fprintf(net_out, "%d ", tmp_instance->fuses[j]); - } - fclose(net_out); - } - - free_instance(tmp_instance, &c); - - if (include_breaking) { - for (unsigned int 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"); - - free_instance(perm_instance, &c); - free_net(network, &c); - - if (save_cluster_dist) { - FILE *cluster_out = fopen(c_filename, "w"); - - for (int i = 0; i < c_dist_size; i++) { - fprintf(cluster_out, "%u ", cluster_size_dist[i]); - } - fprintf(cluster_out, "\n"); - for (int i = 0; i < a_dist_size; i++) { - fprintf(cluster_out, "%u ", avalanche_size_dist[i]); - } - fclose(cluster_out); - - free(c_filename); - free(cluster_size_dist); - free(avalanche_size_dist); - } - - if (save_corr) { - char *corr_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(corr_filename, filename_len, "corr_s_%c_%c_%d_%g.txt", boundc, boundc2, L, - beta); - FILE *corr_file = fopen(corr_filename, "w"); - for (unsigned int i = 0; i < c_dist_size; i++) { - fprintf(corr_file, "%g ", avg_corr[i]); - } - fclose(corr_file); - free(corr_filename); - free(avg_corr); - } - - if (save_conductivity) { - char *cond_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(cond_filename, filename_len, "cond_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); - FILE *cond_file = fopen(cond_filename, "a"); - for (unsigned int i = 0; i < N; i++) { - fprintf(cond_file, "%g ", conductivity[i]); - } - 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_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); - FILE *tough_file = fopen(tough_filename, "a"); - for (unsigned int i = 0; i < N; i++) { - fprintf(tough_file, "%g ", toughness[i]); - } - fclose(tough_file); - free(tough_filename); - free(toughness); - } - - if (save_damage) { - char *hdam_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(hdam_filename, filename_len, "damg_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); - FILE *hdam_file = fopen(hdam_filename, "a"); - for (unsigned int i = 0; i < N; i++) { - fprintf(hdam_file, "%g ", damage[i]); - } - fclose(hdam_file); - free(hdam_filename); - free(damage); - } - - if (include_breaking) { - fclose(break_out); - } - - if (save_crit_stress) { - char *a_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(a_filename, filename_len, "strs_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta); - FILE *a_file = fopen(a_filename, "a"); - for (int i = 0; i < N; i++) { - fprintf(a_file, "%g ", crit_stress[i]); - } - fclose(a_file); - free(a_filename); - free(crit_stress); - } - - CHOL_F(finish)(&c); - - return 0; -} diff --git a/src/homo_voronoi_fracture.c b/src/homo_voronoi_fracture.c deleted file mode 100644 index 85a9c2b..0000000 --- a/src/homo_voronoi_fracture.c +++ /dev/null @@ -1,405 +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; - bool include_breaking, save_cluster_dist, use_voltage_boundaries, use_dual, save_network, - save_crit_stress, save_toughness, save_conductivity, - save_damage; - bound_t boundary; - - filename_len = 100; - - N = 100; - L = 16; - 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_damage = false; - save_conductivity = false; - save_toughness = false; - - uint8_t bound_i; - char boundc2 = 'f'; - - while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrtD")) != -1) { - switch (opt) { - case 'n': - N = atoi(optarg); - break; - case 'L': - L = atoi(optarg); - break; - case 'b': - beta = 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; - 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 '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 '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.txt", boundc, boundc2, L, beta); - 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.dat", boundc, boundc2, L, beta); - snprintf(a_filename, filename_len, "avln_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta); - - 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 *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.dat", boundc, boundc2, L, beta); - - 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 *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->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); - 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]; - } - - 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_conductivity) { - if (!use_voltage_boundaries) { - double *tmp_voltage = get_voltage(tmp_instance, &c); - conductivity[i] = 1/fabs(tmp_voltage[tmp_instance->graph->nv + 1] - tmp_voltage[tmp_instance->graph->nv]); - free(tmp_voltage); - } else { - conductivity[i] = breaking_data->conductivity[max_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_damage) { - damage[max_pos]++; - } - - 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 < network->nv; j++) { - 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++) { - 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->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++) { - 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(network, &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_conductivity) { - char *cond_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(cond_filename, filename_len, "cond_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta); - 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.dat", boundc, boundc2, L, beta); - 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.dat", boundc, boundc2, L, beta); - 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/toy_setup.c b/src/toy_setup.c deleted file mode 100644 index 228acbc..0000000 --- a/src/toy_setup.c +++ /dev/null @@ -1,37 +0,0 @@ - -#include "fracture.h" - -bool break_fuses(unsigned int width, double *fuse_thres, double *field, - bool *fuses) { - assert(fuse_thres != NULL); - assert(field != NULL); - assert(fuses != NULL); - unsigned int size = pow(width, 2); - - for (unsigned int i = 0; i < size; i++) { - if (fuses[i]) { - if (1 / fuse_thres[i] > field[i]) { - fuses[i] = 0; - } - } - } - - return true; -} - -bool gen_toy_field(unsigned int width, double strength, double *field) { - assert(field != NULL); - unsigned int size = pow(width, 2); - - for (unsigned intint i = 0; i < size; i++) { - int x = i / width + 1; - int y = i % width + 1; - if (y < (width + 1) / 2.) - field[i] = fabs((width) / 2. - x) / strength; - else - field[i] = - sqrt(pow((width) / 2. - x, 2) + pow((width) / 2. - y, 2)) / strength; - } - - return true; -} diff --git a/src/update_beta.c b/src/update_beta.c deleted file mode 100644 index 4c1bf65..0000000 --- a/src/update_beta.c +++ /dev/null @@ -1,34 +0,0 @@ - -#include "fracture.h" - -double f(double damage) { - assert(damage <= 1 && damage >= 0); - return sqrt(1 - sqrt(1 - damage)); - // return 0.5 - 0.68182 * (0.5 - damage); -} - -double update_beta(double beta, unsigned int width, const double *stress, - const double *damage, double bound_total) { - - double total = 0; - unsigned int num_totaled = 0; - - for (unsigned int i = 0; i < pow(width, 2); i++) { - unsigned int stress_index = - width / 4 + (width / 4 + (i / width) / 2) * width + (i % width) / 2; - double outer_damage = f(pow(fabs(stress[i]), beta)); - double inner_stress = fabs(2 * stress[stress_index] / bound_total); - - if (outer_damage > 0 && inner_stress > 0 && inner_stress != 1) { - total += log(outer_damage) / log(inner_stress); - num_totaled++; - } - } - - assert(num_totaled > 0); - assert(total == total); - - double new_beta = total / num_totaled; - - return new_beta; -} diff --git a/src/update_boundary.c b/src/update_boundary.c deleted file mode 100644 index 6bc4fca..0000000 --- a/src/update_boundary.c +++ /dev/null @@ -1,86 +0,0 @@ - -#include "fracture.h" - -void update_boundary(net_t *instance, const double *avg_field) { - int size = instance->graph->ne; - int width = sqrt(size); - int num_verts = instance->graph->nv; - double *boundary = (double *)instance->boundary_cond->x; - - boundary[num_verts] = 0; - boundary[num_verts + 1] = 0; - - if (instance->voltage_bound) { - for (int i = 0; i < width / 4; i++) { - int t_ind = (width + 1) * (width / 8) + width / 4 - 1 + i; - double t_val = avg_field[t_ind]; - boundary[2 * i] = t_val; - boundary[2 * i + 1] = t_val; - - int b_ind = num_verts - (width + 1) * (width / 8) - width / 4 - i; - double b_val = avg_field[b_ind]; - boundary[num_verts - 1 - 2 * i] = b_val; - boundary[num_verts - 2 - 2 * i] = b_val; - - int l_ind = (width + 1) * (width / 8) + width / 2 + width / 4 - 1 + - i * (width + 1); - double l_val = avg_field[l_ind]; - boundary[width / 2 + 2 * i * (width + 1)] = l_val; - boundary[width / 2 + (2 * i + 1) * (width + 1)] = l_val; - - int r_ind = (width + 1) * (width / 8) + width - 1 + i * (width + 1); - double r_val = avg_field[r_ind]; - boundary[width + 2 * i * (width + 1)] = r_val; - boundary[width + (2 * i + 1) * (width + 1)] = r_val; - } - } else { - const double *stress = avg_field; - - for (int i = 0; i < width / 4; i++) { - boundary[2 * i] = stress[(width / 4 - 1) * width + width / 4 + 2 * i] + - stress[(width / 4 - 1) * width + width / 4 + 1 + 2 * i]; - boundary[2 * i + 1] = - stress[(width / 4 - 1) * width + width / 4 + 2 * i] + - stress[(width / 4 - 1) * width + width / 4 + 1 + 2 * i]; - - boundary[num_verts - 2 * i - 1] = - -stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 2 * i)] - - stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 1 + 2 * i)]; - boundary[num_verts - 2 * i - 2] = - -stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 2 * i)] - - stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 1 + 2 * i)]; - - boundary[(width + 1) / 2 + 2 * i * (width + 1)] = - stress[((width / 4) + 2 * i) * width + width / 4 - 1] - - stress[((width / 4) + 2 * i + 1) * width + width / 4 - 1]; - boundary[(width + 1) / 2 + (2 * i + 1) * (width + 1)] = - stress[((width / 4) + 2 * i) * width + width / 4 - 1] - - stress[((width / 4) + 2 * i + 1) * width + width / 4 - 1]; - - boundary[(width + 1) / 2 + 2 * i * (width + 1) + width / 2] = - stress[((width / 4) + 2 * i) * width + width / 4 + width / 2] - - stress[((width / 4) + 2 * i + 1) * width + width / 4 + width / 2]; - boundary[(width + 1) / 2 + (2 * i + 1) * (width + 1) + width / 2] = - stress[((width / 4) + 2 * i) * width + width / 4 + width / 2] - - stress[((width / 4) + 2 * i + 1) * width + width / 4 + width / 2]; - } - - double s_total = 0; - double total = 0; - for (int i = 0; i < num_verts; i++) { - if (boundary[i] != boundary[i]) { - boundary[i] = 0; - } else { - total += fabs(boundary[i]); - s_total += boundary[i]; - } - } - - // assert(fabs(s_total) < inf); - boundary[num_verts] -= s_total; - - for (int i = 0; i < num_verts + 1; i++) { - boundary[i] /= total; - } - } -} diff --git a/src/voro_fracture.c b/src/voro_fracture.c deleted file mode 100644 index 237192a..0000000 --- a/src/voro_fracture.c +++ /dev/null @@ -1,492 +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 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; -} |