diff options
Diffstat (limited to 'src/homo_square_fracture.c')
| -rw-r--r-- | src/homo_square_fracture.c | 418 |
1 files changed, 418 insertions, 0 deletions
diff --git a/src/homo_square_fracture.c b/src/homo_square_fracture.c new file mode 100644 index 0000000..a0dd305 --- /dev/null +++ b/src/homo_square_fracture.c @@ -0,0 +1,418 @@ + + +#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; + boundary_type 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; + } + + + fnet *network = ini_square_network(L, boundary, false, &c); + finst *perm_instance = create_instance(network, inf, use_voltage_boundaries, true, &c); + unsigned int c_dist_size = network->num_dual_verts; + unsigned int 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_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->num_edges, network->edge_coords, beta, beta_scaling_flat); + finst *instance = copy_instance(perm_instance, &c); + break_data *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]; + } + + finst *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->network->num_verts + 1] - tmp_voltage[tmp_instance->network->num_verts]); + 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->network->num_edges; + } + + if (save_cluster_dist) { + unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c); + for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; 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->network->num_dual_verts; 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->num_verts; j++) { + fprintf(net_out, "%f %f ", network->vert_coords[2 * j], + tmp_instance->network->vert_coords[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) { + fprintf(net_out, "%u %u ", tmp_instance->network->edges_to_verts[2 * j], + tmp_instance->network->edges_to_verts[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) { + fprintf(net_out, "%f %f ", tmp_instance->network->dual_vert_coords[2 * j], + tmp_instance->network->dual_vert_coords[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) { + fprintf(net_out, "%u %u ", tmp_instance->network->dual_edges_to_verts[2 * j], + tmp_instance->network->dual_edges_to_verts[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (unsigned int j = 0; j < tmp_instance->network->num_edges; 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_fnet(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; +} |