Finished first part of 2D Voronoi intersection algorithm. Now *only* need to...

Finished first part of 2D Voronoi intersection algorithm. Now *only* need to write the second part. And maybe deal with degeneracies...

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 */
 }
 
 /**