diff options
Diffstat (limited to 'src/fracture.c')
| -rw-r--r-- | src/fracture.c | 578 |
1 files changed, 578 insertions, 0 deletions
diff --git a/src/fracture.c b/src/fracture.c new file mode 100644 index 0000000..6798231 --- /dev/null +++ b/src/fracture.c @@ -0,0 +1,578 @@ + + +#include "fracture.h" + +int main(int argc, char *argv[]) { + int opt; + + // defining variables to be (potentially) set by command line flags + unsigned int num, width, filename_len; + double crack_len, crack_width, beta, inf, cutoff; + boundary_type periodic; + bool include_breaking, supplied_bound, save_clusters, voronoi, + side_bounds, voltage_bound, repeat_voronoi, network_out, save_avg_stress, save_avg_ztress, + save_app_stress, save_corr, save_cond, use_crit, use_first, save_damage, save_homo_damage; + char *bound_filename; + + filename_len = 100; + + num = 100; + width = 16; + crack_len = 16; + beta = .3; + inf = 1e-8; + cutoff = 1e-8; + periodic = FREE_BOUND; + include_breaking = false; + supplied_bound = false; + save_clusters = false; + voronoi = false; + side_bounds = false; + voltage_bound = false; + repeat_voronoi = false; + network_out = false; + save_avg_stress = false; + save_avg_ztress = false; + save_app_stress = false; + save_damage = false; + save_homo_damage = false; + save_corr = false; + save_cond = false; + use_crit = true; + use_first = false; + + int periodic_int; + + while ((opt = getopt(argc, argv, "n:w:b:l:f:ocp:vVsrNCaSFtzdH")) != -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 'p': + periodic_int = atoi(optarg); + switch (periodic_int) { + case 0: + periodic = FREE_BOUND; + break; + case 1: + periodic = CYLINDER_BOUND; + break; + case 2: + periodic = TORUS_BOUND; + break; + } + break; + case 'C': + save_avg_stress = true; + use_crit = true; + break; + case 'd': + save_damage = true; + break; + case 'H': + save_homo_damage = true; + break; + case 'z': + save_avg_ztress = true; + use_crit = true; + break; + case 'v': + voronoi = true; + break; + case 'F': + use_first = true; + break; + case 'r': + repeat_voronoi = true; + break; + case 'V': + voltage_bound = true; + break; + case 's': + side_bounds = true; + break; + case 'c': + save_clusters = true; + break; + case 'o': + include_breaking = true; + break; + case 'N': + network_out = true; + break; + case 'a': + save_app_stress = true; + break; + case 'S': + save_corr = true; + break; + case 't': + save_cond = true; + inf = 1; + break; + case 'f': + supplied_bound = true; + bound_filename = optarg; + break; + default: /* '?' */ + exit(EXIT_FAILURE); + } + } + + FILE *break_out; + if (include_breaking) { + char *break_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(break_filename, filename_len, "breaks_%u_%g_%g_%u.txt", width, + crack_len, beta, num); + break_out = fopen(break_filename, "w"); + 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 (voltage_bound) { + (&c)->supernodal = CHOLMOD_SUPERNODAL; + } else { + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + } + + fnet *network; + finst *perm_instance; + unsigned int c_dist_size; + unsigned int a_dist_size; + unsigned int voronoi_num = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2); + + if (voronoi) { + crack_width = 1. / (2 * width); + c_dist_size = 2 * voronoi_num; + a_dist_size = 2 * voronoi_num; + if (repeat_voronoi) { + while ((network = ini_voronoi_network(width, periodic, + genfunc_uniform, &c)) == NULL) + ; + perm_instance = create_instance(network, inf, voltage_bound, false, &c); + gen_crack(perm_instance, crack_len, crack_width, &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); + finish_instance(perm_instance, &c); + c_dist_size = network->num_dual_verts; + a_dist_size = network->num_verts; + } + + // define arrays for saving cluster and avalanche distributions + unsigned int *cluster_size_dist; + unsigned int *avalanche_size_dist; + char *c_filename; + if (save_clusters) { + 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, "cluster_%d_%g_%g.txt", width, crack_len, + 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 *app_stress; + double *app_stress_crit; + if (save_app_stress) { + app_stress = (double *)malloc(num * sizeof(double)); + app_stress_crit = (double *)malloc(num * sizeof(double)); + } + + double *avg_corr; + unsigned int **dists = NULL; + if (save_corr) { + avg_corr = (double *)calloc(2 * voronoi_num, sizeof(double)); + if (!voronoi) + dists = get_dists(network); + } + + double *conductivity; + if (save_cond) { + conductivity = (double *)malloc(num * sizeof(double)); + } + + double *avg_stress; + unsigned int *num_stress; + if (save_avg_stress) { + avg_stress = (double *)calloc(pow(width, 2), sizeof(double)); + num_stress = (unsigned int *)calloc(pow(width, 2), sizeof(unsigned int)); + } + double *avg_ztress; + if (save_avg_ztress) { + avg_ztress = (double *)calloc(pow(width, 2), sizeof(double)); + } + + double *damage; + if (save_damage) { + damage = (double *)calloc(pow(width, 2), sizeof(double)); + } + + double *homo_damage; + if (save_homo_damage) { + homo_damage = (double *)calloc(num, sizeof(double)); + } + + double *square_bound; + unsigned int square_num_verts = 2 * ((width + 1) / 2) * (width / 2 + 1); + if (supplied_bound) { + FILE *bound_file = fopen(bound_filename, "r"); + square_bound = (double *)calloc(square_num_verts, sizeof(double)); + for (unsigned int i = 0; i < square_num_verts; i++) { + double tmp; + fscanf(bound_file, "%lg ", &tmp); + square_bound[i] = tmp; + } + fclose(bound_file); + if (!voronoi) { + double total = 0; + for (unsigned int i = 0; i < square_num_verts; i++) { + ((double *)perm_instance->boundary_cond->x)[i] = square_bound[i]; + total += square_bound[i]; + } + ((double *)perm_instance->boundary_cond->x)[square_num_verts] = 0; + ((double *)perm_instance->boundary_cond->x)[square_num_verts + 1] = 0; + } + } + + printf("\n"); + for (int DUMB = 0; DUMB < num; DUMB++) { + printf("\033[F\033[JFRACTURE: %0*d / %d\n", (int)log10(num) + 1, DUMB + 1, + num); + + break_data *breaking_data = NULL; + while (breaking_data == NULL) { + if (voronoi && !repeat_voronoi) { + while ((network = ini_voronoi_network(width, periodic, + genfunc_uniform, &c)) == NULL) + ; + perm_instance = create_instance(network, inf, voltage_bound, false, &c); + gen_crack(perm_instance, crack_len, crack_width, &c); + finish_instance(perm_instance, &c); + if (supplied_bound) { + voronoi_bound_ini(perm_instance, square_bound, width); + } + } + double *fuse_thres = gen_fuse_thres( + network->num_edges, network->edge_coords, beta, beta_scaling_flat); + finst *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_app_stress) { + app_stress[DUMB] = + 1 / fabs(breaking_data->extern_field[breaking_data->num_broken - 1]); + app_stress_crit[DUMB] = 1 / fabs(breaking_data->extern_field[min_pos]); + } + + finst *tmp_instance = copy_instance(perm_instance, &c); + + int stop_at = breaking_data->num_broken; + if (use_crit) + stop_at = min_pos; + if (use_first) + stop_at = 0; + + unsigned int av_size = 0; + double cur_val = DBL_MAX; + for (unsigned int i = 0; i < min_pos; i++) { + double val = fabs(breaking_data->extern_field[i]); + if (save_clusters) { + if (val > cur_val) { + av_size++; + } + + if (val < cur_val) { + avalanche_size_dist[av_size]++; + av_size = 0; + cur_val = val; + } + } + } + + for (unsigned int i = 0; i < stop_at; i++) { + break_edge(tmp_instance, breaking_data->break_list[i], &c); + } + + if (save_cond) { + double *tmp_voltage = get_voltage(tmp_instance, &c); + conductivity[DUMB] = fabs(tmp_voltage[tmp_instance->network->num_verts + 1] - tmp_voltage[tmp_instance->network->num_verts]); + free(tmp_voltage); + } + + if (save_homo_damage) { + homo_damage[DUMB] = ((double)min_pos) / tmp_instance->network->num_edges; + } + + if (save_avg_stress || save_avg_ztress) { + double *tmp_stress = get_current(tmp_instance, &c); + if (voronoi) { + double *tmp_stress_2 = + bin_values(tmp_instance->network, width, tmp_stress); + free(tmp_stress); + tmp_stress = tmp_stress_2; + } + for (unsigned int i = 0; i < pow(width, 2); i++) { + double sigma = tmp_stress[i]; + if (sigma != 0 && sigma == sigma && save_avg_stress) { + avg_stress[i] += sigma; + num_stress[i]++; + } + if (sigma == sigma && save_avg_ztress) { + avg_ztress[i] += sigma; + } + } + free(tmp_stress); + } + + if (save_damage) { + double *tmp_damage = (double *)calloc(tmp_instance->network->num_edges, sizeof(double)); + for (unsigned int i = 0; i < stop_at; i++) { + tmp_damage[breaking_data->break_list[i]] += 1; + } + if (voronoi) { + double *tmp_damage_2 = bin_values(tmp_instance->network, width, tmp_damage); + free(tmp_damage); + tmp_damage = tmp_damage_2; + } + for (unsigned int i = 0; i < pow(width, 2); i++) { + damage[i] += tmp_damage[i]; + } + free(tmp_damage); + } + + if (save_clusters) { + unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c); + for (unsigned int i = 0; i < network->num_dual_verts; i++) { + cluster_size_dist[i] += tmp_cluster_dist[i]; + } + free(tmp_cluster_dist); + } + + if (save_corr) { + double *tmp_corr = get_corr(tmp_instance, dists, &c); + for (unsigned int i = 0; i < tmp_instance->network->num_dual_verts; i++) { + avg_corr[i] += tmp_corr[i] / num; + } + free(tmp_corr); + } + + if (network_out) { + FILE *net_out = fopen("network.txt", "w"); + for (unsigned int i = 0; i < network->num_verts; i++) { + fprintf(net_out, "%f %f ", network->vert_coords[2 * i], + network->vert_coords[2 * i + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int i = 0; i < network->num_edges; i++) { + fprintf(net_out, "%u %u ", network->edges_to_verts[2 * i], + network->edges_to_verts[2 * i + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int i = 0; i < network->num_dual_verts; i++) { + fprintf(net_out, "%f %f ", network->dual_vert_coords[2 * i], + network->dual_vert_coords[2 * i + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int i = 0; i < network->num_edges; i++) { + fprintf(net_out, "%u %u ", network->dual_edges_to_verts[2 * i], + network->dual_edges_to_verts[2 * i + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int i = 0; i < network->num_edges; i++) { + fprintf(net_out, "%d ", tmp_instance->fuses[i]); + } + fclose(net_out); + } + + free_instance(tmp_instance, &c); + if (voronoi && !repeat_voronoi) { + free_fnet(network, &c); + free_instance(perm_instance, &c); + } + if (include_breaking) { + 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_break_data(breaking_data); + } + + printf("\033[F\033[JFRACTURE: COMPLETE\n"); + + if (save_clusters) { + 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(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_%d_%g_%g.txt", width, crack_len, + beta); + FILE *corr_file = fopen(corr_filename, "w"); + for (unsigned int i = 0; i < 2 * voronoi_num; i++) { + fprintf(corr_file, "%g ", avg_corr[i]); + } + fclose(corr_file); + free(corr_filename); + } + + if (save_cond) { + char *cond_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(cond_filename, filename_len, "conductivity_%d_%g_%g.txt", width, crack_len, beta); + FILE *cond_file = fopen(cond_filename, "a"); + for (unsigned int i = 0; i < num; i++) { + fprintf(cond_file, "%g ", conductivity[i]); + } + fclose(cond_file); + free(cond_filename); + free(conductivity); + } + + if (save_homo_damage) { + char *hdam_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(hdam_filename, filename_len, "damagep_%d_%g.txt", width, beta); + FILE *hdam_file = fopen(hdam_filename, "a"); + for (unsigned int i = 0; i < num; i++) { + fprintf(hdam_file, "%g ", homo_damage[i]); + } + fclose(hdam_file); + free(hdam_filename); + free(homo_damage); + } + + if (!voronoi || repeat_voronoi) { + free_instance(perm_instance, &c); + //free_fnet(network, &c); + } + + if (include_breaking) { + fclose(break_out); + } + + if (save_avg_stress) { + char *stress_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(stress_filename, filename_len, "current_%d_%g_%g_%u.txt", width, + crack_len, beta, num); + FILE *stress_file = fopen(stress_filename, "w"); + for (unsigned int i = 0; i < pow(width, 2); i++) { + if (num_stress[i] != 0) + avg_stress[i] /= num_stress[i]; + fprintf(stress_file, "%g ", avg_stress[i]); + } + fclose(stress_file); + free(stress_filename); + free(avg_stress); + } + if (save_avg_ztress) { + char *ztress_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(ztress_filename, filename_len, "zurrent_%d_%g_%g_%u.txt", width, + crack_len, beta, num); + FILE *stress_file = fopen(ztress_filename, "w"); + for (unsigned int i = 0; i < pow(width, 2); i++) { + fprintf(stress_file, "%g ", avg_ztress[i] / num); + } + fclose(stress_file); + free(ztress_filename); + free(avg_ztress); + } + + if (save_app_stress) { + char *a_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(a_filename, filename_len, "astress_%d_%g_%g.txt", width, crack_len, + beta); + FILE *a_file = fopen(a_filename, "a"); + for (int i = 0; i < num; i++) { + fprintf(a_file, "%g %g\n", app_stress[i], app_stress_crit[i]); + } + fclose(a_file); + free(a_filename); + free(app_stress); + free(app_stress_crit); + } + + if (save_damage) { + char *damage_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(damage_filename, filename_len, "damage_%d_%g_%g_%u.txt", width, + crack_len, beta, num); + FILE *damage_file = fopen(damage_filename, "w"); + for (unsigned int i = 0; i < pow(width, 2); i++) { + fprintf(damage_file, "%g ", damage[i] / num); + } + fclose(damage_file); + free(damage_filename); + free(damage); + + } + + CHOL_F(finish)(&c); + + return 0; +} |