diff --git a/src/hydro/Shadowswift/voronoi2d_algorithm.h b/src/hydro/Shadowswift/voronoi2d_algorithm.h
index c956c731ac84b89d6b5c8a239554f305b2bdae6d..1a30d02e558f5b3810ff0ae25f66a264597a8fd9 100644
--- a/src/hydro/Shadowswift/voronoi2d_algorithm.h
+++ b/src/hydro/Shadowswift/voronoi2d_algorithm.h
@@ -28,11 +28,10 @@
#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!"); \
+/* Check if the number of vertices exceeds the maximal allowed number */
+#define VORONOI_CHECK_SIZE() \
+ if (nvert > VORONOI2D_MAXNUMVERT) { \
+ error("Too many vertices!"); \
}
/* IDs used to keep track of cells neighbouring walls of the simulation box
@@ -187,7 +186,6 @@ __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 */
#ifdef VORONOI_VERBOSE
@@ -200,18 +198,16 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
coordinates */
b1 = test;
- /* move on until we find a vertex that is above the midline */
+ /* move on until we find a vertex that is above or on the midline */
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:
@@ -236,7 +232,10 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
a1 = test;
increment = 1;
} else {
-/* vertex is above midline */
+/* vertex is above or on midline
+ in the case where it is on the midline, we count that as above as well:
+ the vertex will be removed, and a new vertex will be created at the same
+ position */
#ifdef VORONOI_VERBOSE
message("First vertex is above midline (%g %g --> %g)!",
cell->vertices[0][0] + cell->x[0],
@@ -250,14 +249,12 @@ __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.) {
+ 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:
@@ -290,9 +287,9 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
#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
- midline */
+ we start from the vertex above (or on) the midline and search in the
+ direction opposite to the intersected edge direction until we find a vertex
+ below the midline */
/* we make sure we store the test result for the second vertex above the
midline as well, since we need this for intersection point computations
@@ -307,10 +304,9 @@ __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.) {
+ while (test >= 0.) {
/* make sure we always store the most recent test result */
a2 = test;
i += increment;
@@ -322,7 +318,6 @@ __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();
}
index_below2 = i;
@@ -343,6 +338,15 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
error("Only 1 intersected edge found!");
}
+ /* there is exactly one degenerate case we have not addressed yet: the case
+ where index_above1 and index_above2 are the same and are on the midline.
+ In this case we don't want to create 2 new vertices. Instead, we just keep
+ index_above1, which basically means nothing happens at all and we can just
+ return */
+ if (index_above1 == index_above2 && a1 == 0.) {
+ return;
+ }
+
/* 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
@@ -378,6 +382,10 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
cell->vertices[index_above2][1] + cell->x[1], index_above2, a2);
#endif
+ if (b1 == 0. || b2 == 0.) {
+ error("Vertex below midline is on midline!");
+ }
+
/* 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 */
@@ -400,6 +408,7 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
vertices[nvert][0] = cell->vertices[i][0];
vertices[nvert][1] = cell->vertices[i][1];
++nvert;
+ VORONOI_CHECK_SIZE();
++i;
if (i == cell->nvert) {
i = 0;
@@ -411,11 +420,13 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
vertices[nvert][1] = a1 * cell->vertices[index_below1][1] +
b1 * cell->vertices[index_above1][1];
++nvert;
+ VORONOI_CHECK_SIZE();
vertices[nvert][0] = a2 * cell->vertices[index_below2][0] +
b2 * cell->vertices[index_above2][0];
vertices[nvert][1] = a2 * cell->vertices[index_below2][1] +
b2 * cell->vertices[index_above2][1];
++nvert;
+ VORONOI_CHECK_SIZE();
/* overwrite the original vertices */
cell->nvert = nvert;
diff --git a/tests/testVoronoi2D.c b/tests/testVoronoi2D.c
index 47eb40bf0e0a922f75296c6d08e0d49448490bcb..308f83d221fb6419308d0aa620de19c179a44a3d 100644
--- a/tests/testVoronoi2D.c
+++ b/tests/testVoronoi2D.c
@@ -49,7 +49,7 @@ int main() {
assert(cell.centroid[1] == 0.5f);
}
- /* test interaction algorithm */
+ /* test interaction algorithm: normal case */
{
/* create 100 cells with random positions in [0,1]x[0,1] */
struct voronoi_cell cells[TESTVORONOI2D_NUMCELL];
@@ -110,5 +110,64 @@ int main() {
assert(fabs(Atot - 1.0f) < 1.e-6);
}
+ /* test interaction algorithm */
+ {
+ struct voronoi_cell cells[100];
+ double x[2];
+ float dx[2];
+ int i, j;
+ float Atot;
+ struct voronoi_cell *cell_i, *cell_j;
+
+ for (i = 0; i < 10; ++i) {
+ for (j = 0; j < 10; ++j) {
+ x[0] = (i + 0.5f) * 0.1;
+ x[1] = (j + 0.5f) * 0.1;
+ voronoi_cell_init(&cells[10 * i + j], x);
+ }
+ }
+
+ /* interact the cells (with periodic boundaries) */
+ for (i = 0; i < 100; ++i) {
+ cell_i = &cells[i];
+ for (j = 0; j < 100; ++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);
+ }
+ }
+ }
+
+ /* print the cells to the stdout */
+ for (i = 0; i < 100; ++i) {
+ voronoi_print_cell(&cells[i]);
+ voronoi_cell_finalize(&cells[i]);
+ Atot += cells[i].volume;
+ }
+
+ assert(fabs(Atot - 1.0f) < 1.e-6);
+ }
+
return 0;
}