diff --git a/src/hydro/Shadowswift/voronoi2d_algorithm.h b/src/hydro/Shadowswift/voronoi2d_algorithm.h
index a27b652601589fd4cfe4faeffad58adf12114d24..c956c731ac84b89d6b5c8a239554f305b2bdae6d 100644
--- a/src/hydro/Shadowswift/voronoi2d_algorithm.h
+++ b/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) {
diff --git a/tests/testVoronoi2D.c b/tests/testVoronoi2D.c
index 285eb88a917f2ee8605c697225ea92e53f1f3348..47eb40bf0e0a922f75296c6d08e0d49448490bcb 100644
--- a/tests/testVoronoi2D.c
+++ b/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;
}