can generate voronoi networks with regular boundaries

6 files changed, 88 insertions, 66 deletions

diff --git a/makefile_hal b/makefile_hal
index 25ebf48..0b2f6ff 100644
--- a/makefile_hal
+++ b/makefile_hal
@@ -4,7 +4,7 @@ CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-fi
 LDFLAGS = -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lmetis -lgsl -lprofiler -ltcmalloc
 
 OBJ = break_data voronoi_bound_ini bin_values correlations beta_scales randfuncs instance get_dual_clusters coursegrain break_edge graph_components gen_laplacian geometry fracture_network get_current cracking_ini update_factor update_boundary get_file update_beta gen_voltcurmat ini_network free_network fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory get_conductivity
-BIN = fracture current_scaling course_grain_square corr_test homo_voronoi_fracture homo_square_fracture compare_voronoi_fracture
+BIN = current_scaling course_grain_square corr_test homo_voronoi_fracture homo_square_fracture compare_voronoi_fracture
 
 all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
 
diff --git a/src/fortune/main.c b/src/fortune/main.c
index 543f6e4..a1a4372 100644
--- a/src/fortune/main.c
+++ b/src/fortune/main.c
@@ -95,8 +95,8 @@ intptr_t *run_voronoi(unsigned int num, double *lattice, bool periodic, double x
 	double *eff_lattice = lattice;
 
 	if (periodic) {
-		unsigned int eff_num = 9 * num;
-		double *eff_lattice = (double *)malloc(2 * eff_num * sizeof(double));
+		eff_num = 9 * num;
+		eff_lattice = (double *)malloc(2 * eff_num * sizeof(double));
 
 		for (unsigned int i = 0; i < num; i++) {
 			// original sites - our baby boys
@@ -139,7 +139,7 @@ intptr_t *run_voronoi(unsigned int num, double *lattice, bool periodic, double x
 
 	nsites = eff_num;
 	readsites(eff_lattice);
-	free(eff_lattice);
+	if (periodic) free(eff_lattice);
 	next = nextone;
 
 	siteidx = 0;
diff --git a/src/fracture_network.c b/src/fracture_network.c
index 148b08c..fe520fd 100644
--- a/src/fracture_network.c
+++ b/src/fracture_network.c
@@ -39,9 +39,11 @@ break_data *fracture_network(finst *instance, double *fuse_thres,
 		}
 
 		int last_broke = inc_break_fuses(instance, fuse_thres, field, cutoff);
-		if (last_broke > num_edges || last_broke < -1) {
-			printf("%g \n", conductivity);
-			getchar();
+		if (last_broke > num_edges || last_broke < -1 || conductivity < 1e-8) {
+			printf("whoops %u\n\n", breaking_data->num_broken);
+			free(voltages);
+			free(field);
+			break;
 		}
 
 		update_break_data(breaking_data, last_broke, fabs(conductivity * fuse_thres[last_broke] / field[last_broke]), conductivity);
diff --git a/src/homo_voronoi_fracture.c b/src/homo_voronoi_fracture.c
index 65a25f8..26ceed4 100644
--- a/src/homo_voronoi_fracture.c
+++ b/src/homo_voronoi_fracture.c
@@ -7,7 +7,7 @@ int main(int argc, char *argv[]) {
 
 	// defining variables to be (potentially) set by command line flags
 	uint8_t filename_len;
-	uint64_t N;
+	uint32_t N;
 	uint_t L;
 	double beta, inf, cutoff;
 	bool include_breaking, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
@@ -63,6 +63,11 @@ int main(int argc, char *argv[]) {
 						use_voltage_boundaries = true;
 						boundc2 = 't';
 						break;
+					case 3:
+						boundary = EMBEDDED_BOUND;
+						boundc2 = 'e';
+						use_dual = true;
+						break;
 					default:
 						printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n");
 						exit(EXIT_FAILURE);
@@ -124,15 +129,15 @@ int main(int argc, char *argv[]) {
 	}
 
 	// define arrays for saving cluster and avalanche distributions
-	uint64_t *cluster_size_dist;
-	uint64_t *avalanche_size_dist;
+	uint32_t *cluster_size_dist;
+	uint32_t *avalanche_size_dist;
 	char *c_filename;
 	char *a_filename;
 	if (save_cluster_dist) {
 		cluster_size_dist =
-				(uint64_t *)malloc(c_dist_size * sizeof(uint64_t));
+				(uint32_t *)malloc(c_dist_size * sizeof(uint32_t));
 		avalanche_size_dist =
-				(uint64_t *)malloc(a_dist_size * sizeof(uint64_t));
+				(uint32_t *)malloc(a_dist_size * sizeof(uint32_t));
 
 		c_filename = (char *)malloc(filename_len * sizeof(char));
 		a_filename = (char *)malloc(filename_len * sizeof(char));
@@ -143,11 +148,11 @@ int main(int argc, char *argv[]) {
 		FILE *avalanche_out = fopen(a_filename, "rb");
 
 		if (cluster_out != NULL) {
-			fread(cluster_size_dist, sizeof(uint64_t), c_dist_size, cluster_out);
+			fread(cluster_size_dist, sizeof(uint32_t), c_dist_size, cluster_out);
 			fclose(cluster_out);
 		}
 		if (avalanche_out != NULL) {
-			fread(avalanche_size_dist, sizeof(uint64_t), a_dist_size, avalanche_out);
+			fread(avalanche_size_dist, sizeof(uint32_t), a_dist_size, avalanche_out);
 			fclose(avalanche_out);
 		}
 	}
@@ -163,11 +168,11 @@ int main(int argc, char *argv[]) {
 	}
 
 	// define arrays for saving damage distributions
-	uint64_t *damage;
+	uint32_t *damage;
 	char *d_filename;
 	if (save_damage) {
 		damage =
-				(uint64_t *)malloc(a_dist_size * sizeof(uint64_t));
+				(uint32_t *)malloc(a_dist_size * sizeof(uint32_t));
 
 		d_filename = (char *)malloc(filename_len * sizeof(char));
 		snprintf(d_filename, filename_len, "damg_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta);
@@ -175,7 +180,7 @@ int main(int argc, char *argv[]) {
 		FILE *damage_out = fopen(d_filename, "rb");
 
 		if (damage_out != NULL) {
-			fread(damage, sizeof(uint64_t), a_dist_size, damage_out);
+			fread(damage, sizeof(uint32_t), a_dist_size, damage_out);
 			fclose(damage_out);
 		}
 	}
@@ -203,7 +208,7 @@ int main(int argc, char *argv[]) {
 
 
 	printf("\n");
-	for (uint64_t i = 0; i < N; i++) {
+	for (uint32_t i = 0; i < N; i++) {
 		printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1, N);
 
 		fnet *network = ini_voronoi_network(L, boundary, use_dual, genfunc_hyperuniform, &c);
@@ -340,8 +345,8 @@ int main(int argc, char *argv[]) {
 		FILE *cluster_out = fopen(c_filename, "wb");
 		FILE *avalanche_out = fopen(a_filename, "wb");
 
-		fwrite(cluster_size_dist, sizeof(uint64_t), c_dist_size, cluster_out);
-		fwrite(avalanche_size_dist, sizeof(uint64_t), a_dist_size, avalanche_out);
+		fwrite(cluster_size_dist, sizeof(uint32_t), c_dist_size, cluster_out);
+		fwrite(avalanche_size_dist, sizeof(uint32_t), a_dist_size, avalanche_out);
 
 		fclose(cluster_out);
 		fclose(avalanche_out);
@@ -374,7 +379,7 @@ int main(int argc, char *argv[]) {
 
 	if (save_damage) {
 		FILE *hdam_file = fopen(d_filename, "wb");
-		fwrite(damage, sizeof(uint64_t), a_dist_size, hdam_file);
+		fwrite(damage, sizeof(uint32_t), a_dist_size, hdam_file);
 		fclose(hdam_file);
 		free(d_filename);
 		free(damage);
diff --git a/src/ini_network.c b/src/ini_network.c
index edf1539..263b1b2 100644
--- a/src/ini_network.c
+++ b/src/ini_network.c
@@ -23,12 +23,14 @@ double *get_edge_coords(unsigned int num_edges, double *vert_coords,
 												unsigned int *edges_to_verts) {
 	double *output = (double *)malloc(2 * num_edges * sizeof(double));
 
-#pragma omp parallel for
+	#pragma omp parallel for
 	for (unsigned int i = 0; i < num_edges; i++) {
 		unsigned int v1, v2;
 		double v1x, v1y, v2x, v2y, dx, dy;
 		v1 = edges_to_verts[2 * i];
 		v2 = edges_to_verts[2 * i + 1];
+		output[2 * i] = 0;
+		output[2 * i + 1] = 0;
 		v1x = vert_coords[2 * v1];
 		v1y = vert_coords[2 * v1 + 1];
 		v2x = vert_coords[2 * v2];
@@ -270,7 +272,8 @@ unsigned int *get_voro_dual_edges(unsigned int num_edges,
 																	unsigned int *triangles) {
 	unsigned int *dual_edges =
 			(unsigned int *)malloc(2 * num_edges * sizeof(unsigned int));
-#pragma omp parallel for
+	unsigned int place = 0;
+	#pragma omp parallel for
 	for (unsigned int i = 0; i < num_edges; i++) {
 		unsigned int v1, v2;
 		v1 = edges[2 * i];
@@ -285,8 +288,9 @@ unsigned int *get_voro_dual_edges(unsigned int num_edges,
 					t21 = triangles[3 * v2 + k];
 					t22 = triangles[3 * v2 + ((k + 1) % 3)];
 					if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) {
-						dual_edges[2 * i] = t11 < t12 ? t11 : t12;
-						dual_edges[2 * i + 1] = t11 < t12 ? t12 : t11;
+						dual_edges[2 * place] = t11 < t12 ? t11 : t12;
+						dual_edges[2 * place + 1] = t11 < t12 ? t12 : t11;
+						place++;
 						found_match = true;
 						break;
 					}
@@ -361,7 +365,7 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 	// prune the edges of the lattice and assign boundary vertices based on the
 	// desired boundary conditions
 	unsigned int num_bounds;
-	unsigned num_verts;
+	unsigned int num_verts;
 	double *vert_coords;
 	unsigned int *bound_inds;
 	unsigned int *bound_verts;
@@ -379,9 +383,9 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
 			dual_edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
 			unsigned int num_t, num_b, num_l, num_r;
-			bool *bound_t, *bound_b, *bound_l, *bound_r;
+			bool *bound_top, *bound_b, *bound_l, *bound_r;
 			num_t = 0; num_b = 0; num_l = 0; num_r = 0;
-			bound_t = (bool *)calloc(num_verts, sizeof(bool));
+			bound_top = (bool *)calloc(num_verts, sizeof(bool));
 			bound_b = (bool *)calloc(num_verts, sizeof(bool));
 			bound_l = (bool *)calloc(num_verts, sizeof(bool));
 			bound_r = (bool *)calloc(num_verts, sizeof(bool));
@@ -394,15 +398,15 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 				dx = v1x - v2x; dy = v1y - v2y;
 				if (fabs(dy) > 0.5) {
 					if (dy > 0) {
-						if (!bound_t[v1] && !bound_l[v1] && !bound_r[v1]) {
-							bound_t[v1] = true; num_t++;
+						if (!bound_top[v1] && !bound_l[v1] && !bound_r[v1]) {
+							bound_top[v1] = true; num_t++;
 						}
 						if (!bound_b[v2] && !bound_l[v2] && !bound_r[v2]) {
 							bound_b[v2] = true; num_b++;
 						}
 					} else {
-						if (!bound_t[v2] && !bound_l[v2] && !bound_r[v2]) {
-							bound_t[v2] = true; num_t++;
+						if (!bound_top[v2] && !bound_l[v2] && !bound_r[v2]) {
+							bound_top[v2] = true; num_t++;
 						}
 						if (!bound_b[v1] && !bound_l[v1] && !bound_r[v1]) {
 							bound_b[v1] = true; num_b++;
@@ -410,17 +414,17 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 					}
 				} else if (fabs(dx) > 0.5) {
 					if (dx > 0) {
-						if (!bound_r[v1] && !bound_t[v1] && !bound_b[v1]) {
+						if (!bound_r[v1] && !bound_top[v1] && !bound_b[v1]) {
 							bound_r[v1] = true; num_r++;
 						}
-						if (!bound_l[v2] && !bound_t[v2] && !bound_b[v2]) {
+						if (!bound_l[v2] && !bound_top[v2] && !bound_b[v2]) {
 							bound_l[v2] = true; num_l++;
 						}
 					} else {
-						if (!bound_r[v2] && !bound_t[v2] && !bound_b[v2]) {
+						if (!bound_r[v2] && !bound_top[v2] && !bound_b[v2]) {
 							bound_r[v2] = true; num_r++;
 						}
-						if (!bound_l[v1] && !bound_t[v1] && !bound_b[v1]) {
+						if (!bound_l[v1] && !bound_top[v1] && !bound_b[v1]) {
 							bound_l[v1] = true; num_l++;
 						}
 					}
@@ -441,7 +445,7 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			unsigned int pos_t, pos_b, pos_l, pos_r;
 			pos_t = 0; pos_b = 0; pos_l = 0; pos_r = 0;
 			for (unsigned int i = 0; i < num_verts; i++) {
-				if (bound_t[i]) {
+				if (bound_top[i]) {
 					bound_verts[pos_t] = i; pos_t++;
 				} else if (bound_b[i]) {
 					bound_verts[num_t + pos_b] = i; pos_b++;
@@ -451,9 +455,10 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 					bound_verts[num_t + num_b + num_l + pos_r] = i; pos_r++;
 				}
 			}
-			free(bound_l); free(bound_r); free(bound_t); free(bound_b);
+			free(bound_l); free(bound_r); free(bound_top); free(bound_b);
 			free(tmp_edges);
 			free(tmp_dual_edges);
+			num_bounds = 2;
 			network->edges_to_verts_break = edges;
 			network->edges_to_verts = edges;
 			network->num_verts_break = num_verts;
@@ -471,9 +476,9 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
 			dual_edges = (unsigned int *)malloc(2 * tmp_num_edges * sizeof(unsigned int));
 			unsigned int num_t, num_b;
-			bool *bound_t, *bound_b;
+			bool *bound_top, *bound_b;
 			num_t = 0; num_b = 0;
-			bound_t = (bool *)calloc(num_verts, sizeof(bool));
+			bound_top = (bool *)calloc(num_verts, sizeof(bool));
 			bound_b = (bool *)calloc(num_verts, sizeof(bool));
 			for (unsigned int i = 0; i < tmp_num_edges; i++) {
 				unsigned int v1, v2;
@@ -483,15 +488,15 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 				dy = v1y - v2y;
 				if (fabs(dy) > 0.5) {
 					if (dy > 0) {
-						if (!bound_t[v1]) {
-							bound_t[v1] = true; num_t++;
+						if (!bound_top[v1]) {
+							bound_top[v1] = true; num_t++;
 						}
 						if (!bound_b[v2]) {
 							bound_b[v2] = true; num_b++;
 						}
 					} else {
-						if (!bound_t[v2]) {
-							bound_t[v2] = true; num_t++;
+						if (!bound_top[v2]) {
+							bound_top[v2] = true; num_t++;
 						}
 						if (!bound_b[v1]) {
 							bound_b[v1] = true; num_b++;
@@ -512,13 +517,13 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			unsigned int pos_t, pos_b;
 			pos_t = 0; pos_b = 0;
 			for (unsigned int i = 0; i < num_verts; i++) {
-				if (bound_t[i]) {
+				if (bound_top[i]) {
 					bound_verts[pos_t] = i; pos_t++;
 				} else if (bound_b[i]) {
 					bound_verts[num_t + pos_b] = i; pos_b++;
 				}
 			}
-			free(bound_t); free(bound_b);
+			free(bound_top); free(bound_b);
 			free(tmp_edges);
 			free(tmp_dual_edges);
 			network->edges_to_verts_break = edges;
@@ -539,7 +544,7 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 				edges[2*i+1] = tmp_edges[2*i+1];
 			}
 			dual_edges = tmp_dual_edges;
-			bool *bound_t = (bool *)calloc(tmp_num_verts, sizeof(bool));
+			bool *bound_top = (bool *)calloc(tmp_num_verts, sizeof(bool));
 			int *edge_change = (int *)calloc(num_edges, sizeof(int));
 			unsigned int num_t = 0;
 			for (unsigned int i = 0; i < num_edges; i++) {
@@ -552,13 +557,13 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 				if (fabs(dy) > 0.5) {
 					if (dy > 0) {
 						edge_change[i] = 1;
-						if (!bound_t[v1]) {
-							bound_t[v1] = true; num_t++;
+						if (!bound_top[v1]) {
+							bound_top[v1] = true; num_t++;
 						}
 					} else {
 						edge_change[i] = 2;
-						if (!bound_t[v2]) {
-							bound_t[v2] = true; num_t++;
+						if (!bound_top[v2]) {
+							bound_top[v2] = true; num_t++;
 						}
 					}
 				}
@@ -571,7 +576,7 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			for (unsigned int i = 0; i < tmp_num_verts; i++) {
 				vert_coords[2*i] = tmp_vert_coords[2*i];
 				vert_coords[2*i+1] = tmp_vert_coords[2*i+1];
-				if (bound_t[i]) {
+				if (bound_top[i]) {
 					bound_verts[pos_t] = i;
 					vert_coords[2*(tmp_num_verts + pos_t)] = tmp_vert_coords[2*i];
 					vert_coords[2*(tmp_num_verts + pos_t)+1] = tmp_vert_coords[2*i+1];
@@ -589,7 +594,7 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 				}
 			}
 			free(tmp_vert_coords);
-			free(bound_t);
+			free(bound_top);
 			free(edge_change);
 			network->num_verts_break = num_verts;
 			network->num_verts = tmp_num_verts;
@@ -599,17 +604,25 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			break;
 											}
 		case EMBEDDED_BOUND: {
-			num_bounds = 4;
+			num_bounds = 2;
 			bound_inds = (unsigned int *)malloc(5 * sizeof(unsigned int));
 			bound_verts = (unsigned int *)malloc(2 * L * sizeof(unsigned int));
 			for (unsigned int i = 0; i < 5; i++) bound_inds[i] = i * L / 2;
 			for (unsigned int i = 0; i < 2 * L; i++) bound_verts[i] = i;
-			num_edges = tmp_num_edges;
+			unsigned int num_away = 0;
+			for (unsigned int i = 0; i < tmp_num_edges; i++) {
+				if (tmp_dual_edges[2*i] > num || tmp_dual_edges[2*i+1] > num) num_away++;
+			}
+			num_edges = (int)tmp_num_edges - (int)num_away;
 			num_verts = tmp_num_verts;
 			edges = tmp_edges;
 			dual_edges = tmp_dual_edges;
-			num_edges = tmp_num_edges;
 			vert_coords = tmp_vert_coords;
+			network->num_verts_break = num_verts;
+			network->num_verts = num_verts;
+			network->edges_to_verts_break = edges;
+			network->edges_to_verts = edges;
+			network->break_dim = num_verts + num_bounds;
 		}
 	}
 
diff --git a/src/randfuncs.c b/src/randfuncs.c
index 35b9f56..830e078 100644
--- a/src/randfuncs.c
+++ b/src/randfuncs.c
@@ -22,23 +22,25 @@ double *genfunc_hyperuniform(unsigned int L, bound_t boundary, gsl_rng *r, unsig
 	double *lattice = (double *)malloc(2 * (*num) * sizeof(double));
 	double rho = *num;
 	unsigned int to_gen = *num;
+	unsigned int start = 0;
 
 	if (boundary == EMBEDDED_BOUND) {
 		for (unsigned int i = 0; i < L / 2; i++) {
-			lattice[2 * i] = 0;
-			lattice[2 * i + 1] = (2. * i + 1.) / L;
-
-			lattice[L / 2 + 2 * i] = 1;
-			lattice[L / 2 + 2 * i + 1] = (2. * i + 1.) / L;
+			lattice[2 * i + 1] = 0;
+			lattice[2 * i] = (2. * i + 1.) / L;
 
+			lattice[L + 2 * i + 1] = 1;
 			lattice[L + 2 * i] = (2. * i + 1.) / L;
-			lattice[L + 2 * i + 1] = 0;
 
-			lattice[3 * L / 2 + 2 * i] = (2. * i + 1.) / L;
-			lattice[3 * L / 2 + 2 * i + 1] = 1;
+			lattice[2 * L + 2 * i + 1] = (2. * i + 1.) / L;
+			lattice[2 * L + 2 * i] = 0;
+
+			lattice[3 * L + 2 * i + 1] = (2. * i + 1.) / L;
+			lattice[3 * L + 2 * i] = 1;
 		}
 
 		to_gen -= 2 * L;
+		start = 2 * L;
 	}
 
 	for (unsigned int i = 0; i < to_gen; i++) {
@@ -67,8 +69,8 @@ double *genfunc_hyperuniform(unsigned int L, bound_t boundary, gsl_rng *r, unsig
 				}
 			}
 		}
-		lattice[2 * i] = x;
-		lattice[2 * i + 1] = y;
+		lattice[2*start + 2 * i] = x;
+		lattice[2*start + 2 * i + 1] = y;
 	}
 
 	return lattice;