diff options
Diffstat (limited to 'src/course_grain_square.c')
| -rw-r--r-- | src/course_grain_square.c | 168 |
1 files changed, 168 insertions, 0 deletions
diff --git a/src/course_grain_square.c b/src/course_grain_square.c new file mode 100644 index 0000000..9035cf7 --- /dev/null +++ b/src/course_grain_square.c @@ -0,0 +1,168 @@ + +#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; + } + + fnet *network = ini_square_network(width, periodic, false, &c); + fnet *network_p = ini_square_network(width / 2, periodic, false, &c); + finst *perm_instance = create_instance(network, inf, voltage_bound, true, &c); + unsigned int c_dist_size = network->num_dual_verts; + unsigned int c_p_dist_size = network_p->num_dual_verts; + + // 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); + + break_data *breaking_data = NULL; + finst *instance = NULL; + while (breaking_data == NULL) { + double *fuse_thres = gen_fuse_thres( + network->num_edges, network->edge_coords, 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; + } + + finst *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->num_dual_verts; i++) { + cluster_size_dist[i] += tmp_cluster_dist[i]; + } + free(tmp_cluster_dist); + + finst *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->num_dual_verts; 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_fnet(network, &c); + free_fnet(network_p, &c); + + CHOL_F(finish)(&c); + + return 0; +} |