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

4 files changed, 79 insertions, 49 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;
 }
diff --git a/src/fracture.c b/src/fracture.c
index f38568f..718fc05 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -10,7 +10,7 @@ int main(int argc, char *argv[]) {
 	uint32_t N;
 	uint_t L;
 	double beta, inf, cutoff, crack_len;
-	bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
+	bool crack_growth_crit, save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
 			 save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_conductivity,
 			 save_damage, save_damage_field;
 	bound_t boundary;
@@ -32,6 +32,7 @@ int main(int argc, char *argv[]) {
 	cutoff = 1e-9;
 	boundary = FREE_BOUND;
 	lattice = VORONOI_LATTICE;
+	crack_growth_crit = false;
 	save_data = false;
 	save_cluster_dist = false;
 	use_voltage_boundaries = false;
@@ -54,7 +55,7 @@ int main(int argc, char *argv[]) {
 
 	// get commandline options
 
-	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:")) != -1) {
+	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:g")) != -1) {
 		switch (opt) {
 			case 'n':
 				N = atoi(optarg);
@@ -111,6 +112,9 @@ int main(int argc, char *argv[]) {
 						exit(EXIT_FAILURE);
 				}
 				break;
+			case 'g':
+				crack_growth_crit = true;
+				break;
 			case 'd':
 				save_damage = true;
 				break;
@@ -294,16 +298,31 @@ int main(int argc, char *argv[]) {
 			}
 		}
 
-		if (save_crit_stress) crit_stress[i] = data->extern_field[max_pos];
-
-		if (save_conductivity) conductivity[i] = data->conductivity[max_pos];
-
-		if (save_damage) damage[max_pos]++;
-
 		uint_t av_size = 0;
 		double cur_val = 0;
+		uint_t notch_edge;
+
+		if (crack_growth_crit) {
+			for (uint_t i = 0; i < net->graph->ne; i++) {
+				if (net->fuses[i]) {
+					notch_edge = i;
+					break;
+				}
+			}
+		}
+
 		for (uint_t j = 0; j < max_pos; j++) {
-			break_edge(net, data->break_list[j], &c);
+			uint_t next_broken = data->break_list[j];
+
+			bool attached_to_notch = net->dual_marks[next_broken] == net->dual_marks[notch_edge];
+			bool grew_forward = net->graph->ex[2 * next_broken] > crack_len;
+
+			if (crack_growth_crit && attached_to_notch && grew_forward) {
+				max_pos = j;
+				break;
+			}
+
+			break_edge(net, next_broken, &c);
 
 			double val = data->extern_field[j];
 			if (save_cluster_dist) {
@@ -319,6 +338,12 @@ int main(int argc, char *argv[]) {
 			}
 		}
 
+		if (save_crit_stress) crit_stress[i] = data->extern_field[max_pos];
+
+		if (save_conductivity) conductivity[i] = data->conductivity[max_pos];
+
+		if (save_damage) damage[max_pos]++;
+
 		if (save_stress_field || save_voltage_field) {
 			double *tmp_voltages = net_voltages(net, &c);
 			if (save_voltage_field) {
diff --git a/src/fracture.h b/src/fracture.h
index 5e1b3e4..f0a22e8 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -81,6 +81,7 @@ typedef struct {
 	bool voltage_bound;
 	unsigned int num_components;
 	cholmod_sparse *adjacency;
+	cholmod_sparse *dual_adjacency;
 	bool debug_stop;
 } net_t;
 
diff --git a/src/net.c b/src/net.c
index 9f6965b..f3aeda9 100644
--- a/src/net.c
+++ b/src/net.c
@@ -57,8 +57,8 @@ net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, b
 	net->voltage_bound = vb;
 	net->boundary_cond = bound_set(g, vb, notch_len, c);
 
-	if (g->boundary != TORUS_BOUND) net->adjacency = gen_adjacency(net, false, false, 0, c);
-	else net->adjacency = gen_adjacency(net, true, false, 0, c);
+	net->adjacency = gen_adjacency(net, false, false, 0, c);
+	net->dual_adjacency = gen_adjacency(net, true, false, 0, c);
 
 	net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t));
 	assert(net->marks != NULL);
@@ -112,6 +112,7 @@ net_t *net_copy(const net_t *net, cholmod_common *c) {
 	memcpy(net_copy->dual_marks, net->dual_marks, dual_marks_size);
 
 	net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c);
+	net_copy->dual_adjacency = CHOL_F(copy_sparse)(net->dual_adjacency, c);
 	net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c);
 	net_copy->factor = CHOL_F(copy_factor)(net->factor, c);
 
@@ -123,6 +124,7 @@ void net_free(net_t *net, cholmod_common *c) {
 	free(net->thres);
 	CHOL_F(free_dense)(&(net->boundary_cond), c);
 	CHOL_F(free_sparse)(&(net->adjacency), c);
+	CHOL_F(free_sparse)(&(net->dual_adjacency), c);
 	CHOL_F(free_factor)(&(net->factor), c);
 	free(net->marks);
 	free(net->dual_marks);