diff --git a/src/cell_drift.c b/src/cell_drift.c
index 008516fdb61275cb854341251cf7baf5cce41b48..73a7ad1c5cbacf984b7ece2f1bdde648847d5eb9 100644
--- a/src/cell_drift.c
+++ b/src/cell_drift.c
@@ -241,7 +241,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
}
#ifdef SHADOWFAX_SPH
- cell_malloc_delaunay_tesselation(c, &e->s->hs);
+ cell_malloc_delaunay_tessellation(c, &e->s->hs);
#endif
/* Now, get the maximal particle motion from its square */
diff --git a/src/cell_hydro.h b/src/cell_hydro.h
index 566e231cb2e17f74477c2af8c4de05af15e0bd6f..5f53ed8b09b17d0575edc775c70d612d8d9c5a92 100644
--- a/src/cell_hydro.h
+++ b/src/cell_hydro.h
@@ -45,6 +45,8 @@ struct cell_hydro {
struct xpart *xparts;
#ifdef SHADOWFAX_SPH
+ /*! Was Shadowfax functionality enabled for this cell? */
+ /* int shadowfax_enabled;*/
/*! Delaunay tessellation. */
struct delaunay deltess;
/*! Voronoi tessellation. */
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 42cb26e4c37ae6482f9b823b1fea2ccd3a018941..20f4f378c3dcfbc4964d037e2aeca261ba77c481 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -213,6 +213,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.wcount_dh = p->h / (1.1f * hnew - p->h);
return;
}
+ p->density.wcount = 1.0f;
volume = p->cell.volume;
#ifdef SWIFT_DEBUG_CHECKS
diff --git a/src/hydro/Shadowswift/hydro_part.h b/src/hydro/Shadowswift/hydro_part.h
index 343697f3473a9a4e6487ac3318f95271639de8a8..027b0e2ff72d60c3e207a51f7ee180c86b6774dd 100644
--- a/src/hydro/Shadowswift/hydro_part.h
+++ b/src/hydro/Shadowswift/hydro_part.h
@@ -217,6 +217,14 @@ struct part {
/* Voronoi cell. */
struct voronoi_cell cell;
+ /* New Voronoi cell. */
+ struct {
+ /* Flag that keeps track of where this particle was added to the Delaunay
+ tessellation. */
+ int flag;
+
+ } voronoi;
+
} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_SHADOWSWIFT_HYDRO_PART_H */
diff --git a/src/hydro_space.c b/src/hydro_space.c
index 9c5ab9fc619fca98dc22f15124bf4a6a63dcd539..982307474cdf526b0ab68292cb2aab2d02c14939 100644
--- a/src/hydro_space.c
+++ b/src/hydro_space.c
@@ -19,6 +19,7 @@
#include "hydro_space.h"
+#include "error.h"
#include "space.h"
/**
@@ -46,6 +47,11 @@ __attribute__((always_inline)) INLINE void hydro_space_init(
hs->side[1] = s->dim[1];
hs->side[2] = s->dim[2];
}
+
+ message(
+ "Initialised hydro space with anchor [%g, %g, %g] and side [%g, %g, %g].",
+ hs->anchor[0], hs->anchor[1], hs->anchor[2], hs->side[0], hs->side[1],
+ hs->side[2]);
}
#else
void hydro_space_init(struct hydro_space *hs, const struct space *s) {}
diff --git a/src/runner_doiact_functions_hydro.h b/src/runner_doiact_functions_hydro.h
index 4bdbcec05cf8fd9e5d5b06b13737d02886c94164..2d8bf4fbb958ae407b089d9cce6d69e8031ed43d 100644
--- a/src/runner_doiact_functions_hydro.h
+++ b/src/runner_doiact_functions_hydro.h
@@ -26,6 +26,10 @@
#include "runner_doiact_hydro.h"
+#ifdef SHADOWFAX_SPH
+#include "shadowfax/cell_shadowfax.h"
+#endif
+
/**
* @brief Compute the interactions between a cell pair (non-symmetric case).
*
@@ -69,6 +73,10 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
shift[k] = -e->s->dim[k];
}
+#if defined(SHADOWFAX_SPH) && (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ cell_shadowfax_do_pair1_naive(e, ci, cj, -1, shift);
+#endif
+
/* Loop over the parts in ci. */
for (int pid = 0; pid < count_i; pid++) {
@@ -195,6 +203,10 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
shift[k] = -e->s->dim[k];
}
+#if defined(SHADOWFAX_SPH) && (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ cell_shadowfax_do_pair2_naive(e, ci, cj, -1, shift);
+#endif
+
/* Loop over the parts in ci. */
for (int pid = 0; pid < count_i; pid++) {
@@ -585,6 +597,10 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci,
const float a = cosmo->a;
const float H = cosmo->H;
+#if defined(SHADOWFAX_SPH) && (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ cell_shadowfax_do_pair_subset_naive(e, ci, cj, -1, shift);
+#endif
+
/* Loop over the parts_i. */
for (int pid = 0; pid < count; pid++) {
@@ -686,6 +702,11 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
const struct sort_entry *sort_j = cell_get_hydro_sorts(cj, sid);
const float dxj = cj->hydro.dx_max_sort;
+#if defined(SHADOWFAX_SPH) && (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ cell_shadowfax_do_pair_subset(e, ci, parts_i, ind, count, cj, sid, flipped,
+ shift);
+#endif
+
/* Parts are on the left? */
if (!flipped) {
@@ -1044,6 +1065,10 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
const float a = cosmo->a;
const float H = cosmo->H;
+#if defined(SHADOWFAX_SPH) && (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ cell_shadowfax_do_pair1(e, ci, cj, sid, shift);
+#endif
+
if (cell_is_active_hydro(ci, e)) {
/* Loop over the parts in ci. */
@@ -1437,6 +1462,10 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
}
}
+#if defined(SHADOWFAX_SPH) && (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ cell_shadowfax_do_pair2(e, ci, cj, sid, shift);
+#endif
+
/* Loop over *all* the parts in ci starting from the centre until
we are out of range of anything in cj (using the maximal hi). */
for (int pid = count_i - 1;
@@ -1957,6 +1986,10 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
TIMER_TIC;
+#if defined(SHADOWFAX_SPH) && (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ cell_shadowfax_do_self1(e, c);
+#endif
+
struct part *restrict parts = c->hydro.parts;
const int count = c->hydro.count;
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index b98e47a0f31c4c5f235f6f20589597015934fb7e..ba12e952b46c9007de13bcfbdc2c91f2cac6066e 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -41,6 +41,10 @@
#include "timestep_limiter.h"
#include "tracers.h"
+#ifdef SHADOWFAX_SPH
+#include "shadowfax/cell_shadowfax.h"
+#endif
+
/* Import the density loop functions. */
#define FUNCTION density
#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY
@@ -1062,6 +1066,20 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
struct part *p = &parts[pid[i]];
struct xpart *xp = &xparts[pid[i]];
+#ifdef SHADOWFAX_SPH
+ float hnew = delaunay_get_search_radius(
+ &c->hydro.deltess, pid[i] + c->hydro.deltess.vertex_start);
+ if (hnew >= p->h) {
+ p->h *= 1.5f;
+ pid[redo] = pid[i];
+ h_0[redo] = h_0[i];
+ left[redo] = left[i];
+ right[redo] = right[i];
+ redo += 1;
+ }
+ continue;
+#endif
+
#ifdef SWIFT_DEBUG_CHECKS
/* Is this part within the timestep? */
if (!part_is_active(p, e)) error("Ghost applied to inactive particle");
@@ -1410,6 +1428,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Update h_max */
c->hydro.h_max = h_max;
+#ifdef SHADOWFAX_SPH
+ cell_shadowfax_end_density(c);
+#endif
+
/* The ghost may not always be at the top level.
* Therefore we need to update h_max between the super- and top-levels */
if (c->hydro.ghost) {
diff --git a/src/shadowfax/cell_shadowfax.h b/src/shadowfax/cell_shadowfax.h
index 88a55d3cf45d9a525da813a7297fd6961575d56f..aba28d9bfdde5d9c026c5bd67101b3a19c18c884 100644
--- a/src/shadowfax/cell_shadowfax.h
+++ b/src/shadowfax/cell_shadowfax.h
@@ -1,13 +1,397 @@
#ifndef SWIFT_CELL_SHADOWFAX_H
#define SWIFT_CELL_SHADOWFAX_H
+#include "active.h"
#include "cell.h"
#include "delaunay.h"
+#include "voronoi.h"
+
+__attribute__((always_inline)) INLINE static void shadowfax_flag_particle_added(
+ struct part *restrict p, int sid) {
+ p->voronoi.flag |= (1 << sid);
+}
+
+__attribute__((always_inline)) INLINE static int shadowfax_particle_was_added(
+ const struct part *restrict p, int sid) {
+ return (p->voronoi.flag & (1 << sid)) > 0;
+}
__attribute__((always_inline)) INLINE static void
-cell_malloc_delaunay_tesselation(struct cell *c, const struct hydro_space *hs) {
+cell_malloc_delaunay_tessellation(struct cell *c,
+ const struct hydro_space *hs) {
+
+ /* if(c->hydro.shadowfax_enabled){*/
+ delaunay_destroy(&c->hydro.deltess);
+ /* }*/
+ delaunay_init(&c->hydro.deltess, c->loc, c->width, c->hydro.count,
+ 10 * c->hydro.count);
+
+ const int count = c->hydro.count;
+ struct part *restrict parts = c->hydro.parts;
+
+ for (int pd = 0; pd < count; pd++) {
+ /* Get a pointer to the ith particle. */
+ struct part *restrict p = &parts[pd];
+ p->voronoi.flag = 0;
+ }
+
+ /* c->hydro.shadowfax_enabled = 1;*/
+}
+
+__attribute__((always_inline)) INLINE static void cell_shadowfax_do_self1(
+ const struct engine *e, struct cell *restrict c) {
+
+ const int count = c->hydro.count;
+ struct part *restrict parts = c->hydro.parts;
+
+ for (int pd = 0; pd < count; pd++) {
+ /* Get a pointer to the ith particle. */
+ struct part *restrict p = &parts[pd];
+ p->voronoi.flag = 0;
+ }
+
+ /* Loop over the parts in c. */
+ for (int pd = 0; pd < count; pd++) {
+
+ /* Get a pointer to the ith particle. */
+ struct part *restrict p = &parts[pd];
+
+ if (!shadowfax_particle_was_added(p, 0)) {
+ delaunay_add_local_vertex(&c->hydro.deltess, pd, p->x[0], p->x[1]);
+ shadowfax_flag_particle_added(p, 0);
+ }
+ }
+
+ /* delaunay_consolidate(&c->hydro.deltess);*/
+}
+
+__attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair_naive(
+ const struct engine *e, struct cell *restrict ci, struct cell *restrict cj,
+ int sid, const double *shift) {
+
+ if (sid < 0) {
+ error("Doing a naive interaction!");
+ }
+
+ const int count_i = ci->hydro.count;
+ const int count_j = cj->hydro.count;
+ struct part *restrict parts_i = ci->hydro.parts;
+ struct part *restrict parts_j = cj->hydro.parts;
+
+ /* Loop over the parts in cj. */
+ for (int pid = 0; pid < count_i; pid++) {
+
+ /* Get a pointer to the ith particle. */
+ struct part *restrict pi = &parts_i[pid];
+
+ if (!shadowfax_particle_was_added(pi, 1 + sid)) {
+ delaunay_add_new_vertex(&cj->hydro.deltess, pi->x[0] - shift[0],
+ pi->x[1] - shift[1]);
+ shadowfax_flag_particle_added(pi, 1 + sid);
+ }
+ }
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_j; pjd++) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts_j[pjd];
+
+ if (!shadowfax_particle_was_added(pj, 13 + sid)) {
+ delaunay_add_new_vertex(&ci->hydro.deltess, pj->x[0] + shift[0],
+ pj->x[1] + shift[1]);
+ shadowfax_flag_particle_added(pj, 13 + sid);
+ }
+ }
+}
+
+__attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair1_naive(
+ const struct engine *e, struct cell *restrict ci, struct cell *restrict cj,
+ int sid, const double *shift) {
+
+ if (ci == cj) error("Interacting cell with itself!");
+
+ cell_shadowfax_do_pair_naive(e, ci, cj, sid, shift);
+}
+
+__attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair2_naive(
+ const struct engine *e, struct cell *restrict ci, struct cell *restrict cj,
+ int sid, const double *shift) {
+
+ error("Shouldn't be using this function!");
+ if (ci == cj) error("Interacting cell with itself!");
+
+ cell_shadowfax_do_pair_naive(e, ci, cj, sid, shift);
+}
+
+__attribute__((always_inline)) INLINE static void
+cell_shadowfax_do_pair_subset_naive(const struct engine *e,
+ struct cell *restrict ci,
+ struct cell *restrict cj, int sid,
+ const double *shift) {
+
+ if (ci == cj) error("Interacting cell with itself!");
+
+ cell_shadowfax_do_pair_naive(e, ci, cj, sid, shift);
+}
+
+__attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair_subset(
+ const struct engine *e, struct cell *restrict ci,
+ struct part *restrict parts_i, int *restrict ind, int count,
+ struct cell *restrict cj, const int sid, const int flipped,
+ const double *shift) {
+
+ if (ci == cj) error("Interacting cell with itself!");
+
+ const int count_j = cj->hydro.count;
+ struct part *restrict parts_j = cj->hydro.parts;
+
+ /* Pick-out the sorted lists. */
+ const struct sort_entry *sort_j = cell_get_hydro_sorts(cj, sid);
+ const float dxj = cj->hydro.dx_max_sort;
+
+ /* Parts are on the left? */
+ if (!flipped) {
+
+ /* Loop over the parts_i. */
+ for (int pid = 0; pid < count; pid++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct part *restrict pi = &parts_i[ind[pid]];
+ const double pix = pi->x[0] - (shift[0]);
+ const double piy = pi->x[1] - (shift[1]);
+ const double piz = pi->x[2] - (shift[2]);
+ const float hi = pi->h;
+ const float hig2 = hi * hi * kernel_gamma2;
+ const double di = hi * kernel_gamma + dxj + pix * runner_shift[sid][0] +
+ piy * runner_shift[sid][1] + piz * runner_shift[sid][2];
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_j && sort_j[pjd].d < di; pjd++) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts_j[sort_j[pjd].i];
+
+ /* Skip inhibited particles. */
+ if (part_is_inhibited(pj, e)) continue;
+
+ const double pjx = pj->x[0];
+ const double pjy = pj->x[1];
+ const double pjz = pj->x[2];
+
+ /* Compute the pairwise distance. */
+ float dx[3] = {(float)(pix - pjx), (float)(piy - pjy),
+ (float)(piz - pjz)};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+ /* Hit or miss? */
+ if (r2 < hig2) {
+ if (!shadowfax_particle_was_added(pj, 1 + sid)) {
+ delaunay_add_new_vertex(&ci->hydro.deltess, pj->x[0] + shift[0],
+ pj->x[1] + shift[1]);
+ /* shadowfax_flag_particle_added(pj, 1+sid);*/
+ }
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the parts in ci. */
+ }
+
+ /* Parts are on the right. */
+ else {
+
+ /* Loop over the parts_i. */
+ for (int pid = 0; pid < count; pid++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct part *restrict pi = &parts_i[ind[pid]];
+ const double pix = pi->x[0] - (shift[0]);
+ const double piy = pi->x[1] - (shift[1]);
+ const double piz = pi->x[2] - (shift[2]);
+ const float hi = pi->h;
+ const float hig2 = hi * hi * kernel_gamma2;
+ const double di = -hi * kernel_gamma - dxj + pix * runner_shift[sid][0] +
+ piy * runner_shift[sid][1] + piz * runner_shift[sid][2];
+
+ /* Loop over the parts in cj. */
+ for (int pjd = count_j - 1; pjd >= 0 && di < sort_j[pjd].d; pjd--) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts_j[sort_j[pjd].i];
+
+ /* Skip inhibited particles. */
+ if (part_is_inhibited(pj, e)) continue;
+
+ const double pjx = pj->x[0];
+ const double pjy = pj->x[1];
+ const double pjz = pj->x[2];
+
+ /* Compute the pairwise distance. */
+ float dx[3] = {(float)(pix - pjx), (float)(piy - pjy),
+ (float)(piz - pjz)};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+ /* Hit or miss? */
+ if (r2 < hig2) {
+ if (!shadowfax_particle_was_added(pj, 13 + sid)) {
+ delaunay_add_new_vertex(&ci->hydro.deltess, pj->x[0] + shift[0],
+ pj->x[1] + shift[1]);
+ /* shadowfax_flag_particle_added(pj, 13+sid);*/
+ }
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the parts in ci. */
+ }
+}
+
+__attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair1(
+ const struct engine *e, struct cell *restrict ci, struct cell *restrict cj,
+ int sid, const double *shift) {
+
+ if (ci == cj) error("Interacting cell with itself!");
+
+ /* cell_shadowfax_do_pair_naive(e, ci, cj, sid, shift);*/
+ /* return;*/
+
+ /* Get the cutoff shift. */
+ double rshift = 0.0;
+ for (int k = 0; k < 3; k++) rshift += shift[k] * runner_shift[sid][k];
+
+ /* Pick-out the sorted lists. */
+ const struct sort_entry *restrict sort_i = cell_get_hydro_sorts(ci, sid);
+ const struct sort_entry *restrict sort_j = cell_get_hydro_sorts(cj, sid);
+
+ /* Get some other useful values. */
+ const double hi_max = ci->hydro.h_max * kernel_gamma - rshift;
+ const double hj_max = cj->hydro.h_max * kernel_gamma;
+ const int count_i = ci->hydro.count;
+ const int count_j = cj->hydro.count;
+ struct part *restrict parts_i = ci->hydro.parts;
+ struct part *restrict parts_j = cj->hydro.parts;
+ const double di_max = sort_i[count_i - 1].d - rshift;
+ const double dj_min = sort_j[0].d;
+ const float dx_max = (ci->hydro.dx_max_sort + cj->hydro.dx_max_sort);
+
+ if (cell_is_active_hydro(ci, e)) {
+
+ /* Loop over the parts in ci. */
+ for (int pid = count_i - 1;
+ pid >= 0 && sort_i[pid].d + hi_max + dx_max > dj_min; pid--) {
+
+ /* Get a hold of the ith part in ci. */
+ struct part *restrict pi = &parts_i[sort_i[pid].i];
+ const float hi = pi->h;
+
+ /* Skip inactive particles */
+ if (!part_is_active(pi, e)) continue;
+
+ /* Is there anything we need to interact with ? */
+ const double di = sort_i[pid].d + hi * kernel_gamma + dx_max - rshift;
+ if (di < dj_min) continue;
+
+ /* Get some additional information about pi */
+ const float hig2 = hi * hi * kernel_gamma2;
+ const float pix = pi->x[0] - (cj->loc[0] + shift[0]);
+ const float piy = pi->x[1] - (cj->loc[1] + shift[1]);
+ const float piz = pi->x[2] - (cj->loc[2] + shift[2]);
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_j && sort_j[pjd].d < di; pjd++) {
+
+ /* Recover pj */
+ struct part *pj = &parts_j[sort_j[pjd].i];
+
+ /* Skip inhibited particles. */
+ if (part_is_inhibited(pj, e)) continue;
+
+ const float pjx = pj->x[0] - cj->loc[0];
+ const float pjy = pj->x[1] - cj->loc[1];
+ const float pjz = pj->x[2] - cj->loc[2];
+
+ /* Compute the pairwise distance. */
+ float dx[3] = {pix - pjx, piy - pjy, piz - pjz};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+ /* Hit or miss? */
+ if (r2 < hig2) {
+
+ if (!shadowfax_particle_was_added(pj, 13 + sid)) {
+ delaunay_add_new_vertex(&ci->hydro.deltess, pj->x[0] + shift[0],
+ pj->x[1] + shift[1]);
+ shadowfax_flag_particle_added(pj, 13 + sid);
+ }
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the parts in ci. */
+ } /* Cell ci is active */
+
+ if (cell_is_active_hydro(cj, e)) {
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_j && sort_j[pjd].d - hj_max - dx_max < di_max;
+ pjd++) {
+
+ /* Get a hold of the jth part in cj. */
+ struct part *pj = &parts_j[sort_j[pjd].i];
+ const float hj = pj->h;
+
+ /* Skip inactive particles */
+ if (!part_is_active(pj, e)) continue;
+
+ /* Is there anything we need to interact with ? */
+ const double dj = sort_j[pjd].d - hj * kernel_gamma - dx_max + rshift;
+ if (dj - rshift > di_max) continue;
+
+ /* Get some additional information about pj */
+ const float hjg2 = hj * hj * kernel_gamma2;
+ const float pjx = pj->x[0] - cj->loc[0];
+ const float pjy = pj->x[1] - cj->loc[1];
+ const float pjz = pj->x[2] - cj->loc[2];
+
+ /* Loop over the parts in ci. */
+ for (int pid = count_i - 1; pid >= 0 && sort_i[pid].d > dj; pid--) {
+
+ /* Recover pi */
+ struct part *pi = &parts_i[sort_i[pid].i];
+
+ /* Skip inhibited particles. */
+ if (part_is_inhibited(pi, e)) continue;
+
+ const float pix = pi->x[0] - (cj->loc[0] + shift[0]);
+ const float piy = pi->x[1] - (cj->loc[1] + shift[1]);
+ const float piz = pi->x[2] - (cj->loc[2] + shift[2]);
+
+ /* Compute the pairwise distance. */
+ float dx[3] = {pjx - pix, pjy - piy, pjz - piz};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+ /* Hit or miss? */
+ if (r2 < hjg2) {
+
+ if (!shadowfax_particle_was_added(pi, 1 + sid)) {
+ delaunay_add_new_vertex(&cj->hydro.deltess, pi->x[0] - shift[0],
+ pi->x[1] - shift[1]);
+ shadowfax_flag_particle_added(pi, 1 + sid);
+ }
+ }
+ } /* loop over the parts in ci. */
+ } /* loop over the parts in cj. */
+ } /* Cell cj is active */
+}
+
+__attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair2(
+ const struct engine *e, struct cell *restrict ci, struct cell *restrict cj,
+ int sid, const double *shift) {
+
+ error("Shouldn't be using this function!");
+ if (ci == cj) error("Interacting cell with itself!");
+
+ cell_shadowfax_do_pair_naive(e, ci, cj, sid, shift);
+}
+
+__attribute__((always_inline)) INLINE static void cell_shadowfax_end_density(
+ struct cell *restrict c) {
- delaunay_init(&c->hydro.deltess, hs, 100, 100);
+ voronoi_init(&c->hydro.vortess, &c->hydro.deltess);
}
#endif /* SWIFT_CELL_SHADOWFAX_H */
diff --git a/src/shadowfax/delaunay.h b/src/shadowfax/delaunay.h
index b9b5ab841f0530415b1a670548d778bec28f5d2d..9523b13572352a3e7d96ddc68fc8f8a997260eee 100644
--- a/src/shadowfax/delaunay.h
+++ b/src/shadowfax/delaunay.h
@@ -28,8 +28,10 @@
#ifndef SWIFT_DELAUNAY_H
#define SWIFT_DELAUNAY_H
+#include "error.h"
#include "geometry.h"
#include "hydro_space.h"
+#include "memuse.h"
#include "triangle.h"
#include <float.h>
@@ -94,6 +96,9 @@
*/
struct delaunay {
+ /* Activity flag, useful for debugging. */
+ int active;
+
/* Box geometry: used to set up the initial vertices and triangles and to
* convert input coordinates to integer coordinates. */
@@ -231,6 +236,45 @@ static inline unsigned long int delaunay_double_to_int(double d) {
return (u.ull & 0xFFFFFFFFFFFFFllu);
}
+inline static void delaunay_init_vertex(struct delaunay* restrict d, int v,
+ double x, double y) {
+
+ /* store a copy of the vertex coordinates (we should get rid of this for
+ SWIFT) */
+ d->vertices[2 * v] = x;
+ d->vertices[2 * v + 1] = y;
+
+ /* compute the rescaled coordinates. We do this because floating point values
+ in the range [1,2[ all have the same exponent (0), which guarantees that
+ their mantissas form a linear sequence */
+ double rescaled_x = 1. + (x - d->anchor[0]) * d->inverse_side;
+ double rescaled_y = 1. + (y - d->anchor[1]) * d->inverse_side;
+
+ delaunay_assert(rescaled_x >= 1.);
+ delaunay_assert(rescaled_x < 2.);
+ delaunay_assert(rescaled_y >= 1.);
+ delaunay_assert(rescaled_y < 2.);
+
+#ifdef DELAUNAY_NONEXACT
+ /* store a copy of the rescaled coordinates to apply non-exact tests */
+ d->rescaled_vertices[2 * v] = rescaled_x;
+ d->rescaled_vertices[2 * v + 1] = rescaled_y;
+#endif
+
+ /* convert the rescaled coordinates to integer coordinates and store these */
+ d->integer_vertices[2 * v] = delaunay_double_to_int(rescaled_x);
+ d->integer_vertices[2 * v + 1] = delaunay_double_to_int(rescaled_y);
+
+ /* initialise the variables that keep track of the link between vertices and
+ triangles.
+ We use negative values so that we can later detect missing links. */
+ d->vertex_triangles[v] = -1;
+ d->vertex_triangle_index[v] = -1;
+
+ /* initialise the search radii to the largest possible value */
+ d->search_radii[v] = DBL_MAX;
+}
+
/**
* @brief Add a new vertex with the given coordinates.
*
@@ -252,57 +296,28 @@ inline static int delaunay_new_vertex(struct delaunay* restrict d, double x,
if (d->vertex_index == d->vertex_size) {
/* dynamically grow the size of the arrays with a factor 2 */
d->vertex_size <<= 1;
- d->vertices =
- (double*)realloc(d->vertices, d->vertex_size * 2 * sizeof(double));
+ d->vertices = (double*)swift_realloc("delaunay vertices", d->vertices,
+ d->vertex_size * 2 * sizeof(double));
#ifdef DELAUNAY_NONEXACT
- d->rescaled_vertices = (double*)realloc(
- d->rescaled_vertices, d->vertex_size * 2 * sizeof(double));
+ d->rescaled_vertices = (double*)swift_realloc(
+ "delaunay rescaled vertices", d->rescaled_vertices,
+ d->vertex_size * 2 * sizeof(double));
#endif
- d->integer_vertices = (unsigned long int*)realloc(
- d->integer_vertices, d->vertex_size * 2 * sizeof(unsigned long int));
+ d->integer_vertices = (unsigned long int*)swift_realloc(
+ "delaunay integer vertices", d->integer_vertices,
+ d->vertex_size * 2 * sizeof(unsigned long int));
d->vertex_triangles =
- (int*)realloc(d->vertex_triangles, d->vertex_size * sizeof(int));
- d->vertex_triangle_index =
- (int*)realloc(d->vertex_triangle_index, d->vertex_size * sizeof(int));
+ (int*)swift_realloc("delaunay vertex-triangle links",
+ d->vertex_triangles, d->vertex_size * sizeof(int));
+ d->vertex_triangle_index = (int*)swift_realloc(
+ "delaunay vertex-triangle index links", d->vertex_triangle_index,
+ d->vertex_size * sizeof(int));
d->search_radii =
- (double*)realloc(d->search_radii, d->vertex_size * sizeof(double));
+ (double*)swift_realloc("delaunay search radii", d->search_radii,
+ d->vertex_size * sizeof(double));
}
- /* store a copy of the vertex coordinates (we should get rid of this for
- SWIFT) */
- d->vertices[2 * d->vertex_index] = x;
- d->vertices[2 * d->vertex_index + 1] = y;
-
- /* compute the rescaled coordinates. We do this because floating point values
- in the range [1,2[ all have the same exponent (0), which guarantees that
- their mantissas form a linear sequence */
- double rescaled_x = 1. + (x - d->anchor[0]) * d->inverse_side;
- double rescaled_y = 1. + (y - d->anchor[1]) * d->inverse_side;
-
- delaunay_assert(rescaled_x >= 1.);
- delaunay_assert(rescaled_x < 2.);
- delaunay_assert(rescaled_y >= 1.);
- delaunay_assert(rescaled_y < 2.);
-
-#ifdef DELAUNAY_NONEXACT
- /* store a copy of the rescaled coordinates to apply non-exact tests */
- d->rescaled_vertices[2 * d->vertex_index] = rescaled_x;
- d->rescaled_vertices[2 * d->vertex_index + 1] = rescaled_y;
-#endif
-
- /* convert the rescaled coordinates to integer coordinates and store these */
- d->integer_vertices[2 * d->vertex_index] = delaunay_double_to_int(rescaled_x);
- d->integer_vertices[2 * d->vertex_index + 1] =
- delaunay_double_to_int(rescaled_y);
-
- /* initialise the variables that keep track of the link between vertices and
- triangles.
- We use negative values so that we can later detect missing links. */
- d->vertex_triangles[d->vertex_index] = -1;
- d->vertex_triangle_index[d->vertex_index] = -1;
-
- /* initialise the search radii to the largest possible value */
- d->search_radii[d->vertex_index] = DBL_MAX;
+ delaunay_init_vertex(d, d->vertex_index, x, y);
/* return the vertex index and then increase it by 1.
After this operation, vertex_index will correspond to the size of the
@@ -328,8 +343,9 @@ inline static int delaunay_new_triangle(struct delaunay* restrict d) {
/* no: increase the size of the triangle array with a factor 2 and
reallocate it in memory */
d->triangle_size <<= 1;
- d->triangles = (struct triangle*)realloc(
- d->triangles, d->triangle_size * sizeof(struct triangle));
+ d->triangles = (struct triangle*)swift_realloc(
+ "delaunay triangles", d->triangles,
+ d->triangle_size * sizeof(struct triangle));
}
/* return the triangle index and then increase it by 1.
@@ -522,32 +538,49 @@ inline static void delaunay_check_tessellation(struct delaunay* restrict d) {
* @param triangle_size Initial size of the triangle array.
*/
inline static void delaunay_init(struct delaunay* restrict d,
- const struct hydro_space* restrict hs,
- int vertex_size, int triangle_size) {
+ const double* cell_loc,
+ const double* cell_width, int vertex_size,
+ int triangle_size) {
+
+ if (d->active != 0) {
+ error("Delaunay tessellation corruption!");
+ }
+
+ if (vertex_size == 0) {
+ /* Don't bother setting up a Delaunay tessellation for emtpy cells */
+ return;
+ }
+
+ d->active = 1;
/* allocate memory for the vertex arrays */
- d->vertices = (double*)malloc(vertex_size * 2 * sizeof(double));
+ d->vertices = (double*)swift_malloc("delaunay vertices",
+ vertex_size * 2 * sizeof(double));
#ifdef DELAUNAY_NONEXACT
- d->rescaled_vertices = (double*)malloc(vertex_size * 2 * sizeof(double));
+ d->rescaled_vertices = (double*)swift_malloc(
+ "delaunay rescaled vertices", vertex_size * 2 * sizeof(double));
#endif
- d->integer_vertices =
- (unsigned long int*)malloc(vertex_size * 2 * sizeof(unsigned long int));
- d->vertex_triangles = (int*)malloc(vertex_size * sizeof(int));
- d->vertex_triangle_index = (int*)malloc(vertex_size * sizeof(int));
- d->search_radii = (double*)malloc(vertex_size * sizeof(double));
- d->vertex_index = 0;
+ d->integer_vertices = (unsigned long int*)swift_malloc(
+ "delaunay integer vertices", vertex_size * 2 * sizeof(unsigned long int));
+ d->vertex_triangles = (int*)swift_malloc("delaunay vertex-triangle links",
+ vertex_size * sizeof(int));
+ d->vertex_triangle_index = (int*)swift_malloc(
+ "delaunay vertex-triangle index links", vertex_size * sizeof(int));
+ d->search_radii = (double*)swift_malloc("delaunay search radii",
+ vertex_size * sizeof(double));
+ d->vertex_index = vertex_size;
d->vertex_size = vertex_size;
/* allocate memory for the triangle array */
- d->triangles =
- (struct triangle*)malloc(triangle_size * sizeof(struct triangle));
+ d->triangles = (struct triangle*)swift_malloc(
+ "delaunay triangles", triangle_size * sizeof(struct triangle));
d->triangle_index = 0;
d->triangle_size = triangle_size;
/* allocate memory for the queue (note that the queue size of 10 was chosen
arbitrarily, and a proper value should be chosen based on performance
measurements) */
- d->queue = (int*)malloc(10 * sizeof(int));
+ d->queue = (int*)swift_malloc("delaunay queue", 10 * sizeof(int));
d->queue_index = 0;
d->queue_size = 10;
@@ -557,9 +590,9 @@ inline static void delaunay_init(struct delaunay* restrict d,
work for even the very degenerate case of a single vertex that could be
anywhere in the box. Also note that we convert the generally rectangular
box to a square. */
- double box_anchor[2] = {hs->anchor[0] - hs->side[0],
- hs->anchor[1] - hs->side[1]};
- double box_side = 6. * fmax(hs->side[0], hs->side[1]);
+ double box_anchor[2] = {cell_loc[0] - cell_width[0],
+ cell_loc[1] - cell_width[1]};
+ double box_side = 6. * fmax(cell_width[0], cell_width[1]);
/* store the anchor and inverse side_length for the conversion from box
coordinates to rescaled (integer) coordinates */
@@ -615,13 +648,13 @@ inline static void delaunay_init(struct delaunay* restrict d,
d->last_triangle = first_triangle;
- d->vertex_start = d->vertex_index;
+ d->vertex_start = 0;
/* initialise the last vertex index to a negative value to signal that the
Delaunay tessellation was not consolidated yet.
delaunay_update_search_radii() will not work until delaunay_consolidate()
was called. Neither will it be possible to convert the Delaunay
tessellation into a Voronoi grid before this happens. */
- d->vertex_end = -1;
+ d->vertex_end = vertex_size;
/* Perform potential log output and sanity checks */
delaunay_log("Post init check");
@@ -634,17 +667,21 @@ inline static void delaunay_init(struct delaunay* restrict d,
* @param d Delaunay tessellation.
*/
inline static void delaunay_destroy(struct delaunay* restrict d) {
- free(d->vertices);
+ d->active = 0;
+ if (d->vertices != NULL) {
+ swift_free("delaunay vertices", d->vertices);
#ifdef DELAUNAY_NONEXACT
- free(d->rescaled_vertices);
+ swift_free("delaunay rescaled vertices", d->rescaled_vertices);
#endif
- free(d->integer_vertices);
- free(d->vertex_triangles);
- free(d->vertex_triangle_index);
- free(d->search_radii);
- free(d->triangles);
- free(d->queue);
- geometry_destroy(&d->geometry);
+ swift_free("delaunay integer vertices", d->integer_vertices);
+ swift_free("delaunay vertex-triangle links", d->vertex_triangles);
+ swift_free("delaunay vertex-triangle index links",
+ d->vertex_triangle_index);
+ swift_free("delaunay search radii", d->search_radii);
+ swift_free("delaunay triangles", d->triangles);
+ swift_free("delaunay queue", d->queue);
+ geometry_destroy(&d->geometry);
+ }
}
/**
@@ -692,7 +729,7 @@ inline static double delaunay_get_radius(const struct delaunay* restrict d,
* @param d Delaunay tessellation.
*/
inline static void delaunay_consolidate(struct delaunay* restrict d) {
- d->vertex_end = d->vertex_index;
+ d->vertex_end = d->vertex_index - 3;
}
/**
@@ -738,6 +775,23 @@ inline static int delaunay_update_search_radii(struct delaunay* restrict d,
return count;
}
+inline static float delaunay_get_search_radius(struct delaunay* restrict d,
+ int vi) {
+ int t0 = d->vertex_triangles[vi];
+ int vi0 = d->vertex_triangle_index[vi];
+ int vi0p1 = (vi0 + 1) % 3;
+ float r = 2.0f * delaunay_get_radius(d, t0);
+ int t1 = d->triangles[t0].neighbours[vi0p1];
+ int vi1 = d->triangles[t0].index_in_neighbour[vi0p1];
+ while (t1 != t0) {
+ r = max(r, 2.f * delaunay_get_radius(d, t1));
+ int vi1p2 = (vi1 + 2) % 3;
+ vi1 = d->triangles[t1].index_in_neighbour[vi1p2];
+ t1 = d->triangles[t1].neighbours[vi1p2];
+ }
+ return r;
+}
+
/**
* @brief Randomly choose a neighbour from the given two options.
*
@@ -911,14 +965,17 @@ inline static int delaunay_test_point_inside_triangle(
case 26:
/* testi0 is zero */
*t_new = d->triangles[t].neighbours[0];
+ *ngb_index = 0;
return 2;
case 38:
/* testi1 is zero */
*t_new = d->triangles[t].neighbours[1];
+ *ngb_index = 1;
return 2;
case 41:
/* testi2 is zero */
*t_new = d->triangles[t].neighbours[2];
+ *ngb_index = 2;
return 2;
case 42:
/* all tests returned 1 (testsum is 32 + 8 + 2) */
@@ -928,8 +985,14 @@ inline static int delaunay_test_point_inside_triangle(
/* we have ended up in a geometrically impossible scenario. This can
happen if the triangle vertices are colinear, or if one or multiple
vertices coincide. */
- fprintf(stderr, "Impossible case (%i)!\n", testsum);
- abort();
+ return -1;
+ /* fprintf(stderr, "Impossible case (%i)!\n", testsum);*/
+ /* fprintf(stderr, "Vertex positions:\n");*/
+ /* fprintf(stderr, "a: %lu %lu\n", aix, aiy);*/
+ /* fprintf(stderr, "b: %lu %lu\n", bix, biy);*/
+ /* fprintf(stderr, "c: %lu %lu\n", cix, ciy);*/
+ /* fprintf(stderr, "d: %lu %lu\n", dix, diy);*/
+ /* error("Impossible case during Delaunay construction!");*/
}
}
@@ -945,7 +1008,8 @@ inline static void delaunay_enqueue(struct delaunay* restrict d, int t) {
if (d->queue_index == d->queue_size) {
/* there isn't: increase the size of the queue with a factor 2. */
d->queue_size <<= 1;
- d->queue = (int*)realloc(d->queue, d->queue_size * sizeof(int));
+ d->queue = (int*)swift_realloc("delaunay queue", d->queue,
+ d->queue_size * sizeof(int));
}
delaunay_log("Enqueuing triangle %i and vertex 2", t);
@@ -1164,17 +1228,19 @@ inline static void delaunay_check_triangles(struct delaunay* restrict d) {
* vertex should always be the final vertex of any newly created triangle.
*
* @param d Delaunay tessellation.
+ * @param v Vertex index.
* @param x Horizontal coordinate of the new vertex.
* @param y Vertical coordinate of the new vertex.
+ * @return 0 on success, -1 if the vertex already exists.
*/
-inline static void delaunay_add_vertex(struct delaunay* restrict d, double x,
- double y) {
+inline static int delaunay_add_vertex(struct delaunay* restrict d, int v,
+ double x, double y) {
- delaunay_log("Adding vertex with position %g %g", x, y);
+ if (d->active != 1) {
+ error("Trying to add a vertex to an inactive Delaunay tessellation!");
+ }
- /* create the new vertex */
- int v = delaunay_new_vertex(d, x, y);
- delaunay_log("Created new vertex with index %i", v);
+ delaunay_log("Adding vertex with position %g %g", x, y);
/* find the triangle that contains the new vertex. Starting from an initial
guess (the last triangle accessed by the previous insertion), we test if
@@ -1185,10 +1251,16 @@ inline static void delaunay_add_vertex(struct delaunay* restrict d, double x,
int t0, t1, ngb_index;
t0 = d->last_triangle;
int flag = delaunay_test_point_inside_triangle(d, v, t0, &t1, &ngb_index);
+ if (flag == -1) {
+ return -1;
+ }
int count = 0;
while (flag == 0) {
t0 = t1;
flag = delaunay_test_point_inside_triangle(d, v, t0, &t1, &ngb_index);
+ if (flag == -1) {
+ return -1;
+ }
++count;
delaunay_assert(count < d->triangle_index);
}
@@ -1289,10 +1361,10 @@ inline static void delaunay_add_vertex(struct delaunay* restrict d, double x,
int ngb0_2 = d->triangles[t0].neighbours[i0_2];
int ngbi0_2 = d->triangles[t0].index_in_neighbour[i0_2];
- int ngb1_1 = d->triangles[t0].neighbours[i1_1];
- int ngbi1_1 = d->triangles[t0].index_in_neighbour[i1_1];
- int ngb1_2 = d->triangles[t0].neighbours[i1_2];
- int ngbi1_2 = d->triangles[t0].index_in_neighbour[i1_2];
+ int ngb1_1 = d->triangles[t1].neighbours[i1_1];
+ int ngbi1_1 = d->triangles[t1].index_in_neighbour[i1_1];
+ int ngb1_2 = d->triangles[t1].neighbours[i1_2];
+ int ngbi1_2 = d->triangles[t1].index_in_neighbour[i1_2];
int t2 = delaunay_new_triangle(d);
int t3 = delaunay_new_triangle(d);
@@ -1354,6 +1426,46 @@ inline static void delaunay_add_vertex(struct delaunay* restrict d, double x,
delaunay_log("Post vertex %i check", v);
/* perform a consistency test if enabled */
delaunay_check_tessellation(d);
+ return 0;
+}
+
+inline static void delaunay_add_local_vertex(struct delaunay* restrict d, int v,
+ double x, double y) {
+ delaunay_init_vertex(d, v, x, y);
+ if (delaunay_add_vertex(d, v, x, y) != 0) {
+ error("Local vertices cannot be added twice!");
+ }
+}
+
+inline static void delaunay_add_new_vertex(struct delaunay* restrict d,
+ double x, double y) {
+ if (d->active != 1) {
+ error("Trying to add a vertex to an inactive Delaunay tessellation!");
+ }
+
+ /* create the new vertex */
+ int v = delaunay_new_vertex(d, x, y);
+ delaunay_log("Created new vertex with index %i", v);
+ int flag = delaunay_add_vertex(d, v, x, y);
+ if (flag == -1) {
+ /* vertex already exists. Delete the last vertex. */
+ --d->vertex_index;
+ }
+}
+
+inline static void delaunay_write_tessellation(
+ const struct delaunay* restrict d, FILE* file, size_t* offset) {
+
+ for (int i = 0; i < d->vertex_index; ++i) {
+ fprintf(file, "V\t%lu\t%g\t%g\n", *offset + i, d->vertices[2 * i],
+ d->vertices[2 * i + 1]);
+ }
+ for (int i = 3; i < d->triangle_index; ++i) {
+ fprintf(file, "T\t%lu\t%lu\t%lu\n", *offset + d->triangles[i].vertices[0],
+ *offset + d->triangles[i].vertices[1],
+ *offset + d->triangles[i].vertices[2]);
+ }
+ *offset += d->vertex_index;
}
/**
@@ -1376,14 +1488,8 @@ inline static void delaunay_print_tessellation(
FILE* file = fopen(file_name, "w");
- for (int i = 0; i < d->vertex_index; ++i) {
- fprintf(file, "V\t%i\t%g\t%g\n", i, d->vertices[2 * i],
- d->vertices[2 * i + 1]);
- }
- for (int i = 0; i < d->triangle_index; ++i) {
- fprintf(file, "T\t%i\t%i\t%i\n", d->triangles[i].vertices[0],
- d->triangles[i].vertices[1], d->triangles[i].vertices[2]);
- }
+ size_t offset = 0;
+ delaunay_write_tessellation(d, file, &offset);
fclose(file);
}
diff --git a/src/shadowfax/voronoi.h b/src/shadowfax/voronoi.h
index ef8f0fecbb2f5dc743165ba95b08807dc02d4d74..017d6a78d203c66ea3da628297d55d68df3f2261 100644
--- a/src/shadowfax/voronoi.h
+++ b/src/shadowfax/voronoi.h
@@ -364,6 +364,36 @@ static inline void voronoi_check_grid(const struct voronoi *restrict v) {
printf("Total volume: %g\n", V);
}
+static inline void voronoi_write_grid(const struct voronoi *restrict v,
+ FILE *file, size_t *offset) {
+
+ for (int i = 0; i < v->number_of_cells; ++i) {
+ fprintf(file, "G\t%g\t%g\n", v->generators[2 * i],
+ v->generators[2 * i + 1]);
+ fprintf(file, "M\t%g\t%g\n", v->cell_centroid[2 * i],
+ v->cell_centroid[2 * i + 1]);
+ for (int j = 1; j < v->vertex_number[i]; ++j) {
+ int cjm1 = v->vertex_offset[i] + j - 1;
+ int cj = v->vertex_offset[i] + j;
+ fprintf(file, "C\t%lu\t%lu\n", *offset + v->connections[cjm1],
+ *offset + v->connections[cj]);
+ fprintf(file, "F\t%g\t%g\n", v->face_midpoints[2 * cj],
+ v->face_midpoints[2 * cj + 1]);
+ }
+ fprintf(
+ file, "C\t%lu\t%lu\n",
+ *offset + v->connections[v->vertex_offset[i] + v->vertex_number[i] - 1],
+ *offset + v->connections[v->vertex_offset[i]]);
+ fprintf(file, "F\t%g\t%g\n", v->face_midpoints[2 * v->vertex_offset[i]],
+ v->face_midpoints[2 * v->vertex_offset[i] + 1]);
+ }
+ for (int i = 0; i < v->vertex_index; ++i) {
+ fprintf(file, "V\t%g\t%g\n", v->vertices[2 * i], v->vertices[2 * i + 1]);
+ }
+
+ *offset += v->vertex_index;
+}
+
/**
* @brief Print the Voronoi grid to a file with the given name.
*
@@ -385,27 +415,8 @@ static inline void voronoi_print_grid(const struct voronoi *restrict v,
FILE *file = fopen(file_name, "w");
- for (int i = 0; i < v->number_of_cells; ++i) {
- fprintf(file, "G\t%g\t%g\n", v->generators[2 * i],
- v->generators[2 * i + 1]);
- fprintf(file, "M\t%g\t%g\n", v->cell_centroid[2 * i],
- v->cell_centroid[2 * i + 1]);
- for (int j = 1; j < v->vertex_number[i]; ++j) {
- int cjm1 = v->vertex_offset[i] + j - 1;
- int cj = v->vertex_offset[i] + j;
- fprintf(file, "C\t%i\t%i\n", v->connections[cjm1], v->connections[cj]);
- fprintf(file, "F\t%g\t%g\n", v->face_midpoints[2 * cj],
- v->face_midpoints[2 * cj + 1]);
- }
- fprintf(file, "C\t%i\t%i\n",
- v->connections[v->vertex_offset[i] + v->vertex_number[i] - 1],
- v->connections[v->vertex_offset[i]]);
- fprintf(file, "F\t%g\t%g\n", v->face_midpoints[2 * v->vertex_offset[i]],
- v->face_midpoints[2 * v->vertex_offset[i] + 1]);
- }
- for (int i = 0; i < v->vertex_index; ++i) {
- fprintf(file, "V\t%g\t%g\n", v->vertices[2 * i], v->vertices[2 * i + 1]);
- }
+ size_t offset = 0;
+ voronoi_write_grid(v, file, &offset);
fclose(file);
}
diff --git a/src/space_rebuild.c b/src/space_rebuild.c
index 0b68ca2b0cf9cfaea16ba95908fa38cacff71a25..4fb2ebfcb754ce824d575af2a325982ffb88aad3 100644
--- a/src/space_rebuild.c
+++ b/src/space_rebuild.c
@@ -30,6 +30,10 @@
#include "engine.h"
#include "memswap.h"
+#ifdef SHADOWFAX_SPH
+#include "shadowfax/cell_shadowfax.h"
+#endif
+
/*! Expected maximal number of strays received at a rebuild */
extern int space_expected_max_nr_strays;
@@ -955,6 +959,10 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
s->local_cells_with_particles_top[s->nr_local_cells_with_particles] = k;
s->nr_local_cells_with_particles++;
}
+
+#ifdef SHADOWFAX_SPH
+ cell_malloc_delaunay_tessellation(c, &s->hs);
+#endif
}
if (verbose) {
message("Have %d local top-level cells with particles (total=%d)",