Implemented the whole 2D Voronoi algorithm, except for the actual Voronoi cell construction.

examples/main.c

@@ -101,8 +101,8 @@ void print_help_message() {
       "parameter file.\n");
 }
 
-#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D)
-VORONOI3D_DECLARE_GLOBAL_VARIABLES()
+#if defined(SHADOWFAX_SPH)
+VORONOI_DECLARE_GLOBAL_VARIABLES()
 #endif
 
 /**
@@ -439,7 +439,7 @@ int main(int argc, char *argv[]) {
     fflush(stdout);
   }
 
-#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D)
+#if defined(SHADOWFAX_SPH)
   /* set the *global* box dimensions */
   float box_anchor[3], box_side[3];
   if (periodic) {

src/hydro/Shadowswift/voronoi1d_algorithm.h

@@ -27,6 +27,12 @@
 #include "inline.h"
 #include "voronoi1d_cell.h"
 
+#define VORONOI_DECLARE_GLOBAL_VARIABLES()
+
+__attribute__((always_inline)) INLINE static void voronoi_set_box(float *anchor,
+                                                                  float *side) {
+}
+
 /**
  * @brief Initialize a 1D Voronoi cell
  *

src/hydro/Shadowswift/voronoi2d_algorithm.h

@@ -25,8 +25,36 @@
 #include <stdlib.h>
 #include "error.h"
 #include "inline.h"
+#include "minmax.h"
 #include "voronoi2d_cell.h"
 
+#define VORONOI2D_BOX_LEFT 18446744073709551602llu
+#define VORONOI2D_BOX_RIGHT 18446744073709551603llu
+#define VORONOI2D_BOX_TOP 18446744073709551604llu
+#define VORONOI2D_BOX_BOTTOM 18446744073709551605llu
+
+extern float global_voronoi_box_anchor[2];
+extern float global_voronoi_box_side[2];
+
+#define VORONOI_DECLARE_GLOBAL_VARIABLES() \
+  float global_voronoi_box_anchor[2];      \
+  float global_voronoi_box_side[2];
+
+#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]
+
+__attribute__((always_inline)) INLINE static void voronoi_set_box(float *anchor,
+                                                                  float *side) {
+
+  global_voronoi_box_anchor[0] = anchor[0];
+  global_voronoi_box_anchor[1] = anchor[1];
+
+  global_voronoi_box_side[0] = side[0];
+  global_voronoi_box_side[1] = side[1];
+}
+
 /**
  * @brief Initialize a 2D Voronoi cell
  *
@@ -34,7 +62,38 @@
  * @param x Position of the generator of the cell.
  */
 __attribute__((always_inline)) INLINE void voronoi_cell_init(
-    struct voronoi_cell *cell, double *x) {}
+    struct voronoi_cell *cell, double *x) {
+
+  cell->x[0] = x[0];
+  cell->x[1] = x[1];
+
+  cell->nvert = 4;
+
+  cell->vertices[0][0] = VORONOI2D_BOX_ANCHOR_X - cell->x[0];
+  cell->vertices[0][1] = VORONOI2D_BOX_ANCHOR_Y - cell->x[1];
+
+  cell->vertices[1][0] = VORONOI2D_BOX_ANCHOR_X - cell->x[0];
+  cell->vertices[1][1] =
+      VORONOI2D_BOX_ANCHOR_Y + VORONOI2D_BOX_SIDE_Y - cell->x[1];
+
+  cell->vertices[2][0] =
+      VORONOI2D_BOX_ANCHOR_X + VORONOI2D_BOX_SIDE_X - cell->x[0];
+  cell->vertices[2][1] =
+      VORONOI2D_BOX_ANCHOR_Y + VORONOI2D_BOX_SIDE_Y - cell->x[1];
+
+  cell->vertices[3][0] =
+      VORONOI2D_BOX_ANCHOR_X + VORONOI2D_BOX_SIDE_X - cell->x[0];
+  cell->vertices[3][1] = VORONOI2D_BOX_ANCHOR_Y - cell->x[1];
+
+  cell->ngbs[0] = VORONOI2D_BOX_LEFT;
+  cell->ngbs[1] = VORONOI2D_BOX_TOP;
+  cell->ngbs[2] = VORONOI2D_BOX_RIGHT;
+  cell->ngbs[3] = VORONOI2D_BOX_BOTTOM;
+
+  cell->volume = 0.0f;
+  cell->centroid[0] = 0.0f;
+  cell->centroid[1] = 0.0f;
+}
 
 /**
  * @brief Interact a 2D Voronoi cell with a particle with given relative
@@ -57,7 +116,51 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
 __attribute__((always_inline)) INLINE float voronoi_cell_finalize(
     struct voronoi_cell *cell) {
 
-  return 1.0f;
+  int i;
+  float vertices[VORONOI2D_MAXNUMVERT][2];
+  float A, x[2], y[2], r2, r2max;
+
+  /* make a copy of the vertices (they are overwritten when the face midpoints
+     are computed */
+  for (i = 0; i < cell->nvert; ++i) {
+    vertices[i][0] = cell->vertices[i][0];
+    vertices[i][1] = cell->vertices[i][1];
+  }
+
+  r2max = 0.0f;
+  for (i = 0; i < cell->nvert; ++i) {
+    if (i < cell->nvert - 1) {
+      x[0] = vertices[i][0];
+      y[0] = vertices[i][1];
+      x[1] = vertices[i + 1][0];
+      y[1] = vertices[i + 1][1];
+    } else {
+      x[0] = vertices[i][0];
+      y[0] = vertices[i][1];
+      x[1] = vertices[0][0];
+      y[1] = vertices[0][1];
+    }
+    A = x[1] * y[0] - x[0] * y[1];
+    cell->volume += A;
+    cell->centroid[0] += (x[0] + x[1]) * A;
+    cell->centroid[1] += (y[0] + y[1]) * A;
+
+    cell->face_midpoints[i][0] = 0.5f * (x[0] + x[1]) + cell->x[0];
+    cell->face_midpoints[i][1] = 0.5f * (y[0] + y[1]) + cell->x[1];
+
+    r2 = x[0] * x[0] + y[0] * y[0];
+    r2max = max(r2max, r2);
+  }
+
+  cell->volume *= 0.5f;
+  A = 6 * cell->volume;
+  cell->centroid[0] /= A;
+  cell->centroid[1] /= A;
+
+  cell->centroid[0] += cell->x[0];
+  cell->centroid[1] += cell->x[1];
+
+  return 2.0f * sqrtf(r2max);
 }
 
 /**
@@ -72,7 +175,21 @@ __attribute__((always_inline)) INLINE float voronoi_cell_finalize(
 __attribute__((always_inline)) INLINE float voronoi_get_face(
     struct voronoi_cell *cell, unsigned long long ngb, float *midpoint) {
 
-  return 0.0f;
+  /* look up the neighbour */
+  int i = 0;
+  while (i < cell->nvert && cell->ngbs[i] != ngb) {
+    ++i;
+  }
+
+  if (i == cell->nvert) {
+    /* The given cell is not a neighbour. */
+    return 0.0f;
+  }
+
+  midpoint[0] = cell->face_midpoints[i][0];
+  midpoint[1] = cell->face_midpoints[i][1];
+
+  return cell->face_lengths[i];
 }
 
 /**
@@ -82,6 +199,10 @@ __attribute__((always_inline)) INLINE float voronoi_get_face(
  * @param centroid Array to store the centroid in.
  */
 __attribute__((always_inline)) INLINE void voronoi_get_centroid(
-    struct voronoi_cell *cell, float *centroid) {}
+    struct voronoi_cell *cell, float *centroid) {
+
+  centroid[0] = cell->centroid[0];
+  centroid[1] = cell->centroid[1];
+}
 
 #endif  // SWIFT_VORONOIXD_ALGORITHM_H

src/hydro/Shadowswift/voronoi2d_cell.h

@@ -20,6 +20,10 @@
 #ifndef SWIFT_VORONOIXD_CELL_H
 #define SWIFT_VORONOIXD_CELL_H
 
+/* Maximal number of vertices (and neighbours) that can be stored in a
+   voronoi_cell struct. */
+#define VORONOI2D_MAXNUMVERT 100
+
 /* 2D Voronoi cell */
 struct voronoi_cell {
 
@@ -31,6 +35,24 @@ struct voronoi_cell {
 
   /* The centroid of the cell. */
   float centroid[2];
+
+  /* Number of cell vertices (and neighbours). */
+  int nvert;
+
+  /* We only need to store one of these at the same time. */
+  union {
+    /* The relative positions of the vertices of the cell. */
+    float vertices[VORONOI2D_MAXNUMVERT][2];
+
+    /* The midpoints of the faces. */
+    float face_midpoints[VORONOI2D_MAXNUMVERT][2];
+  };
+
+  /* The ids of the neighbouring cells. */
+  unsigned long long ngbs[VORONOI2D_MAXNUMVERT];
+
+  /* The lengths of the faces. */
+  float face_lengths[VORONOI2D_MAXNUMVERT];
 };
 
 #endif  // SWIFT_VORONOIXD_CELL_H

src/hydro/Shadowswift/voronoi3d_algorithm.h

@@ -81,8 +81,8 @@ __attribute__((always_inline)) INLINE int check_counter(int *counter,
 extern float global_voronoi_box_anchor[3];
 extern float global_voronoi_box_side[3];
 
-#define VORONOI3D_DECLARE_GLOBAL_VARIABLES() \
-  float global_voronoi_box_anchor[3];        \
+#define VORONOI_DECLARE_GLOBAL_VARIABLES() \
+  float global_voronoi_box_anchor[3];      \
   float global_voronoi_box_side[3];
 
 /* Bottom front left corner and side lengths of the large box that contains all

tests/benchmarkInteractions.c

@@ -450,8 +450,8 @@ void test_interactions(struct part test_part, struct part *parts, size_t count,
 #endif
 }
 
-#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D)
-VORONOI3D_DECLARE_GLOBAL_VARIABLES()
+#if defined(SHADOWFAX_SPH)
+VORONOI_DECLARE_GLOBAL_VARIABLES()
 #endif
 
 /* And go... */

tests/test125cells.c

@@ -422,8 +422,8 @@ void runner_doself1_density(struct runner *r, struct cell *ci);
 void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
 void runner_doself2_force(struct runner *r, struct cell *ci);
 
-#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D)
-VORONOI3D_DECLARE_GLOBAL_VARIABLES()
+#if defined(SHADOWFAX_SPH)
+VORONOI_DECLARE_GLOBAL_VARIABLES()
 #endif
 
 /* And go... */

tests/test27cells.c

@@ -299,8 +299,8 @@ void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
 void runner_doself1_density(struct runner *r, struct cell *ci);
 void runner_doself1_density_vec(struct runner *r, struct cell *ci);
 
-#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D)
-VORONOI3D_DECLARE_GLOBAL_VARIABLES()
+#if defined(SHADOWFAX_SPH)
+VORONOI_DECLARE_GLOBAL_VARIABLES()
 #endif
 
 /* And go... */

tests/testPair.c

@@ -187,8 +187,8 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
 /* Just a forward declaration... */
 void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
 
-#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D)
-VORONOI3D_DECLARE_GLOBAL_VARIABLES()
+#if defined(SHADOWFAX_SPH)
+VORONOI_DECLARE_GLOBAL_VARIABLES()
 #endif
 
 int main(int argc, char *argv[]) {

tests/testSymmetry.c

@@ -26,8 +26,8 @@
 
 #include "swift.h"
 
-#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D)
-VORONOI3D_DECLARE_GLOBAL_VARIABLES()
+#if defined(SHADOWFAX_SPH)
+VORONOI_DECLARE_GLOBAL_VARIABLES()
 #endif
 
 int main(int argc, char *argv[]) {
@@ -35,7 +35,7 @@ int main(int argc, char *argv[]) {
   /* Choke if need be */
   feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
 
-#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D)
+#if defined(SHADOWFAX_SPH)
   /* Initialize the Voronoi simulation box */
   float box_anchor[3] = {-2.0f, -2.0f, -2.0f};
   float box_side[3] = {6.0f, 6.0f, 6.0f};

tests/testTimeIntegration.c

@@ -22,8 +22,8 @@
 #include <stdlib.h>
 #include <string.h>
 
-#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D)
-VORONOI3D_DECLARE_GLOBAL_VARIABLES()
+#if defined(SHADOWFAX_SPH)
+VORONOI_DECLARE_GLOBAL_VARIABLES()
 #endif
 
 /**

tests/testVoronoi2D.c

@@ -19,4 +19,27 @@
 
 #include "hydro/Shadowswift/voronoi2d_algorithm.h"
 
-int main() { return 0; }
+VORONOI_DECLARE_GLOBAL_VARIABLES()
+
+int main() {
+
+  float anchor[3] = {-0.5f, -0.5f, -0.5f};
+  float side[3] = {2.0f, 2.0f, 2.0f};
+  voronoi_set_box(anchor, side);
+
+  struct voronoi_cell cell;
+  double x[3] = {0.5, 0.5, 0.5};
+
+  voronoi_cell_init(&cell, x);
+
+  float maxradius = voronoi_cell_finalize(&cell);
+
+  assert(maxradius == 2.0f * sqrtf(2.0f));
+
+  assert(cell.volume == 4.0f);
+
+  assert(cell.centroid[0] == 0.5f);
+  assert(cell.centroid[1] == 0.5f);
+
+  return 0;
+}

tests/testVoronoi3D.c

@@ -1192,7 +1192,7 @@ void test_degeneracies() {
 #endif
 }
 
-VORONOI3D_DECLARE_GLOBAL_VARIABLES()
+VORONOI_DECLARE_GLOBAL_VARIABLES()
 
 int main() {