1 files changed, 99 insertions, 0 deletions
diff --git a/src/get_current.c b/src/get_current.c
new file mode 100644
index 0000000..5360054
--- /dev/null
+++ b/src/get_current.c
@@ -0,0 +1,99 @@
+
+#include "fracture.h"
+
+double *get_voltage(const finst *instance, cholmod_common *c) {
+	cholmod_dense *b = instance->boundary_cond;
+	cholmod_factor *factor = instance->factor;
+
+	cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
+
+	if (((double *)x->x)[0] != ((double *)x->x)[0]) {
+		printf("ERROR: GET_VOLTAGE FAILED\n\n");
+		CHOL_F(free_dense)(&x, c);
+		return NULL;
+	}
+
+	double *field = (double *)x->x;
+	x->x = NULL;
+
+	CHOL_F(free_dense)(&x, c);
+
+
+	return field;
+}
+
+double *get_current(const finst *instance, cholmod_common *c) {
+	unsigned int num_edges = instance->network->num_edges;
+	unsigned int num_verts = instance->network->num_verts_break;
+	unsigned int num_gverts = instance->network->break_dim;
+	unsigned int num_bounds = instance->network->num_bounds;
+	cholmod_sparse *voltcurmat = instance->network->voltcurmat;
+
+	double *voltages = get_voltage(instance, c);
+	if (voltages == NULL) {
+		return NULL;
+	}
+	cholmod_dense *x = CHOL_F(allocate_dense)(
+			num_gverts, 1, num_gverts, CHOLMOD_REAL, c);
+	double *tmp_x = x->x;
+	x->x = voltages;
+
+	cholmod_dense *y =
+			CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c);
+
+	double alpha[2] = {1, 0};
+	double beta[2] = {0, 0};
+	CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+
+	double *field = (double *)malloc(num_edges * sizeof(double));
+
+	for (int i = 0; i < num_edges; i++) {
+		if (instance->fuses[i])
+			field[i] = 0;
+		else
+			field[i] = ((double *)y->x)[i];
+	}
+
+	x->x = tmp_x;
+	free(voltages);
+	CHOL_F(free_dense)(&x, c);
+	CHOL_F(free_dense)(&y, c);
+
+	return field;
+}
+
+
+double *get_current_v(const finst *instance, double *voltages, cholmod_common *c) {
+	unsigned int num_edges = instance->network->num_edges;
+	unsigned int num_verts = instance->network->num_verts_break;
+	unsigned int num_gverts = instance->network->break_dim;
+	unsigned int num_bounds = instance->network->num_bounds;
+	cholmod_sparse *voltcurmat = instance->network->voltcurmat;
+
+	cholmod_dense *x = CHOL_F(allocate_dense)(
+			num_gverts, 1, num_gverts, CHOLMOD_REAL, c);
+	double *tmp_x = x->x;
+	x->x = voltages;
+
+	cholmod_dense *y =
+			CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c);
+
+	double alpha[2] = {1, 0};
+	double beta[2] = {0, 0};
+	CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+
+	double *field = (double *)malloc(num_edges * sizeof(double));
+
+	for (int i = 0; i < num_edges; i++) {
+		if (instance->fuses[i])
+			field[i] = 0;
+		else
+			field[i] = ((double *)y->x)[i];
+	}
+
+	x->x = tmp_x;
+	CHOL_F(free_dense)(&x, c);
+	CHOL_F(free_dense)(&y, c);
+
+	return field;
+}