From fd14c5e39d962be94a1f68b0d4cacb7a4aa9c3e7 Mon Sep 17 00:00:00 2001
From: pants <jaron@kent-dobias.com>
Date: Fri, 2 Sep 2016 15:24:34 -0400
Subject: embedded systems with crack fully supported

---
 src/voronoi_bound_ini.c | 98 +++++++++++++------------------------------------
 1 file changed, 26 insertions(+), 72 deletions(-)

(limited to 'src/voronoi_bound_ini.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);
+	}
 }
+
-- 
cgit v1.2.3-70-g09d2