voronoi3d_algorithm.h

/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

#ifndef SWIFT_VORONOIXD_ALGORITHM_H
#define SWIFT_VORONOIXD_ALGORITHM_H

#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "error.h"
#include "inline.h"
#include "voronoi3d_cell.h"

/* For debugging purposes */
//#define LOOP_CHECK 1000

#ifdef LOOP_CHECK
/* We need to do the trickery below to get a unique counter for each call to the
   macro. This only works if the macro is never called twice on the same line.
 */
#define MERGE(a, b) a##b
#define LOOPCOUNTER_NAME(line) MERGE(loopcount, line)

/**
 * @brief Increase the given counter variable and check if it is still valid.
 *
 * @param counter Counter to increase.
 * @param line_number Line number where the while is called.
 * @return 1 if the counter is still valid, 0 otherwise.
 */
__attribute__((always_inline)) INLINE int check_counter(int *counter,
                                                        int line_number) {
  ++(*counter);
  if ((*counter) == LOOP_CHECK) {
    error("Number of iterations reached maximum (=%i) in while on line %i!",
          LOOP_CHECK, line_number);
  }
  return 1;
}

/* safewhile is a wrapper around a while that adds a unique counter variable to
   the loop that is increased by 1 for each time the loop is executed, and
   causes the code to crash if this number exceeds a given value.
   We use this to quickly enable or disable number of iterations checks for a
   large number of while loops */
#define safewhile(condition)          \
  int LOOPCOUNTER_NAME(__LINE__) = 0; \
  while (check_counter(&LOOPCOUNTER_NAME(__LINE__), __LINE__) && (condition))

#else /* LOOP_CHECK */

/* If LOOP_CHECK is not defined, safewhile and while are EXACTLY the same */
#define safewhile(condition) while (condition)

#endif /* LOOP_CHECK */

/* This flag activates a number of expensive geometrical checks that help
   finding bugs. */
//#define VORONOI3D_EXPENSIVE_CHECKS

/* Tolerance parameter used to decide when to use more precise geometric
   criteria */
#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
   we make the strong assumption that there will never be a particle that has
   one of these as particle ID. */
#define VORONOI3D_BOX_FRONT 18446744073709551600llu
#define VORONOI3D_BOX_BACK 18446744073709551601llu
#define VORONOI3D_BOX_TOP 18446744073709551602llu
#define VORONOI3D_BOX_BOTTOM 18446744073709551603llu
#define VORONOI3D_BOX_LEFT 18446744073709551604llu
#define VORONOI3D_BOX_RIGHT 18446744073709551605llu

/*******************************************************************************
 * 3D specific methods
 *
 * Most of these methods are based on the source code of voro++:
 *  http://math.lbl.gov/voro++/
 ******************************************************************************/

/**
 * @brief Print the given cell to the stderr in a format that can be easily
 * plotted using gnuplot.
 *
 * This method prints to the stderr instead of stdout to make it possible to use
 * it right before crashing the code.
 *
 * @param c Voronoi cell to print.
 */
__attribute__((always_inline)) INLINE void voronoi_print_gnuplot_c(
    const struct voronoi_cell *c) {

  int i, j, v;
  const 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];
      if (v < 0) {
        v = -v - 1;
      }
      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
 *
 * @param cell 3D Voronoi cell
 */
__attribute__((always_inline)) INLINE void voronoi_print_cell(
    const struct voronoi_cell *cell) {

  int i, j;

  fprintf(stderr, "x: %g %g %g\n", cell->x[0], cell->x[1], cell->x[2]);
  fprintf(stderr, "nvert: %i\n", cell->nvert);

  for (i = 0; i < cell->nvert; ++i) {
    fprintf(stderr, "%i: %g %g %g (%i)\n", i, cell->vertices[3 * i],
            cell->vertices[3 * i + 1], cell->vertices[3 * i + 2],
            cell->orders[i]);
    for (j = 0; j < cell->orders[i]; ++j) {
      fprintf(stderr, "%i (%i)\n", cell->edges[cell->offsets[i] + j],
              cell->edgeindices[cell->offsets[i] + j]);
    }
  }
  fprintf(stderr, "\n");
}

/**
 * @brief Get the index of the vertex pointed to by the given edge of the given
 * vertex.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of the cell.
 * @param edge Edge of that vertex.
 * @return Index of the vertex on the other side of the edge.
 */
__attribute__((always_inline)) INLINE int voronoi_get_edge(
    const struct voronoi_cell *c, int vertex, int edge) {
  return c->edges[c->offsets[vertex] + edge];
}

/**
 * @brief Get the index of the given edge in the edge list of the vertex on the
 * other side of the edge of the given vertex.
 *
 * Suppose that the given vertex has edges [edge1, edge2, given_edge], and that
 * the vertex on the other side of given_edge has edges [edge1, given_edge,
 * edge2], then this method returns 1.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of the cell.
 * @param edge Edge of that vertex.
 * @return Index of that edge in the edge list of the vertex on the other side
 * of the edge.
 */
__attribute__((always_inline)) INLINE int voronoi_get_edgeindex(
    const struct voronoi_cell *c, int vertex, int edge) {
  return c->edgeindices[c->offsets[vertex] + edge];
}

/**
 * @brief Set the index of the vertex on the other side of the given edge of the
 * given vertex.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of the cell.
 * @param edge Edge of that vertex.
 * @param value Index of the vertex on the other side of that edge.
 */
__attribute__((always_inline)) INLINE void voronoi_set_edge(
    struct voronoi_cell *c, int vertex, int edge, int value) {
  c->edges[c->offsets[vertex] + edge] = value;
}

/**
 * @brief Set the index of the given edge in the edge list of the vertex on the
 * other side of the edge of the given vertex.
 *
 * Suppose that the given vertex has edges [edge1, edge2, given_edge], and we
 * want to tell this method that the vertex on the other side of given_edge has
 * edges [edge1, given_edge, edge2], then we need to pass on a value of 1 to
 * this method.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of that cell.
 * @param edge Edge of that vertex.
 * @param value Index of that edge in the edge list of the vertex on the other
 * side of the edge.
 */
__attribute__((always_inline)) INLINE void voronoi_set_edgeindex(
    struct voronoi_cell *c, int vertex, int edge, int value) {
  c->edgeindices[c->offsets[vertex] + edge] = value;
}

/**
 * @brief Get the neighbour for the given edge of the given vertex.
 *
 * An edge is shared by two faces, and each face has a neighbour. Luckily, each
 * edge also has two endpoints, so we can get away with storing only one
 * neighbour per endpoint of an edge. We have complete freedom in choosing which
 * neighbour to store in which endpoint, but we need to be consistent about it.
 * Here we use the following convention: if we take a vector pointing away from
 * the given vertex along the given edge direction, then we store the neighbour
 * that corresponds to the face to the right if looking to the cell from the
 * outside. This is the face that you encounter first when rotating counter-
 * clockwise around that vector, starting from outside the cell.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of that cell.
 * @param edge Edge of that vertex.
 * @return Index of the neighbour corresponding to that edge and vertex.
 */
__attribute__((always_inline)) INLINE int voronoi_get_ngb(
    const struct voronoi_cell *c, int vertex, int edge) {
  return c->ngbs[c->offsets[vertex] + edge];
}

/**
 * @brief Set the neighbour for the given edge of the given vertex.
 *
 * An edge is shared by two faces, and each face has a neighbour. Luckily, each
 * edge also has two endpoints, so we can get away with storing only one
 * neighbour per endpoint of an edge. We have complete freedom in choosing which
 * neighbour to store in which endpoint, but we need to be consistent about it.
 * Here we use the following convention: if we take a vector pointing away from
 * the given vertex along the given edge direction, then we store the neighbour
 * that corresponds to the face to the right if looking to the cell from the
 * outside. This is the face that you encounter first when rotating counter-
 * clockwise around that vector, starting from outside the cell.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of that cell.
 * @param edge Edge of that vertex.
 * @param value Index of the neighbour corresponding to that edge and vertex.
 */
__attribute__((always_inline)) INLINE void voronoi_set_ngb(
    struct voronoi_cell *c, int vertex, int edge, int value) {
  c->ngbs[c->offsets[vertex] + edge] = value;
}

/**
 * @brief Check if the 3D Voronoi cell is still consistent.
 *
 * A cell is consistent if its edges are consistent, i.e. if edge e of vertex v1
 * points to vertex v2, then v2 should have an edge that points to v1 as well,
 * and then the edge index of vertex v1 should contain the index of that edge
 * in the edge list of v2. We also check if all vertices have orders of at least
 * 3, and if all vertices are actually part of the vertex list.
 * Oh, and we check if the cell actually has vertices.
 *
 * @param cell 3D Voronoi cell to check
 */
__attribute__((always_inline)) INLINE void voronoi_check_cell_consistency(
    const struct voronoi_cell *c) {

  int i, j, e, l, m;

  if (c->nvert < 4) {
    error("Found cell with only %i vertices!", c->nvert);
  }

  for (i = 0; i < c->nvert; ++i) {
    if (c->orders[i] < 3) {
      voronoi_print_cell(c);
      error("Found cell with vertex of order %i!", c->orders[i]);
    }
    for (j = 0; j < c->orders[i]; ++j) {
      e = voronoi_get_edge(c, i, j);
      if (e >= c->nvert) {
        voronoi_print_cell(c);
        error("Found cell with edges that lead to non-existing vertices!");
      }
      if (e < 0) {
        continue;
      }
      l = voronoi_get_edgeindex(c, i, j);
      m = voronoi_get_edge(c, e, l);
      if (m != i) {
        /* voronoi_print_gnuplot_c(c); */
        voronoi_print_cell(c);
        fprintf(stderr, "i: %i, j: %i, e: %i, l: %i, m: %i\n", i, j, e, l, m);
        error("Cell inconsistency!");
      }
    }
  }
}

/**
 * @brief Check if the given vertex is above, below or on the cutting plane
 * defined by the given parameters.
 *
 * @param v Coordinates of a cell vertex, relative w.r.t. the position of the
 * generator of the cell.
 * @param dx Half of the relative distance vector between the position of the
 * generator of the cell and the position of the neighbouring cell that
 * generates the cutting plane, pointing from the generator position to the
 * cutting plane.
 * @param r2 Squared length of dx.
 * @param test Variable to store the result of the geometric test in, which
 * corresponds to the projected distance between the generator and the vertex
 * along dx.
 * @param teststack Stack to store the results of the N last tests in (for
 * debugging purposes only).
 * @param teststack_size Next available field in the teststack, is reset to 0 if
 * the teststack is full (so the N+1th results is overwritten; for debugging
 * purposes only).
 * @return Result of the test: -1 if the vertex is below the cutting plane, +1
 * if it is above, and 0 if it is on the cutting plane.
 */
__attribute__((always_inline)) INLINE int voronoi_test_vertex(
    const float *v, const float *dx, float r2, float *test, float *teststack,
    int *teststack_size) {

  *test = v[0] * dx[0] + v[1] * dx[1] + v[2] * dx[2] - r2;

  teststack[*teststack_size] = *test;
  *teststack_size = *teststack_size + 1;
  if (*teststack_size == 2 * VORONOI3D_MAXNUMVERT) {
    *teststack_size = 0;
  }

  if (*test < -VORONOI3D_TOLERANCE) {
    return -1;
  }
  if (*test > VORONOI3D_TOLERANCE) {
    return 1;
  }
  return 0;
}

/**
 * @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, const double *anchor, const double *side) {

  cell->nvert = 8;

  /* (0, 0, 0) -- 0 */
  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] = 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] = 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] = 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] = 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] = 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] = 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] = 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;
  cell->orders[2] = 3;
  cell->orders[3] = 3;
  cell->orders[4] = 3;
  cell->orders[5] = 3;
  cell->orders[6] = 3;
  cell->orders[7] = 3;

  /* edges are ordered counterclockwise w.r.t. a vector pointing from the
     cell generator to the vertex
     (0, 0, 0) corner */
  cell->offsets[0] = 0;
  cell->edges[0] = 1;
  cell->edges[1] = 2;
  cell->edges[2] = 4;
  cell->edgeindices[0] = 0;
  cell->edgeindices[1] = 2;
  cell->edgeindices[2] = 0;

  /* (0, 0, 1) corner */
  cell->offsets[1] = 3;
  cell->edges[3] = 0;
  cell->edges[4] = 5;
  cell->edges[5] = 3;
  cell->edgeindices[3] = 0;
  cell->edgeindices[4] = 2;
  cell->edgeindices[5] = 1;

  /* (0, 1, 0) corner */
  cell->offsets[2] = 6;
  cell->edges[6] = 3;
  cell->edges[7] = 6;
  cell->edges[8] = 0;
  cell->edgeindices[6] = 0;
  cell->edgeindices[7] = 0;
  cell->edgeindices[8] = 1;

  /* (0, 1, 1) corner */
  cell->offsets[3] = 9;
  cell->edges[9] = 2;
  cell->edges[10] = 1;
  cell->edges[11] = 7;
  cell->edgeindices[9] = 0;
  cell->edgeindices[10] = 2;
  cell->edgeindices[11] = 0;

  /* (1, 0, 0) corner */
  cell->offsets[4] = 12;
  cell->edges[12] = 0;
  cell->edges[13] = 6;
  cell->edges[14] = 5;
  cell->edgeindices[12] = 2;
  cell->edgeindices[13] = 2;
  cell->edgeindices[14] = 0;

  /* (1, 0, 1) corner */
  cell->offsets[5] = 15;
  cell->edges[15] = 4;
  cell->edges[16] = 7;
  cell->edges[17] = 1;
  cell->edgeindices[15] = 2;
  cell->edgeindices[16] = 1;
  cell->edgeindices[17] = 1;

  /* (1, 1, 0) corner */
  cell->offsets[6] = 18;
  cell->edges[18] = 2;
  cell->edges[19] = 7;
  cell->edges[20] = 4;
  cell->edgeindices[18] = 1;
  cell->edgeindices[19] = 2;
  cell->edgeindices[20] = 1;

  /* (1, 1, 1) corner */
  cell->offsets[7] = 21;
  cell->edges[21] = 3;
  cell->edges[22] = 5;
  cell->edges[23] = 6;
  cell->edgeindices[21] = 2;
  cell->edgeindices[22] = 1;
  cell->edgeindices[23] = 1;

  /* ngbs[3*i+j] is the neighbour corresponding to the plane clockwise of
     edge j of vertex i (when going from edge j to vertex i)
     we set them to a ridiculously large value to be able to track faces without
     neighbour */
  cell->ngbs[0] = VORONOI3D_BOX_FRONT;  /* (000) - (001) */
  cell->ngbs[1] = VORONOI3D_BOX_LEFT;   /* (000) - (010) */
  cell->ngbs[2] = VORONOI3D_BOX_BOTTOM; /* (000) - (100) */

  cell->ngbs[3] = VORONOI3D_BOX_LEFT;  /* (001) - (000) */
  cell->ngbs[4] = VORONOI3D_BOX_FRONT; /* (001) - (101) */
  cell->ngbs[5] = VORONOI3D_BOX_TOP;   /* (001) - (011) */

  cell->ngbs[6] = VORONOI3D_BOX_LEFT;   /* (010) - (011) */
  cell->ngbs[7] = VORONOI3D_BOX_BACK;   /* (010) - (110) */
  cell->ngbs[8] = VORONOI3D_BOX_BOTTOM; /* (010) - (000) */

  cell->ngbs[9] = VORONOI3D_BOX_BACK;  /* (011) - (010) */
  cell->ngbs[10] = VORONOI3D_BOX_LEFT; /* (011) - (001) */
  cell->ngbs[11] = VORONOI3D_BOX_TOP;  /* (011) - (111) */

  cell->ngbs[12] = VORONOI3D_BOX_FRONT;  /* (100) - (000) */
  cell->ngbs[13] = VORONOI3D_BOX_BOTTOM; /* (100) - (110) */
  cell->ngbs[14] = VORONOI3D_BOX_RIGHT;  /* (100) - (101) */

  cell->ngbs[15] = VORONOI3D_BOX_FRONT; /* (101) - (100) */
  cell->ngbs[16] = VORONOI3D_BOX_RIGHT; /* (101) - (111) */
  cell->ngbs[17] = VORONOI3D_BOX_TOP;   /* (101) - (001) */

  cell->ngbs[18] = VORONOI3D_BOX_BOTTOM; /* (110) - (010) */
  cell->ngbs[19] = VORONOI3D_BOX_BACK;   /* (110) - (111) */
  cell->ngbs[20] = VORONOI3D_BOX_RIGHT;  /* (110) - (100) */

  cell->ngbs[21] = VORONOI3D_BOX_BACK;  /* (111) - (011) */
  cell->ngbs[22] = VORONOI3D_BOX_TOP;   /* (111) - (101) */
  cell->ngbs[23] = VORONOI3D_BOX_RIGHT; /* (111) - (110) */
}

/**
 * @brief Find an edge of the voronoi_cell that intersects the cutting plane.
 *
 * There is a large number of possible paths through this method, each of which
 * is covered by a separate unit test in testVoronoi3D. Paths have been numbered
 * in the inline comments to help identify them.
 *
 * @param c 3D Voronoi cell.
 * @param dx Vector pointing from pj to the midpoint of the line segment between
 * pi and pj.
 * @param r2 Squared length of dx.
 * @param u Projected distance between the plane and the closest vertex above
 * the plane, along dx.
 * @param up Index of the closest vertex above the plane.
 * @param us Index of the edge of vertex up that intersects the plane.
 * @param uw Result of the last test_vertex call for vertex up.
 * @param l Projected distance between the plane and the closest vertex below
 * the plane, along dx.
 * @param lp Index of the closest vertex below the plane.
 * @param ls Index of the edge of vertex lp that intersects the plane.
 * @param lw Result of the last test_vertex call for vertex lp.
 * @param q Projected distance between the plane and a test vertex, along dx.
 * @param qp Index of the test vertex.
 * @param qs Index of the edge of the test vertex that is connected to up.
 * @param qw Result of the last test_vertex call involving qp.
 * @return A negative value if an error occurred, 0 if the plane does not
 * intersect the cell, 1 if nothing special happened and 2 if we have a
 * complicated setup.
 */
__attribute__((always_inline)) INLINE int voronoi_intersect_find_closest_vertex(
    struct voronoi_cell *c, const float *dx, float r2, float *u, int *up,
    int *us, int *uw, float *l, int *lp, int *ls, int *lw, float *q, int *qp,
    int *qs, int *qw) {

  /* stack to store all vertices that have already been tested (debugging
     only) */
  float teststack[2 * VORONOI3D_MAXNUMVERT];
  /* size of the used part of the stack */
  int teststack_size = 0;
  /* flag signalling a complicated setup */
  int complicated;

  /* test the first vertex: uw = -1 if it is below the plane, 1 if it is above
     0 if it is very close to the plane, and things become complicated... */
  *uw = voronoi_test_vertex(&c->vertices[0], dx, r2, u, teststack,
                            &teststack_size);
  *up = 0;
  complicated = 0;
  if ((*uw) == 0) {

    /* PATH 0 */
    complicated = 1;

  } else {

    /* two options: either the vertex is above or below the plane */

    if ((*uw) == 1) {

      /* PATH 1 */

      /* above: try to find a vertex below
         we test all edges of the current vertex stored in up (vertex 0) until
         we either find one below the plane or closer to the plane */
      *lp = voronoi_get_edge(c, (*up), 0);
      *lw = voronoi_test_vertex(&c->vertices[3 * (*lp)], dx, r2, l, teststack,
                                &teststack_size);
      *us = 1;
      /* Not in while: PATH 1.0 */
      /* somewhere in while: PATH 1.1 */
      /* last valid option of while: PATH 1.2 */
      safewhile((*us) < c->orders[(*up)] && (*l) >= (*u)) {
        *lp = voronoi_get_edge(c, (*up), (*us));
        *lw = voronoi_test_vertex(&c->vertices[3 * (*lp)], dx, r2, l, teststack,
                                  &teststack_size);
        ++(*us);
      }
      /* we increased us too much, correct this */
      --(*us);
      if ((*l) >= (*u)) {
        /* PATH 1.3 */
        /* up is the closest vertex to the plane, but is above the plane
           since the entire cell is convex, up is the closest vertex of all
           vertices of the cell
           this means the entire cell is supposedly above the plane, which is
           impossible */
        message(
            "Cell completely gone! This should not happen. (l >= u, l = %g, u "
            "= %g)",
            (*l), (*u));
        return -1;
      }
      /* we know that lp is closer to the plane or below the plane
         now find the index of the edge up-lp in the edge list of lp */
      *ls = voronoi_get_edgeindex(c, (*up), (*us));

      /* if lp is also above the plane, replace up by lp and repeat the process
         until lp is below the plane */
      safewhile((*lw) == 1) {
        /* PATH 1.4 */
        *u = (*l);
        *up = (*lp);
        *us = 0;
        /* no while: PATH 1.4.0 */
        /* somewhere in while: PATH 1.4.1 */
        /* last valid option of while: PATH 1.4.2 */
        safewhile((*us) < (*ls) && (*l) >= (*u)) {
          *lp = voronoi_get_edge(c, (*up), (*us));
          *lw = voronoi_test_vertex(&c->vertices[3 * (*lp)], dx, r2, l,
                                    teststack, &teststack_size);
          ++(*us);
        }
        if ((*l) >= (*u)) {
          ++(*us);
          /* no while: PATH 1.4.3 */
          /* somewhere in while: PATH 1.4.4 */
          /* last valid option of while: PATH 1.4.5 */
          safewhile((*us) < c->orders[(*up)] && (*l) >= (*u)) {
            *lp = voronoi_get_edge(c, (*up), (*us));
            *lw = voronoi_test_vertex(&c->vertices[3 * (*lp)], dx, r2, l,
                                      teststack, &teststack_size);
            ++(*us);
          }
          if ((*l) >= (*u)) {
            /* PATH 1.4.6 */
            message(
                "Cell completely gone! This should not happen. (l >= u, l = "
                "%g, u = %g)",
                (*l), (*u));
            return -1;
          }
        }
        --(*us);
        *ls = voronoi_get_edgeindex(c, (*up), (*us));
      }
      /* if lp is too close to the plane, replace up by lp and proceed to
         complicated setup */
      if ((*lw) == 0) {
        /* PATH 1.5 */
        *up = (*lp);
        complicated = 1;
      }
    } else { /* if(uw == 1) */

      /* PATH 2 */

      /* below: try to find a vertex above
         we test all edges of the current vertex stored in up (vertex 0) until
         we either find one above the plane or closer to the plane */

      *qp = voronoi_get_edge(c, (*up), 0);
      *qw = voronoi_test_vertex(&c->vertices[3 * (*qp)], dx, r2, q, teststack,
                                &teststack_size);
      *us = 1;
      /* not in while: PATH 2.0 */
      /* somewhere in while: PATH 2.1 */
      /* last valid option of while: PATH 2.2 */
      safewhile((*us) < c->orders[(*up)] && (*u) >= (*q)) {
        *qp = voronoi_get_edge(c, (*up), (*us));
        *qw = voronoi_test_vertex(&c->vertices[3 * (*qp)], dx, r2, q, teststack,
                                  &teststack_size);
        ++(*us);
      }
      if ((*u) >= (*q)) {
        /* PATH 2.3 */
        /* up is the closest vertex to the plane and is below the plane
           since the cell is convex, up is the closest vertex of all vertices of
           the cell
           this means that the entire cell is below the plane
           The cell is unaltered. */
        return 0;
      } else {
        /* the last increase in the loop pushed us too far, correct this */
        --(*us);
      }

      /* repeat the above process until qp is closer or above the plane */
      safewhile((*qw) == -1) {
        /* PATH 2.4 */
        *qs = voronoi_get_edgeindex(c, (*up), (*us));
        *u = (*q);
        *up = (*qp);
        *us = 0;
        /* no while: PATH 2.4.0 */
        /* somewhere in while: PATH 2.4.1 */
        /* last valid option of while: 2.4.2 */
        safewhile((*us) < (*qs) && (*u) >= (*q)) {
          *qp = voronoi_get_edge(c, (*up), (*us));
          *qw = voronoi_test_vertex(&c->vertices[3 * (*qp)], dx, r2, q,
                                    teststack, &teststack_size);
          ++(*us);
        }
        if ((*u) >= (*q)) {
          ++(*us);
          /* no while: PATH 2.4.3 */
          /* somewhere in while: PATH 2.4.4 */
          /* last valid option of while: PATH 2.4.5 */
          safewhile((*us) < c->orders[(*up)] && (*u) >= (*q)) {
            *qp = voronoi_get_edge(c, (*up), (*us));
            *qw = voronoi_test_vertex(&c->vertices[3 * (*qp)], dx, r2, q,
                                      teststack, &teststack_size);
            ++(*us);
          }
          if ((*u) >= (*q)) {
            /* PATH 2.4.6 */
            /* cell unaltered */
            return 0;
          }
        }
        --(*us);
      }
      if ((*qw) == 1) {
        /* qp is above the plane: initialize lp to up and replace up by qp */
        *lp = (*up);
        *ls = (*us);
        *l = (*u);
        *up = (*qp);
        *us = voronoi_get_edgeindex(c, (*lp), (*ls));
        *u = (*q);
      } else {
        /* PATH 2.5 */
        /* too close to call: go to complicated setup */
        *up = (*qp);
        complicated = 1;
      }

    } /* if(uw == 1) */

  } /* if(uw == 0) */

  if (complicated) {
    return 2;
  } else {
    return 1;
  }
}

/**
 * @brief Intersect the given cell with the midplane between the cell generator
 * and a neighbouring cell at the given relative position and with the given ID.
 *
 * This method is the core of the Voronoi algorithm. If anything goes wrong
 * geometrically, it most likely goes wrong somewhere within this method.
 *
 * @param c 3D Voronoi cell.
 * @param odx The original relative distance vector between the cell generator
 * and the intersecting neighbour, as it is passed on to runner_iact_density
 * (remember: odx = pi->x - pj->x).
 * @param ngb ID of the intersecting neighbour (pj->id in runner_iact_density).
 */
__attribute__((always_inline)) INLINE void voronoi_intersect(
    struct voronoi_cell *c, const float *odx, unsigned long long ngb) {

  /* vector pointing from pi to the midpoint of the line segment between pi and
     pj. This corresponds to -0.5*odx */
  float dx[3];
  /* squared norm of dx */
  float r2;
  /* u: distance between the plane and the closest vertex above the plane (up)
     l: distance between the plane and the closest vertex below the plane (lp)
     q: distance between the plane and the vertex that is currently being
     tested (qp) */
  float u = 0.0f, l = 0.0f, q = 0.0f;
  /* up: index of the closest vertex above the plane
     us: index of the edge of vertex up that intersects the plane
     uw: result of the last orientation test involving vertex u
     same naming used for vertex l and vertex q */
  int up = -1, us = -1, uw = -1, lp = -1, ls = -1, lw = -1, qp = -1, qs = -1,
      qw = -1;
  /* auxiliary flag used to capture degeneracies */
  int complicated = -1;

  /* stack to store all vertices that have already been tested (debugging
     only) */
  float teststack[2 * VORONOI3D_MAXNUMVERT];
  /* size of the used part of the stack */
  int teststack_size = 0;

#ifdef VORONOI3D_EXPENSIVE_CHECKS
  voronoi_check_cell_consistency(c);
#endif

  /* initialize dx and r2 */
  dx[0] = -0.5f * odx[0];
  dx[1] = -0.5f * odx[1];
  dx[2] = -0.5f * odx[2];
  r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];

  /* find an intersected edge of the cell */
  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) {
    /* the closest_vertex test only found vertices above the intersecting plane
       this would mean that the entire cell lies above the midplane of the line
       segment connecting a point inside the cell (the generator) and a point
       that could be inside or outside the cell (the neighbour). This is
       geometrically absurd and should NEVER happen. */
    voronoi_print_gnuplot_c(c);
    error("Error while searching intersected edge!");
  }
  if (result == 0) {
    /* no intersection */
    return;
  }
  if (result == 2) {
    complicated = 1;
  } else {
    complicated = 0;
  }

  /* At this point:
      up contains a vertex above the plane
      lp contains a vertex below the plane
      us and ls contain the index of the edge that connects up and lp, this edge
      is intersected by the midplane
      u and l contain the projected distances of up and lp to the midplane,
      along dx
     IF complicated is 1, up contains a vertex that is considered to be on the
     plane. All other variables can be considered to be uninitialized in this
     case. */

  int vindex = -1;
  int visitflags[VORONOI3D_MAXNUMVERT];
  int dstack[2 * VORONOI3D_MAXNUMVERT];
  int dstack_size = 0;
  float r = 0.0f;
  int cs = -1, rp = -1;
  int double_edge = 0;
  int i = -1, j = -1, k = -1;

  /* initialize visitflags */
  for (i = 0; i < VORONOI3D_MAXNUMVERT; ++i) {
    visitflags[i] = 0;
  }

  if (complicated) {

    /* We've entered the complicated setup, which means that somewhere along the
       way we found a vertex that is on or very close to the midplane. The index
       of that vertex is stored in up, all other variables are meaningless at
       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). 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 make sure that up contains the index of a vertex extending beyond the
       plane on exit. */
    dstack[dstack_size] = up;
    ++dstack_size;
    lw = 0;
    j = 0;
    safewhile(j < dstack_size && lw != -1) {
      up = dstack[j];
      for (i = 0; i < c->orders[up]; ++i) {
        lp = voronoi_get_edge(c, up, i);
        lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
                                 &teststack_size);
        if (lw == -1) {
          /* jump out of the for loop */
          break;
        }
        if (lw == 0) {
          /* only add each vertex to the stack once */
          k = 0;
          safewhile(k < dstack_size && dstack[k] != lp) { ++k; }
          if (k == dstack_size) {
            dstack[dstack_size] = lp;
            ++dstack_size;
          }
        }
      }
      ++j;
    }

    /* we increased j after lw was calculated, so only the value of lw should be
       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 midplane!");
    }
    /* reset the delete stack, we need it later on */
    dstack_size = 0;

    /* 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 */

    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);
      lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
                               &teststack_size);
      safewhile(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...
             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 == "
              "c->order[up], i = %d, c->orders[up] = %d, up = %d)\n"
              "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,
                                 &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] && lw == -1) {
        lp = voronoi_get_edge(c, up, j);
        lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
                                 &teststack_size);
        ++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 */