diff --git a/src/shadowfax/cell_shadowfax.h b/src/shadowfax/cell_shadowfax.h index 7d56e72e5db478165e2f5b120f756f79e7a9fb1e..a95fda6e2ad03bbc0c8f82aeb25d45d6c7db0dc8 100644 --- a/src/shadowfax/cell_shadowfax.h +++ b/src/shadowfax/cell_shadowfax.h @@ -45,12 +45,6 @@ __attribute__((always_inline)) INLINE static void cell_shadowfax_do_self1( 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++) { diff --git a/src/shadowfax/delaunay.h b/src/shadowfax/delaunay.h index 2cc11becc7c9796c32ac7eed9f12644d9ebacf41..89060af37fcfce796137c74827400ccf2b3c158a 100644 --- a/src/shadowfax/delaunay.h +++ b/src/shadowfax/delaunay.h @@ -1582,6 +1582,7 @@ inline static void delaunay_add_new_vertex(struct delaunay* restrict d, d->ngb_indices = (int*)swift_realloc( "delaunay ngb indices", d->ngb_indices, d->ngb_size * sizeof(int)); } + delaunay_assert(d->ngb_index == v - d->ngb_offset); d->ngb_cells[d->ngb_index] = cell_index; d->ngb_indices[d->ngb_index] = index_in_cell; ++d->ngb_index; diff --git a/src/shadowfax/voronoi.h b/src/shadowfax/voronoi.h index 6f24712c03ee9990eddb953ccce91a539ef67d1b..2ff017e6abb7d627a4740fe8776160a0bdc20d75 100644 --- a/src/shadowfax/voronoi.h +++ b/src/shadowfax/voronoi.h @@ -32,6 +32,10 @@ #include <string.h> +/*! @brief Store the edges of faces (so that the actual Voronoi grid can be + * reconstructed). */ +#define VORONOI_STORE_CONNECTIONS + /** * @brief Voronoi interface. * @@ -344,8 +348,9 @@ static inline void voronoi_init(struct voronoi *restrict v, /* Get a triangle containing this generator and the index of the generator within that triangle */ - int t0 = d->vertex_triangles[i + d->vertex_start]; - int vi0 = d->vertex_triangle_index[i + d->vertex_start]; + int dverti = i + d->vertex_start; + int t0 = d->vertex_triangles[dverti]; + int vi0 = d->vertex_triangle_index[dverti]; /* Add the first vertex for this cell: the circumcircle midpoint of this triangle */ v->vertex_number[i] = 1; @@ -399,17 +404,17 @@ static inline void voronoi_init(struct voronoi *restrict v, /* the neighbour corresponding to the face is the same vertex that determines the next triangle */ - int other_vertex = d->triangles[t0].vertices[vi1p2]; + int other_vertex = d->triangles[t1].vertices[vi1p2]; if (other_vertex < d->ngb_offset) { - if (other_vertex > i) { + if (other_vertex > dverti) { /* only store pairs once */ - voronoi_add_pair(v, 0, i, other_vertex, bx, by, cx, cy); + voronoi_add_pair(v, 0, dverti, other_vertex, bx, by, cx, cy); } } else { /* no check on other_vertex > i required, since this is always true */ int sid = d->ngb_cells[other_vertex - d->ngb_offset]; int pid = d->ngb_indices[other_vertex - d->ngb_offset]; - voronoi_add_pair(v, sid, i, pid, bx, by, cx, cy); + voronoi_add_pair(v, sid, dverti, pid, bx, by, cx, cy); } vi1 = d->triangles[t1].index_in_neighbour[vi1p2]; @@ -435,9 +440,9 @@ static inline void voronoi_init(struct voronoi *restrict v, v->face_midpoints[2 * c0 + 1] = face_midpoint[1]; if (first_vertex < d->ngb_offset) { - if (first_vertex > i) { + if (first_vertex > dverti) { /* only store pairs once */ - voronoi_add_pair(v, 0, i, first_vertex, bx, by, cx, cy); + voronoi_add_pair(v, 0, dverti, first_vertex, bx, by, cx, cy); } } else { /* no check on other_vertex > i required, since this is always true */ @@ -484,6 +489,29 @@ static inline void voronoi_check_grid(const struct voronoi *restrict v) { printf("Total volume: %g\n", V); } +static inline void voronoi_write_grid_new(const struct voronoi *restrict v, + const char *file_name) { + + 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]); + } + for (int i = 0; i < v->pair_index[0]; ++i) { + struct voronoi_pair *pair = &v->pairs[0][i]; + fprintf(file, "F\t%g\t%g\t%g\t%g\t%i\t%i\n", pair->a[0], pair->a[1], + pair->b[0], pair->b[1], pair->left, pair->right); + } + for (int ngb = 1; ngb < 27; ++ngb) { + for (int i = 0; i < v->pair_index[ngb]; ++i) { + struct voronoi_pair *pair = &v->pairs[ngb][i]; + fprintf(file, "F\t%g\t%g\t%g\t%g\t%i\t%i\n", pair->a[0], pair->a[1], + pair->b[0], pair->b[1], pair->left, -1); + } + } + fclose(file); +} + static inline void voronoi_write_grid(const struct voronoi *restrict v, FILE *file, size_t *offset) {