From a66e8e2d7a1e9d5eb51f4954435513f6ecb8383c Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sat, 14 Jan 2017 11:59:39 -0500
Subject: tried to fix voltage boundary conditions

---
 src/fracture.c      |  8 ++++++++
 src/gen_laplacian.c | 28 ++++++++++++++++++++++++++--
 src/graph_create.c  |  4 ++--
 3 files changed, 36 insertions(+), 4 deletions(-)

diff --git a/src/fracture.c b/src/fracture.c
index bb7701b..af9cd0f 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -325,6 +325,14 @@ 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->bound_inds[1]; j++) {
+			fprintf(testout, "%g ", tmp_voltage[g->bound_verts[j]]-tmp_voltage[g->nv+j]);
+		}
+		fclose(testout);
+		free(tmp_voltage);
+
 		if (save_damage) {
 			uint_t would_break = 0;
 			double *tmp_voltage = net_voltages(net, &c);
diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c
index 42134f0..04f1f9b 100644
--- a/src/gen_laplacian.c
+++ b/src/gen_laplacian.c
@@ -102,11 +102,11 @@ cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
 
 	uint_t num_gedges;
 	if (voltage_bound) {
-		num_gedges = 0;
+		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];
+	if (network->boundary == TORUS_BOUND) num_gedges = 3 * bound_inds[1];
 	uint_t num_gverts = network->break_dim;
 
 	uint_t nnz = num_gverts + 2 * num_gedges;
@@ -147,6 +147,24 @@ cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
 			acoo[num_gverts + 2*i+1] = -1;
 			acoo[vv]++;
 			acoo[network->nv + i]++;
+
+			rowind[num_gverts + 2 * network->bound_inds[1] + 2*i] = network->nv + i;
+			colind[num_gverts + 2 * network->bound_inds[1] + 2*i] = network->nv + network->bound_inds[1];
+			acoo[num_gverts + 2 * network->bound_inds[1] + 2*i] = -1;
+			rowind[num_gverts + 2 * network->bound_inds[1] + 2*i+1] = network->nv + network->bound_inds[1];
+			colind[num_gverts + 2 * network->bound_inds[1] + 2*i+1] = network->nv + i;
+			acoo[num_gverts + 2 * network->bound_inds[1] + 2*i+1] = -1;
+			acoo[network->nv + i]++;
+			acoo[network->nv + network->bound_inds[1]]++;
+
+			rowind[num_gverts + 4 * network->bound_inds[1] + 2*i] = vv;
+			colind[num_gverts + 4 * network->bound_inds[1] + 2*i] = network->nv + network->bound_inds[1];
+			acoo[num_gverts + 4 * network->bound_inds[1] + 2*i] = -1;
+			rowind[num_gverts + 4 * network->bound_inds[1] + 2*i+1] = network->nv + network->bound_inds[1];
+			colind[num_gverts + 4 * network->bound_inds[1] + 2*i+1] = vv;
+			acoo[num_gverts + 4 * network->bound_inds[1] + 2*i+1] = -1;
+			acoo[vv]++;
+			acoo[network->nv + network->bound_inds[1]]++;
 		}
 		acoo[bound_verts[0]]++;
 	}
@@ -158,6 +176,12 @@ cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
 			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[num_gverts + 2 * (bound_inds[i] + j)] = network->nv;
+					colind[num_gverts + 2 * (bound_inds[i] + j)] = bound_verts[bound_inds[i] + j];
+					rowind[num_gverts + 2 * (bound_inds[i] + j) + 1] = bound_verts[bound_inds[i] + j];
+					colind[num_gverts + 2 * (bound_inds[i] + j) + 1] = network->nv;
+					acoo[network->nv]++;
+					acoo[bound_verts[bound_inds[i] + j]]++;
 				}
 			}
 		} else {
diff --git a/src/graph_create.c b/src/graph_create.c
index b071c49..c1f556c 100644
--- a/src/graph_create.c
+++ b/src/graph_create.c
@@ -158,7 +158,7 @@ graph_t *ini_square_network(uint_t width, bound_t boundary, bool side_bounds,
 		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;
+		network->break_dim = network->nv_break + 1;
 	}
 	if (boundary != TORUS_BOUND) {
 		for (uint_t i = 0; i < (width + 1) / 2; i++) {
@@ -598,7 +598,7 @@ graph_t *ini_voro_graph(uint_t L, bound_t boundary, bool use_dual,
 			g->nv = tmp_num_verts;
 			g->ev_break = edges;
 			g->ev = tmp_edges;
-			g->break_dim = num_verts;
+			g->break_dim = num_verts + 1;
 			break;
 											}
 		case EMBEDDED_BOUND: {
-- 
cgit v1.2.3-70-g09d2