diff options
Diffstat (limited to 'src/fracture.c')
| -rw-r--r-- | src/fracture.c | 1017 |
1 files changed, 511 insertions, 506 deletions
diff --git a/src/fracture.c b/src/fracture.c index 73059ff..ede7a24 100644 --- a/src/fracture.c +++ b/src/fracture.c @@ -3,510 +3,515 @@ #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_energy, save_conductivity, - save_damage, save_threshold, save_current_load; - bound_t boundary; - lattice_t lattice; - - - // 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; - lattice = VORONOI_LATTICE; - save_data = 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_energy = false; - save_threshold = false; - save_current_load = false; - - - uint8_t bound_i; - char boundc2 = 'f'; - uint8_t lattice_i; - char lattice_c = 'v'; - char dual_c = 'o'; - - - // get commandline options - - while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TE")) != -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 'q': - lattice_i = atoi(optarg); - switch (lattice_i) { - case 0: - lattice = VORONOI_LATTICE; - lattice_c = 'v'; - break; - case 1: - lattice = DIAGONAL_LATTICE; - lattice_c = 'd'; - break; - case 2: - lattice = VORONOI_HYPERUNIFORM_LATTICE; - lattice_c = 'h'; - break; - case 3: - lattice = TRIANGULAR_LATTICE; - lattice_c = 't'; - break; - case 4: - lattice = SQUARE_LATTICE; - lattice_c = 's'; - break; - default: - printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 (DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n"); - exit(EXIT_FAILURE); - } - break; - case 'd': - save_damage = true; - break; - case 'V': - use_voltage_boundaries = true; - break; - case 'D': - use_dual = true; - dual_c = 'd'; - 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 'r': - save_conductivity = true; - break; - case 'E': - save_energy = true; - break; - case 'T': - save_threshold = true; - break; - case 'C': - save_current_load = 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_%c_%c_%c_%c_%u_%g_%g.txt", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - data_out = fopen(data_filename, "a"); - free(data_filename); - } - - uint_t max_verts, max_edges; - - // these are very liberal estimates - max_verts = 4 * pow(L, 2); - max_edges = 4 * pow(L, 2); - - if (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 *)calloc(max_verts, sizeof(uint32_t)); - avalanche_size_dist = - (uint32_t *)calloc(max_edges, 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_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, 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), max_verts, cluster_out); - fclose(cluster_out); - } - if (avalanche_out != NULL) { - fread(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out); - fclose(avalanche_out); - } - } - - long double *crit_stress; - if (save_crit_stress) { - crit_stress = (long double *)malloc(N * sizeof(long 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 *)calloc(max_edges, sizeof(uint32_t)); - - d_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - - FILE *damage_out = fopen(d_filename, "rb"); - - if (damage_out != NULL) { - fread(damage, sizeof(uint32_t), max_edges, damage_out); - fclose(damage_out); - } - } - - long double *energy; - if (save_energy) { - energy = (long double *)malloc(N * sizeof(long double)); - } - - long double *thresholds; - if (save_threshold) { - thresholds = (long double *)malloc(N * sizeof(long double)); - } - - long double *loads; - if (save_current_load) { - loads = (long double *)malloc(N * sizeof(long 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 = graph_create(lattice, boundary, L, use_dual); - 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; - long double max_val = 0; - - double cond0; - { - double *tmp_voltages = net_voltages(net, &c); - cond0 = net_conductivity(net, tmp_voltages); - free(tmp_voltages); - } - - for (uint_t j = 0; j < data->num_broken; j++) { - long double val = data->extern_field[j]; - - if (val > max_val) { - max_pos = j; - max_val = val; - } - } - - uint_t av_size = 0; - long double cur_val = 0; - - for (uint_t j = 0; j < max_pos; j++) { - uint_t next_broken = data->break_list[j]; - - break_edge(net, next_broken, &c); - - long 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_crit_stress) crit_stress[i] = data->extern_field[max_pos]; - - if (save_conductivity) { - if (max_pos > 0) { - conductivity[i] = data->conductivity[max_pos - 1]; - } else { - conductivity[i] = cond0; - } - } - - if (save_damage) { - uint_t would_break = 0; - double *tmp_voltage = net_voltages(net, &c); - double *tmp_current = net_currents(net, tmp_voltage, &c); - free(tmp_voltage); - for (uint_t j = 0; j < g->ne; j++) { - bool broken = net->fuses[j]; - bool under_thres = net->thres[j] < net->thres[data->break_list[max_pos]]; - bool zero_field = fabs(tmp_current[j]) < cutoff; - if (!broken && under_thres && zero_field) { - break_edge(net, j, &c); - } - } - damage[net->num_broken]++; - free(tmp_current); - } - - if (save_energy) { - long double tmp_energy = 0; - if (max_pos > 0) { - long double sigma1 = data->extern_field[0]; - double cond1 = cond0; - for (uint_t j = 0; j < max_pos - 1; j++) { - long double sigma2 = data->extern_field[j+1]; - double cond2 = data->conductivity[j]; - if (sigma2 > sigma1) { - tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1; - sigma1 = sigma2; - cond1 = cond2; - } - } - } - energy[i] = tmp_energy; - } - - if (save_threshold) { - thresholds[i] = net->thres[data->break_list[max_pos]]; - } - - if (save_current_load) { - loads[i] = data->extern_field[max_pos] / net->thres[data->break_list[max_pos]]; - } - - if (save_data) { - for (uint_t j = 0; j < data->num_broken; j++) { - fprintf(data_out, "%u %Lg %g ", data->break_list[j], - data->extern_field[j], data->conductivity[j]); - } - fprintf(data_out, "\n"); - } - - data_free(data); - 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); - } - - if (save_cluster_dist) { - uint_t *tmp_cluster_dist = get_cluster_dist(net); - for (uint_t j = 0; j < g->dnv; j++) { - cluster_size_dist[j] += tmp_cluster_dist[j]; - } - free(tmp_cluster_dist); - } - - - net_free(net, &c); - graph_free(g); - } - - 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), max_verts, cluster_out); - fwrite(avalanche_size_dist, sizeof(uint32_t), max_edges, 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_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, 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_energy) { - char *tough_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *tough_file = fopen(tough_filename, "ab"); - fwrite(energy, sizeof(long double), N, tough_file); - fclose(tough_file); - free(tough_filename); - free(energy); - } - - if (save_threshold) { - char *thres_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *thres_file = fopen(thres_filename, "ab"); - fwrite(thresholds, sizeof(long double), N, thres_file); - fclose(thres_file); - free(thres_filename); - free(thresholds); - } - - if (save_current_load) { - char *load_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(load_filename, filename_len, "load_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *load_file = fopen(load_filename, "ab"); - fwrite(loads, sizeof(long double), N, load_file); - fclose(load_file); - free(load_filename); - free(loads); - } - - if (save_damage) { - FILE *hdam_file = fopen(d_filename, "wb"); - fwrite(damage, sizeof(uint32_t), max_edges, 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_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *str_file = fopen(str_filename, "ab"); - fwrite(crit_stress, sizeof(long double), N, str_file); - fclose(str_file); - free(str_filename); - free(crit_stress); - } - - CHOL_F(finish)(&c); - - return 0; + 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_energy, save_conductivity, + save_damage, save_threshold, save_current_load; + bound_t boundary; + lattice_t lattice; + + // 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; + lattice = VORONOI_LATTICE; + save_data = 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_energy = false; + save_threshold = false; + save_current_load = false; + + uint8_t bound_i; + char boundc2 = 'f'; + uint8_t lattice_i; + char lattice_c = 'v'; + char dual_c = 'o'; + + // get commandline options + + while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TE")) != -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 'q': + lattice_i = atoi(optarg); + switch (lattice_i) { + case 0: + lattice = VORONOI_LATTICE; + lattice_c = 'v'; + break; + case 1: + lattice = DIAGONAL_LATTICE; + lattice_c = 'd'; + break; + case 2: + lattice = VORONOI_HYPERUNIFORM_LATTICE; + lattice_c = 'h'; + break; + case 3: + lattice = TRIANGULAR_LATTICE; + lattice_c = 't'; + break; + case 4: + lattice = SQUARE_LATTICE; + lattice_c = 's'; + break; + default: + printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 " + "(DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n"); + exit(EXIT_FAILURE); + } + break; + case 'd': + save_damage = true; + break; + case 'V': + use_voltage_boundaries = true; + break; + case 'D': + use_dual = true; + dual_c = 'd'; + 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 'r': + save_conductivity = true; + break; + case 'E': + save_energy = true; + break; + case 'T': + save_threshold = true; + break; + case 'C': + save_current_load = 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_%c_%c_%c_%c_%u_%g_%g.txt", + lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); + data_out = fopen(data_filename, "a"); + free(data_filename); + } + + uint_t max_verts, max_edges; + + // these are very liberal estimates + max_verts = 4 * pow(L, 2); + max_edges = 4 * pow(L, 2); + + if (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 *)calloc(max_verts, sizeof(uint32_t)); + avalanche_size_dist = (uint32_t *)calloc(max_edges, 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_%c_%c_%c_%c_%d_%g_%g.dat", + lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); + snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%d_%g_%g.dat", + lattice_c, dual_c, 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), max_verts, cluster_out); + fclose(cluster_out); + } + if (avalanche_out != NULL) { + fread(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out); + fclose(avalanche_out); + } + } + + long double *crit_stress; + if (save_crit_stress) { + crit_stress = (long double *)malloc(N * sizeof(long 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 *)calloc(max_edges, sizeof(uint32_t)); + + d_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%d_%g_%g.dat", + lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); + + FILE *damage_out = fopen(d_filename, "rb"); + + if (damage_out != NULL) { + fread(damage, sizeof(uint32_t), max_edges, damage_out); + fclose(damage_out); + } + } + + long double *energy; + if (save_energy) { + energy = (long double *)malloc(N * sizeof(long double)); + } + + long double *thresholds; + if (save_threshold) { + thresholds = (long double *)malloc(N * sizeof(long double)); + } + + long double *loads; + if (save_current_load) { + loads = (long double *)malloc(N * sizeof(long 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 = graph_create(lattice, boundary, L, use_dual); + 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; + long double max_val = 0; + + double cond0; + { + double *tmp_voltages = net_voltages(net, &c); + cond0 = net_conductivity(net, tmp_voltages); + free(tmp_voltages); + } + + for (uint_t j = 0; j < data->num_broken; j++) { + long double val = data->extern_field[j]; + + if (val > max_val) { + max_pos = j; + max_val = val; + } + } + + uint_t av_size = 0; + long double cur_val = 0; + + for (uint_t j = 0; j < max_pos; j++) { + uint_t next_broken = data->break_list[j]; + + break_edge(net, next_broken, &c); + + long 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_crit_stress) + crit_stress[i] = data->extern_field[max_pos]; + + if (save_conductivity) { + if (max_pos > 0) { + conductivity[i] = data->conductivity[max_pos - 1]; + } else { + conductivity[i] = cond0; + } + } + + if (save_damage) { + uint_t would_break = 0; + double *tmp_voltage = net_voltages(net, &c); + double *tmp_current = net_currents(net, tmp_voltage, &c); + free(tmp_voltage); + for (uint_t j = 0; j < g->ne; j++) { + bool broken = net->fuses[j]; + bool under_thres = + net->thres[j] < net->thres[data->break_list[max_pos]]; + bool zero_field = fabs(tmp_current[j]) < cutoff; + if (!broken && under_thres && zero_field) { + break_edge(net, j, &c); + } + } + damage[net->num_broken]++; + free(tmp_current); + } + + if (save_energy) { + long double tmp_energy = 0; + if (max_pos > 0) { + long double sigma1 = data->extern_field[0]; + double cond1 = cond0; + for (uint_t j = 0; j < max_pos - 1; j++) { + long double sigma2 = data->extern_field[j + 1]; + double cond2 = data->conductivity[j]; + if (sigma2 > sigma1) { + tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1; + sigma1 = sigma2; + cond1 = cond2; + } + } + } + energy[i] = tmp_energy; + } + + if (save_threshold) { + thresholds[i] = net->thres[data->break_list[max_pos]]; + } + + if (save_current_load) { + loads[i] = + data->extern_field[max_pos] / net->thres[data->break_list[max_pos]]; + } + + if (save_data) { + for (uint_t j = 0; j < data->num_broken; j++) { + fprintf(data_out, "%u %Lg %g ", data->break_list[j], + data->extern_field[j], data->conductivity[j]); + } + fprintf(data_out, "\n"); + } + + data_free(data); + 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); + } + + if (save_cluster_dist) { + uint_t *tmp_cluster_dist = get_cluster_dist(net); + for (uint_t j = 0; j < g->dnv; j++) { + cluster_size_dist[j] += tmp_cluster_dist[j]; + } + free(tmp_cluster_dist); + } + + net_free(net, &c); + graph_free(g); + } + + 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), max_verts, cluster_out); + fwrite(avalanche_size_dist, sizeof(uint32_t), max_edges, 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_%c_%c_%c_%c_%d_%g_%g.dat", + lattice_c, dual_c, 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_energy) { + char *tough_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat", + lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); + FILE *tough_file = fopen(tough_filename, "ab"); + fwrite(energy, sizeof(long double), N, tough_file); + fclose(tough_file); + free(tough_filename); + free(energy); + } + + if (save_threshold) { + char *thres_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat", + lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); + FILE *thres_file = fopen(thres_filename, "ab"); + fwrite(thresholds, sizeof(long double), N, thres_file); + fclose(thres_file); + free(thres_filename); + free(thresholds); + } + + if (save_current_load) { + char *load_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(load_filename, filename_len, "load_%c_%c_%c_%c_%d_%g_%g.dat", + lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); + FILE *load_file = fopen(load_filename, "ab"); + fwrite(loads, sizeof(long double), N, load_file); + fclose(load_file); + free(load_filename); + free(loads); + } + + if (save_damage) { + FILE *hdam_file = fopen(d_filename, "wb"); + fwrite(damage, sizeof(uint32_t), max_edges, 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_%c_%c_%c_%c_%d_%g_%g.dat", + lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); + FILE *str_file = fopen(str_filename, "ab"); + fwrite(crit_stress, sizeof(long double), N, str_file); + fclose(str_file); + free(str_filename); + free(crit_stress); + } + + CHOL_F(finish)(&c); + + return 0; } |