created new fracture program which has full capability (support for variable lattices, boundaries, notch or no). to do: get embedded square lattice working, add flag for constant lattices

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);
+	}
+}
+