diff --git a/examples/main.c b/examples/main.c index 39b2e3a7578a306cee39b3f443a7060261eb7772..41fbd8aba892cc1677c3317fa3579ca6ee5dd7fa 100644 --- a/examples/main.c +++ b/examples/main.c @@ -101,8 +101,8 @@ void print_help_message() { "parameter file.\n"); } -#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D) -VORONOI3D_DECLARE_GLOBAL_VARIABLES() +#if defined(SHADOWFAX_SPH) +VORONOI_DECLARE_GLOBAL_VARIABLES() #endif /** @@ -439,7 +439,7 @@ int main(int argc, char *argv[]) { fflush(stdout); } -#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D) +#if defined(SHADOWFAX_SPH) /* set the *global* box dimensions */ float box_anchor[3], box_side[3]; if (periodic) { diff --git a/src/hydro/Shadowswift/voronoi1d_algorithm.h b/src/hydro/Shadowswift/voronoi1d_algorithm.h index 1bb6111337a8a0771fd31616870a5131729cfaf7..5df2a1fe4b6c3566f16734b524277facddfc49cb 100644 --- a/src/hydro/Shadowswift/voronoi1d_algorithm.h +++ b/src/hydro/Shadowswift/voronoi1d_algorithm.h @@ -27,6 +27,12 @@ #include "inline.h" #include "voronoi1d_cell.h" +#define VORONOI_DECLARE_GLOBAL_VARIABLES() + +__attribute__((always_inline)) INLINE static void voronoi_set_box(float *anchor, + float *side) { +} + /** * @brief Initialize a 1D Voronoi cell * diff --git a/src/hydro/Shadowswift/voronoi2d_algorithm.h b/src/hydro/Shadowswift/voronoi2d_algorithm.h index 437c786f6c002ee85dd66eafb9e9056a8e20cf18..a69d317fa98f83d468b5071812de7d00b152eaab 100644 --- a/src/hydro/Shadowswift/voronoi2d_algorithm.h +++ b/src/hydro/Shadowswift/voronoi2d_algorithm.h @@ -25,8 +25,36 @@ #include <stdlib.h> #include "error.h" #include "inline.h" +#include "minmax.h" #include "voronoi2d_cell.h" +#define VORONOI2D_BOX_LEFT 18446744073709551602llu +#define VORONOI2D_BOX_RIGHT 18446744073709551603llu +#define VORONOI2D_BOX_TOP 18446744073709551604llu +#define VORONOI2D_BOX_BOTTOM 18446744073709551605llu + +extern float global_voronoi_box_anchor[2]; +extern float global_voronoi_box_side[2]; + +#define VORONOI_DECLARE_GLOBAL_VARIABLES() \ + float global_voronoi_box_anchor[2]; \ + float global_voronoi_box_side[2]; + +#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] + +__attribute__((always_inline)) INLINE static void voronoi_set_box(float *anchor, + 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 * @@ -34,7 +62,38 @@ * @param x Position of the generator of the cell. */ __attribute__((always_inline)) INLINE void voronoi_cell_init( - struct voronoi_cell *cell, double *x) {} + struct voronoi_cell *cell, double *x) { + + cell->x[0] = x[0]; + cell->x[1] = x[1]; + + 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[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[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[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->ngbs[0] = VORONOI2D_BOX_LEFT; + cell->ngbs[1] = VORONOI2D_BOX_TOP; + cell->ngbs[2] = VORONOI2D_BOX_RIGHT; + cell->ngbs[3] = VORONOI2D_BOX_BOTTOM; + + cell->volume = 0.0f; + cell->centroid[0] = 0.0f; + cell->centroid[1] = 0.0f; +} /** * @brief Interact a 2D Voronoi cell with a particle with given relative @@ -57,7 +116,51 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact( __attribute__((always_inline)) INLINE float voronoi_cell_finalize( struct voronoi_cell *cell) { - return 1.0f; + int i; + float vertices[VORONOI2D_MAXNUMVERT][2]; + float A, x[2], y[2], r2, r2max; + + /* make a copy of the vertices (they are overwritten when the face midpoints + are computed */ + for (i = 0; i < cell->nvert; ++i) { + vertices[i][0] = cell->vertices[i][0]; + vertices[i][1] = cell->vertices[i][1]; + } + + r2max = 0.0f; + for (i = 0; i < cell->nvert; ++i) { + if (i < cell->nvert - 1) { + x[0] = vertices[i][0]; + y[0] = vertices[i][1]; + x[1] = vertices[i + 1][0]; + y[1] = vertices[i + 1][1]; + } else { + x[0] = vertices[i][0]; + y[0] = vertices[i][1]; + x[1] = vertices[0][0]; + y[1] = vertices[0][1]; + } + A = x[1] * y[0] - x[0] * y[1]; + cell->volume += A; + cell->centroid[0] += (x[0] + x[1]) * A; + cell->centroid[1] += (y[0] + y[1]) * A; + + cell->face_midpoints[i][0] = 0.5f * (x[0] + x[1]) + cell->x[0]; + cell->face_midpoints[i][1] = 0.5f * (y[0] + y[1]) + cell->x[1]; + + r2 = x[0] * x[0] + y[0] * y[0]; + r2max = max(r2max, r2); + } + + cell->volume *= 0.5f; + A = 6 * cell->volume; + cell->centroid[0] /= A; + cell->centroid[1] /= A; + + cell->centroid[0] += cell->x[0]; + cell->centroid[1] += cell->x[1]; + + return 2.0f * sqrtf(r2max); } /** @@ -72,7 +175,21 @@ __attribute__((always_inline)) INLINE float voronoi_cell_finalize( __attribute__((always_inline)) INLINE float voronoi_get_face( struct voronoi_cell *cell, unsigned long long ngb, float *midpoint) { - return 0.0f; + /* look up the neighbour */ + int i = 0; + while (i < cell->nvert && cell->ngbs[i] != ngb) { + ++i; + } + + if (i == cell->nvert) { + /* The given cell is not a neighbour. */ + return 0.0f; + } + + midpoint[0] = cell->face_midpoints[i][0]; + midpoint[1] = cell->face_midpoints[i][1]; + + return cell->face_lengths[i]; } /** @@ -82,6 +199,10 @@ __attribute__((always_inline)) INLINE float voronoi_get_face( * @param centroid Array to store the centroid in. */ __attribute__((always_inline)) INLINE void voronoi_get_centroid( - struct voronoi_cell *cell, float *centroid) {} + struct voronoi_cell *cell, float *centroid) { + + centroid[0] = cell->centroid[0]; + centroid[1] = cell->centroid[1]; +} #endif // SWIFT_VORONOIXD_ALGORITHM_H diff --git a/src/hydro/Shadowswift/voronoi2d_cell.h b/src/hydro/Shadowswift/voronoi2d_cell.h index d67cdaabfae4de30247f1460204dd612ff67e7ca..3c54ea8d0aa9ca1c915da1f245f889ebab2073d3 100644 --- a/src/hydro/Shadowswift/voronoi2d_cell.h +++ b/src/hydro/Shadowswift/voronoi2d_cell.h @@ -20,6 +20,10 @@ #ifndef SWIFT_VORONOIXD_CELL_H #define SWIFT_VORONOIXD_CELL_H +/* Maximal number of vertices (and neighbours) that can be stored in a + voronoi_cell struct. */ +#define VORONOI2D_MAXNUMVERT 100 + /* 2D Voronoi cell */ struct voronoi_cell { @@ -31,6 +35,24 @@ struct voronoi_cell { /* The centroid of the cell. */ float centroid[2]; + + /* Number of cell vertices (and neighbours). */ + int nvert; + + /* We only need to store one of these at the same time. */ + union { + /* The relative positions of the vertices of the cell. */ + float vertices[VORONOI2D_MAXNUMVERT][2]; + + /* The midpoints of the faces. */ + float face_midpoints[VORONOI2D_MAXNUMVERT][2]; + }; + + /* The ids of the neighbouring cells. */ + unsigned long long ngbs[VORONOI2D_MAXNUMVERT]; + + /* The lengths of the faces. */ + float face_lengths[VORONOI2D_MAXNUMVERT]; }; #endif // SWIFT_VORONOIXD_CELL_H diff --git a/src/hydro/Shadowswift/voronoi3d_algorithm.h b/src/hydro/Shadowswift/voronoi3d_algorithm.h index 6399a249fb55832435b2e74854d2269c385ae936..36f609befe1d1cd9835552990de8a3b836e4d58e 100644 --- a/src/hydro/Shadowswift/voronoi3d_algorithm.h +++ b/src/hydro/Shadowswift/voronoi3d_algorithm.h @@ -81,8 +81,8 @@ __attribute__((always_inline)) INLINE int check_counter(int *counter, extern float global_voronoi_box_anchor[3]; extern float global_voronoi_box_side[3]; -#define VORONOI3D_DECLARE_GLOBAL_VARIABLES() \ - float global_voronoi_box_anchor[3]; \ +#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 diff --git a/tests/benchmarkInteractions.c b/tests/benchmarkInteractions.c index da22b37f1df14acfa897ed6837de10bee5a0ce43..14648b0693714d668d35c0795503cc3ccfe9754d 100644 --- a/tests/benchmarkInteractions.c +++ b/tests/benchmarkInteractions.c @@ -450,8 +450,8 @@ void test_interactions(struct part test_part, struct part *parts, size_t count, #endif } -#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D) -VORONOI3D_DECLARE_GLOBAL_VARIABLES() +#if defined(SHADOWFAX_SPH) +VORONOI_DECLARE_GLOBAL_VARIABLES() #endif /* And go... */ diff --git a/tests/test125cells.c b/tests/test125cells.c index 6e8a3fd206142a64ac65a9004143f342677740c0..3ff812593abcbb1c7d9779b571ccf300cc184fe6 100644 --- a/tests/test125cells.c +++ b/tests/test125cells.c @@ -422,8 +422,8 @@ 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) && defined(HYDRO_DIMENSION_3D) -VORONOI3D_DECLARE_GLOBAL_VARIABLES() +#if defined(SHADOWFAX_SPH) +VORONOI_DECLARE_GLOBAL_VARIABLES() #endif /* And go... */ diff --git a/tests/test27cells.c b/tests/test27cells.c index 5cf7d7cdc881745dc6cb10fad3c6584381efdbb6..e5e73ad1e503b71c25bbe95cf5ede9010977c6b6 100644 --- a/tests/test27cells.c +++ b/tests/test27cells.c @@ -299,8 +299,8 @@ 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) && defined(HYDRO_DIMENSION_3D) -VORONOI3D_DECLARE_GLOBAL_VARIABLES() +#if defined(SHADOWFAX_SPH) +VORONOI_DECLARE_GLOBAL_VARIABLES() #endif /* And go... */ diff --git a/tests/testPair.c b/tests/testPair.c index bbcc17e61f1fcc78b7b9ecfad0d3db1857423df9..3719fb195581fade6b49b9c7a07e1060637f371f 100644 --- a/tests/testPair.c +++ b/tests/testPair.c @@ -187,8 +187,8 @@ 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) && defined(HYDRO_DIMENSION_3D) -VORONOI3D_DECLARE_GLOBAL_VARIABLES() +#if defined(SHADOWFAX_SPH) +VORONOI_DECLARE_GLOBAL_VARIABLES() #endif int main(int argc, char *argv[]) { diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c index 850e5aee2b38498d19d5cf37258a473fd6ff5d70..cbe857f3881ee2f0542afbf0cfd62f03185f3928 100644 --- a/tests/testSymmetry.c +++ b/tests/testSymmetry.c @@ -26,8 +26,8 @@ #include "swift.h" -#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D) -VORONOI3D_DECLARE_GLOBAL_VARIABLES() +#if defined(SHADOWFAX_SPH) +VORONOI_DECLARE_GLOBAL_VARIABLES() #endif int main(int argc, char *argv[]) { @@ -35,7 +35,7 @@ int main(int argc, char *argv[]) { /* Choke if need be */ feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); -#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D) +#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}; diff --git a/tests/testTimeIntegration.c b/tests/testTimeIntegration.c index 6ce59eccc62815ca10e04a8837d306c2f88290cf..185db6edde06ec82d0a0f87bbb166f05bd89e8b0 100644 --- a/tests/testTimeIntegration.c +++ b/tests/testTimeIntegration.c @@ -22,8 +22,8 @@ #include <stdlib.h> #include <string.h> -#if defined(SHADOWFAX_SPH) && defined(HYDRO_DIMENSION_3D) -VORONOI3D_DECLARE_GLOBAL_VARIABLES() +#if defined(SHADOWFAX_SPH) +VORONOI_DECLARE_GLOBAL_VARIABLES() #endif /** diff --git a/tests/testVoronoi2D.c b/tests/testVoronoi2D.c index f7605b1ede3ccfcb09b1da80a352c6b912ae7ea4..24210aa33d2c2e7c6e51fa06017ca21821733c53 100644 --- a/tests/testVoronoi2D.c +++ b/tests/testVoronoi2D.c @@ -19,4 +19,27 @@ #include "hydro/Shadowswift/voronoi2d_algorithm.h" -int main() { return 0; } +VORONOI_DECLARE_GLOBAL_VARIABLES() + +int main() { + + float anchor[3] = {-0.5f, -0.5f, -0.5f}; + float side[3] = {2.0f, 2.0f, 2.0f}; + voronoi_set_box(anchor, side); + + struct voronoi_cell cell; + double x[3] = {0.5, 0.5, 0.5}; + + voronoi_cell_init(&cell, x); + + float maxradius = voronoi_cell_finalize(&cell); + + assert(maxradius == 2.0f * sqrtf(2.0f)); + + assert(cell.volume == 4.0f); + + assert(cell.centroid[0] == 0.5f); + assert(cell.centroid[1] == 0.5f); + + return 0; +} diff --git a/tests/testVoronoi3D.c b/tests/testVoronoi3D.c index a1ff127c6b74aad7c9f502c874295f14d1231f95..80d45f7e01f647a3cd7158ac27960e8c9da377c4 100644 --- a/tests/testVoronoi3D.c +++ b/tests/testVoronoi3D.c @@ -1192,7 +1192,7 @@ void test_degeneracies() { #endif } -VORONOI3D_DECLARE_GLOBAL_VARIABLES() +VORONOI_DECLARE_GLOBAL_VARIABLES() int main() {