From 2bb0740b68fdb62d45adc00204b3990ca42ade77 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jkentdobias@g.hmc.edu>
Date: Mon, 22 Aug 2016 10:11:14 -0400
Subject: started repo again without all the data files gunking the works

---
 src/ini_network.c | 620 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 620 insertions(+)
 create mode 100644 src/ini_network.c

(limited to 'src/ini_network.c')

diff --git a/src/ini_network.c b/src/ini_network.c
new file mode 100644
index 0000000..0434d98
--- /dev/null
+++ b/src/ini_network.c
@@ -0,0 +1,620 @@
+
+#include "fracture.h"
+
+unsigned int *get_spanning_edges(unsigned int num_edges, unsigned int *edges_to_verts, double *vert_coords, double cut, unsigned int *n) {
+	unsigned int *spanning_edges = (unsigned int *)malloc(num_edges * sizeof(unsigned int));
+	(*n) = 0;
+	for (unsigned int i = 0; i < num_edges; i++) {
+		unsigned int v1, v2;
+		v1 = edges_to_verts[2 * i];
+		v2 = edges_to_verts[2 * i + 1];
+		double v1y, v2y;
+		v1y = vert_coords[2 * v1 + 1];
+		v2y = vert_coords[2 * v2 + 1];
+		if ((fabs(v1y - v2y) < 0.25) && ((v1y < cut && v2y > cut) || (v1y > cut && v2y < cut))) {
+			spanning_edges[*n] = i;
+			(*n)++;
+		}
+	}
+	return spanning_edges;
+}
+
+double *get_edge_coords(unsigned int num_edges, double *vert_coords,
+												unsigned int *edges_to_verts) {
+	double *output = (double *)malloc(2 * num_edges * sizeof(double));
+
+#pragma omp parallel for
+	for (unsigned int i = 0; i < num_edges; i++) {
+		unsigned int v1, v2;
+		double v1x, v1y, v2x, v2y, dx, dy;
+		v1 = edges_to_verts[2 * i];
+		v2 = edges_to_verts[2 * i + 1];
+		v1x = vert_coords[2 * v1];
+		v1y = vert_coords[2 * v1 + 1];
+		v2x = vert_coords[2 * v2];
+		v2y = vert_coords[2 * v2 + 1];
+		dx = v1x - v2x;
+		dy = v1y - v2y;
+		if (fabs(dx) > 0.5) {
+			if (dx > 0) {
+				v2x = v2x + 1;
+			} else {
+				v1x = v1x + 1;
+			}
+		}
+		if (fabs(dy) > 0.5) {
+			if (dy > 0) {
+				v2y = v2y + 1;
+			} else {
+				v1y = v1y + 1;
+			}
+		}
+		output[2 * i] = (v1x + v2x) / 2;
+		output[2 * i + 1] = (v1y + v2y) / 2;
+	}
+
+	return output;
+}
+
+unsigned int *get_verts_to_edges_ind(unsigned int num_verts,
+																		 unsigned int num_edges,
+																		 const unsigned int *edges_to_verts) {
+	unsigned int *output =
+			(unsigned int *)calloc(num_verts + 1, sizeof(unsigned int));
+	assert(output != NULL);
+
+	for (unsigned int i = 0; i < 2 * num_edges; i++) {
+		if (edges_to_verts[i] < num_verts) {
+			output[edges_to_verts[i] + 1]++;
+		}
+	}
+
+	for (unsigned int i = 0; i < num_verts; i++) {
+		output[i + 1] += output[i];
+	}
+
+	return output;
+}
+
+unsigned int *get_verts_to_edges(unsigned int num_verts, unsigned int num_edges,
+																 const unsigned int *edges_to_verts,
+																 const unsigned int *verts_to_edges_ind) {
+	unsigned int *output = (unsigned int *)calloc(verts_to_edges_ind[num_verts],
+																								sizeof(unsigned int));
+	unsigned int *counts =
+			(unsigned int *)calloc(num_verts, sizeof(unsigned int));
+	for (int i = 0; i < 2 * num_edges; i++) {
+		if (edges_to_verts[i] < num_verts) {
+			output[verts_to_edges_ind[edges_to_verts[i]] +
+						 counts[edges_to_verts[i]]] = i / 2;
+			counts[edges_to_verts[i]]++;
+		}
+	}
+
+	free(counts);
+
+	return output;
+}
+
+fnet *ini_square_network(unsigned int width, boundary_type boundary, bool side_bounds,
+												 cholmod_common *c) {
+	fnet *network = (fnet *)calloc(1, sizeof(fnet));
+
+	network->boundary = boundary;
+	bool periodic = (boundary == CYLINDER_BOUND) || (boundary == TORUS_BOUND) ? true : false;
+
+	network->num_edges = pow(width, 2);
+	if (boundary == CYLINDER_BOUND) {
+		assert(width % 2 == 0);
+		assert(!side_bounds);
+		network->num_verts = (width / 2) * (width + 1);
+		network->num_verts_break = (width / 2) * (width + 1);
+		network->num_dual_verts = (width / 2 + 1) * (width / 2) + pow(width / 2, 2);
+		network->num_bounds = 2;
+		network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) *
+																								 sizeof(unsigned int));
+		network->bound_inds[0] = 0;
+		network->bound_inds[1] = width / 2;
+		network->bound_inds[2] = width;
+		network->bound_verts = (unsigned int *)calloc(width, sizeof(unsigned int));
+		network->break_dim = network->num_verts + network->num_bounds;
+	} else if (boundary == FREE_BOUND) {
+		network->num_verts = 2 * ((width + 1) / 2) * (width / 2 + 1);
+		network->num_verts_break = 2 * ((width + 1) / 2) * (width / 2 + 1);
+		network->num_dual_verts = pow(width / 2 + 1, 2) + pow((width + 1) / 2, 2);
+		if (side_bounds)
+			network->num_bounds = 4;
+		else
+			network->num_bounds = 2;
+		network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) *
+																								 sizeof(unsigned int));
+		for (unsigned int i = 0; i < network->num_bounds + 1; i++) {
+			network->bound_inds[i] = i * ((width + 1) / 2);
+		}
+		network->bound_verts = (unsigned int *)calloc(
+				network->num_bounds * ((width + 1) / 2), sizeof(unsigned int));
+		if (side_bounds) {
+			for (unsigned int i = 0; i < (width + 1) / 2; i++) {
+				network->bound_verts[2 * ((width + 1) / 2) + i] =
+						(width + 1) / 2 + i * (width + 1);
+				if (width % 2) {
+					network->bound_verts[3 * ((width + 1) / 2) + i] =
+							(width + 1) / 2 + i * (width + 1) - 1;
+				} else {
+					network->bound_verts[3 * ((width + 1) / 2) + i] =
+							(width + 1) / 2 + i * (width + 1) + width / 2;
+				}
+			}
+		}
+		network->break_dim = network->num_verts + network->num_bounds;
+	} else if (boundary == TORUS_BOUND) {
+		network->num_verts = (width / 2) * (width + 1) - (width / 2);
+		network->num_verts_break = (width / 2) * (width + 1);
+		network->num_dual_verts = (width / 2 + 1) * (width / 2) + pow(width / 2, 2) - (width / 2);
+		network->num_bounds = 1;
+		network->bound_inds = (unsigned int *)malloc((network->num_bounds + 1) *
+																								 sizeof(unsigned int));
+		network->bound_inds[0] = 0;
+		network->bound_inds[1] = width / 2;
+		network->bound_verts = (unsigned int *)calloc(width / 2, sizeof(unsigned int));
+		network->break_dim = network->num_verts_break;
+	}
+	if (boundary != TORUS_BOUND) {
+		for (unsigned int i = 0; i < (width + 1) / 2; i++) {
+			network->bound_verts[i] = i;
+			network->bound_verts[(width + 1) / 2 + i] = network->num_verts - 1 - i;
+		} 
+	} else {
+		for (unsigned int i = 0; i < width / 2; i++) {
+			network->bound_verts[i] = i;
+		}
+	}
+	network->edges_to_verts_break =
+			(unsigned int *)calloc(2 * network->num_edges, sizeof(unsigned int));
+	network->edges_to_verts =
+			(unsigned int *)calloc(2 * network->num_edges, sizeof(unsigned int));
+	for (unsigned int i = 0; i < network->num_edges; i++) {
+		network->edges_to_verts_break[2 * i] = edge_to_verts(width, periodic, i, 1);
+		network->edges_to_verts_break[2 * i + 1] = edge_to_verts(width, periodic, i, 0);
+		network->edges_to_verts[2 * i] = network->edges_to_verts_break[2 * i] % network->num_verts;
+		network->edges_to_verts[2 * i + 1] = network->edges_to_verts_break[2 * i + 1] % network->num_verts;
+	}
+	network->verts_to_edges_ind =
+			(unsigned int *)calloc(network->num_verts + 1, sizeof(unsigned int));
+	network->verts_to_edges_ind[0] = 0;
+	unsigned int pos1 = 0;
+	for (unsigned int i = 0; i < network->num_verts; i++) {
+		bool in_bound = false;
+		for (unsigned int j = 0; j < network->num_bounds; j++) {
+			for (unsigned int k = 0;
+					 k < network->bound_inds[j + 1] - network->bound_inds[j]; k++) {
+				if (i == network->bound_verts[network->bound_inds[j] + k]) {
+					in_bound = true;
+					break;
+				}
+			}
+		}
+		if (in_bound)
+			pos1 += 2;
+		else
+			pos1 += 4;
+
+		network->verts_to_edges_ind[i + 1] = pos1;
+	}
+	network->verts_to_edges = (unsigned int *)calloc(
+			network->verts_to_edges_ind[network->num_verts], sizeof(unsigned int));
+	unsigned int *vert_counts =
+			(unsigned int *)calloc(network->num_verts, sizeof(unsigned int));
+	for (unsigned int i = 0; i < network->num_edges; i++) {
+		unsigned int v0 = network->edges_to_verts[2 * i];
+		unsigned int v1 = network->edges_to_verts[2 * i + 1];
+		network->verts_to_edges[network->verts_to_edges_ind[v0] + vert_counts[v0]] =
+				i;
+		network->verts_to_edges[network->verts_to_edges_ind[v1] + vert_counts[v1]] =
+				i;
+		vert_counts[v0]++;
+		vert_counts[v1]++;
+	}
+	free(vert_counts);
+
+	network->vert_coords =
+			(double *)malloc(2 * network->num_verts * sizeof(double));
+	for (unsigned int i = 0; i < network->num_verts; i++) {
+		if (!periodic) {
+		network->vert_coords[2 * i] = ((double)((2 * i + 1) / (width + 1)))/width;
+		network->vert_coords[2 * i + 1] = ((double)((2 * i + 1) % (width + 1)))/width;
+		} else {
+		network->vert_coords[2 * i] = ((double)((2 * i + 1) / (width)))/width;
+		network->vert_coords[2 * i + 1] = ((double)((((2 * i + 1) / (width)+1) % 2) + ((2 * i) % (width))))/width;
+		}
+	}
+
+	network->L = width;
+	network->edge_coords = get_edge_coords(
+			network->num_edges, network->vert_coords, network->edges_to_verts);
+
+	network->dual_edges_to_verts =
+			(unsigned int *)malloc(2 * network->num_edges * sizeof(unsigned int));
+	for (unsigned int i = 0; i < network->num_edges; i++) {
+		network->dual_edges_to_verts[2 * i] =
+				dual_edge_to_verts(width, periodic, i, 0) % network->num_dual_verts;
+		network->dual_edges_to_verts[2 * i + 1] =
+				dual_edge_to_verts(width, periodic, i, 1) % network->num_dual_verts;
+	}
+	network->dual_vert_coords =
+			(double *)malloc(2 * network->num_dual_verts * sizeof(double));
+	for (unsigned int i = 0; i < network->num_dual_verts; i++) {
+		network->dual_vert_coords[2 * i] =
+				2*dual_vert_to_coord(width, periodic, i, 0);
+		network->dual_vert_coords[2 * i + 1] =
+				2*dual_vert_to_coord(width, periodic, i, 1);
+	}
+
+	network->voltcurmat = gen_voltcurmat(network->num_edges,
+																			 network->break_dim,
+																			 network->edges_to_verts_break, c);
+
+	network->dual_verts_to_edges_ind =
+			get_verts_to_edges_ind(network->num_dual_verts, network->num_edges,
+														 network->dual_edges_to_verts);
+	network->dual_verts_to_edges = get_verts_to_edges(
+			network->num_dual_verts, network->num_edges, network->dual_edges_to_verts,
+			network->dual_verts_to_edges_ind);
+	network->spanning_edges = get_spanning_edges(network->num_edges, network->edges_to_verts_break, network->vert_coords, 0.51, &(network->num_spanning_edges));
+
+	return network;
+}
+
+unsigned int *get_voro_dual_edges(unsigned int num_edges,
+																	unsigned int num_verts, unsigned int *edges,
+																	unsigned int *triangles) {
+	unsigned int *dual_edges =
+			(unsigned int *)malloc(2 * num_edges * sizeof(unsigned int));
+#pragma omp parallel for
+	for (unsigned int i = 0; i < num_edges; i++) {
+		unsigned int v1, v2;
+		v1 = edges[2 * i];
+		v2 = edges[2 * i + 1];
+		if (v1 < num_verts && v2 < num_verts) {
+			bool found_match = false;
+			for (unsigned int j = 0; j < 3; j++) {
+				for (unsigned int k = 0; k < 3; k++) {
+					unsigned int t11, t12, t21, t22;
+					t11 = triangles[3 * v1 + j];
+					t12 = triangles[3 * v1 + ((j + 1) % 3)];
+					t21 = triangles[3 * v2 + k];
+					t22 = triangles[3 * v2 + ((k + 1) % 3)];
+					if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) {
+						dual_edges[2 * i] = t11 < t12 ? t11 : t12;
+						dual_edges[2 * i + 1] = t11 < t12 ? t12 : t11;
+						found_match = true;
+						break;
+					}
+				}
+				if (found_match)
+					break;
+			}
+		}
+	}
+
+	return dual_edges;
+}
+
+fnet *ini_voronoi_network(unsigned int L, boundary_type boundary,
+													double *(*genfunc)(unsigned int, gsl_rng *, unsigned int *),
+													cholmod_common *c) {
+	fnet *network = (fnet *)calloc(1, sizeof(fnet));
+
+	// generate the dual lattice
+	double *lattice;
+	unsigned int num;
+	{
+		gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+		FILE *rf = fopen("/dev/urandom", "r");
+		unsigned long int seed;
+		fread(&seed, sizeof(unsigned long int), 1, rf);
+		fclose(rf);
+		gsl_rng_set(r, seed);
+		lattice = genfunc(L, r, &num);
+		gsl_rng_free(r);
+	}
+
+	// retrieve a periodic voronoi tesselation of the lattice
+	intptr_t *vout = run_voronoi(num, lattice, 0, 1, 0, 1);
+
+	unsigned int tmp_num_verts = ((unsigned int *)vout[0])[0];
+	unsigned int tmp_num_edges = ((unsigned int *)vout[0])[1];
+	double *tmp_vert_coords = (double *)vout[2];
+	unsigned int *tmp_edges = (unsigned int *)vout[3];
+	unsigned int *tmp_tris = (unsigned int *)vout[5];
+
+	free((void *)vout[0]);
+	free((void *)vout[1]);
+	free((void *)vout[4]);
+	free(vout);
+
+	// get dual edges of the fully periodic graph
+	unsigned int *tmp_dual_edges =
+			get_voro_dual_edges(tmp_num_edges, tmp_num_verts, tmp_edges, tmp_tris);
+
+	// prune the edges of the lattice and assign boundary vertices based on the
+	// desired boundary conditions
+	unsigned int num_bounds;
+	unsigned num_verts;
+	double *vert_coords;
+	unsigned int *bound_inds;
+	unsigned int *bound_verts;
+	unsigned int num_edges;
+	unsigned int *edges;
+	unsigned int *dual_edges;
+	switch (boundary) {
+		case FREE_BOUND: {
+			num_bounds = 4;
+			bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int));
+			bound_inds[0] = 0;
+			vert_coords = tmp_vert_coords;
+			num_verts = tmp_num_verts;
+			num_edges = 0;
+			edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
+			dual_edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
+			unsigned int num_t, num_b, num_l, num_r;
+			bool *bound_t, *bound_b, *bound_l, *bound_r;
+			num_t = 0; num_b = 0; num_l = 0; num_r = 0;
+			bound_t = (bool *)calloc(num_verts, sizeof(bool));
+			bound_b = (bool *)calloc(num_verts, sizeof(bool));
+			bound_l = (bool *)calloc(num_verts, sizeof(bool));
+			bound_r = (bool *)calloc(num_verts, sizeof(bool));
+			for (unsigned int i = 0; i < tmp_num_edges; i++) {
+				unsigned int v1, v2;
+				double v1x, v1y, v2x, v2y, dx, dy;
+				v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
+				v1x = vert_coords[2 * v1]; v1y = vert_coords[2 * v1 + 1];
+				v2x = vert_coords[2 * v2]; v2y = vert_coords[2 * v2 + 1];
+				dx = v1x - v2x; dy = v1y - v2y;
+				if (fabs(dy) > 0.5) {
+					if (dy > 0) {
+						if (!bound_t[v1] && !bound_l[v1] && !bound_r[v1]) {
+							bound_t[v1] = true; num_t++;
+						}
+						if (!bound_b[v2] && !bound_l[v2] && !bound_r[v2]) {
+							bound_b[v2] = true; num_b++;
+						}
+					} else {
+						if (!bound_t[v2] && !bound_l[v2] && !bound_r[v2]) {
+							bound_t[v2] = true; num_t++;
+						}
+						if (!bound_b[v1] && !bound_l[v1] && !bound_r[v1]) {
+							bound_b[v1] = true; num_b++;
+						}
+					}
+				} else if (fabs(dx) > 0.5) {
+					if (dx > 0) {
+						if (!bound_r[v1] && !bound_t[v1] && !bound_b[v1]) {
+							bound_r[v1] = true; num_r++;
+						}
+						if (!bound_l[v2] && !bound_t[v2] && !bound_b[v2]) {
+							bound_l[v2] = true; num_l++;
+						}
+					} else {
+						if (!bound_r[v2] && !bound_t[v2] && !bound_b[v2]) {
+							bound_r[v2] = true; num_r++;
+						}
+						if (!bound_l[v1] && !bound_t[v1] && !bound_b[v1]) {
+							bound_l[v1] = true; num_l++;
+						}
+					}
+				} else {
+					edges[2 * num_edges] = v1 < v2 ? v1 : v2;
+					edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
+					unsigned int d1 = tmp_dual_edges[2 * i];
+					unsigned int d2 = tmp_dual_edges[2 * i + 1];
+					dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2;
+					dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1;
+					num_edges++;
+				}
+			}
+			bound_verts = malloc((num_t + num_b + num_l + num_r) * sizeof(unsigned int));
+			bound_inds[1] = num_t; bound_inds[2] = num_t + num_b;
+			bound_inds[3] = num_l + num_t + num_b;
+			bound_inds[4] = num_t + num_b + num_r + num_l;
+			unsigned int pos_t, pos_b, pos_l, pos_r;
+			pos_t = 0; pos_b = 0; pos_l = 0; pos_r = 0;
+			for (unsigned int i = 0; i < num_verts; i++) {
+				if (bound_t[i]) {
+					bound_verts[pos_t] = i; pos_t++;
+				} else if (bound_b[i]) {
+					bound_verts[num_t + pos_b] = i; pos_b++;
+				} else if (bound_l[i]) {
+					bound_verts[num_t + num_b + pos_l] = i; pos_l++;
+				} else if (bound_r[i]) {
+					bound_verts[num_t + num_b + num_l + pos_r] = i; pos_r++;
+				}
+			}
+			free(bound_l); free(bound_r); free(bound_t); free(bound_b);
+			free(tmp_edges);
+			free(tmp_dual_edges);
+			network->edges_to_verts_break = edges;
+			network->edges_to_verts = edges;
+			network->num_verts_break = num_verts;
+			network->num_verts = num_verts;
+			network->break_dim = num_verts + num_bounds;
+			break;
+										 }
+		case CYLINDER_BOUND: {
+			num_bounds = 2;
+			bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int));
+			bound_inds[0] = 0;
+			vert_coords = tmp_vert_coords;
+			num_verts = tmp_num_verts;
+			num_edges = 0;
+			edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
+			dual_edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
+			unsigned int num_t, num_b;
+			bool *bound_t, *bound_b;
+			num_t = 0; num_b = 0;
+			bound_t = (bool *)calloc(num_verts, sizeof(bool));
+			bound_b = (bool *)calloc(num_verts, sizeof(bool));
+			for (unsigned int i = 0; i < tmp_num_edges; i++) {
+				unsigned int v1, v2;
+				double v1x, v1y, v2x, v2y, dx, dy;
+				v1 = tmp_edges[2 * i]; v2 = tmp_edges[2 * i + 1];
+				v1y = vert_coords[2 * v1 + 1]; v2y = vert_coords[2 * v2 + 1];
+				dy = v1y - v2y;
+				if (fabs(dy) > 0.5) {
+					if (dy > 0) {
+						if (!bound_t[v1]) {
+							bound_t[v1] = true; num_t++;
+						}
+						if (!bound_b[v2]) {
+							bound_b[v2] = true; num_b++;
+						}
+					} else {
+						if (!bound_t[v2]) {
+							bound_t[v2] = true; num_t++;
+						}
+						if (!bound_b[v1]) {
+							bound_b[v1] = true; num_b++;
+						}
+					}
+				} else {
+					edges[2 * num_edges] = v1 < v2 ? v1 : v2;
+					edges[2 * num_edges + 1] = v1 < v2 ? v2 : v1;
+					unsigned int d1 = tmp_dual_edges[2 * i];
+					unsigned int d2 = tmp_dual_edges[2 * i + 1];
+					dual_edges[2 * num_edges] = d1 < d2 ? d1 : d2;
+					dual_edges[2 * num_edges + 1] = d1 < d2 ? d2 : d1;
+					num_edges++;
+				}
+			}
+			bound_verts = malloc((num_t + num_b) * sizeof(unsigned int));
+			bound_inds[1] = num_t; bound_inds[2] = num_t + num_b;
+			unsigned int pos_t, pos_b;
+			pos_t = 0; pos_b = 0;
+			for (unsigned int i = 0; i < num_verts; i++) {
+				if (bound_t[i]) {
+					bound_verts[pos_t] = i; pos_t++;
+				} else if (bound_b[i]) {
+					bound_verts[num_t + pos_b] = i; pos_b++;
+				}
+			}
+			free(bound_t); free(bound_b);
+			free(tmp_edges);
+			free(tmp_dual_edges);
+			network->edges_to_verts_break = edges;
+			network->edges_to_verts = edges;
+			network->num_verts_break = num_verts;
+			network->num_verts = num_verts;
+			network->break_dim = num_verts + num_bounds;
+			break;
+												 }
+		case TORUS_BOUND: {
+			num_bounds = 1;
+			bound_inds = (unsigned int *)malloc((1 + num_bounds) * sizeof(unsigned int));
+			bound_inds[0] = 0;
+			num_edges = tmp_num_edges;
+			edges = (unsigned int *)malloc(2* num_edges*sizeof(unsigned int));
+			for (unsigned int i = 0; i < num_edges; i++) {
+				edges[2*i] = tmp_edges[2*i];
+				edges[2*i+1] = tmp_edges[2*i+1];
+			}
+			dual_edges = tmp_dual_edges;
+			bool *bound_t = (bool *)calloc(tmp_num_verts, sizeof(bool));
+			int *edge_change = (int *)calloc(num_edges, sizeof(int));
+			unsigned int num_t = 0;
+			for (unsigned int i = 0; i < num_edges; i++) {
+				unsigned int v1, v2;
+				double v1x, v1y, v2x, v2y, dx, dy;
+				v1 = edges[2 * i]; v2 = edges[2 * i + 1];
+				v1x = tmp_vert_coords[2 * v1]; v1y = tmp_vert_coords[2 * v1 + 1];
+				v2x = tmp_vert_coords[2 * v2]; v2y = tmp_vert_coords[2 * v2 + 1];
+				dy = v1y - v2y;
+				if (fabs(dy) > 0.5) {
+					if (dy > 0) {
+						edge_change[i] = 1;
+						if (!bound_t[v1]) {
+							bound_t[v1] = true; num_t++;
+						}
+					} else {
+						edge_change[i] = 2;
+						if (!bound_t[v2]) {
+							bound_t[v2] = true; num_t++;
+						}
+					}
+				}
+			}
+			num_verts = tmp_num_verts + num_t;
+			vert_coords = (double *)malloc(2 * num_verts * sizeof(double));
+			bound_verts = malloc(num_t * sizeof(unsigned int));
+			bound_inds[1] = num_t;
+			unsigned int pos_t = 0;
+			for (unsigned int i = 0; i < tmp_num_verts; i++) {
+				vert_coords[2*i] = tmp_vert_coords[2*i];
+				vert_coords[2*i+1] = tmp_vert_coords[2*i+1];
+				if (bound_t[i]) {
+					bound_verts[pos_t] = i;
+					vert_coords[2*(tmp_num_verts + pos_t)] = tmp_vert_coords[2*i];
+					vert_coords[2*(tmp_num_verts + pos_t)+1] = tmp_vert_coords[2*i+1];
+					pos_t++;
+				}
+			}
+			for (unsigned int i = 0; i < num_edges; i++) {
+				if (edge_change[i]) {
+					for (unsigned int j = 0; j < num_t; j++) {
+						if (edges[2*i+(edge_change[i]-1)] == bound_verts[j]) {
+							edges[2*i+(edge_change[i]-1)] = tmp_num_verts + j;
+							break;
+						}
+					}
+				}
+			}
+			free(tmp_vert_coords);
+			free(bound_t);
+			free(edge_change);
+			network->num_verts_break = num_verts;
+			network->num_verts = tmp_num_verts;
+			network->edges_to_verts_break = edges;
+			network->edges_to_verts = tmp_edges;
+			network->break_dim = num_verts;
+			break;
+											}
+	}
+
+	network->boundary = boundary;
+	network->num_bounds = num_bounds;
+	network->bound_inds = bound_inds;
+	network->bound_verts = bound_verts;
+	network->num_edges = num_edges;
+	network->dual_edges_to_verts = dual_edges;
+	network->vert_coords = vert_coords;
+
+	network->verts_to_edges_ind = get_verts_to_edges_ind(
+			network->num_verts, network->num_edges, network->edges_to_verts);
+	network->verts_to_edges =
+			get_verts_to_edges(network->num_verts, network->num_edges,
+												 network->edges_to_verts, network->verts_to_edges_ind);
+
+	network->L = 1;
+
+	network->edge_coords = get_edge_coords(
+			network->num_edges, network->vert_coords, network->edges_to_verts);
+
+	free(tmp_tris);
+
+	network->dual_vert_coords = lattice;
+	network->num_dual_verts = num;
+
+	network->voltcurmat = gen_voltcurmat(network->num_edges,
+																			 network->break_dim,
+																			 network->edges_to_verts_break, c);
+
+	network->dual_verts_to_edges_ind =
+			get_verts_to_edges_ind(network->num_dual_verts, network->num_edges,
+														 network->dual_edges_to_verts);
+	network->dual_verts_to_edges = get_verts_to_edges(
+			network->num_dual_verts, network->num_edges, network->dual_edges_to_verts,
+			network->dual_verts_to_edges_ind);
+
+	network->spanning_edges = get_spanning_edges(network->num_edges, network->edges_to_verts_break, network->vert_coords, 0.5, &(network->num_spanning_edges));
+
+	return network;
+}
-- 
cgit v1.2.3-70-g09d2