1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
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;
}
|
