#include "fracture.h"

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;
}

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)));

	return sqrt(r) * sin(th / 2);
}

void voronoi_bound_ini(net_t *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);
	}
}