fixed voltage and torus conditions, current and free boundaries and broken right now

16 files changed, 647 insertions, 942 deletions

diff --git a/makefile b/makefile
index ed1142d..9ea3865 100644
--- a/makefile
+++ b/makefile
@@ -6,7 +6,7 @@ CPFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-f
 CCPFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -shared -fopenmp -march=native 
 LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler -ltcmalloc
 
-OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand
+OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand frame_create
 OBJP = spheres_poly/box spheres_poly/event spheres_poly/heap spheres_poly/neighbor spheres_poly/sphere
 OBJCCP = spheres_poly/spheres
 BIN = corr_test fracture anal_process big_anal_process cen_anal_process long_anal_process
@@ -32,7 +32,7 @@ ${OBJCCP:%=obj/%.o}: obj/%.o: src/%.cpp ${OBJP:%=obj/%.o}
 	@echo ${CCP} ${CCPFLAGS} ${OBJP:%=obj/%.o} $< -o $@
 	@${CCP} ${CCPFLAGS} ${OBJP:%=obj/%.o} -fuse-ld=gold $< -o $@
 
-bin/%: obj/%.o ${OBJ:%=obj/%.o} ${OBJCCP:%=obj/%.o}
+bin/%: obj/%.o ${OBJ:%=obj/%.o} ${OBJCCP:%=obj/%.o} src/fracture.h
 	@echo ${CC} ${OBJ:%=obj/%.o} ${OBJCCP:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} $< -o $@
 	@${CC} ${OBJ:%=obj/%.o} ${OBJCCP:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} $< -o $@
 
diff --git a/src/bound_set.c b/src/bound_set.c
index f513aa1..65f1e6f 100644
--- a/src/bound_set.c
+++ b/src/bound_set.c
@@ -26,35 +26,73 @@ void bound_set_embedded(double *bound, const graph_t *g, double notch_len) {
 		y3 = (2. * i + 1.) / L - 0.5; x3 = 0.5 - 1.;
 		y4 = (2. * i + 1.) / L - 0.5; x4 = 0.5 - 0.;
 
-		bound[g->bound_verts[g->bound_inds[0] + i]] = u_y(x1, y1);
-		bound[g->bound_verts[g->bound_inds[1] + i]] = u_y(x2, y2);
-		bound[g->bound_verts[g->bound_inds[2] + i]] = u_y(x3, y3);
-		bound[g->bound_verts[g->bound_inds[3] + i]] = u_y(x4, y4);
+		bound[g->b[g->bi[0] + i]] = u_y(x1, y1);
+		bound[g->b[g->bi[1] + i]] = u_y(x2, y2);
+		bound[g->b[g->bi[2] + i]] = u_y(x3, y3);
+		bound[g->b[g->bi[3] + i]] = u_y(x4, y4);
 	}
 }
 
+bool is_in(uint_t len, uint_t *list, uint_t element) {
+	for (uint_t i = 0; i < len; i++) {
+		if (list[i] == element) { 
+			return true;
+		}
+	}
+	return false;
+}
+
 cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_common *c) {
-	cholmod_dense *boundary = CHOL_F(zeros)(g->break_dim, 1, CHOLMOD_REAL, c);
+
+	uint_t dim = g->nv;
+
+	if (vb && g->boundary != TORUS_BOUND) {
+		dim -= g->bi[g->nb];
+	} else if (!vb) {
+		dim += 2;
+	}
+
+	cholmod_dense *boundary = CHOL_F(zeros)(dim, 1, CHOLMOD_REAL, c);
 	double *bound = (double *)boundary->x;
 
 	switch (g->boundary) {
 		case TORUS_BOUND:
-			for (uint_t i = 0; i < g->bound_inds[1]; i++) {
-				bound[g->bound_verts[i]] = 1;
-				bound[g->nv + i] = -1;
+			for (uint_t i = 0; i < g->bi[1]; i++) {
+				uint_t be = g->b[i];
+				uint_t v1 = g->ev[2 * be];
+				uint_t v2 = g->ev[2 * be + 1];
+				double v1y = g->vx[2 * v1 + 1];
+				double v2y = g->vx[2 * v2 + 1];
+
+				uint_t ind = v1y < v2y ? 0 : 1;
+
+				bound[g->ev[2 * be + ind]] += 1;
+				bound[g->ev[2 * be + !ind]] -= 1;
 			}
-			bound[g->bound_verts[0]] = 1;
 			break;
+			/*
 		case EMBEDDED_BOUND:
 			bound_set_embedded(bound, g, notch_len);
 			break;
+			*/
 		default:
 			if (vb) {
-				for (uint_t i = 0; i < g->bound_inds[1]; i++) {
-					bound[g->bound_verts[i]] = 1;
+				for (uint_t i = 0; i < dim; i++) {
+					uint_t v = g->nbi[i];
+					for (uint_t j = 0; j < g->vei[v+1] - g->vei[v]; j++) {
+						uint_t e = g->ve[g->vei[v] + j];
+						uint_t v0 = g->ev[2 * e];
+						uint_t v1 = g->ev[2 * e + 1];
+
+						if (g->bq[v0] || g->bq[v1]) {
+							uint_t vv = v0 == v ? v1 : v0;
+							if (is_in(g->bi[1], g->b, vv)) {
+								bound[i]++;
+							}
+						}
+					}
 				}
 			} else {
-				bound[0] = 1;
 				bound[g->nv] = 1;
 				bound[g->nv + 1] = -1;
 			}
diff --git a/src/break_edge.c b/src/break_edge.c
index 4e8bdcf..2f112c2 100644
--- a/src/break_edge.c
+++ b/src/break_edge.c
@@ -1,104 +1,51 @@
 
 #include "fracture.h"
 
-bool break_edge(net_t *instance, uint_t edge, cholmod_common *c) {
-	instance->fuses[edge] = true;
-	instance->num_broken++;
 
-	if (instance->factor != NULL) {
-		uint_t w1 = instance->graph->ev_break[2 * edge];
-		uint_t w2 = instance->graph->ev_break[2 * edge + 1];
-		factor_update(instance->factor, w1, w2, c);
-	}
-
-	uint_t v1, v2, s1, s2, dv1, dv2, ds1, ds2;
+bool break_edge(net_t *net, uint_t edge, cholmod_common *c) {
+	net->fuses[edge] = true;
+	net->num_broken++;
 
-	v1 = instance->graph->ev[2 * edge];
-	v2 = instance->graph->ev[2 * edge + 1];
-	dv1 = instance->graph->dev[2 * edge];
-	dv2 = instance->graph->dev[2 * edge + 1];
+	double *x = (double *)net->boundary_cond->x;
+	const graph_t *g = net->graph;
 
-	s1 = v1 > v2 ? v1 : v2;
-	s2 = v1 > v2 ? v2 : v1;
-	ds1 = dv1 > dv2 ? dv1 : dv2;
-	ds2 = dv1 > dv2 ? dv2 : dv1;
+	uint_t v1 = net->graph->ev[2 * edge];
+	uint_t v2 = net->graph->ev[2 * edge + 1];
 
-	{
-		int_t *lap_p = (int_t *)instance->adjacency->p;
-		int_t *lap_i = (int_t *)instance->adjacency->i;
-		double *lap_x = (double *)instance->adjacency->x;
+	uint_t s1 = v1 < v2 ? v1 : v2;
+	uint_t s2 = v1 < v2 ? v2 : v1;
 
-		for (int i = 0; i < lap_p[s1 + 1] - lap_p[s1]; i++) {
-			if (lap_i[lap_p[s1] + i] == s2)
-				lap_x[lap_p[s1] + i] = 0;
+	if (net->graph->boundary == TORUS_BOUND) {
+		if (net->graph->bq[edge]) {
+			double v1y = g->vx[2 * v1 + 1];
+			double v2y = g->vx[2 * v2 + 1];
+			uint_t ind = v1y < v2y ? 0 : 1;
+			x[g->ev[2 * edge + ind]] -= 1;
+			x[g->ev[2 * edge + !ind]] += 1;
 		}
-		for (int i = 0; i < lap_p[s2 + 1] - lap_p[s2]; i++) {
-			if (lap_i[lap_p[s2] + i] == s1)
-				lap_x[lap_p[s2] + i] = 0;
+		if (net->factor != NULL) {
+			factor_update(net->factor, s1, s2, c);
 		}
-	}
-
-	int_t old_num_components = instance->num_components;
-
-	instance->num_components = update_components(
-			instance->adjacency, instance->marks, old_num_components, s1, s2, 0);
-
-
-	if (instance->graph->boundary == TORUS_BOUND) {
-		if (instance->dual_marks[dv1] == instance->dual_marks[dv2]) {
-			int **cycles = (int **)malloc(4*instance->graph->ne * sizeof(int *));
-			unsigned int num_cycles = find_cycles(instance->graph->ne, instance->fuses, instance->graph->dev, instance->graph->dvei, instance->graph->dve, 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->graph->dev[2 * ee + !side];
-					v2 = instance->graph->dev[2 * ee + side];
-					v1x = instance->graph->dvx[2 * v1];
-					v1y = instance->graph->dvx[2 * v1 + 1];
-					v2x = instance->graph->dvx[2 * v2];
-					v2y = instance->graph->dvx[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]);
+	} else if (net->voltage_bound) {
+		if (g->bq[v1] || g->bq[v2]) {
+			uint_t vv = g->bq[v1] ? v2 : v1;
+			uint_t uu = v1 == vv ? v2 : v1;
+			if (is_in(g->bi[1], g->b, uu)) {
+				x[g->bni[vv]] -= 1;
+			}
+			if (net->factor != NULL) {
+				factor_update2(net->factor, g->bni[vv], c);
+			}
+		} else {
+			if (net->factor != NULL) {
+				factor_update(net->factor, g->bni[s1], g->bni[s2], c);
 			}
-			free(cycles);
-		}
-	}
-
-	{
-		int_t *lap_p = (int_t *)instance->dual_adjacency->p;
-		int_t *lap_i = (int_t *)instance->dual_adjacency->i;
-		double *lap_x = (double *)instance->dual_adjacency->x;
-
-		for (int i = 0; i < lap_p[ds1 + 1] - lap_p[ds1]; i++) {
-			if (lap_i[lap_p[ds1] + i] == ds2)
-				lap_x[lap_p[ds1] + i] = 1;
 		}
-		for (int i = 0; i < lap_p[ds2 + 1] - lap_p[ds2]; i++) {
-			if (lap_i[lap_p[ds2] + i] == ds1)
-				lap_x[lap_p[ds2] + i] = 1;
+	} else {
+		if (net->factor != NULL) {
+			factor_update(net->factor, s1, s2, c);
 		}
 	}
 
-	set_connected(instance->dual_adjacency, instance->dual_marks, dv1, instance->dual_marks[dv2], -1, 0);
-
 	return true;
 }
diff --git a/src/factor_update.c b/src/factor_update.c
index f27971b..a51705a 100644
--- a/src/factor_update.c
+++ b/src/factor_update.c
@@ -36,3 +36,33 @@ void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common
 	CHOL_F(free_sparse)(&perm_update_mat, c);
 	CHOL_F(free_sparse)(&update_mat, c);
 }
+
+void factor_update2(cholmod_factor *factor, uint_t v, cholmod_common *c) {
+	uint_t n = factor->n;
+
+	cholmod_sparse *update_mat =
+			CHOL_F(allocate_sparse)(n, n, 1, true, true, 0, CHOLMOD_REAL, c);
+
+	int_t *pp = (int_t *)update_mat->p;
+	int_t *ii = (int_t *)update_mat->i;
+	double *xx = (double *)update_mat->x;
+
+	for (uint_t i = 0; i <= v; i++) {
+		pp[i] = 0;
+	}
+
+	for (uint_t i = v + 1; i <= n; i++) {
+		pp[i] = 1;
+	}
+
+	ii[0] = v;
+	xx[0] = 1;
+
+	cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
+			update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
+
+	CHOL_F(updown)(false, perm_update_mat, factor, c);
+
+	CHOL_F(free_sparse)(&perm_update_mat, c);
+	CHOL_F(free_sparse)(&update_mat, c);
+}
diff --git a/src/fracture.c b/src/fracture.c
index 1ed594f..6ad2f26 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -263,13 +263,6 @@ int main(int argc, char *argv[]) {
 
 		graph_t *g = graph_create(lattice, boundary, L, use_dual, &c);
 		net_t *net = net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c);
-		if (net->marks[g->nv] != net->marks[g->nv+1]) {
-			if (save_crit_stress) crit_stress[i] = 0;
-			if (save_conductivity) conductivity[i] = 0;
-			if (save_damage) damage[net->num_broken]++;
-			if (save_energy) energy[i] = 0;
-			if (save_threshold) thresholds[i] = 0;
-		} else {
 		net_t *tmp_net = net_copy(net, &c);
 		data_t *data = net_fracture(tmp_net, &c, cutoff);
 		net_free(tmp_net, &c);
@@ -325,14 +318,6 @@ int main(int argc, char *argv[]) {
 			}
 		}
 
-		FILE *testout = fopen("test.txt", "w");
-		double *tmp_voltage = net_voltages(net, &c);
-		for (uint_t j = 0; j < g->nv; j++) {
-			fprintf(testout, "%g ", tmp_voltage[j]);
-		}
-		fclose(testout);
-		free(tmp_voltage);
-
 		if (save_damage) {
 			uint_t would_break = 0;
 			double *tmp_voltage = net_voltages(net, &c);
@@ -382,7 +367,6 @@ int main(int argc, char *argv[]) {
 		}
 
 		data_free(data);
-		}
 		if (save_network) {
 			FILE *net_out = fopen("network.txt", "w");
 			for (uint_t j = 0; j < g->nv; j++) {
diff --git a/src/fracture.h b/src/fracture.h
index 8f883bc..0a3a687 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -48,26 +48,35 @@ typedef struct {
 	uint_t ne;
 	uint_t nv;
 	uint_t dnv;
-	uint_t nv_break;
-	uint_t num_bounds;
+	uint_t *ev;
+	uint_t *dev;
+	double *vx;
+	double *dvx;
+} frame_t;
+
+typedef struct {
+	uint_t L;
+	bound_t boundary;
+	uint_t ne;
+	uint_t nv;
+	uint_t dnv;
+	uint_t nb;
 	uint_t *ev; 
-	uint_t *ev_break;
+	uint_t *dev;
+	double *vx;
+	double *dvx;
 	uint_t *vei;
 	uint_t *ve;
-	uint_t *bound_inds;
-	uint_t *bound_verts;
-	double *vx;
-	double *ex;
-	uint_t *dev;
 	uint_t *dvei;
 	uint_t *dve;
-	double *dvx;
-	uint_t num_spanning_edges;
+	uint_t *bi;
+	uint_t *b;
+	uint_t *nbi;
+	uint_t *bni;
+	bool *bq;
 	uint_t *spanning_edges;
-	uint_t L;
-	uint_t break_dim;
+	uint_t num_spanning_edges;
 	cholmod_sparse *voltcurmat;
-	bound_t boundary;
 } graph_t;
 
 typedef struct {
@@ -77,14 +86,11 @@ typedef struct {
 	double inf;
 	cholmod_dense *boundary_cond;
 	cholmod_factor *factor;
-	uint_t *marks;
-	uint_t *dual_marks;
 	bool voltage_bound;
-	uint_t num_components;
-	cholmod_sparse *adjacency;
-	cholmod_sparse *dual_adjacency;
-	bool debug_stop;
 	uint_t num_broken;
+	uint_t dim;
+	uint_t nep;
+	uint_t *evp;
 } net_t;
 
 typedef struct {
@@ -101,11 +107,9 @@ int update_components(const cholmod_sparse *laplacian, uint_t *marks,
 
 uint_t *find_components(const cholmod_sparse *laplacian, uint_t skip);
 
-cholmod_sparse *gen_adjacency(const net_t *instance, bool dual, bool breakv,
-															uint_t pad, cholmod_common *c);
+cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c);
 
-cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
-															bool symmetric);
+cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c);
 
 int edge_to_verts(uint_t width, bool periodic, uint_t edge,
 									bool index);
@@ -160,8 +164,8 @@ uint_t *get_clusters(net_t *instance, cholmod_common *c);
 
 uint_t *get_cluster_dist(net_t *instance, cholmod_common *c);
 
-double *genfunc_uniform(uint_t L, bound_t boundary, gsl_rng *r, uint_t *num);
-double *genfunc_hyperuniform(uint_t L, bound_t boundary, gsl_rng *r, uint_t *num);
+double *genfunc_uniform(uint_t num, gsl_rng *r);
+double *genfunc_hyperuniform(uint_t num, gsl_rng *r);
 void randfunc_flat(gsl_rng *r, double *x, double *y);
 void randfunc_gaus(gsl_rng *r, double *x, double *y);
 
@@ -188,3 +192,7 @@ unsigned long int rand_seed();
 long double rand_dist_pow(const gsl_rng *r, double beta);
 
 double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate);
+
+frame_t *frame_create(lattice_t lattice, uint_t L, bool dual);
+void frame_free(frame_t *frame);
+bool is_in(uint_t len, uint_t *list, uint_t element);
diff --git a/src/frame_create.c b/src/frame_create.c
new file mode 100644
index 0000000..89ce69d
--- /dev/null
+++ b/src/frame_create.c
@@ -0,0 +1,167 @@
+
+#include "fracture.h"
+
+uint_t *get_voro_dual_edges(uint_t num_edges,
+																	uint_t num_verts, uint_t *edges,
+																	uint_t *triangles) {
+	uint_t *dual_edges =
+			(uint_t *)malloc(2 * num_edges * sizeof(uint_t));
+	uint_t place = 0;
+	for (uint_t i = 0; i < num_edges; i++) {
+		uint_t v1, v2;
+		v1 = edges[2 * i];
+		v2 = edges[2 * i + 1];
+		if (v1 < num_verts && v2 < num_verts) {
+			bool found_match = false;
+			for (uint_t j = 0; j < 3; j++) {
+				for (uint_t k = 0; k < 3; k++) {
+					uint_t t11, t12, t21, t22;
+					t11 = triangles[3 * v1 + j];
+					t12 = triangles[3 * v1 + ((j + 1) % 3)];
+					t21 = triangles[3 * v2 + k];
+					t22 = triangles[3 * v2 + ((k + 1) % 3)];
+					if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) {
+						dual_edges[2 * place] = t11 < t12 ? t11 : t12;
+						dual_edges[2 * place + 1] = t11 < t12 ? t12 : t11;
+						place++;
+						found_match = true;
+						break;
+					}
+				}
+				if (found_match)
+					break;
+			}
+		}
+	}
+
+	return dual_edges;
+}
+
+frame_t *frame_create_voronoi(uint_t L, bool dual, bool hyperuniform) {
+	double *dvx = NULL;
+	uint_t dnv = 2 * pow(L / 2, 2);
+
+	{
+		gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+		gsl_rng_set(r, rand_seed());
+
+		if (hyperuniform) {
+			dvx = genfunc_hyperuniform(dnv, r);
+		} else {
+			dvx = genfunc_uniform(dnv, r);
+		}
+
+		gsl_rng_free(r);
+	}
+
+	// retrieve a periodic voronoi tesselation of the lattice
+	intptr_t *vout = run_voronoi(dnv, dvx, true, 0, 1, 0, 1);
+
+	uint_t nv = ((uint_t *)vout[0])[0];
+	uint_t ne = ((uint_t *)vout[0])[1];
+	double *vx = (double *)vout[2];
+	uint_t *ev = (uint_t *)vout[3];
+	uint_t *voro_tris = (uint_t *)vout[5];
+
+	free((void *)vout[0]);
+	free((void *)vout[1]);
+	free((void *)vout[4]);
+	free(vout);
+
+	// get dual edges of the fully periodic graph
+	uint_t *dev = get_voro_dual_edges(ne, nv, ev, voro_tris);
+
+	frame_t *frame = (frame_t *)malloc(sizeof(frame_t));
+	frame->ne = ne;
+
+	// when use_dual is specificed, the edge and vertex sets are swapped with the
+	// dual edge and dual vertex sets.  once formally relabelled, everything
+	// works the same way
+	if (dual) {
+		frame->nv = dnv;
+		frame->dnv = nv;
+		frame->ev = dev;
+		frame->dev = ev;
+		frame->vx = dvx;
+		frame->dvx = vx;
+	} else {
+		frame->nv = nv;
+		frame->dnv = dnv;
+		frame->ev = ev;
+		frame->dev = dev;
+		frame->vx = vx;
+		frame->dvx = dvx;
+	}
+
+	return frame;
+}
+
+frame_t *frame_create_square(uint_t L, bool dual) {
+	uint_t nv = 2 * pow(L / 2, 2);
+	uint_t ne = pow(L, 2);
+
+	uint_t *ev = (uint_t *)malloc(2 * ne * sizeof(uint_t));
+	uint_t *dev = (uint_t *)malloc(2 * ne * sizeof(uint_t));
+	double *vx = (double *)malloc(2 * nv * sizeof(double));
+	double *dvx = (double *)malloc(2 * nv * sizeof(double));
+
+	for (uint_t i = 0; i < ne; i++) {
+		uint_t x = i / L;
+		uint_t y = i % L;
+
+		ev[2 * i] = (L * x) / 2 + ((y + x % 2) / 2) % (L / 2);
+		ev[2 * i + 1] = ((L * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2)) % nv;
+
+		dev[2 * i] = (L * x) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2);
+		dev[2 * i + 1] = ((L * (x + 1)) / 2 + ((y + x % 2) / 2) % (L / 2)) % nv;
+	}
+
+	double dx = 1. / L;
+
+	for (uint_t i = 0; i < nv; i++) {
+		vx[2 * i] = dx * ((1 + i / (L / 2)) % 2 + 2 * (i % (L / 2)));
+		vx[2 * i + 1] = dx * (i / (L / 2));
+
+		dvx[2 * i] = dx * ((i / (L / 2)) % 2 + 2 * (i % (L / 2)));
+		dvx[2 * i + 1] = dx * (i / (L / 2));
+	}
+
+	frame_t *frame = (frame_t *)malloc(sizeof(frame_t));
+	frame->ne = ne;
+	frame->nv = nv;
+	frame->dnv = nv;
+
+	if (dual) {
+		frame->ev = dev;
+		frame->dev = ev;
+		frame->vx = dvx;
+		frame->dvx = vx;
+	} else {
+		frame->ev = ev;
+		frame->dev = dev;
+		frame->vx = vx;
+		frame->dvx = dvx;
+	}
+
+	return frame;
+}
+
+frame_t *frame_create(lattice_t lattice, uint_t L, bool dual) {
+	switch (lattice) {
+		case (SQUARE_LATTICE):
+			return frame_create_square(L, dual);
+		case (VORONOI_LATTICE):
+			return frame_create_voronoi(L, dual, false);
+		case (VORONOI_HYPERUNIFORM_LATTICE):
+			return frame_create_voronoi(L, dual, true);
+	}
+}
+
+void frame_free(frame_t *frame) {
+	free(frame->ev);
+	free(frame->dev);
+	free(frame->vx);
+	free(frame->dvx);
+	free(frame);
+}
+
diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c
index 411734c..1cc93a4 100644
--- a/src/gen_laplacian.c
+++ b/src/gen_laplacian.c
@@ -1,118 +1,79 @@
 
 #include "fracture.h"
 
-cholmod_sparse *gen_adjacency(const net_t *instance, bool dual, bool breakv,
-															uint_t pad, cholmod_common *c) {
-	uint_t ne = instance->graph->ne;
-	if (!breakv && instance->graph->boundary != TORUS_BOUND && !dual) {
-		ne += instance->graph->bound_inds[2];
+cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c) {
+	const graph_t *g = net->graph;
+	uint_t nv;
+	uint_t ne;
+	uint_t nre;
+	uint_t *ev;
+
+	if (use_gp) {
+		nv = net->dim;
+		ne = net->nep;
+		nre = (int_t)net->nep - (int_t)net->num_broken;
+		ev = net->evp;
+	} else if (dual) {
+		nv = g->dnv;
+		ne = g->ne;
+		nre = net->num_broken;
+		ev = g->dev;
+	} else {
+		nv = g->nv;
+		ne = g->ne;
+		nre = (int_t)g->ne - (int_t)net->num_broken;
+		ev = g->ev;
 	}
 
-	uint_t nnz = 2 * ne;
+	uint_t nnz = nre;
 
-	uint_t nv;
-	if (dual)
-		nv = instance->graph->dnv;
-	else {
-		if (breakv) nv = instance->graph->nv_break;
-		else nv = instance->graph->nv;
-		if (instance->graph->boundary != TORUS_BOUND && !breakv) {
-			nv += 2;
-		}
-	}
+	cholmod_triplet *t = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
 
-	cholmod_triplet *t =
-			CHOL_F(allocate_triplet)(nv + pad, nv + pad, nnz, 0, CHOLMOD_REAL, c);
 	int_t *ri = (int_t *)t->i;
 	int_t *ci = (int_t *)t->j;
 	double *ai = (double *)t->x;
 
 	t->nnz = nnz;
 
-	uint_t *etv;
-	if (dual)
-		etv = instance->graph->dev;
-	else {
-		if (breakv) etv = instance->graph->ev_break;
-		else etv = instance->graph->ev;
-	}
-
-	uint_t count = 0;
-
-	for (uint_t i = 0; i < instance->graph->ne; i++) {
-		if ((instance->fuses[i] && dual) || (!instance->fuses[i] && !dual)) {
-			uint_t v1 = etv[2 * i];
-			uint_t v2 = etv[2 * i + 1];
+	uint_t a = 0;
 
-			ri[2 * count] = v1;
-			ci[2 * count] = v2;
-			ai[2 * count] = 1;
+	for (uint_t i = 0; i < ne; i++) {
+		if ((net->fuses[i] && dual) || (!net->fuses[i] && !dual)) {
+			uint_t v1 = ev[2 * i];
+			uint_t v2 = ev[2 * i + 1];
 
-			ri[2 * count + 1] = v2;
-			ci[2 * count + 1] = v1;
-			ai[2 * count + 1] = 1;
+			uint_t s1 = v1 < v2 ? v1 : v2;
+			uint_t s2 = v1 < v2 ? v2 : v1;
 
-			count++;
-		} else {
-			uint_t v1 = etv[2 * i];
-			uint_t v2 = etv[2 * i + 1];
+			ri[a] = s2;
+			ci[a] = s1;
+			ai[a] = 1;
 
-			ri[2 * count] = v1;
-			ci[2 * count] = v2;
-			ai[2 * count] = 0;
-
-			ri[2 * count + 1] = v2;
-			ci[2 * count + 1] = v1;
-			ai[2 * count + 1] = 0;
-
-			count++;
+			a++;
 		}
 	}
 
-	if (!breakv && instance->graph->boundary != TORUS_BOUND && !dual) {
-		for (uint_t i = 0; i < 2; i++) {
-			for (uint_t j = 0; j < instance->graph->bound_inds[i+1] - instance->graph->bound_inds[i]; j++) {
-				ri[2*count] = instance->graph->nv + i;
-				ci[2*count] = instance->graph->bound_verts[instance->graph->bound_inds[i] + j];
-				ai[2*count] = 1;
-				ri[2*count+1] = instance->graph->bound_verts[instance->graph->bound_inds[i] + j];
-				ci[2*count+1] = instance->graph->nv + i;
-				ai[2*count+1] = 1;
-				count++;
-			}
-		}
-	}
-
-
 	cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c);
 	CHOL_F(free_triplet)(&t, c);
 
+	if (!symmetric) {
+		cholmod_sparse *tmp_s = CHOL_F(copy)(s, 0, 1, c);
+		CHOL_F(free_sparse)(&s, c);
+		s = tmp_s;
+	}
+
 	return s;
 }
 
-cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
-															bool symmetric) {
-	const graph_t *network = instance->graph;
-	uint_t num_verts = network->nv_break;
-	uint_t num_bounds = network->num_bounds;
-	double inf = instance->inf;
-	bool voltage_bound = instance->voltage_bound;
-	uint_t *bound_inds = network->bound_inds;
-	uint_t *bound_verts = network->bound_verts;
-
-	uint_t num_gedges;
-	if (voltage_bound) {
-		num_gedges = network->bound_inds[num_bounds];
-	} else {
-		num_gedges = network->bound_inds[num_bounds] + (num_bounds - 2) * 2;
-	}
-	if (network->boundary == TORUS_BOUND) num_gedges = bound_inds[1];
-	uint_t num_gverts = network->break_dim;
+cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c) {
+	const graph_t *g = net->graph;
+	uint_t nv = net->dim;
+	uint_t ne = net->nep;
+	uint_t *ev = net->evp;
 
-	uint_t nnz = num_gverts + 2 * num_gedges;
+	uint_t nnz = nv;
 
-	cholmod_triplet *temp_m =
-			CHOL_F(allocate_triplet)(num_gverts, num_gverts, nnz, 0, CHOLMOD_REAL, c);
+	cholmod_triplet *temp_m = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
 
 	int_t *rowind = (int_t *)temp_m->i;
 	int_t *colind = (int_t *)temp_m->j;
@@ -120,126 +81,56 @@ cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
 
 	temp_m->nnz = nnz;
 
-	for (uint_t i = 0; i < num_gverts; i++) {
+	for (uint_t i = 0; i < nv; 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);
-	int_t *ia = (int_t *)adjacency->p;
-	double *aa = (double *)adjacency->x;
+	cholmod_sparse *adjacency = gen_adjacency(net, false, true, true, c);
 
-	for (uint_t i = 0; i < num_verts; i++) {
-		for (uint_t j = 0; j < ia[i + 1] - ia[i]; j++) {
-			acoo[i] += aa[ia[i] + j];
-		}
-	}
+	for (uint_t i = 0; i < ne; i++) {
+		if (!net->fuses[i]){
+			uint_t v1 = ev[2 * i];
+			uint_t v2 = ev[2 * i + 1];
 
-	if (network->boundary == TORUS_BOUND) {
-		for (uint_t i = 0; i < network->bound_inds[1]; i++) {
-			uint_t vv = network->bound_verts[i];
-			rowind[num_gverts + 2*i] = vv;
-			colind[num_gverts + 2*i] = network->nv + i;
-			acoo[num_gverts + 2*i] = -1;
-			rowind[num_gverts + 2*i+1] = network->nv + i;
-			colind[num_gverts + 2*i+1] = vv;
-			acoo[num_gverts + 2*i+1] = -1;
-			acoo[vv]++;
-			acoo[network->nv + i]++;
+			acoo[v1]++;
+			acoo[v2]++;
 		}
-		acoo[bound_verts[0]]++;
 	}
 
-	if (network->boundary != TORUS_BOUND) {
-		if (voltage_bound) {
-			for (uint_t i = 0; i < num_bounds; i++) {
-				for (uint_t j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) {
-					acoo[bound_verts[bound_inds[i] + j]]++;
-
-					rowind[nnz - 1 - 2 * (bound_inds[i] + j)] = bound_verts[bound_inds[i] + j];
-					colind[nnz - 1 - 2 * (bound_inds[i] + j)] = num_gverts - 1;
-					acoo[nnz - 1 - 2 * (bound_inds[i] + j)] = -1;
-					rowind[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = num_gverts - 1;
-					colind[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = bound_verts[bound_inds[i] + j];
-					acoo[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = -1;
-					acoo[num_gverts - 1]++;
-					acoo[bound_verts[bound_inds[i] + j]]++;
-				}
-			}
-		} else {
-			for (uint_t i = 0; i < num_bounds; i++) {
-				if (network->boundary != EMBEDDED_BOUND) acoo[0]++;
-
-				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;
-			}
+	if (net->voltage_bound && g->boundary != TORUS_BOUND) {
+		for (uint_t i = 0; i < net->dim; i++) {
+			uint_t v = g->nbi[i];
+			for (uint_t j = 0; j < g->vei[v+1] - g->vei[v]; j++) {
+				uint_t e = g->ve[g->vei[v] + j];
+				uint_t v0 = g->ev[2 * e];
+				uint_t v1 = g->ev[2 * e + 1];
 
-			uint_t start = num_gverts;
-
-			for (uint_t i = 0; i < 2; i++) {
-				for (uint_t 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;
+				if (g->bq[v0] || g->bq[v1]) {
+					acoo[i]++;
 				}
 			}
-
-			for (uint_t 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;
-			}
 		}
+	} else {
+		acoo[0]++;
 	}
 
-
-	for (uint_t i = 0; i < num_gverts; i++) {
-		if (instance->marks[i] != instance->marks[network->nv])
-			acoo[i]++;
+	for (uint_t i = 0; i < nv; i++) {
+		if (acoo[i] == 0) acoo[i]++;
 	}
 
-	assert(CHOL_F(check_triplet)(temp_m, c));
+	//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));
+	//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);
+	cholmod_sparse *laplacian = CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 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/get_dual_clusters.c b/src/get_dual_clusters.c
index 3a51c38..78bf185 100644
--- a/src/get_dual_clusters.c
+++ b/src/get_dual_clusters.c
@@ -2,7 +2,7 @@
 #include "fracture.h"
 
 unsigned int *get_clusters(net_t *instance, cholmod_common *c) {
-	cholmod_sparse *s_dual = gen_adjacency(instance, true, false, 0, c);
+	cholmod_sparse *s_dual = gen_adjacency(instance, true, false, true, c);
 
 	unsigned int *dual_marks = find_components(s_dual, 0);
 	CHOL_F(free_sparse)(&s_dual, c);
diff --git a/src/graph_create.c b/src/graph_create.c
index c1f556c..09a3aed 100644
--- a/src/graph_create.c
+++ b/src/graph_create.c
@@ -97,586 +97,203 @@ uint_t *get_verts_to_edges(uint_t num_verts, uint_t num_edges,
 	return output;
 }
 
-graph_t *ini_square_network(uint_t width, bound_t boundary, bool side_bounds,
-												 cholmod_common *c) {
-	graph_t *network = (graph_t *)calloc(1, sizeof(graph_t));
-
-	network->boundary = boundary;
-	bool periodic = (boundary == CYLINDER_BOUND) || (boundary == TORUS_BOUND) ? true : false;
-
-	network->ne = pow(width, 2);
-	if (boundary == CYLINDER_BOUND) {
-		assert(width % 2 == 0);
-		assert(!side_bounds);
-		network->nv = (width / 2) * (width + 1);
-		network->nv_break = (width / 2) * (width + 1);
-		network->dnv = (width / 2 + 1) * (width / 2) + pow(width / 2, 2);
-		network->num_bounds = 2;
-		network->bound_inds = (uint_t *)malloc((network->num_bounds + 1) *
-																								 sizeof(uint_t));
-		network->bound_inds[0] = 0;
-		network->bound_inds[1] = width / 2;
-		network->bound_inds[2] = width;
-		network->bound_verts = (uint_t *)calloc(width, sizeof(uint_t));
-		network->break_dim = network->nv + network->num_bounds;
-	} else if (boundary == FREE_BOUND || boundary == EMBEDDED_BOUND) {
-		network->nv = 2 * ((width + 1) / 2) * (width / 2 + 1);
-		network->nv_break = 2 * ((width + 1) / 2) * (width / 2 + 1);
-		network->dnv = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2);
-		if (side_bounds)
-			network->num_bounds = 4;
-		else
-			network->num_bounds = 2;
-		network->bound_inds = (uint_t *)malloc((network->num_bounds + 1) *
-																								 sizeof(uint_t));
-		for (uint_t i = 0; i < network->num_bounds + 1; i++) {
-			network->bound_inds[i] = i * ((width + 1) / 2);
-		}
-		network->bound_verts = (uint_t *)malloc(
-				network->num_bounds * ((width + 1) / 2) * sizeof(uint_t));
-		if (side_bounds) {
-			for (uint_t i = 0; i < (width + 1) / 2; i++) {
-				network->bound_verts[2 * ((width + 1) / 2) + i] =
-						(width + 1) / 2 + i * (width + 1);
-				if (width % 2) {
-					network->bound_verts[3 * ((width + 1) / 2) + i] =
-							(width + 1) / 2 + i * (width + 1) - 1;
-				} else {
-					network->bound_verts[3 * ((width + 1) / 2) + i] =
-							(width + 1) / 2 + i * (width + 1) + width / 2;
-				}
-			}
-		}
-		network->break_dim = network->nv + 2;
-	} else if (boundary == TORUS_BOUND) {
-		network->nv = (width / 2) * (width + 1) - (width / 2);
-		network->nv_break = (width / 2) * (width + 1);
-		network->dnv = (width / 2 + 1) * (width / 2) + pow(width / 2, 2) - (width / 2);
-		network->num_bounds = 1;
-		network->bound_inds = (uint_t *)malloc((network->num_bounds + 1) *
-																								 sizeof(uint_t));
-		network->bound_inds[0] = 0;
-		network->bound_inds[1] = width / 2;
-		network->bound_verts = (uint_t *)calloc(width / 2, sizeof(uint_t));
-		network->break_dim = network->nv_break + 1;
-	}
-	if (boundary != TORUS_BOUND) {
-		for (uint_t i = 0; i < (width + 1) / 2; i++) {
-			network->bound_verts[i] = i;
-			network->bound_verts[(width + 1) / 2 + i] = network->nv - (width + 1) / 2 + i;
-		} 
-	} else {
-		for (uint_t i = 0; i < width / 2; i++) {
-			network->bound_verts[i] = i;
-		}
-	}
-	network->ev_break =
-			(uint_t *)calloc(2 * network->ne, sizeof(uint_t));
-	network->ev =
-			(uint_t *)calloc(2 * network->ne, sizeof(uint_t));
-	for (uint_t i = 0; i < network->ne; i++) {
-		network->ev_break[2 * i] = edge_to_verts(width, periodic, i, 1);
-		network->ev_break[2 * i + 1] = edge_to_verts(width, periodic, i, 0);
-		network->ev[2 * i] = network->ev_break[2 * i] % network->nv;
-		network->ev[2 * i + 1] = network->ev_break[2 * i + 1] % network->nv;
-	}
-	network->vei =
-			(uint_t *)calloc(network->nv + 1, sizeof(uint_t));
-	network->vei[0] = 0;
-	uint_t pos1 = 0;
-	for (uint_t i = 0; i < network->nv; i++) {
-		bool in_bound = false;
-		for (uint_t j = 0; j < network->num_bounds; j++) {
-			for (uint_t 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;
-				}
+uint_t get_cut_edges(uint_t ne, const uint_t *ev, const double *vx, bool both, uint_t *ce) {
+	uint_t nce = 0;
+
+	for (uint_t i = 0; i < ne; i++) {
+		uint_t v1 = ev[2 * i];
+		uint_t v2 = ev[2 * i + 1];
+
+		double v1y = vx[2 * v1 + 1];
+		double v2y = vx[2 * v2 + 1];
+
+		if (fabs(v1y - v2y) > 0.5) {
+			ce[nce] = i;
+			nce++;
+		} else if (both) {
+			double v1x = vx[2 * v1];
+			double v2x = vx[2 * v2];
+
+			if (fabs(v1x - v2x) > 0.5) {
+				ce[nce] = i;
+				nce++;
 			}
 		}
-		if (in_bound)
-			pos1 += 2;
-		else
-			pos1 += 4;
-
-		network->vei[i + 1] = pos1;
-	}
-	network->ve = (uint_t *)calloc(
-			network->vei[network->nv], sizeof(uint_t));
-	uint_t *vert_counts =
-			(uint_t *)calloc(network->nv, sizeof(uint_t));
-	for (uint_t i = 0; i < network->ne; i++) {
-		uint_t v0 = network->ev[2 * i];
-		uint_t v1 = network->ev[2 * i + 1];
-		network->ve[network->vei[v0] + vert_counts[v0]] =
-				i;
-		network->ve[network->vei[v1] + vert_counts[v1]] =
-				i;
-		vert_counts[v0]++;
-		vert_counts[v1]++;
-	}
-	free(vert_counts);
-
-	network->vx =
-			(double *)malloc(2 * network->nv * sizeof(double));
-	for (uint_t i = 0; i < network->nv; i++) {
-		if (!periodic) {
-		network->vx[2 * i + 1] = ((double)((2 * i + 1) / (width + 1)))/width;
-		network->vx[2 * i] = ((double)((2 * i + 1) % (width + 1)))/width;
-		} else {
-		network->vx[2 * i + 1] = ((double)((2 * i + 1) / (width)))/width;
-		network->vx[2 * i] = ((double)((((2 * i + 1) / (width)+1) % 2) + ((2 * i) % (width))))/width;
-		}
 	}
 
-	network->L = width;
-	network->ex = get_edge_coords(
-			network->ne, network->vx, network->ev);
-
-	network->dev =
-			(uint_t *)malloc(2 * network->ne * sizeof(uint_t));
-	for (uint_t i = 0; i < network->ne; i++) {
-		network->dev[2 * i] =
-				dual_edge_to_verts(width, periodic, i, 0) % network->dnv;
-		network->dev[2 * i + 1] =
-				dual_edge_to_verts(width, periodic, i, 1) % network->dnv;
-	}
-	network->dvx =
-			(double *)malloc(2 * network->dnv * sizeof(double));
-	for (uint_t i = 0; i < network->dnv; i++) {
-		network->dvx[2 * i + 1] =
-				dual_vert_to_coord(width, periodic, i, 0) / width;
-		network->dvx[2 * i] =
-				dual_vert_to_coord(width, periodic, i, 1) / width;
-	}
+	return nce;
+}
 
-	network->voltcurmat = gen_voltcurmat(network->ne,
-																			 network->break_dim,
-																			 network->ev_break, c);
+graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c) {
+	graph_t *g = (graph_t *)calloc(1, sizeof(graph_t));
+	frame_t *f = frame_create(lattice, L, dual);
 
-	network->dvei =
-			get_verts_to_edges_ind(network->dnv, network->ne,
-														 network->dev);
-	network->dve = get_verts_to_edges(
-			network->dnv, network->ne, network->dev,
-			network->dvei);
-	network->spanning_edges = get_spanning_edges(network->ne, network->ev, network->vx, 0.51, &(network->num_spanning_edges));
+	g->L = L;
+	g->boundary = bound;
+	g->ne = f->ne;
+
+	if (bound == TORUS_BOUND) {
+		g->nv = f->nv;
+		g->dnv = f->dnv;
+		g->nb = 1;
+
+		g->ev = f->ev;
+		f->ev = NULL;
+		g->dev = f->dev;
+		f->dev = NULL;
+		g->vx = f->vx;
+		f->vx = NULL;
+		g->dvx = f->dvx;
+		f->dvx = NULL;
+
+		// the boundary for the torus consists of a list of edges required to cut
+		// the torus into a cylinder
+		g->bi = (uint_t *)calloc(2, sizeof(uint_t));
+		g->b = (uint_t *)malloc(g->ne * sizeof(uint_t));
+		g->bi[1] = get_cut_edges(g->ne, g->ev, g->vx, false, g->b);
+		g->bq = (bool *)calloc(g->ne, sizeof(bool));
+		for (uint_t i = 0; i < g->bi[1]; i++) {
+			g->bq[g->b[i]] = true;
+		}
+	} else {
+		uint_t *ce = (uint_t *)malloc(g->ne * sizeof(uint_t));
+		uint_t nce = 0;
 
-	return network;
-}
+		if (bound == CYLINDER_BOUND) {
+			g->nb = 2;
+			nce = get_cut_edges(f->ne, f->ev, f->vx, false, ce);
+		} else {
+			g->nb = 4;
+			nce = get_cut_edges(f->ne, f->ev, f->vx, true, ce);
+		}
 
-uint_t *get_voro_dual_edges(uint_t num_edges,
-																	uint_t num_verts, uint_t *edges,
-																	uint_t *triangles) {
-	uint_t *dual_edges =
-			(uint_t *)malloc(2 * num_edges * sizeof(uint_t));
-	uint_t place = 0;
-	for (uint_t i = 0; i < num_edges; i++) {
-		uint_t v1, v2;
-		v1 = edges[2 * i];
-		v2 = edges[2 * i + 1];
-		if (v1 < num_verts && v2 < num_verts) {
-			bool found_match = false;
-			for (uint_t j = 0; j < 3; j++) {
-				for (uint_t k = 0; k < 3; k++) {
-					uint_t t11, t12, t21, t22;
-					t11 = triangles[3 * v1 + j];
-					t12 = triangles[3 * v1 + ((j + 1) % 3)];
-					t21 = triangles[3 * v2 + k];
-					t22 = triangles[3 * v2 + ((k + 1) % 3)];
-					if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) {
-						dual_edges[2 * place] = t11 < t12 ? t11 : t12;
-						dual_edges[2 * place + 1] = t11 < t12 ? t12 : t11;
-						place++;
-						found_match = true;
-						break;
-					}
-				}
-				if (found_match)
-					break;
+		g->nv = f->nv;
+		g->dnv = f->dnv;
+		g->vx = (double *)malloc(2 * (f->nv + nce) * sizeof(double));
+		g->dvx = (double *)malloc(2 * (f->dnv + nce) * sizeof(double));
+		g->ev = f->ev;
+		g->dev = f->dev;
+		f->ev = NULL;
+		f->dev = NULL;
+		memcpy(g->vx, f->vx, 2 * f->nv * sizeof(double));
+		memcpy(g->dvx, f->dvx, 2 * f->dnv * sizeof(double));
+
+		uint_t nbv = 0;
+		uint_t *bv = (uint_t *)calloc(f->nv, sizeof(uint_t));
+		uint_t *dbv = (uint_t *)calloc(f->dnv, sizeof(uint_t));
+		uint_t nside = 0;
+		bool *side = (bool *)calloc(f->nv, sizeof(bool));
+
+		for (uint_t i = 0; i < nce; i++) {
+			uint_t cv1 = g->ev[2 * ce[i]];
+			uint_t cv2 = g->ev[2 * ce[i] + 1];
+			uint_t dv1 = g->dev[2 * ce[i]];
+			uint_t dv2 = g->dev[2 * ce[i] + 1];
+
+			uint_t cin = 1;
+
+			if (bound == FREE_BOUND && (f->vx[2 * cv2] - f->vx[2 * cv1]) > 0.5) {
+				cin = 0; 
 			}
-		}
-	}
 
-	return dual_edges;
-}
+			uint_t vin = f->vx[2 * cv1 + cin] < f->vx[2 * cv2 + cin] ? 0 : 1;
+			uint_t dvin = f->dvx[2 * dv1 + cin] < f->dvx[2 * dv2 + cin] ? 0 : 1;
 
-graph_t *ini_voro_graph(uint_t L, bound_t boundary, bool use_dual,
-													double *(*genfunc)(uint_t, bound_t, gsl_rng *, uint_t *),
-													cholmod_common *c) {
-	graph_t *g = (graph_t *)calloc(1, sizeof(graph_t));
+			if (bv[g->ev[2 * ce[i] + vin]] == 0) {
+				nbv++;
+				if (cin == 0) {
+					nside++;
+					side[g->ev[2 * ce[i] + vin]] = true;
+				}
 
-	// generate the dual lattice
-	double *lattice;
-	uint_t num;
-	{
-		gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
-		FILE *rf = fopen("/dev/urandom", "r");
-		unsigned long int seed;
-		fread(&seed, sizeof(unsigned long int), 1, rf);
-		fclose(rf);
-		gsl_rng_set(r, seed);
-		lattice = genfunc(L, boundary, r, &num);
-		gsl_rng_free(r);
-	}
+				bv[g->ev[2 * ce[i] + vin]] = g->nv;
 
-	// retrieve a periodic voronoi tesselation of the lattice
-	bool run_periodic;
-	if (boundary == EMBEDDED_BOUND) run_periodic = false;
-	else run_periodic = true;
-	intptr_t *vout = run_voronoi(num, lattice, run_periodic, 0, 1, 0, 1);
-
-	uint_t tmp_num_verts = ((uint_t *)vout[0])[0];
-	uint_t tmp_num_edges = ((uint_t *)vout[0])[1];
-	double *tmp_vert_coords = (double *)vout[2];
-	uint_t *tmp_edges = (uint_t *)vout[3];
-	uint_t *tmp_tris = (uint_t *)vout[5];
-
-	free((void *)vout[0]);
-	free((void *)vout[1]);
-	free((void *)vout[4]);
-	free(vout);
-
-	// get dual edges of the fully periodic graph
-	uint_t *tmp_dual_edges =
-			get_voro_dual_edges(tmp_num_edges, tmp_num_verts, tmp_edges, tmp_tris);
-
-	// when use_dual is specificed, the edge and vertex sets are swapped with the
-	// dual edge and dual vertex sets.  once formally relabelled, everything
-	// works the same way
-	if (use_dual) {
-		uint_t *tmp_tmp_dual_edges = tmp_edges;
-		double *tmp_lattice = tmp_vert_coords;
-		uint_t tmp_num = tmp_num_verts;
-
-		tmp_edges = tmp_dual_edges;
-		tmp_dual_edges = tmp_tmp_dual_edges;
-
-		tmp_vert_coords = lattice;
-		lattice = tmp_lattice;
-
-		tmp_num_verts = num;
-		num = tmp_num;
-	}
+				g->vx[2 * g->nv + cin] = 1 + f->vx[2 * g->ev[2 * ce[i] + vin] + cin];
+				g->vx[2 * g->nv + !cin] = f->vx[2 * g->ev[2 * ce[i] + vin] + !cin];
+				g->ev[2 * ce[i] + vin] = g->nv;
 
-	// prune the edges of the lattice and assign boundary vertices based on the
-	// desired boundary conditions
-	uint_t num_bounds;
-	uint_t num_verts;
-	double *vert_coords;
-	uint_t *bound_inds;
-	uint_t *bound_verts;
-	uint_t num_edges;
-	uint_t *edges;
-	uint_t *dual_edges;
-	switch (boundary) {
-		case FREE_BOUND: {
-			num_bounds = 4;
-			bound_inds = (uint_t *)malloc((1 + num_bounds) * sizeof(uint_t));
-			bound_inds[0] = 0;
-			vert_coords = tmp_vert_coords;
-			num_verts = tmp_num_verts;
-			num_edges = 0;
-			edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
-			dual_edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
-			uint_t num_t, num_b, num_l, num_r;
-			bool *bound_top, *bound_b, *bound_l, *bound_r;
-			num_t = 0; num_b = 0; num_l = 0; num_r = 0;
-			bound_top = (bool *)calloc(num_verts, sizeof(bool));
-			bound_b = (bool *)calloc(num_verts, sizeof(bool));
-			bound_l = (bool *)calloc(num_verts, sizeof(bool));
-			bound_r = (bool *)calloc(num_verts, sizeof(bool));
-			for (uint_t i = 0; i < tmp_num_edges; i++) {
-				uint_t v1, v2;
-				double v1x, v1y, v2x, v2y, dx, dy;
-				v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
-				v1x = vert_coords[2 * v1]; v1y = vert_coords[2 * v1 + 1];
-				v2x = vert_coords[2 * v2]; v2y = vert_coords[2 * v2 + 1];
-				dx = v1x - v2x; dy = v1y - v2y;
-				if (fabs(dy) > 0.5) {
-					if (dy > 0) {
-						if (!bound_top[v1] && !bound_l[v1] && !bound_r[v1]) {
-							bound_top[v1] = true; num_t++;
-						}
-						if (!bound_b[v2] && !bound_l[v2] && !bound_r[v2]) {
-							bound_b[v2] = true; num_b++;
-						}
-					} else {
-						if (!bound_top[v2] && !bound_l[v2] && !bound_r[v2]) {
-							bound_top[v2] = true; num_t++;
-						}
-						if (!bound_b[v1] && !bound_l[v1] && !bound_r[v1]) {
-							bound_b[v1] = true; num_b++;
-						}
-					}
-				} else if (fabs(dx) > 0.5) {
-					if (dx > 0) {
-						if (!bound_r[v1] && !bound_top[v1] && !bound_b[v1]) {
-							bound_r[v1] = true; num_r++;
-						}
-						if (!bound_l[v2] && !bound_top[v2] && !bound_b[v2]) {
-							bound_l[v2] = true; num_l++;
-						}
-					} else {
-						if (!bound_r[v2] && !bound_top[v2] && !bound_b[v2]) {
-							bound_r[v2] = true; num_r++;
-						}
-						if (!bound_l[v1] && !bound_top[v1] && !bound_b[v1]) {
-							bound_l[v1] = true; num_l++;
-						}
-					}
-				} else {
-					edges[2 * num_edges] = v1 < v2 ? v1 : v2;
-					edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
-					uint_t d1 = tmp_dual_edges[2 * i];
-					uint_t 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(uint_t));
-			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;
-			uint_t pos_t, pos_b, pos_l, pos_r;
-			pos_t = 0; pos_b = 0; pos_l = 0; pos_r = 0;
-			for (uint_t i = 0; i < num_verts; i++) {
-				if (bound_top[i]) {
-					bound_verts[pos_t] = i; pos_t++;
-				} else if (bound_b[i]) {
-					bound_verts[num_t + pos_b] = i; pos_b++;
-				} else if (bound_l[i]) {
-					bound_verts[num_t + num_b + pos_l] = i; pos_l++;
-				} else if (bound_r[i]) {
-					bound_verts[num_t + num_b + num_l + pos_r] = i; pos_r++;
-				}
+				g->nv++;
+			} else {
+				g->ev[2 * ce[i] + vin] = bv[g->ev[2 * ce[i] + vin]];
 			}
-			free(bound_l); free(bound_r); free(bound_top); free(bound_b);
-			free(tmp_edges);
-			free(tmp_dual_edges);
-			num_bounds = 2;
-			g->ev_break = edges;
-			g->ev = edges;
-			g->nv_break = num_verts;
-			g->nv = num_verts;
-			g->break_dim = num_verts + num_bounds;
-			break;
-										 }
-		case CYLINDER_BOUND: {
-			num_bounds = 2;
-			bound_inds = (uint_t *)malloc((1 + num_bounds) * sizeof(uint_t));
-			bound_inds[0] = 0;
-			vert_coords = tmp_vert_coords;
-			num_verts = tmp_num_verts;
-			num_edges = 0;
-			edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
-			dual_edges = (uint_t *)malloc(2 * tmp_num_edges * sizeof(uint_t));
-			uint_t num_t, num_b;
-			bool *bound_top, *bound_b;
-			num_t = 0; num_b = 0;
-			bound_top = (bool *)calloc(num_verts, sizeof(bool));
-			bound_b = (bool *)calloc(num_verts, sizeof(bool));
-			for (uint_t i = 0; i < tmp_num_edges; i++) {
-				uint_t v1, v2;
-				double v1y, v2y, dy;
-				v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
-				v1y = vert_coords[2 * v1 + 1]; v2y = vert_coords[2 * v2 + 1];
-				dy = v1y - v2y;
-				if (fabs(dy) > 0.5) {
-					if (dy > 0) {
-						if (!bound_top[v1]) {
-							bound_top[v1] = true; num_t++;
-						}
-						if (!bound_b[v2]) {
-							bound_b[v2] = true; num_b++;
-						}
-					} else {
-						if (!bound_top[v2]) {
-							bound_top[v2] = true; num_t++;
-						}
-						if (!bound_b[v1]) {
-							bound_b[v1] = true; num_b++;
-						}
-					}
+			if (dbv[g->dev[2 * ce[i] + dvin]] == 0) {
+				dbv[g->dev[2 * ce[i] + dvin]] = g->dnv;
+
+				if (f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin] < 0.5) {
+					g->dvx[2 * g->dnv + cin] = 1 + f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin];
 				} else {
-					edges[2 * num_edges] = v1 < v2 ? v1 : v2;
-					edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
-					uint_t d1 = tmp_dual_edges[2 * i];
-					uint_t 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(uint_t));
-			bound_inds[1] = num_t; bound_inds[2] = num_t + num_b;
-			uint_t pos_t, pos_b;
-			pos_t = 0; pos_b = 0;
-			for (uint_t i = 0; i < num_verts; i++) {
-				if (bound_top[i]) {
-					bound_verts[pos_t] = i; pos_t++;
-				} else if (bound_b[i]) {
-					bound_verts[num_t + pos_b] = i; pos_b++;
+					g->dvx[2 * g->dnv + cin] = f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin];
 				}
+				g->dvx[2 * g->dnv + !cin] = f->dvx[2 * g->dev[2 * ce[i] + dvin] + !cin];
+				g->dev[2 * ce[i] + dvin] = g->dnv;
+
+				g->dnv++;
+			} else {
+				g->dev[2 * ce[i] + dvin] = dbv[g->dev[2 * ce[i] + dvin]];
 			}
-			free(bound_top); free(bound_b);
-			free(tmp_edges);
-			free(tmp_dual_edges);
-			g->ev_break = edges;
-			g->ev = edges;
-			g->nv_break = num_verts;
-			g->nv = num_verts;
-			g->break_dim = num_verts + num_bounds;
-			break;
-												 }
-		case TORUS_BOUND: {
-			num_bounds = 1;
-			bound_inds = (uint_t *)malloc((1 + num_bounds) * sizeof(uint_t));
-			bound_inds[0] = 0;
-			num_edges = tmp_num_edges;
-			edges = (uint_t *)malloc(2* num_edges*sizeof(uint_t));
-			for (uint_t i = 0; i < num_edges; i++) {
-				edges[2*i] = tmp_edges[2*i];
-				edges[2*i+1] = tmp_edges[2*i+1];
-			}
-			dual_edges = tmp_dual_edges;
-			bool *bound_top = (bool *)calloc(tmp_num_verts, sizeof(bool));
-			int *edge_change = (int *)calloc(num_edges, sizeof(int));
-			uint_t num_t = 0;
-			for (uint_t i = 0; i < num_edges; i++) {
-				uint_t v1, v2;
-				double v1y, v2y, dy;
-				v1 = edges[2 * i]; v2 = edges[2 * i + 1];
-				v1y = tmp_vert_coords[2 * v1 + 1];
-				v2y = tmp_vert_coords[2 * v2 + 1];
-				dy = v1y - v2y;
-				if (fabs(dy) > 0.5) {
-					if (dy > 0) {
-						edge_change[i] = 1;
-						if (!bound_top[v1]) {
-							bound_top[v1] = true; num_t++;
-						}
-					} else {
-						edge_change[i] = 2;
-						if (!bound_top[v2]) {
-							bound_top[v2] = true; num_t++;
-						}
-					}
-				}
-			}
-			num_verts = tmp_num_verts + num_t;
-			vert_coords = (double *)malloc(2 * num_verts * sizeof(double));
-			bound_verts = malloc(num_t * sizeof(uint_t));
-			bound_inds[1] = num_t;
-			uint_t pos_t = 0;
-			for (uint_t i = 0; i < tmp_num_verts; i++) {
-				vert_coords[2*i] = tmp_vert_coords[2*i];
-				vert_coords[2*i+1] = tmp_vert_coords[2*i+1];
-				if (bound_top[i]) {
-					bound_verts[pos_t] = i;
-					vert_coords[2*(tmp_num_verts + pos_t)] = tmp_vert_coords[2*i];
-					vert_coords[2*(tmp_num_verts + pos_t)+1] = tmp_vert_coords[2*i+1];
-					pos_t++;
-				}
-			}
-			for (uint_t i = 0; i < num_edges; i++) {
-				if (edge_change[i]) {
-					for (uint_t j = 0; j < num_t; j++) {
-						if (edges[2*i+(edge_change[i]-1)] == bound_verts[j]) {
-							edges[2*i+(edge_change[i]-1)] = tmp_num_verts + j;
-							break;
-						}
-					}
-				}
-			}
-			free(tmp_vert_coords);
-			free(bound_top);
-			free(edge_change);
-			g->nv_break = num_verts;
-			g->nv = tmp_num_verts;
-			g->ev_break = edges;
-			g->ev = tmp_edges;
-			g->break_dim = num_verts + 1;
-			break;
-											}
-		case EMBEDDED_BOUND: {
-			num_bounds = 4;
-			bound_inds = (uint_t *)malloc(5 * sizeof(uint_t));
-			bound_verts = (uint_t *)malloc(2 * L * sizeof(uint_t));
-			for (uint_t i = 0; i < 5; i++) bound_inds[i] = i * L / 2;
-			for (uint_t i = 0; i < 2 * L; i++) bound_verts[i] = i;
-			uint_t num_away = 0;
-			for (uint_t i = 0; i < tmp_num_edges; i++) {
-				if (tmp_dual_edges[2*i] > num || tmp_dual_edges[2*i+1] > num) num_away++;
-			}
-			num_edges = (int)tmp_num_edges - (int)num_away;
-			num_verts = tmp_num_verts;
-			edges = tmp_edges;
-			dual_edges = tmp_dual_edges;
-			vert_coords = tmp_vert_coords;
-			g->nv_break = num_verts;
-			g->nv = num_verts;
-			g->ev_break = edges;
-			g->ev = edges;
-			g->break_dim = num_verts + 2;
 		}
-	}
 
-	g->boundary = boundary;
-	g->num_bounds = num_bounds;
-	g->bound_inds = bound_inds;
-	g->bound_verts = bound_verts;
-	g->ne = num_edges;
-	g->dev = dual_edges;
-	g->vx = vert_coords;
+		g->bi = (uint_t *)calloc(1 + g->nb, sizeof(uint_t));
+		g->b = (uint_t *)malloc(2 * nbv * sizeof(uint_t));
 
-	g->vei = get_verts_to_edges_ind(g->nv, g->ne, g->ev);
-	g->ve =
-			get_verts_to_edges(g->nv, g->ne,
-												 g->ev, g->vei);
+		g->bi[1] = ((int_t)nbv - (int_t)nside);
+		g->bi[g->nb] = 2 * nbv;
 
-	g->L = L;
+		if (bound == FREE_BOUND) {
+			g->bi[2] = 2 * ((int_t)nbv - (int_t)nside);
+			g->bi[3] = 2 * (int_t)nbv - (int_t)nside;
+		}
 
-	g->ex = get_edge_coords(
-			g->ne, g->vx, g->ev);
+		uint_t n_ins0 = 0;
+		uint_t n_ins1 = 0;
+
+		g->nbi = (uint_t *)malloc(((int_t)g->nv - (int_t)g->bi[g->nb]) * sizeof(uint_t));
+		g->bni = (uint_t *)malloc((g->nv + 1) * sizeof(uint_t));
+		g->bq = (bool *)calloc(g->nv, sizeof(bool));
+		uint_t n_tmp = 0;
+
+		for (uint_t i = 0; i < f->nv; i++) {
+			g->bni[i] = n_tmp;
+			if (bv[i] != 0) {
+				g->bq[i] = true;
+				g->bq[bv[i]] = true;
+				if (side[i]) {
+					g->b[g->bi[2] + n_ins1] = i;
+					g->b[g->bi[3] + n_ins1] = bv[i];
+					n_ins1++;
+				} else {
+					g->b[g->bi[0] + n_ins0] = i;
+					g->b[g->bi[1] + n_ins0] = bv[i];
+					n_ins0++;
+				}
+			} else {
+				g->nbi[n_tmp] = i;
+				n_tmp++;
+			}
+		}
 
-	free(tmp_tris);
+		for (uint_t i = 0; i < g->nv - f->nv + 1; i++) {
+			g->bni[f->nv + i] = n_tmp;
+		}
 
-	g->dvx = lattice;
-	g->dnv = num;
+		free(bv);
+		free(dbv);
+		free(side);
+	}
 
-	g->voltcurmat = gen_voltcurmat(g->ne,
-																			 g->break_dim,
-																			 g->ev_break, c);
+	g->vei = get_verts_to_edges_ind(g->nv, g->ne, g->ev);
+	g->ve = get_verts_to_edges(g->nv, g->ne, g->ev, g->vei);
+	g->dvei = get_verts_to_edges_ind(g->dnv, g->ne, g->dev);
+	g->dve = get_verts_to_edges(g->dnv, g->ne, g->dev, g->dvei);
 
-	g->dvei =
-			get_verts_to_edges_ind(g->dnv, g->ne,
-														 g->dev);
-	g->dve = get_verts_to_edges(
-			g->dnv, g->ne, g->dev,
-			g->dvei);
+	g->voltcurmat = gen_voltcurmat(g->ne, g->nv, g->ev, c);
+	uint_t num_spanning_edges;
+	g->spanning_edges = get_spanning_edges(g->ne, g->ev, g->vx, 0.5, &num_spanning_edges);
+	g->num_spanning_edges = num_spanning_edges;
 
-	g->spanning_edges = get_spanning_edges(g->ne, g->ev_break, g->vx, 0.5, &(g->num_spanning_edges));
+	frame_free(f);
 
 	return g;
 }
 
-graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c) {
-	bool side_bounds;
-	switch (lattice) {
-		case (VORONOI_LATTICE):
-			return ini_voro_graph(L, bound, dual, genfunc_uniform, c);
-		case (SQUARE_LATTICE):
-			if (bound == EMBEDDED_BOUND) side_bounds = true;
-			else side_bounds = false;
-			return ini_square_network(L, bound, side_bounds, c);
-		case (VORONOI_HYPERUNIFORM_LATTICE):
-			return ini_voro_graph(L, bound, dual, genfunc_hyperuniform, c);
-		default:
-			printf("lattice type unsupported\n");
-			exit(EXIT_FAILURE);
-	}
-}
 
diff --git a/src/graph_free.c b/src/graph_free.c
index e9c55d7..ee30c99 100644
--- a/src/graph_free.c
+++ b/src/graph_free.c
@@ -3,20 +3,17 @@
 
 void graph_free(graph_t *g, cholmod_common *c) {
 	free(g->ev);
-	if (g->ev_break != g->ev) {
-		free(g->ev_break);
-	}
 	free(g->vei);
 	free(g->ve);
-	free(g->bound_inds);
-	free(g->bound_verts);
+	free(g->bi);
+	free(g->b);
 	free(g->vx);
-	free(g->ex);
-	free(g->dev);
 	free(g->dvx);
+	free(g->dev);
 	free(g->dvei);
 	free(g->dve);
-	free(g->spanning_edges);
+	free(g->nbi);
+	free(g->bq);
 	CHOL_F(free_sparse)(&(g->voltcurmat), c);
 	free(g);
 }
diff --git a/src/graph_genfunc.c b/src/graph_genfunc.c
index 8dcd0df..d396206 100644
--- a/src/graph_genfunc.c
+++ b/src/graph_genfunc.c
@@ -1,26 +1,18 @@
 
 #include "fracture.h"
 
-double *genfunc_uniform(unsigned int L, bound_t boundary, gsl_rng *r, unsigned int *num) {
-	*num = 2 * pow(L / 2, 2);
+double *genfunc_uniform(uint_t num, gsl_rng *r) {
+	double *lattice = (double *)malloc(2 * num * sizeof(double));
 
-	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);
+	for (uint_t i = 0; i < 2 * num; i++) {
+		lattice[i] = gsl_rng_uniform(r);
 	}
 
 	return lattice;
 }
 
-double g(double rho, double dist) {
-	return 1 - gsl_sf_exp(-M_PI * rho * dist);
-}
-
-double *genfunc_hyperuniform(unsigned int L, bound_t boundary, gsl_rng *r, unsigned int *num) {
-	*num = 2 * pow(L / 2, 2);
-
-	double *lattice = spheres(*num, 0.8, 0.5, 0.9, 100, 100000);
+double *genfunc_hyperuniform(uint_t num, gsl_rng *r) {
+	double *lattice = spheres(num, 0.8, 0.5, 0.9, 100, 100000);
 
 	return lattice;
 }
diff --git a/src/net.c b/src/net.c
index b61c9ea..a1047e2 100644
--- a/src/net.c
+++ b/src/net.c
@@ -50,35 +50,53 @@ net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, b
 
 	net->graph = g;
 	net->num_broken = 0;
+
 	net->fuses = (bool *)calloc(g->ne, sizeof(bool));
 	assert(net->fuses != NULL);
+
 	net->thres = get_thres(g->ne, beta);
 	net->inf = inf;
 
+	net->dim = g->nv;
+
+	if (g->boundary == TORUS_BOUND) {
+		net->nep = g->ne;
+		net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
+		memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
+	} else {
+		if (vb) {
+			net->dim -= g->bi[g->nb];
+			net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
+			net->nep = 0;
+			for (uint_t i = 0; i < g->ne; i++) {
+				if (!(g->bq[g->ev[2*i]] || g->bq[g->ev[2*i+1]])) {
+					net->evp[2*net->nep] = g->bni[g->ev[2*i]];
+					net->evp[2*net->nep + 1] = g->bni[g->ev[2*i + 1]];
+					net->nep++;
+				}
+			}
+		} else {
+			net->dim += 2;
+			net->evp = (uint_t *)malloc(2 * (g->ne + g->bi[2]) * sizeof(uint_t));
+			memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
+			net->nep = g->ne + g->bi[2];
+
+			for (uint_t i = 0; i < 2; i++) {
+				for (uint_t j = 0; j < g->bi[i+1] - g->bi[i]; j++){ 
+					net->evp[2 * (g->ne + g->bi[i] + j)] = g->b[g->bi[i] + j];
+					net->evp[2 * (g->ne + g->bi[i] + j) + 1] = g->nv + i;
+				}
+			}
+		}
+	}
+
 	net->voltage_bound = vb;
 	net->boundary_cond = bound_set(g, vb, notch_len, c);
 
-	net->adjacency = gen_adjacency(net, false, false, 0, c);
-	net->dual_adjacency = gen_adjacency(net, true, false, 0, c);
-
-	net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t));
-	assert(net->marks != NULL);
-
-	net->dual_marks = (uint_t *)malloc((net->graph->dnv) * sizeof(uint_t));
-	assert(net->dual_marks != NULL);
-
-	for (uint_t i = 0; i < (net->graph->break_dim); i++) {
-		net->marks[i] = 1;
-	}
-	for (uint_t i = 0; i < (net->graph->dnv); i++) {
-		net->dual_marks[i] = i+1;
-	}
-	net->num_components = 1;
-
 	net_notch(net, notch_len, c);
 
 	{
-		cholmod_sparse *laplacian = gen_laplacian(net, c, true);
+		cholmod_sparse *laplacian = gen_laplacian(net, c);
 		net->factor = CHOL_F(analyze)(laplacian, c);
 		CHOL_F(factorize)(laplacian, net->factor, c);
 		CHOL_F(free_sparse)(&laplacian, c);
@@ -102,18 +120,11 @@ net_t *net_copy(const net_t *net, cholmod_common *c) {
 	assert(net_copy->thres != NULL);
 	memcpy(net_copy->thres, net->thres, thres_size);
 
-	size_t marks_size = (net->graph->break_dim) * sizeof(uint_t);
-	net_copy->marks = (uint_t *)malloc(marks_size);
-	assert(net_copy->marks != NULL);
-	memcpy(net_copy->marks, net->marks, marks_size);
-
-	size_t dual_marks_size = (net->graph->dnv) * sizeof(uint_t);
-	net_copy->dual_marks = (uint_t *)malloc(dual_marks_size);
-	assert(net_copy->dual_marks != NULL);
-	memcpy(net_copy->dual_marks, net->dual_marks, dual_marks_size);
+	size_t evp_size = 2 * net->nep * sizeof(uint_t);
+	net_copy->evp = (uint_t *)malloc(thres_size);
+	assert(net_copy->evp != NULL);
+	memcpy(net_copy->evp, net->evp, evp_size);
 
-	net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c);
-	net_copy->dual_adjacency = CHOL_F(copy_sparse)(net->dual_adjacency, c);
 	net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c);
 	net_copy->factor = CHOL_F(copy_factor)(net->factor, c);
 
@@ -124,10 +135,7 @@ void net_free(net_t *net, cholmod_common *c) {
 	free(net->fuses);
 	free(net->thres);
 	CHOL_F(free_dense)(&(net->boundary_cond), c);
-	CHOL_F(free_sparse)(&(net->adjacency), c);
-	CHOL_F(free_sparse)(&(net->dual_adjacency), c);
 	CHOL_F(free_factor)(&(net->factor), c);
-	free(net->marks);
-	free(net->dual_marks);
+	free(net->evp);
 	free(net);
 }
diff --git a/src/net_currents.c b/src/net_currents.c
index 431818f..9bb2874 100644
--- a/src/net_currents.c
+++ b/src/net_currents.c
@@ -3,7 +3,7 @@
 
 double *net_currents(const net_t *net, const double *voltages, cholmod_common *c) {
 	uint_t ne = net->graph->ne;
-	uint_t dim = net->graph->break_dim;
+	uint_t dim = net->graph->nv;
 	cholmod_sparse *voltcurmat = net->graph->voltcurmat;
 
 	cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c);
@@ -27,6 +27,22 @@ double *net_currents(const net_t *net, const double *voltages, cholmod_common *c
 		}
 	}
 
+	if (net->graph->boundary == TORUS_BOUND) {
+		for (uint_t i = 0; i < net->graph->bi[1]; i++) {
+			uint_t e = net->graph->b[i];
+			uint_t v1 = net->graph->ev[2 * e];
+			uint_t v2 = net->graph->ev[2 * e + 1];
+			double v1y = net->graph->vx[2 * v1 + 1];
+			double v2y = net->graph->vx[2 * v2 + 1];
+
+			if (v1y > v2y) {
+				currents[e] += 1;
+			} else {
+				currents[e] -= 1;
+			}
+		}
+	}
+
 	x->x = tmp_x;
 	CHOL_F(free_dense)(&x, c);
 	CHOL_F(free_dense)(&y, c);
diff --git a/src/net_fracture.c b/src/net_fracture.c
index 71e8c4b..e7f18fc 100644
--- a/src/net_fracture.c
+++ b/src/net_fracture.c
@@ -38,7 +38,7 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
 
 		double conductivity = net_conductivity(net, voltages);
 
-		if (conductivity < 1e-12 && net->graph->boundary == TORUS_BOUND) {
+		if (conductivity < cutoff) {
 			free(voltages);
 			free(currents);
 			break;
@@ -60,14 +60,6 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
 		free(currents);
 
 		break_edge(net, last_broke, c);
-
-		if (net->num_components > 1 && net->graph->boundary == TORUS_BOUND) {
-			break;
-		}
-
-		if (net->marks[net->graph->nv] != net->marks[net->graph->nv + 1] && net->graph->boundary != TORUS_BOUND) {
-			break;
-		}
 	}
 
 	return data;
diff --git a/src/net_voltages.c b/src/net_voltages.c
index dedf5b2..c3537a5 100644
--- a/src/net_voltages.c
+++ b/src/net_voltages.c
@@ -12,11 +12,29 @@ double *net_voltages(const net_t *net, cholmod_common *c) {
 		exit(EXIT_FAILURE);
 	}
 
-	double *voltages = (double *)x->x;
+	double *t_voltages = (double *)x->x;
 	x->x = NULL;
-
 	CHOL_F(free_dense)(&x, c);
 
-	return voltages;
+	graph_t *g = net->graph;
+
+	if (g->boundary == TORUS_BOUND) {
+		return t_voltages;
+	} else if (net->voltage_bound) {
+		double *voltages = (double *)malloc(g->nv * sizeof(double));
+		for (uint_t i = 0; i < g->nv - g->bi[g->nb]; i++) {
+			voltages[g->nbi[i]] = t_voltages[i];
+		}
+		for (uint_t i = 0; i < 2; i++) {
+			for (uint_t j = 0; j < g->bi[i + 1] - g->bi[i]; j++) {
+				voltages[g->b[g->bi[i] + j]] = 1 - i;
+			}
+		}
+
+		free(t_voltages);
+		return voltages;
+	} else {
+		return t_voltages;
+	}
 }