From 6590154ae3e4ee97e5e1a2792f9f2ebf716ed251 Mon Sep 17 00:00:00 2001
From: pants <jaron@kent-dobias.com>
Date: Wed, 7 Sep 2016 17:30:19 -0400
Subject: 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

---
 makefile           |   4 +-
 makefile_hal       |   4 +-
 src/fracture.c     | 819 ++++++++++++++++++++++++-----------------------------
 src/fracture.h     |   6 +
 src/graph_create.c | 680 ++++++++++++++++++++++++++++++++++++++++++++
 src/ini_network.c  | 666 -------------------------------------------
 6 files changed, 1066 insertions(+), 1113 deletions(-)
 create mode 100644 src/graph_create.c
 delete mode 100644 src/ini_network.c

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/graph_create.c b/src/graph_create.c
new file mode 100644
index 0000000..2746310
--- /dev/null
+++ b/src/graph_create.c
@@ -0,0 +1,680 @@
+
+#include "fracture.h"
+
+unsigned int *get_spanning_edges(unsigned int num_edges, unsigned int *edges_to_verts, double *vert_coords, double cut, unsigned int *n) {
+	unsigned int *spanning_edges = (unsigned int *)malloc(num_edges * sizeof(unsigned int));
+	(*n) = 0;
+	for (unsigned int i = 0; i < num_edges; i++) {
+		unsigned int v1, v2;
+		v1 = edges_to_verts[2 * i];
+		v2 = edges_to_verts[2 * i + 1];
+		double v1y, v2y;
+		v1y = vert_coords[2 * v1 + 1];
+		v2y = vert_coords[2 * v2 + 1];
+		if ((fabs(v1y - v2y) < 0.25) && ((v1y < cut && v2y > cut) || (v1y > cut && v2y < cut))) {
+			spanning_edges[*n] = i;
+			(*n)++;
+		}
+	}
+	return spanning_edges;
+}
+
+double *get_edge_coords(unsigned int num_edges, double *vert_coords,
+												unsigned int *edges_to_verts) {
+	double *output = (double *)malloc(2 * num_edges * sizeof(double));
+
+	#pragma omp parallel for
+	for (unsigned int i = 0; i < num_edges; i++) {
+		unsigned int v1, v2;
+		double v1x, v1y, v2x, v2y, dx, dy;
+		v1 = edges_to_verts[2 * i];
+		v2 = edges_to_verts[2 * i + 1];
+		output[2 * i] = 0;
+		output[2 * i + 1] = 0;
+		v1x = vert_coords[2 * v1];
+		v1y = vert_coords[2 * v1 + 1];
+		v2x = vert_coords[2 * v2];
+		v2y = vert_coords[2 * v2 + 1];
+		dx = v1x - v2x;
+		dy = v1y - v2y;
+		if (fabs(dx) > 0.5) {
+			if (dx > 0) {
+				v2x = v2x + 1;
+			} else {
+				v1x = v1x + 1;
+			}
+		}
+		if (fabs(dy) > 0.5) {
+			if (dy > 0) {
+				v2y = v2y + 1;
+			} else {
+				v1y = v1y + 1;
+			}
+		}
+		output[2 * i] = (v1x + v2x) / 2;
+		output[2 * i + 1] = (v1y + v2y) / 2;
+	}
+
+	return output;
+}
+
+unsigned int *get_verts_to_edges_ind(unsigned int num_verts,
+																		 unsigned int num_edges,
+																		 const unsigned int *edges_to_verts) {
+	unsigned int *output =
+			(unsigned int *)calloc(num_verts + 1, sizeof(unsigned int));
+	assert(output != NULL);
+
+	for (unsigned int i = 0; i < 2 * num_edges; i++) {
+		if (edges_to_verts[i] < num_verts) {
+			output[edges_to_verts[i] + 1]++;
+		}
+	}
+
+	for (unsigned int i = 0; i < num_verts; i++) {
+		output[i + 1] += output[i];
+	}
+
+	return output;
+}
+
+unsigned int *get_verts_to_edges(unsigned int num_verts, unsigned int num_edges,
+																 const unsigned int *edges_to_verts,
+																 const unsigned int *verts_to_edges_ind) {
+	unsigned int *output = (unsigned int *)calloc(verts_to_edges_ind[num_verts],
+																								sizeof(unsigned int));
+	unsigned int *counts =
+			(unsigned int *)calloc(num_verts, sizeof(unsigned int));
+	for (int i = 0; i < 2 * num_edges; i++) {
+		if (edges_to_verts[i] < num_verts) {
+			output[verts_to_edges_ind[edges_to_verts[i]] +
+						 counts[edges_to_verts[i]]] = i / 2;
+			counts[edges_to_verts[i]]++;
+		}
+	}
+
+	free(counts);
+
+	return output;
+}
+
+graph_t *ini_square_network(unsigned int width, bound_t boundary, bool side_bounds,
+												 cholmod_common *c) {
+	graph_t *network = (graph_t *)calloc(1, sizeof(graph_t));
+
+	network->boundary = boundary;
+	bool periodic = (boundary == CYLINDER_BOUND) || (boundary == TORUS_BOUND) ? true : false;
+
+	network->ne = pow(width, 2);
+	if (boundary == CYLINDER_BOUND) {
+		assert(width % 2 == 0);
+		assert(!side_bounds);
+		network->nv = (width / 2) * (width + 1);
+		network->nv_break = (width / 2) * (width + 1);
+		network->dnv = (width / 2 + 1) * (width / 2) + pow(width / 2, 2);
+		network->num_bounds = 2;
+		network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) *
+																								 sizeof(unsigned int));
+		network->bound_inds[0] = 0;
+		network->bound_inds[1] = width / 2;
+		network->bound_inds[2] = width;
+		network->bound_verts = (unsigned int *)calloc(width, sizeof(unsigned int));
+		network->break_dim = network->nv + network->num_bounds;
+	} else if (boundary == FREE_BOUND) {
+		network->nv = 2 * ((width + 1) / 2) * (width / 2 + 1);
+		network->nv_break = 2 * ((width + 1) / 2) * (width / 2 + 1);
+		network->dnv = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2);
+		if (side_bounds)
+			network->num_bounds = 4;
+		else
+			network->num_bounds = 2;
+		network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) *
+																								 sizeof(unsigned int));
+		for (unsigned int i = 0; i < network->num_bounds + 1; i++) {
+			network->bound_inds[i] = i * ((width + 1) / 2);
+		}
+		network->bound_verts = (unsigned int *)calloc(
+				network->num_bounds * ((width + 1) / 2), sizeof(unsigned int));
+		if (side_bounds) {
+			for (unsigned int i = 0; i < (width + 1) / 2; i++) {
+				network->bound_verts[2 * ((width + 1) / 2) + i] =
+						(width + 1) / 2 + i * (width + 1);
+				if (width % 2) {
+					network->bound_verts[3 * ((width + 1) / 2) + i] =
+							(width + 1) / 2 + i * (width + 1) - 1;
+				} else {
+					network->bound_verts[3 * ((width + 1) / 2) + i] =
+							(width + 1) / 2 + i * (width + 1) + width / 2;
+				}
+			}
+		}
+		network->break_dim = network->nv + network->num_bounds;
+	} else if (boundary == TORUS_BOUND) {
+		network->nv = (width / 2) * (width + 1) - (width / 2);
+		network->nv_break = (width / 2) * (width + 1);
+		network->dnv = (width / 2 + 1) * (width / 2) + pow(width / 2, 2) - (width / 2);
+		network->num_bounds = 1;
+		network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) *
+																								 sizeof(unsigned int));
+		network->bound_inds[0] = 0;
+		network->bound_inds[1] = width / 2;
+		network->bound_verts = (unsigned int *)calloc(width / 2, sizeof(unsigned int));
+		network->break_dim = network->nv_break;
+	}
+	if (boundary != TORUS_BOUND) {
+		for (unsigned int i = 0; i < (width + 1) / 2; i++) {
+			network->bound_verts[i] = i;
+			network->bound_verts[(width + 1) / 2 + i] = network->nv - 1 - i;
+		} 
+	} else {
+		for (unsigned int i = 0; i < width / 2; i++) {
+			network->bound_verts[i] = i;
+		}
+	}
+	network->ev_break =
+			(unsigned int *)calloc(2 * network->ne, sizeof(unsigned int));
+	network->ev =
+			(unsigned int *)calloc(2 * network->ne, sizeof(unsigned int));
+	for (unsigned int i = 0; i < network->ne; i++) {
+		network->ev_break[2 * i] = edge_to_verts(width, periodic, i, 1);
+		network->ev_break[2 * i + 1] = edge_to_verts(width, periodic, i, 0);
+		network->ev[2 * i] = network->ev_break[2 * i] % network->nv;
+		network->ev[2 * i + 1] = network->ev_break[2 * i + 1] % network->nv;
+	}
+	network->vei =
+			(unsigned int *)calloc(network->nv + 1, sizeof(unsigned int));
+	network->vei[0] = 0;
+	unsigned int pos1 = 0;
+	for (unsigned int i = 0; i < network->nv; i++) {
+		bool in_bound = false;
+		for (unsigned int j = 0; j < network->num_bounds; j++) {
+			for (unsigned int k = 0;
+					 k < network->bound_inds[j + 1] - network->bound_inds[j]; k++) {
+				if (i == network->bound_verts[network->bound_inds[j] + k]) {
+					in_bound = true;
+					break;
+				}
+			}
+		}
+		if (in_bound)
+			pos1 += 2;
+		else
+			pos1 += 4;
+
+		network->vei[i + 1] = pos1;
+	}
+	network->ve = (unsigned int *)calloc(
+			network->vei[network->nv], sizeof(unsigned int));
+	unsigned int *vert_counts =
+			(unsigned int *)calloc(network->nv, sizeof(unsigned int));
+	for (unsigned int i = 0; i < network->ne; i++) {
+		unsigned int v0 = network->ev[2 * i];
+		unsigned int v1 = network->ev[2 * i + 1];
+		network->ve[network->vei[v0] + vert_counts[v0]] =
+				i;
+		network->ve[network->vei[v1] + vert_counts[v1]] =
+				i;
+		vert_counts[v0]++;
+		vert_counts[v1]++;
+	}
+	free(vert_counts);
+
+	network->vx =
+			(double *)malloc(2 * network->nv * sizeof(double));
+	for (unsigned int i = 0; i < network->nv; i++) {
+		if (!periodic) {
+		network->vx[2 * i] = ((double)((2 * i + 1) / (width + 1)))/width;
+		network->vx[2 * i + 1] = ((double)((2 * i + 1) % (width + 1)))/width;
+		} else {
+		network->vx[2 * i] = ((double)((2 * i + 1) / (width)))/width;
+		network->vx[2 * i + 1] = ((double)((((2 * i + 1) / (width)+1) % 2) + ((2 * i) % (width))))/width;
+		}
+	}
+
+	network->L = width;
+	network->ex = get_edge_coords(
+			network->ne, network->vx, network->ev);
+
+	network->dev =
+			(unsigned int *)malloc(2 * network->ne * sizeof(unsigned int));
+	for (unsigned int i = 0; i < network->ne; i++) {
+		network->dev[2 * i] =
+				dual_edge_to_verts(width, periodic, i, 0) % network->dnv;
+		network->dev[2 * i + 1] =
+				dual_edge_to_verts(width, periodic, i, 1) % network->dnv;
+	}
+	network->dvx =
+			(double *)malloc(2 * network->dnv * sizeof(double));
+	for (unsigned int i = 0; i < network->dnv; i++) {
+		network->dvx[2 * i] =
+				2*dual_vert_to_coord(width, periodic, i, 0);
+		network->dvx[2 * i + 1] =
+				2*dual_vert_to_coord(width, periodic, i, 1);
+	}
+
+	network->voltcurmat = gen_voltcurmat(network->ne,
+																			 network->break_dim,
+																			 network->ev_break, c);
+
+	network->dvei =
+			get_verts_to_edges_ind(network->dnv, network->ne,
+														 network->dev);
+	network->dve = get_verts_to_edges(
+			network->dnv, network->ne, network->dev,
+			network->dvei);
+	network->spanning_edges = get_spanning_edges(network->ne, network->ev_break, network->vx, 0.51, &(network->num_spanning_edges));
+
+	return network;
+}
+
+unsigned int *get_voro_dual_edges(unsigned int num_edges,
+																	unsigned int num_verts, unsigned int *edges,
+																	unsigned int *triangles) {
+	unsigned int *dual_edges =
+			(unsigned int *)malloc(2 * num_edges * sizeof(unsigned int));
+	unsigned int place = 0;
+	#pragma omp parallel for
+	for (unsigned int i = 0; i < num_edges; i++) {
+		unsigned int v1, v2;
+		v1 = edges[2 * i];
+		v2 = edges[2 * i + 1];
+		if (v1 < num_verts && v2 < num_verts) {
+			bool found_match = false;
+			for (unsigned int j = 0; j < 3; j++) {
+				for (unsigned int k = 0; k < 3; k++) {
+					unsigned int t11, t12, t21, t22;
+					t11 = triangles[3 * v1 + j];
+					t12 = triangles[3 * v1 + ((j + 1) % 3)];
+					t21 = triangles[3 * v2 + k];
+					t22 = triangles[3 * v2 + ((k + 1) % 3)];
+					if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) {
+						dual_edges[2 * place] = t11 < t12 ? t11 : t12;
+						dual_edges[2 * place + 1] = t11 < t12 ? t12 : t11;
+						place++;
+						found_match = true;
+						break;
+					}
+				}
+				if (found_match)
+					break;
+			}
+		}
+	}
+
+	return dual_edges;
+}
+
+graph_t *ini_voro_graph(unsigned int L, bound_t boundary, bool use_dual,
+													double *(*genfunc)(unsigned int, bound_t, gsl_rng *, unsigned int *),
+													cholmod_common *c) {
+	graph_t *g = (graph_t *)calloc(1, sizeof(graph_t));
+
+	// generate the dual lattice
+	double *lattice;
+	unsigned int num;
+	{
+		gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+		FILE *rf = fopen("/dev/urandom", "r");
+		unsigned long int seed;
+		fread(&seed, sizeof(unsigned long int), 1, rf);
+		fclose(rf);
+		gsl_rng_set(r, seed);
+		lattice = genfunc(L, boundary, r, &num);
+		gsl_rng_free(r);
+	}
+
+	// retrieve a periodic voronoi tesselation of the lattice
+	bool run_periodic;
+	if (boundary == EMBEDDED_BOUND) run_periodic = false;
+	else run_periodic = true;
+	intptr_t *vout = run_voronoi(num, lattice, run_periodic, 0, 1, 0, 1);
+
+	unsigned int tmp_num_verts = ((unsigned int *)vout[0])[0];
+	unsigned int tmp_num_edges = ((unsigned int *)vout[0])[1];
+	double *tmp_vert_coords = (double *)vout[2];
+	unsigned int *tmp_edges = (unsigned int *)vout[3];
+	unsigned int *tmp_tris = (unsigned int *)vout[5];
+
+	free((void *)vout[0]);
+	free((void *)vout[1]);
+	free((void *)vout[4]);
+	free(vout);
+
+	// get dual edges of the fully periodic graph
+	unsigned int *tmp_dual_edges =
+			get_voro_dual_edges(tmp_num_edges, tmp_num_verts, tmp_edges, tmp_tris);
+
+	// when use_dual is specificed, the edge and vertex sets are swapped with the
+	// dual edge and dual vertex sets.  once formally relabelled, everything
+	// works the same way
+	if (use_dual) {
+		unsigned int *tmp_tmp_dual_edges = tmp_edges;
+		double *tmp_lattice = tmp_vert_coords;
+		unsigned int tmp_num = tmp_num_verts;
+
+		tmp_edges = tmp_dual_edges;
+		tmp_dual_edges = tmp_tmp_dual_edges;
+
+		tmp_vert_coords = lattice;
+		lattice = tmp_lattice;
+
+		tmp_num_verts = num;
+		num = tmp_num;
+	}
+
+	// prune the edges of the lattice and assign boundary vertices based on the
+	// desired boundary conditions
+	unsigned int num_bounds;
+	unsigned int num_verts;
+	double *vert_coords;
+	unsigned int *bound_inds;
+	unsigned int *bound_verts;
+	unsigned int num_edges;
+	unsigned int *edges;
+	unsigned int *dual_edges;
+	switch (boundary) {
+		case FREE_BOUND: {
+			num_bounds = 4;
+			bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int));
+			bound_inds[0] = 0;
+			vert_coords = tmp_vert_coords;
+			num_verts = tmp_num_verts;
+			num_edges = 0;
+			edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
+			dual_edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
+			unsigned int num_t, num_b, num_l, num_r;
+			bool *bound_top, *bound_b, *bound_l, *bound_r;
+			num_t = 0; num_b = 0; num_l = 0; num_r = 0;
+			bound_top = (bool *)calloc(num_verts, sizeof(bool));
+			bound_b = (bool *)calloc(num_verts, sizeof(bool));
+			bound_l = (bool *)calloc(num_verts, sizeof(bool));
+			bound_r = (bool *)calloc(num_verts, sizeof(bool));
+			for (unsigned int i = 0; i < tmp_num_edges; i++) {
+				unsigned int v1, v2;
+				double v1x, v1y, v2x, v2y, dx, dy;
+				v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
+				v1x = vert_coords[2 * v1]; v1y = vert_coords[2 * v1 + 1];
+				v2x = vert_coords[2 * v2]; v2y = vert_coords[2 * v2 + 1];
+				dx = v1x - v2x; dy = v1y - v2y;
+				if (fabs(dy) > 0.5) {
+					if (dy > 0) {
+						if (!bound_top[v1] && !bound_l[v1] && !bound_r[v1]) {
+							bound_top[v1] = true; num_t++;
+						}
+						if (!bound_b[v2] && !bound_l[v2] && !bound_r[v2]) {
+							bound_b[v2] = true; num_b++;
+						}
+					} else {
+						if (!bound_top[v2] && !bound_l[v2] && !bound_r[v2]) {
+							bound_top[v2] = true; num_t++;
+						}
+						if (!bound_b[v1] && !bound_l[v1] && !bound_r[v1]) {
+							bound_b[v1] = true; num_b++;
+						}
+					}
+				} else if (fabs(dx) > 0.5) {
+					if (dx > 0) {
+						if (!bound_r[v1] && !bound_top[v1] && !bound_b[v1]) {
+							bound_r[v1] = true; num_r++;
+						}
+						if (!bound_l[v2] && !bound_top[v2] && !bound_b[v2]) {
+							bound_l[v2] = true; num_l++;
+						}
+					} else {
+						if (!bound_r[v2] && !bound_top[v2] && !bound_b[v2]) {
+							bound_r[v2] = true; num_r++;
+						}
+						if (!bound_l[v1] && !bound_top[v1] && !bound_b[v1]) {
+							bound_l[v1] = true; num_l++;
+						}
+					}
+				} else {
+					edges[2 * num_edges] = v1 < v2 ? v1 : v2;
+					edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
+					unsigned int d1 = tmp_dual_edges[2 * i];
+					unsigned int d2 = tmp_dual_edges[2 * i + 1];
+					dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2;
+					dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1;
+					num_edges++;
+				}
+			}
+			bound_verts = malloc((num_t + num_b + num_l + num_r) * sizeof(unsigned int));
+			bound_inds[1] = num_t; bound_inds[2] = num_t + num_b;
+			bound_inds[3] = num_l + num_t + num_b;
+			bound_inds[4] = num_t + num_b + num_r + num_l;
+			unsigned int pos_t, pos_b, pos_l, pos_r;
+			pos_t = 0; pos_b = 0; pos_l = 0; pos_r = 0;
+			for (unsigned int i = 0; i < num_verts; i++) {
+				if (bound_top[i]) {
+					bound_verts[pos_t] = i; pos_t++;
+				} else if (bound_b[i]) {
+					bound_verts[num_t + pos_b] = i; pos_b++;
+				} else if (bound_l[i]) {
+					bound_verts[num_t + num_b + pos_l] = i; pos_l++;
+				} else if (bound_r[i]) {
+					bound_verts[num_t + num_b + num_l + pos_r] = i; pos_r++;
+				}
+			}
+			free(bound_l); free(bound_r); free(bound_top); free(bound_b);
+			free(tmp_edges);
+			free(tmp_dual_edges);
+			num_bounds = 2;
+			g->ev_break = edges;
+			g->ev = edges;
+			g->nv_break = num_verts;
+			g->nv = num_verts;
+			g->break_dim = num_verts + num_bounds;
+			break;
+										 }
+		case CYLINDER_BOUND: {
+			num_bounds = 2;
+			bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int));
+			bound_inds[0] = 0;
+			vert_coords = tmp_vert_coords;
+			num_verts = tmp_num_verts;
+			num_edges = 0;
+			edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
+			dual_edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
+			unsigned int num_t, num_b;
+			bool *bound_top, *bound_b;
+			num_t = 0; num_b = 0;
+			bound_top = (bool *)calloc(num_verts, sizeof(bool));
+			bound_b = (bool *)calloc(num_verts, sizeof(bool));
+			for (unsigned int i = 0; i < tmp_num_edges; i++) {
+				unsigned int v1, v2;
+				double v1x, v1y, v2x, v2y, dx, dy;
+				v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
+				v1y = vert_coords[2 * v1 + 1]; v2y = vert_coords[2 * v2 + 1];
+				dy = v1y - v2y;
+				if (fabs(dy) > 0.5) {
+					if (dy > 0) {
+						if (!bound_top[v1]) {
+							bound_top[v1] = true; num_t++;
+						}
+						if (!bound_b[v2]) {
+							bound_b[v2] = true; num_b++;
+						}
+					} else {
+						if (!bound_top[v2]) {
+							bound_top[v2] = true; num_t++;
+						}
+						if (!bound_b[v1]) {
+							bound_b[v1] = true; num_b++;
+						}
+					}
+				} else {
+					edges[2 * num_edges] = v1 < v2 ? v1 : v2;
+					edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
+					unsigned int d1 = tmp_dual_edges[2 * i];
+					unsigned int d2 = tmp_dual_edges[2 * i + 1];
+					dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2;
+					dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1;
+					num_edges++;
+				}
+			}
+			bound_verts = malloc((num_t + num_b) * sizeof(unsigned int));
+			bound_inds[1] = num_t; bound_inds[2] = num_t + num_b;
+			unsigned int pos_t, pos_b;
+			pos_t = 0; pos_b = 0;
+			for (unsigned int i = 0; i < num_verts; i++) {
+				if (bound_top[i]) {
+					bound_verts[pos_t] = i; pos_t++;
+				} else if (bound_b[i]) {
+					bound_verts[num_t + pos_b] = i; pos_b++;
+				}
+			}
+			free(bound_top); free(bound_b);
+			free(tmp_edges);
+			free(tmp_dual_edges);
+			g->ev_break = edges;
+			g->ev = edges;
+			g->nv_break = num_verts;
+			g->nv = num_verts;
+			g->break_dim = num_verts + num_bounds;
+			break;
+												 }
+		case TORUS_BOUND: {
+			num_bounds = 1;
+			bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int));
+			bound_inds[0] = 0;
+			num_edges = tmp_num_edges;
+			edges = (unsigned int *)malloc(2* num_edges*sizeof(unsigned int));
+			for (unsigned int i = 0; i < num_edges; i++) {
+				edges[2*i] = tmp_edges[2*i];
+				edges[2*i+1] = tmp_edges[2*i+1];
+			}
+			dual_edges = tmp_dual_edges;
+			bool *bound_top = (bool *)calloc(tmp_num_verts, sizeof(bool));
+			int *edge_change = (int *)calloc(num_edges, sizeof(int));
+			unsigned int num_t = 0;
+			for (unsigned int i = 0; i < num_edges; i++) {
+				unsigned int v1, v2;
+				double v1x, v1y, v2x, v2y, dx, dy;
+				v1 = edges[2 * i]; v2 = edges[2 * i + 1];
+				v1x = tmp_vert_coords[2 * v1]; v1y = tmp_vert_coords[2 * v1 + 1];
+				v2x = tmp_vert_coords[2 * v2]; v2y = tmp_vert_coords[2 * v2 + 1];
+				dy = v1y - v2y;
+				if (fabs(dy) > 0.5) {
+					if (dy > 0) {
+						edge_change[i] = 1;
+						if (!bound_top[v1]) {
+							bound_top[v1] = true; num_t++;
+						}
+					} else {
+						edge_change[i] = 2;
+						if (!bound_top[v2]) {
+							bound_top[v2] = true; num_t++;
+						}
+					}
+				}
+			}
+			num_verts = tmp_num_verts + num_t;
+			vert_coords = (double *)malloc(2 * num_verts * sizeof(double));
+			bound_verts = malloc(num_t * sizeof(unsigned int));
+			bound_inds[1] = num_t;
+			unsigned int pos_t = 0;
+			for (unsigned int i = 0; i < tmp_num_verts; i++) {
+				vert_coords[2*i] = tmp_vert_coords[2*i];
+				vert_coords[2*i+1] = tmp_vert_coords[2*i+1];
+				if (bound_top[i]) {
+					bound_verts[pos_t] = i;
+					vert_coords[2*(tmp_num_verts + pos_t)] = tmp_vert_coords[2*i];
+					vert_coords[2*(tmp_num_verts + pos_t)+1] = tmp_vert_coords[2*i+1];
+					pos_t++;
+				}
+			}
+			for (unsigned int i = 0; i < num_edges; i++) {
+				if (edge_change[i]) {
+					for (unsigned int j = 0; j < num_t; j++) {
+						if (edges[2*i+(edge_change[i]-1)] == bound_verts[j]) {
+							edges[2*i+(edge_change[i]-1)] = tmp_num_verts + j;
+							break;
+						}
+					}
+				}
+			}
+			free(tmp_vert_coords);
+			free(bound_top);
+			free(edge_change);
+			g->nv_break = num_verts;
+			g->nv = tmp_num_verts;
+			g->ev_break = edges;
+			g->ev = tmp_edges;
+			g->break_dim = num_verts;
+			break;
+											}
+		case EMBEDDED_BOUND: {
+			num_bounds = 4;
+			bound_inds = (unsigned int *)malloc(5 * sizeof(unsigned int));
+			bound_verts = (unsigned int *)malloc(2 * L * sizeof(unsigned int));
+			for (unsigned int i = 0; i < 5; i++) bound_inds[i] = i * L / 2;
+			for (unsigned int i = 0; i < 2 * L; i++) bound_verts[i] = i;
+			unsigned int num_away = 0;
+			for (unsigned int i = 0; i < tmp_num_edges; i++) {
+				if (tmp_dual_edges[2*i] > num || tmp_dual_edges[2*i+1] > num) num_away++;
+			}
+			num_edges = (int)tmp_num_edges - (int)num_away;
+			num_verts = tmp_num_verts;
+			edges = tmp_edges;
+			dual_edges = tmp_dual_edges;
+			vert_coords = tmp_vert_coords;
+			g->nv_break = num_verts;
+			g->nv = num_verts;
+			g->ev_break = edges;
+			g->ev = edges;
+			g->break_dim = num_verts + 2;
+		}
+	}
+
+	g->boundary = boundary;
+	g->num_bounds = num_bounds;
+	g->bound_inds = bound_inds;
+	g->bound_verts = bound_verts;
+	g->ne = num_edges;
+	g->dev = dual_edges;
+	g->vx = vert_coords;
+
+	g->vei = get_verts_to_edges_ind(g->nv, g->ne, g->ev);
+	g->ve =
+			get_verts_to_edges(g->nv, g->ne,
+												 g->ev, g->vei);
+
+	g->L = L;
+
+	g->ex = get_edge_coords(
+			g->ne, g->vx, g->ev);
+
+	free(tmp_tris);
+
+	g->dvx = lattice;
+	g->dnv = num;
+
+	g->voltcurmat = gen_voltcurmat(g->ne,
+																			 g->break_dim,
+																			 g->ev_break, c);
+
+	g->dvei =
+			get_verts_to_edges_ind(g->dnv, g->ne,
+														 g->dev);
+	g->dve = get_verts_to_edges(
+			g->dnv, g->ne, g->dev,
+			g->dvei);
+
+	g->spanning_edges = get_spanning_edges(g->ne, g->ev_break, g->vx, 0.5, &(g->num_spanning_edges));
+
+	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);
+	}
+}
+
diff --git a/src/ini_network.c b/src/ini_network.c
deleted file mode 100644
index cfb8d53..0000000
--- a/src/ini_network.c
+++ /dev/null
@@ -1,666 +0,0 @@
-
-#include "fracture.h"
-
-unsigned int *get_spanning_edges(unsigned int num_edges, unsigned int *edges_to_verts, double *vert_coords, double cut, unsigned int *n) {
-	unsigned int *spanning_edges = (unsigned int *)malloc(num_edges * sizeof(unsigned int));
-	(*n) = 0;
-	for (unsigned int i = 0; i < num_edges; i++) {
-		unsigned int v1, v2;
-		v1 = edges_to_verts[2 * i];
-		v2 = edges_to_verts[2 * i + 1];
-		double v1y, v2y;
-		v1y = vert_coords[2 * v1 + 1];
-		v2y = vert_coords[2 * v2 + 1];
-		if ((fabs(v1y - v2y) < 0.25) && ((v1y < cut && v2y > cut) || (v1y > cut && v2y < cut))) {
-			spanning_edges[*n] = i;
-			(*n)++;
-		}
-	}
-	return spanning_edges;
-}
-
-double *get_edge_coords(unsigned int num_edges, double *vert_coords,
-												unsigned int *edges_to_verts) {
-	double *output = (double *)malloc(2 * num_edges * sizeof(double));
-
-	#pragma omp parallel for
-	for (unsigned int i = 0; i < num_edges; i++) {
-		unsigned int v1, v2;
-		double v1x, v1y, v2x, v2y, dx, dy;
-		v1 = edges_to_verts[2 * i];
-		v2 = edges_to_verts[2 * i + 1];
-		output[2 * i] = 0;
-		output[2 * i + 1] = 0;
-		v1x = vert_coords[2 * v1];
-		v1y = vert_coords[2 * v1 + 1];
-		v2x = vert_coords[2 * v2];
-		v2y = vert_coords[2 * v2 + 1];
-		dx = v1x - v2x;
-		dy = v1y - v2y;
-		if (fabs(dx) > 0.5) {
-			if (dx > 0) {
-				v2x = v2x + 1;
-			} else {
-				v1x = v1x + 1;
-			}
-		}
-		if (fabs(dy) > 0.5) {
-			if (dy > 0) {
-				v2y = v2y + 1;
-			} else {
-				v1y = v1y + 1;
-			}
-		}
-		output[2 * i] = (v1x + v2x) / 2;
-		output[2 * i + 1] = (v1y + v2y) / 2;
-	}
-
-	return output;
-}
-
-unsigned int *get_verts_to_edges_ind(unsigned int num_verts,
-																		 unsigned int num_edges,
-																		 const unsigned int *edges_to_verts) {
-	unsigned int *output =
-			(unsigned int *)calloc(num_verts + 1, sizeof(unsigned int));
-	assert(output != NULL);
-
-	for (unsigned int i = 0; i < 2 * num_edges; i++) {
-		if (edges_to_verts[i] < num_verts) {
-			output[edges_to_verts[i] + 1]++;
-		}
-	}
-
-	for (unsigned int i = 0; i < num_verts; i++) {
-		output[i + 1] += output[i];
-	}
-
-	return output;
-}
-
-unsigned int *get_verts_to_edges(unsigned int num_verts, unsigned int num_edges,
-																 const unsigned int *edges_to_verts,
-																 const unsigned int *verts_to_edges_ind) {
-	unsigned int *output = (unsigned int *)calloc(verts_to_edges_ind[num_verts],
-																								sizeof(unsigned int));
-	unsigned int *counts =
-			(unsigned int *)calloc(num_verts, sizeof(unsigned int));
-	for (int i = 0; i < 2 * num_edges; i++) {
-		if (edges_to_verts[i] < num_verts) {
-			output[verts_to_edges_ind[edges_to_verts[i]] +
-						 counts[edges_to_verts[i]]] = i / 2;
-			counts[edges_to_verts[i]]++;
-		}
-	}
-
-	free(counts);
-
-	return output;
-}
-
-graph_t *ini_square_network(unsigned int width, bound_t boundary, bool side_bounds,
-												 cholmod_common *c) {
-	graph_t *network = (graph_t *)calloc(1, sizeof(graph_t));
-
-	network->boundary = boundary;
-	bool periodic = (boundary == CYLINDER_BOUND) || (boundary == TORUS_BOUND) ? true : false;
-
-	network->ne = pow(width, 2);
-	if (boundary == CYLINDER_BOUND) {
-		assert(width % 2 == 0);
-		assert(!side_bounds);
-		network->nv = (width / 2) * (width + 1);
-		network->nv_break = (width / 2) * (width + 1);
-		network->dnv = (width / 2 + 1) * (width / 2) + pow(width / 2, 2);
-		network->num_bounds = 2;
-		network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) *
-																								 sizeof(unsigned int));
-		network->bound_inds[0] = 0;
-		network->bound_inds[1] = width / 2;
-		network->bound_inds[2] = width;
-		network->bound_verts = (unsigned int *)calloc(width, sizeof(unsigned int));
-		network->break_dim = network->nv + network->num_bounds;
-	} else if (boundary == FREE_BOUND) {
-		network->nv = 2 * ((width + 1) / 2) * (width / 2 + 1);
-		network->nv_break = 2 * ((width + 1) / 2) * (width / 2 + 1);
-		network->dnv = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2);
-		if (side_bounds)
-			network->num_bounds = 4;
-		else
-			network->num_bounds = 2;
-		network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) *
-																								 sizeof(unsigned int));
-		for (unsigned int i = 0; i < network->num_bounds + 1; i++) {
-			network->bound_inds[i] = i * ((width + 1) / 2);
-		}
-		network->bound_verts = (unsigned int *)calloc(
-				network->num_bounds * ((width + 1) / 2), sizeof(unsigned int));
-		if (side_bounds) {
-			for (unsigned int i = 0; i < (width + 1) / 2; i++) {
-				network->bound_verts[2 * ((width + 1) / 2) + i] =
-						(width + 1) / 2 + i * (width + 1);
-				if (width % 2) {
-					network->bound_verts[3 * ((width + 1) / 2) + i] =
-							(width + 1) / 2 + i * (width + 1) - 1;
-				} else {
-					network->bound_verts[3 * ((width + 1) / 2) + i] =
-							(width + 1) / 2 + i * (width + 1) + width / 2;
-				}
-			}
-		}
-		network->break_dim = network->nv + network->num_bounds;
-	} else if (boundary == TORUS_BOUND) {
-		network->nv = (width / 2) * (width + 1) - (width / 2);
-		network->nv_break = (width / 2) * (width + 1);
-		network->dnv = (width / 2 + 1) * (width / 2) + pow(width / 2, 2) - (width / 2);
-		network->num_bounds = 1;
-		network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) *
-																								 sizeof(unsigned int));
-		network->bound_inds[0] = 0;
-		network->bound_inds[1] = width / 2;
-		network->bound_verts = (unsigned int *)calloc(width / 2, sizeof(unsigned int));
-		network->break_dim = network->nv_break;
-	}
-	if (boundary != TORUS_BOUND) {
-		for (unsigned int i = 0; i < (width + 1) / 2; i++) {
-			network->bound_verts[i] = i;
-			network->bound_verts[(width + 1) / 2 + i] = network->nv - 1 - i;
-		} 
-	} else {
-		for (unsigned int i = 0; i < width / 2; i++) {
-			network->bound_verts[i] = i;
-		}
-	}
-	network->ev_break =
-			(unsigned int *)calloc(2 * network->ne, sizeof(unsigned int));
-	network->ev =
-			(unsigned int *)calloc(2 * network->ne, sizeof(unsigned int));
-	for (unsigned int i = 0; i < network->ne; i++) {
-		network->ev_break[2 * i] = edge_to_verts(width, periodic, i, 1);
-		network->ev_break[2 * i + 1] = edge_to_verts(width, periodic, i, 0);
-		network->ev[2 * i] = network->ev_break[2 * i] % network->nv;
-		network->ev[2 * i + 1] = network->ev_break[2 * i + 1] % network->nv;
-	}
-	network->vei =
-			(unsigned int *)calloc(network->nv + 1, sizeof(unsigned int));
-	network->vei[0] = 0;
-	unsigned int pos1 = 0;
-	for (unsigned int i = 0; i < network->nv; i++) {
-		bool in_bound = false;
-		for (unsigned int j = 0; j < network->num_bounds; j++) {
-			for (unsigned int k = 0;
-					 k < network->bound_inds[j + 1] - network->bound_inds[j]; k++) {
-				if (i == network->bound_verts[network->bound_inds[j] + k]) {
-					in_bound = true;
-					break;
-				}
-			}
-		}
-		if (in_bound)
-			pos1 += 2;
-		else
-			pos1 += 4;
-
-		network->vei[i + 1] = pos1;
-	}
-	network->ve = (unsigned int *)calloc(
-			network->vei[network->nv], sizeof(unsigned int));
-	unsigned int *vert_counts =
-			(unsigned int *)calloc(network->nv, sizeof(unsigned int));
-	for (unsigned int i = 0; i < network->ne; i++) {
-		unsigned int v0 = network->ev[2 * i];
-		unsigned int v1 = network->ev[2 * i + 1];
-		network->ve[network->vei[v0] + vert_counts[v0]] =
-				i;
-		network->ve[network->vei[v1] + vert_counts[v1]] =
-				i;
-		vert_counts[v0]++;
-		vert_counts[v1]++;
-	}
-	free(vert_counts);
-
-	network->vx =
-			(double *)malloc(2 * network->nv * sizeof(double));
-	for (unsigned int i = 0; i < network->nv; i++) {
-		if (!periodic) {
-		network->vx[2 * i] = ((double)((2 * i + 1) / (width + 1)))/width;
-		network->vx[2 * i + 1] = ((double)((2 * i + 1) % (width + 1)))/width;
-		} else {
-		network->vx[2 * i] = ((double)((2 * i + 1) / (width)))/width;
-		network->vx[2 * i + 1] = ((double)((((2 * i + 1) / (width)+1) % 2) + ((2 * i) % (width))))/width;
-		}
-	}
-
-	network->L = width;
-	network->ex = get_edge_coords(
-			network->ne, network->vx, network->ev);
-
-	network->dev =
-			(unsigned int *)malloc(2 * network->ne * sizeof(unsigned int));
-	for (unsigned int i = 0; i < network->ne; i++) {
-		network->dev[2 * i] =
-				dual_edge_to_verts(width, periodic, i, 0) % network->dnv;
-		network->dev[2 * i + 1] =
-				dual_edge_to_verts(width, periodic, i, 1) % network->dnv;
-	}
-	network->dvx =
-			(double *)malloc(2 * network->dnv * sizeof(double));
-	for (unsigned int i = 0; i < network->dnv; i++) {
-		network->dvx[2 * i] =
-				2*dual_vert_to_coord(width, periodic, i, 0);
-		network->dvx[2 * i + 1] =
-				2*dual_vert_to_coord(width, periodic, i, 1);
-	}
-
-	network->voltcurmat = gen_voltcurmat(network->ne,
-																			 network->break_dim,
-																			 network->ev_break, c);
-
-	network->dvei =
-			get_verts_to_edges_ind(network->dnv, network->ne,
-														 network->dev);
-	network->dve = get_verts_to_edges(
-			network->dnv, network->ne, network->dev,
-			network->dvei);
-	network->spanning_edges = get_spanning_edges(network->ne, network->ev_break, network->vx, 0.51, &(network->num_spanning_edges));
-
-	return network;
-}
-
-unsigned int *get_voro_dual_edges(unsigned int num_edges,
-																	unsigned int num_verts, unsigned int *edges,
-																	unsigned int *triangles) {
-	unsigned int *dual_edges =
-			(unsigned int *)malloc(2 * num_edges * sizeof(unsigned int));
-	unsigned int place = 0;
-	#pragma omp parallel for
-	for (unsigned int i = 0; i < num_edges; i++) {
-		unsigned int v1, v2;
-		v1 = edges[2 * i];
-		v2 = edges[2 * i + 1];
-		if (v1 < num_verts && v2 < num_verts) {
-			bool found_match = false;
-			for (unsigned int j = 0; j < 3; j++) {
-				for (unsigned int k = 0; k < 3; k++) {
-					unsigned int t11, t12, t21, t22;
-					t11 = triangles[3 * v1 + j];
-					t12 = triangles[3 * v1 + ((j + 1) % 3)];
-					t21 = triangles[3 * v2 + k];
-					t22 = triangles[3 * v2 + ((k + 1) % 3)];
-					if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) {
-						dual_edges[2 * place] = t11 < t12 ? t11 : t12;
-						dual_edges[2 * place + 1] = t11 < t12 ? t12 : t11;
-						place++;
-						found_match = true;
-						break;
-					}
-				}
-				if (found_match)
-					break;
-			}
-		}
-	}
-
-	return dual_edges;
-}
-
-graph_t *ini_voro_graph(unsigned int L, bound_t boundary, bool use_dual,
-													double *(*genfunc)(unsigned int, bound_t, gsl_rng *, unsigned int *),
-													cholmod_common *c) {
-	graph_t *g = (graph_t *)calloc(1, sizeof(graph_t));
-
-	// generate the dual lattice
-	double *lattice;
-	unsigned int num;
-	{
-		gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
-		FILE *rf = fopen("/dev/urandom", "r");
-		unsigned long int seed;
-		fread(&seed, sizeof(unsigned long int), 1, rf);
-		fclose(rf);
-		gsl_rng_set(r, seed);
-		lattice = genfunc(L, boundary, r, &num);
-		gsl_rng_free(r);
-	}
-
-	// retrieve a periodic voronoi tesselation of the lattice
-	bool run_periodic;
-	if (boundary == EMBEDDED_BOUND) run_periodic = false;
-	else run_periodic = true;
-	intptr_t *vout = run_voronoi(num, lattice, run_periodic, 0, 1, 0, 1);
-
-	unsigned int tmp_num_verts = ((unsigned int *)vout[0])[0];
-	unsigned int tmp_num_edges = ((unsigned int *)vout[0])[1];
-	double *tmp_vert_coords = (double *)vout[2];
-	unsigned int *tmp_edges = (unsigned int *)vout[3];
-	unsigned int *tmp_tris = (unsigned int *)vout[5];
-
-	free((void *)vout[0]);
-	free((void *)vout[1]);
-	free((void *)vout[4]);
-	free(vout);
-
-	// get dual edges of the fully periodic graph
-	unsigned int *tmp_dual_edges =
-			get_voro_dual_edges(tmp_num_edges, tmp_num_verts, tmp_edges, tmp_tris);
-
-	// when use_dual is specificed, the edge and vertex sets are swapped with the
-	// dual edge and dual vertex sets.  once formally relabelled, everything
-	// works the same way
-	if (use_dual) {
-		unsigned int *tmp_tmp_dual_edges = tmp_edges;
-		double *tmp_lattice = tmp_vert_coords;
-		unsigned int tmp_num = tmp_num_verts;
-
-		tmp_edges = tmp_dual_edges;
-		tmp_dual_edges = tmp_tmp_dual_edges;
-
-		tmp_vert_coords = lattice;
-		lattice = tmp_lattice;
-
-		tmp_num_verts = num;
-		num = tmp_num;
-	}
-
-	// prune the edges of the lattice and assign boundary vertices based on the
-	// desired boundary conditions
-	unsigned int num_bounds;
-	unsigned int num_verts;
-	double *vert_coords;
-	unsigned int *bound_inds;
-	unsigned int *bound_verts;
-	unsigned int num_edges;
-	unsigned int *edges;
-	unsigned int *dual_edges;
-	switch (boundary) {
-		case FREE_BOUND: {
-			num_bounds = 4;
-			bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int));
-			bound_inds[0] = 0;
-			vert_coords = tmp_vert_coords;
-			num_verts = tmp_num_verts;
-			num_edges = 0;
-			edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
-			dual_edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
-			unsigned int num_t, num_b, num_l, num_r;
-			bool *bound_top, *bound_b, *bound_l, *bound_r;
-			num_t = 0; num_b = 0; num_l = 0; num_r = 0;
-			bound_top = (bool *)calloc(num_verts, sizeof(bool));
-			bound_b = (bool *)calloc(num_verts, sizeof(bool));
-			bound_l = (bool *)calloc(num_verts, sizeof(bool));
-			bound_r = (bool *)calloc(num_verts, sizeof(bool));
-			for (unsigned int i = 0; i < tmp_num_edges; i++) {
-				unsigned int v1, v2;
-				double v1x, v1y, v2x, v2y, dx, dy;
-				v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
-				v1x = vert_coords[2 * v1]; v1y = vert_coords[2 * v1 + 1];
-				v2x = vert_coords[2 * v2]; v2y = vert_coords[2 * v2 + 1];
-				dx = v1x - v2x; dy = v1y - v2y;
-				if (fabs(dy) > 0.5) {
-					if (dy > 0) {
-						if (!bound_top[v1] && !bound_l[v1] && !bound_r[v1]) {
-							bound_top[v1] = true; num_t++;
-						}
-						if (!bound_b[v2] && !bound_l[v2] && !bound_r[v2]) {
-							bound_b[v2] = true; num_b++;
-						}
-					} else {
-						if (!bound_top[v2] && !bound_l[v2] && !bound_r[v2]) {
-							bound_top[v2] = true; num_t++;
-						}
-						if (!bound_b[v1] && !bound_l[v1] && !bound_r[v1]) {
-							bound_b[v1] = true; num_b++;
-						}
-					}
-				} else if (fabs(dx) > 0.5) {
-					if (dx > 0) {
-						if (!bound_r[v1] && !bound_top[v1] && !bound_b[v1]) {
-							bound_r[v1] = true; num_r++;
-						}
-						if (!bound_l[v2] && !bound_top[v2] && !bound_b[v2]) {
-							bound_l[v2] = true; num_l++;
-						}
-					} else {
-						if (!bound_r[v2] && !bound_top[v2] && !bound_b[v2]) {
-							bound_r[v2] = true; num_r++;
-						}
-						if (!bound_l[v1] && !bound_top[v1] && !bound_b[v1]) {
-							bound_l[v1] = true; num_l++;
-						}
-					}
-				} else {
-					edges[2 * num_edges] = v1 < v2 ? v1 : v2;
-					edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
-					unsigned int d1 = tmp_dual_edges[2 * i];
-					unsigned int d2 = tmp_dual_edges[2 * i + 1];
-					dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2;
-					dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1;
-					num_edges++;
-				}
-			}
-			bound_verts = malloc((num_t + num_b + num_l + num_r) * sizeof(unsigned int));
-			bound_inds[1] = num_t; bound_inds[2] = num_t + num_b;
-			bound_inds[3] = num_l + num_t + num_b;
-			bound_inds[4] = num_t + num_b + num_r + num_l;
-			unsigned int pos_t, pos_b, pos_l, pos_r;
-			pos_t = 0; pos_b = 0; pos_l = 0; pos_r = 0;
-			for (unsigned int i = 0; i < num_verts; i++) {
-				if (bound_top[i]) {
-					bound_verts[pos_t] = i; pos_t++;
-				} else if (bound_b[i]) {
-					bound_verts[num_t + pos_b] = i; pos_b++;
-				} else if (bound_l[i]) {
-					bound_verts[num_t + num_b + pos_l] = i; pos_l++;
-				} else if (bound_r[i]) {
-					bound_verts[num_t + num_b + num_l + pos_r] = i; pos_r++;
-				}
-			}
-			free(bound_l); free(bound_r); free(bound_top); free(bound_b);
-			free(tmp_edges);
-			free(tmp_dual_edges);
-			num_bounds = 2;
-			g->ev_break = edges;
-			g->ev = edges;
-			g->nv_break = num_verts;
-			g->nv = num_verts;
-			g->break_dim = num_verts + num_bounds;
-			break;
-										 }
-		case CYLINDER_BOUND: {
-			num_bounds = 2;
-			bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int));
-			bound_inds[0] = 0;
-			vert_coords = tmp_vert_coords;
-			num_verts = tmp_num_verts;
-			num_edges = 0;
-			edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
-			dual_edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
-			unsigned int num_t, num_b;
-			bool *bound_top, *bound_b;
-			num_t = 0; num_b = 0;
-			bound_top = (bool *)calloc(num_verts, sizeof(bool));
-			bound_b = (bool *)calloc(num_verts, sizeof(bool));
-			for (unsigned int i = 0; i < tmp_num_edges; i++) {
-				unsigned int v1, v2;
-				double v1x, v1y, v2x, v2y, dx, dy;
-				v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
-				v1y = vert_coords[2 * v1 + 1]; v2y = vert_coords[2 * v2 + 1];
-				dy = v1y - v2y;
-				if (fabs(dy) > 0.5) {
-					if (dy > 0) {
-						if (!bound_top[v1]) {
-							bound_top[v1] = true; num_t++;
-						}
-						if (!bound_b[v2]) {
-							bound_b[v2] = true; num_b++;
-						}
-					} else {
-						if (!bound_top[v2]) {
-							bound_top[v2] = true; num_t++;
-						}
-						if (!bound_b[v1]) {
-							bound_b[v1] = true; num_b++;
-						}
-					}
-				} else {
-					edges[2 * num_edges] = v1 < v2 ? v1 : v2;
-					edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
-					unsigned int d1 = tmp_dual_edges[2 * i];
-					unsigned int d2 = tmp_dual_edges[2 * i + 1];
-					dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2;
-					dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1;
-					num_edges++;
-				}
-			}
-			bound_verts = malloc((num_t + num_b) * sizeof(unsigned int));
-			bound_inds[1] = num_t; bound_inds[2] = num_t + num_b;
-			unsigned int pos_t, pos_b;
-			pos_t = 0; pos_b = 0;
-			for (unsigned int i = 0; i < num_verts; i++) {
-				if (bound_top[i]) {
-					bound_verts[pos_t] = i; pos_t++;
-				} else if (bound_b[i]) {
-					bound_verts[num_t + pos_b] = i; pos_b++;
-				}
-			}
-			free(bound_top); free(bound_b);
-			free(tmp_edges);
-			free(tmp_dual_edges);
-			g->ev_break = edges;
-			g->ev = edges;
-			g->nv_break = num_verts;
-			g->nv = num_verts;
-			g->break_dim = num_verts + num_bounds;
-			break;
-												 }
-		case TORUS_BOUND: {
-			num_bounds = 1;
-			bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int));
-			bound_inds[0] = 0;
-			num_edges = tmp_num_edges;
-			edges = (unsigned int *)malloc(2* num_edges*sizeof(unsigned int));
-			for (unsigned int i = 0; i < num_edges; i++) {
-				edges[2*i] = tmp_edges[2*i];
-				edges[2*i+1] = tmp_edges[2*i+1];
-			}
-			dual_edges = tmp_dual_edges;
-			bool *bound_top = (bool *)calloc(tmp_num_verts, sizeof(bool));
-			int *edge_change = (int *)calloc(num_edges, sizeof(int));
-			unsigned int num_t = 0;
-			for (unsigned int i = 0; i < num_edges; i++) {
-				unsigned int v1, v2;
-				double v1x, v1y, v2x, v2y, dx, dy;
-				v1 = edges[2 * i]; v2 = edges[2 * i + 1];
-				v1x = tmp_vert_coords[2 * v1]; v1y = tmp_vert_coords[2 * v1 + 1];
-				v2x = tmp_vert_coords[2 * v2]; v2y = tmp_vert_coords[2 * v2 + 1];
-				dy = v1y - v2y;
-				if (fabs(dy) > 0.5) {
-					if (dy > 0) {
-						edge_change[i] = 1;
-						if (!bound_top[v1]) {
-							bound_top[v1] = true; num_t++;
-						}
-					} else {
-						edge_change[i] = 2;
-						if (!bound_top[v2]) {
-							bound_top[v2] = true; num_t++;
-						}
-					}
-				}
-			}
-			num_verts = tmp_num_verts + num_t;
-			vert_coords = (double *)malloc(2 * num_verts * sizeof(double));
-			bound_verts = malloc(num_t * sizeof(unsigned int));
-			bound_inds[1] = num_t;
-			unsigned int pos_t = 0;
-			for (unsigned int i = 0; i < tmp_num_verts; i++) {
-				vert_coords[2*i] = tmp_vert_coords[2*i];
-				vert_coords[2*i+1] = tmp_vert_coords[2*i+1];
-				if (bound_top[i]) {
-					bound_verts[pos_t] = i;
-					vert_coords[2*(tmp_num_verts + pos_t)] = tmp_vert_coords[2*i];
-					vert_coords[2*(tmp_num_verts + pos_t)+1] = tmp_vert_coords[2*i+1];
-					pos_t++;
-				}
-			}
-			for (unsigned int i = 0; i < num_edges; i++) {
-				if (edge_change[i]) {
-					for (unsigned int j = 0; j < num_t; j++) {
-						if (edges[2*i+(edge_change[i]-1)] == bound_verts[j]) {
-							edges[2*i+(edge_change[i]-1)] = tmp_num_verts + j;
-							break;
-						}
-					}
-				}
-			}
-			free(tmp_vert_coords);
-			free(bound_top);
-			free(edge_change);
-			g->nv_break = num_verts;
-			g->nv = tmp_num_verts;
-			g->ev_break = edges;
-			g->ev = tmp_edges;
-			g->break_dim = num_verts;
-			break;
-											}
-		case EMBEDDED_BOUND: {
-			num_bounds = 4;
-			bound_inds = (unsigned int *)malloc(5 * sizeof(unsigned int));
-			bound_verts = (unsigned int *)malloc(2 * L * sizeof(unsigned int));
-			for (unsigned int i = 0; i < 5; i++) bound_inds[i] = i * L / 2;
-			for (unsigned int i = 0; i < 2 * L; i++) bound_verts[i] = i;
-			unsigned int num_away = 0;
-			for (unsigned int i = 0; i < tmp_num_edges; i++) {
-				if (tmp_dual_edges[2*i] > num || tmp_dual_edges[2*i+1] > num) num_away++;
-			}
-			num_edges = (int)tmp_num_edges - (int)num_away;
-			num_verts = tmp_num_verts;
-			edges = tmp_edges;
-			dual_edges = tmp_dual_edges;
-			vert_coords = tmp_vert_coords;
-			g->nv_break = num_verts;
-			g->nv = num_verts;
-			g->ev_break = edges;
-			g->ev = edges;
-			g->break_dim = num_verts + 2;
-		}
-	}
-
-	g->boundary = boundary;
-	g->num_bounds = num_bounds;
-	g->bound_inds = bound_inds;
-	g->bound_verts = bound_verts;
-	g->ne = num_edges;
-	g->dev = dual_edges;
-	g->vx = vert_coords;
-
-	g->vei = get_verts_to_edges_ind(g->nv, g->ne, g->ev);
-	g->ve =
-			get_verts_to_edges(g->nv, g->ne,
-												 g->ev, g->vei);
-
-	g->L = L;
-
-	g->ex = get_edge_coords(
-			g->ne, g->vx, g->ev);
-
-	free(tmp_tris);
-
-	g->dvx = lattice;
-	g->dnv = num;
-
-	g->voltcurmat = gen_voltcurmat(g->ne,
-																			 g->break_dim,
-																			 g->ev_break, c);
-
-	g->dvei =
-			get_verts_to_edges_ind(g->dnv, g->ne,
-														 g->dev);
-	g->dve = get_verts_to_edges(
-			g->dnv, g->ne, g->dev,
-			g->dvei);
-
-	g->spanning_edges = get_spanning_edges(g->ne, g->ev_break, g->vx, 0.5, &(g->num_spanning_edges));
-
-	return g;
-}
-- 
cgit v1.2.3-70-g09d2