summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--makefile4
-rw-r--r--makefile_hal4
-rw-r--r--src/fracture.c819
-rw-r--r--src/fracture.h6
-rw-r--r--src/graph_create.c (renamed from src/ini_network.c)14
5 files changed, 400 insertions, 447 deletions
diff --git a/makefile b/makefile
index 14954e7..7b7b3e4 100644
--- a/makefile
+++ b/makefile
@@ -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);
+ }
+}
+