removed percolation thing, changed definition of damage to make percolation exponents work

9 files changed, 45 insertions, 158 deletions

diff --git a/src/anal_process.c b/src/anal_process.c
index c434297..0ebec84 100644
--- a/src/anal_process.c
+++ b/src/anal_process.c
@@ -47,7 +47,7 @@ int main(int argc, char *argv[]) {
 		char *fn = argv[1 + i];
 		char *buff = (char *)malloc(20 * sizeof(char));
 		uint_t pos = 0; uint_t c = 0;
-		while (c < 6) {
+		while (c < 5) {
 			if (fn[pos] == '_') {
 				c++;
 			}
diff --git a/src/big_anal_process.c b/src/big_anal_process.c
index f255f20..0f7106f 100644
--- a/src/big_anal_process.c
+++ b/src/big_anal_process.c
@@ -63,7 +63,7 @@ int main(int argc, char *argv[]) {
 		char *fn = argv[1 + i];
 		char *buff = (char *)malloc(20 * sizeof(char));
 		uint_t pos = 0; uint_t c = 0;
-		while (c < 6) {
+		while (c < 5) {
 			if (fn[pos] == '_') {
 				c++;
 			}
diff --git a/src/break_edge.c b/src/break_edge.c
index 7fe7161..4e8bdcf 100644
--- a/src/break_edge.c
+++ b/src/break_edge.c
@@ -1,22 +1,14 @@
 
 #include "fracture.h"
 
-bool break_edge(net_t *instance, uint_t edge, cholmod_common *c, bool refactor) {
+bool break_edge(net_t *instance, uint_t edge, cholmod_common *c) {
 	instance->fuses[edge] = true;
 	instance->num_broken++;
 
 	if (instance->factor != NULL) {
-		if (refactor) {
-			cholmod_sparse *laplacian = gen_laplacian(instance, c, true);
-			CHOL_F(free_factor)(&instance->factor, c);
-			instance->factor = CHOL_F(analyze)(laplacian, c);
-			CHOL_F(factorize)(laplacian, instance->factor, c);
-			CHOL_F(free_sparse)(&laplacian, c);
-		} else {
-			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);
-		}
+		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);
 	}
 
 	uint_t v1, v2, s1, s2, dv1, dv2, ds1, ds2;
diff --git a/src/cen_anal_process.c b/src/cen_anal_process.c
index 7aad0e2..98f68ac 100644
--- a/src/cen_anal_process.c
+++ b/src/cen_anal_process.c
@@ -67,7 +67,7 @@ int main(int argc, char *argv[]) {
 		char *fn = argv[1 + i];
 		char *buff = (char *)malloc(20 * sizeof(char));
 		uint_t pos = 0; uint_t c = 0;
-		while (c < 6) {
+		while (c < 5) {
 			if (fn[pos] == '_') {
 				c++;
 			}
diff --git a/src/corr_test.c b/src/corr_test.c
index 5011ecc..4dc7efc 100644
--- a/src/corr_test.c
+++ b/src/corr_test.c
@@ -8,7 +8,7 @@ int main() {
 	unsigned int width = 64;
 
 	graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 128, false, &c);
-	net_t *instance = net_create(network, 1e-14, 3, 0, true, false, &c);
+	net_t *instance = net_create(network, 1e-14, 3, 0, true, &c);
 	net_fracture(instance, &c, 1e-10, 40);
 
 	double *corr = get_corr(instance, NULL, &c);
diff --git a/src/fracture.c b/src/fracture.c
index 5256eac..08961ed 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -8,11 +8,11 @@ int main(int argc, char *argv[]) {
 	// defining variables to be (potentially) set by command line flags
 	uint8_t filename_len;
 	uint32_t N;
-	uint_t L, refactor_every;
+	uint_t L;
 	double beta, inf, cutoff, crack_len;
-	bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, use_perc, save_network,
-			 save_crit_stress, save_stress_field, save_voltage_field, save_energy, save_conductivity,
-			 save_damage, save_damage_field, save_threshold;
+	bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
+			 save_crit_stress, save_energy, save_conductivity,
+			 save_damage, save_threshold;
 	bound_t boundary;
 	lattice_t lattice;
 
@@ -25,7 +25,6 @@ int main(int argc, char *argv[]) {
 	// set default values
 
 	N = 100;
-	refactor_every = UINT_MAX;
 	L = 16;
 	crack_len = 0.;
 	beta = .3;
@@ -39,14 +38,10 @@ int main(int argc, char *argv[]) {
 	use_dual = false;
 	save_network = false;
 	save_crit_stress = false;
-	save_stress_field = false;
-	save_voltage_field = false;
 	save_damage = false;
-	save_damage_field = false;
 	save_conductivity = false;
 	save_energy = false;
 	save_threshold = false;
-	use_perc = false;
 
 
 	uint8_t bound_i;
@@ -54,19 +49,15 @@ int main(int argc, char *argv[]) {
 	uint8_t lattice_i;
 	char lattice_c = 'v';
 	char dual_c = 'o';
-	char perc_c = 't';
 
 
 	// get commandline options
 
-	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDSvel:TER:p")) != -1) {
+	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TE")) != -1) {
 		switch (opt) {
 			case 'n':
 				N = atoi(optarg);
 				break;
-			case 'R':
-				refactor_every = atoi(optarg);
-				break;
 			case 'L':
 				L = atoi(optarg);
 				break;
@@ -76,10 +67,6 @@ int main(int argc, char *argv[]) {
 			case 'l':
 				crack_len = atof(optarg);
 				break;
-			case 'p':
-				use_perc = true;
-				perc_c = 'b';
-				break;
 			case 'B':
 				bound_i = atoi(optarg);
 				switch (bound_i) {
@@ -126,9 +113,6 @@ int main(int argc, char *argv[]) {
 			case 'd':
 				save_damage = true;
 				break;
-			case 'e':
-				save_damage_field = true;
-				break;
 			case 'V':
 				use_voltage_boundaries = true;
 				break;
@@ -148,12 +132,6 @@ int main(int argc, char *argv[]) {
 			case 's':
 				save_crit_stress = true;
 				break;
-			case 'S':
-				save_stress_field = true;
-				break;
-			case 'v':
-				save_voltage_field = true;
-				break;
 			case 'r':
 				save_conductivity = true;
 				break;
@@ -176,7 +154,7 @@ int main(int argc, char *argv[]) {
 	FILE *data_out;
 	if (save_data) {
 		char *data_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(data_filename, filename_len, "data_%c_%c_%c_%c_%c_%u_%g_%g.txt", perc_c, lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		snprintf(data_filename, filename_len, "data_%c_%c_%c_%c_%u_%g_%g.txt", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
 		data_out = fopen(data_filename, "a");
 		free(data_filename);
 	}
@@ -204,8 +182,8 @@ int main(int argc, char *argv[]) {
 
 		c_filename = (char *)malloc(filename_len * sizeof(char));
 		a_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(c_filename, filename_len, "cstr_%c_%c_%c_%c_%c_%d_%g_%g.dat", perc_c, lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%c_%d_%g_%g.dat", perc_c, lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		snprintf(c_filename, filename_len, "cstr_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
 
 		FILE *cluster_out = fopen(c_filename, "rb");
 		FILE *avalanche_out = fopen(a_filename, "rb");
@@ -225,24 +203,6 @@ int main(int argc, char *argv[]) {
 		crit_stress = (long double *)malloc(N * sizeof(long double));
 	}
 
-	double *stress_field;
-	unsigned int stress_pos = 0;
-	if (save_stress_field) {
-		stress_field = (double *)malloc(3 * N * max_verts * sizeof(double));
-	}
-
-	double *voltage_field;
-	unsigned int voltage_pos = 0;
-	if (save_voltage_field) {
-		voltage_field = (double *)malloc(3 * N * max_verts * sizeof(double));
-	}
-
-	double *damage_field;
-	unsigned int damage_pos = 0;
-	if (save_damage_field) {
-		damage_field = (double *)malloc(2 * N * max_verts * sizeof(double));
-	}
-
 	double *conductivity;
 	if (save_conductivity) {
 		conductivity = (double *)malloc(N * sizeof(double));
@@ -256,7 +216,7 @@ int main(int argc, char *argv[]) {
 				(uint32_t *)calloc(max_edges, sizeof(uint32_t));
 
 		d_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%c_%d_%g_%g.dat", perc_c, lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
 
 		FILE *damage_out = fopen(d_filename, "rb");
 
@@ -298,7 +258,7 @@ int main(int argc, char *argv[]) {
 		printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1, N);
 
 		graph_t *g = graph_create(lattice, boundary, L, use_dual, &c);
-		net_t *net = net_create(g, inf, beta, crack_len, use_voltage_boundaries, use_perc, &c);
+		net_t *net = net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c);
 		if (net->marks[g->nv] != net->marks[g->nv+1]) {
 			if (save_crit_stress) crit_stress[i] = 0;
 			if (save_conductivity) conductivity[i] = 0;
@@ -307,7 +267,7 @@ int main(int argc, char *argv[]) {
 			if (save_threshold) thresholds[i] = 0;
 		} else {
 		net_t *tmp_net = net_copy(net, &c);
-		data_t *data = net_fracture(tmp_net, &c, cutoff, refactor_every);
+		data_t *data = net_fracture(tmp_net, &c, cutoff);
 		net_free(tmp_net, &c);
 
 		uint_t max_pos = 0;
@@ -335,8 +295,7 @@ int main(int argc, char *argv[]) {
 		for (uint_t j = 0; j < max_pos; j++) {
 			uint_t next_broken = data->break_list[j];
 
-			bool refactor = ((j + 1) % refactor_every) == 0;
-			break_edge(net, next_broken, &c, refactor);
+			break_edge(net, next_broken, &c);
 
 			long double val = data->extern_field[j];
 			if (save_cluster_dist) {
@@ -362,37 +321,19 @@ int main(int argc, char *argv[]) {
 			}
 		}
 
-		if (save_damage) damage[net->num_broken]++;
-
-		if (save_stress_field || save_voltage_field) {
-			double *tmp_voltages = net_voltages(net, &c);
-			if (save_voltage_field) {
-				for (uint_t j = 0; j < g->nv; j++) {
-					voltage_field[3 * voltage_pos] = g->vx[2 * j];
-					voltage_field[3 * voltage_pos + 1] = g->vx[2 * j + 1];
-					voltage_field[3 * voltage_pos + 2] = tmp_voltages[j];
-					voltage_pos++;
-				}
-			}
-			if (save_stress_field) {
-				double *tmp_currents = net_currents(net, tmp_voltages, &c);
-				for (uint_t j = 0; j < g->ne; j++) {
-					stress_field[3 * stress_pos] = g->ex[2 * j];
-					stress_field[3 * stress_pos + 1] = g->ex[2 * j + 1];
-					stress_field[3 * stress_pos + 2] = tmp_currents[j];
-					stress_pos++;
-				}
-				free(tmp_currents);
-			}
-			free(tmp_voltages);
-		}
-
-		if (save_damage_field) {
-			for (uint_t j = 0; j < max_pos; j++) {
-				damage_field[2 * damage_pos] = g->ex[2 * data->break_list[j]];
-				damage_field[2 * damage_pos + 1] = g->ex[2 * data->break_list[j] + 1];
-				damage_pos++;
+		if (save_damage) {
+			uint_t would_break = 0;
+			double *tmp_voltage = net_voltages(net, &c);
+			double *tmp_current = net_currents(net, tmp_voltage, &c);
+			free(tmp_voltage);
+			for (uint_t j = 0; j < g->ne; j++) {
+				bool broken = net->fuses[j];
+				bool under_thres = net->thres[j] < net->thres[data->break_list[max_pos]];
+				bool zero_field = tmp_current[j] < inf;
+				if (!broken && under_thres && zero_field) would_break++;
 			}
+			printf("%u %g\n\n", would_break, ((double)(would_break + net->num_broken)) / pow(L,2));
+			damage[net->num_broken + would_break]++;
 		}
 
 		if (save_energy) {
@@ -485,39 +426,9 @@ int main(int argc, char *argv[]) {
 		free(avalanche_size_dist);
 	}
 
-	if (save_voltage_field) {
-		char *vfld_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(vfld_filename, filename_len, "vfld_%c_%c_%c_%c_%c_%d_%g_%g.dat", perc_c, lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		FILE *vfld_file = fopen(vfld_filename, "ab");
-		fwrite(voltage_field, sizeof(double), 3 * voltage_pos, vfld_file);
-		fclose(vfld_file);
-		free(vfld_filename);
-		free(voltage_field);
-	}
-
-	if (save_stress_field) {
-		char *cfld_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(cfld_filename, filename_len, "cfld_%c_%c_%c_%c_%c_%d_%g_%g.dat", perc_c, lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		FILE *cfld_file = fopen(cfld_filename, "ab");
-		fwrite(stress_field, sizeof(double), 3 * stress_pos, cfld_file);
-		fclose(cfld_file);
-		free(cfld_filename);
-		free(stress_field);
-	}
-
-	if (save_damage_field) {
-		char *dfld_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(dfld_filename, filename_len, "dfld_%c_%c_%c_%c_%c_%d_%g_%g.dat", perc_c, lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		FILE *dfld_file = fopen(dfld_filename, "ab");
-		fwrite(damage_field, sizeof(double), 2 * damage_pos, dfld_file);
-		fclose(dfld_file);
-		free(dfld_filename);
-		free(damage_field);
-	}
-
 	if (save_conductivity) {
 		char *cond_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(cond_filename, filename_len, "cond_%c_%c_%c_%c_%c_%d_%g_%g.dat", perc_c, lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		snprintf(cond_filename, filename_len, "cond_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
 		FILE *cond_file = fopen(cond_filename, "ab");
 		fwrite(conductivity, sizeof(double), N, cond_file);
 		fclose(cond_file);
@@ -527,7 +438,7 @@ int main(int argc, char *argv[]) {
 
 	if (save_energy) {
 		char *tough_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%c_%d_%g_%g.dat", perc_c, lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
 		FILE *tough_file = fopen(tough_filename, "ab");
 		fwrite(energy, sizeof(long double), N, tough_file);
 		fclose(tough_file);
@@ -537,7 +448,7 @@ int main(int argc, char *argv[]) {
 
 	if (save_threshold) {
 		char *thres_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%c_%d_%g_%g.dat", perc_c, lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
 		FILE *thres_file = fopen(thres_filename, "ab");
 		fwrite(thresholds, sizeof(long double), N, thres_file);
 		fclose(thres_file);
@@ -559,7 +470,7 @@ int main(int argc, char *argv[]) {
 
 	if (save_crit_stress) {
 		char *str_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(str_filename, filename_len, "strs_%c_%c_%c_%c_%c_%d_%g_%g.dat", perc_c, lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		snprintf(str_filename, filename_len, "strs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
 		FILE *str_file = fopen(str_filename, "ab");
 		fwrite(crit_stress, sizeof(long double), N, str_file);
 		fclose(str_file);
diff --git a/src/fracture.h b/src/fracture.h
index e02bcef..f4f7334 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -118,7 +118,7 @@ double dual_vert_to_coord(uint_t width, bool periodic, uint_t vert,
 void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c);
 
 void net_notch(net_t *net, double notch_len, cholmod_common *c);
-data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refactor_every);
+data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff);
 double *net_voltages(const net_t *net, cholmod_common *c);
 double *net_currents(const net_t *net, const double *voltages, cholmod_common *c);
 double net_conductivity(const net_t *net, const double *voltages);
@@ -143,13 +143,13 @@ graph_t *ini_square_network(uint_t width, bound_t boundary, bool side_bounds,
 void graph_free(graph_t *network, cholmod_common *c);
 void net_free(net_t *instance, cholmod_common *c);
 
-net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, bool perc, cholmod_common *c);
+net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c);
 
 graph_t *ini_voro_graph(uint_t L, bound_t boundary, bool use_dual,
 													double *(*genfunc)(uint_t, bound_t, gsl_rng *, uint_t *),
 													cholmod_common *c);
 
-bool break_edge(net_t *instance, uint_t edge, cholmod_common *c, bool refactor);
+bool break_edge(net_t *instance, uint_t edge, cholmod_common *c);
 
 void finish_instance(net_t *instance, cholmod_common *c);
 
diff --git a/src/net.c b/src/net.c
index 1be5613..21a2139 100644
--- a/src/net.c
+++ b/src/net.c
@@ -1,7 +1,7 @@
 
 #include "fracture.h"
 
-long double *get_thres(uint_t ne, double beta, bool perc) {
+long double *get_thres(uint_t ne, double beta) {
 	long double *thres = (long double *)malloc(ne * sizeof(long double));
 	assert(thres != NULL);
 
@@ -9,8 +9,7 @@ long double *get_thres(uint_t ne, double beta, bool perc) {
 	gsl_rng_set(r, rand_seed());
 
 	for (uint_t i = 0; i < ne; i++) {
-		if (perc) thres[i] = 1;
-		else thres[i] = rand_dist_pow(r, beta);
+		thres[i] = rand_dist_pow(r, beta);
 	}
 
 	gsl_rng_free(r);
@@ -45,19 +44,7 @@ void net_notch(net_t *net, double notch_len, cholmod_common *c) {
 	}
 }
 
-void net_perc(net_t *net, double p, cholmod_common *c) {
-	gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN);
-	gsl_rng_set(r, rand_seed());
-
-	for (uint_t i = 0; i < net->graph->ne; i++) {
-		double x = gsl_rng_uniform_pos(r);
-		if (!net->fuses[i] && p > x) break_edge(net, i, c, false);
-	}
-
-	gsl_rng_free(r);
-}
-
-net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, bool perc, cholmod_common *c) {
+net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c) {
 	net_t *net = (net_t *)calloc(1, sizeof(net_t));
 	assert(net != NULL);
 
@@ -65,7 +52,7 @@ net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, b
 	net->num_broken = 0;
 	net->fuses = (bool *)calloc(g->ne, sizeof(bool));
 	assert(net->fuses != NULL);
-	net->thres = get_thres(g->ne, beta, perc);
+	net->thres = get_thres(g->ne, beta);
 	net->inf = inf;
 
 	net->voltage_bound = vb;
@@ -90,8 +77,6 @@ net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, b
 
 	net_notch(net, notch_len, c);
 
-	if (perc) net_perc(net, beta, c);
-
 	{
 		cholmod_sparse *laplacian = gen_laplacian(net, c, true);
 		net->factor = CHOL_F(analyze)(laplacian, c);
diff --git a/src/net_fracture.c b/src/net_fracture.c
index b5ab524..71e8c4b 100644
--- a/src/net_fracture.c
+++ b/src/net_fracture.c
@@ -27,7 +27,7 @@ uint_t get_next_broken(net_t *net, double *currents, double cutoff) {
 	return max_pos;
 }
 
-data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refactor_every) {
+data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
 	data_t *data = data_create(net->graph->ne);
 
 	uint_t n = 0;
@@ -59,8 +59,7 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff, uint_t refact
 		free(voltages);
 		free(currents);
 
-		bool refactor = (n % refactor_every) == 0;
-		break_edge(net, last_broke, c, refactor);
+		break_edge(net, last_broke, c);
 
 		if (net->num_components > 1 && net->graph->boundary == TORUS_BOUND) {
 			break;