diff options
Diffstat (limited to 'src/voro_fracture.c')
| -rw-r--r-- | src/voro_fracture.c | 492 |
1 files changed, 492 insertions, 0 deletions
diff --git a/src/voro_fracture.c b/src/voro_fracture.c new file mode 100644 index 0000000..237192a --- /dev/null +++ b/src/voro_fracture.c @@ -0,0 +1,492 @@ + + +#include "fracture.h" + +int main(int argc, char *argv[]) { + int opt; + + // defining variables to be (potentially) set by command line flags + uint8_t filename_len; + uint32_t N; + uint_t L; + double beta, inf, cutoff, crack_len; + bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network, + save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_conductivity, + save_damage, save_damage_field; + bound_t boundary; + + + // assume filenames will be less than 100 characters + + filename_len = 100; + + + // set default values + + N = 100; + L = 16; + crack_len = 0.; + beta = .3; + inf = 1e10; + cutoff = 1e-9; + boundary = FREE_BOUND; + save_data = false; + save_cluster_dist = false; + use_voltage_boundaries = false; + use_dual = false; + save_network = false; + save_crit_stress = false; + save_stress_field = false; + save_voltage_field = false; + save_damage = false; + save_damage_field = false; + save_conductivity = false; + save_toughness = false; + + + uint8_t bound_i; + char boundc2 = 'f'; + + + // get commandline options + + while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrtDSvel:")) != -1) { + switch (opt) { + case 'n': + N = atoi(optarg); + break; + case 'L': + L = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + case 'l': + crack_len = atof(optarg); + break; + case 'B': + bound_i = atoi(optarg); + switch (bound_i) { + case 0: + boundary = FREE_BOUND; + boundc2 = 'f'; + break; + case 1: + boundary = CYLINDER_BOUND; + boundc2 = 'c'; + break; + case 2: + boundary = TORUS_BOUND; + use_voltage_boundaries = true; + boundc2 = 't'; + break; + case 3: + boundary = EMBEDDED_BOUND; + boundc2 = 'e'; + use_dual = true; + use_voltage_boundaries = true; + break; + default: + printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n"); + exit(EXIT_FAILURE); + } + break; + case 'd': + save_damage = true; + break; + case 'e': + save_damage_field = true; + break; + case 'V': + use_voltage_boundaries = true; + break; + case 'D': + use_dual = true; + break; + case 'c': + save_cluster_dist = true; + break; + case 'o': + save_data = true; + break; + case 'N': + save_network = true; + break; + case 's': + save_crit_stress = true; + break; + case 'S': + save_stress_field = true; + break; + case 'v': + save_voltage_field = true; + break; + case 'r': + save_conductivity = true; + //inf = 1; + break; + case 't': + save_toughness = true; + break; + default: /* '?' */ + exit(EXIT_FAILURE); + } + } + + + char boundc; + if (use_voltage_boundaries) boundc = 'v'; + else boundc = 'c'; + + FILE *data_out; + if (save_data) { + char *data_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(data_filename, filename_len, "data_v_%c_%c_%u_%g_%g.txt", boundc, boundc2, L, beta, crack_len); + data_out = fopen(data_filename, "a"); + free(data_filename); + } + + uint_t voronoi_max_verts, c_dist_size, a_dist_size; + + voronoi_max_verts = 4 * pow(L, 2); + c_dist_size = voronoi_max_verts; + a_dist_size = voronoi_max_verts; + + if (voronoi_max_verts > CINT_MAX) { + exit(EXIT_FAILURE); + } + + // define arrays for saving cluster and avalanche distributions + uint32_t *cluster_size_dist; + uint32_t *avalanche_size_dist; + char *c_filename; + char *a_filename; + if (save_cluster_dist) { + cluster_size_dist = + (uint32_t *)malloc(c_dist_size * sizeof(uint32_t)); + avalanche_size_dist = + (uint32_t *)malloc(a_dist_size * sizeof(uint32_t)); + + c_filename = (char *)malloc(filename_len * sizeof(char)); + a_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(c_filename, filename_len, "cstr_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + snprintf(a_filename, filename_len, "avln_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + + FILE *cluster_out = fopen(c_filename, "rb"); + FILE *avalanche_out = fopen(a_filename, "rb"); + + if (cluster_out != NULL) { + fread(cluster_size_dist, sizeof(uint32_t), c_dist_size, cluster_out); + fclose(cluster_out); + } + if (avalanche_out != NULL) { + fread(avalanche_size_dist, sizeof(uint32_t), a_dist_size, avalanche_out); + fclose(avalanche_out); + } + } + + double *crit_stress; + if (save_crit_stress) { + crit_stress = (double *)malloc(N * sizeof(double)); + } + + double *stress_field; + unsigned int stress_pos = 0; + if (save_stress_field) { + stress_field = (double *)malloc(3 * N * voronoi_max_verts * sizeof(double)); + } + + double *voltage_field; + unsigned int voltage_pos = 0; + if (save_voltage_field) { + voltage_field = (double *)malloc(3 * N * voronoi_max_verts * sizeof(double)); + } + + double *damage_field; + unsigned int damage_pos = 0; + if (save_damage_field) { + damage_field = (double *)malloc(2 * N * voronoi_max_verts * sizeof(double)); + } + + double *conductivity; + if (save_conductivity) { + conductivity = (double *)malloc(N * sizeof(double)); + } + + // define arrays for saving damage distributions + uint32_t *damage; + char *d_filename; + if (save_damage) { + damage = + (uint32_t *)malloc(a_dist_size * sizeof(uint32_t)); + + d_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(d_filename, filename_len, "damg_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + + FILE *damage_out = fopen(d_filename, "rb"); + + if (damage_out != NULL) { + fread(damage, sizeof(uint32_t), a_dist_size, damage_out); + fclose(damage_out); + } + } + + double *toughness; + if (save_toughness) { + toughness = (double *)malloc(N * sizeof(double)); + } + + + // start cholmod + cholmod_common c; + CHOL_F(start)(&c); + + /* if we use voltage boundary conditions, the laplacian matrix is positive + * definite and we can use a supernodal LL decomposition. otherwise we need + * to use the simplicial LDL decomposition + */ + if (use_voltage_boundaries) { + //(&c)->supernodal = CHOLMOD_SUPERNODAL; + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + } else { + (&c)->supernodal = CHOLMOD_SIMPLICIAL; + } + + + printf("\n"); + for (uint32_t i = 0; i < N; i++) { + printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1, N); + + graph_t *g = ini_voro_graph(L, boundary, use_dual, genfunc_hyperuniform, &c); + net_t *net = net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c); + net_t *tmp_net = net_copy(net, &c); + data_t *data = net_fracture(tmp_net, &c, cutoff); + net_free(tmp_net, &c); + + uint_t max_pos = 0; + double max_val = 0; + + for (uint_t j = 0; j < data->num_broken; j++) { + double val = data->extern_field[j]; + + if (val > max_val) { + max_pos = j; + max_val = val; + } + } + + if (save_crit_stress) crit_stress[i] = data->extern_field[max_pos]; + + if (save_conductivity) conductivity[i] = data->conductivity[max_pos]; + + if (save_damage) damage[max_pos]++; + + uint_t av_size = 0; + double cur_val = 0; + for (uint_t j = 0; j < max_pos; j++) { + break_edge(net, data->break_list[j], &c); + + double val = data->extern_field[j]; + if (save_cluster_dist) { + if (val < cur_val) { + av_size++; + } + + if (val > cur_val) { + avalanche_size_dist[av_size]++; + av_size = 0; + cur_val = val; + } + } + } + + if (save_stress_field || save_voltage_field) { + double *tmp_voltages = get_voltage(net, &c); + if (save_voltage_field) { + for (uint_t j = 0; j < g->nv; j++) { + voltage_field[3 * voltage_pos] = g->vx[2 * j]; + voltage_field[3 * voltage_pos + 1] = g->vx[2 * j + 1]; + voltage_field[3 * voltage_pos + 2] = tmp_voltages[j]; + voltage_pos++; + } + } + if (save_stress_field) { + double *tmp_currents = get_current_v(net, tmp_voltages, &c); + for (uint_t j = 0; j < g->ne; j++) { + stress_field[3 * stress_pos] = g->ex[2 * j]; + stress_field[3 * stress_pos + 1] = g->ex[2 * j + 1]; + stress_field[3 * stress_pos + 2] = tmp_currents[j]; + stress_pos++; + } + free(tmp_currents); + } + free(tmp_voltages); + } + + if (save_damage_field) { + for (uint_t j = 0; j < max_pos; j++) { + damage_field[2 * damage_pos] = g->ex[2 * data->break_list[j]]; + damage_field[2 * damage_pos + 1] = g->ex[2 * data->break_list[j] + 1]; + damage_pos++; + } + } + + if (save_toughness) { + double tmp_toughness = 0; + if (max_pos > 0) { + double sigma1 = data->extern_field[0]; + double epsilon1 = sigma1 / data->conductivity[0]; + for (uint_t j = 0; j < max_pos - 1; j++) { + double sigma2 = data->extern_field[j+1]; + double epsilon2 = sigma2 / data->conductivity[j+1]; + if (epsilon2 > epsilon1) { + tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2; + sigma1 = sigma2; epsilon1 = epsilon2; + } + } + } + toughness[i] = tmp_toughness; + } + + if (save_cluster_dist) { + uint_t *tmp_cluster_dist = get_cluster_dist(net, &c); + for (uint_t j = 0; j < g->dnv; j++) { + cluster_size_dist[j] += tmp_cluster_dist[j]; + } + free(tmp_cluster_dist); + } + + if (save_network) { + FILE *net_out = fopen("network.txt", "w"); + for (uint_t j = 0; j < g->nv; j++) { + fprintf(net_out, "%f %f ", g->vx[2 * j], + g->vx[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (uint_t j = 0; j < g->ne; j++) { + fprintf(net_out, "%u %u ", g->ev[2 * j], g->ev[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (uint_t j = 0; j < g->dnv; j++) { + fprintf(net_out, "%f %f ", g->dvx[2 * j], + g->dvx[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (uint_t j = 0; j < g->ne; j++) { + fprintf(net_out, "%u %u ", g->dev[2 * j], g->dev[2 * j + 1]); + } + fprintf(net_out, "\n"); + for (uint_t j = 0; j < g->ne; j++) { + fprintf(net_out, "%d ", net->fuses[j]); + } + fclose(net_out); + } + + net_free(net, &c); + graph_free(g, &c); + + if (save_data) { + for (uint_t j = 0; j < data->num_broken; j++) { + fprintf(data_out, "%u %g %g ", data->break_list[j], + data->extern_field[j], data->conductivity[j]); + } + fprintf(data_out, "\n"); + } + + free_break_data(data); + } + + printf("\033[F\033[JFRACTURE: COMPLETE\n"); + + if (save_cluster_dist) { + FILE *cluster_out = fopen(c_filename, "wb"); + FILE *avalanche_out = fopen(a_filename, "wb"); + + fwrite(cluster_size_dist, sizeof(uint32_t), c_dist_size, cluster_out); + fwrite(avalanche_size_dist, sizeof(uint32_t), a_dist_size, avalanche_out); + + fclose(cluster_out); + fclose(avalanche_out); + + free(c_filename); + free(a_filename); + free(cluster_size_dist); + free(avalanche_size_dist); + } + + if (save_voltage_field) { + char *vfld_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(vfld_filename, filename_len, "vfld_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + FILE *vfld_file = fopen(vfld_filename, "ab"); + fwrite(voltage_field, sizeof(double), 3 * voltage_pos, vfld_file); + fclose(vfld_file); + free(vfld_filename); + free(voltage_field); + } + + if (save_stress_field) { + char *cfld_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(cfld_filename, filename_len, "cfld_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + FILE *cfld_file = fopen(cfld_filename, "ab"); + fwrite(stress_field, sizeof(double), 3 * stress_pos, cfld_file); + fclose(cfld_file); + free(cfld_filename); + free(stress_field); + } + + if (save_damage_field) { + char *dfld_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(dfld_filename, filename_len, "dfld_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + FILE *dfld_file = fopen(dfld_filename, "ab"); + fwrite(damage_field, sizeof(double), 2 * damage_pos, dfld_file); + fclose(dfld_file); + free(dfld_filename); + free(damage_field); + } + + if (save_conductivity) { + char *cond_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(cond_filename, filename_len, "cond_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + FILE *cond_file = fopen(cond_filename, "ab"); + fwrite(conductivity, sizeof(double), N, cond_file); + fclose(cond_file); + free(cond_filename); + free(conductivity); + } + + if (save_toughness) { + char *tough_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(tough_filename, filename_len, "tuff_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + FILE *tough_file = fopen(tough_filename, "ab"); + fwrite(toughness, sizeof(double), N, tough_file); + fclose(tough_file); + free(tough_filename); + free(toughness); + } + + if (save_damage) { + FILE *hdam_file = fopen(d_filename, "wb"); + fwrite(damage, sizeof(uint32_t), a_dist_size, hdam_file); + fclose(hdam_file); + free(d_filename); + free(damage); + } + + if (save_data) { + fclose(data_out); + } + + if (save_crit_stress) { + char *str_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(str_filename, filename_len, "strs_v_%c_%c_%d_%g_%g.dat", boundc, boundc2, L, beta, crack_len); + FILE *str_file = fopen(str_filename, "ab"); + fwrite(crit_stress, sizeof(double), N, str_file); + fclose(str_file); + free(str_filename); + free(crit_stress); + } + + CHOL_F(finish)(&c); + + return 0; +} |