Fixed a bug in 2D Voronoi algorithm, added a lot of debugging output...

Fixed a bug in 2D Voronoi algorithm, added a lot of debugging output possibilities. 2D algorithm now works for 100 uniform random generators in a periodic box.

src/hydro/Shadowswift/voronoi2d_algorithm.h

@@ -28,6 +28,13 @@
 #include "minmax.h"
 #include "voronoi2d_cell.h"
 
+/* Check if the test variable contains a value that is small enough to be
+   affected by round off error */
+#define VORONOI_CHECK_TEST()     \
+  if (fabs(test) < 1.e-10) {     \
+    error("Too close to call!"); \
+  }
+
 /* 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
@@ -180,8 +187,14 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
   /* start testing a random vertex: the first one */
   test = cell->vertices[0][0] * half_dx[0] + cell->vertices[0][1] * half_dx[1] -
          r2;
+  VORONOI_CHECK_TEST();
   if (test < 0.) {
-    /* vertex is below midline */
+/* vertex is below midline */
+#ifdef VORONOI_VERBOSE
+    message("First vertex is below midline (%g %g --> %g)!",
+            cell->vertices[0][0] + cell->x[0],
+            cell->vertices[0][1] + cell->x[1], test);
+#endif
 
     /* store the test result; we might need it to compute the intersection
        coordinates */
@@ -191,12 +204,14 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
     i = 1;
     test = cell->vertices[i][0] * half_dx[0] +
            cell->vertices[i][1] * half_dx[1] - r2;
+    VORONOI_CHECK_TEST();
     while (i < cell->nvert && test < 0.) {
       /* make sure we always store the latest test result */
       b1 = test;
       ++i;
       test = cell->vertices[i][0] * half_dx[0] +
              cell->vertices[i][1] * half_dx[1] - r2;
+      VORONOI_CHECK_TEST();
     }
 
     /* loop finished, there are two possibilities:
@@ -204,7 +219,10 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
           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 */
+/* the given neighbour is not an actual neighbour: exit the routine */
+#ifdef VORONOI_VERBOSE
+      message("Not a neighbour!");
+#endif
       return;
     }
 
@@ -218,7 +236,12 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
     a1 = test;
     increment = 1;
   } else {
-    /* vertex is above midline */
+/* vertex is above midline */
+#ifdef VORONOI_VERBOSE
+    message("First vertex is above midline (%g %g --> %g)!",
+            cell->vertices[0][0] + cell->x[0],
+            cell->vertices[0][1] + cell->x[1], test);
+#endif
 
     /* store the test result */
     a1 = test;
@@ -227,12 +250,14 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
     i = 1;
     test = cell->vertices[i][0] * half_dx[0] +
            cell->vertices[i][1] * half_dx[1] - r2;
+    VORONOI_CHECK_TEST();
     while (i < cell->nvert && test > 0.) {
       /* make sure we always store the most recent test result */
       a1 = test;
       ++i;
       test = cell->vertices[i][0] * half_dx[0] +
              cell->vertices[i][1] * half_dx[1] - r2;
+      VORONOI_CHECK_TEST();
     }
 
     /* loop finished, there are two possibilities:
@@ -255,6 +280,15 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
     b1 = test;
   }
 
+#ifdef VORONOI_VERBOSE
+  message("First intersected edge: %g %g --> %g %g (%i --> %i)",
+          cell->vertices[index_below1][0] + cell->x[0],
+          cell->vertices[index_below1][1] + cell->x[1],
+          cell->vertices[index_above1][0] + cell->x[0],
+          cell->vertices[index_above1][1] + cell->x[1], index_below1,
+          index_above1);
+#endif
+
   /* 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
@@ -273,6 +307,7 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
   }
   test = cell->vertices[i][0] * half_dx[0] + cell->vertices[i][1] * half_dx[1] -
          r2;
+  VORONOI_CHECK_TEST();
   /* this loop can never deadlock, as we know there is at least 1 vertex below
      the midline */
   while (test > 0.) {
@@ -287,11 +322,7 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
     }
     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!");
+    VORONOI_CHECK_TEST();
   }
 
   index_below2 = i;
@@ -305,6 +336,13 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
   /* we also store the test result for the second vertex below the midline */
   b2 = test;
 
+  if (index_above1 == index_above2 && index_below1 == index_below2) {
+    /* There can be only 1 vertex above or below the midline, but we need 2
+       intersected edges, so if the vertices above the midline are the same, the
+       ones below need to be different and vice versa */
+    error("Only 1 intersected edge found!");
+  }
+
   /* to make the code below more clear, we make sure index_above1 always holds
      the first vertex to remove, and index_above2 the last one, in clockwise
      order
@@ -325,6 +363,21 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
     a2 = test;
   }
 
+#ifdef VORONOI_VERBOSE
+  message("First vertex below: %g %g (%i, %g)",
+          cell->vertices[index_below1][0] + cell->x[0],
+          cell->vertices[index_below1][1] + cell->x[1], index_below1, b1);
+  message("First vertex above: %g %g (%i, %g)",
+          cell->vertices[index_above1][0] + cell->x[0],
+          cell->vertices[index_above1][1] + cell->x[1], index_above1, a1);
+  message("Second vertex below: %g %g (%i, %g)",
+          cell->vertices[index_below2][0] + cell->x[0],
+          cell->vertices[index_below2][1] + cell->x[1], index_below2, b2);
+  message("Second vertex above: %g %g (%i, %g)",
+          cell->vertices[index_above2][0] + cell->x[0],
+          cell->vertices[index_above2][1] + cell->x[1], index_above2, a2);
+#endif
+
   /* convert the test results (which correspond to the projected distance
      between the vertex and the midline) to the fractions of the intersected
      edges above and below the midline */
@@ -484,6 +537,11 @@ __attribute__((always_inline)) INLINE void voronoi_print_cell(
     struct voronoi_cell *cell) {
 
   int i, ip1;
+
+  /* print cell generator */
+  printf("%g %g\n\n", cell->x[0], cell->x[1]);
+
+  /* print cell vertices */
   for (i = 0; i < cell->nvert; ++i) {
     ip1 = i + 1;
     if (ip1 == cell->nvert) {

tests/testVoronoi2D.c

@@ -19,35 +19,96 @@
 
 #include "hydro/Shadowswift/voronoi2d_algorithm.h"
 
+/* Number of cells used to test the 2D interaction algorithm */
+#define TESTVORONOI2D_NUMCELL 100
+
+/* declare the global variables needed to store the simulation box size */
 VORONOI_DECLARE_GLOBAL_VARIABLES()
 
 int main() {
 
+  /* initialize simulation box */
   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};
+  /* test initialization and finalization algorithms */
+  {
+    struct voronoi_cell cell;
+    double x[3] = {0.5, 0.5, 0.5};
+
+    voronoi_cell_init(&cell, x);
+
+    float maxradius = voronoi_cell_finalize(&cell);
 
-  voronoi_cell_init(&cell, x);
+    assert(maxradius == 2.0f * sqrtf(2.0f));
 
-  float maxradius = voronoi_cell_finalize(&cell);
+    assert(cell.volume == 4.0f);
 
-  assert(maxradius == 2.0f * sqrtf(2.0f));
+    assert(cell.centroid[0] == 0.5f);
+    assert(cell.centroid[1] == 0.5f);
+  }
 
-  assert(cell.volume == 4.0f);
+  /* test interaction algorithm */
+  {
+    /* create 100 cells with random positions in [0,1]x[0,1] */
+    struct voronoi_cell cells[TESTVORONOI2D_NUMCELL];
+    double x[2];
+    float dx[2];
+    int i, j;
+    float Atot;
+    struct voronoi_cell *cell_i, *cell_j;
 
-  assert(cell.centroid[0] == 0.5f);
-  assert(cell.centroid[1] == 0.5f);
+    for (i = 0; i < TESTVORONOI2D_NUMCELL; ++i) {
+      x[0] = ((double)rand()) / ((double)RAND_MAX);
+      x[1] = ((double)rand()) / ((double)RAND_MAX);
+      voronoi_cell_init(&cells[i], x);
+#ifdef VORONOI_VERBOSE
+      message("cell[%i]: %g %g", i, x[0], x[1]);
+#endif
+    }
 
-  /* reinitialize cell */
-  voronoi_cell_init(&cell, x);
+    /* interact the cells (with periodic boundaries) */
+    for (i = 0; i < TESTVORONOI2D_NUMCELL; ++i) {
+      cell_i = &cells[i];
+      for (j = 0; j < TESTVORONOI2D_NUMCELL; ++j) {
+        if (i != j) {
+          cell_j = &cells[j];
+          dx[0] = cell_i->x[0] - cell_j->x[0];
+          dx[1] = cell_i->x[1] - cell_j->x[1];
+          /* periodic boundaries */
+          if (dx[0] >= 0.5f) {
+            dx[0] -= 1.0f;
+          }
+          if (dx[0] < -0.5f) {
+            dx[0] += 1.0f;
+          }
+          if (dx[1] >= 0.5f) {
+            dx[1] -= 1.0f;
+          }
+          if (dx[1] < -0.5f) {
+            dx[1] += 1.0f;
+          }
+#ifdef VORONOI_VERBOSE
+          message("Cell %i before:", i);
+          voronoi_print_cell(&cells[i]);
+          message("Interacting cell %i with cell %i (%g %g, %g %g", i, j,
+                  cells[i].x[0], cells[i].x[1], cells[j].x[0], cells[j].x[1]);
+#endif
+          voronoi_cell_interact(cell_i, dx, j);
+        }
+      }
+    }
 
-  float dx[2] = {0.25, 0.25};
-  voronoi_cell_interact(&cell, dx, 0);
+    /* print the cells to the stdout */
+    for (i = 0; i < TESTVORONOI2D_NUMCELL; ++i) {
+      voronoi_print_cell(&cells[i]);
+      voronoi_cell_finalize(&cells[i]);
+      Atot += cells[i].volume;
+    }
 
-  voronoi_print_cell(&cell);
+    assert(fabs(Atot - 1.0f) < 1.e-6);
+  }
 
   return 0;
 }