diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index b8376c3d7a843651870e0d6afaa2508c98a5d56b..e2327cfff6159cec89092535687f9b8eb5931fe5 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -20,6 +20,7 @@
 #include <float.h>
 #include "adiabatic_index.h"
 #include "approx_math.h"
+#include "equation_of_state.h"
 #include "hydro_gradients.h"
 #include "voronoi_algorithm.h"
 
@@ -437,3 +438,35 @@ __attribute__((always_inline)) INLINE static float hydro_get_density(
 
   return p->primitives.rho;
 }
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed internal
+ * energy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param u The new internal energy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
+    struct part* restrict p, float u) {
+
+  p->conserved.energy = u * p->conserved.mass;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed entropy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param S The new entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_entropy(
+    struct part* restrict p, float S) {
+
+  p->conserved.energy = gas_internal_energy_from_entropy(p->primitives.rho, S) *
+                        p->conserved.mass;
+}
diff --git a/src/hydro/Shadowswift/voronoi3d_algorithm.h b/src/hydro/Shadowswift/voronoi3d_algorithm.h
index 47ec42656922ac1701ccdfb30770973a88959693..2f1208245ab7600d7acd6e3e3a9b7feab4b48876 100644
--- a/src/hydro/Shadowswift/voronoi3d_algorithm.h
+++ b/src/hydro/Shadowswift/voronoi3d_algorithm.h
@@ -30,7 +30,7 @@
 
 /* Tolerance parameter used to decide when to use more precise geometric
    criteria */
-#define VORONOI3D_TOLERANCE 1.e-5f
+#define VORONOI3D_TOLERANCE 1.e-6f
 
 /* Box boundary flags used to signal cells neighbouring the box boundary
    These values correspond to the top range of possible 64-bit integers, and
@@ -48,12 +48,68 @@
 /* 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 -1.0f
-#define VORONOI3D_BOX_ANCHOR_Y -1.0f
-#define VORONOI3D_BOX_ANCHOR_Z -1.0f
-#define VORONOI3D_BOX_SIDE_X 3.0f
-#define VORONOI3D_BOX_SIDE_Y 3.0f
-#define VORONOI3D_BOX_SIDE_Z 3.0f
+#define VORONOI3D_BOX_ANCHOR_X -2.0f
+#define VORONOI3D_BOX_ANCHOR_Y -2.0f
+#define VORONOI3D_BOX_ANCHOR_Z -2.0f
+#define VORONOI3D_BOX_SIDE_X 6.0f
+#define VORONOI3D_BOX_SIDE_Y 6.0f
+#define VORONOI3D_BOX_SIDE_Z 6.0f
+
+__attribute__((always_inline)) INLINE static float voronoi_get_box_volume() {
+  return VORONOI3D_BOX_SIDE_X * VORONOI3D_BOX_SIDE_Y * VORONOI3D_BOX_SIDE_Z;
+}
+
+__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;
+}
+
+__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
@@ -62,6 +118,26 @@
  *  http://math.lbl.gov/voro++/
  ******************************************************************************/
 
+__attribute__((always_inline)) INLINE void voronoi_print_gnuplot_c(
+    struct voronoi_cell *c) {
+
+  int i, j, v;
+  double *x = c->x;
+
+  fprintf(stderr, "%g\t%g\t%g\n\n", x[0], x[1], x[2]);
+
+  for (i = 0; i < c->nvert; i++) {
+    for (j = 0; j < c->orders[i]; j++) {
+      v = c->edges[c->offsets[i] + j];
+      fprintf(stderr, "%g\t%g\t%g\n", c->vertices[3 * i + 0] + x[0],
+              c->vertices[3 * i + 1] + x[1], c->vertices[3 * i + 2] + x[2]);
+      fprintf(stderr, "%g\t%g\t%g\n\n", c->vertices[3 * v + 0] + x[0],
+              c->vertices[3 * v + 1] + x[1], c->vertices[3 * v + 2] + x[2]);
+    }
+  }
+  fprintf(stderr, "\n");
+}
+
 /**
  * @brief Print the contents of a 3D Voronoi cell
  *
@@ -611,6 +687,7 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
   int result = voronoi_intersect_find_closest_vertex(
       c, dx, r2, &u, &up, &us, &uw, &l, &lp, &ls, &lw, &q, &qp, &qs, &qw);
   if (result < 0) {
+    voronoi_print_gnuplot_c(c);
     error("Error while searching intersected edge!");
   }
   if (!result) {
@@ -625,7 +702,7 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
 
   int vindex = -1;
   int visitflags[VORONOI3D_MAXNUMVERT];
-  int dstack[VORONOI3D_MAXNUMVERT];
+  int dstack[2 * VORONOI3D_MAXNUMVERT];
   int dstack_size = 1;
   float r = 0.0f;
   int cs = -1, rp = -1;
@@ -639,12 +716,19 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
 
   if (complicated) {
 
+    // the search routine detected a vertex very close to the plane
+    // the index of this vertex is stored in up
+    // we proceed by checking the edges of this vertex
+
     lp = voronoi_get_edge(c, up, 0);
     lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
                              &teststack_size);
 
+    // the first edge can be below, above or on the plane
     if (lw != -1) {
 
+      // above or on the plane: we try to find one below the plane
+
       rp = lw;
       i = 1;
       lp = voronoi_get_edge(c, up, i);
@@ -653,11 +737,17 @@ __attribute__((always_inline)) INLINE void voronoi_intersect(
       while (lw != -1) {
         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...
+          voronoi_print_gnuplot_c(c);
           error(
               "Cell completely gone! This should not happen. (i == "
               "c->order[up], i = %d, c->orders[up] = %d, up = %d)\n"
-              "dx: [%g %g %g]",
-              i, c->orders[up], up, dx[0], dx[1], dx[2]);
+              "dx: [%g %g %g]\nv[up]: [%g %g %g]\nx: [%g %g %g]",
+              i, c->orders[up], up, dx[0], dx[1], dx[2], c->vertices[3 * up],
+              c->vertices[3 * up + 1], c->vertices[3 * up + 2], c->x[0],
+              c->x[1], c->x[2]);
         }
         lp = voronoi_get_edge(c, up, i);
         lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
diff --git a/tests/test125cells.c b/tests/test125cells.c
index d24d4d52e544dac103012b0d4fa91837d0c6290c..d9c04fdedec395ad03314d2732b49633663b1465 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -197,7 +197,7 @@ void reset_particles(struct cell *c, enum velocity_field vel,
     set_energy_state(p, press, size, density);
 
 #if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
-    float volume = p->mass / density;
+    float volume = p->conserved.mass / density;
 #if defined(GIZMO_SPH)
     p->geometry.volume = volume;
 #else
@@ -263,7 +263,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
         part->x[2] = offset[2] + size * (z + 0.5) / (float)n;
         part->h = size * h / (float)n;
 
-#ifdef GIZMO_SPH
+#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
         part->conserved.mass = density * volume / count;
 #else
         part->mass = density * volume / count;
@@ -279,7 +279,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
         hydro_first_init_part(part, xpart);
 
 #if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
-        float volume = part->mass / density;
+        float volume = part->conserved.mass / density;
 #ifdef GIZMO_SPH
         part->geometry.volume = volume;
 #else
diff --git a/tests/testVoronoi3D.c b/tests/testVoronoi3D.c
index 54aa5087327851671ca643978b6226f749ee0daa..b289d22e5668d130ce8f960f612ed305af24ac60 100644
--- a/tests/testVoronoi3D.c
+++ b/tests/testVoronoi3D.c
@@ -69,10 +69,12 @@ void test_calculate_cell() {
   voronoi_calculate_faces(&cell);
 
   /* Update these values if you ever change to another large cube! */
-  assert(cell.volume == 27.0f);
-  assert(cell.centroid[0] = 0.5f);
-  assert(cell.centroid[1] = 0.5f);
-  assert(cell.centroid[2] = 0.5f);
+  assert(cell.volume == voronoi_get_box_volume());
+  float box_centroid[3];
+  voronoi_get_box_centroid(box_centroid);
+  assert(cell.centroid[0] = box_centroid[0]);
+  assert(cell.centroid[1] = box_centroid[1]);
+  assert(cell.centroid[2] = box_centroid[2]);
 
   /* Check cell neighbours. */
   assert(cell.nface == 6);
@@ -84,35 +86,42 @@ void test_calculate_cell() {
   assert(cell.ngbs[5] == VORONOI3D_BOX_RIGHT);
 
   /* Check cell faces */
-  assert(cell.face_areas[0] == 9.0f);
-  assert(cell.face_midpoints[0][0] == 0.5f);
-  assert(cell.face_midpoints[0][1] == -1.0f);
-  assert(cell.face_midpoints[0][2] == 0.5f);
+  float face_midpoint[3], face_area;
+  face_area = voronoi_get_box_face(VORONOI3D_BOX_FRONT, face_midpoint);
+  assert(cell.face_areas[0] == face_area);
+  assert(cell.face_midpoints[0][0] == face_midpoint[0]);
+  assert(cell.face_midpoints[0][1] == face_midpoint[1]);
+  assert(cell.face_midpoints[0][2] == face_midpoint[2]);
 
-  assert(cell.face_areas[1] == 9.0f);
-  assert(cell.face_midpoints[1][0] == -1.0f);
-  assert(cell.face_midpoints[1][1] == 0.5f);
-  assert(cell.face_midpoints[1][2] == 0.5f);
+  face_area = voronoi_get_box_face(VORONOI3D_BOX_LEFT, face_midpoint);
+  assert(cell.face_areas[1] == face_area);
+  assert(cell.face_midpoints[1][0] == face_midpoint[0]);
+  assert(cell.face_midpoints[1][1] == face_midpoint[1]);
+  assert(cell.face_midpoints[1][2] == face_midpoint[2]);
 
-  assert(cell.face_areas[2] == 9.0f);
-  assert(cell.face_midpoints[2][0] == 0.5f);
-  assert(cell.face_midpoints[2][1] == 0.5f);
-  assert(cell.face_midpoints[2][2] == -1.0f);
+  face_area = voronoi_get_box_face(VORONOI3D_BOX_BOTTOM, face_midpoint);
+  assert(cell.face_areas[2] == face_area);
+  assert(cell.face_midpoints[2][0] == face_midpoint[0]);
+  assert(cell.face_midpoints[2][1] == face_midpoint[1]);
+  assert(cell.face_midpoints[2][2] == face_midpoint[2]);
 
-  assert(cell.face_areas[3] == 9.0f);
-  assert(cell.face_midpoints[3][0] == 0.5f);
-  assert(cell.face_midpoints[3][1] == 0.5f);
-  assert(cell.face_midpoints[3][2] == 2.0f);
+  face_area = voronoi_get_box_face(VORONOI3D_BOX_TOP, face_midpoint);
+  assert(cell.face_areas[3] == face_area);
+  assert(cell.face_midpoints[3][0] == face_midpoint[0]);
+  assert(cell.face_midpoints[3][1] == face_midpoint[1]);
+  assert(cell.face_midpoints[3][2] == face_midpoint[2]);
 
-  assert(cell.face_areas[4] == 9.0f);
-  assert(cell.face_midpoints[4][0] == 0.5f);
-  assert(cell.face_midpoints[4][1] == 2.0f);
-  assert(cell.face_midpoints[4][2] == 0.5f);
+  face_area = voronoi_get_box_face(VORONOI3D_BOX_BACK, face_midpoint);
+  assert(cell.face_areas[4] == face_area);
+  assert(cell.face_midpoints[4][0] == face_midpoint[0]);
+  assert(cell.face_midpoints[4][1] == face_midpoint[1]);
+  assert(cell.face_midpoints[4][2] == face_midpoint[2]);
 
-  assert(cell.face_areas[5] == 9.0f);
-  assert(cell.face_midpoints[5][0] == 2.0f);
-  assert(cell.face_midpoints[5][1] == 0.5f);
-  assert(cell.face_midpoints[5][2] == 0.5f);
+  face_area = voronoi_get_box_face(VORONOI3D_BOX_RIGHT, face_midpoint);
+  assert(cell.face_areas[5] == face_area);
+  assert(cell.face_midpoints[5][0] == face_midpoint[0]);
+  assert(cell.face_midpoints[5][1] == face_midpoint[1]);
+  assert(cell.face_midpoints[5][2] == face_midpoint[2]);
 }
 
 void test_paths() {
@@ -1104,6 +1113,7 @@ void test_paths() {
   }
 }
 
+#ifdef SHADOWSWIFT
 void set_coordinates(struct part *p, double x, double y, double z,
                      unsigned int id) {
   p->x[0] = x;
@@ -1112,8 +1122,10 @@ void set_coordinates(struct part *p, double x, double y, double z,
   p->id = id;
   voronoi_cell_init(&p->cell, p->x);
 }
+#endif
 
 void test_degeneracies() {
+#ifdef SHADOWSWIFT
   int idx = 0;
   /* make a small cube */
   struct part particles[100];
@@ -1177,6 +1189,7 @@ void test_degeneracies() {
     dx[2] = particles[0].x[2] - particles[i].x[2];
     voronoi_cell_interact(&particles[0].cell, dx, particles[i].id);
   }
+#endif
 }
 
 int main() {
@@ -1217,7 +1230,7 @@ int main() {
   /* Finalize cell and check results */
   voronoi_cell_finalize(&cell);
 
-  if (cell.volume != 0.125f) {
+  if (fabs(cell.volume - 0.125f) > 1.e-5) {
     error("Wrong volume: %g!", cell.volume);
   }
   if (fabs(cell.centroid[0] - 0.5f) > 1.e-5f ||
@@ -1231,7 +1244,7 @@ int main() {
   float A, midpoint[3];
   A = voronoi_get_face(&cell, 1, midpoint);
   if (A) {
-    if (A != 0.25f) {
+    if (fabs(A - 0.25f) > 1.e-5) {
       error("Wrong surface area: %g!", A);
     }
     if (fabs(midpoint[0] - 0.25f) > 1.e-5 || fabs(midpoint[1] - 0.5f) > 1.e-5 ||