diff options
| -rw-r--r-- | makefile | 4 | ||||
| -rw-r--r-- | makefile_hal | 4 | ||||
| -rw-r--r-- | src/fracture.c | 819 | ||||
| -rw-r--r-- | src/fracture.h | 6 | ||||
| -rw-r--r-- | src/graph_create.c (renamed from src/ini_network.c) | 14 |
5 files changed, 400 insertions, 447 deletions
@@ -3,8 +3,8 @@ CC = clang CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -march=native #-fopenmp LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler -ltcmalloc -OBJ = break_data bound_set bin_values correlations beta_scales randfuncs net get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry net_fracture get_current update_factor update_boundary get_file update_beta gen_voltcurmat ini_network free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity net_notch -BIN = corr_test voro_fracture +OBJ = break_data bound_set bin_values correlations beta_scales randfuncs net get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry net_fracture get_current update_factor update_boundary get_file update_beta gen_voltcurmat graph_create free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity net_notch +BIN = corr_test voro_fracture fracture all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} diff --git a/makefile_hal b/makefile_hal index f2bc954..be9bccc 100644 --- a/makefile_hal +++ b/makefile_hal @@ -3,8 +3,8 @@ CC = clang CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto #-fopenmp LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler -ltcmalloc -OBJ = break_data bound_set bin_values correlations beta_scales randfuncs net get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry net_fracture get_current update_factor update_boundary get_file update_beta gen_voltcurmat ini_network free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity net_notch -BIN = corr_test voro_fracture +OBJ = break_data bound_set bin_values correlations beta_scales randfuncs net get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry net_fracture get_current update_factor update_boundary get_file update_beta gen_voltcurmat graph_create free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity net_notch +BIN = corr_test voro_fracture fracture all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} diff --git a/src/fracture.c b/src/fracture.c index 6053146..e7abaec 100644 --- a/src/fracture.c +++ b/src/fracture.c @@ -6,337 +6,312 @@ 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; + 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; + lattice_t lattice; + + + // assume filenames will be less than 100 characters filename_len = 100; - num = 100; - width = 16; - crack_len = 16; + + // set default values + + N = 100; + L = 16; + crack_len = 0.; 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; + 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_stress_field = false; + save_voltage_field = false; save_damage = false; - save_homo_damage = false; - save_corr = false; - save_cond = false; - use_crit = true; - use_first = false; + save_damage_field = false; + save_conductivity = false; + save_toughness = false; + + + uint8_t bound_i; + char boundc2 = 'f'; + uint8_t lattice_i; + char lattice_c = 'v'; + - int periodic_int; + // get commandline options - while ((opt = getopt(argc, argv, "n:w:b:l:f:ocp:vVsrNCaSFtzdH")) != -1) { + while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:")) != -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); + 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 = SQUARE_LATTICE; + lattice_c = 's'; + break; + default: + printf("lattice specifier must be 0 (VORONOI_LATTICE) or 1 (SQUARE_LATTICE).\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; + break; + case 't': + save_toughness = true; + 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); + char boundc; + if (use_voltage_boundaries) boundc = 'v'; + else boundc = '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; + FILE *data_out; + if (save_data) { + char *data_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(data_filename, filename_len, "data_%c_%c_%c_%u_%g_%g.txt", lattice_c, boundc, boundc2, L, beta, crack_len); + data_out = fopen(data_filename, "a"); + free(data_filename); } - graph_t *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, &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, &c); - finish_instance(perm_instance, &c); - c_dist_size = network->dnv; - a_dist_size = network->nv; + 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 - unsigned int *cluster_size_dist; - unsigned int *avalanche_size_dist; + uint32_t *cluster_size_dist; + uint32_t *avalanche_size_dist; char *c_filename; - if (save_clusters) { + char *a_filename; + if (save_cluster_dist) { cluster_size_dist = - (unsigned int *)calloc(c_dist_size, sizeof(unsigned int)); + (uint32_t *)malloc(max_verts * sizeof(uint32_t)); avalanche_size_dist = - (unsigned int *)calloc(a_dist_size, sizeof(unsigned int)); + (uint32_t *)malloc(max_edges * sizeof(uint32_t)); c_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(c_filename, filename_len, "cluster_%d_%g_%g.txt", width, crack_len, - beta); + a_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(c_filename, filename_len, "cstr_%c_%c_%c_%d_%g_%g.dat", lattice_c, boundc, boundc2, L, beta, crack_len); + snprintf(a_filename, filename_len, "avln_%c_%c_%c_%d_%g_%g.dat", lattice_c, boundc, boundc2, L, beta, crack_len); - FILE *cluster_out = fopen(c_filename, "r"); + FILE *cluster_out = fopen(c_filename, "rb"); + FILE *avalanche_out = fopen(a_filename, "rb"); 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; - } + 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); + } } - 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 *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(2 * voronoi_num, sizeof(double)); - if (!voronoi) - dists = get_dists(network); + double *stress_field; + unsigned int stress_pos = 0; + if (save_stress_field) { + stress_field = (double *)malloc(3 * N * max_verts * sizeof(double)); } - double *conductivity; - if (save_cond) { - conductivity = (double *)malloc(num * sizeof(double)); + double *voltage_field; + unsigned int voltage_pos = 0; + if (save_voltage_field) { + voltage_field = (double *)malloc(3 * N * max_verts * 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 *damage_field; + unsigned int damage_pos = 0; + if (save_damage_field) { + damage_field = (double *)malloc(2 * N * max_verts * sizeof(double)); } - double *avg_ztress; - if (save_avg_ztress) { - avg_ztress = (double *)calloc(pow(width, 2), sizeof(double)); + + double *conductivity; + if (save_conductivity) { + conductivity = (double *)malloc(N * sizeof(double)); } - double *damage; + // define arrays for saving damage distributions + uint32_t *damage; + char *d_filename; if (save_damage) { - damage = (double *)calloc(pow(width, 2), sizeof(double)); + damage = + (uint32_t *)malloc(max_edges * sizeof(uint32_t)); + + d_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(d_filename, filename_len, "damg_%c_%c_%c_%d_%g_%g.dat", lattice_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); + } } - double *homo_damage; - if (save_homo_damage) { - homo_damage = (double *)calloc(num, sizeof(double)); + double *toughness; + if (save_toughness) { + toughness = (double *)malloc(N * 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; - } + + // 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 (int DUMB = 0; DUMB < num; DUMB++) { - printf("\033[F\033[JFRACTURE: %0*d / %d\n", (int)log10(num) + 1, DUMB + 1, - num); - - data_t *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, &c); - finish_instance(perm_instance, &c); - if (supplied_bound) { - voronoi_bound_ini(perm_instance, square_bound, width); - } - } - double *fuse_thres = gen_fuse_thres( - network->ne, network->ex, 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); - } + 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); - unsigned int min_pos = 0; - double min_val = DBL_MAX; + graph_t *g = graph_create(lattice, boundary, L, use_dual, &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); - for (unsigned int j = 0; j < breaking_data->num_broken; j++) { - double val = fabs(breaking_data->extern_field[j]); + 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 < min_val) { - min_pos = j; - min_val = val; + if (val > max_val) { + max_pos = j; + max_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]); - } + if (save_crit_stress) crit_stress[i] = data->extern_field[max_pos]; - finst *tmp_instance = copy_instance(perm_instance, &c); + if (save_conductivity) conductivity[i] = data->conductivity[max_pos]; - int stop_at = breaking_data->num_broken; - if (use_crit) - stop_at = min_pos; - if (use_first) - stop_at = 0; + if (save_damage) damage[max_pos]++; - 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) { + 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) { + if (val > cur_val) { avalanche_size_dist[av_size]++; av_size = 0; cur_val = val; @@ -344,232 +319,190 @@ int main(int argc, char *argv[]) { } } - 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->graph->nv + 1] - tmp_voltage[tmp_instance->graph->nv]); - free(tmp_voltage); - } - - if (save_homo_damage) { - homo_damage[DUMB] = ((double)min_pos) / tmp_instance->graph->ne; - } - - 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->graph, 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 (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 (sigma == sigma && save_avg_ztress) { - avg_ztress[i] += sigma; + } + 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_stress); + free(tmp_voltages); } - if (save_damage) { - double *tmp_damage = (double *)calloc(tmp_instance->graph->ne, 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->graph, 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]; + 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++; } - free(tmp_damage); } - if (save_clusters) { - unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c); - for (unsigned int i = 0; i < network->dnv; i++) { - cluster_size_dist[i] += tmp_cluster_dist[i]; + 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; + } + } } - free(tmp_cluster_dist); + toughness[i] = tmp_toughness; } - if (save_corr) { - double *tmp_corr = get_corr(tmp_instance, dists, &c); - for (unsigned int i = 0; i < tmp_instance->graph->dnv; i++) { - avg_corr[i] += tmp_corr[i] / num; + 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_corr); + free(tmp_cluster_dist); } - if (network_out) { + if (save_network) { FILE *net_out = fopen("network.txt", "w"); - for (unsigned int i = 0; i < network->nv; i++) { - fprintf(net_out, "%f %f ", network->vx[2 * i], - network->vx[2 * i + 1]); + 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 (unsigned int i = 0; i < network->ne; i++) { - fprintf(net_out, "%u %u ", network->ev[2 * i], - network->ev[2 * i + 1]); + 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 (unsigned int i = 0; i < network->dnv; i++) { - fprintf(net_out, "%f %f ", network->dvx[2 * i], - network->dvx[2 * i + 1]); + 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 (unsigned int i = 0; i < network->ne; i++) { - fprintf(net_out, "%u %u ", network->dev[2 * i], - network->dev[2 * i + 1]); + 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 (unsigned int i = 0; i < network->ne; i++) { - fprintf(net_out, "%d ", tmp_instance->fuses[i]); + for (uint_t j = 0; j < g->ne; j++) { + fprintf(net_out, "%d ", net->fuses[j]); } fclose(net_out); } - free_instance(tmp_instance, &c); - if (voronoi && !repeat_voronoi) { - free_net(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]); + 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(break_out, "\n"); + fprintf(data_out, "\n"); } - free_break_data(breaking_data); + + free_break_data(data); } printf("\033[F\033[JFRACTURE: COMPLETE\n"); - if (save_clusters) { - FILE *cluster_out = fopen(c_filename, "w"); + 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); - 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); + fclose(avalanche_out); + free(c_filename); + free(a_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_%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_voltage_field) { + char *vfld_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(vfld_filename, filename_len, "vfld_%c_%c_%c_%d_%g_%g.dat", lattice_c, 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_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_stress_field) { + char *cfld_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(cfld_filename, filename_len, "cfld_%c_%c_%c_%d_%g_%g.dat", lattice_c, 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_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 (save_damage_field) { + char *dfld_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(dfld_filename, filename_len, "dfld_%c_%c_%c_%d_%g_%g.dat", lattice_c, 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 (!voronoi || repeat_voronoi) { - free_instance(perm_instance, &c); - //free_net(network, &c); + if (save_conductivity) { + char *cond_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(cond_filename, filename_len, "cond_%c_%c_%c_%d_%g_%g.dat", lattice_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 (include_breaking) { - fclose(break_out); + if (save_toughness) { + char *tough_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(tough_filename, filename_len, "tuff_%c_%c_%c_%d_%g_%g.dat", lattice_c, 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_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_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_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_data) { + fclose(data_out); } - 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); - + if (save_crit_stress) { + char *str_filename = (char *)malloc(filename_len * sizeof(char)); + snprintf(str_filename, filename_len, "strs_%c_%c_%c_%d_%g_%g.dat", lattice_c, 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); diff --git a/src/fracture.h b/src/fracture.h index 17ab8fd..cae1a49 100644 --- a/src/fracture.h +++ b/src/fracture.h @@ -29,6 +29,11 @@ #define CINT_MAX INT_MAX #define CHOL_F(x) cholmod_##x +typedef enum lattice_t { + VORONOI_LATTICE, + SQUARE_LATTICE +} lattice_t; + typedef enum bound_t { FREE_BOUND, CYLINDER_BOUND, @@ -177,3 +182,4 @@ void update_break_data(data_t *data, unsigned int last_broke, double strength, d double get_conductivity(net_t *inst, double *current, cholmod_common *c); +graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c); diff --git a/src/ini_network.c b/src/graph_create.c index cfb8d53..2746310 100644 --- a/src/ini_network.c +++ b/src/graph_create.c @@ -664,3 +664,17 @@ graph_t *ini_voro_graph(unsigned int L, bound_t boundary, bool use_dual, return g; } + + +graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c) { + switch (lattice) { + case (VORONOI_LATTICE): + return ini_voro_graph(L, bound, dual, genfunc_hyperuniform, c); + case (SQUARE_LATTICE): + return ini_square_network(L, bound, true, c); + default: + printf("lattice type unsupported\n"); + exit(EXIT_FAILURE); + } +} + |