From 1e1fdfc2e3892667bccaf317a01defd8832041c7 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 16 Jan 2017 01:31:10 -0500
Subject: fixed voltage and torus conditions, current and free boundaries and
 broken right now

---
 src/bound_set.c | 62 ++++++++++++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 50 insertions(+), 12 deletions(-)

(limited to 'src/bound_set.c')

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;
 			}
-- 
cgit v1.2.3-70-g09d2