added support for controlling disorder by random initial bond breaking

6 files changed, 66 insertions, 34 deletions

diff --git a/src/break_edge.c b/src/break_edge.c
index 538b88d..7fe7161 100644
--- a/src/break_edge.c
+++ b/src/break_edge.c
@@ -3,6 +3,7 @@
 
 bool break_edge(net_t *instance, uint_t edge, cholmod_common *c, bool refactor) {
 	instance->fuses[edge] = true;
+	instance->num_broken++;
 
 	if (instance->factor != NULL) {
 		if (refactor) {
diff --git a/src/corr_test.c b/src/corr_test.c
index 4dc7efc..5011ecc 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, &c);
+	net_t *instance = net_create(network, 1e-14, 3, 0, true, false, &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 59bd424..6050182 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, refactor_every;
 	double beta, inf, cutoff, crack_len;
-	bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
+	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;
 	bound_t boundary;
@@ -44,8 +44,9 @@ int main(int argc, char *argv[]) {
 	save_damage = false;
 	save_damage_field = false;
 	save_conductivity = false;
-	save_energy = true;
+	save_energy = false;
 	save_threshold = false;
+	use_perc = false;
 
 
 	uint8_t bound_i;
@@ -53,11 +54,12 @@ 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:")) != -1) {
+	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDSvel:TER:p")) != -1) {
 		switch (opt) {
 			case 'n':
 				N = atoi(optarg);
@@ -74,6 +76,10 @@ 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) {
@@ -170,7 +176,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_%u_%g_%g.txt", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		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);
 		data_out = fopen(data_filename, "a");
 		free(data_filename);
 	}
@@ -198,8 +204,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_%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);
+		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);
 
 		FILE *cluster_out = fopen(c_filename, "rb");
 		FILE *avalanche_out = fopen(a_filename, "rb");
@@ -250,7 +256,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_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		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);
 
 		FILE *damage_out = fopen(d_filename, "rb");
 
@@ -292,7 +298,14 @@ 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, &c);
+		net_t *net = net_create(g, inf, beta, crack_len, use_voltage_boundaries, use_perc, &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;
+			if (save_damage) damage[net->num_broken]++;
+			if (save_energy) energy[i] = 0;
+			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);
 		net_free(tmp_net, &c);
@@ -349,7 +362,7 @@ int main(int argc, char *argv[]) {
 			}
 		}
 
-		if (save_damage) damage[max_pos]++;
+		if (save_damage) damage[net->num_broken]++;
 
 		if (save_stress_field || save_voltage_field) {
 			double *tmp_voltages = net_voltages(net, &c);
@@ -412,6 +425,17 @@ int main(int argc, char *argv[]) {
 			free(tmp_cluster_dist);
 		}
 
+
+		if (save_data) {
+			for (uint_t j = 0; j < data->num_broken; j++) {
+				fprintf(data_out, "%u %Lg %g ", data->break_list[j],
+								data->extern_field[j], data->conductivity[j]);
+			}
+			fprintf(data_out, "\n");
+		}
+
+		data_free(data);
+		}
 		if (save_network) {
 			FILE *net_out = fopen("network.txt", "w");
 			for (uint_t j = 0; j < g->nv; j++) {
@@ -440,16 +464,6 @@ int main(int argc, char *argv[]) {
 
 		net_free(net, &c);
 		graph_free(g, &c);
-
-		if (save_data) {
-			for (uint_t j = 0; j < data->num_broken; j++) {
-				fprintf(data_out, "%u %Lg %g ", data->break_list[j],
-								data->extern_field[j], data->conductivity[j]);
-			}
-			fprintf(data_out, "\n");
-		}
-
-		data_free(data);
 	}
 
 	printf("\033[F\033[JFRACTURE: COMPLETE\n");
@@ -472,7 +486,7 @@ int main(int argc, char *argv[]) {
 
 	if (save_voltage_field) {
 		char *vfld_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(vfld_filename, filename_len, "vfld_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		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);
@@ -482,7 +496,7 @@ int main(int argc, char *argv[]) {
 
 	if (save_stress_field) {
 		char *cfld_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(cfld_filename, filename_len, "cfld_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		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);
@@ -492,7 +506,7 @@ int main(int argc, char *argv[]) {
 
 	if (save_damage_field) {
 		char *dfld_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(dfld_filename, filename_len, "dfld_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		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);
@@ -502,7 +516,7 @@ int main(int argc, char *argv[]) {
 
 	if (save_conductivity) {
 		char *cond_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(cond_filename, filename_len, "cond_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		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);
 		FILE *cond_file = fopen(cond_filename, "ab");
 		fwrite(conductivity, sizeof(double), N, cond_file);
 		fclose(cond_file);
@@ -512,7 +526,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_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		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);
 		FILE *tough_file = fopen(tough_filename, "ab");
 		fwrite(energy, sizeof(long double), N, tough_file);
 		fclose(tough_file);
@@ -522,7 +536,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_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		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);
 		FILE *thres_file = fopen(thres_filename, "ab");
 		fwrite(thresholds, sizeof(long double), N, thres_file);
 		fclose(thres_file);
@@ -544,7 +558,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_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		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);
 		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 556e1e5..e02bcef 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -83,6 +83,7 @@ typedef struct {
 	cholmod_sparse *adjacency;
 	cholmod_sparse *dual_adjacency;
 	bool debug_stop;
+	uint_t num_broken;
 } net_t;
 
 typedef struct {
@@ -142,7 +143,7 @@ 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, 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);
 
 graph_t *ini_voro_graph(uint_t L, bound_t boundary, bool use_dual,
 													double *(*genfunc)(uint_t, bound_t, gsl_rng *, uint_t *),
diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c
index 73416e4..42134f0 100644
--- a/src/gen_laplacian.c
+++ b/src/gen_laplacian.c
@@ -208,8 +208,8 @@ cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c,
 
 
 	for (uint_t i = 0; i < num_gverts; i++) {
-		if (acoo[i] == 0)
-			acoo[i] = 1;
+		if (instance->marks[i] != instance->marks[network->nv])
+			acoo[i]++;
 	}
 
 	assert(CHOL_F(check_triplet)(temp_m, c));
diff --git a/src/net.c b/src/net.c
index 444e76c..1be5613 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) {
+long double *get_thres(uint_t ne, double beta, bool perc) {
 	long double *thres = (long double *)malloc(ne * sizeof(long double));
 	assert(thres != NULL);
 
@@ -9,7 +9,8 @@ long double *get_thres(uint_t ne, double beta) {
 	gsl_rng_set(r, rand_seed());
 
 	for (uint_t i = 0; i < ne; i++) {
-		thres[i] = rand_dist_pow(r, beta);
+		if (perc) thres[i] = 1;
+		else thres[i] = rand_dist_pow(r, beta);
 	}
 
 	gsl_rng_free(r);
@@ -44,14 +45,27 @@ void net_notch(net_t *net, double notch_len, cholmod_common *c) {
 	}
 }
 
-net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, 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 = (net_t *)calloc(1, sizeof(net_t));
 	assert(net != NULL);
 
 	net->graph = g;
+	net->num_broken = 0;
 	net->fuses = (bool *)calloc(g->ne, sizeof(bool));
 	assert(net->fuses != NULL);
-	net->thres = get_thres(g->ne, beta);
+	net->thres = get_thres(g->ne, beta, perc);
 	net->inf = inf;
 
 	net->voltage_bound = vb;
@@ -76,6 +90,8 @@ 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);