1 files changed, 86 insertions, 0 deletions
diff --git a/src/update_boundary.c b/src/update_boundary.c
new file mode 100644
index 0000000..3d2a086
--- /dev/null
+++ b/src/update_boundary.c
@@ -0,0 +1,86 @@
+
+#include "fracture.h"
+
+void update_boundary(finst *instance, const double *avg_field) {
+	int size = instance->network->num_edges;
+	int width = sqrt(size);
+	int num_verts = instance->network->num_verts;
+	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;
+		}
+	}
+}