Skip to content
Snippets Groups Projects
Commit 77e62725 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Added support for degeneracies to 2D Voronoi algorithm.

parent 51d1c406
No related branches found
No related tags found
1 merge request!3211D and 2D moving mesh algorithm
Hide whitespace changes
Inline Side-by-side
src/hydro/Shadowswift/voronoi2d_algorithm.h
+ 30
− 19
View file @ 77e62725
......@@ -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;
......
tests/testVoronoi2D.c
+ 60
− 1
View file @ 77e62725
......@@ -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;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment