#include "fracture.h"

double *get_voltage(const net_t *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]) {
		for (uint_t i = 0; i < instance->graph->ne; i++) { 
			printf("%d ", instance->fuses[i]);
		}
		printf("\n");
		printf("GET_VOLTAGE: value is NaN\n");
		exit(EXIT_FAILURE);
	}

	double *field = (double *)x->x;
	x->x = NULL;

	CHOL_F(free_dense)(&x, c);


	return field;
}

double *get_current(const net_t *instance, cholmod_common *c) {
	unsigned int num_edges = instance->graph->ne;
	unsigned int num_verts = instance->graph->nv_break;
	unsigned int num_gverts = instance->graph->break_dim;
	unsigned int num_bounds = instance->graph->num_bounds;
	cholmod_sparse *voltcurmat = instance->graph->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 net_t *instance, double *voltages, cholmod_common *c) {
	unsigned int num_edges = instance->graph->ne;
	unsigned int num_verts = instance->graph->nv_break;
	unsigned int num_gverts = instance->graph->break_dim;
	unsigned int num_bounds = instance->graph->num_bounds;
	cholmod_sparse *voltcurmat = instance->graph->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;
}