From 901b9f16494f37890be17ef4bb66e6efc6873340 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 10 Feb 2017 12:18:11 -0500 Subject: changed code to rely on jst --- lib/net_currents.c | 51 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 lib/net_currents.c (limited to 'lib/net_currents.c') diff --git a/lib/net_currents.c b/lib/net_currents.c new file mode 100644 index 0000000..a078336 --- /dev/null +++ b/lib/net_currents.c @@ -0,0 +1,51 @@ + +#include "fracture.h" + +double *net_currents(const net_t *net, const double *voltages, cholmod_common *c) { + uint_t ne = net->graph->ne; + uint_t dim = net->graph->nv; + cholmod_sparse *voltcurmat = net->voltcurmat; + + cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c); + + double *tmp_x = x->x; + x->x = (void *)voltages; + + cholmod_dense *y = CHOL_F(allocate_dense)(ne, 1, ne, CHOLMOD_REAL, c); + + double alpha[2] = {1, 0}; + double beta[2] = {0, 0}; + CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); + + double *currents = (double *)malloc(ne * sizeof(double)); + + for (int i = 0; i < ne; i++) { + if (net->fuses[i]) { + currents[i] = 0; + } else { + currents[i] = ((double *)y->x)[i]; + } + } + + if (net->graph->boundary == TORUS_BOUND) { + for (uint_t i = 0; i < net->graph->bi[1]; i++) { + uint_t e = net->graph->b[i]; + uint_t v1 = net->graph->ev[2 * e]; + uint_t v2 = net->graph->ev[2 * e + 1]; + double v1y = net->graph->vx[2 * v1 + 1]; + double v2y = net->graph->vx[2 * v2 + 1]; + + if (v1y > v2y) { + currents[e] += 1; + } else { + currents[e] -= 1; + } + } + } + + x->x = tmp_x; + CHOL_F(free_dense)(&x, c); + CHOL_F(free_dense)(&y, c); + + return currents; +} -- cgit v1.2.3-70-g09d2