diff --git a/src/Makefile.am b/src/Makefile.am index 0fdc51571424e2a9ad91d1d7c714d91446702161..fae2e297888f9c5288e92ed01b2934ef2c6dccce 100644 --- a/src/Makefile.am +++ b/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 \ diff --git a/src/engine.c b/src/engine.c index 354fa42af9c955f09b64ea2bb6af84e010f3537b..6fa6a9158958bd1e75ed666bd752e72aa1051224 100644 --- a/src/engine.c +++ b/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); diff --git a/src/poisson_disc.c b/src/poisson_disc.c index 36070f5ffb5a3e5e08a851add544fb11d8ef7615..7baf84c1a2921b5274d63b0eaf8b5797f24bfd93 100644 --- a/src/poisson_disc.c +++ b/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; } diff --git a/src/poisson_disc.h b/src/poisson_disc.h new file mode 100644 index 0000000000000000000000000000000000000000..e7c87e2470e86cbd35b79ccf23f24ad04b190a3a --- /dev/null +++ b/src/poisson_disc.h @@ -0,0 +1,25 @@ +/******************************************************************************* + * 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 */