more major refactoring

25 files changed, 465 insertions, 519 deletions

diff --git a/src/bin_values.c b/src/bin_values.c
index 8cde548..70009c1 100644
--- a/src/bin_values.c
+++ b/src/bin_values.c
@@ -6,8 +6,8 @@ double *bin_values(graph_t *network, unsigned int width, double *values) {
 	unsigned int *num_binned = calloc(pow(width, 2), sizeof(unsigned int));
 	for (unsigned int i = 0; i < network->ne; 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;
+			unsigned int x = ((int)(network->ex[2 * i] * width)) % width;
+			unsigned int y = ((int)(network->ex[2 * i + 1] * width)) % width;
 			binned[width * x + y] += fabs(values[i]);
 			num_binned[width * x + y]++;
 		}
diff --git a/src/bound_set.c b/src/bound_set.c
new file mode 100644
index 0000000..9c67f82
--- /dev/null
+++ b/src/bound_set.c
@@ -0,0 +1,62 @@
+
+#include "fracture.h"
+
+double th_p(double x, double y, double th) {
+	if (x >= 0 && y >= 0) return th;
+	else if (x <  0 && y >= 0) return M_PI - th;
+	else if (x <  0 && y <  0) return th - M_PI;
+	else return -th;
+
+}
+
+double u_y(double x, double y) {
+	double r = sqrt(pow(x, 2) + pow(y, 2));
+	double th = th_p(x, y, atan(fabs(y / x)));
+
+	return sqrt(r) * sin(th / 2);
+}
+
+void bound_set_embedded(double *bound, uint_t L, double notch_len) {
+	for (uint_t i = 0; i < L / 2; i++) {
+		double x1, y1, x2, y2, x3, y3, x4, y4;
+		x1 = (2. * i + 1.) / L - notch_len; y1 = 0.5 - 1.;
+		x2 = (2. * i + 1.) / L - notch_len; y2 = 0.5 - 0.;
+		y3 = (2. * i + 1.) / L - 0.5; x3 = 0.5 - 1.;
+		y4 = (2. * i + 1.) / L - 0.5; x4 = 0.5 - 0.;
+
+		bound[i] = u_y(x1, y1);
+		bound[L / 2 + i] = u_y(x2, y2);
+		bound[L + i] = u_y(x3, y3);
+		bound[3 * L / 2 + i] = u_y(x4, y4);
+	}
+}
+
+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);
+	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;
+			}
+			bound[g->bound_verts[0]] = 1;
+			break;
+		case EMBEDDED_BOUND:
+			bound_set_embedded(bound, g->L, notch_len);
+			break;
+		default:
+			if (vb) {
+				for (uint_t i = 0; i < g->bound_inds[1]; i++) {
+					bound[g->bound_verts[i]] = 1;
+				}
+			} else {
+				bound[0] = 1;
+				bound[g->nv] = 1;
+				bound[g->nv + 1] = -1;
+			}
+	}
+
+	return boundary;
+}
diff --git a/src/break_edge.c b/src/break_edge.c
index 7b94160..4e1559b 100644
--- a/src/break_edge.c
+++ b/src/break_edge.c
@@ -3,7 +3,6 @@
 
 bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) {
 	instance->fuses[edge] = true;
-	instance->num_remaining_edges--;
 
 	unsigned int v1 = instance->graph->ev_break[2 * edge];
 	unsigned int v2 = instance->graph->ev_break[2 * edge + 1];
@@ -56,10 +55,10 @@ bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) {
 					double v1x, v1y, v2x, v2y;
 					v1 = instance->graph->dev[2 * ee + !side];
 					v2 = instance->graph->dev[2 * ee + side];
-					v1x = instance->graph->dual_vert_coords[2 * v1];
-					v1y = instance->graph->dual_vert_coords[2 * v1 + 1];
-					v2x = instance->graph->dual_vert_coords[2 * v2];
-					v2y = instance->graph->dual_vert_coords[2 * v2 + 1];
+					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) {
diff --git a/src/compare_voronoi_fracture.c b/src/compare_voronoi_fracture.c
index 2edcae2..91fdcea 100644
--- a/src/compare_voronoi_fracture.c
+++ b/src/compare_voronoi_fracture.c
@@ -66,10 +66,10 @@ int main(int argc, char *argv[]) {
 	(&c)->supernodal = CHOLMOD_SIMPLICIAL;
 
 
-	graph_t *network = ini_voronoi_network(L, false, boundary, genfunc_hyperuniform, &c);
+	graph_t *network = ini_voro_graph(L, false, boundary, genfunc_hyperuniform, &c);
 	net_t *perm_voltage_instance = create_instance(network, inf, true, true, &c);
 	net_t *perm_current_instance = create_instance(network, inf, false, true, &c);
-	double *fuse_thres = gen_fuse_thres(network->ne, network->edge_coords, beta, beta_scaling_flat);
+	double *fuse_thres = gen_fuse_thres(network->ne, network->ex, beta, beta_scaling_flat);
 	net_t *voltage_instance = copy_instance(perm_voltage_instance, &c);
 	net_t *current_instance = copy_instance(perm_current_instance, &c);
 	data_t *breaking_data_voltage = fracture_network(voltage_instance, fuse_thres, &c, cutoff);
@@ -82,8 +82,8 @@ int main(int argc, char *argv[]) {
 
 	FILE *net_out = fopen("network.txt", "w");
 	for (unsigned int j = 0; j < network->nv; j++) {
-		fprintf(net_out, "%f %f ", network->vert_coords[2 * j],
-						network->vert_coords[2 * j + 1]);
+		fprintf(net_out, "%f %f ", network->vx[2 * j],
+						network->vx[2 * j + 1]);
 	}
 	fprintf(net_out, "\n");
 	for (unsigned int j = 0; j < network->ne; j++) {
@@ -92,8 +92,8 @@ int main(int argc, char *argv[]) {
 	}
 	fprintf(net_out, "\n");
 	for (unsigned int j = 0; j < network->dnv; j++) {
-		fprintf(net_out, "%f %f ", network->dual_vert_coords[2 * j],
-						network->dual_vert_coords[2 * j + 1]);
+		fprintf(net_out, "%f %f ", network->dvx[2 * j],
+						network->dvx[2 * j + 1]);
 	}
 	fprintf(net_out, "\n");
 	for (unsigned int j = 0; j < network->ne; j++) {
diff --git a/src/corr_test.c b/src/corr_test.c
index e0de2d8..a8eaff5 100644
--- a/src/corr_test.c
+++ b/src/corr_test.c
@@ -9,10 +9,8 @@ int main() {
 	unsigned int n = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2);
 
 	graph_t *network = ini_square_network(width, true, false, &c);
-	net_t *instance = create_instance(network, 1e-14, true, true, &c);
-	double *fuse_thres = gen_fuse_thres(network->ne, network->edge_coords,
-																			0.001, beta_scaling_flat);
-	fracture_network(instance, fuse_thres, &c, 1e-10);
+	net_t *instance = net_create(network, 1e-14, 0.001, 0, true, &c);
+	net_fracture(instance, &c, 1e-10);
 
 	double *corr = get_corr(instance, NULL, &c);
 
@@ -21,8 +19,8 @@ int main() {
 	}
 	printf("\n");
 
-	free_instance(instance, &c);
-	free_net(network, &c);
+	net_free(instance, &c);
+	graph_free(network, &c);
 
 	CHOL_F(finish)(&c);
 
diff --git a/src/correlations.c b/src/correlations.c
index 4fc44ac..7fedfbf 100644
--- a/src/correlations.c
+++ b/src/correlations.c
@@ -298,7 +298,7 @@ double *get_corr(net_t *instance, unsigned int **dists, cholmod_common *c) {
 /*
 double *get_space_corr(net_t *instance, cholmod_common *c) {
 	unsigned int nv = instance->graph->dnv;
-	double *vc = instance->graph->dual_vert_coords;
+	double *vc = instance->graph->dvx;
 	unsigned int ne = instance->graph->ne;
 	unsigned int *ev = instance->graph->dev;
 	double *corr = calloc(nv, sizeof(unsigned int));
diff --git a/src/course_grain_square.c b/src/course_grain_square.c
index 9c9ca6e..6587e83 100644
--- a/src/course_grain_square.c
+++ b/src/course_grain_square.c
@@ -84,7 +84,7 @@ int main(int argc, char *argv[]) {
 		net_t *instance = NULL;
 		while (breaking_data == NULL) {
 			double *fuse_thres = gen_fuse_thres(
-					network->ne, network->edge_coords, beta, beta_scaling_flat);
+					network->ne, network->ex, beta, beta_scaling_flat);
 			instance = copy_instance(perm_instance, &c);
 			breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
 			free_instance(instance, &c);
diff --git a/src/coursegrain.c b/src/coursegrain.c
index 1e0d5a7..3e9db1a 100644
--- a/src/coursegrain.c
+++ b/src/coursegrain.c
@@ -5,7 +5,7 @@ net_t *coursegrain_square(net_t *instance, graph_t *network_p, cholmod_common *c
 	unsigned int width = sqrt(instance->graph->ne);
 	assert(width % 4 == 0);
 
-	net_t *instance_p = create_instance(network_p, instance->inf,
+	net_t *instance_p = create_instance(network_p, instance->inf, 1,
 																			instance->voltage_bound, true, c);
 
 	unsigned int width_p = width / 2;
diff --git a/src/cracking_ini.c b/src/cracking_ini.c
deleted file mode 100644
index 93e5765..0000000
--- a/src/cracking_ini.c
+++ /dev/null
@@ -1,77 +0,0 @@
-
-#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;
-}
-
-void gen_voro_crack(net_t *instance, double crack_len, cholmod_common *c) {
-	for (uint_t i = 0; i < instance->graph->ne; i++) {
-		uint_t v1, v2;
-		double v1x, v1y, v2x, v2y, dx, dy;
-
-		v1 = instance->graph->ev[2 * i];
-		v2 = instance->graph->ev[2 * i + 1];
-
-		v1x = instance->graph->vert_coords[2 * v1];
-		v1y = instance->graph->vert_coords[2 * v1 + 1];
-		v2x = instance->graph->vert_coords[2 * v2];
-		v2y = instance->graph->vert_coords[2 * v2 + 1];
-
-		dx = v1x - v2x;
-		dy = v1y - v2y;
-
-		if (((v1y > 0.5 && v2y < 0.5) || (v1y < 0.5 && v2y > 0.5)) && fabs(dy) < 0.5) {
-			if (v1x + dx / dy * (v1y - 0.5) <= crack_len) {
-				break_edge(instance, i, c);
-			}
-		}
-	}
-}
-
-bool gen_crack(net_t *instance, double crack_len, double crack_width,
-							 cholmod_common *c) {
-	assert(instance != NULL);
-	bool *fuses = instance->fuses;
-	assert(fuses != NULL);
-	const graph_t *network = instance->graph;
-	unsigned int num_edges = network->ne;
-	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
index f750cb1..ef36b0e 100644
--- a/src/current_scaling.c
+++ b/src/current_scaling.c
@@ -90,13 +90,13 @@ int main(int argc, char *argv[]) {
 	graph_t *network = ini_square_network(width, false, true, &c);
 	net_t *perm_instance =
 			create_instance(network, inf, voltage_bound, false, &c);
-	gen_crack(perm_instance, crack_len, 1, &c);
+	gen_crack(perm_instance, crack_len, &c);
 	finish_instance(perm_instance, &c);
 
 	if (voltage_bound) {
 		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
 		net_t *tmp_instance = create_instance(network, inf, false, false, &c);
-		gen_crack(tmp_instance, crack_len, 1, &c);
+		gen_crack(tmp_instance, crack_len, &c);
 		finish_instance(tmp_instance, &c);
 		double *voltage = get_voltage(tmp_instance, &c);
 
@@ -151,7 +151,7 @@ int main(int argc, char *argv[]) {
 			data_t *breaking_data = NULL;
 			while (breaking_data == NULL) {
 				double *fuse_thres = gen_fuse_thres(
-						network->ne, network->edge_coords, beta, beta_scaling_flat);
+						network->ne, network->ex, beta, beta_scaling_flat);
 				net_t *instance = copy_instance(perm_instance, &c);
 				breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
 				free_instance(instance, &c);
diff --git a/src/fracture.c b/src/fracture.c
index 23bce03..6053146 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -166,14 +166,14 @@ int main(int argc, char *argv[]) {
 																						genfunc_uniform, &c)) == NULL)
 				;
 			perm_instance = create_instance(network, inf, voltage_bound, false, &c);
-			gen_crack(perm_instance, crack_len, crack_width, &c);
+			gen_crack(perm_instance, crack_len, &c);
 			finish_instance(perm_instance, &c);
 		}
 	} else {
 		network = ini_square_network(width, periodic, side_bounds, &c);
 		crack_width = 1;
 		perm_instance = create_instance(network, inf, voltage_bound, false, &c);
-		gen_crack(perm_instance, crack_len, crack_width, &c);
+		gen_crack(perm_instance, crack_len, &c);
 		finish_instance(perm_instance, &c);
 		c_dist_size = network->dnv;
 		a_dist_size = network->nv;
@@ -286,14 +286,14 @@ int main(int argc, char *argv[]) {
 																							genfunc_uniform, &c)) == NULL)
 					;
 				perm_instance = create_instance(network, inf, voltage_bound, false, &c);
-				gen_crack(perm_instance, crack_len, crack_width, &c);
+				gen_crack(perm_instance, crack_len, &c);
 				finish_instance(perm_instance, &c);
 				if (supplied_bound) {
 					voronoi_bound_ini(perm_instance, square_bound, width);
 				}
 			}
 			double *fuse_thres = gen_fuse_thres(
-					network->ne, network->edge_coords, beta, beta_scaling_flat);
+					network->ne, network->ex, beta, beta_scaling_flat);
 			finst *instance = copy_instance(perm_instance, &c);
 			breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
 			free_instance(instance, &c);
@@ -414,8 +414,8 @@ int main(int argc, char *argv[]) {
 		if (network_out) {
 			FILE *net_out = fopen("network.txt", "w");
 			for (unsigned int i = 0; i < network->nv; i++) {
-				fprintf(net_out, "%f %f ", network->vert_coords[2 * i],
-								network->vert_coords[2 * i + 1]);
+				fprintf(net_out, "%f %f ", network->vx[2 * i],
+								network->vx[2 * i + 1]);
 			}
 			fprintf(net_out, "\n");
 			for (unsigned int i = 0; i < network->ne; i++) {
@@ -424,8 +424,8 @@ int main(int argc, char *argv[]) {
 			}
 			fprintf(net_out, "\n");
 			for (unsigned int i = 0; i < network->dnv; i++) {
-				fprintf(net_out, "%f %f ", network->dual_vert_coords[2 * i],
-								network->dual_vert_coords[2 * i + 1]);
+				fprintf(net_out, "%f %f ", network->dvx[2 * i],
+								network->dvx[2 * i + 1]);
 			}
 			fprintf(net_out, "\n");
 			for (unsigned int i = 0; i < network->ne; i++) {
diff --git a/src/fracture.h b/src/fracture.h
index a9973ea..17ab8fd 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -37,38 +37,35 @@ typedef enum bound_t {
 } bound_t;
 
 typedef struct {
-	uint_t ne;          // number of edges
-	uint_t nv;          // number of vertices
-	uint_t dnv;         // number of dual vertices
+	uint_t ne;
+	uint_t nv;
+	uint_t dnv;
 	uint_t nv_break;
 	uint_t num_bounds;
-	// 2*ne length array.  edge i connects vertices ev[2*i], ev[2*i+1] < nv
 	uint_t *ev; 
 	uint_t *ev_break;
-	// nv+1 length array.  vertex i has degree vei[i+1]-vei[i], vei[0] = 0
 	uint_t *vei;
-	// vei[nv] length array.  vertex i is connected to edges ve[vei[i]+j], 0 <= j < vei[i+1]-vei[i]
 	uint_t *ve;
 	uint_t *bound_inds;
 	uint_t *bound_verts;
-	double *vert_coords;
-	double *edge_coords;
+	double *vx;
+	double *ex;
 	uint_t *dev;
 	uint_t *dvei;
 	uint_t *dve;
-	double *dual_vert_coords;
+	double *dvx;
 	uint_t num_spanning_edges;
 	uint_t *spanning_edges;
-	double L;
+	uint_t L;
 	uint_t break_dim;
 	cholmod_sparse *voltcurmat;
 	bound_t boundary;
 } graph_t;
 
 typedef struct {
-	graph_t *graph;
-	unsigned int num_remaining_edges;
+	const graph_t *graph;
 	bool *fuses;
+	double *thres;
 	double inf;
 	cholmod_dense *boundary_cond;
 	cholmod_factor *factor;
@@ -112,18 +109,13 @@ double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert,
 bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2,
 									 cholmod_common *c);
 
-data_t *fracture_network(net_t *instance, double *fuse_thres,
-														 cholmod_common *c, double cutoff);
+data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff);
 
 double *get_current(const net_t *instance, cholmod_common *c);
 double *get_current_v(const net_t *instance, double *voltages, cholmod_common *c);
 double *get_voltage(const net_t *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(net_t *instance, double crack_len, double crack_width,
-							 cholmod_common *c);
+void net_notch(net_t *net, double notch_len, cholmod_common *c);
 
 void update_boundary(net_t *instance, const double *avg_field);
 
@@ -137,18 +129,17 @@ double update_beta(double beta, unsigned int width, const double *stress,
 cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts,
 															 unsigned int *edges_to_verts, cholmod_common *c);
 
-net_t *copy_instance(const net_t *instance, cholmod_common *c);
+net_t *net_copy(const net_t *net, cholmod_common *c);
 
 graph_t *ini_square_network(unsigned int width, bound_t boundary, bool side_bounds,
 												 cholmod_common *c);
 
-void free_net(graph_t *network, cholmod_common *c);
-void free_instance(net_t *instance, cholmod_common *c);
+void graph_free(graph_t *network, cholmod_common *c);
+void net_free(net_t *instance, cholmod_common *c);
 
-net_t *create_instance(graph_t *network, double inf, bool voltage_bound,
-											 bool startnow, cholmod_common *c);
+net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c);
 
-graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
+graph_t *ini_voro_graph(unsigned int L, bound_t boundary, bool use_dual,
 													double *(*genfunc)(unsigned int, bound_t, gsl_rng *, unsigned int *),
 													cholmod_common *c);
 
@@ -178,7 +169,7 @@ double *get_corr(net_t *instance, unsigned int **dists, cholmod_common *c);
 
 double *bin_values(graph_t *network, unsigned int width, double *values);
 
-void voronoi_bound_ini(net_t *instance, uint_t L, double crack_len);
+cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_common *c);
 
 data_t *alloc_break_data(unsigned int num_edges);
 void free_break_data(data_t *data);
@@ -186,4 +177,3 @@ void update_break_data(data_t *data, unsigned int last_broke, double strength, d
 
 double get_conductivity(net_t *inst, double *current, cholmod_common *c);
 
-void gen_voro_crack(net_t *instance, double crack_len, cholmod_common *c);
diff --git a/src/fracture_network.c b/src/fracture_network.c
deleted file mode 100644
index 428b092..0000000
--- a/src/fracture_network.c
+++ /dev/null
@@ -1,67 +0,0 @@
-
-#include "fracture.h"
-
-int inc_break_fuses(net_t *instance, double *thres, double *field,
-														 double cutoff) {
-	unsigned int size = (instance->graph)->ne;
-
-	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;
-}
-
-data_t *fracture_network(net_t *instance, double *fuse_thres,
-														 cholmod_common *c, double cutoff) {
-	unsigned int num_edges = instance->graph->ne;
-	unsigned int num_verts = instance->graph->nv;
-
-	data_t *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 || conductivity < 1e-8) {
-			printf("whoops %u\n\n", breaking_data->num_broken);
-			free(voltages);
-			free(field);
-			break;
-		}
-
-		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 && instance->graph->boundary == TORUS_BOUND) {
-			break;
-		}
-
-		if (instance->marks[num_verts] != instance->marks[num_verts + 1] && instance->graph->boundary != TORUS_BOUND) {
-			break;
-		}
-	}
-
-	return breaking_data;
-}
-
diff --git a/src/free_network.c b/src/free_network.c
index 5fd3290..2001479 100644
--- a/src/free_network.c
+++ b/src/free_network.c
@@ -1,7 +1,7 @@
 
 #include "fracture.h"
 
-void free_net(graph_t *network, cholmod_common *c) {
+void graph_free(graph_t *network, cholmod_common *c) {
 	free(network->ev);
 	if (network->ev_break != network->ev) {
 		free(network->ev_break);
@@ -10,10 +10,10 @@ void free_net(graph_t *network, cholmod_common *c) {
 	free(network->ve);
 	free(network->bound_inds);
 	free(network->bound_verts);
-	free(network->vert_coords);
-	free(network->edge_coords);
+	free(network->vx);
+	free(network->ex);
 	free(network->dev);
-	free(network->dual_vert_coords);
+	free(network->dvx);
 	free(network->dvei);
 	free(network->dve);
 	free(network->spanning_edges);
diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c
index 6034a82..4196937 100644
--- a/src/gen_laplacian.c
+++ b/src/gen_laplacian.c
@@ -3,14 +3,9 @@
 
 cholmod_sparse *gen_adjacency(const net_t *instance, bool dual, bool breakv,
 															unsigned int pad, cholmod_common *c) {
-	unsigned int ne;
-	if (dual)
-		ne = ((int)instance->graph->ne);
-	else {
-		ne = instance->num_remaining_edges;
-		if (!breakv && instance->graph->boundary != TORUS_BOUND) {
-			ne += instance->graph->bound_inds[2];
-		}
+	unsigned int ne = instance->graph->ne;
+	if (!breakv && instance->graph->boundary != TORUS_BOUND && !dual) {
+		ne += instance->graph->bound_inds[2];
 	}
 
 	unsigned int nnz = 2 * ne;
@@ -58,7 +53,7 @@ cholmod_sparse *gen_adjacency(const net_t *instance, bool dual, bool breakv,
 			ai[2 * count + 1] = 1;
 
 			count++;
-		} else if (dual) {
+		} else {
 			unsigned int v1 = etv[2 * i];
 			unsigned int v2 = etv[2 * i + 1];
 
@@ -99,7 +94,7 @@ cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
 															bool symmetric) {
 	const graph_t *network = instance->graph;
 	unsigned int num_verts = network->nv_break;
-	double *vert_coords = network->vert_coords;
+	double *vert_coords = network->vx;
 	unsigned int num_bounds = network->num_bounds;
 	double inf = instance->inf;
 	bool voltage_bound = instance->voltage_bound;
diff --git a/src/get_conductivity.c b/src/get_conductivity.c
index 793987e..8c4d228 100644
--- a/src/get_conductivity.c
+++ b/src/get_conductivity.c
@@ -9,8 +9,8 @@ double get_conductivity(net_t *inst, double *voltage, cholmod_common *c) {
 			if (!inst->fuses[e]) {
 				unsigned int v1 = inst->graph->ev[2*e];
 				unsigned int v2 = inst->graph->ev[2*e+1];
-				double v1y = inst->graph->vert_coords[2 * v1 + 1];
-				double v2y = inst->graph->vert_coords[2 * v2 + 1];
+				double v1y = inst->graph->vx[2 * v1 + 1];
+				double v2y = inst->graph->vx[2 * v2 + 1];
 				unsigned int s1 = v1y < v2y ? v1 : v2;
 				unsigned int s2 = v1y < v2y ? v2 : v1;
 				tot_cur += voltage[s1] - voltage[s2];
diff --git a/src/homo_square_fracture.c b/src/homo_square_fracture.c
index e3e8ad3..b301136 100644
--- a/src/homo_square_fracture.c
+++ b/src/homo_square_fracture.c
@@ -193,7 +193,7 @@ int main(int argc, char *argv[]) {
 	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->ne, network->edge_coords, beta, beta_scaling_flat);
+		double *fuse_thres = gen_fuse_thres(network->ne, network->ex, beta, beta_scaling_flat);
 		net_t *instance = copy_instance(perm_instance, &c);
 		data_t *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
 		free_instance(instance, &c);
@@ -287,8 +287,8 @@ int main(int argc, char *argv[]) {
 		if (save_network) {
 			FILE *net_out = fopen("network.txt", "w");
 			for (unsigned int j = 0; j < network->nv; j++) {
-				fprintf(net_out, "%f %f ", network->vert_coords[2 * j],
-								tmp_instance->graph->vert_coords[2 * j + 1]);
+				fprintf(net_out, "%f %f ", network->vx[2 * j],
+								tmp_instance->graph->vx[2 * j + 1]);
 			}
 			fprintf(net_out, "\n");
 			for (unsigned int j = 0; j < tmp_instance->graph->ne; j++) {
@@ -297,8 +297,8 @@ int main(int argc, char *argv[]) {
 			}
 			fprintf(net_out, "\n");
 			for (unsigned int j = 0; j < tmp_instance->graph->dnv; j++) {
-				fprintf(net_out, "%f %f ", tmp_instance->graph->dual_vert_coords[2 * j],
-								tmp_instance->graph->dual_vert_coords[2 * j + 1]);
+				fprintf(net_out, "%f %f ", tmp_instance->graph->dvx[2 * j],
+								tmp_instance->graph->dvx[2 * j + 1]);
 			}
 			fprintf(net_out, "\n");
 			for (unsigned int j = 0; j < tmp_instance->graph->ne; j++) {
diff --git a/src/homo_voronoi_fracture.c b/src/homo_voronoi_fracture.c
index 6171480..85a9c2b 100644
--- a/src/homo_voronoi_fracture.c
+++ b/src/homo_voronoi_fracture.c
@@ -211,9 +211,9 @@ int main(int argc, char *argv[]) {
 	for (uint32_t i = 0; i < N; i++) {
 		printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1, N);
 
-		graph_t *network = ini_voronoi_network(L, boundary, use_dual, genfunc_hyperuniform, &c);
+		graph_t *network = ini_voro_graph(L, boundary, use_dual, genfunc_hyperuniform, &c);
 		net_t *perm_instance = create_instance(network, inf, use_voltage_boundaries, true, &c);
-		double *fuse_thres = gen_fuse_thres(network->ne, network->edge_coords, beta, beta_scaling_flat);
+		double *fuse_thres = gen_fuse_thres(network->ne, network->ex, beta, beta_scaling_flat);
 		net_t *instance = copy_instance(perm_instance, &c);
 		data_t *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
 		free_instance(instance, &c);
@@ -299,8 +299,8 @@ int main(int argc, char *argv[]) {
 		if (save_network) {
 			FILE *net_out = fopen("network.txt", "w");
 			for (uint_t j = 0; j < network->nv; j++) {
-				fprintf(net_out, "%f %f ", network->vert_coords[2 * j],
-								tmp_instance->graph->vert_coords[2 * j + 1]);
+				fprintf(net_out, "%f %f ", network->vx[2 * j],
+								tmp_instance->graph->vx[2 * j + 1]);
 			}
 			fprintf(net_out, "\n");
 			for (uint_t j = 0; j < tmp_instance->graph->ne; j++) {
@@ -309,8 +309,8 @@ int main(int argc, char *argv[]) {
 			}
 			fprintf(net_out, "\n");
 			for (uint_t j = 0; j < tmp_instance->graph->dnv; j++) {
-				fprintf(net_out, "%f %f ", tmp_instance->graph->dual_vert_coords[2 * j],
-								tmp_instance->graph->dual_vert_coords[2 * j + 1]);
+				fprintf(net_out, "%f %f ", tmp_instance->graph->dvx[2 * j],
+								tmp_instance->graph->dvx[2 * j + 1]);
 			}
 			fprintf(net_out, "\n");
 			for (uint_t j = 0; j < tmp_instance->graph->ne; j++) {
diff --git a/src/ini_network.c b/src/ini_network.c
index ebfea76..cfb8d53 100644
--- a/src/ini_network.c
+++ b/src/ini_network.c
@@ -219,21 +219,21 @@ graph_t *ini_square_network(unsigned int width, bound_t boundary, bool side_boun
 	}
 	free(vert_counts);
 
-	network->vert_coords =
+	network->vx =
 			(double *)malloc(2 * network->nv * sizeof(double));
 	for (unsigned int i = 0; i < network->nv; 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;
+		network->vx[2 * i] = ((double)((2 * i + 1) / (width + 1)))/width;
+		network->vx[2 * i + 1] = ((double)((2 * i + 1) % (width + 1)))/width;
 		} else {
-		network->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->vx[2 * i] = ((double)((2 * i + 1) / (width)))/width;
+		network->vx[2 * i + 1] = ((double)((((2 * i + 1) / (width)+1) % 2) + ((2 * i) % (width))))/width;
 		}
 	}
 
 	network->L = width;
-	network->edge_coords = get_edge_coords(
-			network->ne, network->vert_coords, network->ev);
+	network->ex = get_edge_coords(
+			network->ne, network->vx, network->ev);
 
 	network->dev =
 			(unsigned int *)malloc(2 * network->ne * sizeof(unsigned int));
@@ -243,12 +243,12 @@ graph_t *ini_square_network(unsigned int width, bound_t boundary, bool side_boun
 		network->dev[2 * i + 1] =
 				dual_edge_to_verts(width, periodic, i, 1) % network->dnv;
 	}
-	network->dual_vert_coords =
+	network->dvx =
 			(double *)malloc(2 * network->dnv * sizeof(double));
 	for (unsigned int i = 0; i < network->dnv; i++) {
-		network->dual_vert_coords[2 * i] =
+		network->dvx[2 * i] =
 				2*dual_vert_to_coord(width, periodic, i, 0);
-		network->dual_vert_coords[2 * i + 1] =
+		network->dvx[2 * i + 1] =
 				2*dual_vert_to_coord(width, periodic, i, 1);
 	}
 
@@ -262,7 +262,7 @@ graph_t *ini_square_network(unsigned int width, bound_t boundary, bool side_boun
 	network->dve = get_verts_to_edges(
 			network->dnv, network->ne, network->dev,
 			network->dvei);
-	network->spanning_edges = get_spanning_edges(network->ne, network->ev_break, network->vert_coords, 0.51, &(network->num_spanning_edges));
+	network->spanning_edges = get_spanning_edges(network->ne, network->ev_break, network->vx, 0.51, &(network->num_spanning_edges));
 
 	return network;
 }
@@ -304,10 +304,10 @@ unsigned int *get_voro_dual_edges(unsigned int num_edges,
 	return dual_edges;
 }
 
-graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
+graph_t *ini_voro_graph(unsigned int L, bound_t boundary, bool use_dual,
 													double *(*genfunc)(unsigned int, bound_t, gsl_rng *, unsigned int *),
 													cholmod_common *c) {
-	graph_t *network = (graph_t *)calloc(1, sizeof(graph_t));
+	graph_t *g = (graph_t *)calloc(1, sizeof(graph_t));
 
 	// generate the dual lattice
 	double *lattice;
@@ -459,11 +459,11 @@ graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			free(tmp_edges);
 			free(tmp_dual_edges);
 			num_bounds = 2;
-			network->ev_break = edges;
-			network->ev = edges;
-			network->nv_break = num_verts;
-			network->nv = num_verts;
-			network->break_dim = num_verts + num_bounds;
+			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: {
@@ -526,11 +526,11 @@ graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			free(bound_top); free(bound_b);
 			free(tmp_edges);
 			free(tmp_dual_edges);
-			network->ev_break = edges;
-			network->ev = edges;
-			network->nv_break = num_verts;
-			network->nv = num_verts;
-			network->break_dim = num_verts + num_bounds;
+			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: {
@@ -596,11 +596,11 @@ graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			free(tmp_vert_coords);
 			free(bound_top);
 			free(edge_change);
-			network->nv_break = num_verts;
-			network->nv = tmp_num_verts;
-			network->ev_break = edges;
-			network->ev = tmp_edges;
-			network->break_dim = num_verts;
+			g->nv_break = num_verts;
+			g->nv = tmp_num_verts;
+			g->ev_break = edges;
+			g->ev = tmp_edges;
+			g->break_dim = num_verts;
 			break;
 											}
 		case EMBEDDED_BOUND: {
@@ -618,50 +618,49 @@ graph_t *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			edges = tmp_edges;
 			dual_edges = tmp_dual_edges;
 			vert_coords = tmp_vert_coords;
-			network->nv_break = num_verts;
-			network->nv = num_verts;
-			network->ev_break = edges;
-			network->ev = edges;
-			network->break_dim = num_verts + 2;
+			g->nv_break = num_verts;
+			g->nv = num_verts;
+			g->ev_break = edges;
+			g->ev = edges;
+			g->break_dim = num_verts + 2;
 		}
 	}
 
-	network->boundary = boundary;
-	network->num_bounds = num_bounds;
-	network->bound_inds = bound_inds;
-	network->bound_verts = bound_verts;
-	network->ne = num_edges;
-	network->dev = dual_edges;
-	network->vert_coords = vert_coords;
+	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;
 
-	network->vei = get_verts_to_edges_ind(
-			network->nv, network->ne, network->ev);
-	network->ve =
-			get_verts_to_edges(network->nv, network->ne,
-												 network->ev, network->vei);
+	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);
 
-	network->L = 1;
+	g->L = L;
 
-	network->edge_coords = get_edge_coords(
-			network->ne, network->vert_coords, network->ev);
+	g->ex = get_edge_coords(
+			g->ne, g->vx, g->ev);
 
 	free(tmp_tris);
 
-	network->dual_vert_coords = lattice;
-	network->dnv = num;
+	g->dvx = lattice;
+	g->dnv = num;
 
-	network->voltcurmat = gen_voltcurmat(network->ne,
-																			 network->break_dim,
-																			 network->ev_break, c);
+	g->voltcurmat = gen_voltcurmat(g->ne,
+																			 g->break_dim,
+																			 g->ev_break, c);
 
-	network->dvei =
-			get_verts_to_edges_ind(network->dnv, network->ne,
-														 network->dev);
-	network->dve = get_verts_to_edges(
-			network->dnv, network->ne, network->dev,
-			network->dvei);
+	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);
 
-	network->spanning_edges = get_spanning_edges(network->ne, network->ev_break, network->vert_coords, 0.5, &(network->num_spanning_edges));
+	g->spanning_edges = get_spanning_edges(g->ne, g->ev_break, g->vx, 0.5, &(g->num_spanning_edges));
 
-	return network;
+	return g;
 }
diff --git a/src/instance.c b/src/instance.c
deleted file mode 100644
index bb1ac8c..0000000
--- a/src/instance.c
+++ /dev/null
@@ -1,124 +0,0 @@
-
-#include "fracture.h"
-
-net_t *create_instance(graph_t *network, double inf, bool voltage_bound,
-											 bool startnow, cholmod_common *c) {
-	net_t *instance = (net_t *)calloc(1, sizeof(net_t));
-	assert(instance != NULL);
-
-	instance->graph = network;
-	instance->num_remaining_edges = network->ne;
-	instance->fuses = (bool *)calloc(network->ne, 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) {
-		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->nv + i] = -1;
-		}
-		((double *)instance->boundary_cond->x)[network->bound_verts[0]] = 1;
-	} else if (network->boundary == EMBEDDED_BOUND) {
-		// do nothing
-		for (unsigned int i = 0; i < network->bound_inds[1]; i++) {
-			((double *)instance->boundary_cond->x)[network->bound_verts[i]] =1;
-		}
-	} else {
-		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->nv] = 1;
-			((double *)instance->boundary_cond->x)[network->nv + 1] = -1;
-		}
-	}
-
-	if (network->boundary != TORUS_BOUND) instance->adjacency = gen_adjacency(instance, false, false, 0, c);
-	else instance->adjacency = gen_adjacency(instance, true, 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->graph->break_dim) *
-			sizeof(unsigned int));
-	instance->dual_marks = (unsigned int *)malloc(
-			(instance->graph->dnv) *
-			sizeof(unsigned int));
-	assert(instance->marks != NULL);
-
-	for (unsigned int i = 0;
-			 i < (instance->graph->break_dim);
-			 i++) {
-		instance->marks[i] = 1;
-	}
-	for (unsigned int i = 0;
-			 i < (instance->graph->dnv);
-			 i++) {
-		instance->dual_marks[i] = i+1;
-	}
-	instance->num_components = 1;
-
-	return instance;
-}
-
-void finish_instance(net_t *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);
-}
-
-net_t *copy_instance(const net_t *instance, cholmod_common *c) {
-	net_t *instance_copy = (net_t *)calloc(1, sizeof(net_t));
-	memcpy(instance_copy, instance, sizeof(net_t));
-
-	size_t fuses_size = (instance->graph)->ne * sizeof(bool);
-	instance_copy->fuses = (bool *)malloc(fuses_size);
-	memcpy(instance_copy->fuses, instance->fuses, fuses_size);
-
-	size_t marks_size =
-			(instance->graph->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->graph->dnv) *
-			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(net_t *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 net_t *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/net.c b/src/net.c
new file mode 100644
index 0000000..d090a9a
--- /dev/null
+++ b/src/net.c
@@ -0,0 +1,121 @@
+
+#include "fracture.h"
+
+double *get_thres(uint_t ne, double beta) {
+	assert(beta > 0);
+
+	double *thres = (double *)malloc(ne * sizeof(double));
+	assert(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 (uint_t i = 0; i < ne; i++) {
+		while ((thres[i] = exp(log(gsl_ran_flat(r, 0, 1)) / beta)) == 0.0);
+	}
+
+	gsl_rng_free(r);
+
+	return thres;
+}
+
+
+net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c) {
+	net_t *net = (net_t *)calloc(1, sizeof(net_t));
+	assert(net != NULL);
+
+	net->graph = g;
+	net->fuses = (bool *)calloc(g->ne, sizeof(bool));
+	assert(net->fuses != NULL);
+	net->thres = get_thres(g->ne, beta);
+	net->inf = inf;
+
+	net->voltage_bound = vb;
+	net->boundary_cond = bound_set(g, vb, notch_len, c);
+
+	if (g->boundary != TORUS_BOUND) net->adjacency = gen_adjacency(net, false, false, 0, c);
+	else net->adjacency = gen_adjacency(net, true, false, 0, c);
+
+	net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t));
+	net->dual_marks = (uint_t *)malloc((net->graph->dnv) * sizeof(uint_t));
+	assert(net->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);
+		net->factor = CHOL_F(analyze)(laplacian, c);
+		CHOL_F(factorize)(laplacian, net->factor, c);
+		CHOL_F(free_sparse)(&laplacian, c);
+	}
+
+	return net;
+}
+
+net_t *net_copy(const net_t *net, cholmod_common *c) {
+	net_t *net_copy = (net_t *)calloc(1, sizeof(net_t));
+	assert(net_copy != NULL);
+	memcpy(net_copy, net, sizeof(net_t));
+
+	size_t fuses_size = (net->graph)->ne * sizeof(bool);
+	net_copy->fuses = (bool *)malloc(fuses_size);
+	assert(net_copy->fuses != NULL);
+	memcpy(net_copy->fuses, net->fuses, fuses_size);
+
+	size_t thres_size = (net->graph)->ne * sizeof(double);
+	net_copy->thres = (double *)malloc(thres_size);
+	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);
+
+	net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c);
+	net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c);
+	net_copy->factor = CHOL_F(copy_factor)(net->factor, c);
+
+	return net_copy;
+}
+
+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_factor)(&(net->factor), c);
+	free(net->marks);
+	free(net->dual_marks);
+	free(net);
+}
+
+bool check_instance(const net_t *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/net_fracture.c b/src/net_fracture.c
new file mode 100644
index 0000000..b7fa61d
--- /dev/null
+++ b/src/net_fracture.c
@@ -0,0 +1,66 @@
+
+#include "fracture.h"
+
+uint_t get_next_broken(net_t *net, double *currents, double cutoff) {
+	uint_t max_pos = UINT_MAX;
+	double max_val = 0;
+
+	for (uint_t i = 0; i < net->graph->ne; i++) {
+		double current = fabs(currents[i]);
+		bool broken = net->fuses[i];
+
+		if (!broken && current > cutoff) {
+			double val = current / net->thres[i];
+
+			if (val > max_val) {
+				max_val = val;
+				max_pos = i;
+			}
+		}
+	}
+
+	if (max_pos == UINT_MAX) {
+		printf("GET_NEXT_BROKEN: currents is zero or NaN, no max_val found\n");
+		exit(EXIT_FAILURE);
+	}
+
+	return max_pos;
+}
+
+
+data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
+	data_t *data = alloc_break_data(net->graph->ne);
+
+	while (true) {
+		double *voltages = get_voltage(net, c);
+		double *currents = get_current_v(net, voltages, c);
+
+		double conductivity = get_conductivity(net, voltages, c);
+
+		if (conductivity < 1e-12 && net->graph->boundary == TORUS_BOUND) {
+			free(voltages);
+			free(currents);
+			break;
+		}
+
+		uint_t last_broke = get_next_broken(net, currents, cutoff);
+
+		update_break_data(data, last_broke, fabs(conductivity * (net->thres)[last_broke] / currents[last_broke]), conductivity);
+
+		free(voltages);
+		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_notch.c b/src/net_notch.c
new file mode 100644
index 0000000..608d0b5
--- /dev/null
+++ b/src/net_notch.c
@@ -0,0 +1,29 @@
+
+#include "fracture.h"
+
+void net_notch(net_t *net, double notch_len, cholmod_common *c) {
+	for (uint_t i = 0; i < net->graph->ne; i++) {
+		uint_t v1, v2;
+		double v1x, v1y, v2x, v2y, dx, dy;
+		bool crosses_center, not_wrapping, correct_length;
+
+		v1 = net->graph->ev[2 * i];
+		v2 = net->graph->ev[2 * i + 1];
+
+		v1x = net->graph->vx[2 * v1];
+		v1y = net->graph->vx[2 * v1 + 1];
+		v2x = net->graph->vx[2 * v2];
+		v2y = net->graph->vx[2 * v2 + 1];
+
+		dx = v1x - v2x;
+		dy = v1y - v2y;
+
+		crosses_center = (v1y >= 0.5 && v2y <= 0.5) || (v1y <= 0.5 && v2y >= 0.5);
+		not_wrapping = fabs(dy) < 0.5;
+		correct_length = v1x + dx / dy * (v1y - 0.5) <= notch_len;
+
+		if (crosses_center && not_wrapping && correct_length) {
+			break_edge(net, i, c);
+		}
+	}
+}
diff --git a/src/cracked_voronoi_fracture.c b/src/voro_fracture.c
index 93eb112..237192a 100644
--- a/src/cracked_voronoi_fracture.c
+++ b/src/voro_fracture.c
@@ -10,7 +10,7 @@ int main(int argc, char *argv[]) {
 	uint32_t N;
 	uint_t L;
 	double beta, inf, cutoff, crack_len;
-	bool include_breaking, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
+	bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
 			 save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_conductivity,
 			 save_damage, save_damage_field;
 	bound_t boundary;
@@ -30,7 +30,7 @@ int main(int argc, char *argv[]) {
 	inf = 1e10;
 	cutoff = 1e-9;
 	boundary = FREE_BOUND;
-	include_breaking = false;
+	save_data = false;
 	save_cluster_dist = false;
 	use_voltage_boundaries = false;
 	use_dual = false;
@@ -107,7 +107,7 @@ int main(int argc, char *argv[]) {
 				save_cluster_dist = true;
 				break;
 			case 'o':
-				include_breaking = true;
+				save_data = true;
 				break;
 			case 'N':
 				save_network = true;
@@ -138,12 +138,12 @@ int main(int argc, char *argv[]) {
 	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_%g.txt", boundc, boundc2, L, beta, crack_len);
-		break_out = fopen(break_filename, "a");
-		free(break_filename);
+	FILE *data_out;
+	if (save_data) {
+		char *data_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(data_filename, filename_len, "data_v_%c_%c_%u_%g_%g.txt", boundc, boundc2, L, beta, crack_len);
+		data_out = fopen(data_filename, "a");
+		free(data_filename);
 	}
 
 	uint_t voronoi_max_verts, c_dist_size, a_dist_size;
@@ -257,24 +257,17 @@ int main(int argc, char *argv[]) {
 	for (uint32_t i = 0; i < N; i++) {
 		printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1, N);
 
-		graph_t *graph = ini_voronoi_network(L, boundary, use_dual, genfunc_hyperuniform, &c);
-		net_t *perm_instance = create_instance(graph, inf, use_voltage_boundaries, false, &c);
-		if (boundary == EMBEDDED_BOUND) {
-			voronoi_bound_ini(perm_instance, L, crack_len);
-		}
-		gen_voro_crack(perm_instance, crack_len, &c);
-		finish_instance(perm_instance, &c);
-		double *fuse_thres = gen_fuse_thres(graph->ne, graph->edge_coords, beta, beta_scaling_flat);
-		net_t *instance = copy_instance(perm_instance, &c);
-		data_t *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
-		free_instance(instance, &c);
-		free(fuse_thres);
+		graph_t *g = ini_voro_graph(L, boundary, use_dual, genfunc_hyperuniform, &c);
+		net_t *net = net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c);
+		net_t *tmp_net = net_copy(net, &c);
+		data_t *data = net_fracture(tmp_net, &c, cutoff);
+		net_free(tmp_net, &c);
 
 		uint_t max_pos = 0;
 		double max_val = 0;
 
-		for (uint_t j = 0; j < breaking_data->num_broken; j++) {
-			double val = breaking_data->extern_field[j];
+		for (uint_t j = 0; j < data->num_broken; j++) {
+			double val = data->extern_field[j];
 
 			if (val > max_val) {
 				max_pos = j;
@@ -282,20 +275,18 @@ int main(int argc, char *argv[]) {
 			}
 		}
 
-		if (save_crit_stress) crit_stress[i] = breaking_data->extern_field[max_pos];
+		if (save_crit_stress) crit_stress[i] = data->extern_field[max_pos];
 
-		if (save_conductivity) conductivity[i] = breaking_data->conductivity[max_pos];
+		if (save_conductivity) conductivity[i] = data->conductivity[max_pos];
 
 		if (save_damage) damage[max_pos]++;
 
-		net_t *tmp_instance = copy_instance(perm_instance, &c);
-
 		uint_t av_size = 0;
 		double cur_val = 0;
 		for (uint_t j = 0; j < max_pos; j++) {
-			break_edge(tmp_instance, breaking_data->break_list[j], &c);
+			break_edge(net, data->break_list[j], &c);
 
-			double val = breaking_data->extern_field[j];
+			double val = data->extern_field[j];
 			if (save_cluster_dist) {
 				if (val < cur_val) {
 					av_size++;
@@ -310,20 +301,20 @@ int main(int argc, char *argv[]) {
 		}
 
 		if (save_stress_field || save_voltage_field) {
-			double *tmp_voltages = get_voltage(tmp_instance, &c);
+			double *tmp_voltages = get_voltage(net, &c);
 			if (save_voltage_field) {
-				for (uint_t j = 0; j < tmp_instance->graph->nv; j++) {
-					voltage_field[3 * voltage_pos] = tmp_instance->graph->vert_coords[2 * j];
-					voltage_field[3 * voltage_pos + 1] = tmp_instance->graph->vert_coords[2 * j + 1];
+				for (uint_t j = 0; j < g->nv; j++) {
+					voltage_field[3 * voltage_pos] = g->vx[2 * j];
+					voltage_field[3 * voltage_pos + 1] = g->vx[2 * j + 1];
 					voltage_field[3 * voltage_pos + 2] = tmp_voltages[j];
 					voltage_pos++;
 				}
 			}
 			if (save_stress_field) {
-				double *tmp_currents = get_current_v(tmp_instance, tmp_voltages, &c);
-				for (uint_t j = 0; j < tmp_instance->graph->ne; j++) {
-					stress_field[3 * stress_pos] = tmp_instance->graph->edge_coords[2 * j];
-					stress_field[3 * stress_pos + 1] = tmp_instance->graph->edge_coords[2 * j + 1];
+				double *tmp_currents = get_current_v(net, tmp_voltages, &c);
+				for (uint_t j = 0; j < g->ne; j++) {
+					stress_field[3 * stress_pos] = g->ex[2 * j];
+					stress_field[3 * stress_pos + 1] = g->ex[2 * j + 1];
 					stress_field[3 * stress_pos + 2] = tmp_currents[j];
 					stress_pos++;
 				}
@@ -334,8 +325,8 @@ int main(int argc, char *argv[]) {
 
 		if (save_damage_field) {
 			for (uint_t j = 0; j < max_pos; j++) {
-				damage_field[2 * damage_pos] = tmp_instance->graph->edge_coords[2 * breaking_data->break_list[j]];
-				damage_field[2 * damage_pos + 1] = tmp_instance->graph->edge_coords[2 * breaking_data->break_list[j] + 1];
+				damage_field[2 * damage_pos] = g->ex[2 * data->break_list[j]];
+				damage_field[2 * damage_pos + 1] = g->ex[2 * data->break_list[j] + 1];
 				damage_pos++;
 			}
 		}
@@ -343,11 +334,11 @@ int main(int argc, char *argv[]) {
 		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];
+				double sigma1 = data->extern_field[0];
+				double epsilon1 = sigma1 / data->conductivity[0];
 				for (uint_t j = 0; j < max_pos - 1; j++) {
-					double sigma2 = breaking_data->extern_field[j+1];
-					double epsilon2 = sigma2 / breaking_data->conductivity[j+1];
+					double sigma2 = data->extern_field[j+1];
+					double epsilon2 = sigma2 / data->conductivity[j+1];
 					if (epsilon2 > epsilon1) {
 						tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2;
 						sigma1 = sigma2; epsilon1 = epsilon2;
@@ -358,8 +349,8 @@ int main(int argc, char *argv[]) {
 		}
 
 		if (save_cluster_dist) {
-			uint_t *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c);
-			for (uint_t j = 0; j < tmp_instance->graph->dnv; j++) {
+			uint_t *tmp_cluster_dist = get_cluster_dist(net, &c);
+			for (uint_t j = 0; j < g->dnv; j++) {
 				cluster_size_dist[j] += tmp_cluster_dist[j];
 			}
 			free(tmp_cluster_dist);
@@ -367,45 +358,42 @@ int main(int argc, char *argv[]) {
 
 		if (save_network) {
 			FILE *net_out = fopen("network.txt", "w");
-			for (uint_t j = 0; j < graph->nv; j++) {
-				fprintf(net_out, "%f %f ", graph->vert_coords[2 * j],
-								tmp_instance->graph->vert_coords[2 * j + 1]);
+			for (uint_t j = 0; j < g->nv; j++) {
+				fprintf(net_out, "%f %f ", g->vx[2 * j],
+						g->vx[2 * j + 1]);
 			}
 			fprintf(net_out, "\n");
-			for (uint_t j = 0; j < tmp_instance->graph->ne; j++) {
-				fprintf(net_out, "%u %u ", tmp_instance->graph->ev[2 * j],
-								tmp_instance->graph->ev[2 * j + 1]);
+			for (uint_t j = 0; j < g->ne; j++) {
+				fprintf(net_out, "%u %u ", g->ev[2 * j], g->ev[2 * j + 1]);
 			}
 			fprintf(net_out, "\n");
-			for (uint_t j = 0; j < tmp_instance->graph->dnv; j++) {
-				fprintf(net_out, "%f %f ", tmp_instance->graph->dual_vert_coords[2 * j],
-								tmp_instance->graph->dual_vert_coords[2 * j + 1]);
+			for (uint_t j = 0; j < g->dnv; j++) {
+				fprintf(net_out, "%f %f ", g->dvx[2 * j],
+						g->dvx[2 * j + 1]);
 			}
 			fprintf(net_out, "\n");
-			for (uint_t j = 0; j < tmp_instance->graph->ne; j++) {
-				fprintf(net_out, "%u %u ", tmp_instance->graph->dev[2 * j],
-								tmp_instance->graph->dev[2 * j + 1]);
+			for (uint_t j = 0; j < g->ne; j++) {
+				fprintf(net_out, "%u %u ", g->dev[2 * j], g->dev[2 * j + 1]);
 			}
 			fprintf(net_out, "\n");
-			for (uint_t j = 0; j < tmp_instance->graph->ne; j++) {
-				fprintf(net_out, "%d ", tmp_instance->fuses[j]);
+			for (uint_t j = 0; j < g->ne; j++) {
+				fprintf(net_out, "%d ", net->fuses[j]);
 			}
 			fclose(net_out);
 		}
 
-		free_instance(tmp_instance, &c);
-		free_instance(perm_instance, &c);
-		free_net(graph, &c);
+		net_free(net, &c);
+		graph_free(g, &c);
 
-		if (include_breaking) {
-			for (uint_t 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]);
+		if (save_data) {
+			for (uint_t j = 0; j < data->num_broken; j++) {
+				fprintf(data_out, "%u %g %g ", data->break_list[j],
+								data->extern_field[j], data->conductivity[j]);
 			}
-			fprintf(break_out, "\n");
+			fprintf(data_out, "\n");
 		}
 
-		free_break_data(breaking_data);
+		free_break_data(data);
 	}
 
 	printf("\033[F\033[JFRACTURE: COMPLETE\n");
@@ -484,8 +472,8 @@ int main(int argc, char *argv[]) {
 		free(damage);
 	}
 
-	if (include_breaking) {
-		fclose(break_out);
+	if (save_data) {
+		fclose(data_out);
 	}
 
 	if (save_crit_stress) {
diff --git a/src/voronoi_bound_ini.c b/src/voronoi_bound_ini.c
deleted file mode 100644
index 0b65ef5..0000000
--- a/src/voronoi_bound_ini.c
+++ /dev/null
@@ -1,33 +0,0 @@
-
-#include "fracture.h"
-
-double th_p(double x, double y, double th) {
-	if (x >= 0 && y >= 0) return th;
-	if (x <  0 && y >= 0) return M_PI - th;
-	if (x <  0 && y <  0) return th - M_PI;
-	if (x >= 0 && y <  0) return -th;
-}
-
-double u_y(double x, double y) {
-	double r = sqrt(pow(x, 2) + pow(y, 2));
-	double th = th_p(x, y, atan(fabs(y / x)));
-
-	return sqrt(r) * sin(th / 2);
-}
-
-void voronoi_bound_ini(net_t *instance, uint_t L, double crack_len) {
-	double *bound = (double *)instance->boundary_cond->x;
-	for (uint_t i = 0; i < L / 2; i++) {
-		double x1, y1, x2, y2, x3, y3, x4, y4;
-		x1 = (2. * i + 1.) / L - crack_len; y1 = 0.5 - 1.;
-		x2 = (2. * i + 1.) / L - crack_len; y2 = 0.5 - 0.;
-		y3 = (2. * i + 1.) / L - 0.5; x3 = 0.5 - 1.;
-		y4 = (2. * i + 1.) / L - 0.5; x4 = 0.5 - 0.;
-
-		bound[i] = u_y(x1, y1);
-		bound[L / 2 + i] = u_y(x2, y2);
-		bound[L + i] = u_y(x3, y3);
-		bound[3 * L / 2 + i] = u_y(x4, y4);
-	}
-}
-