Commit e07b2c04 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Added actual grid test to testVoronoi3D; degenerate case crashes. Started...

Added actual grid test to testVoronoi3D; degenerate case crashes. Started documenting 3D Voronoi algorithm to refresh my memory about how it works.
parent 0cf4064e
......@@ -56,12 +56,21 @@ __attribute__((always_inline)) INLINE int check_counter(int *counter,
return 1;
}
/* safewhile is a wrapper around a while that adds a unique counter variable to
the loop that is increased by 1 for each time the loop is executed, and
causes the code to crash if this number exceeds a given value.
We use this to quickly enable or disable number of iterations checks for a
large number of while loops */
#define safewhile(condition) \
int LOOPCOUNTER_NAME(__LINE__) = 0; \
while (check_counter(&LOOPCOUNTER_NAME(__LINE__), __LINE__) && (condition))
#else
#else // LOOP_CHECK
/* If LOOP_CHECK is not defined, safewhile and while are EXACTLY the same */
#define safewhile(condition) while (condition)
#endif
#endif // LOOP_CHECK
/* Tolerance parameter used to decide when to use more precise geometric
criteria */
......@@ -78,9 +87,12 @@ __attribute__((always_inline)) INLINE int check_counter(int *counter,
#define VORONOI3D_BOX_LEFT 18446744073709551604llu
#define VORONOI3D_BOX_RIGHT 18446744073709551605llu
/* Global (!) variables used to store the extents of the simulation box */
extern float global_voronoi_box_anchor[3];
extern float global_voronoi_box_side[3];
/* Macro used to declare the global variables in files that include this header
file */
#define VORONOI_DECLARE_GLOBAL_VARIABLES() \
float global_voronoi_box_anchor[3]; \
float global_voronoi_box_side[3];
......@@ -97,6 +109,12 @@ extern float global_voronoi_box_side[3];
#define VORONOI3D_BOX_SIDE_Y global_voronoi_box_side[1]
#define VORONOI3D_BOX_SIDE_Z global_voronoi_box_side[2]
/**
* @brief Store the extents of the simulation box in the global variables.
*
* @param anchor Corner of the simulation box with the lowest coordinate values.
* @param side Side lengths of the simulation box.
*/
__attribute__((always_inline)) INLINE static void voronoi_set_box(float *anchor,
float *side) {
global_voronoi_box_anchor[0] = anchor[0];
......@@ -107,10 +125,28 @@ __attribute__((always_inline)) INLINE static void voronoi_set_box(float *anchor,
global_voronoi_box_side[2] = side[2];
}
/*******************************************************************************
* 3D specific simulation box related methods
******************************************************************************/
/**
* @brief Get the volume of the simulation box stored in the global variables.
*
* This method is only used for debugging purposes.
*
* @return Volume of the simulation box as it is stored in the global variables.
*/
__attribute__((always_inline)) INLINE static float voronoi_get_box_volume() {
return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Y * VORONOI3D_BOX_SIDE_Z;
}
/**
* @brief Get the centroid of the simulation box stored in the global variables.
*
* This method is only used for debugging purposes.
*
* @param box_centroid Array to store the resulting coordinates in.
*/
__attribute__((always_inline)) INLINE static void voronoi_get_box_centroid(
float *box_centroid) {
box_centroid[0] = 0.5f * VORONOI3D_BOX_SIDE_X + VORONOI3D_BOX_ANCHOR_X;
......@@ -118,6 +154,17 @@ __attribute__((always_inline)) INLINE static void voronoi_get_box_centroid(
box_centroid[2] = 0.5f * VORONOI3D_BOX_SIDE_Z + VORONOI3D_BOX_ANCHOR_Z;
}
/**
* @brief Check if the given vertex (of the cell generated by the generator with
* the given position) is inside the simulation box.
*
* This method does not return a value, but causes the code to crash if the
* given vertex is found to be outside the simulation box.
*
* @param v Coordinates of a vertex, relative to the position of the cell
* generator.
* @param x Position of the cell generator.
*/
__attribute__((always_inline)) INLINE void voronoi_box_test_inside(float *v,
double *x) {
float vpos[3];
......@@ -139,6 +186,18 @@ __attribute__((always_inline)) INLINE void voronoi_box_test_inside(float *v,
}
}
/**
* @brief Get the surface area and coordinates of the face midpoint for the
* face of the simulation box with the given unique ID.
*
* This method is only used for debugging purposes.
*
* @param id Unique ID of one of the 6 faces of the simulation box.
* @param face_midpoint Array to store the coordinates of the requested
* midpoint in.
* @return Surface area of the face, or 0 if the given ID does not correspond to
* a known face of the simulation box.
*/
__attribute__((always_inline)) INLINE static float voronoi_get_box_face(
unsigned long long id, float *face_midpoint) {
......@@ -191,6 +250,15 @@ __attribute__((always_inline)) INLINE static float voronoi_get_box_face(
* http://math.lbl.gov/voro++/
******************************************************************************/
/**
* @brief Print the given cell to the stderr in a format that can be easily
* plotted using gnuplot.
*
* This method prints to the stderr instead of stdout to make it possible to use
* it right before crashing the code.
*
* @param c Voronoi cell to print.
*/
__attribute__((always_inline)) INLINE void voronoi_print_gnuplot_c(
struct voronoi_cell *c) {
......@@ -241,7 +309,12 @@ __attribute__((always_inline)) INLINE void voronoi_print_cell(
/**
* @brief Get the index of the vertex pointed to by the given edge of the given
* vertex
* vertex.
*
* @param c 3D Voronoi cell.
* @param vertex Index of a vertex of the cell.
* @param edge Edge of that vertex.
* @return Index of the vertex on the other side of the edge.
*/
__attribute__((always_inline)) INLINE int voronoi_get_edge(
struct voronoi_cell *c, int vertex, int edge) {
......@@ -249,8 +322,18 @@ __attribute__((always_inline)) INLINE int voronoi_get_edge(
}
/**
* @brief Get the index of vertex in the edge list of the vertex pointed to by
* the given edge of the given vertex
* @brief Get the index of the given edge in the edge list of the vertex on the
* other side of the edge of the given vertex.
*
* Suppose that the given vertex has edges [edge1, edge2, given_edge], and that
* the vertex on the other side of given_edge has edges [edge1, given_edge,
* edge2], then this method returns 1.
*
* @param c 3D Voronoi cell.
* @param vertex Index of a vertex of the cell.
* @param edge Edge of that vertex.
* @return Index of that edge in the edge list of the vertex on the other side
* of the edge.
*/
__attribute__((always_inline)) INLINE int voronoi_get_edgeindex(
struct voronoi_cell *c, int vertex, int edge) {
......@@ -258,8 +341,13 @@ __attribute__((always_inline)) INLINE int voronoi_get_edgeindex(
}
/**
* @brief Set the index of the vertex pointed to by the given edge of the given
* vertex
* @brief Set the index of the vertex on the other side of the given edge of the
* given vertex.
*
* @param c 3D Voronoi cell.
* @param vertex Index of a vertex of the cell.
* @param edge Edge of that vertex.
* @param value Index of the vertex on the other side of that edge.
*/
__attribute__((always_inline)) INLINE void voronoi_set_edge(
struct voronoi_cell *c, int vertex, int edge, int value) {
......@@ -267,8 +355,19 @@ __attribute__((always_inline)) INLINE void voronoi_set_edge(
}
/**
* @brief Set the index of vertex in the edge list of the vertex pointed to by
* the given edge of the given vertex
* @brief Set the index of the given edge in the edge list of the vertex on the
* other side of the edge of the given vertex.
*
* Suppose that the given vertex has edges [edge1, edge2, given_edge], and we
* want to tell this method that the vertex on the other side of given_edge has
* edges [edge1, given_edge, edge2], then we need to pass on a value of 1 to
* this method.
*
* @param c 3D Voronoi cell.
* @param vertex Index of a vertex of that cell.
* @param edge Edge of that vertex.
* @param value Index of that edge in the edge list of the vertex on the other
* side of the edge.
*/
__attribute__((always_inline)) INLINE void voronoi_set_edgeindex(
struct voronoi_cell *c, int vertex, int edge, int value) {
......@@ -276,7 +375,22 @@ __attribute__((always_inline)) INLINE void voronoi_set_edgeindex(
}
/**
* @brief Get the neighbour for the given edge of the given vertex
* @brief Get the neighbour for the given edge of the given vertex.
*
* An edge is shared by two faces, and each face has a neighbour. Luckily, each
* edge also has two endpoints, so we can get away with storing only one
* neighbour per endpoint of an edge. We have complete freedom in choosing which
* neighbour to store in which endpoint, but we need to be consistent about it.
* Here we use the following convention: if we take a vector pointing away from
* the given vertex along the given edge direction, then we store the neighbour
* that corresponds to the face to the right if looking to the cell from the
* outside. This is the face that you encounter first when rotating counter-
* clockwise around that vector, starting from outside the cell.
*
* @param c 3D Voronoi cell.
* @param vertex Index of a vertex of that cell.
* @param edge Edge of that vertex.
* @return Index of the neighbour corresponding to that edge and vertex.
*/
__attribute__((always_inline)) INLINE int voronoi_get_ngb(
struct voronoi_cell *c, int vertex, int edge) {
......@@ -284,7 +398,22 @@ __attribute__((always_inline)) INLINE int voronoi_get_ngb(
}
/**
* @brief Set the neighbour for the given edge of the given vertex
* @brief Set the neighbour for the given edge of the given vertex.
*
* An edge is shared by two faces, and each face has a neighbour. Luckily, each
* edge also has two endpoints, so we can get away with storing only one
* neighbour per endpoint of an edge. We have complete freedom in choosing which
* neighbour to store in which endpoint, but we need to be consistent about it.
* Here we use the following convention: if we take a vector pointing away from
* the given vertex along the given edge direction, then we store the neighbour
* that corresponds to the face to the right if looking to the cell from the
* outside. This is the face that you encounter first when rotating counter-
* clockwise around that vector, starting from outside the cell.
*
* @param c 3D Voronoi cell.
* @param vertex Index of a vertex of that cell.
* @param edge Edge of that vertex.
* @param value Index of the neighbour corresponding to that edge and vertex.
*/
__attribute__((always_inline)) INLINE void voronoi_set_ngb(
struct voronoi_cell *c, int vertex, int edge, int value) {
......@@ -292,7 +421,12 @@ __attribute__((always_inline)) INLINE void voronoi_set_ngb(
}
/**
* @brief Check if the 3D Voronoi cell is still consistent
* @brief Check if the 3D Voronoi cell is still consistent.
*
* A cell is consistent if its edges are consistent, i.e. if edge e of vertex v1
* points to vertex v2, then v2 should have an edge that points to v1 as well,
* and then the edge index of vertex v1 should contain the index of that edge
* in the edge list of v2.
*
* @param cell 3D Voronoi cell to check
*/
......@@ -321,6 +455,25 @@ __attribute__((always_inline)) INLINE void voronoi_check_cell_consistency(
/**
* @brief Check if the given vertex is above, below or on the cutting plane
* defined by the given parameters.
*
* @param v Coordinates of a cell vertex, relative w.r.t. the position of the
* generator of the cell.
* @param dx Half of the relative distance vector between the position of the
* generator of the cell and the position of the neighbouring cell that
* generates the cutting plane, pointing from the generator position to the
* cutting plane.
* @param r2 Squared length of dx.
* @param test Variable to store the result of the geometric test in, which
* corresponds to the projected distance between the generator and the vertex
* along dx.
* @param teststack Stack to store the results of the N last tests in (for
* debugging purposes only).
* @param teststack_size Next available field in the teststack, is reset to 0 if
* the teststack is full (so the N+1th results is overwritten; for debugging
* purposes only).
* @return Result of the test: -1 if the vertex is below the cutting plane, +1
* if it is above, and 0 if it is on the cutting plane.
*/
__attribute__((always_inline)) INLINE int voronoi_test_vertex(
float *v, float *dx, float r2, float *test, float *teststack,
......@@ -344,8 +497,9 @@ __attribute__((always_inline)) INLINE int voronoi_test_vertex(
}
/**
* @brief Initialize the cell as a cube with side 2*h, centered on the particle
* position
* @brief Initialize the cell as a cube that spans the entire simulation box.
*
* @param c 3D Voronoi cell to initialize.
*/
__attribute__((always_inline)) INLINE void voronoi_initialize(
struct voronoi_cell *cell) {
......@@ -525,9 +679,9 @@ __attribute__((always_inline)) INLINE void voronoi_initialize(
}
/**
* @brief Find an edge of the voronoi_cell that intersects the cutting plane
* @brief Find an edge of the voronoi_cell that intersects the cutting plane.
*
* @param c voronoi_cell
* @param c 3D Voronoi cell.
* @param dx vector pointing from the midpoint of the line segment between pi
* and pj to pj
* @param r2 Squared norm of dx
......
......@@ -107,6 +107,7 @@ int main() {
}
}
Atot = 0.0f;
/* print the cells to the stdout */
for (i = 0; i < TESTVORONOI2D_NUMCELL; ++i) {
printf("Cell %i:\n", i);
......@@ -184,6 +185,7 @@ int main() {
}
}
Atot = 0.0f;
/* print the cells to the stdout */
for (i = 0; i < 100; ++i) {
printf("Cell %i:\n", i);
......
......@@ -22,6 +22,8 @@
#include "hydro/Shadowswift/voronoi3d_algorithm.h"
#include "part.h"
#define TESTVORONOI3D_NUMCELL 100
/**
* @brief Check if voronoi_volume_tetrahedron() works
*/
......@@ -1209,90 +1211,183 @@ int main() {
/* Test the different paths */
test_paths();
/* Create a Voronoi cell */
double x[3] = {0.5f, 0.5f, 0.5f};
struct voronoi_cell cell;
voronoi_cell_init(&cell, x);
/* Interact with neighbours */
float x0[3] = {0.5f, 0.0f, 0.0f};
float x1[3] = {-0.5f, 0.0f, 0.0f};
float x2[3] = {0.0f, 0.5f, 0.0f};
float x3[3] = {0.0f, -0.5f, 0.0f};
float x4[3] = {0.0f, 0.0f, 0.5f};
float x5[3] = {0.0f, 0.0f, -0.5f};
voronoi_cell_interact(&cell, x0, 1);
voronoi_cell_interact(&cell, x1, 2);
voronoi_cell_interact(&cell, x2, 3);
voronoi_cell_interact(&cell, x3, 4);
voronoi_cell_interact(&cell, x4, 5);
voronoi_cell_interact(&cell, x5, 6);
float expected_midpoints[6][3], expected_areas[6];
expected_areas[0] = 0.25f;
expected_midpoints[0][0] = 0.25f;
expected_midpoints[0][1] = 0.5f;
expected_midpoints[0][2] = 0.5f;
expected_areas[1] = 0.25f;
expected_midpoints[1][0] = 0.75f;
expected_midpoints[1][1] = 0.5f;
expected_midpoints[1][2] = 0.5f;
expected_areas[2] = 0.25f;
expected_midpoints[2][0] = 0.5f;
expected_midpoints[2][1] = 0.25f;
expected_midpoints[2][2] = 0.5f;
expected_areas[3] = 0.25f;
expected_midpoints[3][0] = 0.5f;
expected_midpoints[3][1] = 0.75f;
expected_midpoints[3][2] = 0.5f;
expected_areas[4] = 0.25f;
expected_midpoints[4][0] = 0.5f;
expected_midpoints[4][1] = 0.5f;
expected_midpoints[4][2] = 0.25f;
expected_areas[5] = 0.25f;
expected_midpoints[5][0] = 0.5f;
expected_midpoints[5][1] = 0.5f;
expected_midpoints[5][2] = 0.75f;
/* Interact with some more neighbours to check if they are properly ignored */
float xE0[3] = {0.6f, 0.0f, 0.1f};
float xE1[3] = {-0.7f, 0.2f, 0.04f};
voronoi_cell_interact(&cell, xE0, 7);
voronoi_cell_interact(&cell, xE1, 8);
/* Finalize cell and check results */
voronoi_cell_finalize(&cell);
if (fabs(cell.volume - 0.125f) > 1.e-5) {
error("Wrong volume: %g!", cell.volume);
/* Test the interaction and geometry algorithms */
{
/* Create a Voronoi cell */
double x[3] = {0.5f, 0.5f, 0.5f};
struct voronoi_cell cell;
voronoi_cell_init(&cell, x);
/* Interact with neighbours */
float x0[3] = {0.5f, 0.0f, 0.0f};
float x1[3] = {-0.5f, 0.0f, 0.0f};
float x2[3] = {0.0f, 0.5f, 0.0f};
float x3[3] = {0.0f, -0.5f, 0.0f};
float x4[3] = {0.0f, 0.0f, 0.5f};
float x5[3] = {0.0f, 0.0f, -0.5f};
voronoi_cell_interact(&cell, x0, 1);
voronoi_cell_interact(&cell, x1, 2);
voronoi_cell_interact(&cell, x2, 3);
voronoi_cell_interact(&cell, x3, 4);
voronoi_cell_interact(&cell, x4, 5);
voronoi_cell_interact(&cell, x5, 6);
float expected_midpoints[6][3], expected_areas[6];
expected_areas[0] = 0.25f;
expected_midpoints[0][0] = 0.25f;
expected_midpoints[0][1] = 0.5f;
expected_midpoints[0][2] = 0.5f;
expected_areas[1] = 0.25f;
expected_midpoints[1][0] = 0.75f;
expected_midpoints[1][1] = 0.5f;
expected_midpoints[1][2] = 0.5f;
expected_areas[2] = 0.25f;
expected_midpoints[2][0] = 0.5f;
expected_midpoints[2][1] = 0.25f;
expected_midpoints[2][2] = 0.5f;
expected_areas[3] = 0.25f;
expected_midpoints[3][0] = 0.5f;
expected_midpoints[3][1] = 0.75f;
expected_midpoints[3][2] = 0.5f;
expected_areas[4] = 0.25f;
expected_midpoints[4][0] = 0.5f;
expected_midpoints[4][1] = 0.5f;
expected_midpoints[4][2] = 0.25f;
expected_areas[5] = 0.25f;
expected_midpoints[5][0] = 0.5f;
expected_midpoints[5][1] = 0.5f;
expected_midpoints[5][2] = 0.75f;
/* Interact with some more neighbours to check if they are properly
ignored */
float xE0[3] = {0.6f, 0.0f, 0.1f};
float xE1[3] = {-0.7f, 0.2f, 0.04f};
voronoi_cell_interact(&cell, xE0, 7);
voronoi_cell_interact(&cell, xE1, 8);
/* Finalize cell and check results */
voronoi_cell_finalize(&cell);
if (fabs(cell.volume - 0.125f) > 1.e-5) {
error("Wrong volume: %g!", cell.volume);
}
if (fabs(cell.centroid[0] - 0.5f) > 1.e-5f ||
fabs(cell.centroid[1] - 0.5f) > 1.e-5f ||
fabs(cell.centroid[2] - 0.5f) > 1.e-5f) {
error("Wrong centroid: %g %g %g!", cell.centroid[0], cell.centroid[1],
cell.centroid[2]);
}
/* Check faces. */
float A, midpoint[3];
for (int i = 0; i < 6; ++i) {
A = voronoi_get_face(&cell, i + 1, midpoint);
if (A) {
if (fabs(A - expected_areas[i]) > 1.e-5) {
error("Wrong surface area: %g!", A);
}
if (fabs(midpoint[0] - expected_midpoints[i][0] + cell.x[0]) > 1.e-5 ||
fabs(midpoint[1] - expected_midpoints[i][1] + cell.x[1]) > 1.e-5 ||
fabs(midpoint[2] - expected_midpoints[i][2] + cell.x[2]) > 1.e-5) {
error("Wrong face midpoint: %g %g %g (should be %g %g %g)!",
midpoint[0], midpoint[1], midpoint[2], expected_midpoints[i][0],
expected_midpoints[i][1], expected_midpoints[i][2]);
}
} else {
error("Neighbour %i not found!", i);
}
}
}
if (fabs(cell.centroid[0] - 0.5f) > 1.e-5f ||
fabs(cell.centroid[1] - 0.5f) > 1.e-5f ||
fabs(cell.centroid[2] - 0.5f) > 1.e-5f) {
error("Wrong centroid: %g %g %g!", cell.centroid[0], cell.centroid[1],
cell.centroid[2]);
/* Test degenerate cases */
test_degeneracies();
/* Construct a small random grid */
{
int i, j;
double x[3];
float dx[3];
float Vtot;
struct voronoi_cell cells[TESTVORONOI3D_NUMCELL];
struct voronoi_cell *cell_i, *cell_j;
/* initialize cells with random generator locations */
for (i = 0; i < TESTVORONOI3D_NUMCELL; ++i) {
x[0] = ((double)rand()) / ((double)RAND_MAX);
x[1] = ((double)rand()) / ((double)RAND_MAX);
x[2] = ((double)rand()) / ((double)RAND_MAX);
voronoi_cell_init(&cells[i], x);
}
/* interact the cells */
for (i = 0; i < TESTVORONOI3D_NUMCELL; ++i) {
cell_i = &cells[i];
for (j = 0; j < TESTVORONOI3D_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];
dx[2] = cell_i->x[2] - cell_j->x[2];
voronoi_cell_interact(cell_i, dx, j);
}
}
}
Vtot = 0.0f;
/* print the cells to the stdout */
for (i = 0; i < TESTVORONOI3D_NUMCELL; ++i) {
/* voronoi_print_gnuplot_c(&cells[i]);*/
voronoi_cell_finalize(&cells[i]);
Vtot += cells[i].volume;
}
assert(fabs(Vtot - 1.0f) < 1.e-6);
}
/* Check faces. */
float A, midpoint[3];
for (int i = 0; i < 6; ++i) {
A = voronoi_get_face(&cell, i + 1, midpoint);
if (A) {
if (fabs(A - expected_areas[i]) > 1.e-5) {
error("Wrong surface area: %g!", A);
/* Construct a small Cartesian grid full of degeneracies */
{
int i, j, k;
double x[3];
float dx[3];
float Vtot;
struct voronoi_cell cells[1000];
struct voronoi_cell *cell_i, *cell_j;
/* initialize cells with random generator locations */
for (i = 0; i < 10; ++i) {
for (j = 0; j < 10; ++j) {
for (k = 0; k < 10; ++k) {
x[0] = (i + 0.5f) * 0.1;
x[1] = (j + 0.5f) * 0.1;
x[2] = (k + 0.5f) * 0.1;
voronoi_cell_init(&cells[i], x);
}
}
if (fabs(midpoint[0] - expected_midpoints[i][0] + cell.x[0]) > 1.e-5 ||
fabs(midpoint[1] - expected_midpoints[i][1] + cell.x[1]) > 1.e-5 ||
fabs(midpoint[2] - expected_midpoints[i][2] + cell.x[2]) > 1.e-5) {
error("Wrong face midpoint: %g %g %g (should be %g %g %g)!",
midpoint[0], midpoint[1], midpoint[2], expected_midpoints[i][0],
expected_midpoints[i][1], expected_midpoints[i][2]);
}
/* interact the cells */
for (i = 0; i < TESTVORONOI3D_NUMCELL; ++i) {
cell_i = &cells[i];
for (j = 0; j < TESTVORONOI3D_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];
dx[2] = cell_i->x[2] - cell_j->x[2];
voronoi_cell_interact(cell_i, dx, j);
}
}
} else {
error("Neighbour %i not found!", i);
}
}
test_degeneracies();
Vtot = 0.0f;
/* print the cells to the stdout */
for (i = 0; i < TESTVORONOI3D_NUMCELL; ++i) {
voronoi_print_gnuplot_c(&cells[i]);
voronoi_cell_finalize(&cells[i]);
Vtot += cells[i].volume;
}
assert(fabs(Vtot - 1.0f) < 1.e-6);
}
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment