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

1 files changed, 72 insertions, 181 deletions

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++;
-		}
-	}
-
-	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++;
-			}
+			a++;
 		}
 	}
 
-
 	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;
+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 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;
+	uint_t nnz = nv;
 
-	uint_t nnz = num_gverts + 2 * num_gedges;
-
-	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]]++;
+	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];
 
-					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]]++;
+				if (g->bq[v0] || g->bq[v1]) {
+					acoo[i]++;
 				}
 			}
-		} 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;
-			}
-
-			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;
-				}
-			}
-
-			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;