1 files changed, 137 insertions, 0 deletions
diff --git a/lib/gen_laplacian.c b/lib/gen_laplacian.c
new file mode 100644
index 0000000..1cc93a4
--- /dev/null
+++ b/lib/gen_laplacian.c
@@ -0,0 +1,137 @@
+
+#include "fracture.h"
+
+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 = nre;
+
+	cholmod_triplet *t = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, 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 a = 0;
+
+	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];
+
+			uint_t s1 = v1 < v2 ? v1 : v2;
+			uint_t s2 = v1 < v2 ? v2 : v1;
+
+			ri[a] = s2;
+			ci[a] = s1;
+			ai[a] = 1;
+
+			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 *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 nnz = nv;
+
+	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;
+	double *acoo = (double *)temp_m->x;
+
+	temp_m->nnz = nnz;
+
+	for (uint_t i = 0; i < nv; i++) {
+		rowind[i] = i;
+		colind[i] = i;
+		acoo[i] = 0;
+	}
+
+	cholmod_sparse *adjacency = gen_adjacency(net, false, true, true, c);
+
+	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];
+
+			acoo[v1]++;
+			acoo[v2]++;
+		}
+	}
+
+	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];
+
+				if (g->bq[v0] || g->bq[v1]) {
+					acoo[i]++;
+				}
+			}
+		}
+	} else {
+		acoo[0]++;
+	}
+
+	for (uint_t i = 0; i < nv; i++) {
+		if (acoo[i] == 0) acoo[i]++;
+	}
+
+	//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));
+
+	double alpha[2] = {1, 0};
+	double beta[2] = {-1, 0};
+	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_triplet)(&temp_m, c);
+
+	return laplacian;
+}