From a66c747990cc1155050d93c2083d7fa24141d248 Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@gmail.com>
Date: Fri, 17 Mar 2017 17:19:19 +0000
Subject: [PATCH] Replaced global variables by a new hydro_space struct that is
part of space. This struct is then passed on to hydro_init_part, and used to
initialize the Voronoi cells. For non SHADOWFAX_SPH schemes, this struct is
empty. Replaced random_uniform declarations in some unit tests with the
tools.h version. Added some more comments to 3D Voronoi code; it is now
almost completely documented (meaning I almost completely understand it now).
---
examples/main.c | 32 --
src/Makefile.am | 6 +-
src/hydro/Default/hydro.h | 6 +-
src/hydro/Gadget2/hydro.h | 6 +-
src/hydro/Gizmo/hydro.h | 4 +-
src/hydro/Minimal/hydro.h | 6 +-
src/hydro/PressureEntropy/hydro.h | 6 +-
src/hydro/Shadowswift/hydro.h | 6 +-
src/hydro/Shadowswift/voronoi1d_algorithm.h | 12 +-
src/hydro/Shadowswift/voronoi2d_algorithm.h | 59 +--
src/hydro/Shadowswift/voronoi3d_algorithm.h | 461 +++++++++-----------
src/hydro/Shadowswift/voronoi3d_cell.h | 18 +-
src/hydro_properties.c | 6 +
src/runner.c | 6 +-
src/space.c | 2 +
src/space.h | 4 +
src/tools.c | 4 +-
tests/benchmarkInteractions.c | 4 -
tests/test125cells.c | 4 -
tests/test27cells.c | 14 +-
tests/testMatrixInversion.c | 5 +-
tests/testPair.c | 6 +-
tests/testRiemannExact.c | 5 +-
tests/testRiemannHLLC.c | 5 +-
tests/testRiemannTRRS.c | 5 +-
tests/testSymmetry.c | 14 +-
tests/testTimeIntegration.c | 4 -
tests/testVoronoi1D.c | 5 +-
tests/testVoronoi2D.c | 19 +-
tests/testVoronoi3D.c | 129 +++++-
30 files changed, 427 insertions(+), 436 deletions(-)
diff --git a/examples/main.c b/examples/main.c
index 41fbd8aba8..8d4289c872 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 29c1102184..a29cf5521f 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 a614d08c30..051c22f46b 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 cc7b422ccb..747c81a8e6 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 643489912a..60ff8ccee0 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 56078a8256..8f216a550a 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 20238896f1..4c4868cd37 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 82db9f4b18..0568d47ee7 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 31045895be..74cc5f1dbf 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 a799057ed6..ad4d565fbf 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 318c83bbd5..961f6e91ff 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 db7bf95b5a..ef43eff174 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 46785b4b2d..818c1b6349 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 792f03f9b4..8240f77763 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 00e5b87ad6..5899cc396b 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 e3886364bf..55ebdad91d 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 89ac286fb4..73684c8266 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 14648b0693..0f5b3d2eb2 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 3ff812593a..21fc3f5407 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 e5e73ad1e5..2390ae8aea 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 d28bb49b03..9a45cd52d6 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 3719fb1955..c2533b63b9 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 2a6b1ab80d..82b12449f1 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 59ccdb40a7..6bdf1192a6 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 fbdd723e85..4a0eac0be2 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 cbe857f388..73c5708a6a 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 185db6edde..42a3d224f4 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 dfc0d77b35..d16a36d944 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 fa6c150092..8e2c875ee1 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 c4a688758b..b4f219a413 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);
}
}
}
--
GitLab