First attempt at integration

src/Makefile.am

@@ -35,13 +35,13 @@ endif
 # List required headers
 include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
-    common_io.h single_io.h multipole.h map.h tools.h
+    common_io.h single_io.h multipole.h map.h tools.h poisson_disc.h
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
     units.c common_io.c single_io.c multipole.c version.c map.c \
-    kernel.c tools.c
+    kernel.c tools.c poisson_disc.c
 
 # Include files for distribution, not installation.
 noinst_HEADERS = atomic.h cycle.h error.h inline.h kernel.h vector.h \

src/engine.c

@@ -48,6 +48,7 @@
 #include "debug.h"
 #include "error.h"
 #include "timers.h"
+#include "poisson_disc.h"
 
 #ifdef LEGACY_GADGET2_SPH
 #include "runner_iact_legacy.h"
@@ -2071,24 +2072,27 @@ void engine_makeproxies(struct engine *e) {
 
 void engine_split(struct engine *e, int *grid) {
 
-  int j, k;
-  int ind[3];
+  //int j, k;
+  //int ind[3];
   struct space *s = e->s;
-  struct cell *c;
+  //struct cell *c;
 
   /* If we've got the wrong number of nodes, fail. */
-  if (e->nr_nodes != grid[0] * grid[1] * grid[2])
-    error("Grid size does not match number of nodes.");
+  //if (e->nr_nodes != grid[0] * grid[1] * grid[2])
+  //   error("Grid size does not match number of nodes.");
 
   /* Run through the cells and set their nodeID. */
   // message("s->dim = [%e,%e,%e]", s->dim[0], s->dim[1], s->dim[2]);
-  for (k = 0; k < s->nr_cells; k++) {
-    c = &s->cells[k];
-    for (j = 0; j < 3; j++) ind[j] = c->loc[j] / s->dim[j] * grid[j];
-    c->nodeID = ind[0] + grid[0] * (ind[1] + grid[1] * ind[2]);
-    // message("cell at [%e,%e,%e]: ind = [%i,%i,%i], nodeID = %i", c->loc[0],
-    // c->loc[1], c->loc[2], ind[0], ind[1], ind[2], c->nodeID);
-  }
+  //for (k = 0; k < s->nr_cells; k++) {
+  //  c = &s->cells[k];
+  //  for (j = 0; j < 3; j++) ind[j] = c->loc[j] / s->dim[j] * grid[j];
+  //  c->nodeID = ind[0] + grid[0] * (ind[1] + grid[1] * ind[2]);
+  //  // message("cell at [%e,%e,%e]: ind = [%i,%i,%i], nodeID = %i", c->loc[0],
+  //  // c->loc[1], c->loc[2], ind[0], ind[1], ind[2], c->nodeID);
+  //}
+
+  if (! poisson_split(s, e->nr_nodes) )
+    error("Failed to partition cells");
 
   /* Make the proxies. */
   engine_makeproxies(e);

src/poisson_disc.c

@@ -42,6 +42,7 @@
 
 /* Local headers. */
 #include "error.h"
+#include "poisson_disc.h"
 
 /* Useful defines. */
 #define MAX(a, b) ((a) > (b) ? (a) : (b))
@@ -211,7 +212,8 @@ void poisson_free(struct grid_data *grid) {
  * @param k number of attempts to pick an unmarked position (30 is a good
  * choice).
  */
-void poisson_disc(struct grid_data *grid, int dims[3], float radius, int k) {
+static void poisson_disc(struct grid_data *grid, int dims[3], float radius, 
+                         int k) {
 
   grid->radius = radius;
   grid->cell_size = radius / sqrtf(3.0f);
@@ -293,53 +295,54 @@ void poisson_disc(struct grid_data *grid, int dims[3], float radius, int k) {
   }
 }
 
-int main(int argc, char *argv[]) {
-  int dims[3];
-  float radius = 0.0f;
-  int k = 30;
-  struct grid_data grid;
-
-  /* Expected sample size. */
-  int N = 90;
-
-  /* Dimensions. */
-  dims[0] = 300;
-  dims[1] = 300;
-  dims[2] = 300;
+/**
+ * @brief Apply poisson disc partitioning to a cell structure.
+ *
+ * Generates a poisson disc sample for the 3D space of the cells and assigns
+ * each cells nodeID to the nearest sample point index, thus partitioning 
+ * the space into contiguous regions.
+ *
+ * @param s the space containing the cells to split into partitions.
+ * @param nparts number of parititions to generate.
+ * @result true if the partitioning succeeded.
+ */
+int poisson_split(struct space *s, int nparts) {
 
+  /* Randomize randoms. */
   srand(time(NULL));
 
   /*  Pick radius for expected sample size. */
-  radius = pow((dims[0] * dims[1] * dims[2]) / (N), 0.3333);
-  printf("# Radius = %f\n", radius);
+  float radius = pow((s->cdim[0]*s->cdim[1]*s->cdim[2])/(1.5*nparts), 0.3333);
+  message("# Radius = %f\n", radius);
 
-  /*  Sample is stocastic, so we may need to ask more than one to get the
+  /*  Sample is stocastic, so we may need to ask more than once to get the
    *  number of samples we require as a minimum. */
+  struct grid_data grid;
   grid.samplelen = 0;
-  while (grid.samplelen < N) {
-    printf("# Sampling...\n");
-    poisson_disc(&grid, dims, radius, k);
-    printf("# Samples = %d\n", grid.samplelen);
+  while (grid.samplelen < nparts) {
+    message("# Sampling...\n");
+    poisson_disc(&grid, s->cdim, radius, 30);
+    message("# Samples = %d\n", grid.samplelen);
   }
 
   for (int i = 0; i < grid.samplelen; i++) {
-    printf("# %f %f %f\n", grid.samplelist[i].x[0], grid.samplelist[i].x[1],
+    message("# %f %f %f\n", grid.samplelist[i].x[0], grid.samplelist[i].x[1],
            grid.samplelist[i].x[2]);
   }
 
   /*  Partition the space. Slow .... */
-
-  unsigned long int counts[N];
-  for (int i = 0; i < N; i++) {
+  unsigned long int counts[nparts];
+  for (int i = 0; i < nparts; i++) {
     counts[i] = 0;
   }
 
-  for (int i = 0; i < dims[0]; i++) {
-    for (int j = 0; j < dims[1]; j++) {
-      for (int k = 0; k < dims[2]; k++) {
+  int n = 0;
+  for (int i = 0; i < s->cdim[0]; i++) {
+    for (int j = 0; j < s->cdim[1]; j++) {
+      for (int k = 0; k < s->cdim[2]; k++) {
         int select = -1;
         float rsqmax = FLT_MAX;
-        for (int l = 0; l < N; l++) {
+        for (int l = 0; l < nparts; l++) {
           float dx = grid.samplelist[l].x[0] - (i + 0.5);
           float dy = grid.samplelist[l].x[1] - (j + 0.5);
           float dz = grid.samplelist[l].x[2] - (k + 0.5);
@@ -349,21 +352,22 @@ int main(int argc, char *argv[]) {
             select = l;
           }
         }
+        s->cells[n].nodeID = select;
         counts[select]++;
-        printf("%f %f %f %d\n", i + 0.5, j + 0.5, k + 0.5, select);
+        message("%f %f %f %d\n", i + 0.5, j + 0.5, k + 0.5, select);
       }
     }
   }
 
-  printf("# Counts:\n");
+  message("# Counts:\n");
   unsigned long int total = 0;
-  for (int i = 0; i < N; i++) {
-    printf("#  %d %ld\n", i, counts[i]);
+  for (int i = 0; i < nparts; i++) {
+    message("#  %d %ld\n", i, counts[i]);
     total += counts[i];
   }
-  printf("# total = %ld\n", total);
+  message("# total = %ld\n", total);
 
   poisson_free( &grid );
 
-  return 0;
+  return 1;
 }

src/poisson_disc.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ *
+ * 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_POISSON_DISC_H
+#define SWIFT_POISSON_DISC_H
+
+#include "space.h"
+int poisson_split(struct space *s, int nparts);
+
+#endif /* SWIFT_POISSON_DISC_H */