diff --git a/src/hydro/Shadowswift/voronoi2d_algorithm.h b/src/hydro/Shadowswift/voronoi2d_algorithm.h
index d02f690cdc7357ce795c3f836c9519de759ee0e3..717c929f6ab02e4072f209bd0b1a8794c394e249 100644
--- a/src/hydro/Shadowswift/voronoi2d_algorithm.h
+++ b/src/hydro/Shadowswift/voronoi2d_algorithm.h
@@ -28,23 +28,39 @@
#include "minmax.h"
#include "voronoi2d_cell.h"
+/* IDs used to keep track of cells neighbouring walls of the simulation box
+ This will only work if these IDs are never used for actual particles (which
+ in practice means you want to have less than 2^63-4 (~9e18) particles in your
+ simulation) */
#define VORONOI2D_BOX_LEFT 18446744073709551602llu
#define VORONOI2D_BOX_RIGHT 18446744073709551603llu
#define VORONOI2D_BOX_TOP 18446744073709551604llu
#define VORONOI2D_BOX_BOTTOM 18446744073709551605llu
+/* Global variables used to store the extent of the simulation box */
extern float global_voronoi_box_anchor[2];
extern float global_voronoi_box_side[2];
+/* Macro used to declare the global variables in any piece of code that uses
+ this header file */
#define VORONOI_DECLARE_GLOBAL_VARIABLES() \
float global_voronoi_box_anchor[2]; \
float global_voronoi_box_side[2];
+/* Macros that make it easier to change the box extents that are used throughout
+ the code */
#define VORONOI2D_BOX_ANCHOR_X global_voronoi_box_anchor[0]
#define VORONOI2D_BOX_ANCHOR_Y global_voronoi_box_anchor[1]
#define VORONOI2D_BOX_SIDE_X global_voronoi_box_side[0]
#define VORONOI2D_BOX_SIDE_Y global_voronoi_box_side[1]
+/**
+ * @brief Set the values of the global variables that store the simulation box
+ * extents.
+ *
+ * @param anchor Corner of the box with the lowest coordinate values.
+ * @param side Side lengths of the box.
+ */
__attribute__((always_inline)) INLINE static void voronoi_set_box(float *anchor,
float *side) {
@@ -64,9 +80,12 @@ __attribute__((always_inline)) INLINE static void voronoi_set_box(float *anchor,
__attribute__((always_inline)) INLINE void voronoi_cell_init(
struct voronoi_cell *cell, double *x) {
+ /* Set the position of the generator of the cell (for reference) */
cell->x[0] = x[0];
cell->x[1] = x[1];
+ /* Initialize the cell as a box with the same extents as the simulation box
+ (note: all vertex coordinates are relative w.r.t. the cell generator) */
cell->nvert = 4;
cell->vertices[0][0] = VORONOI2D_BOX_ANCHOR_X - cell->x[0];
@@ -85,11 +104,16 @@ __attribute__((always_inline)) INLINE void voronoi_cell_init(
VORONOI2D_BOX_ANCHOR_X + VORONOI2D_BOX_SIDE_X - cell->x[0];
cell->vertices[3][1] = VORONOI2D_BOX_ANCHOR_Y - cell->x[1];
+ /* The neighbours are ordered such that neighbour i shares the face in between
+ vertices i and i+1 (with last vertex + 1 = first vertex)
+ We walk around the cell in clockwise direction */
cell->ngbs[0] = VORONOI2D_BOX_LEFT;
cell->ngbs[1] = VORONOI2D_BOX_TOP;
cell->ngbs[2] = VORONOI2D_BOX_RIGHT;
cell->ngbs[3] = VORONOI2D_BOX_BOTTOM;
+ /* These variables are initialized to zero, we will compute them after the
+ neighbour iteration has finished */
cell->volume = 0.0f;
cell->centroid[0] = 0.0f;
cell->centroid[1] = 0.0f;
@@ -109,6 +133,24 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
float half_dx[2];
float r2;
+ float test;
+ int i;
+ int index_above1, index_above2;
+ int index_below1, index_below2;
+ int increment;
+
+ /* The process of cutting the current cell with the midline of the generator
+ and the given relative neighbour position proceeds in two steps:
+ - first we need to locate an edge of the current cell that is intersected
+ by this midline. Such an edge does not necessarily exist; in this case
+ the given neighbour is not an actual neighbour of this cell
+ - Once we have an intersected edge, we create a new edge starting at the
+ intersection point. We follow the edges connected to the intersected
+ edge until we find another intersected edge, and use its intersection
+ point as end point of the new edge. */
+
+ /* First, we set up some variables that are used to check if a vertex is above
+ or below the midplane. */
/* we need a vector with half the size of the vector joining generator and
neighbour, pointing to the neighbour */
@@ -118,7 +160,130 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
/* we need the squared length of this vector */
r2 = half_dx[0] * half_dx[0] + half_dx[1] * half_dx[1];
- (void)r2;
+ /* a vertex v = (vx, vy) is above the midline if
+ vx*half_dx[0] + vy*half_dx[1] > r2
+ i.e., if the length of the projected vertex position is longer than the
+ length of the vector pointing to the closest point on the midline (both
+ vectors originate at the position of the generator)
+ the vertex is below the midline if the projected position vector is shorter
+ if the projected position vector has the same length, the vertex is on the
+ midline */
+
+ /* start testing a random vertex: the first one */
+ test = cell->vertices[0][0] * half_dx[0] + cell->vertices[0][1] * half_dx[1] -
+ r2;
+ if (test < 0.) {
+ /* vertex is below midline */
+
+ /* move on until we find a vertex that is above the midline */
+ i = 1;
+ test = cell->vertices[i][0] * half_dx[0] +
+ cell->vertices[i][1] * half_dx[1] - r2;
+ while (i < cell->nvert && test < 0.) {
+ ++i;
+ test = cell->vertices[i][0] * half_dx[0] +
+ cell->vertices[i][1] * half_dx[1] - r2;
+ }
+
+ /* loop finished, there are two possibilities:
+ - i == cell->nvert, all vertices lie below the midline and the given
+ neighbour is not an actual neighbour of this cell
+ - test >= 0., we found a vertex above (or on) the midline */
+ if (i == cell->nvert) {
+ /* the given neighbour is not an actual neighbour: exit the routine */
+ return;
+ }
+
+ /* we have found an intersected edge: i-1 -> i
+ we store the index of the vertex above the midline, and set the increment
+ for the consecutive vertex search to +1 */
+ index_below1 = i - 1;
+ index_above1 = i;
+ increment = 1;
+ } else {
+ /* vertex is above midline */
+
+ /* move on until we find a vertex that is below the midline */
+ i = 1;
+ test = cell->vertices[i][0] * half_dx[0] +
+ cell->vertices[i][1] * half_dx[1] - r2;
+ while (i < cell->nvert && test > 0.) {
+ ++i;
+ test = cell->vertices[i][0] * half_dx[0] +
+ cell->vertices[i][1] * half_dx[1] - r2;
+ }
+
+ /* loop finished, there are two possibilities:
+ - i == cell->nvert, all vertices lie above the midline. This should
+ never happen.
+ - test <= 0., we found a vertex below (or on) the midline */
+ if (i == cell->nvert) {
+ /* fatal error! */
+ error("Could not find a vertex below the midline!");
+ }
+
+ /* we have found an intersected edge: i-1 -> i
+ we store the index of the vertex above the midline, and set the increment
+ for the consecutive vertex search to -1 */
+ index_below1 = i;
+ index_above1 = i - 1;
+ increment = -1;
+ }
+
+ /* now we need to find the second intersected edge
+ we start from the vertex above the midline and search in the direction
+ opposite to the intersected edge direction until we find a vertex below the
+ midline */
+ i = index_above1 + increment;
+ if (i < 0) {
+ i = cell->nvert - 1;
+ }
+ if (i == cell->nvert) {
+ i = 0;
+ }
+ test = cell->vertices[i][0] * half_dx[0] + cell->vertices[i][1] * half_dx[1] -
+ r2;
+ /* this loop can never deadlock, as we know there is at least 1 vertex below
+ the midline */
+ while (test > 0.) {
+ i += increment;
+ if (i < 0) {
+ i = cell->nvert - 1;
+ }
+ if (i == cell->nvert) {
+ i = 0;
+ }
+ test = cell->vertices[i][0] * half_dx[0] +
+ cell->vertices[i][1] * half_dx[1] - r2;
+ }
+
+ if (i == index_below1) {
+ /* we only found 1 intersected edge. This is impossible. */
+ error("Only 1 intersected edge found!");
+ }
+
+ index_below2 = i;
+ index_above2 = i - increment;
+ if (index_above2 < 0) {
+ index_above2 = cell->nvert - 1;
+ }
+ if (index_above2 == cell->nvert) {
+ index_above2 = 0;
+ }
+
+ if (increment < 0) {
+ /* interchange index_below and above 1 and 2 */
+ i = index_below1;
+ index_below1 = index_below2;
+ index_below2 = i;
+ i = index_above1;
+ index_above1 = index_above2;
+ index_above2 = i;
+ }
+
+ /* index_above1 now holds the first vertex to remove, and index_above2 the
+ last, in clockwise order
+ they can be the same */
}
/**