added preliminary support for gathering statistics at crack growth, not critical stress

1 files changed, 40 insertions, 38 deletions

diff --git a/src/break_edge.c b/src/break_edge.c
index 54eaf34..01570d8 100644
--- a/src/break_edge.c
+++ b/src/break_edge.c
@@ -1,46 +1,50 @@
 
 #include "fracture.h"
 
-bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) {
+bool break_edge(net_t *instance, uint_t edge, cholmod_common *c) {
 	instance->fuses[edge] = true;
 
-	unsigned int v1 = instance->graph->ev_break[2 * edge];
-	unsigned int v2 = instance->graph->ev_break[2 * edge + 1];
+	if (instance->factor != NULL) {
+		uint_t w1 = instance->graph->ev_break[2 * edge];
+		uint_t w2 = instance->graph->ev_break[2 * edge + 1];
+		factor_update(instance->factor, w1, w2, c);
+	}
 
-	if (instance->factor != NULL) factor_update(instance->factor, v1, v2, c);
+	uint_t v1, v2, s1, s2, dv1, dv2, ds1, ds2;
 
-	if (instance->graph->boundary != TORUS_BOUND) {
-		unsigned int w1 = instance->graph->ev[2 * edge];
-		unsigned int w2 = instance->graph->ev[2 * edge + 1];
+	v1 = instance->graph->ev[2 * edge];
+	v2 = instance->graph->ev[2 * edge + 1];
+	dv1 = instance->graph->dev[2 * edge];
+	dv2 = instance->graph->dev[2 * edge + 1];
 
-		unsigned int tw1 = w1 > w2 ? w1 : w2;
-		unsigned int tw2 = w1 > w2 ? w2 : w1;
+	s1 = v1 > v2 ? v1 : v2;
+	s2 = v1 > v2 ? v2 : v1;
+	ds1 = dv1 > dv2 ? dv1 : dv2;
+	ds2 = dv1 > dv2 ? dv2 : dv1;
 
+	{
 		int_t *lap_p = (int_t *)instance->adjacency->p;
 		int_t *lap_i = (int_t *)instance->adjacency->i;
 		double *lap_x = (double *)instance->adjacency->x;
 
-		for (int i = 0; i < lap_p[tw1 + 1] - lap_p[tw1]; i++) {
-			if (lap_i[lap_p[tw1] + i] == tw2)
-				lap_x[lap_p[tw1] + i] = 0;
+		for (int i = 0; i < lap_p[s1 + 1] - lap_p[s1]; i++) {
+			if (lap_i[lap_p[s1] + i] == s2)
+				lap_x[lap_p[s1] + i] = 0;
 		}
-		for (int i = 0; i < lap_p[tw2 + 1] - lap_p[tw2]; i++) {
-			if (lap_i[lap_p[tw2] + i] == tw1)
-				lap_x[lap_p[tw2] + i] = 0;
+		for (int i = 0; i < lap_p[s2 + 1] - lap_p[s2]; i++) {
+			if (lap_i[lap_p[s2] + i] == s1)
+				lap_x[lap_p[s2] + i] = 0;
 		}
+	}
 
-		int old_num_components = instance->num_components;
+	int_t old_num_components = instance->num_components;
 
-		instance->num_components = update_components(
-				instance->adjacency, instance->marks, old_num_components, (int)tw1,
-				(int)tw2, 0);
-	}
+	instance->num_components = update_components(
+			instance->adjacency, instance->marks, old_num_components, s1, s2, 0);
 
-	if (instance->graph->boundary == TORUS_BOUND) {
-		unsigned int dw1 = instance->graph->dev[2 * edge];
-		unsigned int dw2 = instance->graph->dev[2 * edge + 1];
 
-		if (instance->dual_marks[dw1] == instance->dual_marks[dw2]) {
+	if (instance->graph->boundary == TORUS_BOUND) {
+		if (instance->dual_marks[dv1] == instance->dual_marks[dv2]) {
 			int **cycles = (int **)malloc(4*instance->graph->ne * sizeof(int *));
 			unsigned int num_cycles = find_cycles(instance->graph->ne, instance->fuses, instance->graph->dev, instance->graph->dvei, instance->graph->dve, cycles);
 
@@ -76,26 +80,24 @@ bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) {
 			}
 			free(cycles);
 		}
+	}
 
-		unsigned int tw1 = dw1 > dw2 ? dw1 : dw2;
-		unsigned int tw2 = dw1 > dw2 ? dw2 : dw1;
-
-		int_t *lap_p = (int_t *)instance->adjacency->p;
-		int_t *lap_i = (int_t *)instance->adjacency->i;
-		double *lap_x = (double *)instance->adjacency->x;
+	{
+		int_t *lap_p = (int_t *)instance->dual_adjacency->p;
+		int_t *lap_i = (int_t *)instance->dual_adjacency->i;
+		double *lap_x = (double *)instance->dual_adjacency->x;
 
-		for (int i = 0; i < lap_p[tw1 + 1] - lap_p[tw1]; i++) {
-			if (lap_i[lap_p[tw1] + i] == tw2)
-				lap_x[lap_p[tw1] + i] = 1;
+		for (int i = 0; i < lap_p[ds1 + 1] - lap_p[ds1]; i++) {
+			if (lap_i[lap_p[ds1] + i] == ds2)
+				lap_x[lap_p[ds1] + i] = 1;
 		}
-		for (int i = 0; i < lap_p[tw2 + 1] - lap_p[tw2]; i++) {
-			if (lap_i[lap_p[tw2] + i] == tw1)
-				lap_x[lap_p[tw2] + i] = 1;
+		for (int i = 0; i < lap_p[ds2 + 1] - lap_p[ds2]; i++) {
+			if (lap_i[lap_p[ds2] + i] == ds1)
+				lap_x[lap_p[ds2] + i] = 1;
 		}
-
-		set_connected(instance->adjacency, instance->dual_marks, dw1, instance->dual_marks[dw2], -1, 0);
 	}
 
+	set_connected(instance->dual_adjacency, instance->dual_marks, dv1, instance->dual_marks[dv2], -1, 0);
 
 	return true;
 }