embedded systems with crack fully supported

9 files changed, 504 insertions, 88 deletions

diff --git a/src/break_edge.c b/src/break_edge.c
index 0ba1ad7..72ccf3e 100644
--- a/src/break_edge.c
+++ b/src/break_edge.c
@@ -8,7 +8,7 @@ bool break_edge(finst *instance, unsigned int edge, cholmod_common *c) {
 	unsigned int v1 = instance->network->edges_to_verts_break[2 * edge];
 	unsigned int v2 = instance->network->edges_to_verts_break[2 * edge + 1];
 
-	update_factor(instance->factor, v1, v2, c);
+	if (instance->factor != NULL) update_factor(instance->factor, v1, v2, c);
 
 	if (instance->network->boundary != TORUS_BOUND) {
 		unsigned int w1 = instance->network->edges_to_verts[2 * edge];
diff --git a/src/cracked_voronoi_fracture.c b/src/cracked_voronoi_fracture.c
new file mode 100644
index 0000000..532c3c3
--- /dev/null
+++ b/src/cracked_voronoi_fracture.c
@@ -0,0 +1,427 @@
+
+
+#include "fracture.h"
+
+int main(int argc, char *argv[]) {
+	int opt;
+
+	// defining variables to be (potentially) set by command line flags
+	uint8_t filename_len;
+	uint32_t N;
+	uint_t L;
+	double beta, inf, cutoff, crack_len;
+	bool include_breaking, 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;
+
+	filename_len = 100;
+
+	N = 100;
+	L = 16;
+	crack_len = 0.;
+	beta = .3;
+	inf = 1e10;
+	cutoff = 1e-9;
+	boundary = FREE_BOUND;
+	include_breaking = false;
+	save_cluster_dist = false;
+	use_voltage_boundaries = false;
+	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_toughness = false;
+
+	uint8_t bound_i;
+	char boundc2 = 'f';
+
+	while ((opt = getopt(argc, argv, "n:L:b:B:dVcoNsCrtDSvel:")) != -1) {
+		switch (opt) {
+			case 'n':
+				N = atoi(optarg);
+				break;
+			case 'L':
+				L = atoi(optarg);
+				break;
+			case 'b':
+				beta = atof(optarg);
+				break;
+			case 'l':
+				crack_len = atof(optarg);
+				break;
+			case 'B':
+				bound_i = atoi(optarg);
+				switch (bound_i) {
+					case 0:
+						boundary = FREE_BOUND;
+						boundc2 = 'f';
+						break;
+					case 1:
+						boundary = CYLINDER_BOUND;
+						boundc2 = 'c';
+						break;
+					case 2:
+						boundary = TORUS_BOUND;
+						use_voltage_boundaries = true;
+						boundc2 = 't';
+						break;
+					case 3:
+						boundary = EMBEDDED_BOUND;
+						boundc2 = 'e';
+						use_dual = true;
+						use_voltage_boundaries = true;
+						break;
+					default:
+						printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n");
+						exit(EXIT_FAILURE);
+				}
+				break;
+			case 'd':
+				save_damage = true;
+				break;
+			case 'e':
+				save_damage_field = true;
+				break;
+			case 'V':
+				use_voltage_boundaries = true;
+				break;
+			case 'D':
+				use_dual = true;
+				break;
+			case 'c':
+				save_cluster_dist = true;
+				break;
+			case 'o':
+				include_breaking = true;
+				break;
+			case 'N':
+				save_network = true;
+				break;
+			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;
+				//inf = 1;
+				break;
+			case 't':
+				save_toughness = true;
+				break;
+			default: /* '?' */
+				exit(EXIT_FAILURE);
+		}
+	}
+
+	char boundc;
+	if (use_voltage_boundaries) boundc = 'v';
+	else boundc = 'c';
+
+	FILE *break_out;
+	if (include_breaking) {
+		char *break_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(break_filename, filename_len, "breaks_v_%c_%c_%u_%g.txt", boundc, boundc2, L, beta);
+		break_out = fopen(break_filename, "a");
+		free(break_filename);
+	}
+
+	uint_t voronoi_max_verts, c_dist_size, a_dist_size;
+
+	voronoi_max_verts = 4 * pow(L, 2);
+	c_dist_size = voronoi_max_verts;
+	a_dist_size = voronoi_max_verts;
+
+	if (voronoi_max_verts > CINT_MAX) {
+		exit(EXIT_FAILURE);
+	}
+
+	// define arrays for saving cluster and avalanche distributions
+	uint32_t *cluster_size_dist;
+	uint32_t *avalanche_size_dist;
+	char *c_filename;
+	char *a_filename;
+	if (save_cluster_dist) {
+		cluster_size_dist =
+				(uint32_t *)malloc(c_dist_size * sizeof(uint32_t));
+		avalanche_size_dist =
+				(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));
+		snprintf(c_filename, filename_len, "cstr_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta);
+		snprintf(a_filename, filename_len, "avln_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta);
+
+		FILE *cluster_out = fopen(c_filename, "rb");
+		FILE *avalanche_out = fopen(a_filename, "rb");
+
+		if (cluster_out != NULL) {
+			fread(cluster_size_dist, sizeof(uint32_t), c_dist_size, cluster_out);
+			fclose(cluster_out);
+		}
+		if (avalanche_out != NULL) {
+			fread(avalanche_size_dist, sizeof(uint32_t), a_dist_size, avalanche_out);
+			fclose(avalanche_out);
+		}
+	}
+
+	double *crit_stress;
+	if (save_crit_stress) {
+		crit_stress = (double *)malloc(N * sizeof(double));
+	}
+
+	double *conductivity;
+	if (save_conductivity) {
+		conductivity = (double *)malloc(N * sizeof(double));
+	}
+
+	// define arrays for saving damage distributions
+	uint32_t *damage;
+	char *d_filename;
+	if (save_damage) {
+		damage =
+				(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);
+
+		FILE *damage_out = fopen(d_filename, "rb");
+
+		if (damage_out != NULL) {
+			fread(damage, sizeof(uint32_t), a_dist_size, damage_out);
+			fclose(damage_out);
+		}
+	}
+
+	double *toughness;
+	if (save_toughness) {
+		toughness = (double *)malloc(N * sizeof(double));
+	}
+
+
+	// start cholmod
+	cholmod_common c;
+	CHOL_F(start)(&c);
+
+	/* if we use voltage boundary conditions, the laplacian matrix is positive
+	 * definite and we can use a supernodal LL decomposition.  otherwise we need
+	 * to use the simplicial LDL decomposition
+	 */
+	if (use_voltage_boundaries) {
+		//(&c)->supernodal = CHOLMOD_SUPERNODAL;
+		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
+	} else {
+		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
+	}
+
+
+	printf("\n");
+	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);
+		finst *perm_instance = create_instance(network, inf, use_voltage_boundaries, false, &c);
+		if (boundary == EMBEDDED_BOUND) {
+			voronoi_bound_ini(perm_instance, L, crack_len);
+		}
+		gen_voro_crack(perm_instance, crack_len, &c);
+		finish_instance(perm_instance, &c);
+		double *fuse_thres = gen_fuse_thres(network->num_edges, network->edge_coords, beta, beta_scaling_flat);
+		finst *instance = copy_instance(perm_instance, &c);
+		break_data *breaking_data = fracture_network(instance, fuse_thres, &c, cutoff);
+		free_instance(instance, &c);
+		free(fuse_thres);
+
+		uint_t max_pos = 0;
+		double max_val = 0;
+
+		for (uint_t j = 0; j < breaking_data->num_broken; j++) {
+			double val = breaking_data->extern_field[j];
+
+			if (val > max_val) {
+				max_pos = j;
+				max_val = val;
+			}
+		}
+
+		if (save_crit_stress) {
+			crit_stress[i] =
+					breaking_data->extern_field[max_pos];
+		}
+
+		finst *tmp_instance = copy_instance(perm_instance, &c);
+
+		uint_t av_size = 0;
+		double cur_val = 0;
+		for (uint_t j = 0; j < max_pos; j++) {
+			break_edge(tmp_instance, breaking_data->break_list[j], &c);
+
+			double val = breaking_data->extern_field[j];
+			if (save_cluster_dist) {
+				if (val < cur_val) {
+					av_size++;
+				}
+
+				if (val > cur_val) {
+					avalanche_size_dist[av_size]++;
+					av_size = 0;
+					cur_val = val;
+				}
+			}
+		}
+
+		if (save_conductivity) {
+			if (!use_voltage_boundaries) {
+				double *tmp_voltage = get_voltage(tmp_instance, &c);
+				conductivity[i] = 1/fabs(tmp_voltage[tmp_instance->network->num_verts + 1] - tmp_voltage[tmp_instance->network->num_verts]);
+				free(tmp_voltage);
+			} else {
+				conductivity[i] = breaking_data->conductivity[max_pos];
+			}
+		}
+
+		if (save_toughness) {
+			double tmp_toughness = 0;
+			if (max_pos > 0) {
+				double sigma1 = breaking_data->extern_field[0];
+				double epsilon1 = sigma1 / breaking_data->conductivity[0];
+				for (uint_t j = 0; j < max_pos - 1; j++) {
+					double sigma2 = breaking_data->extern_field[j+1];
+					double epsilon2 = sigma2 / breaking_data->conductivity[j+1];
+					if (epsilon2 > epsilon1) {
+						tmp_toughness += (sigma1 + sigma2) * (epsilon2 - epsilon1) / 2;
+						sigma1 = sigma2; epsilon1 = epsilon2;
+					}
+				}
+			}
+			toughness[i] = tmp_toughness;
+		}
+
+		if (save_damage) {
+			damage[max_pos]++;
+		}
+
+		if (save_cluster_dist) {
+			uint_t *tmp_cluster_dist = get_cluster_dist(tmp_instance, &c);
+			for (uint_t j = 0; j < tmp_instance->network->num_dual_verts; j++) {
+				cluster_size_dist[j] += tmp_cluster_dist[j];
+			}
+			free(tmp_cluster_dist);
+		}
+
+		if (save_network) {
+			FILE *net_out = fopen("network.txt", "w");
+			for (uint_t j = 0; j < network->num_verts; j++) {
+				fprintf(net_out, "%f %f ", network->vert_coords[2 * j],
+								tmp_instance->network->vert_coords[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (uint_t j = 0; j < tmp_instance->network->num_edges; j++) {
+				fprintf(net_out, "%u %u ", tmp_instance->network->edges_to_verts[2 * j],
+								tmp_instance->network->edges_to_verts[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (uint_t j = 0; j < tmp_instance->network->num_dual_verts; j++) {
+				fprintf(net_out, "%f %f ", tmp_instance->network->dual_vert_coords[2 * j],
+								tmp_instance->network->dual_vert_coords[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (uint_t j = 0; j < tmp_instance->network->num_edges; j++) {
+				fprintf(net_out, "%u %u ", tmp_instance->network->dual_edges_to_verts[2 * j],
+								tmp_instance->network->dual_edges_to_verts[2 * j + 1]);
+			}
+			fprintf(net_out, "\n");
+			for (uint_t j = 0; j < tmp_instance->network->num_edges; j++) {
+				fprintf(net_out, "%d ", tmp_instance->fuses[j]);
+			}
+			fclose(net_out);
+		}
+
+		free_instance(tmp_instance, &c);
+		free_instance(perm_instance, &c);
+		free_fnet(network, &c);
+
+		if (include_breaking) {
+			for (uint_t j = 0; j < breaking_data->num_broken; j++) {
+				fprintf(break_out, "%u %g %g ", breaking_data->break_list[j],
+								breaking_data->extern_field[j], breaking_data->conductivity[j]);
+			}
+			fprintf(break_out, "\n");
+		}
+
+		free_break_data(breaking_data);
+	}
+
+	printf("\033[F\033[JFRACTURE: COMPLETE\n");
+
+	if (save_cluster_dist) {
+		FILE *cluster_out = fopen(c_filename, "wb");
+		FILE *avalanche_out = fopen(a_filename, "wb");
+
+		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);
+
+		free(c_filename);
+		free(a_filename);
+		free(cluster_size_dist);
+		free(avalanche_size_dist);
+	}
+
+	if (save_conductivity) {
+		char *cond_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(cond_filename, filename_len, "cond_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta);
+		FILE *cond_file = fopen(cond_filename, "ab");
+		fwrite(conductivity, sizeof(double), N, cond_file);
+		fclose(cond_file);
+		free(cond_filename);
+		free(conductivity);
+	}
+
+	if (save_toughness) {
+		char *tough_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(tough_filename, filename_len, "tuff_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta);
+		FILE *tough_file = fopen(tough_filename, "ab");
+		fwrite(toughness, sizeof(double), N, tough_file);
+		fclose(tough_file);
+		free(tough_filename);
+		free(toughness);
+	}
+
+	if (save_damage) {
+		FILE *hdam_file = fopen(d_filename, "wb");
+		fwrite(damage, sizeof(uint32_t), a_dist_size, hdam_file);
+		fclose(hdam_file);
+		free(d_filename);
+		free(damage);
+	}
+
+	if (include_breaking) {
+		fclose(break_out);
+	}
+
+	if (save_crit_stress) {
+		char *str_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(str_filename, filename_len, "strs_v_%c_%c_%d_%g.dat", boundc, boundc2, L, beta);
+		FILE *str_file = fopen(str_filename, "ab");
+		fwrite(crit_stress, sizeof(double), N, str_file);
+		fclose(str_file);
+		free(str_filename);
+		free(crit_stress);
+	}
+
+	CHOL_F(finish)(&c);
+
+	return 0;
+}
diff --git a/src/cracking_ini.c b/src/cracking_ini.c
index 9ab4a9d..988af5b 100644
--- a/src/cracking_ini.c
+++ b/src/cracking_ini.c
@@ -31,6 +31,30 @@ double *gen_fuse_thres(unsigned int num_edges, double *edge_coords, double beta,
 	return fuse_thres;
 }
 
+void gen_voro_crack(finst *instance, double crack_len, cholmod_common *c) {
+	for (uint_t i = 0; i < instance->network->num_edges; i++) {
+		uint_t v1, v2;
+		double v1x, v1y, v2x, v2y, dx, dy;
+
+		v1 = instance->network->edges_to_verts[2 * i];
+		v2 = instance->network->edges_to_verts[2 * i + 1];
+
+		v1x = instance->network->vert_coords[2 * v1];
+		v1y = instance->network->vert_coords[2 * v1 + 1];
+		v2x = instance->network->vert_coords[2 * v2];
+		v2y = instance->network->vert_coords[2 * v2 + 1];
+
+		dx = v1x - v2x;
+		dy = v1y - v2y;
+
+		if (((v1y > 0.5 && v2y < 0.5) || (v1y < 0.5 && v2y > 0.5)) && fabs(dy) < 0.5) {
+			if (v1x + dx / dy * (v1y - 0.5) <= crack_len) {
+				break_edge(instance, i, c);
+			}
+		}
+	}
+}
+
 bool gen_crack(finst *instance, double crack_len, double crack_width,
 							 cholmod_common *c) {
 	assert(instance != NULL);
diff --git a/src/fracture.h b/src/fracture.h
index 91b2bfc..93e2ed1 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -175,11 +175,12 @@ double *get_corr(finst *instance, unsigned int **dists, cholmod_common *c);
 
 double *bin_values(fnet *network, unsigned int width, double *values);
 
-void voronoi_bound_ini(finst *instance, double *square_bound,
-											 unsigned int width);
+void voronoi_bound_ini(finst *instance, uint_t L, double crack_len);
 
 break_data *alloc_break_data(unsigned int num_edges);
 void free_break_data(break_data *data);
 void update_break_data(break_data *data, unsigned int last_broke, double strength, double conductivity);
 
 double get_conductivity(finst *inst, double *current, cholmod_common *c);
+
+void gen_voro_crack(finst *instance, double crack_len, cholmod_common *c);
diff --git a/src/fracture_network.c b/src/fracture_network.c
index fe520fd..3f06104 100644
--- a/src/fracture_network.c
+++ b/src/fracture_network.c
@@ -53,7 +53,11 @@ break_data *fracture_network(finst *instance, double *fuse_thres,
 
 		break_edge(instance, last_broke, c);
 
-		if (instance->num_components > 1) {
+		if (instance->num_components > 1 && instance->network->boundary == TORUS_BOUND) {
+			break;
+		}
+
+		if (instance->marks[num_verts] != instance->marks[num_verts + 1] && instance->network->boundary != TORUS_BOUND) {
 			break;
 		}
 	}
diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c
index d8f0c6d..97a2c9d 100644
--- a/src/gen_laplacian.c
+++ b/src/gen_laplacian.c
@@ -159,14 +159,15 @@ cholmod_sparse *gen_laplacian(const finst *instance, cholmod_common *c,
 	}
 
 	if (network->boundary != TORUS_BOUND) {
+		if (network->boundary != EMBEDDED_BOUND) acoo[0]++;
+
 		if (voltage_bound) {
-			for (unsigned int i = 0; i < 2; i++) {
+			for (unsigned int i = 0; i < num_bounds; i++) {
 				for (unsigned int j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) {
 					acoo[bound_verts[bound_inds[i] + j]]++;
 				}
 			}
 		} else {
-			acoo[0]++;
 			for (unsigned int i = 0; i < num_bounds; i++) {
 				rowind[num_verts + i] = num_verts + i;
 				colind[num_verts + i] = num_verts + i;
@@ -181,7 +182,7 @@ cholmod_sparse *gen_laplacian(const finst *instance, cholmod_common *c,
 
 			unsigned int start = num_gverts;
 
-			for (unsigned int i = 0; i < num_bounds; i++) {
+			for (unsigned int i = 0; i < 2; i++) {
 				for (unsigned int j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) {
 					rowind[start + bound_inds[i] + j] = bound_verts[bound_inds[i] + j];
 					colind[start + bound_inds[i] + j] = num_verts + i;
diff --git a/src/ini_network.c b/src/ini_network.c
index 263b1b2..d80f43b 100644
--- a/src/ini_network.c
+++ b/src/ini_network.c
@@ -604,7 +604,7 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			break;
 											}
 		case EMBEDDED_BOUND: {
-			num_bounds = 2;
+			num_bounds = 4;
 			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;
@@ -622,7 +622,7 @@ fnet *ini_voronoi_network(unsigned int L, bound_t boundary, bool use_dual,
 			network->num_verts = num_verts;
 			network->edges_to_verts_break = edges;
 			network->edges_to_verts = edges;
-			network->break_dim = num_verts + num_bounds;
+			network->break_dim = num_verts + 2;
 		}
 	}
 
diff --git a/src/instance.c b/src/instance.c
index 098ed92..cc33591 100644
--- a/src/instance.c
+++ b/src/instance.c
@@ -14,7 +14,18 @@ finst *create_instance(fnet *network, double inf, bool voltage_bound,
 	instance->voltage_bound = voltage_bound;
 	instance->boundary_cond = CHOL_F(zeros)(
 			network->break_dim, 1, CHOLMOD_REAL, c);
-	if (network->boundary != TORUS_BOUND) {
+	if (network->boundary == TORUS_BOUND) {
+		for (unsigned int i = 0; i < network->bound_inds[1]; i++) {
+			((double *)instance->boundary_cond->x)[network->bound_verts[i]] = 1;
+			((double *)instance->boundary_cond->x)[network->num_verts + i] = -1;
+		}
+		((double *)instance->boundary_cond->x)[network->bound_verts[0]] = 1;
+	} else if (network->boundary == EMBEDDED_BOUND) {
+		// do nothing
+		for (unsigned int i = 0; i < network->bound_inds[1]; i++) {
+			((double *)instance->boundary_cond->x)[network->bound_verts[i]] =1;
+		}
+	} else {
 		if (voltage_bound) {
 			for (unsigned int i = 0; i < network->bound_inds[1]; i++) {
 				((double *)instance->boundary_cond->x)[network->bound_verts[i]] =1;
@@ -24,12 +35,6 @@ finst *create_instance(fnet *network, double inf, bool voltage_bound,
 			((double *)instance->boundary_cond->x)[network->num_verts] = 1;
 			((double *)instance->boundary_cond->x)[network->num_verts + 1] = -1;
 		}
-	} else {
-		for (unsigned int i = 0; i < network->bound_inds[1]; i++) {
-			((double *)instance->boundary_cond->x)[network->bound_verts[i]] = 1;
-			((double *)instance->boundary_cond->x)[network->num_verts + i] = -1;
-		}
-		((double *)instance->boundary_cond->x)[network->bound_verts[0]] = 1;
 	}
 
 	if (network->boundary != TORUS_BOUND) instance->adjacency = gen_adjacency(instance, false, false, 0, c);
diff --git a/src/voronoi_bound_ini.c b/src/voronoi_bound_ini.c
index 38f41cc..d201093 100644
--- a/src/voronoi_bound_ini.c
+++ b/src/voronoi_bound_ini.c
@@ -1,79 +1,33 @@
 
 #include "fracture.h"
 
-void voronoi_bound_ini(finst *instance, double *square_bound,
-											 unsigned int width) {
-	unsigned int square_num_verts = 2 * ((width + 1) / 2) * (width / 2 + 1);
-	double total1 = 0;
-	double total2 = 0;
-	double total3 = 0;
-	double total4 = 0;
-	for (unsigned int i = 0; i < instance->network->num_bounds; i++) {
-		for (unsigned int j = 0; j < instance->network->bound_inds[i + 1] -
-																		 instance->network->bound_inds[i];
-				 j++) {
-			double x = instance->network
-										 ->vert_coords[2 *
-																	 instance->network->bound_verts
-																			 [instance->network->bound_inds[i] + j]];
-			double y =
-					instance->network
-							->vert_coords[2 *
-																instance->network->bound_verts
-																		[instance->network->bound_inds[i] + j] +
-														1];
+double th_p(double x, double y, double th) {
+	if (x >= 0 && y >= 0) return th;
+	if (x <  0 && y >= 0) return M_PI - th;
+	if (x <  0 && y <  0) return th - M_PI;
+	if (x >= 0 && y <  0) return -th;
+}
 
-			unsigned int vw = (width + 1) / 2;
-			unsigned int xp = (unsigned int)(x * vw);
-			unsigned int yp = (unsigned int)(y * vw);
+double u_y(double x, double y) {
+	double r = sqrt(pow(x, 2) + pow(y, 2));
+	double th = th_p(x, y, atan(fabs(y / x)));
 
-			if (i == 0) {
-				((double *)instance->boundary_cond
-						 ->x)[instance->network
-											->bound_verts[instance->network->bound_inds[i] + j]] =
-						square_bound[xp];
-				total1 += square_bound[xp];
-			}
-			if (i == 1) {
-				((double *)instance->boundary_cond
-						 ->x)[instance->network
-											->bound_verts[instance->network->bound_inds[i] + j]] =
-						square_bound[square_num_verts - vw - 1 + xp];
-				total2 += square_bound[square_num_verts - vw - 1 + xp];
-			}
-			if (i == 2) {
-				((double *)instance->boundary_cond
-						 ->x)[instance->network
-											->bound_verts[instance->network->bound_inds[i] + j]] =
-						square_bound[(width + 1) / 2 + yp * (width + 1)];
-				total3 += square_bound[(width + 1) / 2 + yp * (width + 1)];
-			}
-			if (i == 3) {
-				((double *)instance->boundary_cond
-						 ->x)[instance->network
-											->bound_verts[instance->network->bound_inds[i] + j]] =
-						square_bound[(width + 1) / 2 + yp * (width + 1) + width / 2];
-				total4 += square_bound[(width + 1) / 2 + yp * (width + 1) + width / 2];
-			}
-		}
-	}
+	return sqrt(r) * sin(th / 2);
+}
 
-	((double *)instance->boundary_cond->x)[instance->network->num_verts] =
-			- total1 / 2 - total2 / 2;
-	((double *)instance->boundary_cond->x)[instance->network->num_verts + 1] =
-			- total1 / 2 - total2 / 2;
-	((double *)instance->boundary_cond->x)[instance->network->num_verts + 2] =
-			-total3;
-	((double *)instance->boundary_cond->x)[instance->network->num_verts + 3] =
-			-total4;
-	/*
-	((double *)instance->boundary_cond->x)[instance->network->num_verts] =
-			0;
-	((double *)instance->boundary_cond->x)[instance->network->num_verts + 1] =
-			0;
-	((double *)instance->boundary_cond->x)[instance->network->num_verts + 2] =
-			0;
-	((double *)instance->boundary_cond->x)[instance->network->num_verts + 3] =
-			0;
-			*/
+void voronoi_bound_ini(finst *instance, uint_t L, double crack_len) {
+	double *bound = (double *)instance->boundary_cond->x;
+	for (uint_t i = 0; i < L / 2; i++) {
+		double x1, y1, x2, y2, x3, y3, x4, y4;
+		x1 = (2. * i + 1.) / L - crack_len; y1 = 0.5 - 1.;
+		x2 = (2. * i + 1.) / L - crack_len; y2 = 0.5 - 0.;
+		y3 = (2. * i + 1.) / L - 0.5; x3 = 0.5 - 1.;
+		y4 = (2. * i + 1.) / L - 0.5; x4 = 0.5 - 0.;
+
+		bound[i] = u_y(x1, y1);
+		bound[L / 2 + i] = u_y(x2, y2);
+		bound[L + i] = u_y(x3, y3);
+		bound[3 * L / 2 + i] = u_y(x4, y4);
+	}
 }
+