#include "fracture.h"

void update_boundary(net_t *instance, const double *avg_field) {
	int size = instance->graph->ne;
	int width = sqrt(size);
	int num_verts = instance->graph->nv;
	double *boundary = (double *)instance->boundary_cond->x;

	boundary[num_verts] = 0;
	boundary[num_verts + 1] = 0;

	if (instance->voltage_bound) {
		for (int i = 0; i < width / 4; i++) {
			int t_ind = (width + 1) * (width / 8) + width / 4 - 1 + i;
			double t_val = avg_field[t_ind];
			boundary[2 * i] = t_val;
			boundary[2 * i + 1] = t_val;

			int b_ind = num_verts - (width + 1) * (width / 8) - width / 4 - i;
			double b_val = avg_field[b_ind];
			boundary[num_verts - 1 - 2 * i] = b_val;
			boundary[num_verts - 2 - 2 * i] = b_val;

			int l_ind = (width + 1) * (width / 8) + width / 2 + width / 4 - 1 +
									i * (width + 1);
			double l_val = avg_field[l_ind];
			boundary[width / 2 + 2 * i * (width + 1)] = l_val;
			boundary[width / 2 + (2 * i + 1) * (width + 1)] = l_val;

			int r_ind = (width + 1) * (width / 8) + width - 1 + i * (width + 1);
			double r_val = avg_field[r_ind];
			boundary[width + 2 * i * (width + 1)] = r_val;
			boundary[width + (2 * i + 1) * (width + 1)] = r_val;
		}
	} else {
		const double *stress = avg_field;

		for (int i = 0; i < width / 4; i++) {
			boundary[2 * i] = stress[(width / 4 - 1) * width + width / 4 + 2 * i] +
												stress[(width / 4 - 1) * width + width / 4 + 1 + 2 * i];
			boundary[2 * i + 1] =
					stress[(width / 4 - 1) * width + width / 4 + 2 * i] +
					stress[(width / 4 - 1) * width + width / 4 + 1 + 2 * i];

			boundary[num_verts - 2 * i - 1] =
					-stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 2 * i)] -
					stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 1 + 2 * i)];
			boundary[num_verts - 2 * i - 2] =
					-stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 2 * i)] -
					stress[size - 1 - ((width / 4 - 1) * width + width / 4 + 1 + 2 * i)];

			boundary[(width + 1) / 2 + 2 * i * (width + 1)] =
					stress[((width / 4) + 2 * i) * width + width / 4 - 1] -
					stress[((width / 4) + 2 * i + 1) * width + width / 4 - 1];
			boundary[(width + 1) / 2 + (2 * i + 1) * (width + 1)] =
					stress[((width / 4) + 2 * i) * width + width / 4 - 1] -
					stress[((width / 4) + 2 * i + 1) * width + width / 4 - 1];

			boundary[(width + 1) / 2 + 2 * i * (width + 1) + width / 2] =
					stress[((width / 4) + 2 * i) * width + width / 4 + width / 2] -
					stress[((width / 4) + 2 * i + 1) * width + width / 4 + width / 2];
			boundary[(width + 1) / 2 + (2 * i + 1) * (width + 1) + width / 2] =
					stress[((width / 4) + 2 * i) * width + width / 4 + width / 2] -
					stress[((width / 4) + 2 * i + 1) * width + width / 4 + width / 2];
		}

		double s_total = 0;
		double total = 0;
		for (int i = 0; i < num_verts; i++) {
			if (boundary[i] != boundary[i]) {
				boundary[i] = 0;
			} else {
				total += fabs(boundary[i]);
				s_total += boundary[i];
			}
		}

		// assert(fabs(s_total) < inf);
		boundary[num_verts] -= s_total;

		for (int i = 0; i < num_verts + 1; i++) {
			boundary[i] /= total;
		}
	}
}