1 files changed, 66 insertions, 0 deletions
diff --git a/src/net_fracture.c b/src/net_fracture.c
new file mode 100644
index 0000000..b7fa61d
--- /dev/null
+++ b/src/net_fracture.c
@@ -0,0 +1,66 @@
+
+#include "fracture.h"
+
+uint_t get_next_broken(net_t *net, double *currents, double cutoff) {
+	uint_t max_pos = UINT_MAX;
+	double max_val = 0;
+
+	for (uint_t i = 0; i < net->graph->ne; i++) {
+		double current = fabs(currents[i]);
+		bool broken = net->fuses[i];
+
+		if (!broken && current > cutoff) {
+			double val = current / net->thres[i];
+
+			if (val > max_val) {
+				max_val = val;
+				max_pos = i;
+			}
+		}
+	}
+
+	if (max_pos == UINT_MAX) {
+		printf("GET_NEXT_BROKEN: currents is zero or NaN, no max_val found\n");
+		exit(EXIT_FAILURE);
+	}
+
+	return max_pos;
+}
+
+
+data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
+	data_t *data = alloc_break_data(net->graph->ne);
+
+	while (true) {
+		double *voltages = get_voltage(net, c);
+		double *currents = get_current_v(net, voltages, c);
+
+		double conductivity = get_conductivity(net, voltages, c);
+
+		if (conductivity < 1e-12 && net->graph->boundary == TORUS_BOUND) {
+			free(voltages);
+			free(currents);
+			break;
+		}
+
+		uint_t last_broke = get_next_broken(net, currents, cutoff);
+
+		update_break_data(data, last_broke, fabs(conductivity * (net->thres)[last_broke] / currents[last_broke]), conductivity);
+
+		free(voltages);
+		free(currents);
+
+		break_edge(net, last_broke, c);
+
+		if (net->num_components > 1 && net->graph->boundary == TORUS_BOUND) {
+			break;
+		}
+
+		if (net->marks[net->graph->nv] != net->marks[net->graph->nv + 1] && net->graph->boundary != TORUS_BOUND) {
+			break;
+		}
+	}
+
+	return data;
+}
+