started repo again without all the data files gunking the works

46 files changed, 5826 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..a8276ab
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+bin/*
+obj/*
+data/*
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..2c1ec50
--- /dev/null
+++ b/README.md
@@ -0,0 +1,4 @@
+fracture
+========
+
+tools for simulating fracture in fuse networks
\ No newline at end of file
diff --git a/makefile b/makefile
new file mode 100644
index 0000000..eb9fc4d
--- /dev/null
+++ b/makefile
@@ -0,0 +1,31 @@
+
+CC = clang
+CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -march=native #-fopenmp 
+LDFLAGS = -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lmetis -lgsl -lprofiler -ltcmalloc
+
+OBJ = break_data voronoi_bound_ini bin_values correlations beta_scales randfuncs instance get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry fracture_network get_current cracking_ini 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
+BIN = fracture current_scaling course_grain_square corr_test homo_voronoi_fracture homo_square_fracture compare_voronoi_fracture
+
+all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+
+opt:
+	@echo build options:
+	@echo "CC       = ${CC}"
+	@echo "CFLAGS   = ${CFLAGS}"
+	@echo "LDFLAGS  = ${LDFLAGS}"
+
+obj/%.o: src/%.c
+	@echo CC -c -o $@
+	@${CC} -c -o $@ $< ${CFLAGS}
+
+bin/%: obj/%.o ${OBJ:%=obj/%.o}
+	@echo CC -o $@
+	@${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS}
+
+clean:
+	@echo cleaning:
+	@echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+	@rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+
+.PHONY: all clean
+
diff --git a/makefile_ept b/makefile_ept
new file mode 100644
index 0000000..826fde8
--- /dev/null
+++ b/makefile_ept
@@ -0,0 +1,31 @@
+
+CC = clang
+CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -I/usr/local/cuda/include -I/home/jaron/.local/include -I/usr/lib/gcc/x86_64-linux-gnu/5/include -march=native -DGPU_BLAS -fopenmp=libiomp5 
+LDFLAGS = -L/usr/local/cuda/lib64 -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lmetis -lgsl -lprofiler -ltcmalloc
+
+OBJ = break_data voronoi_bound_ini bin_values correlations beta_scales randfuncs instance get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry fracture_network get_current cracking_ini 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
+BIN = fracture current_scaling course_grain_square corr_test homo_voronoi_fracture
+
+all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+
+opt:
+	@echo build options:
+	@echo "CC       = ${CC}"
+	@echo "CFLAGS   = ${CFLAGS}"
+	@echo "LDFLAGS  = ${LDFLAGS}"
+
+obj/%.o: src/%.c
+	@echo CC -c -o $@
+	@${CC} -c -o $@ $< ${CFLAGS}
+
+bin/%: obj/%.o ${OBJ:%=obj/%.o}
+	@echo CC -o $@
+	@${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS}
+
+clean:
+	@echo cleaning:
+	@echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+	@rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+
+.PHONY: all clean
+
diff --git a/makefile_hal b/makefile_hal
new file mode 100644
index 0000000..73cb7a8
--- /dev/null
+++ b/makefile_hal
@@ -0,0 +1,31 @@
+
+CC = clang
+CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -fopenmp
+LDFLAGS = -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lmetis -lgsl -lprofiler -ltcmalloc
+
+OBJ = break_data voronoi_bound_ini bin_values correlations beta_scales randfuncs instance get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry fracture_network get_current cracking_ini 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
+BIN = fracture current_scaling course_grain_square corr_test homo_voronoi_fracture homo_square_fracture compare_voronoi_fracture
+
+all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+
+opt:
+	@echo build options:
+	@echo "CC       = ${CC}"
+	@echo "CFLAGS   = ${CFLAGS}"
+	@echo "LDFLAGS  = ${LDFLAGS}"
+
+obj/%.o: src/%.c
+	@echo CC -c -o $@
+	@${CC} -c -o $@ $< ${CFLAGS}
+
+bin/%: obj/%.o ${OBJ:%=obj/%.o}
+	@echo CC -o $@
+	@${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS}
+
+clean:
+	@echo cleaning:
+	@echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+	@rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+
+.PHONY: all clean
+
diff --git a/src/beta_scales.c b/src/beta_scales.c
new file mode 100644
index 0000000..8cbc5d1
--- /dev/null
+++ b/src/beta_scales.c
@@ -0,0 +1,21 @@
+
+#include "fracture.h"
+
+double beta_scaling_flat(double beta, double x, double y) { return beta; }
+
+double beta_mag(double beta) {
+	double aa = -7.52579;
+	double bb = 9.63706;
+	double cc = 0.692515;
+	double dd = -2.47638;
+
+	return gsl_sf_exp(aa + bb * gsl_sf_erf(cc * (gsl_sf_log(beta) - dd)));
+}
+
+double beta_scaling_gaus(double beta, double x, double y) {
+	double sigma = 0.25;
+	double nu_f = 1.56;
+	double betap = beta * gsl_sf_exp((gsl_pow_2(x - 0.5) + gsl_pow_2(y - 0.5)) /
+																	 (2 * sigma * 2 * nu_f));
+	return betap;
+}
diff --git a/src/bin_values.c b/src/bin_values.c
new file mode 100644
index 0000000..2cb9cf6
--- /dev/null
+++ b/src/bin_values.c
@@ -0,0 +1,23 @@
+
+#include "fracture.h"
+
+double *bin_values(fnet *network, unsigned int width, double *values) {
+	double *binned = calloc(pow(width, 2), sizeof(double));
+	unsigned int *num_binned = calloc(pow(width, 2), sizeof(unsigned int));
+	for (unsigned int i = 0; i < network->num_edges; i++) {
+		if (values[i] != 0) {
+			unsigned int x = ((int)(network->edge_coords[2 * i] * width)) % width;
+			unsigned int y = ((int)(network->edge_coords[2 * i + 1] * width)) % width;
+			binned[width * x + y] += fabs(values[i]);
+			num_binned[width * x + y]++;
+		}
+	}
+
+	for (unsigned int i = 0; i < pow(width, 2); i++) {
+		if (num_binned[i] != 0) {
+			binned[i] /= num_binned[i];
+		}
+	}
+
+	return binned;
+}
diff --git a/src/break_data.c b/src/break_data.c
new file mode 100644
index 0000000..0ce6730
--- /dev/null
+++ b/src/break_data.c
@@ -0,0 +1,34 @@
+
+#include "fracture.h"
+
+break_data *alloc_break_data(unsigned int num_edges) {
+	break_data *data = malloc(1 * sizeof(break_data)); assert(data != NULL);
+
+	data->num_broken = 0;
+
+	data->break_list = (unsigned int *)malloc(num_edges * sizeof(unsigned int));
+	assert(data->break_list != NULL);
+
+	data->extern_field = (double *)malloc(num_edges * sizeof(double));
+	assert(data->extern_field != NULL);
+
+	data->conductivity = (double *)malloc(num_edges * sizeof(double));
+	assert(data->conductivity != NULL);
+
+	return data;
+}
+
+void free_break_data(break_data *data) {
+	free(data->break_list);
+	free(data->extern_field);
+	free(data->conductivity);
+	free(data);
+}
+
+void update_break_data(break_data *data, unsigned int last_broke, double strength, double conductivity) {
+	data->break_list[data->num_broken] = last_broke;
+	data->extern_field[data->num_broken] = strength;
+	data->conductivity[data->num_broken] = conductivity;
+	data->num_broken++;
+}
+
diff --git a/src/break_edge.c b/src/break_edge.c
new file mode 100644
index 0000000..0bd0277
--- /dev/null
+++ b/src/break_edge.c
@@ -0,0 +1,89 @@
+
+#include "fracture.h"
+
+bool break_edge(finst *instance, unsigned int edge, cholmod_common *c) {
+	instance->fuses[edge] = true;
+	instance->num_remaining_edges--;
+
+	unsigned int v1 = instance->network->edges_to_verts_break[2 * edge];
+	unsigned int v2 = instance->network->edges_to_verts_break[2 * edge + 1];
+
+	update_factor(instance->factor, v1, v2, c);
+
+
+	unsigned int w1 = instance->network->edges_to_verts[2 * edge];
+	unsigned int w2 = instance->network->edges_to_verts[2 * edge + 1];
+	unsigned int dw1 = instance->network->dual_edges_to_verts[2 * edge];
+	unsigned int dw2 = instance->network->dual_edges_to_verts[2 * edge + 1];
+
+	if (instance->network->boundary != TORUS_BOUND) {
+
+		unsigned int tw1 = w1 > w2 ? w1 : w2;
+		unsigned int tw2 = w1 > w2 ? w2 : w1;
+
+		CHOL_INT *lap_p = (CHOL_INT *)instance->adjacency->p;
+		CHOL_INT *lap_i = (CHOL_INT *)instance->adjacency->i;
+		double *lap_x = (double *)instance->adjacency->x;
+
+		for (int i = 0; i < lap_p[tw1 + 1] - lap_p[tw1]; i++) {
+			if (lap_i[lap_p[tw1] + i] == tw2)
+				lap_x[lap_p[tw1] + i] = 0;
+		}
+		for (int i = 0; i < lap_p[tw2 + 1] - lap_p[tw2]; i++) {
+			if (lap_i[lap_p[tw2] + i] == tw1)
+				lap_x[lap_p[tw2] + i] = 0;
+		}
+
+		int old_num_components = instance->num_components;
+
+		instance->num_components = update_components(
+				instance->adjacency, instance->marks, old_num_components, (int)tw1,
+				(int)tw2, 0);
+	}
+
+	if (instance->network->boundary == TORUS_BOUND) {
+		if (instance->dual_marks[dw1] == instance->dual_marks[dw2]) {
+			int **cycles = (int **)malloc(2*instance->network->num_edges * sizeof(int *));
+			unsigned int num_cycles = find_cycles(instance->network->num_edges, instance->fuses, instance->network->dual_edges_to_verts, instance->network->dual_verts_to_edges_ind, instance->network->dual_verts_to_edges, cycles);
+
+			for (unsigned int i = 0; i < num_cycles; i++) {
+				int x_num_crossings = 0;
+				int y_num_crossings = 0;
+				int ee; unsigned int j = 0;
+				while ((ee = cycles[2*i][j]) >= 0) {
+					int side = cycles[2*i+1][j];
+					j++;
+					unsigned int v1, v2;
+					double v1x, v1y, v2x, v2y;
+					v1 = instance->network->dual_edges_to_verts[2 * ee + !side];
+					v2 = instance->network->dual_edges_to_verts[2 * ee + side];
+					v1x = instance->network->dual_vert_coords[2 * v1];
+					v1y = instance->network->dual_vert_coords[2 * v1 + 1];
+					v2x = instance->network->dual_vert_coords[2 * v2];
+					v2y = instance->network->dual_vert_coords[2 * v2 + 1];
+					double dx = v1x - v2x;
+					double dy = v1y - v2y;
+					if (((v1x > 0.5 && v2x < 0.5) || (v1x < 0.5 && v2x > 0.5)) && fabs(dx) < 0.5) {
+						x_num_crossings += dx > 0 ? 1 : -1;
+					}
+					if (((v1y > 0.5 && v2y < 0.5) || (v1y < 0.5 && v2y > 0.5)) && fabs(dy) < 0.5) {
+						y_num_crossings += dy > 0 ? 1 : -1;
+					}
+				}
+				if ((abs(y_num_crossings) == 0 && abs(x_num_crossings) > 0) || (abs(y_num_crossings) > 0 && abs(x_num_crossings) > 0 && num_cycles > 1)) {
+					instance->num_components = 2;
+				}
+				free(cycles[2*i]);
+				free(cycles[2*i+1]);
+			}
+			free(cycles);
+		}
+
+		unsigned int *tmp_marks = get_clusters(instance, c);
+		free(instance->dual_marks);
+		instance->dual_marks = tmp_marks;
+	}
+
+
+	return true;
+}
diff --git a/src/compare_voronoi_fracture.c b/src/compare_voronoi_fracture.c
new file mode 100644
index 0000000..5578987
--- /dev/null
+++ b/src/compare_voronoi_fracture.c
@@ -0,0 +1,124 @@
+
+
+#include "fracture.h"
+
+int main(int argc, char *argv[]) {
+	int opt;
+
+	// defining variables to be (potentially) set by command line flags
+	unsigned int N, L, filename_len;
+	double beta, inf, cutoff;
+	boundary_type boundary;
+
+	filename_len = 100;
+
+	N = 100;
+	L = 16;
+	beta = .3;
+	inf = 1e10;
+	cutoff = 1e-10;
+	boundary = FREE_BOUND;
+
+	int boundary_int;
+	char boundc2 = 'f';
+
+	while ((opt = getopt(argc, argv, "n:L:b:B:dcoNsCrt")) != -1) {
+		switch (opt) {
+			case 'n':
+				N = atoi(optarg);
+				break;
+			case 'L':
+				L = atoi(optarg);
+				break;
+			case 'b':
+				beta = atof(optarg);
+				break;
+			case 'B':
+				boundary_int = atoi(optarg);
+				switch (boundary_int) {
+					case 0:
+						boundary = FREE_BOUND;
+						boundc2 = 'f';
+						break;
+					case 1:
+						boundary = CYLINDER_BOUND;
+						boundc2 = 'c';
+						break;
+					default:
+						printf("boundary specifier must be 0 (FREE_BOUND) or 1 (CYLINDER_BOUND).\n");
+				}
+				break;
+			default: /* '?' */
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	char *break_filename = (char *)malloc(filename_len * sizeof(char));
+	snprintf(break_filename, filename_len, "breaks_v_vc_%c_%u_%g.txt", boundc2, L, beta);
+	FILE *break_out = fopen(break_filename, "a");
+	free(break_filename);
+
+
+	// start cholmod
+	cholmod_common c;
+	CHOL_F(start)(&c);
+
+	(&c)->supernodal = CHOLMOD_SIMPLICIAL;
+
+
+	fnet *network = ini_voronoi_network(L, boundary, genfunc_hyperuniform, &c);
+	finst *perm_voltage_instance = create_instance(network, inf, true, true, &c);
+	finst *perm_current_instance = create_instance(network, inf, false, true, &c);
+	double *fuse_thres = gen_fuse_thres(network->num_edges, network->edge_coords, beta, beta_scaling_flat);
+	finst *voltage_instance = copy_instance(perm_voltage_instance, &c);
+	finst *current_instance = copy_instance(perm_current_instance, &c);
+	break_data *breaking_data_voltage = fracture_network(voltage_instance, fuse_thres, &c, cutoff);
+	break_data *breaking_data_current = fracture_network(current_instance, fuse_thres, &c, cutoff);
+	free_instance(voltage_instance, &c);
+	free_instance(current_instance, &c);
+	free_instance(perm_voltage_instance, &c);
+	free_instance(perm_current_instance, &c);
+	free(fuse_thres);
+
+	FILE *net_out = fopen("network.txt", "w");
+	for (unsigned int j = 0; j < network->num_verts; j++) {
+		fprintf(net_out, "%f %f ", network->vert_coords[2 * j],
+						network->vert_coords[2 * j + 1]);
+	}
+	fprintf(net_out, "\n");
+	for (unsigned int j = 0; j < network->num_edges; j++) {
+		fprintf(net_out, "%u %u ", network->edges_to_verts[2 * j],
+						network->edges_to_verts[2 * j + 1]);
+	}
+	fprintf(net_out, "\n");
+	for (unsigned int j = 0; j < network->num_dual_verts; j++) {
+		fprintf(net_out, "%f %f ", network->dual_vert_coords[2 * j],
+						network->dual_vert_coords[2 * j + 1]);
+	}
+	fprintf(net_out, "\n");
+	for (unsigned int j = 0; j < network->num_edges; j++) {
+		fprintf(net_out, "%u %u ", network->dual_edges_to_verts[2 * j],
+						network->dual_edges_to_verts[2 * j + 1]);
+	}
+
+	free_fnet(network, &c);
+
+	for (unsigned int j = 0; j < breaking_data_voltage->num_broken; j++) {
+		fprintf(break_out, "%u %g %g ", breaking_data_voltage->break_list[j],
+						breaking_data_voltage->extern_field[j], breaking_data_voltage->conductivity[j]);
+	}
+	fprintf(break_out, "\n");
+	for (unsigned int j = 0; j < breaking_data_current->num_broken; j++) {
+		fprintf(break_out, "%u %g %g ", breaking_data_current->break_list[j],
+						breaking_data_current->extern_field[j], breaking_data_current->conductivity[j]);
+	}
+	fprintf(break_out, "\n");
+
+	free_break_data(breaking_data_voltage);
+	free_break_data(breaking_data_current);
+	fclose(break_out);
+
+	CHOL_F(finish)(&c);
+
+	return 0;
+}
diff --git a/src/corr_test.c b/src/corr_test.c
new file mode 100644
index 0000000..2c93b8e
--- /dev/null
+++ b/src/corr_test.c
@@ -0,0 +1,30 @@
+
+#include "fracture.h"
+
+int main() {
+	cholmod_common c;
+	CHOL_F(start)(&c);
+
+	unsigned int width = 64;
+	unsigned int n = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2);
+
+	fnet *network = ini_square_network(width, true, false, &c);
+	finst *instance = create_instance(network, 1e-14, true, true, &c);
+	double *fuse_thres = gen_fuse_thres(network->num_edges, network->edge_coords,
+																			0.001, beta_scaling_flat);
+	fracture_network(instance, fuse_thres, &c, 1e-10);
+
+	double *corr = get_corr(instance, NULL, &c);
+
+	for (int i = 0; i < 2 * width; i++) {
+		printf("%f ", corr[i]);
+	}
+	printf("\n");
+
+	free_instance(instance, &c);
+	free_fnet(network, &c);
+
+	CHOL_F(finish)(&c);
+
+	return 0;
+}
diff --git a/src/correlations.c b/src/correlations.c
new file mode 100644
index 0000000..bd05daf
--- /dev/null
+++ b/src/correlations.c
@@ -0,0 +1,345 @@
+
+#include "fracture.h"
+
+struct fibn {
+	unsigned int value;
+	unsigned int priority;
+	unsigned int rank;
+	bool marked;
+	struct fibn *first_child;
+	struct fibn *parent;
+	struct fibn *prev;
+	struct fibn *next;
+};
+
+typedef struct {
+	struct fibn *min;
+	unsigned int rank;
+	unsigned int trees;
+	struct fibn **vertex_to_node;
+} fibh;
+
+unsigned int fib_findmin(fibh *heap) { return heap->min->value; }
+
+void ll_setparent(struct fibn *ll, struct fibn *parent) {
+	struct fibn *curnode = ll->next;
+	ll->parent = parent;
+	while (curnode != ll) {
+		curnode->parent = parent;
+		curnode = curnode->next;
+	}
+}
+
+struct fibn *ll_merge(struct fibn *ll1, struct fibn *ll2) {
+	if (ll1 == NULL)
+		return ll2;
+	if (ll2 == NULL)
+		return ll1;
+
+	// link the beginning of list one to the end of list two and vice versa
+	struct fibn *ll1_beg = ll1;
+	struct fibn *ll1_end = ll1->prev;
+	struct fibn *ll2_beg = ll2;
+	struct fibn *ll2_end = ll2->prev;
+
+	ll1_beg->prev = ll2_end;
+	ll1_end->next = ll2_beg;
+	ll2_beg->prev = ll1_end;
+	ll2_end->next = ll1_beg;
+
+	return ll1;
+}
+
+struct fibn *ll_delroot(struct fibn *ll) {
+	if (ll == NULL)
+		return NULL;
+
+	if (ll->next == ll) {
+		return NULL;
+	}
+
+	struct fibn *ll_beg = ll->next;
+	struct fibn *ll_end = ll->prev;
+
+	ll_beg->prev = ll_end;
+	ll_end->next = ll_beg;
+
+	ll->next = ll;
+	ll->prev = ll;
+
+	return ll_beg;
+}
+
+void fib_insert(fibh *heap, unsigned int value, unsigned int priority) {
+	struct fibn *newnode = calloc(1, sizeof(struct fibn));
+	newnode->value = value;
+	newnode->priority = priority;
+	newnode->next = newnode;
+	newnode->prev = newnode;
+
+	heap->min = ll_merge(heap->min, newnode);
+
+	if (priority < heap->min->priority) {
+		heap->min = newnode;
+	}
+
+	heap->vertex_to_node[value] = newnode;
+
+	heap->trees++;
+}
+
+void fib_deletemin(fibh *heap) {
+	unsigned int min_rank = heap->min->rank;
+	struct fibn *min_children = heap->min->first_child;
+
+	struct fibn *trees = ll_delroot(heap->min);
+	heap->vertex_to_node[heap->min->value] = NULL;
+	free(heap->min);
+
+	if (min_children != NULL)
+		ll_setparent(min_children, NULL);
+	trees = ll_merge(trees, min_children);
+
+	heap->trees += min_rank - 1;
+
+	if (trees == NULL) {
+		// min had no children and was only tree, return empty heap
+		heap->min = NULL;
+		heap->rank = 0;
+	} else {
+		// min had children or there were other trees
+		unsigned int min_val = UINT_MAX;
+		struct fibn *min_ptr = NULL;
+		unsigned int max_rank = 0;
+		struct fibn *curnode = trees;
+		for (unsigned int i = 0; i < heap->trees; i++) {
+			if (curnode->priority < min_val) {
+				min_ptr = curnode;
+				min_val = curnode->priority;
+			}
+			if (curnode->rank > max_rank)
+				max_rank = curnode->rank;
+			curnode = curnode->next;
+		}
+
+		if (min_ptr == NULL)
+			min_ptr = trees;
+		heap->min = min_ptr;
+		heap->rank = max_rank;
+
+		struct fibn **rankslots = calloc(max_rank + 1, sizeof(struct fibn *));
+		curnode = heap->min;
+		while (curnode != rankslots[curnode->rank]) {
+			if (rankslots[curnode->rank] == NULL) {
+				rankslots[curnode->rank] = curnode;
+				curnode = curnode->next;
+			} else {
+				struct fibn *oldnode = rankslots[curnode->rank];
+				rankslots[curnode->rank] = NULL;
+				struct fibn *smaller =
+						curnode->priority < oldnode->priority || curnode == heap->min
+								? curnode
+								: oldnode;
+				struct fibn *larger = smaller == curnode ? oldnode : curnode;
+				ll_delroot(larger);
+				ll_setparent(larger, smaller);
+				struct fibn *smaller_children = smaller->first_child;
+				smaller->first_child = ll_merge(smaller_children, larger);
+				heap->trees--;
+				smaller->rank++;
+				if (smaller->rank > heap->rank) {
+					heap->rank = smaller->rank;
+					rankslots =
+							realloc(rankslots, (heap->rank + 1) * sizeof(struct fibn *));
+					rankslots[heap->rank] = NULL;
+				}
+				curnode = smaller;
+			}
+		}
+		free(rankslots);
+	}
+}
+
+void fib_decreasekey(fibh *heap, unsigned int value,
+										 unsigned int new_priority) {
+	struct fibn *node = heap->vertex_to_node[value];
+	if (node != NULL) {
+		node->priority = new_priority;
+		if (node->parent != NULL) {
+			if (node->priority < node->parent->priority) {
+
+				struct fibn *curnode = node;
+				curnode->marked = true;
+				while (curnode->parent != NULL && curnode->marked) {
+					struct fibn *oldparent = curnode->parent;
+					oldparent->rank--;
+					oldparent->first_child = ll_delroot(curnode);
+					ll_setparent(curnode, NULL);
+					ll_merge(heap->min, curnode);
+					heap->trees++;
+					if (curnode->marked) {
+						curnode->marked = false;
+					}
+
+					if (oldparent->marked) {
+						curnode = oldparent;
+					} else {
+						oldparent->marked = true;
+						break;
+					}
+				}
+			}
+		}
+
+		if (new_priority < heap->min->priority) {
+			heap->min = node;
+		}
+	}
+}
+
+unsigned int *dijkstra(fnet *network, unsigned int source) {
+	unsigned int nv = network->num_dual_verts;
+	unsigned int *vei = network->dual_verts_to_edges_ind;
+	unsigned int *ve = network->dual_verts_to_edges;
+	unsigned int *ev = network->dual_edges_to_verts;
+
+	unsigned int *dist = (unsigned int *)calloc(nv, sizeof(unsigned int));
+	fibh *Q = (fibh *)calloc(1, sizeof(fibh));
+	Q->vertex_to_node = (struct fibn **)calloc(nv, sizeof(struct fibn *));
+
+	for (unsigned int i = 0; i < nv; i++) {
+		if (i != source) {
+			dist[i] = UINT_MAX;
+		}
+
+		fib_insert(Q, i, dist[i]);
+	}
+
+	while (Q->min != NULL) {
+		unsigned int u = fib_findmin(Q);
+		fib_deletemin(Q);
+
+		for (unsigned int i = 0; i < vei[u + 1] - vei[u]; i++) {
+			unsigned int e = ve[vei[u] + i];
+			unsigned int v = ev[2 * e] == u ? ev[2 * e + 1] : ev[2 * e];
+			unsigned int alt = dist[u] + 1;
+			if (alt < dist[v]) {
+				dist[v] = alt;
+				fib_decreasekey(Q, v, alt);
+			}
+		}
+	}
+
+	free(Q->vertex_to_node);
+	free(Q);
+
+	return dist;
+}
+
+unsigned int **get_dists(fnet *network) {
+	unsigned int nv = network->num_dual_verts;
+	unsigned int **dists = (unsigned int **)malloc(nv * sizeof(unsigned int *));
+
+#pragma omp parallel for
+	for (unsigned int i = 0; i < nv; i++) {
+		dists[i] = dijkstra(network, i);
+	}
+
+	return dists;
+}
+
+double *get_corr(finst *instance, unsigned int **dists, cholmod_common *c) {
+	unsigned int nv = instance->network->num_dual_verts;
+	unsigned int ne = instance->network->num_edges;
+	unsigned int *ev = instance->network->dual_edges_to_verts;
+	bool nulldists = false;
+	if (dists == NULL) {
+		dists = get_dists(instance->network);
+		nulldists = true;
+	}
+	double *corr = calloc(nv, sizeof(double));
+	unsigned int *marks = get_clusters(instance, c);
+	unsigned int *numat = calloc(nv, sizeof(unsigned int));
+
+	for (unsigned int i = 0; i < ne; i++) {
+		unsigned int v1 = ev[2 * i];
+		unsigned int v2 = ev[2 * i + 1];
+		for (unsigned int j = 0; j < ne; j++) {
+			unsigned int v3 = ev[2 * j];
+			unsigned int v4 = ev[2 * j + 1];
+			unsigned int dist1 = dists[v1][v3];
+			unsigned int dist2 = dists[v1][v4];
+			unsigned int dist3 = dists[v2][v3];
+			unsigned int dist4 = dists[v2][v4];
+			unsigned int dist = (dist1 + dist2 + dist3 + dist4) / 4;
+			corr[dist] += instance->fuses[i] && instance->fuses[j];
+			numat[dist]++;
+		}
+	}
+
+	for (unsigned int i = 0; i < nv; i++) {
+		if (numat[i] > 0) {
+			corr[i] /= numat[i];
+		}
+	}
+
+	if (nulldists) {
+		for (int i = 0; i < nv; i++) {
+			free(dists[i]);
+		}
+		free(dists);
+	}
+
+	free(marks);
+
+	return corr;
+}
+
+/*
+double *get_space_corr(finst *instance, cholmod_common *c) {
+	unsigned int nv = instance->network->num_dual_verts;
+	double *vc = instance->network->dual_vert_coords;
+	unsigned int ne = instance->network->num_edges;
+	unsigned int *ev = instance->network->dual_edges_to_verts;
+	double *corr = calloc(nv, sizeof(unsigned int));
+	unsigned int *numat = calloc(nv, sizeof(unsigned int));
+
+	for (unsigned int i = 0; i < ne; i++) {
+		unsigned int v1 = ev[2 * i];
+		unsigned int v2 = ev[2 * i + 1];
+		double v1x = vc[2 * v1];
+		double v1y = vc[2 * v1 + 1];
+		double v2x = vc[2 * v2];
+		double v2y = vc[2 * v2 + 1];
+		double e1x = (v1x + v2x) / 2;
+		double e1y = (v1y + v2y) / 2;
+		for (unsigned int j = 0; j < ne; j++) {
+			unsigned int v3 = ev[2 * j];
+			unsigned int v4 = ev[2 * j + 1];
+			double v1x = vc[2 * v1];
+			double v1y = vc[2 * v1 + 1];
+			double v2x = vc[2 * v2];
+			double v2y = vc[2 * v2 + 1];
+			double e1x = (v1x + v2x) / 2;
+			double e1y = (v1y + v2y) / 2;
+			corr[dist] += instance->fuses[i] && instance->fuses[j];
+			numat[dist]++;
+		}
+	}
+
+	for (unsigned int i = 0; i < nv; i++) {
+		if (numat[i] > 0) {
+			corr[i] /= numat[i];
+		}
+	}
+
+	for (int i = 0; i < nv; i++) {
+		free(dists[i]);
+	}
+
+	free(marks);
+	free(dists);
+
+	return corr;
+}
+*/
diff --git a/src/course_grain_square.c b/src/course_grain_square.c
new file mode 100644
index 0000000..9035cf7
--- /dev/null
+++ b/src/course_grain_square.c
@@ -0,0 +1,168 @@
+
+#include "fracture.h"
+
+int main(int argc, char *argv[]) {
+	int opt;
+
+	// defining variables to be (potentially) set by command line flags
+	int num = 100;
+	int width = 16;
+	double beta = .3;
+	bool save_clusters = false;
+	bool voltage_bound = false;
+	bool output_break_data = false;
+
+	while ((opt = getopt(argc, argv, "n:w:b:cVo")) != -1) {
+		switch (opt) {
+		case 'n':
+			num = atoi(optarg);
+			break;
+		case 'w':
+			width = atoi(optarg);
+			break;
+		case 'b':
+			beta = atof(optarg);
+			break;
+		case 'c':
+			save_clusters = true;
+			break;
+		case 'V':
+			voltage_bound = true;
+			break;
+		case 'o':
+			output_break_data = true;
+			break;
+		default: /* '?' */
+			exit(EXIT_FAILURE);
+		}
+	}
+
+	FILE *break_out;
+	if (output_break_data) {
+		char *break_filename = (char *)malloc(100 * sizeof(char));
+		snprintf(break_filename, 100, "breaks_%d_%.3f_%d.txt", width, beta, num);
+		break_out = fopen(break_filename, "w");
+		free(break_filename);
+	}
+
+	bool periodic = true;
+	double inf = 1;
+	double cutoff = 1e-10;
+
+	// 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 (voltage_bound) {
+		(&c)->supernodal = CHOLMOD_SUPERNODAL;
+	} else {
+		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
+	}
+
+	fnet *network = ini_square_network(width, periodic, false, &c);
+	fnet *network_p = ini_square_network(width / 2, periodic, false, &c);
+	finst *perm_instance = create_instance(network, inf, voltage_bound, true, &c);
+	unsigned int c_dist_size = network->num_dual_verts;
+	unsigned int c_p_dist_size = network_p->num_dual_verts;
+
+	// define arrays for saving cluster and avalanche distributions
+	unsigned int *cluster_size_dist =
+			(unsigned int *)calloc(c_dist_size, sizeof(unsigned int));
+	unsigned int *cluster_p_size_dist =
+			(unsigned int *)calloc(c_p_dist_size, sizeof(unsigned int));
+
+	printf("\n");
+	for (int DUMB = 0; DUMB < num; DUMB++) {
+		printf("\033[F\033[JCOURSEGRAIN_SQUARE: %0*d / %d\n", (int)log10(num) + 1,
+					 DUMB + 1, num);
+
+		break_data *breaking_data = NULL;
+		finst *instance = NULL;
+		while (breaking_data == NULL) {
+			double *fuse_thres = gen_fuse_thres(
+					network->num_edges, network->edge_coords, beta, beta_scaling_flat);
+			instance = copy_instance(perm_instance, &c);
+			breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
+			free_instance(instance, &c);
+			free(fuse_thres);
+		}
+
+		unsigned int min_pos = 0;
+		double min_val = DBL_MAX;
+
+		for (unsigned int j = 0; j < breaking_data->num_broken; j++) {
+			double val = fabs(breaking_data->extern_field[j]);
+			if (val < min_val) {
+				min_pos = j;
+				min_val = val;
+			}
+			if (val > 10 * min_val)
+				break;
+		}
+
+		finst *tmp_instance = copy_instance(perm_instance, &c);
+
+		for (unsigned int i = 0; i < breaking_data->num_broken; i++) {
+			break_edge(tmp_instance, breaking_data->break_list[i], &c);
+		}
+
+		unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c);
+		for (unsigned int i = 0; i < network->num_dual_verts; i++) {
+			cluster_size_dist[i] += tmp_cluster_dist[i];
+		}
+		free(tmp_cluster_dist);
+
+		finst *instance_p = coursegrain_square(tmp_instance, network_p, &c);
+
+		unsigned int *tmp_cluster_p_dist = get_cluster_dist(instance_p, &c);
+		for (unsigned int i = 0; i < network_p->num_dual_verts; i++) {
+			cluster_p_size_dist[i] += tmp_cluster_p_dist[i];
+		}
+		free(tmp_cluster_p_dist);
+
+		free_instance(tmp_instance, &c);
+		free_instance(instance_p, &c);
+		if (output_break_data) {
+			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]);
+			}
+			fprintf(break_out, "\n");
+		}
+		free(breaking_data->break_list);
+		free(breaking_data->extern_field);
+		free(breaking_data);
+	}
+
+	printf("\033[F\033[JCURRENT_SCALING: COMPLETE");
+
+	if (save_clusters) {
+		FILE *cluster_out = get_file("cluster", width, 0, beta, 1, 1, num, false);
+		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 < c_p_dist_size; i++) {
+			fprintf(cluster_out, "%u ", cluster_p_size_dist[i]);
+		}
+		fclose(cluster_out);
+	}
+
+	if (output_break_data) {
+		fclose(break_out);
+	}
+
+	free(cluster_size_dist);
+	free(cluster_p_size_dist);
+	free_instance(perm_instance, &c);
+	free_fnet(network, &c);
+	free_fnet(network_p, &c);
+
+	CHOL_F(finish)(&c);
+
+	return 0;
+}
diff --git a/src/coursegrain.c b/src/coursegrain.c
new file mode 100644
index 0000000..bd84e1c
--- /dev/null
+++ b/src/coursegrain.c
@@ -0,0 +1,40 @@
+
+#include "fracture.h"
+
+finst *coursegrain_square(finst *instance, fnet *network_p, cholmod_common *c) {
+	unsigned int width = sqrt(instance->network->num_edges);
+	assert(width % 4 == 0);
+
+	finst *instance_p = create_instance(network_p, instance->inf,
+																			instance->voltage_bound, true, c);
+
+	unsigned int width_p = width / 2;
+	bool *fuses = instance->fuses;
+
+	for (unsigned int i = 0; i < network_p->num_edges; i++) {
+		int xp = i / width_p;
+		int yp = i % width_p;
+		unsigned int x1 = 2 * xp;
+		unsigned int y1 = (2 * yp - 1) % width;
+		unsigned int x2 = 2 * xp;
+		unsigned int y2 = 2 * yp;
+		unsigned int x3 = 2 * xp + 1;
+		unsigned int y3 = (2 * yp - 1) % width;
+		unsigned int x4 = 2 * xp + 1;
+		unsigned int y4 = 2 * yp;
+		bool f1 = fuses[width * x1 + y1];
+		bool f2 = fuses[width * x2 + y2];
+		bool f3 = fuses[width * x3 + y3];
+		bool f4 = fuses[width * x4 + y4];
+
+		if ((f1 && f2) || (f3 && f4)) {
+			//		instance_p->fuses[i] = true;
+			//		instance_p->num_remaining_edges--;
+			break_edge(instance_p, i, c);
+		}
+	}
+
+	// fin_instance(instance_p, c);
+
+	return instance_p;
+}
diff --git a/src/cracking_ini.c b/src/cracking_ini.c
new file mode 100644
index 0000000..9ab4a9d
--- /dev/null
+++ b/src/cracking_ini.c
@@ -0,0 +1,53 @@
+
+#include "fracture.h"
+
+double *gen_fuse_thres(unsigned int num_edges, double *edge_coords, double beta,
+											 double (*beta_scale)(double, double, double)) {
+	unsigned int size = num_edges;
+	assert(beta > 0);
+
+	double *fuse_thres = (double *)malloc(size * sizeof(double));
+	assert(fuse_thres != NULL);
+
+	gsl_set_error_handler_off();
+
+	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);
+
+	for (unsigned int i = 0; i < size; i++) {
+		double x = edge_coords[2 * i];
+		double y = edge_coords[2 * i + 1];
+		while ((fuse_thres[i] = gsl_sf_exp(gsl_sf_log(gsl_ran_flat(r, 0, 1)) /
+																			 beta_scale(beta, x, y))) == 0.0)
+			;
+	}
+
+	gsl_rng_free(r);
+
+	return fuse_thres;
+}
+
+bool gen_crack(finst *instance, double crack_len, double crack_width,
+							 cholmod_common *c) {
+	assert(instance != NULL);
+	bool *fuses = instance->fuses;
+	assert(fuses != NULL);
+	const fnet *network = instance->network;
+	unsigned int num_edges = network->num_edges;
+	double *edge_coords = network->edge_coords;
+
+	for (unsigned int j = 0; j < num_edges; j++) {
+		if (edge_coords[2 * j + 1] < crack_len &&
+				fabs(edge_coords[2 * j] - network->L / 2) < crack_width) {
+			instance->fuses[j] = true;
+			instance->num_remaining_edges--;
+		} else
+			fuses[j] = false;
+	}
+
+	return true;
+}
diff --git a/src/current_scaling.c b/src/current_scaling.c
new file mode 100644
index 0000000..6837f6b
--- /dev/null
+++ b/src/current_scaling.c
@@ -0,0 +1,300 @@
+
+#include "fracture.h"
+
+int main(int argc, char *argv[]) {
+	int opt;
+
+	// defining variables to be (potentially) set by command line flags
+	int iter = 1;
+	int num = 100;
+	int width = 16;
+	double crack_len = 8;
+	double beta = .3;
+	double inf = 1e-8;
+	double cutoff = 1e-8;
+	bool beta_shift = false;
+	bool supplied_bound = false;
+	bool ash_beta = false;
+	char *bound_file;
+	bool voltage_bound = false;
+	bool use_first = false;
+	bool save_stress = false;
+	bool save_bound = false;
+	bool save_damage = false;
+	bool save_strength = false;
+
+	while ((opt = getopt(argc, argv, "n:w:b:l:i:Bf:aVFsSed")) != -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 'i':
+			iter = atoi(optarg);
+			break;
+		case 'B':
+			beta_shift = true;
+			break;
+		case 'a':
+			ash_beta = true;
+			break;
+		case 'V':
+			voltage_bound = true;
+			break;
+		case 'F':
+			use_first = true;
+			break;
+		case 's':
+			save_stress = true;
+			break;
+		case 'd':
+			save_damage = true;
+			break;
+		case 'e':
+			save_bound = true;
+			break;
+		case 'S':
+			save_strength = true;
+			break;
+		case 'f':
+			supplied_bound = true;
+			bound_file = optarg;
+			break;
+		default: /* '?' */
+			exit(EXIT_FAILURE);
+		}
+	}
+
+	// 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 (voltage_bound) {
+		(&c)->supernodal = CHOLMOD_SUPERNODAL;
+	} else {
+		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
+	}
+
+	fnet *network = ini_square_network(width, false, true, &c);
+	finst *perm_instance =
+			create_instance(network, inf, voltage_bound, false, &c);
+	gen_crack(perm_instance, crack_len, 1, &c);
+	finish_instance(perm_instance, &c);
+
+	if (voltage_bound) {
+		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
+		finst *tmp_instance = create_instance(network, inf, false, false, &c);
+		gen_crack(tmp_instance, crack_len, 1, &c);
+		finish_instance(tmp_instance, &c);
+		double *voltage = get_voltage(tmp_instance, &c);
+
+		for (int i = 0; i < network->num_bounds; i++) {
+			for (int j = 0; j < network->bound_inds[i + 1] - network->bound_inds[i];
+					 j++) {
+				((double *)perm_instance->boundary_cond
+						 ->x)[network->bound_verts[network->bound_inds[i] + j]] =
+						voltage[network->bound_verts[network->bound_inds[i] + j]];
+			}
+		}
+		(&c)->supernodal = CHOLMOD_SUPERNODAL;
+	}
+
+	if (supplied_bound) {
+		FILE *bound_f = fopen(bound_file, "r");
+		for (int i = 0; i < network->num_verts; i++) {
+			double tmp;
+			fscanf(bound_f, "%lg ", &tmp);
+			((double *)perm_instance->boundary_cond->x)[i] = tmp;
+		}
+
+		((double *)perm_instance->boundary_cond->x)[network->num_verts] = 0;
+		((double *)perm_instance->boundary_cond->x)[network->num_verts + 1] = 0;
+		fclose(bound_f);
+	}
+
+	printf("\n");
+	for (int DUMB2 = 0; DUMB2 < iter; DUMB2++) {
+
+		double *strength;
+		if (save_strength) {
+			strength = (double *)malloc(num * sizeof(double));
+		}
+
+		double *damage;
+		if (save_damage) {
+			damage = (double *)calloc(network->num_edges, sizeof(double));
+		}
+		double *avg_current = (double *)calloc(network->num_edges, sizeof(double));
+		unsigned int *num_current_skipped =
+				(unsigned int *)calloc(network->num_edges, sizeof(unsigned int));
+		double *avg_voltage = (double *)calloc(network->num_verts, sizeof(double));
+		unsigned int *num_voltage_skipped =
+				(unsigned int *)calloc(network->num_verts, sizeof(unsigned int));
+
+		for (int DUMB = 0; DUMB < num; DUMB++) {
+			printf("\033[F\033[JCURRENT_SCALING: ITERATION %0*d: %0*d / %d\n",
+						 (int)log10(iter) + 1, DUMB2 + 1, (int)log10(num) + 1, DUMB + 1,
+						 num);
+
+			break_data *breaking_data = NULL;
+			while (breaking_data == NULL) {
+				double *fuse_thres = gen_fuse_thres(
+						network->num_edges, network->edge_coords, 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);
+			}
+
+			unsigned int min_pos = 0;
+			double min_val = DBL_MAX;
+
+			for (unsigned int j = 0; j < breaking_data->num_broken; j++) {
+				double val = fabs(breaking_data->extern_field[j]);
+				if (val < min_val) {
+					min_pos = j;
+					min_val = val;
+				}
+			}
+
+			if (save_strength) {
+				strength[DUMB] = fabs(breaking_data->extern_field[min_pos]);
+			}
+
+			finst *tmp_instance = copy_instance(perm_instance, &c);
+
+			unsigned int until = min_pos;
+			if (use_first)
+				until = 1;
+			for (unsigned int i = 0; i < until; i++) {
+				break_edge(tmp_instance, breaking_data->break_list[i], &c);
+				if (save_damage) {
+					damage[breaking_data->break_list[i]] += 1. / num;
+				}
+			}
+
+			double *voltage = get_voltage(tmp_instance, &c);
+			double *current = get_current(tmp_instance, &c);
+
+			for (unsigned int i = 0; i < network->num_edges; i++) {
+				avg_current[i] += current[i];
+				if (current[i] == 0)
+					num_current_skipped[i]++;
+			}
+
+			for (unsigned int i = 0; i < network->num_verts; i++) {
+				if (tmp_instance->marks[i] == tmp_instance->marks[network->num_verts]) {
+					avg_voltage[i] += voltage[i];
+				} else {
+					num_voltage_skipped[i]++;
+				}
+			}
+
+			free(current);
+			free(voltage);
+			free_instance(tmp_instance, &c);
+			free(breaking_data->break_list);
+			free(breaking_data->extern_field);
+			free(breaking_data);
+		}
+
+		for (int i = 0; i < network->num_edges; i++) {
+			if (num_current_skipped[i] < num) {
+				avg_current[i] /= num - num_current_skipped[i];
+			}
+		}
+
+		for (int i = 0; i < network->num_verts; i++) {
+			if (num_voltage_skipped[i] < num) {
+				avg_voltage[i] /= num - num_voltage_skipped[i];
+			}
+		}
+
+		double *avg_field;
+		if (voltage_bound)
+			avg_field = avg_voltage;
+		else
+			avg_field = avg_current;
+
+		update_boundary(perm_instance, avg_field);
+
+		if (save_stress) {
+			char *c_filename = (char *)malloc(100 * sizeof(char));
+			snprintf(c_filename, 100, "current_%d_%g_%d_%g.txt", width, crack_len,
+							 iter, beta);
+			FILE *outfile = fopen(c_filename, "w");
+			for (int i = 0; i < network->num_edges; i++) {
+				fprintf(outfile, "%g ", avg_current[i]);
+			}
+			fclose(outfile);
+			free(c_filename);
+		}
+
+		if (save_damage) {
+			char *c_filename = (char *)malloc(100 * sizeof(char));
+			snprintf(c_filename, 100, "damage_%d_%g_%d_%g.txt", width, crack_len,
+							 iter, beta);
+			FILE *outfile = fopen(c_filename, "w");
+			for (int i = 0; i < network->num_edges; i++) {
+				fprintf(outfile, "%g ", damage[i]);
+			}
+			fclose(outfile);
+			free(c_filename);
+		}
+
+		if (save_strength) {
+			char *s_filename = (char *)malloc(100 * sizeof(char));
+			snprintf(s_filename, 100, "strength_%d_%g_%d_%g.txt", width, crack_len, iter, beta);
+			FILE *f = fopen(s_filename, "a");
+			for (int i = 0; i < num; i++) {
+				fprintf(f, "%g ", strength[i]);
+			}
+			fclose(f);
+			free(s_filename);
+		}
+
+		if (save_bound) {
+			char *b_filename = (char *)malloc(100 * sizeof(char));
+			snprintf(b_filename, 100, "bounds_%d_%g_%d_%g.txt", width, crack_len,
+							 iter, beta);
+			FILE *outfile = fopen(b_filename, "w");
+			for (int i = 0; i < network->num_verts; i++) {
+				fprintf(outfile, "%g ", ((double *)perm_instance->boundary_cond->x)[i]);
+			}
+			fclose(outfile);
+			free(b_filename);
+		}
+
+		free(avg_current);
+		free(avg_voltage);
+		if (save_damage) free(damage);
+		free(num_current_skipped);
+		free(num_voltage_skipped);
+		if (save_strength) {
+			free(strength);
+		}
+
+		printf(
+				"\033[F\033[JCURRENT_SCALING: ITERATION %0*d COMPLETE, BETA = %.2f\n\n",
+				(int)log10(iter) + 1, DUMB2 + 1, beta);
+	}
+
+	free_instance(perm_instance, &c);
+	free_fnet(network, &c);
+
+	CHOL_F(finish)(&c);
+
+	return 0;
+}
diff --git a/src/fortune/defs.h b/src/fortune/defs.h
new file mode 100644
index 0000000..b186581
--- /dev/null
+++ b/src/fortune/defs.h
@@ -0,0 +1,134 @@
+#ifndef NULL
+#define NULL 0
+#endif
+#define DELETED -2
+#include <stdlib.h>
+#include <limits.h>
+#include <stdbool.h>
+
+extern int triangulate, sorted, plot, debug;
+
+struct Freenode {
+	struct Freenode *nextfree;
+};
+struct Freelist {
+	struct Freenode *head;
+	int nodesize;
+};
+char *getfree();
+char *myalloc();
+
+extern double xmin, xmax, ymin, ymax, deltax, deltay;
+
+struct Point {
+	double x, y;
+};
+
+/* structure used both for sites and for vertices */
+struct Site {
+	struct Point coord;
+	int sitenbr;
+	int refcnt;
+};
+
+extern struct Site *sites;
+extern int nsites;
+extern int siteidx;
+extern int sqrt_nsites;
+extern int nvertices;
+extern struct Freelist sfl;
+extern struct Site *bottomsite;
+
+char **alloclist;
+int allocnum;
+
+struct Edge {
+	double a, b, c;
+	struct Site *ep[2];
+	struct Site *reg[2];
+	int edgenbr;
+};
+#define le 0
+#define re 1
+extern int nedges;
+extern struct Freelist efl;
+
+int has_endpoint(), right_of();
+struct Site *intersect();
+double dist();
+struct Point PQ_min();
+struct Halfedge *PQextractmin();
+struct Edge *bisect();
+
+struct Halfedge {
+	struct Halfedge *ELleft, *ELright;
+	struct Edge *ELedge;
+	int ELrefcnt;
+	char ELpm;
+	struct Site *vertex;
+	double ystar;
+	struct Halfedge *PQnext;
+};
+
+extern struct Freelist hfl;
+extern struct Halfedge *ELleftend, *ELrightend;
+extern int ELhashsize;
+extern struct Halfedge **ELhash;
+struct Halfedge *HEcreate(), *ELleft(), *ELright(), *ELleftbnd();
+struct Site *leftreg(), *rightreg();
+
+extern int PQhashsize;
+extern struct Halfedge *PQhash;
+extern struct Halfedge *PQfind();
+extern int PQcount;
+extern int PQmin;
+int PQempty();
+
+extern double *site_list;
+extern double *vert_list;
+extern double *line_list;
+extern unsigned int *edge_list;
+extern unsigned int *dual_list;
+
+extern int vert_count;
+extern int edge_count;
+extern int line_count;
+extern int dual_count;
+extern int site_count;
+
+void makefree(struct Freenode *curr, struct Freelist *fl);
+void voronoi(int triangulate, struct Site *(*nextsite)());
+char *getfree(struct Freelist *fl);
+char *myalloc(unsigned n);
+void freeinit(struct Freelist *fl, int size);
+void ELinitialize();
+struct Halfedge *HEcreate(struct Edge *e, int pm);
+void ELinsert(struct Halfedge *lb, struct Halfedge *new);
+struct Halfedge *ELgethash(int b);
+struct Halfedge *ELleftbnd(struct Point *p);
+void ELdelete(struct Halfedge *he);
+struct Halfedge *ELright(struct Halfedge *he);
+struct Halfedge *ELleft(struct Halfedge *he);
+struct Site *leftreg(struct Halfedge *he);
+struct Site *rightreg(struct Halfedge *he);
+void geominit();
+struct Edge *bisect(struct Site *s1, struct Site *s2);
+struct Site *intersect(struct Halfedge *el1, struct Halfedge *el2);
+int right_of(struct Halfedge *el, struct Point *p);
+void endpoint(struct Edge *e, int lr, struct Site *s);
+double dist(struct Site *s, struct Site *t);
+void makevertex(struct Site *v);
+void deref(struct Site *v);
+void ref(struct Site *v);
+void PQinsert(struct Halfedge *he, struct Site *v, double offset);
+void PQdelete(struct Halfedge *he);
+int PQbucket(struct Halfedge *he);
+int PQempty();
+struct Point PQ_min();
+struct Halfedge *PQextractmin();
+void PQinitialize();
+void out_bisector(struct Edge *e);
+void out_ep(struct Edge *e);
+void out_vertex(struct Site *v);
+void out_site(struct Site *s);
+void out_triple(struct Site *s1, struct Site *s2, struct Site *s3);
diff --git a/src/fortune/edgelist.c b/src/fortune/edgelist.c
new file mode 100644
index 0000000..dfd85a7
--- /dev/null
+++ b/src/fortune/edgelist.c
@@ -0,0 +1,136 @@
+#
+#include "defs.h"
+int nedges;
+struct Freelist efl;
+struct Freelist hfl;
+struct Halfedge *ELleftend, *ELrightend;
+int ELhashsize;
+struct Halfedge **ELhash;
+
+int ntry, totalsearch;
+
+void ELinitialize() {
+	int i;
+	freeinit(&hfl, sizeof **ELhash);
+	ELhashsize = 2 * sqrt_nsites;
+	ELhash = (struct Halfedge **)myalloc(sizeof *ELhash * ELhashsize);
+	for (i = 0; i < ELhashsize; i += 1)
+		ELhash[i] = (struct Halfedge *)NULL;
+	ELleftend = HEcreate((struct Edge *)NULL, 0);
+	ELrightend = HEcreate((struct Edge *)NULL, 0);
+	ELleftend->ELleft = (struct Halfedge *)NULL;
+	ELleftend->ELright = ELrightend;
+	ELrightend->ELleft = ELleftend;
+	ELrightend->ELright = (struct Halfedge *)NULL;
+	ELhash[0] = ELleftend;
+	ELhash[ELhashsize - 1] = ELrightend;
+}
+
+struct Halfedge *HEcreate(e, pm) struct Edge *e;
+int pm;
+{
+	struct Halfedge *answer;
+	answer = (struct Halfedge *)getfree(&hfl);
+	answer->ELedge = e;
+	answer->ELpm = pm;
+	answer->PQnext = (struct Halfedge *)NULL;
+	answer->vertex = (struct Site *)NULL;
+	answer->ELrefcnt = 0;
+	return (answer);
+}
+
+void ELinsert(struct Halfedge *lb, struct Halfedge *new) {
+	new->ELleft = lb;
+	new->ELright = lb->ELright;
+	(lb->ELright)->ELleft = new;
+	lb->ELright = new;
+}
+
+/* Get entry from hash table, pruning any deleted nodes */
+struct Halfedge *ELgethash(b) int b;
+{
+	struct Halfedge *he;
+
+	if (b < 0 || b >= ELhashsize)
+		return ((struct Halfedge *)NULL);
+	he = ELhash[b];
+	if (he == (struct Halfedge *)NULL || he->ELedge != (struct Edge *)DELETED)
+		return (he);
+
+	/* Hash table points to deleted half edge.  Patch as necessary. */
+	ELhash[b] = (struct Halfedge *)NULL;
+	if ((he->ELrefcnt -= 1) == 0)
+		makefree((struct Freenode *)he, &hfl);
+	return ((struct Halfedge *)NULL);
+}
+
+struct Halfedge *ELleftbnd(p) struct Point *p;
+{
+	int i, bucket;
+	struct Halfedge *he;
+
+	/* Use hash table to get close to desired halfedge */
+	bucket = (p->x - xmin) / deltax * ELhashsize;
+	if (bucket < 0)
+		bucket = 0;
+	if (bucket >= ELhashsize)
+		bucket = ELhashsize - 1;
+	he = ELgethash(bucket);
+	if (he == (struct Halfedge *)NULL) {
+		for (i = 1; 1; i += 1) {
+			if ((he = ELgethash(bucket - i)) != (struct Halfedge *)NULL)
+				break;
+			if ((he = ELgethash(bucket + i)) != (struct Halfedge *)NULL)
+				break;
+		};
+		totalsearch += i;
+	};
+	ntry += 1;
+	/* Now search linear list of halfedges for the corect one */
+	if (he == ELleftend || (he != ELrightend && right_of(he, p))) {
+		do {
+			he = he->ELright;
+		} while (he != ELrightend && right_of(he, p));
+		he = he->ELleft;
+	} else
+		do {
+			he = he->ELleft;
+		} while (he != ELleftend && !right_of(he, p));
+
+	/* Update hash table and reference counts */
+	if (bucket > 0 && bucket < ELhashsize - 1) {
+		if (ELhash[bucket] != (struct Halfedge *)NULL)
+			ELhash[bucket]->ELrefcnt -= 1;
+		ELhash[bucket] = he;
+		ELhash[bucket]->ELrefcnt += 1;
+	};
+	return (he);
+}
+
+/* This delete routine can't reclaim node, since pointers from hash
+	 table may be present.   */
+void ELdelete(struct Halfedge *he) {
+	(he->ELleft)->ELright = he->ELright;
+	(he->ELright)->ELleft = he->ELleft;
+	he->ELedge = (struct Edge *)DELETED;
+}
+
+struct Halfedge *ELright(he) struct Halfedge *he;
+{ return (he->ELright); }
+
+struct Halfedge *ELleft(he) struct Halfedge *he;
+{ return (he->ELleft); }
+
+struct Site *leftreg(he) struct Halfedge *he;
+{
+	if (he->ELedge == (struct Edge *)NULL)
+		return (bottomsite);
+	return (he->ELpm == le ? he->ELedge->reg[le] : he->ELedge->reg[re]);
+}
+
+struct Site *rightreg(he) struct Halfedge *he;
+{
+	if (he->ELedge == (struct Edge *)NULL)
+		return (bottomsite);
+	return (he->ELpm == le ? he->ELedge->reg[re] : he->ELedge->reg[le]);
+}
diff --git a/src/fortune/geometry.c b/src/fortune/geometry.c
new file mode 100644
index 0000000..131cacf
--- /dev/null
+++ b/src/fortune/geometry.c
@@ -0,0 +1,184 @@
+#
+#include "defs.h"
+#include <math.h>
+
+void geominit() {
+	struct Edge e;
+	double sn;
+
+	freeinit(&efl, sizeof e);
+	nvertices = 0;
+	nedges = 0;
+	sn = nsites + 4;
+	sqrt_nsites = sqrt(sn);
+	deltay = ymax - ymin;
+	deltax = xmax - xmin;
+}
+
+struct Edge *bisect(s1, s2) struct Site *s1, *s2;
+{
+	double dx, dy, adx, ady;
+	struct Edge *newedge;
+
+	newedge = (struct Edge *)getfree(&efl);
+
+	newedge->reg[0] = s1;
+	newedge->reg[1] = s2;
+	ref(s1);
+	ref(s2);
+	newedge->ep[0] = (struct Site *)NULL;
+	newedge->ep[1] = (struct Site *)NULL;
+
+	dx = s2->coord.x - s1->coord.x;
+	dy = s2->coord.y - s1->coord.y;
+	adx = dx > 0 ? dx : -dx;
+	ady = dy > 0 ? dy : -dy;
+	newedge->c = s1->coord.x * dx + s1->coord.y * dy + (dx * dx + dy * dy) * 0.5;
+	if (adx > ady) {
+		newedge->a = 1.0;
+		newedge->b = dy / dx;
+		newedge->c /= dx;
+	} else {
+		newedge->b = 1.0;
+		newedge->a = dx / dy;
+		newedge->c /= dy;
+	};
+
+	newedge->edgenbr = nedges;
+	out_bisector(newedge);
+	nedges += 1;
+	return (newedge);
+}
+
+struct Site *intersect(el1, el2) struct Halfedge *el1, *el2;
+{
+	struct Edge *e1, *e2, *e;
+	struct Halfedge *el;
+	double d, xint, yint;
+	int right_of_site;
+	struct Site *v;
+
+	e1 = el1->ELedge;
+	e2 = el2->ELedge;
+	if (e1 == (struct Edge *)NULL || e2 == (struct Edge *)NULL)
+		return ((struct Site *)NULL);
+	if (e1->reg[1] == e2->reg[1])
+		return ((struct Site *)NULL);
+
+	d = e1->a * e2->b - e1->b * e2->a;
+	/* printf("intersect: d=%g\n", d); */
+	if (-1.0e-20 < d && d < 1.0e-20) {
+		return ((struct Site *)NULL);
+	};
+
+	xint = (e1->c * e2->b - e2->c * e1->b) / d;
+	yint = (e2->c * e1->a - e1->c * e2->a) / d;
+
+	if ((e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
+			(e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
+			 e1->reg[1]->coord.x < e2->reg[1]->coord.x)) {
+		el = el1;
+		e = e1;
+	} else {
+		el = el2;
+		e = e2;
+	};
+	right_of_site = xint >= e->reg[1]->coord.x;
+	if ((right_of_site && el->ELpm == le) || (!right_of_site && el->ELpm == re))
+		return ((struct Site *)NULL);
+
+	v = (struct Site *)getfree(&sfl);
+	v->refcnt = 0;
+	v->coord.x = xint;
+	v->coord.y = yint;
+	return (v);
+}
+
+/* returns 1 if p is to right of halfedge e */
+int right_of(el, p) struct Halfedge *el;
+struct Point *p;
+{
+	struct Edge *e;
+	struct Site *topsite;
+	int right_of_site, above, fast;
+	double dxp, dyp, dxs, t1, t2, t3, yl;
+
+	e = el->ELedge;
+	topsite = e->reg[1];
+	right_of_site = p->x > topsite->coord.x;
+	if (right_of_site && el->ELpm == le)
+		return (1);
+	if (!right_of_site && el->ELpm == re)
+		return (0);
+
+	if (e->a == 1.0) {
+		dyp = p->y - topsite->coord.y;
+		dxp = p->x - topsite->coord.x;
+		fast = 0;
+		if ((!right_of_site && e->b < 0.0) | (right_of_site && e->b >= 0.0)) {
+			above = dyp >= e->b * dxp;
+			fast = above;
+		} else {
+			above = p->x + p->y * e->b > e->c;
+			if (e->b < 0.0)
+				above = !above;
+			if (!above)
+				fast = 1;
+		};
+		if (!fast) {
+			dxs = topsite->coord.x - (e->reg[0])->coord.x;
+			above = e->b * (dxp * dxp - dyp * dyp) <
+							dxs * dyp * (1.0 + 2.0 * dxp / dxs + e->b * e->b);
+			if (e->b < 0.0)
+				above = !above;
+		};
+	} else /*e->b==1.0 */
+			{
+		yl = e->c - e->a * p->x;
+		t1 = p->y - yl;
+		t2 = p->x - topsite->coord.x;
+		t3 = yl - topsite->coord.y;
+		above = t1 * t1 > t2 * t2 + t3 * t3;
+	};
+	return (el->ELpm == le ? above : !above);
+}
+
+void endpoint(e, lr, s) struct Edge *e;
+int lr;
+struct Site *s;
+{
+	e->ep[lr] = s;
+	ref(s);
+	if (e->ep[re - lr] == (struct Site *)NULL) {
+	} else {
+		out_ep(e);
+		deref(e->reg[le]);
+		deref(e->reg[re]);
+		makefree((struct Freenode *)e, &efl);
+	}
+}
+
+double dist(s, t) struct Site *s, *t;
+{
+	double dx, dy;
+	dx = s->coord.x - t->coord.x;
+	dy = s->coord.y - t->coord.y;
+	return (sqrt(dx * dx + dy * dy));
+}
+
+void makevertex(v) struct Site *v;
+{
+	v->sitenbr = nvertices;
+	nvertices += 1;
+	out_vertex(v);
+}
+
+void deref(v) struct Site *v;
+{
+	v->refcnt -= 1;
+	if (v->refcnt == 0)
+		makefree((struct Freenode *)v, &sfl);
+}
+
+void ref(v) struct Site *v;
+{ v->refcnt += 1; }
diff --git a/src/fortune/heap.c b/src/fortune/heap.c
new file mode 100644
index 0000000..2ebbd58
--- /dev/null
+++ b/src/fortune/heap.c
@@ -0,0 +1,94 @@
+#
+#include "defs.h"
+int PQhashsize;
+struct Halfedge *PQhash;
+struct Halfedge *PQfind();
+int PQcount;
+int PQmin;
+int PQempty();
+
+void PQinsert(he, v, offset) struct Halfedge *he;
+struct Site *v;
+double offset;
+{
+	struct Halfedge *last, *next;
+	he->vertex = v;
+	ref(v);
+	he->ystar = v->coord.y + offset;
+	last = &PQhash[PQbucket(he)];
+	while ((next = last->PQnext) != (struct Halfedge *)NULL &&
+				 (he->ystar > next->ystar ||
+					(he->ystar == next->ystar && v->coord.x > next->vertex->coord.x))) {
+		last = next;
+	};
+	he->PQnext = last->PQnext;
+	last->PQnext = he;
+	PQcount += 1;
+}
+
+void PQdelete(he) struct Halfedge *he;
+{
+	struct Halfedge *last;
+
+	if (he->vertex != (struct Site *)NULL) {
+		last = &PQhash[PQbucket(he)];
+		while (last->PQnext != he)
+			last = last->PQnext;
+		last->PQnext = he->PQnext;
+		PQcount -= 1;
+		deref(he->vertex);
+		he->vertex = (struct Site *)NULL;
+	};
+}
+
+int PQbucket(he) struct Halfedge *he;
+{
+	int bucket;
+
+	if (he->ystar < ymin)
+		bucket = 0;
+	else if (he->ystar >= ymax)
+		bucket = PQhashsize - 1;
+	else
+		bucket = (he->ystar - ymin) / deltay * PQhashsize;
+	if (bucket < 0)
+		bucket = 0;
+	if (bucket >= PQhashsize)
+		bucket = PQhashsize - 1;
+	if (bucket < PQmin)
+		PQmin = bucket;
+	return (bucket);
+}
+
+int PQempty() { return (PQcount == 0); }
+
+struct Point PQ_min() {
+	struct Point answer;
+
+	while (PQhash[PQmin].PQnext == (struct Halfedge *)NULL) {
+		PQmin += 1;
+	};
+	answer.x = PQhash[PQmin].PQnext->vertex->coord.x;
+	answer.y = PQhash[PQmin].PQnext->ystar;
+	return (answer);
+}
+
+struct Halfedge *PQextractmin() {
+	struct Halfedge *curr;
+	curr = PQhash[PQmin].PQnext;
+	PQhash[PQmin].PQnext = curr->PQnext;
+	PQcount -= 1;
+	return (curr);
+}
+
+void PQinitialize() {
+	int i;
+//	struct Point *s;
+
+	PQcount = 0;
+	PQmin = 0;
+	PQhashsize = 4 * sqrt_nsites;
+	PQhash = (struct Halfedge *)myalloc(PQhashsize * sizeof *PQhash);
+	for (i = 0; i < PQhashsize; i += 1)
+		PQhash[i].PQnext = (struct Halfedge *)NULL;
+}
diff --git a/src/fortune/main.c b/src/fortune/main.c
new file mode 100644
index 0000000..7571945
--- /dev/null
+++ b/src/fortune/main.c
@@ -0,0 +1,284 @@
+#
+#include <stdio.h>
+#include <stdint.h>
+#include "defs.h"
+#include <math.h>
+struct Site *nextone();
+
+int triangulate, sorted, plot, debug;
+double xmin, xmax, ymin, ymax, deltax, deltay;
+int vert_count, edge_count, line_count, dual_count, site_count;
+double *site_list, *vert_list, *line_list;
+unsigned int *edge_list, *dual_list;
+
+/* sort sites on y, then x, coord */
+int scomp(s1, s2) struct Point *s1, *s2;
+{
+	if (s1->y < s2->y)
+		return (-1);
+	if (s1->y > s2->y)
+		return (1);
+	if (s1->x < s2->x)
+		return (-1);
+	if (s1->x > s2->x)
+		return (1);
+	return (0);
+}
+
+/* read all sites, sort, and compute xmin, xmax, ymin, ymax */
+void readsites(double *lattice) {
+
+	sites = (struct Site *)myalloc(nsites * sizeof *sites);
+	for (unsigned int j = 0; j < nsites; j++) {
+		sites[j].coord.x = lattice[2 * j];
+		sites[j].coord.y = lattice[2 * j + 1];
+		sites[j].sitenbr = j;
+		sites[j].refcnt = 0;
+	}
+	qsort(sites, nsites, sizeof(struct Site), scomp);
+	xmin = sites[0].coord.x;
+	xmax = sites[0].coord.x;
+	for (unsigned int i = 1; i < nsites; i += 1) {
+		if (sites[i].coord.x < xmin)
+			xmin = sites[i].coord.x;
+		if (sites[i].coord.x > xmax)
+			xmax = sites[i].coord.x;
+	};
+	ymin = sites[0].coord.y;
+	ymax = sites[nsites - 1].coord.y;
+}
+
+unsigned int delete_duplicates(unsigned int ne, unsigned int *etv) {
+	unsigned int ndup = 0;
+	bool *duplicates = (bool *)calloc(ne, sizeof(bool));
+	for (unsigned int i = 0; i < ne; i++) {
+		unsigned int v1 = etv[2 * i];
+		unsigned int v2 = etv[2 * i + 1];
+		for (unsigned int j = 0; j < i; j++) {
+			unsigned int w1 = etv[2 * j];
+			unsigned int w2 = etv[2 * j + 1];
+			if (w1 == v1 && w2 == v2) {
+				duplicates[j] = true;
+				ndup++;
+				break;
+			}
+		}
+	}
+	unsigned int newsize = (int)ne - (int)ndup;
+	unsigned int count = 0;
+	for (unsigned int i = 0; i < ne; i++) {
+		if (!duplicates[i]) {
+			etv[2 * count] = etv[2 * i];
+			etv[2 * count + 1] = etv[2 * i + 1];
+			count++;
+		}
+	}
+	free(duplicates);
+	return newsize;
+}
+
+intptr_t *run_voronoi(unsigned int num, double *lattice, double xmin, double xmax, double ymin, double ymax)
+{
+	struct Site *(*next)();
+
+	sorted = 0;
+	triangulate = 0;
+	plot = 0;
+	debug = 0;
+
+	alloclist = (char **)malloc(9 * num * sizeof(char *));
+	allocnum = 0;
+
+	freeinit(&sfl, sizeof *sites);
+
+	unsigned int eff_num = 9 * num;
+	double *eff_lattice = (double *)malloc(2 * eff_num * sizeof(double));
+
+	for (unsigned int i = 0; i < num; i++) {
+		// original sites - our baby boys
+		eff_lattice[2*i] = lattice[2*i];
+		eff_lattice[2*i+1] = lattice[2*i+1];
+
+		// sites to the right
+		eff_lattice[2*(num+i)+1] = lattice[2*i+1] - xmin + xmax;
+		eff_lattice[2*(num+i)] = lattice[2*i];
+
+		// sites to the left
+		eff_lattice[2*(2*num+i)+1] = lattice[2*i+1] - xmax + xmin;
+		eff_lattice[2*(2*num+i)] = lattice[2*i];
+
+		// sites to the top
+		eff_lattice[2*(3*num+i)+1] = lattice[2*i+1];
+		eff_lattice[2*(3*num+i)] = lattice[2*i] - ymin + ymax;
+
+		// sites to the bottom
+		eff_lattice[2*(4*num+i)+1] = lattice[2*i+1];
+		eff_lattice[2*(4*num+i)] = lattice[2*i] - ymax + ymin;
+
+		// sites to the upper right
+		eff_lattice[2*(5*num+i)+1] = lattice[2*i+1] - xmin + xmax;
+		eff_lattice[2*(5*num+i)] = lattice[2*i] - ymin + ymax;
+
+		// sites to the upper left
+		eff_lattice[2*(6*num+i)+1] = lattice[2*i+1] - xmax + xmin;
+		eff_lattice[2*(6*num+i)] = lattice[2*i] - ymin + ymax;
+
+		// sites to the lower left
+		eff_lattice[2*(7*num+i)+1] = lattice[2*i+1] - xmax + xmin;
+		eff_lattice[2*(7*num+i)] = lattice[2*i] + ymin - ymax;
+
+		// sites to the lower right
+		eff_lattice[2*(8*num+i)+1] = lattice[2*i+1] + xmax - xmin;
+		eff_lattice[2*(8*num+i)] = lattice[2*i] + ymin - ymax;
+	}
+
+	nsites = eff_num;
+	readsites(eff_lattice);
+	free(eff_lattice);
+	next = nextone;
+
+	siteidx = 0;
+	geominit();
+
+	site_list = malloc(2 * nsites * sizeof(double));
+	vert_list = NULL;
+	line_list = NULL;
+	edge_list = NULL;
+	dual_list = NULL;
+	site_count = 0;
+	vert_count = 0;
+	edge_count = 0;
+	line_count = 0;
+	dual_count = 0;
+
+	voronoi(triangulate, next);
+
+	unsigned int real_vert_count = vert_count;
+	unsigned int real_edge_count = edge_count;
+	unsigned int real_dual_count = dual_count;
+	double *real_vert_list = vert_list;
+	unsigned int *real_edge_list = edge_list;
+	unsigned int *real_dual_list = dual_list;
+
+	{
+		real_vert_count = 0;
+		real_edge_count = 0;
+		real_dual_count = 0;
+		real_vert_list = (double *)malloc(2 * vert_count * sizeof(double));
+		real_edge_list = (unsigned int *)malloc(2 * edge_count * sizeof(unsigned int));
+		real_dual_list = (unsigned int *)malloc(3 * dual_count * sizeof(unsigned int));
+		unsigned int *triangle_analyzed = (unsigned int *)malloc(4 * dual_count * sizeof(unsigned int));
+		unsigned int q = 0;
+		int *new_vert_ind = (int *)malloc(vert_count * sizeof(int));
+
+		for (unsigned int i = 0; i < vert_count; i++) {
+			unsigned int t1, t2, t3;
+			t1 = dual_list[3*i]; t2 = dual_list[3*i+1]; t3 = dual_list[3*i+2];
+			if (t1 >= num && t2 >= num && t3 >= num) {
+				new_vert_ind[i] = -1;
+			}
+			else if (t1 < num && t2 < num && t3 < num) {
+				real_vert_list[2*real_vert_count] = vert_list[2*i];
+				real_vert_list[2*real_vert_count+1] = vert_list[2*i+1];
+				real_dual_list[3*real_vert_count] = t1;
+				real_dual_list[3*real_vert_count+1] = t2;
+				real_dual_list[3*real_vert_count+2] = t3;
+				new_vert_ind[i] = real_vert_count;
+				real_vert_count++;
+			}
+			else {
+				unsigned int tt1, tt2, tt3;
+				tt1 = t1 % num; tt2 = t2 % num; tt3 = t3 % num;
+				unsigned int s1, s2, s3;
+				s1 = tt1 < tt2 ? (tt1 < tt3 ? tt1 : tt3) : (tt2 < tt3 ? tt2 : tt3);
+				s3 = tt1 > tt2 ? (tt1 > tt3 ? tt1 : tt3) : (tt2 > tt3 ? tt2 : tt3);
+				s2 = tt1 < tt2 ? (tt1 > tt3 ? tt1 : (tt2 < tt3 ? tt2 : tt3)) : (tt1 < tt3 ? tt1 : (tt2 < tt3 ? tt3: tt2));
+
+				bool matched = false;
+				unsigned int ass_vert;
+
+				for (unsigned int j = 0; j < q; j++) {
+					if (s1 == triangle_analyzed[4*j] && s2 == triangle_analyzed[4*j+1] && s3 == triangle_analyzed[4*j+2]) {
+						matched = true;
+						ass_vert = triangle_analyzed[4*j+3];
+						break;
+					}
+				}
+
+				if (matched) {
+					new_vert_ind[i] = ass_vert;
+				} else {
+					triangle_analyzed[4*q] = s1;
+					triangle_analyzed[4*q+1] = s2;
+					triangle_analyzed[4*q+2] = s3;
+					triangle_analyzed[4*q+3] = real_vert_count;
+					q++;
+					real_vert_list[2*real_vert_count+1] = fmod(2*xmax+vert_list[2*i+1], xmax);
+					real_vert_list[2*real_vert_count] = fmod(2*ymax+vert_list[2*i], ymax);
+					real_dual_list[3*real_vert_count] = t1 % num;
+					real_dual_list[3*real_vert_count+1] = t2 % num;
+					real_dual_list[3*real_vert_count+2] = t3 % num;
+					new_vert_ind[i] = real_vert_count;
+					real_vert_count++;
+				}
+			}
+		}
+		for (unsigned int i = 0; i < edge_count; i++) {
+			unsigned int v1, v2;
+			v1 = edge_list[2 * i];
+			v2 = edge_list[2 * i+1];
+			if (v1 != UINT_MAX && v2 != UINT_MAX) {
+				if (new_vert_ind[v1] >= 0 && new_vert_ind[v2] >=0) {
+					unsigned int nv1 = new_vert_ind[v1];
+					unsigned int nv2 = new_vert_ind[v2];
+					real_edge_list[2*real_edge_count] = nv1 < nv2 ? nv1 : nv2;
+					real_edge_list[2*real_edge_count+1] = nv1 < nv2 ? nv2 : nv1;
+					real_edge_count++;
+				}
+			}
+		}
+		unsigned int new_edge_count = delete_duplicates(real_edge_count, real_edge_list);
+		real_edge_count = new_edge_count;
+		free(triangle_analyzed);
+		free(new_vert_ind);
+		free(dual_list);
+		free(edge_list);
+		free(vert_list);
+	}
+
+
+	intptr_t *output = (intptr_t *)malloc(6 * sizeof(intptr_t));
+	output[0] = (intptr_t)malloc(4 * sizeof(unsigned int));
+	((unsigned int *)output[0])[0] = real_vert_count;
+	((unsigned int *)output[0])[1] = real_edge_count;
+	((unsigned int *)output[0])[2] = line_count;
+	((unsigned int *)output[0])[3] = real_vert_count;
+	output[1] = (intptr_t)site_list;
+	output[2] = (intptr_t)real_vert_list;
+	output[3] = (intptr_t)real_edge_list;
+	output[4] = (intptr_t)line_list;
+	output[5] = (intptr_t)real_dual_list;
+
+	free(ELhash);
+	free(PQhash);
+	free(sites);
+	for (unsigned int i = 0; i < allocnum; i++) {
+		free(alloclist[i]);
+	}
+	free(alloclist);
+
+	return output;
+}
+
+/* return a single in-storage site */
+struct Site *nextone() {
+	for (; siteidx < nsites; siteidx += 1) {
+		if (siteidx == 0 || sites[siteidx].coord.x != sites[siteidx - 1].coord.x ||
+				sites[siteidx].coord.y != sites[siteidx - 1].coord.y) {
+			siteidx += 1;
+			return (&sites[siteidx - 1]);
+		};
+	};
+	return ((struct Site *)NULL);
+}
+
diff --git a/src/fortune/memory.c b/src/fortune/memory.c
new file mode 100644
index 0000000..807d1b0
--- /dev/null
+++ b/src/fortune/memory.c
@@ -0,0 +1,44 @@
+#
+#include "defs.h"
+#include <stdio.h>
+
+void freeinit(struct Freelist *fl, int size) {
+	fl->head = (struct Freenode *)NULL;
+	fl->nodesize = size;
+}
+
+char *getfree(fl) struct Freelist *fl;
+{
+	int i;
+	struct Freenode *t;
+	if (fl->head == (struct Freenode *)NULL) {
+		t = (struct Freenode *)myalloc(sqrt_nsites * fl->nodesize);
+		alloclist[allocnum] = t;
+		allocnum++;
+		for (i = 0; i < sqrt_nsites; i += 1)
+			makefree((struct Freenode *)((char *)t + i * fl->nodesize), fl);
+	};
+	t = fl->head;
+	fl->head = (fl->head)->nextfree;
+	return ((char *)t);
+}
+
+void makefree(curr, fl) struct Freenode *curr;struct Freelist *fl;
+{
+	curr->nextfree = fl->head;
+	fl->head = curr;
+}
+
+int total_alloc;
+char *myalloc(n) unsigned n;
+{
+	char *t;
+	if ((t = malloc(n)) == (char *)0) {
+		fprintf(stderr,
+						"Insufficient memory processing site %d (%d bytes in use)\n",
+						siteidx, total_alloc);
+		exit(1);
+	};
+	total_alloc += n;
+	return (t);
+}
diff --git a/src/fortune/output.c b/src/fortune/output.c
new file mode 100644
index 0000000..d496feb
--- /dev/null
+++ b/src/fortune/output.c
@@ -0,0 +1,46 @@
+#
+#include "defs.h"
+#include <stdio.h>
+double pxmin, pxmax, pymin, pymax, cradius;
+
+void out_bisector(e) struct Edge *e;
+{
+	if (line_count % nsites == 0) line_list = realloc(line_list, 3 * (nsites + line_count) * sizeof(double));
+	line_list[3*line_count] = e->a;
+	line_list[3*line_count+1] = e->b;
+	line_list[3*line_count+2] = e->c;
+	line_count++;
+}
+
+void out_ep(e) struct Edge *e;
+{
+	if (edge_count % nsites == 0) edge_list = realloc(edge_list, 2 * (nsites + edge_count) * sizeof(unsigned int));
+	edge_list[2 * edge_count] = (e->ep[le] != (struct Site *)NULL ? e->ep[le]->sitenbr : UINT_MAX);
+	edge_list[2 * edge_count + 1] = (e->ep[re] != (struct Site *)NULL ? e->ep[re]->sitenbr : UINT_MAX);
+	edge_count++;
+}
+
+void out_vertex(v) struct Site *v;
+{
+	if (vert_count % nsites == 0) vert_list = realloc(vert_list, 2 * (nsites + vert_count) * sizeof(double));
+	vert_list[2 * vert_count] = v->coord.x;
+	vert_list[2 * vert_count + 1] = v->coord.y;
+	vert_count++;
+}
+
+void out_site(s) struct Site *s;
+{
+	site_list[2 * site_count] = s->coord.x;
+	site_list[2 * site_count + 1] = s->coord.y;
+	site_count++;
+}
+
+void out_triple(s1, s2, s3) struct Site *s1, *s2, *s3;
+{
+	if (dual_count % nsites == 0) dual_list = realloc(dual_list, 3 * (nsites + dual_count) * sizeof(unsigned int));
+	dual_list[dual_count * 3] = s1->sitenbr;
+	dual_list[dual_count * 3 + 1] = s2->sitenbr;
+	dual_list[dual_count * 3 + 2] = s3->sitenbr;
+	dual_count++;
+}
+
diff --git a/src/fortune/voronoi.c b/src/fortune/voronoi.c
new file mode 100644
index 0000000..dc9945f
--- /dev/null
+++ b/src/fortune/voronoi.c
@@ -0,0 +1,103 @@
+#
+#include "defs.h"
+
+struct Site *sites;
+int nsites;
+int siteidx;
+int sqrt_nsites;
+int nvertices;
+struct Freelist sfl;
+struct Site *bottomsite;
+
+/* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
+	 deltax, deltay (can all be estimates).
+	 Performance suffers if they are wrong; better to make nsites,
+	 deltax, and deltay too big than too small.  (?) */
+
+void voronoi(triangulate, nextsite) int triangulate;
+struct Site *(*nextsite)();
+{
+	struct Site *newsite, *bot, *top, *temp, *p;
+	struct Site *v;
+	struct Point newintstar;
+	int pm;
+	struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
+	struct Edge *e;
+
+	PQinitialize();
+	bottomsite = (*nextsite)();
+	out_site(bottomsite);
+	ELinitialize();
+
+	newsite = (*nextsite)();
+	while (1) {
+		if (!PQempty())
+			newintstar = PQ_min();
+
+		if (newsite != (struct Site *)NULL &&
+				(PQempty() || newsite->coord.y < newintstar.y ||
+				 (newsite->coord.y == newintstar.y &&
+					newsite->coord.x < newintstar.x))) {/* new site is smallest */
+			out_site(newsite);
+			lbnd = ELleftbnd(&(newsite->coord));
+			rbnd = ELright(lbnd);
+			bot = rightreg(lbnd);
+			e = bisect(bot, newsite);
+			bisector = HEcreate(e, le);
+			ELinsert(lbnd, bisector);
+			if ((p = intersect(lbnd, bisector)) != (struct Site *)NULL) {
+				PQdelete(lbnd);
+				PQinsert(lbnd, p, dist(p, newsite));
+			};
+			lbnd = bisector;
+			bisector = HEcreate(e, re);
+			ELinsert(lbnd, bisector);
+			if ((p = intersect(bisector, rbnd)) != (struct Site *)NULL) {
+				PQinsert(bisector, p, dist(p, newsite));
+			};
+			newsite = (*nextsite)();
+		} else if (!PQempty())
+				/* intersection is smallest */
+				{
+			lbnd = PQextractmin();
+			llbnd = ELleft(lbnd);
+			rbnd = ELright(lbnd);
+			rrbnd = ELright(rbnd);
+			bot = leftreg(lbnd);
+			top = rightreg(rbnd);
+			out_triple(bot, top, rightreg(lbnd));
+			v = lbnd->vertex;
+			makevertex(v);
+			endpoint(lbnd->ELedge, lbnd->ELpm, v);
+			endpoint(rbnd->ELedge, rbnd->ELpm, v);
+			ELdelete(lbnd);
+			PQdelete(rbnd);
+			ELdelete(rbnd);
+			pm = le;
+			if (bot->coord.y > top->coord.y) {
+				temp = bot;
+				bot = top;
+				top = temp;
+				pm = re;
+			}
+			e = bisect(bot, top);
+			bisector = HEcreate(e, pm);
+			ELinsert(llbnd, bisector);
+			endpoint(e, re - pm, v);
+			deref(v);
+			if ((p = intersect(llbnd, bisector)) != (struct Site *)NULL) {
+				PQdelete(llbnd);
+				PQinsert(llbnd, p, dist(p, bot));
+			};
+			if ((p = intersect(bisector, rrbnd)) != (struct Site *)NULL) {
+				PQinsert(bisector, p, dist(p, bot));
+			};
+		} else
+			break;
+	};
+
+	for (lbnd = ELright(ELleftend); lbnd != ELrightend; lbnd = ELright(lbnd)) {
+		e = lbnd->ELedge;
+		out_ep(e);
+	};
+}
diff --git a/src/fracture.c b/src/fracture.c
new file mode 100644
index 0000000..6798231
--- /dev/null
+++ b/src/fracture.c
@@ -0,0 +1,578 @@
+
+
+#include "fracture.h"
+
+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;
+
+	filename_len = 100;
+
+	num = 100;
+	width = 16;
+	crack_len = 16;
+	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;
+	save_damage = false;
+	save_homo_damage = false;
+	save_corr = false;
+	save_cond = false;
+	use_crit = true;
+	use_first = false;
+
+	int periodic_int;
+
+	while ((opt = getopt(argc, argv, "n:w:b:l:f:ocp:vVsrNCaSFtzdH")) != -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);
+		}
+	}
+
+	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);
+
+	/* 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;
+	}
+
+	fnet *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, crack_width, &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, crack_width, &c);
+		finish_instance(perm_instance, &c);
+		c_dist_size = network->num_dual_verts;
+		a_dist_size = network->num_verts;
+	}
+
+	// define arrays for saving cluster and avalanche distributions
+	unsigned int *cluster_size_dist;
+	unsigned int *avalanche_size_dist;
+	char *c_filename;
+	if (save_clusters) {
+		cluster_size_dist =
+				(unsigned int *)calloc(c_dist_size, sizeof(unsigned int));
+		avalanche_size_dist =
+				(unsigned int *)calloc(a_dist_size, sizeof(unsigned int));
+
+		c_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(c_filename, filename_len, "cluster_%d_%g_%g.txt", width, crack_len,
+						 beta);
+
+		FILE *cluster_out = fopen(c_filename, "r");
+
+		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;
+			}
+			fclose(cluster_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 *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 *conductivity;
+	if (save_cond) {
+		conductivity = (double *)malloc(num * 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 *avg_ztress;
+	if (save_avg_ztress) {
+		avg_ztress = (double *)calloc(pow(width, 2), sizeof(double));
+	}
+
+	double *damage;
+	if (save_damage) {
+		damage = (double *)calloc(pow(width, 2), sizeof(double));
+	}
+
+	double *homo_damage;
+	if (save_homo_damage) {
+		homo_damage = (double *)calloc(num, 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;
+		}
+	}
+
+	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);
+
+		break_data *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, crack_width, &c);
+				finish_instance(perm_instance, &c);
+				if (supplied_bound) {
+					voronoi_bound_ini(perm_instance, square_bound, width);
+				}
+			}
+			double *fuse_thres = gen_fuse_thres(
+					network->num_edges, network->edge_coords, 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);
+		}
+
+		unsigned int min_pos = 0;
+		double min_val = DBL_MAX;
+
+		for (unsigned int j = 0; j < breaking_data->num_broken; j++) {
+			double val = fabs(breaking_data->extern_field[j]);
+
+
+			if (val < min_val) {
+				min_pos = j;
+				min_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]);
+		}
+
+		finst *tmp_instance = copy_instance(perm_instance, &c);
+
+		int stop_at = breaking_data->num_broken;
+		if (use_crit)
+			stop_at = min_pos;
+		if (use_first)
+			stop_at = 0;
+
+		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) {
+					av_size++;
+				}
+
+				if (val < cur_val) {
+					avalanche_size_dist[av_size]++;
+					av_size = 0;
+					cur_val = val;
+				}
+			}
+		}
+
+		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->network->num_verts + 1] - tmp_voltage[tmp_instance->network->num_verts]);
+			free(tmp_voltage);
+		}
+
+		if (save_homo_damage) {
+			homo_damage[DUMB] = ((double)min_pos) / tmp_instance->network->num_edges;
+		}
+
+		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->network, 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 (sigma == sigma && save_avg_ztress) {
+					avg_ztress[i] += sigma;
+				}
+			}
+			free(tmp_stress);
+		}
+
+		if (save_damage) {
+			double *tmp_damage = (double *)calloc(tmp_instance->network->num_edges, 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->network, 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];
+			}
+			free(tmp_damage);
+		}
+
+		if (save_clusters) {
+			unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c);
+			for (unsigned int i = 0; i < network->num_dual_verts; i++) {
+				cluster_size_dist[i] += tmp_cluster_dist[i];
+			}
+			free(tmp_cluster_dist);
+		}
+
+		if (save_corr) {
+			double *tmp_corr = get_corr(tmp_instance, dists, &c);
+			for (unsigned int i = 0; i < tmp_instance->network->num_dual_verts; i++) {
+				avg_corr[i] += tmp_corr[i] / num;
+			}
+			free(tmp_corr);
+		}
+
+		if (network_out) {
+			FILE *net_out = fopen("network.txt", "w");
+			for (unsigned int i = 0; i < network->num_verts; i++) {
+				fprintf(net_out, "%f %f ", network->vert_coords[2 * i],
+								network->vert_coords[2 * i + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int i = 0; i < network->num_edges; i++) {
+				fprintf(net_out, "%u %u ", network->edges_to_verts[2 * i],
+								network->edges_to_verts[2 * i + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int i = 0; i < network->num_dual_verts; i++) {
+				fprintf(net_out, "%f %f ", network->dual_vert_coords[2 * i],
+								network->dual_vert_coords[2 * i + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int i = 0; i < network->num_edges; i++) {
+				fprintf(net_out, "%u %u ", network->dual_edges_to_verts[2 * i],
+								network->dual_edges_to_verts[2 * i + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int i = 0; i < network->num_edges; i++) {
+				fprintf(net_out, "%d ", tmp_instance->fuses[i]);
+			}
+			fclose(net_out);
+		}
+
+		free_instance(tmp_instance, &c);
+		if (voronoi && !repeat_voronoi) {
+			free_fnet(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]);
+			}
+			fprintf(break_out, "\n");
+		}
+		free_break_data(breaking_data);
+	}
+
+	printf("\033[F\033[JFRACTURE: COMPLETE\n");
+
+	if (save_clusters) {
+		FILE *cluster_out = fopen(c_filename, "w");
+
+		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);
+
+		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_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_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 (!voronoi || repeat_voronoi) {
+		free_instance(perm_instance, &c);
+		//free_fnet(network, &c);
+	}
+
+	if (include_breaking) {
+		fclose(break_out);
+	}
+
+	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_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_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);
+
+	}
+
+	CHOL_F(finish)(&c);
+
+	return 0;
+}
diff --git a/src/fracture.h b/src/fracture.h
new file mode 100644
index 0000000..4ee7050
--- /dev/null
+++ b/src/fracture.h
@@ -0,0 +1,182 @@
+
+#pragma once
+
+#include <assert.h>
+#include <cholmod.h>
+#include <float.h>
+#include <getopt.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_sf_erf.h>
+#include <gsl/gsl_sf_exp.h>
+#include <gsl/gsl_sf_log.h>
+#include <inttypes.h>
+#include <math.h>
+#include <omp.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <unistd.h>
+
+// these two defs allow me to switch to long int cholmod in a sitch
+#define CHOL_INT long int
+#define CHOL_F(x) cholmod_l_##x
+
+typedef enum boundary_type {
+	FREE_BOUND,
+	CYLINDER_BOUND,
+	TORUS_BOUND
+} boundary_type;
+
+typedef struct {
+	unsigned int num_edges;
+	unsigned int num_verts;
+	unsigned int num_verts_break;
+	unsigned int num_bounds;
+	unsigned int *edges_to_verts;
+	unsigned int *edges_to_verts_break;
+	unsigned int *verts_to_edges_ind;
+	unsigned int *verts_to_edges;
+	unsigned int *bound_inds;
+	unsigned int *bound_verts;
+	double *vert_coords;
+	double *edge_coords;
+	unsigned int *dual_edges_to_verts;
+	unsigned int *dual_verts_to_edges_ind;
+	unsigned int *dual_verts_to_edges;
+	double *dual_vert_coords;
+	unsigned int num_spanning_edges;
+	unsigned int *spanning_edges;
+	double L;
+	unsigned int num_dual_verts;
+	unsigned int break_dim;
+	cholmod_sparse *voltcurmat;
+	boundary_type boundary;
+} fnet;
+
+typedef struct {
+	fnet *network;
+	unsigned int num_remaining_edges;
+	bool *fuses;
+	double inf;
+	cholmod_dense *boundary_cond;
+	cholmod_factor *factor;
+	unsigned int *marks;
+	unsigned int *dual_marks;
+	bool voltage_bound;
+	unsigned int num_components;
+	cholmod_sparse *adjacency;
+	bool debug_stop;
+} finst;
+
+typedef struct {
+	unsigned int num_broken;
+	unsigned int *break_list;
+	double *conductivity;
+	double *extern_field;
+} break_data;
+
+intptr_t *run_voronoi(unsigned int num_coords, double *coords, double xmin, double xmax, double ymin, double ymax);
+
+int update_components(const cholmod_sparse *laplacian, unsigned int *marks,
+											int old_num_components, int v1, int v2, int exclude);
+
+unsigned int *find_components(const cholmod_sparse *laplacian, unsigned int skip);
+
+cholmod_sparse *gen_adjacency(const finst *instance, bool dual, bool breakv,
+															unsigned int pad, cholmod_common *c);
+
+cholmod_sparse *gen_laplacian(const finst *instance, cholmod_common *c,
+															bool symmetric);
+
+int edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
+									bool index);
+
+int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
+											 bool index);
+
+double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert,
+													bool index);
+
+bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2,
+									 cholmod_common *c);
+
+break_data *fracture_network(finst *instance, double *fuse_thres,
+														 cholmod_common *c, double cutoff);
+
+double *get_current(const finst *instance, cholmod_common *c);
+double *get_current_v(const finst *instance, double *voltages, cholmod_common *c);
+double *get_voltage(const finst *instance, cholmod_common *c);
+
+double *gen_fuse_thres(unsigned int num_edges, double *edge_coords, double beta,
+											 double (*beta_scale)(double, double, double));
+
+bool gen_crack(finst *instance, double crack_len, double crack_width,
+							 cholmod_common *c);
+
+void update_boundary(finst *instance, const double *avg_field);
+
+FILE *get_file(const char *prefix, unsigned int width, unsigned int crack,
+							 double beta, unsigned int iter, unsigned int num_iter,
+							 unsigned int num, bool read);
+
+double update_beta(double beta, unsigned int width, const double *stress,
+									 const double *damage, double bound_total);
+
+cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts,
+															 unsigned int *edges_to_verts, cholmod_common *c);
+
+finst *copy_instance(const finst *instance, cholmod_common *c);
+
+fnet *ini_square_network(unsigned int width, boundary_type boundary, bool side_bounds,
+												 cholmod_common *c);
+
+void free_fnet(fnet *network, cholmod_common *c);
+void free_instance(finst *instance, cholmod_common *c);
+
+finst *create_instance(fnet *network, double inf, bool voltage_bound,
+											 bool startnow, cholmod_common *c);
+
+fnet *ini_voronoi_network(unsigned int L, boundary_type boundary,
+													double *(*genfunc)(unsigned int, gsl_rng *, unsigned int *),
+													cholmod_common *c);
+
+bool check_instance(const finst *instance, cholmod_common *c);
+
+bool break_edge(finst *instance, unsigned int edge, cholmod_common *c);
+
+void finish_instance(finst *instance, cholmod_common *c);
+
+finst *coursegrain_square(finst *instance, fnet *network_p, cholmod_common *c);
+
+unsigned int *get_clusters(finst *instance, cholmod_common *c);
+
+unsigned int *get_cluster_dist(finst *instance, cholmod_common *c);
+
+double *genfunc_uniform(unsigned int L, gsl_rng *r, unsigned int *num);
+double *genfunc_hyperuniform(unsigned int L, gsl_rng *r, unsigned int *num);
+void randfunc_flat(gsl_rng *r, double *x, double *y);
+void randfunc_gaus(gsl_rng *r, double *x, double *y);
+double beta_scaling_flat(double beta, double x, double y);
+double beta_scaling_gaus(double beta, double x, double y);
+double beta_mag(double beta);
+
+unsigned int *dijkstra(fnet *network, unsigned int source);
+unsigned int **get_dists(fnet *network);
+double *get_corr(finst *instance, unsigned int **dists, cholmod_common *c);
+
+double *bin_values(fnet *network, unsigned int width, double *values);
+
+void voronoi_bound_ini(finst *instance, double *square_bound,
+											 unsigned int width);
+
+break_data *alloc_break_data(unsigned int num_edges);
+void free_break_data(break_data *data);
+void update_break_data(break_data *data, unsigned int last_broke, double strength, double conductivity);
+
+double get_conductivity(finst *inst, double *current, cholmod_common *c);
diff --git a/src/fracture_network.c b/src/fracture_network.c
new file mode 100644
index 0000000..148b08c
--- /dev/null
+++ b/src/fracture_network.c
@@ -0,0 +1,61 @@
+
+#include "fracture.h"
+
+int inc_break_fuses(finst *instance, double *thres, double *field,
+														 double cutoff) {
+	unsigned int size = (instance->network)->num_edges;
+
+	int min_pos = -1;
+	long double min_val = -1;
+
+	for (unsigned int i = 0; i < size; i++) {
+		if (!instance->fuses[i] && fabs(field[i]) > cutoff) {
+			double val = fabs(field[i]) / thres[i];
+			if (min_val < val) {
+				min_val = val; min_pos = i;
+			}
+		}
+	}
+
+	return min_pos;
+}
+
+break_data *fracture_network(finst *instance, double *fuse_thres,
+														 cholmod_common *c, double cutoff) {
+	unsigned int num_edges = instance->network->num_edges;
+	unsigned int num_verts = instance->network->num_verts;
+
+	break_data *breaking_data = alloc_break_data(num_edges);
+
+	while (true) {
+		double *voltages = get_voltage(instance, c);
+		double *field = get_current_v(instance, voltages, c);
+
+		double conductivity = get_conductivity(instance, voltages, c);
+		if (conductivity < 1e-12 && instance->voltage_bound) {
+			free(voltages);
+			free(field);
+			break;
+		}
+
+		int last_broke = inc_break_fuses(instance, fuse_thres, field, cutoff);
+		if (last_broke > num_edges || last_broke < -1) {
+			printf("%g \n", conductivity);
+			getchar();
+		}
+
+		update_break_data(breaking_data, last_broke, fabs(conductivity * fuse_thres[last_broke] / field[last_broke]), conductivity);
+
+		free(voltages);
+		free(field);
+
+		break_edge(instance, last_broke, c);
+
+		if (instance->num_components > 1) {
+			break;
+		}
+	}
+
+	return breaking_data;
+}
+
diff --git a/src/free_network.c b/src/free_network.c
new file mode 100644
index 0000000..651d3ae
--- /dev/null
+++ b/src/free_network.c
@@ -0,0 +1,22 @@
+
+#include "fracture.h"
+
+void free_fnet(fnet *network, cholmod_common *c) {
+	free(network->edges_to_verts);
+	if (network->edges_to_verts_break != network->edges_to_verts) {
+		free(network->edges_to_verts_break);
+	}
+	free(network->verts_to_edges_ind);
+	free(network->verts_to_edges);
+	free(network->bound_inds);
+	free(network->bound_verts);
+	free(network->vert_coords);
+	free(network->edge_coords);
+	free(network->dual_edges_to_verts);
+	free(network->dual_vert_coords);
+	free(network->dual_verts_to_edges_ind);
+	free(network->dual_verts_to_edges);
+	free(network->spanning_edges);
+	CHOL_F(free_sparse)(&(network->voltcurmat), c);
+	free(network);
+}
diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c
new file mode 100644
index 0000000..6f4263a
--- /dev/null
+++ b/src/gen_laplacian.c
@@ -0,0 +1,231 @@
+
+#include "fracture.h"
+
+cholmod_sparse *gen_adjacency(const finst *instance, bool dual, bool breakv,
+															unsigned int pad, cholmod_common *c) {
+	unsigned int ne;
+	if (dual)
+		ne = ((int)instance->network->num_edges -
+				(int)instance->num_remaining_edges);
+	else {
+		ne = instance->num_remaining_edges;
+		if (!breakv && instance->network->boundary != TORUS_BOUND) {
+			ne += instance->network->bound_inds[2];
+		}
+	}
+
+	unsigned int nnz = 2 * ne;
+
+	unsigned int nv;
+	if (dual)
+		nv = instance->network->num_dual_verts;
+	else {
+		if (breakv) nv = instance->network->num_verts_break;
+		else nv = instance->network->num_verts;
+		if (instance->network->boundary != TORUS_BOUND && !breakv) {
+			nv += 2;
+		}
+	}
+
+	cholmod_triplet *t =
+			CHOL_F(allocate_triplet)(nv + pad, nv + pad, nnz, 0, CHOLMOD_REAL, c);
+	CHOL_INT *ri = (CHOL_INT *)t->i;
+	CHOL_INT *ci = (CHOL_INT *)t->j;
+	double *ai = (double *)t->x;
+
+	t->nnz = nnz;
+
+	unsigned int *etv;
+	if (dual)
+		etv = instance->network->dual_edges_to_verts;
+	else {
+		if (breakv) etv = instance->network->edges_to_verts_break;
+		else etv = instance->network->edges_to_verts;
+	}
+
+	unsigned int count = 0;
+
+	for (unsigned int i = 0; i < instance->network->num_edges; i++) {
+		if ((instance->fuses[i] && dual) || (!instance->fuses[i] && !dual)) {
+			unsigned int v1 = etv[2 * i];
+			unsigned int v2 = etv[2 * i + 1];
+
+			ri[2 * count] = v1;
+			ci[2 * count] = v2;
+			ai[2 * count] = 1;
+
+			ri[2 * count + 1] = v2;
+			ci[2 * count + 1] = v1;
+			ai[2 * count + 1] = 1;
+
+			count++;
+		}
+	}
+
+	if (!breakv && instance->network->boundary != TORUS_BOUND && !dual) {
+		for (unsigned int i = 0; i < 2; i++) {
+			for (unsigned int j = 0; j < instance->network->bound_inds[i+1] - instance->network->bound_inds[i]; j++) {
+				ri[2*count] = instance->network->num_verts + i;
+				ci[2*count] = instance->network->bound_verts[instance->network->bound_inds[i] + j];
+				ai[2*count] = 1;
+				ri[2*count+1] = instance->network->bound_verts[instance->network->bound_inds[i] + j];
+				ci[2*count+1] = instance->network->num_verts + i;
+				ai[2*count+1] = 1;
+				count++;
+			}
+		}
+	}
+
+
+	cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c);
+	CHOL_F(free_triplet)(&t, c);
+
+	return s;
+}
+
+cholmod_sparse *gen_laplacian(const finst *instance, cholmod_common *c,
+															bool symmetric) {
+	const fnet *network = instance->network;
+	unsigned int num_verts = network->num_verts_break;
+	double *vert_coords = network->vert_coords;
+	unsigned int num_bounds = network->num_bounds;
+	double inf = instance->inf;
+	bool voltage_bound = instance->voltage_bound;
+	unsigned int *bound_inds = network->bound_inds;
+	unsigned int *bound_verts = network->bound_verts;
+
+	unsigned int num_gedges;
+	if (voltage_bound) {
+		num_gedges = 0;
+	} else {
+		num_gedges = network->bound_inds[num_bounds] + (num_bounds - 2) * 2;
+	}
+	if (network->boundary == TORUS_BOUND) num_gedges = bound_inds[1];
+	unsigned int num_gverts = network->break_dim;
+
+	unsigned int nnz = num_gverts + 2 * num_gedges;
+
+	cholmod_triplet *temp_m =
+			CHOL_F(allocate_triplet)(num_gverts, num_gverts, nnz, 0, CHOLMOD_REAL, c);
+
+	CHOL_INT *rowind = (CHOL_INT *)temp_m->i;
+	CHOL_INT *colind = (CHOL_INT *)temp_m->j;
+	double *acoo = (double *)temp_m->x;
+
+	temp_m->nnz = nnz;
+
+	for (unsigned int i = 0; i < num_gverts; i++) {
+		rowind[i] = i;
+		colind[i] = i;
+		acoo[i] = 0;
+	}
+
+	cholmod_sparse *adjacency = gen_adjacency(instance, false, true, (int)num_gverts - (int)num_verts, c);
+	CHOL_INT *ia = (CHOL_INT *)adjacency->p;
+	CHOL_INT *ja = (CHOL_INT *)adjacency->i;
+	double *aa = (double *)adjacency->x;
+
+	for (unsigned int i = 0; i < num_verts; i++) {
+		for (unsigned int j = 0; j < ia[i + 1] - ia[i]; j++) {
+			acoo[i] += aa[ia[i] + j];
+		}
+	}
+
+	if (network->boundary == TORUS_BOUND) {
+		for (unsigned int i = 0; i < network->bound_inds[1]; i++) {
+			unsigned int vv = network->bound_verts[i];
+			rowind[num_gverts + 2*i] = vv;
+			colind[num_gverts + 2*i] = network->num_verts + i;
+			acoo[num_gverts + 2*i] = -1;
+			rowind[num_gverts + 2*i+1] = network->num_verts + i;
+			colind[num_gverts + 2*i+1] = vv;
+			acoo[num_gverts + 2*i+1] = -1;
+			acoo[vv]++;
+			acoo[network->num_verts + i]++;
+		}
+		acoo[bound_verts[0]]++;
+	}
+
+	if (network->boundary != TORUS_BOUND) {
+		if (voltage_bound) {
+			for (unsigned int i = 0; i < 2; i++) {
+				for (unsigned int j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) {
+					acoo[bound_verts[bound_inds[i] + j]]++;
+				}
+			}
+		} else {
+			acoo[0]++;
+			for (unsigned int i = 0; i < num_bounds; i++) {
+				rowind[num_verts + i] = num_verts + i;
+				colind[num_verts + i] = num_verts + i;
+				if (i < 2)
+					acoo[num_verts + i] = (network->bound_inds[i + 1] -
+																 network->bound_inds[i] + network->num_bounds - 2) *
+																inf;
+				else
+					acoo[num_verts + i] =
+							(network->bound_inds[i + 1] - network->bound_inds[i] + 2) * inf;
+			}
+
+			unsigned int start = num_gverts;
+
+			for (unsigned int i = 0; i < num_bounds; i++) {
+				for (unsigned int j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) {
+					rowind[start + bound_inds[i] + j] = bound_verts[bound_inds[i] + j];
+					colind[start + bound_inds[i] + j] = num_verts + i;
+					acoo[start + bound_inds[i] + j] = -inf;
+					acoo[bound_verts[bound_inds[i] + j]] += inf;
+					colind[start + bound_inds[num_bounds] + bound_inds[i] + j] =
+							bound_verts[bound_inds[i] + j];
+					rowind[start + bound_inds[num_bounds] + bound_inds[i] + j] =
+							num_verts + i;
+					acoo[start + bound_inds[num_bounds] + bound_inds[i] + j] = -inf;
+				}
+			}
+
+			for (unsigned int i = 0; i < num_bounds - 2; i++) {
+				rowind[nnz - 2 * i - 1] = num_verts;
+				colind[nnz - 2 * i - 1] = num_verts + 2 + i;
+				acoo[nnz - 2 * i - 1] = -inf;
+				rowind[nnz - 2 * i - 2] = num_verts + 1;
+				colind[nnz - 2 * i - 2] = num_verts + 2 + i;
+				acoo[nnz - 2 * i - 2] = -inf;
+				colind[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = num_verts;
+				rowind[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = num_verts + 2 + i;
+				acoo[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = -inf;
+				colind[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = num_verts + 1;
+				rowind[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = num_verts + 2 + i;
+				acoo[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = -inf;
+			}
+		}
+	}
+
+
+	for (unsigned int i = 0; i < num_gverts; i++) {
+		if (acoo[i] == 0)
+			acoo[i] = 1;
+	}
+
+	assert(CHOL_F(check_triplet)(temp_m, c));
+
+	cholmod_sparse *t_out = CHOL_F(triplet_to_sparse)(temp_m, nnz, c);
+	assert(CHOL_F(check_sparse)(t_out, c));
+
+	double alpha[2] = {1, 0};
+	double beta[2] = {-1, 0};
+	cholmod_sparse *t_laplacian =
+			CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c);
+
+	cholmod_sparse *laplacian;
+	if (symmetric)
+		laplacian = CHOL_F(copy)(t_laplacian, 1, 1, c);
+	else
+		laplacian = CHOL_F(copy)(t_laplacian, 0, 1, c);
+
+	CHOL_F(free_sparse)(&t_out, c);
+	CHOL_F(free_sparse)(&adjacency, c);
+	CHOL_F(free_sparse)(&t_laplacian, c);
+	CHOL_F(free_triplet)(&temp_m, c);
+
+	return laplacian;
+}
diff --git a/src/gen_voltcurmat.c b/src/gen_voltcurmat.c
new file mode 100644
index 0000000..ab16150
--- /dev/null
+++ b/src/gen_voltcurmat.c
@@ -0,0 +1,36 @@
+
+#include "fracture.h"
+
+cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts,
+															 unsigned int *edges_to_verts,
+															 cholmod_common *c) {
+
+	cholmod_triplet *t_m = CHOL_F(allocate_triplet)(
+			num_edges, num_verts, num_edges * 2, 0, CHOLMOD_REAL, c);
+	assert(t_m != NULL);
+
+	CHOL_INT *rowind = (CHOL_INT *)t_m->i;
+	CHOL_INT *colind = (CHOL_INT *)t_m->j;
+	double *acoo = (double *)t_m->x;
+
+	for (unsigned int i = 0; i < num_edges; i++) {
+		unsigned int v1 = edges_to_verts[2 * i];
+		unsigned int v2 = edges_to_verts[2 * i + 1];
+		rowind[2 * i] = i;
+		rowind[2 * i + 1] = i;
+		colind[2 * i] = v1;
+		colind[2 * i + 1] = v2;
+		acoo[2 * i] = 1;
+		acoo[2 * i + 1] = -1;
+	}
+
+	t_m->nnz = num_edges * 2;
+
+	assert(CHOL_F(check_triplet)(t_m, c));
+
+	cholmod_sparse *m = CHOL_F(triplet_to_sparse)(t_m, num_edges * 2, c);
+
+	CHOL_F(free_triplet)(&t_m, c);
+
+	return m;
+}
diff --git a/src/geometry.c b/src/geometry.c
new file mode 100644
index 0000000..83ef0bf
--- /dev/null
+++ b/src/geometry.c
@@ -0,0 +1,55 @@
+
+#include "fracture.h"
+
+int edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
+									bool index) {
+	assert(edge < pow(width, 2));
+
+	int x = edge / width + 1;
+	int y = edge % width + 1;
+
+	if (periodic) {
+		return (((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) - 1) %
+								width +
+						(x - index) * width) /
+					 2;
+	} else {
+		return ((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) +
+						(x - index) * (width + 1) - 1) /
+					 2;
+	}
+}
+
+int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
+											 bool index) {
+	assert(edge < pow(width, 2));
+
+	int x = edge / width + 1;
+	int y = edge % width + 1;
+
+	if (periodic) {
+		return (((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) - 1) %
+								width +
+						(x - index) * width) /
+					 2;
+	} else {
+		return ((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) +
+						(x - index) * (width + 1) - 1) /
+					 2;
+	}
+}
+
+double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert,
+													bool index) {
+	if (periodic) {
+		if (index)
+			return vert % (width / 2) + (vert / (width / 2)) % 2;
+		else
+			return vert / (width / 2);
+	} else {
+		if (index)
+			return (2 * vert) % (width + 1);
+		else
+			return (2 * vert) / (width + 1);
+	}
+}
diff --git a/src/get_conductivity.c b/src/get_conductivity.c
new file mode 100644
index 0000000..7ef08a8
--- /dev/null
+++ b/src/get_conductivity.c
@@ -0,0 +1,24 @@
+
+#include "fracture.h"
+
+double get_conductivity(finst *inst, double *voltage, cholmod_common *c) {
+	if (inst->voltage_bound) {
+		double tot_cur = 0;
+		for (unsigned int i = 0; i < inst->network->num_spanning_edges; i++) {
+			unsigned int e = inst->network->spanning_edges[i];
+			if (!inst->fuses[e]) {
+				unsigned int v1 = inst->network->edges_to_verts[2*e];
+				unsigned int v2 = inst->network->edges_to_verts[2*e+1];
+				double v1y = inst->network->vert_coords[2 * v1 + 1];
+				double v2y = inst->network->vert_coords[2 * v2 + 1];
+				unsigned int s1 = v1y < v2y ? v1 : v2;
+				unsigned int s2 = v1y < v2y ? v2 : v1;
+				tot_cur += voltage[s1] - voltage[s2];
+			}
+		}
+
+		return fabs(tot_cur);
+	} else {
+		return 1 / fabs(voltage[inst->network->num_verts] - voltage[inst->network->num_verts + 1]);
+	}
+}
diff --git a/src/get_current.c b/src/get_current.c
new file mode 100644
index 0000000..5360054
--- /dev/null
+++ b/src/get_current.c
@@ -0,0 +1,99 @@
+
+#include "fracture.h"
+
+double *get_voltage(const finst *instance, cholmod_common *c) {
+	cholmod_dense *b = instance->boundary_cond;
+	cholmod_factor *factor = instance->factor;
+
+	cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
+
+	if (((double *)x->x)[0] != ((double *)x->x)[0]) {
+		printf("ERROR: GET_VOLTAGE FAILED\n\n");
+		CHOL_F(free_dense)(&x, c);
+		return NULL;
+	}
+
+	double *field = (double *)x->x;
+	x->x = NULL;
+
+	CHOL_F(free_dense)(&x, c);
+
+
+	return field;
+}
+
+double *get_current(const finst *instance, cholmod_common *c) {
+	unsigned int num_edges = instance->network->num_edges;
+	unsigned int num_verts = instance->network->num_verts_break;
+	unsigned int num_gverts = instance->network->break_dim;
+	unsigned int num_bounds = instance->network->num_bounds;
+	cholmod_sparse *voltcurmat = instance->network->voltcurmat;
+
+	double *voltages = get_voltage(instance, c);
+	if (voltages == NULL) {
+		return NULL;
+	}
+	cholmod_dense *x = CHOL_F(allocate_dense)(
+			num_gverts, 1, num_gverts, CHOLMOD_REAL, c);
+	double *tmp_x = x->x;
+	x->x = voltages;
+
+	cholmod_dense *y =
+			CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c);
+
+	double alpha[2] = {1, 0};
+	double beta[2] = {0, 0};
+	CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+
+	double *field = (double *)malloc(num_edges * sizeof(double));
+
+	for (int i = 0; i < num_edges; i++) {
+		if (instance->fuses[i])
+			field[i] = 0;
+		else
+			field[i] = ((double *)y->x)[i];
+	}
+
+	x->x = tmp_x;
+	free(voltages);
+	CHOL_F(free_dense)(&x, c);
+	CHOL_F(free_dense)(&y, c);
+
+	return field;
+}
+
+
+double *get_current_v(const finst *instance, double *voltages, cholmod_common *c) {
+	unsigned int num_edges = instance->network->num_edges;
+	unsigned int num_verts = instance->network->num_verts_break;
+	unsigned int num_gverts = instance->network->break_dim;
+	unsigned int num_bounds = instance->network->num_bounds;
+	cholmod_sparse *voltcurmat = instance->network->voltcurmat;
+
+	cholmod_dense *x = CHOL_F(allocate_dense)(
+			num_gverts, 1, num_gverts, CHOLMOD_REAL, c);
+	double *tmp_x = x->x;
+	x->x = voltages;
+
+	cholmod_dense *y =
+			CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c);
+
+	double alpha[2] = {1, 0};
+	double beta[2] = {0, 0};
+	CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+
+	double *field = (double *)malloc(num_edges * sizeof(double));
+
+	for (int i = 0; i < num_edges; i++) {
+		if (instance->fuses[i])
+			field[i] = 0;
+		else
+			field[i] = ((double *)y->x)[i];
+	}
+
+	x->x = tmp_x;
+	CHOL_F(free_dense)(&x, c);
+	CHOL_F(free_dense)(&y, c);
+
+	return field;
+}
diff --git a/src/get_dual_clusters.c b/src/get_dual_clusters.c
new file mode 100644
index 0000000..488793f
--- /dev/null
+++ b/src/get_dual_clusters.c
@@ -0,0 +1,36 @@
+
+#include "fracture.h"
+
+unsigned int *get_clusters(finst *instance, cholmod_common *c) {
+	cholmod_sparse *s_dual = gen_adjacency(instance, true, false, 0, c);
+
+	unsigned int *dual_marks = find_components(s_dual, 0);
+	CHOL_F(free_sparse)(&s_dual, c);
+
+	return dual_marks;
+}
+
+unsigned int *get_cluster_dist(finst *instance, cholmod_common *c) {
+	unsigned int *clusters = get_clusters(instance, c);
+	unsigned int *cluster_dist = (unsigned int *)calloc(
+			instance->network->num_dual_verts, sizeof(unsigned int));
+
+	unsigned int cur_mark = 0;
+	while (true) {
+		cur_mark++;
+		unsigned int num_in_cluster = 0;
+		for (unsigned int i = 0; i < instance->network->num_dual_verts; i++) {
+			if (clusters[i] == cur_mark)
+				num_in_cluster++;
+		}
+
+		if (num_in_cluster == 0)
+			break;
+
+		cluster_dist[num_in_cluster - 1]++;
+	}
+
+	free(clusters);
+
+	return cluster_dist;
+}
diff --git a/src/get_file.c b/src/get_file.c
new file mode 100644
index 0000000..9141631
--- /dev/null
+++ b/src/get_file.c
@@ -0,0 +1,38 @@
+
+#include "fracture.h"
+
+FILE *get_file(const char *prefix, unsigned int width, unsigned int crack,
+							 double beta, unsigned int iter, unsigned int num_iter,
+							 unsigned int num, bool read) {
+	int prefix_len = strlen(prefix);
+	int width_len = 1 + (int)log10(width);
+	int crack_len;
+	if (crack != 0) {
+		crack_len = 1 + (int)log10(crack);
+	} else {
+		crack_len = 1;
+	}
+	int beta_len;
+	if (beta > 1) {
+		beta_len = 1 + (int)log10(beta) + 3;
+	} else {
+		beta_len = 4;
+	}
+	int iter_len = 1 + (int)log10(num_iter);
+	int num_len = 1 + (int)log10(num);
+	int num_underscores = 5;
+
+	int len = prefix_len + width_len + crack_len + beta_len + iter_len + num_len +
+						num_underscores + 4;
+
+	char filename[len + 1];
+	snprintf(filename, sizeof(filename), "%s_%u_%u_%.2f_%0*u_%u.txt", prefix,
+					 width, crack, beta, iter_len, iter + 1, num);
+
+	char *mode = "w";
+	if (read) {
+		mode = "r";
+	}
+	FILE *file = fopen(filename, mode);
+	return file;
+}
diff --git a/src/graph_components.c b/src/graph_components.c
new file mode 100644
index 0000000..036f2fe
--- /dev/null
+++ b/src/graph_components.c
@@ -0,0 +1,185 @@
+
+#include "fracture.h"
+
+bool set_connected(const cholmod_sparse *laplacian, unsigned int *marks,
+									 int vertex, int label, int stop_at, int exclude) {
+	unsigned int num_verts = (CHOL_INT)laplacian->nrow;
+
+	CHOL_INT *ia = (CHOL_INT *)laplacian->p;
+	CHOL_INT *ja = (CHOL_INT *)laplacian->i;
+	double *a = (double *)laplacian->x;
+
+	bool stopped = false;
+
+	int *queue = (int *)malloc((num_verts) * sizeof(int));
+	assert(queue != NULL);
+
+	int queue_pos = 0;
+
+	int *component_list = (int *)malloc((num_verts) * sizeof(int));
+	assert(component_list != NULL);
+
+	int list_pos = 0;
+	queue[0] = vertex;
+	component_list[0] = vertex;
+
+	unsigned int old_label = marks[vertex];
+	marks[vertex] = label;
+
+	while (queue_pos >= 0) {
+		int next_vertex = queue[queue_pos];
+		queue_pos--;
+		for (int j = 0; j < ia[next_vertex + 1] - ia[next_vertex]; j++) {
+			if (ja[ia[next_vertex] + j] == stop_at && a[ia[next_vertex] + j] != 0) {
+				// if we run into the other vertex, the components haven't changed
+				stopped = true; queue_pos = -1;
+				break;
+			}
+			if (marks[ja[ia[next_vertex] + j]] != label &&
+					a[ia[next_vertex] + j] != 0 &&
+					ja[ia[next_vertex] + j] < num_verts - exclude) {
+				marks[ja[ia[next_vertex] + j]] = label;
+				queue_pos++;
+				queue[queue_pos] = ja[ia[next_vertex] + j];
+				list_pos++;
+				component_list[list_pos] = ja[ia[next_vertex] + j];
+			}
+		}
+	}
+
+	if (stopped) {
+		for (unsigned int i = 0; i <= list_pos; i++) {
+			marks[component_list[i]] = old_label;
+		}
+	}
+
+	free(queue);
+	free(component_list);
+
+	if (stopped)
+		return false;
+	else
+		return true;
+}
+
+int update_components(const cholmod_sparse *laplacian, unsigned int *marks,
+											int old_num_components, int v1, int v2, int exclude) {
+	assert(laplacian != NULL);
+	assert(marks != NULL);
+
+	if (set_connected(laplacian, marks, v1, old_num_components + 1, v2,
+										exclude)) {
+		return old_num_components + 1;
+	} else
+		return old_num_components;
+}
+
+unsigned int *find_components(const cholmod_sparse *laplacian, unsigned int skip) {
+	size_t num_verts = laplacian->nrow;
+	unsigned int *marks = (unsigned int *)calloc(num_verts, sizeof(unsigned int));
+
+	int num_components = 0;
+
+	for (int i = 0; i < num_verts; i++) {
+		if (marks[i] == 0) {
+			num_components++;
+			set_connected(laplacian, marks, i, num_components, -1, skip);
+		}
+	}
+
+	return marks;
+}
+
+unsigned int find_cycles(unsigned int num_edges, const bool *fuses, const unsigned int *ev, const unsigned
+		int *vei, const unsigned int *ve, int **cycles) {
+	unsigned int num_cycles = 0;
+	unsigned int *elist = (unsigned int *)malloc(num_edges * sizeof(unsigned int));
+	unsigned int *side = (unsigned int *)calloc(num_edges, sizeof(unsigned int));
+	unsigned int *num = (unsigned int *)malloc(num_edges * sizeof(unsigned int));
+	for (unsigned int i = 0; i < num_edges; i++) {
+		if (fuses[i]) {
+		bool in_cycle = false;
+		for (unsigned int j = 0; j < num_cycles; j++) {
+			unsigned int k = 0;
+			int e;
+			while ((e = cycles[j][k]) >= 0) {
+				if (e == i) {
+					in_cycle = true;
+					break;
+				}
+				k++;
+			}
+			if (in_cycle) {
+				break;
+			}
+		}
+		if (!in_cycle) {
+			unsigned int pos = 0;
+			unsigned int cur_e = i;
+			bool *included = (bool *)calloc(num_edges, sizeof(bool));
+			included[cur_e] = true;
+			elist[0] = cur_e;
+			side[0] = 1;
+			num[0] = 0;
+			while (true) {
+				unsigned int cv = ev[2 * cur_e + side[pos]];
+				unsigned int new_e = cur_e;
+				unsigned int j;
+				for (j = num[pos]; j < vei[cv+1] - vei[cv]; j++) {
+					new_e = ve[vei[cv] + j];
+					if (new_e != cur_e && fuses[new_e]) break;
+				}
+				if (new_e != cur_e && fuses[new_e]) {
+					if (new_e == elist[0] && cv == ev[2 * elist[0] + !side[0]]) {
+						pos++;
+						cycles[2*num_cycles] = (int *)malloc((pos + 1) * sizeof(int));
+						cycles[2*num_cycles+1] = (int *)malloc((pos + 1) * sizeof(int));
+						for (unsigned int i = 0; i < pos; i++) {
+							cycles[2*num_cycles][i] = elist[i];
+							cycles[2*num_cycles+1][i] = side[i];
+						}
+						cycles[2*num_cycles][pos] = -1;
+						cycles[2*num_cycles+1][pos] = -1;
+						num_cycles++;
+						pos--;
+						num[pos]++;
+						included[pos] = false;
+						cur_e = elist[pos];
+					} else if (included[new_e]) {
+						if (pos == 0) break;
+						included[cur_e] = false;
+						pos--;
+						num[pos]++;
+						cur_e = elist[pos];
+					} else {
+						num[pos] = j;
+						pos++;
+						elist[pos] = new_e;
+						included[new_e] = true;
+						unsigned int nv1 = ev[2 * new_e];
+						unsigned int nv2 = ev[2 * new_e + 1];
+						unsigned int v = nv1 == cv ? nv2 : nv1;
+						side[pos] = v == nv1 ? 0 : 1;
+						num[pos] = 0;
+						cur_e = new_e;
+					}
+				} else {
+					if (pos == 0) break;
+					included[cur_e] = false;
+					pos--;
+					num[pos]++;
+					cur_e = elist[pos];
+				}
+			}
+			free(included);
+		}
+		}
+	}
+
+	free(elist);
+	free(side);
+	free(num);
+
+	return num_cycles;
+}
+
diff --git a/src/homo_square_fracture.c b/src/homo_square_fracture.c
new file mode 100644
index 0000000..a0dd305
--- /dev/null
+++ b/src/homo_square_fracture.c
@@ -0,0 +1,418 @@
+
+
+#include "fracture.h"
+
+int main(int argc, char *argv[]) {
+	int opt;
+
+	// defining variables to be (potentially) set by command line flags
+	unsigned int N, L, filename_len;
+	double beta, inf, cutoff;
+	boundary_type boundary;
+	bool include_breaking, save_cluster_dist, use_voltage_boundaries, save_network,
+			 save_crit_stress, save_toughness, save_corr, save_conductivity,
+			 save_damage;
+
+	filename_len = 100;
+
+	N = 100;
+	L = 16;
+	beta = .3;
+	inf = 1e-8;
+	cutoff = 1e-10;
+	boundary = FREE_BOUND;
+	include_breaking = false;
+	save_cluster_dist = false;
+	use_voltage_boundaries = false;
+	save_network = false;
+	save_crit_stress = false;
+	save_damage = false;
+	save_corr = false;
+	save_conductivity = false;
+	save_toughness = false;
+
+	int boundary_int;
+	char boundc2 = 'f';
+
+	while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrt")) != -1) {
+		switch (opt) {
+			case 'n':
+				N = atoi(optarg);
+				break;
+			case 'L':
+				L = atoi(optarg);
+				break;
+			case 'b':
+				beta = atof(optarg);
+				break;
+			case 'B':
+				boundary_int = atoi(optarg);
+				switch (boundary_int) {
+					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;
+					default:
+						printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n");
+				}
+				break;
+			case 'd':
+				save_damage = true;
+				break;
+			case 'V':
+				use_voltage_boundaries = true;
+				break;
+			case 'c':
+				save_cluster_dist = true;
+				break;
+			case 'o':
+				include_breaking = true;
+				break;
+			case 'N':
+				save_network = true;
+				break;
+			case 's':
+				save_crit_stress = true;
+				break;
+			case 'C':
+				save_corr = true;
+				break;
+			case 'r':
+				save_conductivity = true;
+				inf = 1;
+				break;
+			case 't':
+				save_toughness = true;
+				break;
+			default: /* '?' */
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	char boundc;
+	if (use_voltage_boundaries) boundc = 'v';
+	else boundc = 'c';
+
+	FILE *break_out;
+	if (include_breaking) {
+		char *break_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(break_filename, filename_len, "breaks_s_%c_%c_%u_%g.txt", boundc, boundc2, L, beta);
+		break_out = fopen(break_filename, "a");
+		free(break_filename);
+	}
+
+	// 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;
+	}
+
+
+	fnet *network = ini_square_network(L, boundary, false, &c);
+	finst *perm_instance = create_instance(network, inf, use_voltage_boundaries, true, &c);
+	unsigned int c_dist_size = network->num_dual_verts;
+	unsigned int a_dist_size = network->num_verts;
+
+	// define arrays for saving cluster and avalanche distributions
+	unsigned int *cluster_size_dist;
+	unsigned int *avalanche_size_dist;
+	char *c_filename;
+	if (save_cluster_dist) {
+		cluster_size_dist =
+				(unsigned int *)calloc(c_dist_size, sizeof(unsigned int));
+		avalanche_size_dist =
+				(unsigned int *)calloc(a_dist_size, sizeof(unsigned int));
+
+		c_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(c_filename, filename_len, "cstr_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
+
+		FILE *cluster_out = fopen(c_filename, "r");
+
+		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;
+			}
+			fclose(cluster_out);
+		}
+	}
+
+	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(c_dist_size, sizeof(double));
+	}
+
+	double *conductivity;
+	if (save_conductivity) {
+		conductivity = (double *)malloc(N * sizeof(double));
+	}
+
+	double *damage;
+	if (save_damage) {
+		damage = (double *)malloc(N * sizeof(double));
+	}
+
+	double *toughness;
+	if (save_toughness) {
+		toughness = (double *)malloc(N * sizeof(double));
+	}
+
+
+	printf("\n");
+	for (unsigned int i = 0; i < N; i++) {
+		printf("\033[F\033[JFRACTURE: %0*d / %d\n", (int)log10(N) + 1, i + 1, N);
+
+		double *fuse_thres = gen_fuse_thres(network->num_edges, network->edge_coords, beta, beta_scaling_flat);
+		finst *instance = copy_instance(perm_instance, &c);
+		break_data *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
+		free_instance(instance, &c);
+		free(fuse_thres);
+
+		unsigned int max_pos = 0;
+		double max_val = 0;
+
+		for (unsigned int j = 0; j < breaking_data->num_broken; j++) {
+			double val = breaking_data->extern_field[j];
+
+			if (val > max_val) {
+				max_pos = j;
+				max_val = val;
+			}
+		}
+
+		if (save_crit_stress) {
+			crit_stress[i] =
+					breaking_data->extern_field[max_pos];
+		}
+
+		finst *tmp_instance = copy_instance(perm_instance, &c);
+
+		unsigned int av_size = 0;
+		double cur_val = DBL_MAX;
+		for (unsigned int j = 0; j < max_pos; j++) {
+			break_edge(tmp_instance, breaking_data->break_list[j], &c);
+
+			double val = fabs(breaking_data->extern_field[j]);
+			if (save_cluster_dist) {
+				if (val < cur_val) {
+					av_size++;
+				}
+
+				if (val > cur_val) {
+					avalanche_size_dist[av_size]++;
+					av_size = 0;
+					cur_val = val;
+				}
+			}
+		}
+
+		if (save_conductivity) {
+			if (!use_voltage_boundaries) {
+				double *tmp_voltage = get_voltage(tmp_instance, &c);
+				conductivity[i] = 1/fabs(tmp_voltage[tmp_instance->network->num_verts + 1] - tmp_voltage[tmp_instance->network->num_verts]);
+				free(tmp_voltage);
+			} else {
+				conductivity[i] = breaking_data->conductivity[max_pos];
+			}
+		}
+
+		if (save_toughness) {
+			double tmp_toughness = 0;
+			if (max_pos > 0) {
+				double sigma1 = breaking_data->extern_field[0];
+				double epsilon1 = sigma1 / breaking_data->conductivity[0];
+				for (unsigned int j = 0; j < max_pos - 1; j++) {
+					double sigma2 = breaking_data->extern_field[j+1];
+					double epsilon2 = sigma2 / breaking_data->conductivity[j+1];
+					if (epsilon2 > epsilon1) {
+						tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2;
+						sigma1 = sigma2; epsilon1 = epsilon2;
+					}
+				}
+			}
+			toughness[i] = tmp_toughness;
+		}
+
+		if (save_damage) {
+			damage[i] = ((double)max_pos) / tmp_instance->network->num_edges;
+		}
+
+		if (save_cluster_dist) {
+			unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c);
+			for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) {
+				cluster_size_dist[j] += tmp_cluster_dist[j];
+			}
+			free(tmp_cluster_dist);
+		}
+
+		if (save_corr) {
+			double *tmp_corr = get_corr(tmp_instance, dists, &c);
+			for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) {
+				avg_corr[i] += tmp_corr[j] / N;
+			}
+			free(tmp_corr);
+		}
+
+		if (save_network) {
+			FILE *net_out = fopen("network.txt", "w");
+			for (unsigned int j = 0; j < network->num_verts; j++) {
+				fprintf(net_out, "%f %f ", network->vert_coords[2 * j],
+								tmp_instance->network->vert_coords[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) {
+				fprintf(net_out, "%u %u ", tmp_instance->network->edges_to_verts[2 * j],
+								tmp_instance->network->edges_to_verts[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) {
+				fprintf(net_out, "%f %f ", tmp_instance->network->dual_vert_coords[2 * j],
+								tmp_instance->network->dual_vert_coords[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) {
+				fprintf(net_out, "%u %u ", tmp_instance->network->dual_edges_to_verts[2 * j],
+								tmp_instance->network->dual_edges_to_verts[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) {
+				fprintf(net_out, "%d ", tmp_instance->fuses[j]);
+			}
+			fclose(net_out);
+		}
+
+		free_instance(tmp_instance, &c);
+
+		if (include_breaking) {
+			for (unsigned int j = 0; j < breaking_data->num_broken; j++) {
+				fprintf(break_out, "%u %g %g ", breaking_data->break_list[j],
+								breaking_data->extern_field[j], breaking_data->conductivity[j]);
+			}
+			fprintf(break_out, "\n");
+		}
+
+		free_break_data(breaking_data);
+	}
+
+	printf("\033[F\033[JFRACTURE: COMPLETE\n");
+
+	free_instance(perm_instance, &c);
+	free_fnet(network, &c);
+
+	if (save_cluster_dist) {
+		FILE *cluster_out = fopen(c_filename, "w");
+
+		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);
+
+		free(c_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_s_%c_%c_%d_%g.txt", boundc, boundc2, L,
+						 beta);
+		FILE *corr_file = fopen(corr_filename, "w");
+		for (unsigned int i = 0; i < c_dist_size; i++) {
+			fprintf(corr_file, "%g ", avg_corr[i]);
+		}
+		fclose(corr_file);
+		free(corr_filename);
+		free(avg_corr);
+	}
+
+	if (save_conductivity) {
+		char *cond_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(cond_filename, filename_len, "cond_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
+		FILE *cond_file = fopen(cond_filename, "a");
+		for (unsigned int i = 0; i < N; i++) {
+			fprintf(cond_file, "%g ", conductivity[i]);
+		}
+		fclose(cond_file);
+		free(cond_filename);
+		free(conductivity);
+	}
+
+	if (save_toughness) {
+		char *tough_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(tough_filename, filename_len, "tuff_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
+		FILE *tough_file = fopen(tough_filename, "a");
+		for (unsigned int i = 0; i < N; i++) {
+			fprintf(tough_file, "%g ", toughness[i]);
+		}
+		fclose(tough_file);
+		free(tough_filename);
+		free(toughness);
+	}
+
+	if (save_damage) {
+		char *hdam_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(hdam_filename, filename_len, "damg_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
+		FILE *hdam_file = fopen(hdam_filename, "a");
+		for (unsigned int i = 0; i < N; i++) {
+			fprintf(hdam_file, "%g ", damage[i]);
+		}
+		fclose(hdam_file);
+		free(hdam_filename);
+		free(damage);
+	}
+
+	if (include_breaking) {
+		fclose(break_out);
+	}
+
+	if (save_crit_stress) {
+		char *a_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(a_filename, filename_len, "strs_s_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
+		FILE *a_file = fopen(a_filename, "a");
+		for (int i = 0; i < N; i++) {
+			fprintf(a_file, "%g ", crit_stress[i]);
+		}
+		fclose(a_file);
+		free(a_filename);
+		free(crit_stress);
+	}
+
+	CHOL_F(finish)(&c);
+
+	return 0;
+}
diff --git a/src/homo_voronoi_fracture.c b/src/homo_voronoi_fracture.c
new file mode 100644
index 0000000..1d9d0e6
--- /dev/null
+++ b/src/homo_voronoi_fracture.c
@@ -0,0 +1,420 @@
+
+
+#include "fracture.h"
+
+int main(int argc, char *argv[]) {
+	int opt;
+
+	// defining variables to be (potentially) set by command line flags
+	unsigned int N, L, filename_len;
+	double beta, inf, cutoff;
+	boundary_type boundary;
+	bool include_breaking, save_cluster_dist, use_voltage_boundaries, save_network,
+			 save_crit_stress, save_toughness, save_corr, save_conductivity,
+			 save_damage;
+
+	filename_len = 100;
+
+	N = 100;
+	L = 16;
+	beta = .3;
+	inf = 1e10;
+	cutoff = 1e-9;
+	boundary = FREE_BOUND;
+	include_breaking = false;
+	save_cluster_dist = false;
+	use_voltage_boundaries = false;
+	save_network = false;
+	save_crit_stress = false;
+	save_damage = false;
+	save_corr = false;
+	save_conductivity = false;
+	save_toughness = false;
+
+	int boundary_int;
+	char boundc2 = 'f';
+
+	while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrt")) != -1) {
+		switch (opt) {
+			case 'n':
+				N = atoi(optarg);
+				break;
+			case 'L':
+				L = atoi(optarg);
+				break;
+			case 'b':
+				beta = atof(optarg);
+				break;
+			case 'B':
+				boundary_int = atoi(optarg);
+				switch (boundary_int) {
+					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;
+					default:
+						printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n");
+				}
+				break;
+			case 'd':
+				save_damage = true;
+				break;
+			case 'V':
+				use_voltage_boundaries = true;
+				break;
+			case 'c':
+				save_cluster_dist = true;
+				break;
+			case 'o':
+				include_breaking = true;
+				break;
+			case 'N':
+				save_network = true;
+				break;
+			case 's':
+				save_crit_stress = true;
+				break;
+			case 'C':
+				save_corr = true;
+				break;
+			case 'r':
+				save_conductivity = true;
+				//inf = 1;
+				break;
+			case 't':
+				save_toughness = true;
+				break;
+			default: /* '?' */
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	char boundc;
+	if (use_voltage_boundaries) boundc = 'v';
+	else boundc = 'c';
+
+	FILE *break_out;
+	if (include_breaking) {
+		char *break_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(break_filename, filename_len, "breaks_v_%c_%c_%u_%g.txt", boundc, boundc2, L, beta);
+		break_out = fopen(break_filename, "a");
+		free(break_filename);
+	}
+
+	unsigned int voronoi_max_verts = 4 * pow(L, 2);
+	unsigned int c_dist_size = voronoi_max_verts;
+	unsigned int a_dist_size = voronoi_max_verts;
+
+	// define arrays for saving cluster and avalanche distributions
+	unsigned int *cluster_size_dist;
+	unsigned int *avalanche_size_dist;
+	char *c_filename;
+	if (save_cluster_dist) {
+		cluster_size_dist =
+				(unsigned int *)calloc(c_dist_size, sizeof(unsigned int));
+		avalanche_size_dist =
+				(unsigned int *)calloc(a_dist_size, sizeof(unsigned int));
+
+		c_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(c_filename, filename_len, "cstr_v_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
+
+		FILE *cluster_out = fopen(c_filename, "r");
+
+		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;
+			}
+			fclose(cluster_out);
+		}
+	}
+
+	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(voronoi_max_verts, sizeof(double));
+	}
+
+	double *conductivity;
+	if (save_conductivity) {
+		conductivity = (double *)malloc(N * sizeof(double));
+	}
+
+	double *damage;
+	if (save_damage) {
+		damage = (double *)malloc(N * sizeof(double));
+	}
+
+	double *toughness;
+	if (save_toughness) {
+		toughness = (double *)malloc(N * sizeof(double));
+	}
+
+
+	// start cholmod
+	cholmod_common c;
+	CHOL_F(start)(&c);
+
+	/* if we use voltage boundary conditions, the laplacian matrix is positive
+	 * definite and we can use a supernodal LL decomposition.  otherwise we need
+	 * to use the simplicial LDL decomposition
+	 */
+	if (use_voltage_boundaries) {
+		//(&c)->supernodal = CHOLMOD_SUPERNODAL;
+		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
+	} else {
+		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
+	}
+
+
+	printf("\n");
+	for (unsigned int i = 0; i < N; i++) {
+		printf("\033[F\033[JFRACTURE: %0*d / %d\n", (int)log10(N) + 1, i + 1, N);
+
+		fnet *network = ini_voronoi_network(L, boundary, genfunc_hyperuniform, &c);
+		finst *perm_instance = create_instance(network, inf, use_voltage_boundaries, true, &c);
+		double *fuse_thres = gen_fuse_thres(network->num_edges, network->edge_coords, beta, beta_scaling_flat);
+		finst *instance = copy_instance(perm_instance, &c);
+		break_data *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
+		bool debug_stop = instance->debug_stop;
+		free_instance(instance, &c);
+		free(fuse_thres);
+
+		unsigned int max_pos = 0;
+		double max_val = 0;
+
+		for (unsigned int j = 0; j < breaking_data->num_broken; j++) {
+			double val = breaking_data->extern_field[j];
+
+			if (val > max_val) {
+				max_pos = j;
+				max_val = val;
+			}
+		}
+
+		if (save_crit_stress) {
+			crit_stress[i] =
+					breaking_data->extern_field[max_pos];
+		}
+
+		finst *tmp_instance = copy_instance(perm_instance, &c);
+
+		unsigned int av_size = 0;
+		double cur_val = 0;
+		for (unsigned int j = 0; j < max_pos; j++) {
+			break_edge(tmp_instance, breaking_data->break_list[j], &c);
+
+			double val = breaking_data->extern_field[j];
+			if (save_cluster_dist) {
+				if (val < cur_val) {
+					av_size++;
+				}
+
+				if (val > cur_val) {
+					avalanche_size_dist[av_size]++;
+					av_size = 0;
+					cur_val = val;
+				}
+			}
+		}
+
+		if (save_conductivity) {
+			if (!use_voltage_boundaries) {
+				double *tmp_voltage = get_voltage(tmp_instance, &c);
+				conductivity[i] = 1/fabs(tmp_voltage[tmp_instance->network->num_verts + 1] - tmp_voltage[tmp_instance->network->num_verts]);
+				free(tmp_voltage);
+			} else {
+				conductivity[i] = breaking_data->conductivity[max_pos];
+			}
+		}
+
+		if (save_toughness) {
+			double tmp_toughness = 0;
+			if (max_pos > 0) {
+				double sigma1 = breaking_data->extern_field[0];
+				double epsilon1 = sigma1 / breaking_data->conductivity[0];
+				for (unsigned int j = 0; j < max_pos - 1; j++) {
+					double sigma2 = breaking_data->extern_field[j+1];
+					double epsilon2 = sigma2 / breaking_data->conductivity[j+1];
+					if (epsilon2 > epsilon1) {
+						tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2;
+						sigma1 = sigma2; epsilon1 = epsilon2;
+					}
+				}
+			}
+			toughness[i] = tmp_toughness;
+		}
+
+		if (save_damage) {
+			damage[i] = ((double)max_pos) / tmp_instance->network->num_edges;
+		}
+
+		if (save_cluster_dist) {
+			unsigned int *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c);
+			for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) {
+				cluster_size_dist[j] += tmp_cluster_dist[j];
+			}
+			free(tmp_cluster_dist);
+		}
+
+		if (save_corr) {
+			double *tmp_corr = get_corr(tmp_instance, dists, &c);
+			for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) {
+				avg_corr[i] += tmp_corr[j] / N;
+			}
+			free(tmp_corr);
+		}
+
+		if (save_network) {
+			FILE *net_out = fopen("network.txt", "w");
+			for (unsigned int j = 0; j < network->num_verts; j++) {
+				fprintf(net_out, "%f %f ", network->vert_coords[2 * j],
+								tmp_instance->network->vert_coords[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) {
+				fprintf(net_out, "%u %u ", tmp_instance->network->edges_to_verts[2 * j],
+								tmp_instance->network->edges_to_verts[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int j = 0; j < tmp_instance->network->num_dual_verts; j++) {
+				fprintf(net_out, "%f %f ", tmp_instance->network->dual_vert_coords[2 * j],
+								tmp_instance->network->dual_vert_coords[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) {
+				fprintf(net_out, "%u %u ", tmp_instance->network->dual_edges_to_verts[2 * j],
+								tmp_instance->network->dual_edges_to_verts[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (unsigned int j = 0; j < tmp_instance->network->num_edges; j++) {
+				fprintf(net_out, "%d ", tmp_instance->fuses[j]);
+			}
+			fclose(net_out);
+		}
+
+		free_instance(tmp_instance, &c);
+		free_instance(perm_instance, &c);
+		free_fnet(network, &c);
+
+		if (include_breaking) {
+			for (unsigned int j = 0; j < breaking_data->num_broken; j++) {
+				fprintf(break_out, "%u %g %g ", breaking_data->break_list[j],
+								breaking_data->extern_field[j], breaking_data->conductivity[j]);
+			}
+			fprintf(break_out, "\n");
+		}
+
+		free_break_data(breaking_data);
+		if (debug_stop) break;
+	}
+
+	printf("\033[F\033[JFRACTURE: COMPLETE\n");
+
+	if (save_cluster_dist) {
+		FILE *cluster_out = fopen(c_filename, "w");
+
+		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);
+
+		free(c_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_v_%c_%c_%d_%g.txt", boundc, boundc2, L,
+						 beta);
+		FILE *corr_file = fopen(corr_filename, "w");
+		for (unsigned int i = 0; i < voronoi_max_verts; i++) {
+			fprintf(corr_file, "%g ", avg_corr[i]);
+		}
+		fclose(corr_file);
+		free(corr_filename);
+		free(avg_corr);
+	}
+
+	if (save_conductivity) {
+		char *cond_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(cond_filename, filename_len, "cond_v_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
+		FILE *cond_file = fopen(cond_filename, "a");
+		for (unsigned int i = 0; i < N; i++) {
+			fprintf(cond_file, "%g ", conductivity[i]);
+		}
+		fclose(cond_file);
+		free(cond_filename);
+		free(conductivity);
+	}
+
+	if (save_toughness) {
+		char *tough_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(tough_filename, filename_len, "tuff_v_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
+		FILE *tough_file = fopen(tough_filename, "a");
+		for (unsigned int i = 0; i < N; i++) {
+			fprintf(tough_file, "%g ", toughness[i]);
+		}
+		fclose(tough_file);
+		free(tough_filename);
+		free(toughness);
+	}
+
+	if (save_damage) {
+		char *hdam_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(hdam_filename, filename_len, "damg_v_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
+		FILE *hdam_file = fopen(hdam_filename, "a");
+		for (unsigned int i = 0; i < N; i++) {
+			fprintf(hdam_file, "%g ", damage[i]);
+		}
+		fclose(hdam_file);
+		free(hdam_filename);
+		free(damage);
+	}
+
+	if (include_breaking) {
+		fclose(break_out);
+	}
+
+	if (save_crit_stress) {
+		char *a_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(a_filename, filename_len, "strs_v_%c_%c_%d_%g.txt", boundc, boundc2, L, beta);
+		FILE *a_file = fopen(a_filename, "a");
+		for (int i = 0; i < N; i++) {
+			fprintf(a_file, "%g ", crit_stress[i]);
+		}
+		fclose(a_file);
+		free(a_filename);
+		free(crit_stress);
+	}
+
+	CHOL_F(finish)(&c);
+
+	return 0;
+}
diff --git a/src/ini_network.c b/src/ini_network.c
new file mode 100644
index 0000000..0434d98
--- /dev/null
+++ b/src/ini_network.c
@@ -0,0 +1,620 @@
+
+#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];
+		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;
+}
+
+fnet *ini_square_network(unsigned int width, boundary_type boundary, bool side_bounds,
+												 cholmod_common *c) {
+	fnet *network = (fnet *)calloc(1, sizeof(fnet));
+
+	network->boundary = boundary;
+	bool periodic = (boundary == CYLINDER_BOUND) || (boundary == TORUS_BOUND) ? true : false;
+
+	network->num_edges = pow(width, 2);
+	if (boundary == CYLINDER_BOUND) {
+		assert(width % 2 == 0);
+		assert(!side_bounds);
+		network->num_verts = (width / 2) * (width + 1);
+		network->num_verts_break = (width / 2) * (width + 1);
+		network->num_dual_verts = (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->num_verts + network->num_bounds;
+	} else if (boundary == FREE_BOUND) {
+		network->num_verts = 2 * ((width + 1) / 2) * (width / 2 + 1);
+		network->num_verts_break = 2 * ((width + 1) / 2) * (width / 2 + 1);
+		network->num_dual_verts = 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->num_verts + network->num_bounds;
+	} else if (boundary == TORUS_BOUND) {
+		network->num_verts = (width / 2) * (width + 1) - (width / 2);
+		network->num_verts_break = (width / 2) * (width + 1);
+		network->num_dual_verts = (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->num_verts_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->num_verts - 1 - i;
+		} 
+	} else {
+		for (unsigned int i = 0; i < width / 2; i++) {
+			network->bound_verts[i] = i;
+		}
+	}
+	network->edges_to_verts_break =
+			(unsigned int *)calloc(2 * network->num_edges, sizeof(unsigned int));
+	network->edges_to_verts =
+			(unsigned int *)calloc(2 * network->num_edges, sizeof(unsigned int));
+	for (unsigned int i = 0; i < network->num_edges; i++) {
+		network->edges_to_verts_break[2 * i] = edge_to_verts(width, periodic, i, 1);
+		network->edges_to_verts_break[2 * i + 1] = edge_to_verts(width, periodic, i, 0);
+		network->edges_to_verts[2 * i] = network->edges_to_verts_break[2 * i] % network->num_verts;
+		network->edges_to_verts[2 * i + 1] = network->edges_to_verts_break[2 * i + 1] % network->num_verts;
+	}
+	network->verts_to_edges_ind =
+			(unsigned int *)calloc(network->num_verts + 1, sizeof(unsigned int));
+	network->verts_to_edges_ind[0] = 0;
+	unsigned int pos1 = 0;
+	for (unsigned int i = 0; i < network->num_verts; 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->verts_to_edges_ind[i + 1] = pos1;
+	}
+	network->verts_to_edges = (unsigned int *)calloc(
+			network->verts_to_edges_ind[network->num_verts], sizeof(unsigned int));
+	unsigned int *vert_counts =
+			(unsigned int *)calloc(network->num_verts, sizeof(unsigned int));
+	for (unsigned int i = 0; i < network->num_edges; i++) {
+		unsigned int v0 = network->edges_to_verts[2 * i];
+		unsigned int v1 = network->edges_to_verts[2 * i + 1];
+		network->verts_to_edges[network->verts_to_edges_ind[v0] + vert_counts[v0]] =
+				i;
+		network->verts_to_edges[network->verts_to_edges_ind[v1] + vert_counts[v1]] =
+				i;
+		vert_counts[v0]++;
+		vert_counts[v1]++;
+	}
+	free(vert_counts);
+
+	network->vert_coords =
+			(double *)malloc(2 * network->num_verts * sizeof(double));
+	for (unsigned int i = 0; i < network->num_verts; i++) {
+		if (!periodic) {
+		network->vert_coords[2 * i] = ((double)((2 * i + 1) / (width + 1)))/width;
+		network->vert_coords[2 * i + 1] = ((double)((2 * i + 1) % (width + 1)))/width;
+		} else {
+		network->vert_coords[2 * i] = ((double)((2 * i + 1) / (width)))/width;
+		network->vert_coords[2 * i + 1] = ((double)((((2 * i + 1) / (width)+1) % 2) + ((2 * i) % (width))))/width;
+		}
+	}
+
+	network->L = width;
+	network->edge_coords = get_edge_coords(
+			network->num_edges, network->vert_coords, network->edges_to_verts);
+
+	network->dual_edges_to_verts =
+			(unsigned int *)malloc(2 * network->num_edges * sizeof(unsigned int));
+	for (unsigned int i = 0; i < network->num_edges; i++) {
+		network->dual_edges_to_verts[2 * i] =
+				dual_edge_to_verts(width, periodic, i, 0) % network->num_dual_verts;
+		network->dual_edges_to_verts[2 * i + 1] =
+				dual_edge_to_verts(width, periodic, i, 1) % network->num_dual_verts;
+	}
+	network->dual_vert_coords =
+			(double *)malloc(2 * network->num_dual_verts * sizeof(double));
+	for (unsigned int i = 0; i < network->num_dual_verts; i++) {
+		network->dual_vert_coords[2 * i] =
+				2*dual_vert_to_coord(width, periodic, i, 0);
+		network->dual_vert_coords[2 * i + 1] =
+				2*dual_vert_to_coord(width, periodic, i, 1);
+	}
+
+	network->voltcurmat = gen_voltcurmat(network->num_edges,
+																			 network->break_dim,
+																			 network->edges_to_verts_break, c);
+
+	network->dual_verts_to_edges_ind =
+			get_verts_to_edges_ind(network->num_dual_verts, network->num_edges,
+														 network->dual_edges_to_verts);
+	network->dual_verts_to_edges = get_verts_to_edges(
+			network->num_dual_verts, network->num_edges, network->dual_edges_to_verts,
+			network->dual_verts_to_edges_ind);
+	network->spanning_edges = get_spanning_edges(network->num_edges, network->edges_to_verts_break, network->vert_coords, 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));
+#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 * i] = t11 < t12 ? t11 : t12;
+						dual_edges[2 * i + 1] = t11 < t12 ? t12 : t11;
+						found_match = true;
+						break;
+					}
+				}
+				if (found_match)
+					break;
+			}
+		}
+	}
+
+	return dual_edges;
+}
+
+fnet *ini_voronoi_network(unsigned int L, boundary_type boundary,
+													double *(*genfunc)(unsigned int, gsl_rng *, unsigned int *),
+													cholmod_common *c) {
+	fnet *network = (fnet *)calloc(1, sizeof(fnet));
+
+	// 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, r, &num);
+		gsl_rng_free(r);
+	}
+
+	// retrieve a periodic voronoi tesselation of the lattice
+	intptr_t *vout = run_voronoi(num, lattice, 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);
+
+	// prune the edges of the lattice and assign boundary vertices based on the
+	// desired boundary conditions
+	unsigned int num_bounds;
+	unsigned 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_t, *bound_b, *bound_l, *bound_r;
+			num_t = 0; num_b = 0; num_l = 0; num_r = 0;
+			bound_t = (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_t[v1] && !bound_l[v1] && !bound_r[v1]) {
+							bound_t[v1] = true; num_t++;
+						}
+						if (!bound_b[v2] && !bound_l[v2] && !bound_r[v2]) {
+							bound_b[v2] = true; num_b++;
+						}
+					} else {
+						if (!bound_t[v2] && !bound_l[v2] && !bound_r[v2]) {
+							bound_t[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_t[v1] && !bound_b[v1]) {
+							bound_r[v1] = true; num_r++;
+						}
+						if (!bound_l[v2] && !bound_t[v2] && !bound_b[v2]) {
+							bound_l[v2] = true; num_l++;
+						}
+					} else {
+						if (!bound_r[v2] && !bound_t[v2] && !bound_b[v2]) {
+							bound_r[v2] = true; num_r++;
+						}
+						if (!bound_l[v1] && !bound_t[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_t[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_t); free(bound_b);
+			free(tmp_edges);
+			free(tmp_dual_edges);
+			network->edges_to_verts_break = edges;
+			network->edges_to_verts = edges;
+			network->num_verts_break = num_verts;
+			network->num_verts = num_verts;
+			network->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_t, *bound_b;
+			num_t = 0; num_b = 0;
+			bound_t = (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_t[v1]) {
+							bound_t[v1] = true; num_t++;
+						}
+						if (!bound_b[v2]) {
+							bound_b[v2] = true; num_b++;
+						}
+					} else {
+						if (!bound_t[v2]) {
+							bound_t[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_t[i]) {
+					bound_verts[pos_t] = i; pos_t++;
+				} else if (bound_b[i]) {
+					bound_verts[num_t + pos_b] = i; pos_b++;
+				}
+			}
+			free(bound_t); free(bound_b);
+			free(tmp_edges);
+			free(tmp_dual_edges);
+			network->edges_to_verts_break = edges;
+			network->edges_to_verts = edges;
+			network->num_verts_break = num_verts;
+			network->num_verts = num_verts;
+			network->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_t = (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_t[v1]) {
+							bound_t[v1] = true; num_t++;
+						}
+					} else {
+						edge_change[i] = 2;
+						if (!bound_t[v2]) {
+							bound_t[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_t[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_t);
+			free(edge_change);
+			network->num_verts_break = num_verts;
+			network->num_verts = tmp_num_verts;
+			network->edges_to_verts_break = edges;
+			network->edges_to_verts = tmp_edges;
+			network->break_dim = num_verts;
+			break;
+											}
+	}
+
+	network->boundary = boundary;
+	network->num_bounds = num_bounds;
+	network->bound_inds = bound_inds;
+	network->bound_verts = bound_verts;
+	network->num_edges = num_edges;
+	network->dual_edges_to_verts = dual_edges;
+	network->vert_coords = vert_coords;
+
+	network->verts_to_edges_ind = get_verts_to_edges_ind(
+			network->num_verts, network->num_edges, network->edges_to_verts);
+	network->verts_to_edges =
+			get_verts_to_edges(network->num_verts, network->num_edges,
+												 network->edges_to_verts, network->verts_to_edges_ind);
+
+	network->L = 1;
+
+	network->edge_coords = get_edge_coords(
+			network->num_edges, network->vert_coords, network->edges_to_verts);
+
+	free(tmp_tris);
+
+	network->dual_vert_coords = lattice;
+	network->num_dual_verts = num;
+
+	network->voltcurmat = gen_voltcurmat(network->num_edges,
+																			 network->break_dim,
+																			 network->edges_to_verts_break, c);
+
+	network->dual_verts_to_edges_ind =
+			get_verts_to_edges_ind(network->num_dual_verts, network->num_edges,
+														 network->dual_edges_to_verts);
+	network->dual_verts_to_edges = get_verts_to_edges(
+			network->num_dual_verts, network->num_edges, network->dual_edges_to_verts,
+			network->dual_verts_to_edges_ind);
+
+	network->spanning_edges = get_spanning_edges(network->num_edges, network->edges_to_verts_break, network->vert_coords, 0.5, &(network->num_spanning_edges));
+
+	return network;
+}
diff --git a/src/instance.c b/src/instance.c
new file mode 100644
index 0000000..5b3d8fb
--- /dev/null
+++ b/src/instance.c
@@ -0,0 +1,118 @@
+
+#include "fracture.h"
+
+finst *create_instance(fnet *network, double inf, bool voltage_bound,
+											 bool startnow, cholmod_common *c) {
+	finst *instance = (finst *)calloc(1, sizeof(finst));
+	assert(instance != NULL);
+
+	instance->network = network;
+	instance->num_remaining_edges = network->num_edges;
+	instance->fuses = (bool *)calloc(network->num_edges, sizeof(bool));
+	assert(instance->fuses != NULL);
+	instance->inf = inf;
+	instance->voltage_bound = voltage_bound;
+	instance->boundary_cond = CHOL_F(zeros)(
+			network->break_dim, 1, CHOLMOD_REAL, c);
+	if (network->boundary != TORUS_BOUND) {
+		if (voltage_bound) {
+			for (unsigned int i = 0; i < network->bound_inds[1]; i++) {
+				((double *)instance->boundary_cond->x)[network->bound_verts[i]] =1;
+			}
+		} else {
+			((double *)instance->boundary_cond->x)[0] = 1;
+			((double *)instance->boundary_cond->x)[network->num_verts] = 1;
+			((double *)instance->boundary_cond->x)[network->num_verts + 1] = -1;
+		}
+	} else {
+		for (unsigned int i = 0; i < network->bound_inds[1]; i++) {
+			((double *)instance->boundary_cond->x)[network->bound_verts[i]] = 1;
+			((double *)instance->boundary_cond->x)[network->num_verts + i] = -1;
+		}
+		((double *)instance->boundary_cond->x)[network->bound_verts[0]] = 1;
+	}
+
+	instance->adjacency = gen_adjacency(instance, false, false, 0, c);
+
+	if (startnow) {
+		cholmod_sparse *laplacian = gen_laplacian(instance, c, true);
+		instance->factor = CHOL_F(analyze)(laplacian, c);
+		CHOL_F(factorize)(laplacian, instance->factor, c);
+		CHOL_F(free_sparse)(&laplacian, c);
+	}
+
+	instance->marks = (unsigned int *)malloc(
+			(instance->network->break_dim) *
+			sizeof(unsigned int));
+	instance->dual_marks = (unsigned int *)malloc(
+			(instance->network->num_dual_verts) *
+			sizeof(unsigned int));
+	assert(instance->marks != NULL);
+
+	for (unsigned int i = 0;
+			 i < (instance->network->break_dim);
+			 i++) {
+		instance->marks[i] = 1;
+	}
+	for (unsigned int i = 0;
+			 i < (instance->network->num_dual_verts);
+			 i++) {
+		instance->dual_marks[i] = i+1;
+	}
+	instance->num_components = 1;
+
+	return instance;
+}
+
+void finish_instance(finst *instance, cholmod_common *c) {
+	cholmod_sparse *laplacian = gen_laplacian(instance, c, true);
+	instance->factor = CHOL_F(analyze)(laplacian, c);
+	CHOL_F(factorize)(laplacian, instance->factor, c);
+	CHOL_F(free_sparse)(&laplacian, c);
+}
+
+finst *copy_instance(const finst *instance, cholmod_common *c) {
+	finst *instance_copy = (finst *)calloc(1, sizeof(finst));
+	memcpy(instance_copy, instance, sizeof(finst));
+
+	size_t fuses_size = (instance->network)->num_edges * sizeof(bool);
+	instance_copy->fuses = (bool *)malloc(fuses_size);
+	memcpy(instance_copy->fuses, instance->fuses, fuses_size);
+
+	size_t marks_size =
+			(instance->network->break_dim) *
+			sizeof(unsigned int);
+	instance_copy->marks = (unsigned int *)malloc(marks_size);
+	memcpy(instance_copy->marks, instance->marks, marks_size);
+
+	size_t dual_marks_size =
+			(instance->network->num_dual_verts) *
+			sizeof(unsigned int);
+	instance_copy->dual_marks = (unsigned int *)malloc(dual_marks_size);
+	memcpy(instance_copy->dual_marks, instance->dual_marks, dual_marks_size);
+
+	instance_copy->adjacency = CHOL_F(copy_sparse)(instance->adjacency, c);
+	instance_copy->boundary_cond = CHOL_F(copy_dense)(instance->boundary_cond, c);
+	instance_copy->factor = CHOL_F(copy_factor)(instance->factor, c);
+
+	return instance_copy;
+}
+
+void free_instance(finst *instance, cholmod_common *c) {
+	free(instance->fuses);
+	CHOL_F(free_dense)(&(instance->boundary_cond), c);
+	CHOL_F(free_sparse)(&(instance->adjacency), c);
+	CHOL_F(free_factor)(&(instance->factor), c);
+	free(instance->marks);
+	free(instance->dual_marks);
+	free(instance);
+}
+
+bool check_instance(const finst *instance, cholmod_common *c) {
+	assert(instance != NULL);
+	assert(instance->fuses != NULL);
+	assert(CHOL_F(check_dense)(instance->boundary_cond, c));
+	assert(CHOL_F(check_factor)(instance->factor, c));
+
+	return true;
+}
diff --git a/src/randfuncs.c b/src/randfuncs.c
new file mode 100644
index 0000000..260c047
--- /dev/null
+++ b/src/randfuncs.c
@@ -0,0 +1,71 @@
+
+#include "fracture.h"
+
+double *genfunc_uniform(unsigned int L, gsl_rng *r, unsigned int *num) {
+	*num = pow(L / 2 + 1, 2) + pow((L + 1) / 2, 2);
+
+	double *lattice = (double *)malloc(2 * (*num) * sizeof(double));
+	for (unsigned int i = 0; i < (*num); i++) {
+		lattice[2*i] = gsl_ran_flat(r, 0, 1);
+		lattice[2*i+1] = gsl_ran_flat(r, 0, 1);
+	}
+
+	return lattice;
+}
+
+double *genfunc_hyperuniform(unsigned int L, gsl_rng *r, unsigned int *num) {
+	*num = pow(L / 2 + 1, 2) + pow((L + 1) / 2, 2);
+
+	// necessary to prevent crashing when underflow occurs
+	gsl_set_error_handler_off();
+
+	double *lattice = (double *)malloc(2 * (*num) * sizeof(double));
+	double rho = *num;
+	for (unsigned int i = 0; i < (*num); i++) {
+		bool reject = true;
+		double x, y;
+		while(reject) {
+			x = gsl_ran_flat(r, 0, 1);
+			y = gsl_ran_flat(r, 0, 1);
+			reject = false;
+			for (unsigned int j = 0; j < i; j++) {
+				double *ds = (double *)malloc(5 * sizeof(double));
+				ds[0] = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1],2);
+				ds[1] = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1],2);
+				ds[2] = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1],2);
+				ds[3] = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] + 1,2);
+				ds[4] = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] - 1,2);
+				unsigned int min_pos = 0; double min_val = 100;
+				for (unsigned int k = 0; k < 5; k++) {
+					if (min_val > ds[k]) {
+						min_pos = k; min_val = ds[k];
+					}
+				}
+				if (1-gsl_sf_exp(-M_PI * rho * min_val) < gsl_ran_flat(r, 0, 1)) {
+					reject = true;
+					break;
+				}
+			}
+		}
+		lattice[2 * i] = x;
+		lattice[2 * i + 1] = y;
+	}
+
+	return lattice;
+}
+
+void randfunc_flat(gsl_rng *r, double *x, double *y) {
+	*x = gsl_ran_flat(r, 0, 1);
+	*y = gsl_ran_flat(r, 0, 1);
+}
+
+void randfunc_gaus(gsl_rng *r, double *x, double *y) {
+	*x = 100;
+	*y = 100;
+	double sigma = 0.25;
+	while (fabs(*x) > 0.5 || fabs(*y) > 0.5) {
+		gsl_ran_bivariate_gaussian(r, sigma, sigma, 0, x, y);
+	}
+	*x += 0.5;
+	*y += 0.5;
+}
diff --git a/src/toy_setup.c b/src/toy_setup.c
new file mode 100644
index 0000000..228acbc
--- /dev/null
+++ b/src/toy_setup.c
@@ -0,0 +1,37 @@
+
+#include "fracture.h"
+
+bool break_fuses(unsigned int width, double *fuse_thres, double *field,
+								 bool *fuses) {
+	assert(fuse_thres != NULL);
+	assert(field != NULL);
+	assert(fuses != NULL);
+	unsigned int size = pow(width, 2);
+
+	for (unsigned int i = 0; i < size; i++) {
+		if (fuses[i]) {
+			if (1 / fuse_thres[i] > field[i]) {
+				fuses[i] = 0;
+			}
+		}
+	}
+
+	return true;
+}
+
+bool gen_toy_field(unsigned int width, double strength, double *field) {
+	assert(field != NULL);
+	unsigned int size = pow(width, 2);
+
+	for (unsigned intint i = 0; i < size; i++) {
+		int x = i / width + 1;
+		int y = i % width + 1;
+		if (y < (width + 1) / 2.)
+			field[i] = fabs((width) / 2. - x) / strength;
+		else
+			field[i] =
+					sqrt(pow((width) / 2. - x, 2) + pow((width) / 2. - y, 2)) / strength;
+	}
+
+	return true;
+}
diff --git a/src/update_beta.c b/src/update_beta.c
new file mode 100644
index 0000000..4c1bf65
--- /dev/null
+++ b/src/update_beta.c
@@ -0,0 +1,34 @@
+
+#include "fracture.h"
+
+double f(double damage) {
+	assert(damage <= 1 && damage >= 0);
+	return sqrt(1 - sqrt(1 - damage));
+	//	return 0.5 - 0.68182 * (0.5 - damage);
+}
+
+double update_beta(double beta, unsigned int width, const double *stress,
+									 const double *damage, double bound_total) {
+
+	double total = 0;
+	unsigned int num_totaled = 0;
+
+	for (unsigned int i = 0; i < pow(width, 2); i++) {
+		unsigned int stress_index =
+				width / 4 + (width / 4 + (i / width) / 2) * width + (i % width) / 2;
+		double outer_damage = f(pow(fabs(stress[i]), beta));
+		double inner_stress = fabs(2 * stress[stress_index] / bound_total);
+
+		if (outer_damage > 0 && inner_stress > 0 && inner_stress != 1) {
+			total += log(outer_damage) / log(inner_stress);
+			num_totaled++;
+		}
+	}
+
+	assert(num_totaled > 0);
+	assert(total == total);
+
+	double new_beta = total / num_totaled;
+
+	return new_beta;
+}
diff --git a/src/update_boundary.c b/src/update_boundary.c
new file mode 100644
index 0000000..3d2a086
--- /dev/null
+++ b/src/update_boundary.c
@@ -0,0 +1,86 @@
+
+#include "fracture.h"
+
+void update_boundary(finst *instance, const double *avg_field) {
+	int size = instance->network->num_edges;
+	int width = sqrt(size);
+	int num_verts = instance->network->num_verts;
+	double *boundary = (double *)instance->boundary_cond->x;
+
+	boundary[num_verts] = 0;
+	boundary[num_verts + 1] = 0;
+
+	if (instance->voltage_bound) {
+		for (int i = 0; i < width / 4; i++) {
+			int t_ind = (width + 1) * (width / 8) + width / 4 - 1 + i;
+			double t_val = avg_field[t_ind];
+			boundary[2 * i] = t_val;
+			boundary[2 * i + 1] = t_val;
+
+			int b_ind = num_verts - (width + 1) * (width / 8) - width / 4 - i;
+			double b_val = avg_field[b_ind];
+			boundary[num_verts - 1 - 2 * i] = b_val;
+			boundary[num_verts - 2 - 2 * i] = b_val;
+
+			int l_ind = (width + 1) * (width / 8) + width / 2 + width / 4 - 1 +
+									i * (width + 1);
+			double l_val = avg_field[l_ind];
+			boundary[width / 2 + 2 * i * (width + 1)] = l_val;
+			boundary[width / 2 + (2 * i + 1) * (width + 1)] = l_val;
+
+			int r_ind = (width + 1) * (width / 8) + width - 1 + i * (width + 1);
+			double r_val = avg_field[r_ind];
+			boundary[width + 2 * i * (width + 1)] = r_val;
+			boundary[width + (2 * i + 1) * (width + 1)] = r_val;
+		}
+	} else {
+		const double *stress = avg_field;
+
+		for (int i = 0; i < width / 4; i++) {
+			boundary[2 * i] = stress[(width / 4 - 1) * width + width / 4 + 2 * i] +
+												stress[(width / 4 - 1) * width + width / 4 + 1 + 2 * i];
+			boundary[2 * i + 1] =
+					stress[(width / 4 - 1) * width + width / 4 + 2 * i] +
+					stress[(width / 4 - 1) * width + width / 4 + 1 + 2 * i];
+
+			boundary[num_verts - 2 * i - 1] =
+					-stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 2 * i)] -
+					stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 1 + 2 * i)];
+			boundary[num_verts - 2 * i - 2] =
+					-stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 2 * i)] -
+					stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 1 + 2 * i)];
+
+			boundary[(width + 1) / 2 + 2 * i * (width + 1)] =
+					stress[((width / 4) + 2 * i) * width + width / 4 - 1] -
+					stress[((width / 4) + 2 * i + 1) * width + width / 4 - 1];
+			boundary[(width + 1) / 2 + (2 * i + 1) * (width + 1)] =
+					stress[((width / 4) + 2 * i) * width + width / 4 - 1] -
+					stress[((width / 4) + 2 * i + 1) * width + width / 4 - 1];
+
+			boundary[(width + 1) / 2 + 2 * i * (width + 1) + width / 2] =
+					stress[((width / 4) + 2 * i) * width + width / 4 + width / 2] -
+					stress[((width / 4) + 2 * i + 1) * width + width / 4 + width / 2];
+			boundary[(width + 1) / 2 + (2 * i + 1) * (width + 1) + width / 2] =
+					stress[((width / 4) + 2 * i) * width + width / 4 + width / 2] -
+					stress[((width / 4) + 2 * i + 1) * width + width / 4 + width / 2];
+		}
+
+		double s_total = 0;
+		double total = 0;
+		for (int i = 0; i < num_verts; i++) {
+			if (boundary[i] != boundary[i]) {
+				boundary[i] = 0;
+			} else {
+				total += fabs(boundary[i]);
+				s_total += boundary[i];
+			}
+		}
+
+		// assert(fabs(s_total) < inf);
+		boundary[num_verts] -= s_total;
+
+		for (int i = 0; i < num_verts + 1; i++) {
+			boundary[i] /= total;
+		}
+	}
+}
diff --git a/src/update_factor.c b/src/update_factor.c
new file mode 100644
index 0000000..6b8e05e
--- /dev/null
+++ b/src/update_factor.c
@@ -0,0 +1,44 @@
+
+#include "fracture.h"
+
+bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2,
+									 cholmod_common *c) {
+	int n = factor->n;
+	assert(v1 < n);
+	assert(v2 < n);
+
+	cholmod_sparse *update_mat =
+			CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
+
+	unsigned int v3, v4;
+	v3 = v1 < v2 ? v1 : v2;
+	v4 = v1 > v2 ? v1 : v2;
+
+	for (int i = 0; i < n; i++) {
+		if (i <= v3)
+			((CHOL_INT *)update_mat->p)[i] = 0;
+		else
+			((CHOL_INT *)update_mat->p)[i] = 2;
+	}
+	((CHOL_INT *)update_mat->p)[n] = 2;
+	((CHOL_INT *)update_mat->i)[0] = v3;
+	((CHOL_INT *)update_mat->i)[1] = v4;
+	((double *)update_mat->x)[0] = 1;
+	((double *)update_mat->x)[1] = -1;
+
+	// assert(CHOL_F(check_sparse)(update_mat, c));
+
+	cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
+			update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
+
+	// assert(CHOL_F(check_sparse)(perm_update_mat, c));
+
+	CHOL_F(updown)(false, perm_update_mat, factor, c);
+
+	// assert(CHOL_F(check_factor)(factor, c));
+
+	CHOL_F(free_sparse)(&perm_update_mat, c);
+	CHOL_F(free_sparse)(&update_mat, c);
+
+	return true;
+}
diff --git a/src/voronoi_bound_ini.c b/src/voronoi_bound_ini.c
new file mode 100644
index 0000000..38f41cc
--- /dev/null
+++ b/src/voronoi_bound_ini.c
@@ -0,0 +1,79 @@
+
+#include "fracture.h"
+
+void voronoi_bound_ini(finst *instance, double *square_bound,
+											 unsigned int width) {
+	unsigned int square_num_verts = 2 * ((width + 1) / 2) * (width / 2 + 1);
+	double total1 = 0;
+	double total2 = 0;
+	double total3 = 0;
+	double total4 = 0;
+	for (unsigned int i = 0; i < instance->network->num_bounds; i++) {
+		for (unsigned int j = 0; j < instance->network->bound_inds[i + 1] -
+																		 instance->network->bound_inds[i];
+				 j++) {
+			double x = instance->network
+										 ->vert_coords[2 *
+																	 instance->network->bound_verts
+																			 [instance->network->bound_inds[i] + j]];
+			double y =
+					instance->network
+							->vert_coords[2 *
+																instance->network->bound_verts
+																		[instance->network->bound_inds[i] + j] +
+														1];
+
+			unsigned int vw = (width + 1) / 2;
+			unsigned int xp = (unsigned int)(x * vw);
+			unsigned int yp = (unsigned int)(y * vw);
+
+			if (i == 0) {
+				((double *)instance->boundary_cond
+						 ->x)[instance->network
+											->bound_verts[instance->network->bound_inds[i] + j]] =
+						square_bound[xp];
+				total1 += square_bound[xp];
+			}
+			if (i == 1) {
+				((double *)instance->boundary_cond
+						 ->x)[instance->network
+											->bound_verts[instance->network->bound_inds[i] + j]] =
+						square_bound[square_num_verts - vw - 1 + xp];
+				total2 += square_bound[square_num_verts - vw - 1 + xp];
+			}
+			if (i == 2) {
+				((double *)instance->boundary_cond
+						 ->x)[instance->network
+											->bound_verts[instance->network->bound_inds[i] + j]] =
+						square_bound[(width + 1) / 2 + yp * (width + 1)];
+				total3 += square_bound[(width + 1) / 2 + yp * (width + 1)];
+			}
+			if (i == 3) {
+				((double *)instance->boundary_cond
+						 ->x)[instance->network
+											->bound_verts[instance->network->bound_inds[i] + j]] =
+						square_bound[(width + 1) / 2 + yp * (width + 1) + width / 2];
+				total4 += square_bound[(width + 1) / 2 + yp * (width + 1) + width / 2];
+			}
+		}
+	}
+
+	((double *)instance->boundary_cond->x)[instance->network->num_verts] =
+			- total1 / 2 - total2 / 2;
+	((double *)instance->boundary_cond->x)[instance->network->num_verts + 1] =
+			- total1 / 2 - total2 / 2;
+	((double *)instance->boundary_cond->x)[instance->network->num_verts + 2] =
+			-total3;
+	((double *)instance->boundary_cond->x)[instance->network->num_verts + 3] =
+			-total4;
+	/*
+	((double *)instance->boundary_cond->x)[instance->network->num_verts] =
+			0;
+	((double *)instance->boundary_cond->x)[instance->network->num_verts + 1] =
+			0;
+	((double *)instance->boundary_cond->x)[instance->network->num_verts + 2] =
+			0;
+	((double *)instance->boundary_cond->x)[instance->network->num_verts + 3] =
+			0;
+			*/
+}