diff --git a/examples/main.c b/examples/main.c
index 41fbd8aba892cc1677c3317fa3579ca6ee5dd7fa..8d4289c872ce34ba626cacffca4097c0093fe0c3 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -101,10 +101,6 @@ void print_help_message() {
"parameter file.\n");
}
-#if defined(SHADOWFAX_SPH)
-VORONOI_DECLARE_GLOBAL_VARIABLES()
-#endif
-
/**
* @brief Main routine that loads a few particles and generates some output.
*
@@ -390,13 +386,6 @@ int main(int argc, char *argv[]) {
struct gravity_props gravity_properties;
if (with_self_gravity) gravity_props_init(&gravity_properties, params);
-#if defined(SHADOWFAX_SPH)
- /* Override the variables governing the density iteration
- (see src/hydro/Shadowswift/hydro.h for full explanation) */
- hydro_properties.target_neighbours = 1.0f;
- hydro_properties.delta_neighbours = 0.0f;
-#endif
-
/* Read particles and space information from (GADGET) ICs */
char ICfileName[200] = "";
parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
@@ -439,27 +428,6 @@ int main(int argc, char *argv[]) {
fflush(stdout);
}
-#if defined(SHADOWFAX_SPH)
- /* set the *global* box dimensions */
- float box_anchor[3], box_side[3];
- if (periodic) {
- box_anchor[0] = -0.5f * dim[0];
- box_anchor[1] = -0.5f * dim[1];
- box_anchor[2] = -0.5f * dim[2];
- box_side[0] = 2.0f * dim[0];
- box_side[1] = 2.0f * dim[1];
- box_side[2] = 2.0f * dim[2];
- } else {
- box_anchor[0] = 0.0f;
- box_anchor[1] = 0.0f;
- box_anchor[2] = 0.0f;
- box_side[0] = dim[0];
- box_side[1] = dim[1];
- box_side[2] = dim[2];
- }
- voronoi_set_box(box_anchor, box_side);
-#endif
-
#ifdef SWIFT_DEBUG_CHECKS
/* Check once and for all that we don't have unwanted links */
if (!with_stars) {
diff --git a/src/Makefile.am b/src/Makefile.am
index 29c110218416cc4bece0d766516eeb0a97fe4810..a29cf5521f2b25cbcc91ef85e277d1340ee69dd5 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -45,7 +45,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
physical_constants.h physical_constants_cgs.h potential.h version.h \
hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
- dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h
+ dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
+ hydro_space.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -55,7 +56,8 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
physical_constants.c potential.c hydro_properties.c \
runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
- part_type.c xmf.c gravity_properties.c gravity.c
+ part_type.c xmf.c gravity_properties.c gravity.c \
+ hydro_space.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index a614d08c30b21f9e7d422bf6b6a09d10d2e89799..051c22f46b3ecdff5d3de910e0f75409b0e78f02 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -22,6 +22,7 @@
#include "adiabatic_index.h"
#include "approx_math.h"
#include "equation_of_state.h"
+#include "hydro_space.h"
#include "minmax.h"
#include <float.h>
@@ -165,9 +166,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* the variaous density tasks
*
* @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
- struct part *restrict p) {
+ struct part *restrict p, const struct hydro_space *hs) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
@@ -400,7 +402,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->u_full = p->u;
hydro_reset_acceleration(p);
- hydro_init_part(p);
+ hydro_init_part(p, NULL);
}
#endif /* SWIFT_DEFAULT_HYDRO_H */
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index cc7b422ccbe7c678969df5779a4d4a054c65528e..747c81a8e64c18a06b04160cfab326a3521c5901 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -36,6 +36,7 @@
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
+#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
@@ -169,9 +170,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* the variaous density tasks
*
* @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
- struct part *restrict p) {
+ struct part *restrict p, const struct hydro_space *hs) {
p->rho = 0.f;
p->density.wcount = 0.f;
@@ -456,7 +458,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->entropy_full = p->entropy;
hydro_reset_acceleration(p);
- hydro_init_part(p);
+ hydro_init_part(p, NULL);
}
#endif /* SWIFT_GADGET2_HYDRO_H */
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 643489912a6c6b1db921e73b508910cc670d49ae..60ff8ccee0e1a1e9c3477a10293f8981ff9b837e 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -23,6 +23,7 @@
#include "approx_math.h"
#include "equation_of_state.h"
#include "hydro_gradients.h"
+#include "hydro_space.h"
#include "minmax.h"
#include "riemann.h"
@@ -145,9 +146,10 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
* Simply makes sure all necessary variables are initialized to zero.
*
* @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
- struct part* p) {
+ struct part* p, const struct hydro_space* hs) {
p->density.wcount = 0.0f;
p->density.wcount_dh = 0.0f;
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 56078a82569fb0bc30347d5c01831e9eecfd48f4..8f216a550ae061d01a594ff23d57575e754f85dc 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -38,6 +38,7 @@
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
+#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
@@ -183,9 +184,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* density sub-structure of a particle get zeroed in here.
*
* @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
- struct part *restrict p) {
+ struct part *restrict p, const struct hydro_space *hs) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
@@ -429,7 +431,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->u_full = p->u;
hydro_reset_acceleration(p);
- hydro_init_part(p);
+ hydro_init_part(p, NULL);
}
#endif /* SWIFT_MINIMAL_HYDRO_H */
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 20238896f1458d0abebacca4865968a3a671c886..4c4868cd3703e5ec5466d4878749a61284b19344 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -36,6 +36,7 @@
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
+#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
@@ -169,9 +170,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* the variaous density tasks
*
* @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
- struct part *restrict p) {
+ struct part *restrict p, const struct hydro_space *hs) {
p->rho = 0.f;
p->rho_bar = 0.f;
@@ -474,7 +476,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->v_full[2] = p->v[2];
hydro_reset_acceleration(p);
- hydro_init_part(p);
+ hydro_init_part(p, NULL);
}
#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 82db9f4b18be7dcdc9ad735fa9696d36b14ca5d1..0568d47ee7ed33c59790cbca943cccbf1ceda58f 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -22,6 +22,7 @@
#include "approx_math.h"
#include "equation_of_state.h"
#include "hydro_gradients.h"
+#include "hydro_space.h"
#include "voronoi_algorithm.h"
/**
@@ -128,14 +129,15 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
* Initializes the Voronoi cell.
*
* @param p The particle to act upon
+ * @param hs #hydro_space containing extra information about the space.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
- struct part* p) {
+ struct part* p, const struct hydro_space* hs) {
p->density.wcount = 0.0f;
p->density.wcount_dh = 0.0f;
- voronoi_cell_init(&p->cell, p->x);
+ voronoi_cell_init(&p->cell, p->x, hs->anchor, hs->side);
/* Set the active flag to active. */
p->force.active = 1;
diff --git a/src/hydro/Shadowswift/voronoi1d_algorithm.h b/src/hydro/Shadowswift/voronoi1d_algorithm.h
index 31045895be31327e1e968f953cc66c137bfb3384..74cc5f1dbf3a2d72df55ce73de0321b5493193a7 100644
--- a/src/hydro/Shadowswift/voronoi1d_algorithm.h
+++ b/src/hydro/Shadowswift/voronoi1d_algorithm.h
@@ -20,15 +20,12 @@
#ifndef SWIFT_VORONOIXD_ALGORITHM_H
#define SWIFT_VORONOIXD_ALGORITHM_H
-#include <float.h>
#include <math.h>
#include <stdlib.h>
#include "error.h"
#include "inline.h"
#include "voronoi1d_cell.h"
-#define VORONOI_DECLARE_GLOBAL_VARIABLES()
-
/**
* @brief Store the extents of the simulation box in the global variables.
*
@@ -47,12 +44,15 @@ __attribute__((always_inline)) INLINE static void voronoi_set_box(
*
* @param cell 1D Voronoi cell to initialize.
* @param x Position of the generator of the cell.
+ * @param anchor Anchor of the simulation box.
+ * @param side Side lengths of the simulation box.
*/
__attribute__((always_inline)) INLINE void voronoi_cell_init(
- struct voronoi_cell *cell, const double *x) {
+ struct voronoi_cell *cell, const double *x, const double *anchor,
+ const double *side) {
cell->x = x[0];
- cell->xL = -DBL_MAX;
- cell->xR = DBL_MAX;
+ cell->xL = anchor[0] - cell->x;
+ cell->xR = anchor[0] + side[0] - cell->x;
cell->idL = 0;
cell->idR = 0;
cell->volume = 0.0f;
diff --git a/src/hydro/Shadowswift/voronoi2d_algorithm.h b/src/hydro/Shadowswift/voronoi2d_algorithm.h
index a799057ed6d56ac8c960b576105094a7bb30664b..ad4d565fbf817da7b7116a0c619a93b5aee1917c 100644
--- a/src/hydro/Shadowswift/voronoi2d_algorithm.h
+++ b/src/hydro/Shadowswift/voronoi2d_algorithm.h
@@ -43,48 +43,17 @@
#define VORONOI2D_BOX_TOP 18446744073709551604llu
#define VORONOI2D_BOX_BOTTOM 18446744073709551605llu
-/* Global variables used to store the extent of the simulation box */
-extern float global_voronoi_box_anchor[2];
-extern float global_voronoi_box_side[2];
-
-/* Macro used to declare the global variables in any piece of code that uses
- this header file */
-#define VORONOI_DECLARE_GLOBAL_VARIABLES() \
- float global_voronoi_box_anchor[2]; \
- float global_voronoi_box_side[2];
-
-/* Macros that make it easier to change the box extents that are used throughout
- the code */
-#define VORONOI2D_BOX_ANCHOR_X global_voronoi_box_anchor[0]
-#define VORONOI2D_BOX_ANCHOR_Y global_voronoi_box_anchor[1]
-#define VORONOI2D_BOX_SIDE_X global_voronoi_box_side[0]
-#define VORONOI2D_BOX_SIDE_Y global_voronoi_box_side[1]
-
-/**
- * @brief Set the values of the global variables that store the simulation box
- * extents.
- *
- * @param anchor Corner of the box with the lowest coordinate values.
- * @param side Side lengths of the box.
- */
-__attribute__((always_inline)) INLINE static void voronoi_set_box(
- const float *anchor, const float *side) {
-
- global_voronoi_box_anchor[0] = anchor[0];
- global_voronoi_box_anchor[1] = anchor[1];
-
- global_voronoi_box_side[0] = side[0];
- global_voronoi_box_side[1] = side[1];
-}
-
/**
* @brief Initialize a 2D Voronoi cell.
*
* @param cell 2D Voronoi cell to initialize.
* @param x Position of the generator of the cell.
+ * @param anchor Anchor of the simulation box containing all particles.
+ * @param side Side lengths of the simulation box containing all particles.
*/
__attribute__((always_inline)) INLINE void voronoi_cell_init(
- struct voronoi_cell *cell, const double *x) {
+ struct voronoi_cell *cell, const double *x, const double *anchor,
+ const double *side) {
/* Set the position of the generator of the cell (for reference) */
cell->x[0] = x[0];
@@ -94,21 +63,17 @@ __attribute__((always_inline)) INLINE void voronoi_cell_init(
(note: all vertex coordinates are relative w.r.t. the cell generator) */
cell->nvert = 4;
- cell->vertices[0][0] = VORONOI2D_BOX_ANCHOR_X - cell->x[0];
- cell->vertices[0][1] = VORONOI2D_BOX_ANCHOR_Y - cell->x[1];
+ cell->vertices[0][0] = anchor[0] - cell->x[0];
+ cell->vertices[0][1] = anchor[1] - cell->x[1];
- cell->vertices[1][0] = VORONOI2D_BOX_ANCHOR_X - cell->x[0];
- cell->vertices[1][1] =
- VORONOI2D_BOX_ANCHOR_Y + VORONOI2D_BOX_SIDE_Y - cell->x[1];
+ cell->vertices[1][0] = anchor[0] - cell->x[0];
+ cell->vertices[1][1] = anchor[1] + side[1] - cell->x[1];
- cell->vertices[2][0] =
- VORONOI2D_BOX_ANCHOR_X + VORONOI2D_BOX_SIDE_X - cell->x[0];
- cell->vertices[2][1] =
- VORONOI2D_BOX_ANCHOR_Y + VORONOI2D_BOX_SIDE_Y - cell->x[1];
+ cell->vertices[2][0] = anchor[0] + side[0] - cell->x[0];
+ cell->vertices[2][1] = anchor[1] + side[1] - cell->x[1];
- cell->vertices[3][0] =
- VORONOI2D_BOX_ANCHOR_X + VORONOI2D_BOX_SIDE_X - cell->x[0];
- cell->vertices[3][1] = VORONOI2D_BOX_ANCHOR_Y - cell->x[1];
+ cell->vertices[3][0] = anchor[0] + side[0] - cell->x[0];
+ cell->vertices[3][1] = anchor[1] - cell->x[1];
/* The neighbours are ordered such that neighbour i shares the face in between
vertices i and i+1 (with last vertex + 1 = first vertex)
diff --git a/src/hydro/Shadowswift/voronoi3d_algorithm.h b/src/hydro/Shadowswift/voronoi3d_algorithm.h
index 318c83bbd56535fd0d6f313323fbd97a7f846857..961f6e91ffd24a4d1336b985db004ffde562b8b6 100644
--- a/src/hydro/Shadowswift/voronoi3d_algorithm.h
+++ b/src/hydro/Shadowswift/voronoi3d_algorithm.h
@@ -91,162 +91,6 @@ __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];
-
-/* Bottom front left corner and side lengths of the large box that contains all
- particles and is used as initial cell at the start of the construction */
-/* We should make sure that this box is either so large a particle can never
- fall outside (by using FLT_MAX if that works), or is initialized to be larger
- than the (periodic) simulation box */
-#define VORONOI3D_BOX_ANCHOR_X global_voronoi_box_anchor[0]
-#define VORONOI3D_BOX_ANCHOR_Y global_voronoi_box_anchor[1]
-#define VORONOI3D_BOX_ANCHOR_Z global_voronoi_box_anchor[2]
-#define VORONOI3D_BOX_SIDE_X global_voronoi_box_side[0]
-#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(
- const float *anchor, const float *side) {
- global_voronoi_box_anchor[0] = anchor[0];
- global_voronoi_box_anchor[1] = anchor[1];
- global_voronoi_box_anchor[2] = anchor[2];
- global_voronoi_box_side[0] = side[0];
- global_voronoi_box_side[1] = side[1];
- 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;
- box_centroid[1] = 0.5f * VORONOI3D_BOX_SIDE_Y + VORONOI3D_BOX_ANCHOR_Y;
- 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(
- const float *v, const double *x) {
- float vpos[3];
- vpos[0] = x[0] + v[0];
- vpos[1] = x[1] + v[1];
- vpos[2] = x[2] + v[2];
- int check = (vpos[0] > VORONOI3D_BOX_ANCHOR_X - 1.e-5);
- check &= (vpos[0] < VORONOI3D_BOX_ANCHOR_X + VORONOI3D_BOX_SIDE_X + 1.e-5);
- check &= (vpos[1] > VORONOI3D_BOX_ANCHOR_Y - 1.e-5);
- check &= (vpos[1] < VORONOI3D_BOX_ANCHOR_Y + VORONOI3D_BOX_SIDE_Y + 1.e-5);
- check &= (vpos[2] > VORONOI3D_BOX_ANCHOR_Z - 1.e-5);
- check &= (vpos[2] < VORONOI3D_BOX_ANCHOR_Z + VORONOI3D_BOX_SIDE_Z + 1.e-5);
-
- if (!check) {
- error("New vertex outside box! (%g %g %g, [%g %g %g], [%g %g %g])", vpos[0],
- vpos[1], vpos[2], VORONOI3D_BOX_ANCHOR_X, VORONOI3D_BOX_ANCHOR_Y,
- VORONOI3D_BOX_ANCHOR_Z, VORONOI3D_BOX_SIDE_X, VORONOI3D_BOX_SIDE_Y,
- VORONOI3D_BOX_SIDE_Z);
- }
-}
-
-/**
- * @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) {
-
- if (id == VORONOI3D_BOX_FRONT) {
- face_midpoint[0] = 0.5f * VORONOI3D_BOX_SIDE_X + VORONOI3D_BOX_ANCHOR_X;
- face_midpoint[1] = VORONOI3D_BOX_ANCHOR_Y;
- face_midpoint[2] = 0.5f * VORONOI3D_BOX_SIDE_Z + VORONOI3D_BOX_ANCHOR_Z;
- return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Z;
- }
- if (id == VORONOI3D_BOX_BACK) {
- face_midpoint[0] = 0.5f * VORONOI3D_BOX_SIDE_X + VORONOI3D_BOX_ANCHOR_X;
- face_midpoint[1] = VORONOI3D_BOX_ANCHOR_Y + VORONOI3D_BOX_SIDE_Y;
- face_midpoint[2] = 0.5f * VORONOI3D_BOX_SIDE_Z + VORONOI3D_BOX_ANCHOR_Z;
- return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Z;
- }
-
- if (id == VORONOI3D_BOX_BOTTOM) {
- face_midpoint[0] = 0.5f * VORONOI3D_BOX_SIDE_X + VORONOI3D_BOX_ANCHOR_X;
- face_midpoint[1] = 0.5f * VORONOI3D_BOX_SIDE_Y + VORONOI3D_BOX_ANCHOR_Y;
- face_midpoint[2] = VORONOI3D_BOX_ANCHOR_Z;
- return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Y;
- }
- if (id == VORONOI3D_BOX_TOP) {
- face_midpoint[0] = 0.5f * VORONOI3D_BOX_SIDE_X + VORONOI3D_BOX_ANCHOR_X;
- face_midpoint[1] = 0.5f * VORONOI3D_BOX_SIDE_Y + VORONOI3D_BOX_ANCHOR_Y;
- face_midpoint[2] = VORONOI3D_BOX_ANCHOR_Z + VORONOI3D_BOX_SIDE_Z;
- return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Y;
- }
-
- if (id == VORONOI3D_BOX_LEFT) {
- face_midpoint[0] = VORONOI3D_BOX_ANCHOR_X;
- face_midpoint[1] = 0.5f * VORONOI3D_BOX_SIDE_Y + VORONOI3D_BOX_ANCHOR_Y;
- face_midpoint[2] = 0.5f * VORONOI3D_BOX_SIDE_Z + VORONOI3D_BOX_ANCHOR_Z;
- return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Y;
- }
- if (id == VORONOI3D_BOX_RIGHT) {
- face_midpoint[0] = VORONOI3D_BOX_ANCHOR_X + VORONOI3D_BOX_SIDE_X;
- face_midpoint[1] = 0.5f * VORONOI3D_BOX_SIDE_Y + VORONOI3D_BOX_ANCHOR_Y;
- face_midpoint[2] = 0.5f * VORONOI3D_BOX_SIDE_Z + VORONOI3D_BOX_ANCHOR_Z;
- return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Y;
- }
-
- return 0.0f;
-}
-
/*******************************************************************************
* 3D specific methods
*
@@ -518,63 +362,53 @@ __attribute__((always_inline)) INLINE int voronoi_test_vertex(
* @brief Initialize the cell as a cube that spans the entire simulation box.
*
* @param c 3D Voronoi cell to initialize.
+ * @param anchor Anchor of the simulation box.
+ * @param side Side lengths of the simulation box.
*/
__attribute__((always_inline)) INLINE void voronoi_initialize(
- struct voronoi_cell *cell) {
+ struct voronoi_cell *cell, const double *anchor, const double *side) {
cell->nvert = 8;
/* (0, 0, 0) -- 0 */
- cell->vertices[0] = VORONOI3D_BOX_ANCHOR_X - cell->x[0];
- cell->vertices[1] = VORONOI3D_BOX_ANCHOR_Y - cell->x[1];
- cell->vertices[2] = VORONOI3D_BOX_ANCHOR_Z - cell->x[2];
+ cell->vertices[0] = anchor[0] - cell->x[0];
+ cell->vertices[1] = anchor[1] - cell->x[1];
+ cell->vertices[2] = anchor[2] - cell->x[2];
/* (0, 0, 1)-- 1 */
- cell->vertices[3] = VORONOI3D_BOX_ANCHOR_X - cell->x[0];
- cell->vertices[4] = VORONOI3D_BOX_ANCHOR_Y - cell->x[1];
- cell->vertices[5] =
- VORONOI3D_BOX_ANCHOR_Z + VORONOI3D_BOX_SIDE_Z - cell->x[2];
+ cell->vertices[3] = anchor[0] - cell->x[0];
+ cell->vertices[4] = anchor[1] - cell->x[1];
+ cell->vertices[5] = anchor[2] + side[2] - cell->x[2];
/* (0, 1, 0) -- 2 */
- cell->vertices[6] = VORONOI3D_BOX_ANCHOR_X - cell->x[0];
- cell->vertices[7] =
- VORONOI3D_BOX_ANCHOR_Y + VORONOI3D_BOX_SIDE_Y - cell->x[1];
- cell->vertices[8] = VORONOI3D_BOX_ANCHOR_Z - cell->x[2];
+ cell->vertices[6] = anchor[0] - cell->x[0];
+ cell->vertices[7] = anchor[1] + side[1] - cell->x[1];
+ cell->vertices[8] = anchor[2] - cell->x[2];
/* (0, 1, 1) -- 3 */
- cell->vertices[9] = VORONOI3D_BOX_ANCHOR_X - cell->x[0];
- cell->vertices[10] =
- VORONOI3D_BOX_ANCHOR_Y + VORONOI3D_BOX_SIDE_Y - cell->x[1];
- cell->vertices[11] =
- VORONOI3D_BOX_ANCHOR_Z + VORONOI3D_BOX_SIDE_Z - cell->x[2];
+ cell->vertices[9] = anchor[0] - cell->x[0];
+ cell->vertices[10] = anchor[1] + side[1] - cell->x[1];
+ cell->vertices[11] = anchor[2] + side[2] - cell->x[2];
/* (1, 0, 0) -- 4 */
- cell->vertices[12] =
- VORONOI3D_BOX_ANCHOR_X + VORONOI3D_BOX_SIDE_X - cell->x[0];
- cell->vertices[13] = VORONOI3D_BOX_ANCHOR_Y - cell->x[1];
- cell->vertices[14] = VORONOI3D_BOX_ANCHOR_Z - cell->x[2];
+ cell->vertices[12] = anchor[0] + side[0] - cell->x[0];
+ cell->vertices[13] = anchor[1] - cell->x[1];
+ cell->vertices[14] = anchor[2] - cell->x[2];
/* (1, 0, 1) -- 5 */
- cell->vertices[15] =
- VORONOI3D_BOX_ANCHOR_X + VORONOI3D_BOX_SIDE_X - cell->x[0];
- cell->vertices[16] = VORONOI3D_BOX_ANCHOR_Y - cell->x[1];
- cell->vertices[17] =
- VORONOI3D_BOX_ANCHOR_Z + VORONOI3D_BOX_SIDE_Z - cell->x[2];
+ cell->vertices[15] = anchor[0] + side[0] - cell->x[0];
+ cell->vertices[16] = anchor[1] - cell->x[1];
+ cell->vertices[17] = anchor[2] + side[2] - cell->x[2];
/* (1, 1, 0) -- 6 */
- cell->vertices[18] =
- VORONOI3D_BOX_ANCHOR_X + VORONOI3D_BOX_SIDE_X - cell->x[0];
- cell->vertices[19] =
- VORONOI3D_BOX_ANCHOR_Y + VORONOI3D_BOX_SIDE_Y - cell->x[1];
- cell->vertices[20] = VORONOI3D_BOX_ANCHOR_Z - cell->x[2];
+ cell->vertices[18] = anchor[0] + side[0] - cell->x[0];
+ cell->vertices[19] = anchor[1] + side[1] - cell->x[1];
+ cell->vertices[20] = anchor[2] - cell->x[2];
/* (1, 1, 1) -- 7 */
- cell->vertices[21] =
- VORONOI3D_BOX_ANCHOR_X + VORONOI3D_BOX_SIDE_X - cell->x[0];
- cell->vertices[22] =
- VORONOI3D_BOX_ANCHOR_Y + VORONOI3D_BOX_SIDE_Y - cell->x[1];
- cell->vertices[23] =
- VORONOI3D_BOX_ANCHOR_Z + VORONOI3D_BOX_SIDE_Z - cell->x[2];
+ cell->vertices[21] = anchor[0] + side[0] - cell->x[0];
+ cell->vertices[22] = anchor[1] + side[1] - cell->x[1];
+ cell->vertices[23] = anchor[2] + side[2] - cell->x[2];
cell->orders[0] = 3;
cell->orders[1] = 3;
@@ -1039,13 +873,14 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
this point. */
/* first of all, we need to find a vertex which has edges that extend below
- the plane (since the remainder of our algorithm depends on that)
+ the plane (since the remainder of our algorithm depends on that). This is
+ not necessarily the case: in principle a vertex can only have edges that
+ extend inside or above the plane.
we create a stack of vertices to test (we use dstack for this), and add
vertex up. For each vertex on the stack, we then traverse its edges. If
the edge extends above the plane, we ignore it. If it extends below, we
stop. If the edge lies in the plane, we add the vertex on the other end
- to the stack (we also add vertices above the plane, although that should
- not be strictly necessary)
+ to the stack.
We make sure that up contains the index of a vertex extending beyond the
plane on exit. */
dstack[dstack_size] = up;
@@ -1061,7 +896,8 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
if (lw == -1) {
/* jump out of the for loop */
break;
- } else {
+ }
+ if (lw == 0) {
/* only add each vertex to the stack once */
k = 0;
safewhile(k < dstack_size && dstack[k] != lp) { ++k; }
@@ -1075,15 +911,24 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
}
/* we increased j after lw was calculated, so only the value of lw should be
- used to determine whether or not the loop was succesfull */
+ used to determine whether or not the loop was successful */
if (lw != -1) {
+ /* we did not find an edge that extends below the plane. There are two
+ possible reasons for this: either all vertices of the cell lie above
+ or inside the midplane of the segment connecting a point inside the
+ cell (the generator) with a point inside or outside the cell (the
+ neighbour). This is geometrically absurd.
+ Another reason might be that somehow all vertices in the midplane only
+ have edges that extend outwards. This is contradictory to the fact that
+ a Voronoi cell is convex, and therefore also unacceptable.
+ We conclude that we should NEVER end up here. */
voronoi_print_cell(c);
- error("Unable to find a vertex below the plane!");
+ error("Unable to find a vertex below the midplane!");
}
/* reset the delete stack, we need it later on */
dstack_size = 0;
- /* the search routine detected a vertex very close to the plane
+ /* the search routine detected a vertex very close to or in the midplane
the index of this vertex is stored in up
we proceed by checking the edges of this vertex */
@@ -1102,11 +947,13 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
&teststack_size);
safewhile(lw != -1) {
- i++;
+ ++i;
if (i == c->orders[up]) {
/* none of the edges of up is below the plane. Since the cell is
supposed to be convex, this means the entire cell is above or on
- the plane. This should not happen... */
+ the plane. This should not happen...
+ Furthermore, we should really NEVER end up here, as in this case
+ an error should already have be thrown above. */
voronoi_print_gnuplot_c(c);
error(
"Cell completely gone! This should not happen. (i == "
@@ -1121,21 +968,41 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
&teststack_size);
}
+ /* lp, l and lw now contain values corresponding to an edge below the
+ plane
+ rp contains the result of test_vertex for the first edge of up, for
+ reference */
+
+ /* we go on to the next edge of up, and see if we can find an edge that
+ does not extend below the plane */
+
j = i + 1;
- safewhile(j < c->orders[up]) {
+ safewhile(j < c->orders[up] && lw == -1) {
lp = voronoi_get_edge(c, up, j);
lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
&teststack_size);
- if (lw != -1) {
- break;
- }
j++;
}
+ if (lw != -1) {
+ /* the last iteration increased j by 1 too many, correct this */
+ --j;
+ }
+
+ /* j-i now contains the number of edges below the plane. We will replace
+ up by a new vertex of order this number + 2 (since 2 new edges will be
+ created inside the plane)
+ however, we do not do this if there is exactly one edge that lies in
+ the plane, and all other edges lie below, because in this case we can
+ just keep vertex up as is */
+
if (j == c->orders[up] && i == 1 && rp == 0) {
+ /* keep the order of up, and flag this event for later reference */
k = c->orders[up];
double_edge = 1;
} else {
+ /* general case: keep all edges below the plane, and create 2 new ones
+ in the plane */
k = j - i + 2;
}
@@ -1151,19 +1018,20 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
error("Too many edges!");
}
- k = 1;
-
visitflags[vindex] = 0;
+ /* the new vertex adopts the coordinates of the old vertex */
c->vertices[3 * vindex + 0] = c->vertices[3 * up + 0];
c->vertices[3 * vindex + 1] = c->vertices[3 * up + 1];
c->vertices[3 * vindex + 2] = c->vertices[3 * up + 2];
- voronoi_box_test_inside(&c->vertices[3 * vindex], c->x);
-
+ /* us contains the index of the last edge NOT below the plane
+ note that i is at least 1, so there is no need to wrap in this case */
us = i - 1;
- if (i < 0) {
- i = c->orders[up] - 1;
- }
+
+ /* copy all edges of up below the plane into the new vertex, starting from
+ edge 1 (edge 0 is reserved to connect to a newly created vertex
+ below) */
+ k = 1;
safewhile(i < j) {
qp = voronoi_get_edge(c, up, i);
qs = voronoi_get_edgeindex(c, up, i);
@@ -1172,10 +1040,13 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
voronoi_set_edgeindex(c, vindex, k, qs);
voronoi_set_edge(c, qp, qs, vindex);
voronoi_set_edgeindex(c, qp, qs, k);
+ /* disconnect up, since this vertex will be removed */
voronoi_set_edge(c, up, i, -1);
- i++;
- k++;
+ ++i;
+ ++k;
}
+
+ /* store the index of the first edge not below the plane */
if (i == c->orders[up]) {
qs = 0;
} else {
@@ -1183,6 +1054,9 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
}
} else { /* if(lw != -1) */
+ /* the first edge lies below the plane, try to find one that does not */
+
+ /* we first do a reverse search */
i = c->orders[up] - 1;
lp = voronoi_get_edge(c, up, i);
lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
@@ -1190,7 +1064,7 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
safewhile(lw == -1) {
i--;
if (i == 0) {
- /* cell unaltered */
+ /* No edge above or in the plane found: the cell is unaltered */
return;
}
lp = voronoi_get_edge(c, up, i);
@@ -1198,6 +1072,7 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
&teststack_size);
}
+ /* now we do a forward search */
j = 1;
qp = voronoi_get_edge(c, up, j);
qw = voronoi_test_vertex(&c->vertices[3 * qp], dx, r2, &q, teststack,
@@ -1209,10 +1084,20 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
&teststack_size);
}
+ /* at this point j contains the index of the first edge not below the
+ plane, i the index of the last edge not below the plane
+ we use this to compute the number of edges below the plane. up is
+ replaced by a new vertex that has that number + 2 edges (since 2 new
+ edges are created inside the plane). We again capture the special event
+ where there is only one edge not below the plane, which lies inside the
+ plane. In this case up is copied as is. */
+
if (i == j && qw == 0) {
+ /* we keep up as is, and flag this event */
double_edge = 1;
k = c->orders[up];
} else {
+ /* (c->orders[up]-1 - i) + j is the number of edges below the plane */
k = c->orders[up] - i + j + 1;
}
@@ -1227,17 +1112,23 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
if (c->offsets[vindex] + k >= VORONOI3D_MAXNUMEDGE) {
error("Too many edges!");
}
- k = 1;
visitflags[vindex] = 0;
+ /* the new vertex is just a copy of vertex up */
c->vertices[3 * vindex + 0] = c->vertices[3 * up + 0];
c->vertices[3 * vindex + 1] = c->vertices[3 * up + 1];
c->vertices[3 * vindex + 2] = c->vertices[3 * up + 2];
- voronoi_box_test_inside(&c->vertices[3 * vindex], c->x);
-
+ /* as above, us stores the index of the last edge below the plane */
us = i;
- i++;
+
+ /* copy all edges below the plane into the new vertex, starting from edge
+ 1 (edge 0 will be connected to a newly created vertex below)
+ We have to do this in two steps: first we copy the high index edges of
+ up, then the low index ones (since the edges below the plane are not a
+ continuous block of indices in this case) */
+ k = 1;
+ ++i;
safewhile(i < c->orders[up]) {
qp = voronoi_get_edge(c, up, i);
qs = voronoi_get_edgeindex(c, up, i);
@@ -1246,9 +1137,10 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
voronoi_set_edgeindex(c, vindex, k, qs);
voronoi_set_edge(c, qp, qs, vindex);
voronoi_set_edgeindex(c, qp, qs, k);
+ /* disconnect up, it will be removed */
voronoi_set_edge(c, up, i, -1);
- i++;
- k++;
+ ++i;
+ ++k;
}
i = 0;
safewhile(i < j) {
@@ -1260,28 +1152,51 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
voronoi_set_edge(c, qp, qs, vindex);
voronoi_set_edgeindex(c, qp, qs, k);
voronoi_set_edge(c, up, i, -1);
- i++;
- k++;
+ ++i;
+ ++k;
}
+ /* qs stores the index of the first edge not below the plane */
qs = j;
}
+ /* at this point, we have created a new vertex that contains all edge of up
+ below the plane, and two dangling edges: 0 and k
+ Furthermore, us stores the index of the last edge not below the plane,
+ qs the index of the first edge not below the plane */
+
+ /* now set the neighbours for the dangling edge(s) */
if (!double_edge) {
+ /* the last edge has the same neighbour as the first edge not below the
+ plane */
voronoi_set_ngb(c, vindex, k, voronoi_get_ngb(c, up, qs));
+ /* the first edge has the new neighbour as neighbour */
voronoi_set_ngb(c, vindex, 0, ngb);
} else {
+ /* up is copied as is, so we also copy its last remaining neighbour */
voronoi_set_ngb(c, vindex, 0, voronoi_get_ngb(c, up, qs));
}
+ /* add up to the delete stack */
dstack[dstack_size] = up;
++dstack_size;
+ /* make sure the variables below have the same meaning as they would have
+ if we had the non complicated setup:
+ cs contains the index of the last dangling edge of the new vertex
+ qp and q correspond to the last vertex that has been deleted
+ qs corresponds to the first edge not below the plane
+ up and us correspond to the last edge not below the plane, i.e. the edge
+ that will be the last one to connect to the new vertex
+ note that the value of i is ignored below, it is just used to temporary
+ store the new value of up */
cs = k;
qp = up;
q = u;
i = voronoi_get_edge(c, up, us);
us = voronoi_get_edgeindex(c, up, us);
up = i;
+ /* we store the index of the newly created vertex in the visitflags of the
+ last deleted vertex */
visitflags[qp] = -vindex;
} else { /* if(complicated) */
@@ -1325,10 +1240,6 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
c->vertices[3 * vindex + 2] =
c->vertices[3 * lp + 2] * r + c->vertices[3 * up + 2] * l;
-#ifdef VORONOI3D_EXPENSIVE_CHECKS
- voronoi_box_test_inside(&c->vertices[3 * vindex], c->x);
-#endif
-
/* add vertex up to the delete stack */
dstack[dstack_size] = up;
++dstack_size;
@@ -1368,36 +1279,57 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
} /* if(complicated) */
/* at this point:
- qp and up have the same value, and correspond to a vertex above the plane
- that has been deleted
- qs contains the index of the edge of qp that is next in line to be tested;
+ qp corresponds to the last vertex that has been deleted
+ up corresponds to the last vertex that should be used to connect a new
+ vertex to the newly created vertex above. In the normal case, qp and up
+ are the same vertex, but qp and up can be different if the newly created
+ vertex lies in the midplane
+ qs contains the index of the edge of qp that is next in line to be tested:
the edge that comes after the intersected edge that was deleted above
+ us corresponds to the edge of up that was connected to the vertex that is
+ now connected to the newly created vertex above
q contains the projected distance between qp and the midplane, along dx
cs contains the index of the last dangling edge of the last vertex that
was created above; we still need to connect this edge to a vertex below */
+ /* we have found one intersected edge (or at least an edge that lies inside
+ the midplane) and created one new vertex that lies in the midplane, with
+ dangling edges. We now try to find other intersected edges and create other
+ new vertices that will be connected to the first new vertex. */
+
int cp = -1;
int iqs = -1;
int new_double_edge = -1;
+ /* cp and rp both contain the index of the last vertex that was created
+ cp will be updated if we add more vertices, rp will be kept, as we need it
+ to link the last new vertex to the first new vertex in the end */
cp = vindex;
rp = vindex;
+ /* we traverse connections of the first removed vertex, until we arrive at an
+ edge that links to this vertex (or its equivalent in the degenerate
+ case) */
safewhile(qp != up || qs != us) {
+ /* test the next edge of qp */
lp = voronoi_get_edge(c, qp, qs);
lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
&teststack_size);
if (lw == 0) {
- k = 1;
+ /* degenerate case: next vertex lies inside the plane */
+
+ k = 2;
if (double_edge) {
- k = 0;
+ k = 1;
}
+ /* store the vertex and edge on the other side of the edge in qp and qs */
qs = voronoi_get_edgeindex(c, qp, qs);
qp = lp;
- iqs = qs;
- k++;
- qs++;
+ /* move on to the next edge of qp and keep the original edge for
+ reference */
+ iqs = qs;
+ ++qs;
if (qs == c->orders[qp]) {
qs = 0;
}
@@ -1405,8 +1337,8 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
&teststack_size);
safewhile(lw == -1) {
- k++;
- qs++;
+ ++k;
+ ++qs;
if (qs == c->orders[qp]) {
qs = 0;
}
@@ -1505,8 +1437,6 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
c->vertices[3 * vindex + 1] = c->vertices[3 * qp + 1];
c->vertices[3 * vindex + 2] = c->vertices[3 * qp + 2];
- voronoi_box_test_inside(&c->vertices[3 * vindex], c->x);
-
visitflags[qp] = -vindex;
dstack[dstack_size] = qp;
dstack_size++;
@@ -1558,7 +1488,15 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
double_edge = new_double_edge;
} else { /* if(lw == 0) */
+ /* normal case: next vertex lies below or above the plane */
+
if (lw == 1) {
+
+ /* vertex lies above the plane */
+
+ /* we just delete the vertex and continue with the next edge of this
+ vertex */
+
qs = voronoi_get_edgeindex(c, qp, qs) + 1;
if (qs == c->orders[lp]) {
qs = 0;
@@ -1569,6 +1507,11 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
dstack_size++;
} else {
+ /* vertex lies below the plane */
+
+ /* we have found our next intersected edge: create a new vertex and link
+ it to the other vertices */
+
if (q == l) {
error("Upper and lower vertex are the same!");
}
@@ -1624,8 +1567,10 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
c->vertices[3 * vindex + 2] =
c->vertices[3 * lp + 2] * r + c->vertices[3 * qp + 2] * l;
- voronoi_box_test_inside(&c->vertices[3 * vindex], c->x);
-
+ /* link the edges:
+ the first edge is connected to the last edge of the previous new
+ vertex. The last edge will be connected to the next new vertex, and
+ is left open for the moment */
ls = voronoi_get_edgeindex(c, qp, qs);
voronoi_set_edge(c, vindex, 0, cp);
voronoi_set_edge(c, vindex, 1, lp);
@@ -1641,10 +1586,14 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
voronoi_set_ngb(c, vindex, 1, voronoi_get_ngb(c, qp, qs));
voronoi_set_ngb(c, vindex, 2, voronoi_get_ngb(c, lp, ls));
- qs++;
+ /* continue with the next edge of qp (the last vertex above the
+ midplane */
+ ++qs;
if (qs == c->orders[qp]) {
qs = 0;
}
+ /* store the last newly created vertex and its dangling edge for the
+ next iteration */
cp = vindex;
cs = 2;
} /* if(lw == 1) */
@@ -1653,11 +1602,21 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
} /* while() */
+ /* we finished adding new vertices. Now connect the last dangling edge of the
+ last newly created vertex to the first dangling edge of the first newly
+ created vertex */
voronoi_set_edge(c, cp, cs, rp);
voronoi_set_edge(c, rp, 0, cp);
voronoi_set_edgeindex(c, cp, cs, 0);
voronoi_set_edgeindex(c, rp, 0, cs);
+ /* now remove the vertices in the delete stack */
+
+ /* the algorithm above did not necessarily visit all vertices above the plane.
+ here we scan for vertices that are linked to vertices that are to be
+ removed and add them to the delete stack if necessary
+ this only works because we made sure that all deleted vertices no longer
+ have edges that connect them to vertices that need to stay */
for (i = 0; i < dstack_size; i++) {
for (j = 0; j < c->orders[dstack[i]]; j++) {
if (voronoi_get_edge(c, dstack[i], j) >= 0) {
@@ -1902,11 +1861,6 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
/* Update the cell values. */
voronoi3d_cell_copy(&new_cell, c);
- /* Final checks */
- for (i = 0; i < c->nvert; ++i) {
- voronoi_box_test_inside(&c->vertices[3 * i], c->x);
- }
-
#ifdef VORONOI3D_EXPENSIVE_CHECKS
voronoi_check_cell_consistency(c);
#endif
@@ -1980,10 +1934,6 @@ __attribute__((always_inline)) INLINE void voronoi_calculate_cell(
float tcentroid[3];
float tvol;
- for (i = 0; i < cell->nvert; ++i) {
- voronoi_box_test_inside(&cell->vertices[3 * i], cell->x);
- }
-
/* we need to calculate the volume of the tetrahedra formed by the first
vertex and the triangles that make up the other faces
since we do not store faces explicitly, this means keeping track of the
@@ -2186,15 +2136,18 @@ __attribute__((always_inline)) INLINE void voronoi_calculate_faces(
*
* @param cell 3D Voronoi cell to initialize.
* @param x Position of the generator of the cell.
+ * @param anchor Anchor of the simulation box.
+ * @param side Side lengths of the simulation box.
*/
__attribute__((always_inline)) INLINE void voronoi_cell_init(
- struct voronoi_cell *cell, const double *x) {
+ struct voronoi_cell *cell, const double *x, const double *anchor,
+ const double *side) {
cell->x[0] = x[0];
cell->x[1] = x[1];
cell->x[2] = x[2];
- voronoi_initialize(cell);
+ voronoi_initialize(cell, anchor, side);
cell->volume = 0.0f;
cell->centroid[0] = 0.0f;
diff --git a/src/hydro/Shadowswift/voronoi3d_cell.h b/src/hydro/Shadowswift/voronoi3d_cell.h
index db7bf95b5abe96419062296d55714d86ccc11160..ef43eff1745f48219af14aec2455aaa5e5b0d47a 100644
--- a/src/hydro/Shadowswift/voronoi3d_cell.h
+++ b/src/hydro/Shadowswift/voronoi3d_cell.h
@@ -48,23 +48,29 @@ struct voronoi_cell {
float vertices[3 * VORONOI3D_MAXNUMVERT];
/* Number of edges for every vertex. */
- int orders[VORONOI3D_MAXNUMVERT];
+ char orders[VORONOI3D_MAXNUMVERT];
/* Offsets of the edges, edgeindices and neighbours corresponding to a
particular vertex in the internal arrays */
int offsets[VORONOI3D_MAXNUMVERT];
- /* Edge information. */
+ /* Edge information. Edges are ordered counterclockwise w.r.t. a vector
+ pointing from the cell generator to the vertex. */
int edges[VORONOI3D_MAXNUMEDGE];
/* Additional edge information. */
- int edgeindices[VORONOI3D_MAXNUMEDGE];
-
- /* Neighbour information. */
+ char edgeindices[VORONOI3D_MAXNUMEDGE];
+
+ /* Neighbour information. This field is used differently depending on where we
+ are in the algorithm. During cell construction, it contains, for every edge
+ of every vertex, the index of the neighbour that generates the face
+ counterclockwise of the edge w.r.t. a vector pointing from the vertex along
+ the edge. After cell finalization, it contains a neighbour for every face,
+ in the same order as the face_areas and face_midpoints arrays. */
unsigned long long ngbs[VORONOI3D_MAXNUMEDGE];
/* Number of faces of the cell. */
- int nface;
+ unsigned char nface;
/* Surface areas of the cell faces. */
float face_areas[VORONOI3D_MAXFACE];
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 46785b4b2d5b958f6db3bd9813d139575217d6fe..818c1b6349192ed73b28cd4c3ae771f89a3754cd 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -44,6 +44,12 @@ void hydro_props_init(struct hydro_props *p,
p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
+#ifdef SHADOWFAX_SPH
+ /* change the meaning of target_neighbours and delta_neighbours */
+ p->target_neighbours = 1.0f;
+ p->delta_neighbours = 0.0f;
+#endif
+
/* Maximal smoothing length */
p->h_max = parser_get_opt_param_float(params, "SPH:h_max",
hydro_props_default_h_max);
diff --git a/src/runner.c b/src/runner.c
index 792f03f9b454d2afe0116f9bda194bb00ce79543..8240f77763d73d3ebc1e7a402194054b5c9ea7d5 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -489,6 +489,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
const int count = c->count;
const int gcount = c->gcount;
const struct engine *e = r->e;
+ const struct space *s = e->s;
TIMER_TIC;
@@ -514,7 +515,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
if (part_is_active(p, e)) {
/* Get ready for a density calculation */
- hydro_init_part(p);
+ hydro_init_part(p, &s->hs);
}
}
@@ -596,6 +597,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
const struct engine *e = r->e;
+ const struct space *s = e->s;
const float hydro_h_max = e->hydro_properties->h_max;
const float target_wcount = e->hydro_properties->target_neighbours;
const float max_wcount =
@@ -676,7 +678,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
redo += 1;
/* Re-initialise everything */
- hydro_init_part(p);
+ hydro_init_part(p, &s->hs);
/* Off we go ! */
continue;
diff --git a/src/space.c b/src/space.c
index 00e5b87ad6d3587b02d14d514c9788f8639c5dd5..5899cc396b84cd3aa92aab10a687caaa8d6d81e9 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2692,6 +2692,8 @@ void space_init(struct space *s, const struct swift_params *params,
bzero(s->xparts, Npart * sizeof(struct xpart));
}
+ hydro_space_init(&s->hs, s);
+
/* Set the particles in a state where they are ready for a run */
space_init_parts(s);
space_init_xparts(s);
diff --git a/src/space.h b/src/space.h
index e3886364bfa7b2f85e724993e5fe3c2d30663bd2..55ebdad91db4de7f7baa422b2e00b7dde438d8ca 100644
--- a/src/space.h
+++ b/src/space.h
@@ -31,6 +31,7 @@
/* Includes. */
#include "cell.h"
+#include "hydro_space.h"
#include "lock.h"
#include "parser.h"
#include "part.h"
@@ -69,6 +70,9 @@ struct space {
/*! Is the space periodic? */
int periodic;
+ /*! Extra space information needed for some hydro schemes. */
+ struct hydro_space hs;
+
/*! Are we doing gravity? */
int gravity;
diff --git a/src/tools.c b/src/tools.c
index 89ac286fb435c01b361bdea66e62dd2d7f41ee24..73684c82662870d368f7dd360c84635654f06434 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -144,7 +144,7 @@ void pairs_single_density(double *dim, long long int pid,
p = parts[k];
printf("pairs_single: part[%i].id == %lli.\n", k, pid);
- hydro_init_part(&p);
+ hydro_init_part(&p, NULL);
/* Loop over all particle pairs. */
for (k = 0; k < N; k++) {
@@ -459,7 +459,7 @@ void engine_single_density(double *dim, long long int pid,
p = parts[k];
/* Clear accumulators. */
- hydro_init_part(&p);
+ hydro_init_part(&p, NULL);
/* Loop over all particle pairs (force). */
for (k = 0; k < N; k++) {
diff --git a/tests/benchmarkInteractions.c b/tests/benchmarkInteractions.c
index 14648b0693714d668d35c0795503cc3ccfe9754d..0f5b3d2eb294c13e3035885b511c702e6f0cd540 100644
--- a/tests/benchmarkInteractions.c
+++ b/tests/benchmarkInteractions.c
@@ -450,10 +450,6 @@ void test_interactions(struct part test_part, struct part *parts, size_t count,
#endif
}
-#if defined(SHADOWFAX_SPH)
-VORONOI_DECLARE_GLOBAL_VARIABLES()
-#endif
-
/* And go... */
int main(int argc, char *argv[]) {
size_t runs = 10000;
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 3ff812593abcbb1c7d9779b571ccf300cc184fe6..21fc3f5407f90d9b75485d08bf1077e0a20e88b2 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -422,10 +422,6 @@ void runner_doself1_density(struct runner *r, struct cell *ci);
void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself2_force(struct runner *r, struct cell *ci);
-#if defined(SHADOWFAX_SPH)
-VORONOI_DECLARE_GLOBAL_VARIABLES()
-#endif
-
/* And go... */
int main(int argc, char *argv[]) {
diff --git a/tests/test27cells.c b/tests/test27cells.c
index e5e73ad1e503b71c25bbe95cf5ede9010977c6b6..2390ae8aeac465c49831984c0d24817b50d757b3 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -119,7 +119,13 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
part->conserved.mass = density * volume / count;
- voronoi_cell_init(&part->cell, part->x);
+
+#ifdef SHADOWFAX_SPH
+ double anchor[3] = {0., 0., 0.};
+ double side[3] = {1., 1., 1.};
+ voronoi_cell_init(&part->cell, part->x, anchor, side);
+#endif
+
#else
part->mass = density * volume / count;
#endif
@@ -177,7 +183,7 @@ void clean_up(struct cell *ci) {
*/
void zero_particle_fields(struct cell *c) {
for (int pid = 0; pid < c->count; pid++) {
- hydro_init_part(&c->parts[pid]);
+ hydro_init_part(&c->parts[pid], NULL);
}
}
@@ -299,10 +305,6 @@ void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself1_density(struct runner *r, struct cell *ci);
void runner_doself1_density_vec(struct runner *r, struct cell *ci);
-#if defined(SHADOWFAX_SPH)
-VORONOI_DECLARE_GLOBAL_VARIABLES()
-#endif
-
/* And go... */
int main(int argc, char *argv[]) {
diff --git a/tests/testMatrixInversion.c b/tests/testMatrixInversion.c
index d28bb49b03ef0890807479197def58590c625f2d..9a45cd52d6f5d3ec96cc6d3f34fd683971f4cf19 100644
--- a/tests/testMatrixInversion.c
+++ b/tests/testMatrixInversion.c
@@ -22,10 +22,7 @@
#include "const.h"
#include "dimension.h"
#include "error.h"
-
-double random_uniform(double a, double b) {
- return (rand() / (double)RAND_MAX) * (b - a) + a;
-}
+#include "tools.h"
void setup_matrix(float A[3][3]) {
A[0][0] = random_uniform(-1.0, 1.0);
diff --git a/tests/testPair.c b/tests/testPair.c
index 3719fb195581fade6b49b9c7a07e1060637f371f..c2533b63b902e3bdc7e7cae6fcbcf50c87dee4af 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -116,7 +116,7 @@ void clean_up(struct cell *ci) {
*/
void zero_particle_fields(struct cell *c) {
for (int pid = 0; pid < c->count; pid++) {
- hydro_init_part(&c->parts[pid]);
+ hydro_init_part(&c->parts[pid], NULL);
}
}
@@ -187,10 +187,6 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
/* Just a forward declaration... */
void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
-#if defined(SHADOWFAX_SPH)
-VORONOI_DECLARE_GLOBAL_VARIABLES()
-#endif
-
int main(int argc, char *argv[]) {
size_t particles = 0, runs = 0, volume, type = 0;
double offset[3] = {0, 0, 0}, h = 1.1255, size = 1., rho = 1.;
diff --git a/tests/testRiemannExact.c b/tests/testRiemannExact.c
index 2a6b1ab80d889dae6f7008cc4044be0f3be43478..82b12449f1b199133de5a74fe7b68b5c386c9cf5 100644
--- a/tests/testRiemannExact.c
+++ b/tests/testRiemannExact.c
@@ -20,10 +20,7 @@
#include <string.h>
#include "error.h"
#include "riemann/riemann_exact.h"
-
-double random_uniform(double a, double b) {
- return (rand() / (double)RAND_MAX) * (b - a) + a;
-}
+#include "tools.h"
int opposite(float a, float b) {
if ((a - b)) {
diff --git a/tests/testRiemannHLLC.c b/tests/testRiemannHLLC.c
index 59ccdb40a76c973305d5009cfe61c47af558f9fc..6bdf1192a6da8482d562895027d761f73ecc71de 100644
--- a/tests/testRiemannHLLC.c
+++ b/tests/testRiemannHLLC.c
@@ -20,10 +20,7 @@
#include <string.h>
#include "error.h"
#include "riemann/riemann_hllc.h"
-
-double random_uniform(double a, double b) {
- return (rand() / (double)RAND_MAX) * (b - a) + a;
-}
+#include "tools.h"
int consistent_with_zero(float val) { return fabs(val) < 1.e-4; }
diff --git a/tests/testRiemannTRRS.c b/tests/testRiemannTRRS.c
index fbdd723e850e87f497c549531706313da476b801..4a0eac0be23581e175d2c0e599b786fd4508b14a 100644
--- a/tests/testRiemannTRRS.c
+++ b/tests/testRiemannTRRS.c
@@ -20,10 +20,7 @@
#include <string.h>
#include "error.h"
#include "riemann/riemann_trrs.h"
-
-double random_uniform(double a, double b) {
- return (rand() / (double)RAND_MAX) * (b - a) + a;
-}
+#include "tools.h"
int opposite(float a, float b) {
if ((a - b)) {
diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c
index cbe857f3881ee2f0542afbf0cfd62f03185f3928..73c5708a6add174b88f26cc716a716fa2ad81709 100644
--- a/tests/testSymmetry.c
+++ b/tests/testSymmetry.c
@@ -26,10 +26,6 @@
#include "swift.h"
-#if defined(SHADOWFAX_SPH)
-VORONOI_DECLARE_GLOBAL_VARIABLES()
-#endif
-
int main(int argc, char *argv[]) {
/* Choke if need be */
@@ -37,9 +33,9 @@ int main(int argc, char *argv[]) {
#if defined(SHADOWFAX_SPH)
/* Initialize the Voronoi simulation box */
- float box_anchor[3] = {-2.0f, -2.0f, -2.0f};
- float box_side[3] = {6.0f, 6.0f, 6.0f};
- voronoi_set_box(box_anchor, box_side);
+ double box_anchor[3] = {-2.0f, -2.0f, -2.0f};
+ double box_side[3] = {6.0f, 6.0f, 6.0f};
+/* voronoi_set_box(box_anchor, box_side);*/
#endif
/* Create two random particles (don't do this at home !) */
@@ -106,8 +102,8 @@ int main(int argc, char *argv[]) {
pj.force.dt = 0.001;
#ifdef SHADOWFAX_SPH
- voronoi_cell_init(&pi.cell, pi.x);
- voronoi_cell_init(&pj.cell, pj.x);
+ voronoi_cell_init(&pi.cell, pi.x, box_anchor, box_side);
+ voronoi_cell_init(&pj.cell, pj.x, box_anchor, box_side);
#endif
#endif
diff --git a/tests/testTimeIntegration.c b/tests/testTimeIntegration.c
index 185db6edde06ec82d0a0f87bbb166f05bd89e8b0..42a3d224f43d580e512119edc55051bd22719a3b 100644
--- a/tests/testTimeIntegration.c
+++ b/tests/testTimeIntegration.c
@@ -22,10 +22,6 @@
#include <stdlib.h>
#include <string.h>
-#if defined(SHADOWFAX_SPH)
-VORONOI_DECLARE_GLOBAL_VARIABLES()
-#endif
-
/**
* @brief Test the kick-drift-kick leapfrog integration
* via a Sun-Earth simulation
diff --git a/tests/testVoronoi1D.c b/tests/testVoronoi1D.c
index dfc0d77b35d83976cfeeb4aed68f9c69d5ed395d..d16a36d9449d7bfdb2c74408efad61b219b1d7e3 100644
--- a/tests/testVoronoi1D.c
+++ b/tests/testVoronoi1D.c
@@ -21,10 +21,13 @@
int main() {
+ double box_anchor[1] = {-0.5};
+ double box_side[1] = {2.};
+
/* Create a Voronoi cell */
double x[1] = {0.5f};
struct voronoi_cell cell;
- voronoi_cell_init(&cell, x);
+ voronoi_cell_init(&cell, x, box_anchor, box_side);
/* Interact with a left and right neighbour */
float xL[1] = {0.5f};
diff --git a/tests/testVoronoi2D.c b/tests/testVoronoi2D.c
index fa6c1500928095c1269c989b87e1195d7866b5b0..8e2c875ee1b41828a0b39aa896f02dff0ccbac88 100644
--- a/tests/testVoronoi2D.c
+++ b/tests/testVoronoi2D.c
@@ -18,19 +18,16 @@
******************************************************************************/
#include "hydro/Shadowswift/voronoi2d_algorithm.h"
+#include "tools.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);
+ double anchor[3] = {-0.5f, -0.5f, -0.5f};
+ double side[3] = {2.0f, 2.0f, 2.0f};
/* test initialization and finalization algorithms */
{
@@ -39,7 +36,7 @@ int main() {
struct voronoi_cell cell;
double x[3] = {0.5, 0.5, 0.5};
- voronoi_cell_init(&cell, x);
+ voronoi_cell_init(&cell, x, anchor, side);
float maxradius = voronoi_cell_finalize(&cell);
@@ -67,9 +64,9 @@ int main() {
struct voronoi_cell *cell_i, *cell_j;
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);
+ x[0] = random_uniform(0., 1.);
+ x[1] = random_uniform(0., 1.);
+ voronoi_cell_init(&cells[i], x, anchor, side);
#ifdef VORONOI_VERBOSE
message("cell[%i]: %g %g", i, x[0], x[1]);
#endif
@@ -149,7 +146,7 @@ int main() {
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);
+ voronoi_cell_init(&cells[10 * i + j], x, anchor, side);
}
}
diff --git a/tests/testVoronoi3D.c b/tests/testVoronoi3D.c
index c4a688758b49e07386189f0e8494b9cd787a4cbc..b4f219a41368bb3ce4e8111ae44c43e7fa1f7441 100644
--- a/tests/testVoronoi3D.c
+++ b/tests/testVoronoi3D.c
@@ -21,6 +21,7 @@
#include "error.h"
#include "hydro/Shadowswift/voronoi3d_algorithm.h"
#include "part.h"
+#include "tools.h"
/* Number of random generators to use in the first grid build test */
#define TESTVORONOI3D_NUMCELL_RANDOM 100
@@ -36,6 +37,95 @@
(TESTVORONOI3D_NUMCELL_CARTESIAN_1D * TESTVORONOI3D_NUMCELL_CARTESIAN_1D * \
TESTVORONOI3D_NUMCELL_CARTESIAN_1D)
+/* Bottom front left corner and side lengths of the large box that contains all
+ particles and is used as initial cell at the start of the construction */
+#define VORONOI3D_BOX_ANCHOR_X 0.0f
+#define VORONOI3D_BOX_ANCHOR_Y 0.0f
+#define VORONOI3D_BOX_ANCHOR_Z 0.0f
+#define VORONOI3D_BOX_SIDE_X 1.0f
+#define VORONOI3D_BOX_SIDE_Y 1.0f
+#define VORONOI3D_BOX_SIDE_Z 1.0f
+
+/**
+ * @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.
+ */
+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.
+ */
+void voronoi_get_box_centroid(float *box_centroid) {
+ box_centroid[0] = 0.5f * VORONOI3D_BOX_SIDE_X + VORONOI3D_BOX_ANCHOR_X;
+ box_centroid[1] = 0.5f * VORONOI3D_BOX_SIDE_Y + VORONOI3D_BOX_ANCHOR_Y;
+ box_centroid[2] = 0.5f * VORONOI3D_BOX_SIDE_Z + VORONOI3D_BOX_ANCHOR_Z;
+}
+
+/**
+ * @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.
+ */
+float voronoi_get_box_face(unsigned long long id, float *face_midpoint) {
+
+ if (id == VORONOI3D_BOX_FRONT) {
+ face_midpoint[0] = 0.5f * VORONOI3D_BOX_SIDE_X + VORONOI3D_BOX_ANCHOR_X;
+ face_midpoint[1] = VORONOI3D_BOX_ANCHOR_Y;
+ face_midpoint[2] = 0.5f * VORONOI3D_BOX_SIDE_Z + VORONOI3D_BOX_ANCHOR_Z;
+ return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Z;
+ }
+ if (id == VORONOI3D_BOX_BACK) {
+ face_midpoint[0] = 0.5f * VORONOI3D_BOX_SIDE_X + VORONOI3D_BOX_ANCHOR_X;
+ face_midpoint[1] = VORONOI3D_BOX_ANCHOR_Y + VORONOI3D_BOX_SIDE_Y;
+ face_midpoint[2] = 0.5f * VORONOI3D_BOX_SIDE_Z + VORONOI3D_BOX_ANCHOR_Z;
+ return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Z;
+ }
+
+ if (id == VORONOI3D_BOX_BOTTOM) {
+ face_midpoint[0] = 0.5f * VORONOI3D_BOX_SIDE_X + VORONOI3D_BOX_ANCHOR_X;
+ face_midpoint[1] = 0.5f * VORONOI3D_BOX_SIDE_Y + VORONOI3D_BOX_ANCHOR_Y;
+ face_midpoint[2] = VORONOI3D_BOX_ANCHOR_Z;
+ return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Y;
+ }
+ if (id == VORONOI3D_BOX_TOP) {
+ face_midpoint[0] = 0.5f * VORONOI3D_BOX_SIDE_X + VORONOI3D_BOX_ANCHOR_X;
+ face_midpoint[1] = 0.5f * VORONOI3D_BOX_SIDE_Y + VORONOI3D_BOX_ANCHOR_Y;
+ face_midpoint[2] = VORONOI3D_BOX_ANCHOR_Z + VORONOI3D_BOX_SIDE_Z;
+ return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Y;
+ }
+
+ if (id == VORONOI3D_BOX_LEFT) {
+ face_midpoint[0] = VORONOI3D_BOX_ANCHOR_X;
+ face_midpoint[1] = 0.5f * VORONOI3D_BOX_SIDE_Y + VORONOI3D_BOX_ANCHOR_Y;
+ face_midpoint[2] = 0.5f * VORONOI3D_BOX_SIDE_Z + VORONOI3D_BOX_ANCHOR_Z;
+ return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Y;
+ }
+ if (id == VORONOI3D_BOX_RIGHT) {
+ face_midpoint[0] = VORONOI3D_BOX_ANCHOR_X + VORONOI3D_BOX_SIDE_X;
+ face_midpoint[1] = 0.5f * VORONOI3D_BOX_SIDE_Y + VORONOI3D_BOX_ANCHOR_Y;
+ face_midpoint[2] = 0.5f * VORONOI3D_BOX_SIDE_Z + VORONOI3D_BOX_ANCHOR_Z;
+ return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Y;
+ }
+
+ return 0.0f;
+}
+
/**
* @brief Check if voronoi_volume_tetrahedron() works
*/
@@ -69,6 +159,12 @@ void test_voronoi_centroid_tetrahedron() {
* @brief Check if voronoi_calculate_cell() works
*/
void test_calculate_cell() {
+
+ double box_anchor[3] = {VORONOI3D_BOX_ANCHOR_X, VORONOI3D_BOX_ANCHOR_Y,
+ VORONOI3D_BOX_ANCHOR_Z};
+ double box_side[3] = {VORONOI3D_BOX_SIDE_X, VORONOI3D_BOX_SIDE_Y,
+ VORONOI3D_BOX_SIDE_Z};
+
struct voronoi_cell cell;
cell.x[0] = 0.5f;
@@ -76,7 +172,7 @@ void test_calculate_cell() {
cell.x[2] = 0.5f;
/* Initialize the cell to a large cube. */
- voronoi_initialize(&cell);
+ voronoi_initialize(&cell, box_anchor, box_side);
/* Calculate the volume and centroid of the large cube. */
voronoi_calculate_cell(&cell);
/* Calculate the faces. */
@@ -1130,11 +1226,17 @@ void test_paths() {
#ifdef SHADOWFAX_SPH
void set_coordinates(struct part *p, double x, double y, double z,
unsigned int id) {
+
+ double box_anchor[3] = {VORONOI3D_BOX_ANCHOR_X, VORONOI3D_BOX_ANCHOR_Y,
+ VORONOI3D_BOX_ANCHOR_Z};
+ double box_side[3] = {VORONOI3D_BOX_SIDE_X, VORONOI3D_BOX_SIDE_Y,
+ VORONOI3D_BOX_SIDE_Z};
+
p->x[0] = x;
p->x[1] = y;
p->x[2] = z;
p->id = id;
- voronoi_cell_init(&p->cell, p->x);
+ voronoi_cell_init(&p->cell, p->x, box_anchor, box_side);
}
#endif
@@ -1206,14 +1308,13 @@ void test_degeneracies() {
#endif
}
-VORONOI_DECLARE_GLOBAL_VARIABLES()
-
int main() {
/* Set the all enclosing simulation box dimensions */
- float box_anchor[3] = {0.0f, 0.0f, 0.0f};
- float box_side[3] = {1.0f, 1.0f, 1.0f};
- voronoi_set_box(box_anchor, box_side);
+ double box_anchor[3] = {VORONOI3D_BOX_ANCHOR_X, VORONOI3D_BOX_ANCHOR_Y,
+ VORONOI3D_BOX_ANCHOR_Z};
+ double box_side[3] = {VORONOI3D_BOX_SIDE_X, VORONOI3D_BOX_SIDE_Y,
+ VORONOI3D_BOX_SIDE_Z};
/* Check basic Voronoi cell functions */
test_voronoi_volume_tetrahedron();
@@ -1228,7 +1329,7 @@ int main() {
/* Create a Voronoi cell */
double x[3] = {0.5f, 0.5f, 0.5f};
struct voronoi_cell cell;
- voronoi_cell_init(&cell, x);
+ voronoi_cell_init(&cell, x, box_anchor, box_side);
/* Interact with neighbours */
float x0[3] = {0.5f, 0.0f, 0.0f};
@@ -1326,10 +1427,10 @@ int main() {
/* initialize cells with random generator locations */
for (i = 0; i < TESTVORONOI3D_NUMCELL_RANDOM; ++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);
+ x[0] = random_uniform(0., 1.);
+ x[1] = random_uniform(0., 1.);
+ x[2] = random_uniform(0., 1.);
+ voronoi_cell_init(&cells[i], x, box_anchor, box_side);
}
/* interact the cells */
@@ -1370,7 +1471,7 @@ int main() {
struct voronoi_cell cells[TESTVORONOI3D_NUMCELL_CARTESIAN_3D];
struct voronoi_cell *cell_i, *cell_j;
- /* initialize cells with random generator locations */
+ /* initialize cells with Cartesian generator locations */
for (i = 0; i < TESTVORONOI3D_NUMCELL_CARTESIAN_1D; ++i) {
for (j = 0; j < TESTVORONOI3D_NUMCELL_CARTESIAN_1D; ++j) {
for (k = 0; k < TESTVORONOI3D_NUMCELL_CARTESIAN_1D; ++k) {
@@ -1380,7 +1481,7 @@ int main() {
voronoi_cell_init(&cells[TESTVORONOI3D_NUMCELL_CARTESIAN_1D *
TESTVORONOI3D_NUMCELL_CARTESIAN_1D * i +
TESTVORONOI3D_NUMCELL_CARTESIAN_1D * j + k],
- x);
+ x, box_anchor, box_side);
}
}
}