/*******************************************************************************
* This file is part of SWIFT.
* Copyright (C) 2017 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 .
*
******************************************************************************/
#include
/* Local headers. */
#include "hydro/Shadowswift/voronoi2d_algorithm.h"
#include "tools.h"
/* Number of cells used to test the 2D interaction algorithm */
#define TESTVORONOI2D_NUMCELL 100
int main(int argc, char *argv[]) {
/* initialize simulation box */
double anchor[3] = {-0.5f, -0.5f, -0.5f};
double side[3] = {2.0f, 2.0f, 2.0f};
/* test initialization and finalization algorithms */
{
message("Testing initialization and finalization algorithm...");
struct voronoi_cell cell;
double x[3] = {0.5, 0.5, 0.5};
voronoi_cell_init(&cell, x, anchor, side);
float maxradius = voronoi_cell_finalize(&cell);
assert(maxradius == 2.0f * sqrtf(2.0f));
assert(cell.volume == 4.0f);
assert(cell.centroid[0] == 0.5f);
assert(cell.centroid[1] == 0.5f);
message("Done.");
}
/* test interaction algorithm: normal case */
{
message("Testing %i cell grid with random positions...",
TESTVORONOI2D_NUMCELL);
/* create 100 cells with random positions in [0,1]x[0,1] */
struct voronoi_cell cells[TESTVORONOI2D_NUMCELL];
double x[2];
float dx[2];
int i, j;
float Atot;
struct voronoi_cell *cell_i, *cell_j;
for (i = 0; i < TESTVORONOI2D_NUMCELL; ++i) {
x[0] = random_uniform(0., 1.);
x[1] = random_uniform(0., 1.);
voronoi_cell_init(&cells[i], x, anchor, side);
#ifdef VORONOI_VERBOSE
message("cell[%i]: %g %g", i, x[0], x[1]);
#endif
}
/* interact the cells (with periodic boundaries) */
for (i = 0; i < TESTVORONOI2D_NUMCELL; ++i) {
cell_i = &cells[i];
for (j = 0; j < TESTVORONOI2D_NUMCELL; ++j) {
if (i != j) {
cell_j = &cells[j];
dx[0] = cell_i->x[0] - cell_j->x[0];
dx[1] = cell_i->x[1] - cell_j->x[1];
/* periodic boundaries */
if (dx[0] >= 0.5f) {
dx[0] -= 1.0f;
}
if (dx[0] < -0.5f) {
dx[0] += 1.0f;
}
if (dx[1] >= 0.5f) {
dx[1] -= 1.0f;
}
if (dx[1] < -0.5f) {
dx[1] += 1.0f;
}
#ifdef VORONOI_VERBOSE
message("Cell %i before:", i);
voronoi_print_cell(&cells[i]);
message("Interacting cell %i with cell %i (%g %g, %g %g", i, j,
cells[i].x[0], cells[i].x[1], cells[j].x[0], cells[j].x[1]);
#endif
voronoi_cell_interact(cell_i, dx, j);
}
}
}
Atot = 0.0f;
/* print the cells to the stdout */
for (i = 0; i < TESTVORONOI2D_NUMCELL; ++i) {
#ifdef VORONOI_VERBOSE
printf("Cell %i:\n", i);
voronoi_print_cell(&cells[i]);
#endif
voronoi_cell_finalize(&cells[i]);
Atot += cells[i].volume;
}
/* Check the total surface area */
assert(fabs(Atot - 1.0f) < 1.e-6);
/* Check the neighbour relations for an arbitrary cell: cell 44
We plotted the grid and manually found the correct neighbours and their
order. */
assert(cells[44].nvert == 7);
assert(cells[44].ngbs[0] == 26);
assert(cells[44].ngbs[1] == 38);
assert(cells[44].ngbs[2] == 3);
assert(cells[44].ngbs[3] == 33);
assert(cells[44].ngbs[4] == 5);
assert(cells[44].ngbs[5] == 90);
assert(cells[44].ngbs[6] == 4);
message("Done.");
}
/* test interaction algorithm */
{
message("Testing 100 cell grid with Cartesian mesh positions...");
struct voronoi_cell cells[100];
double x[2];
float dx[2];
int i, j;
float Atot;
struct voronoi_cell *cell_i, *cell_j;
for (i = 0; i < 10; ++i) {
for (j = 0; j < 10; ++j) {
x[0] = (i + 0.5f) * 0.1;
x[1] = (j + 0.5f) * 0.1;
voronoi_cell_init(&cells[10 * i + j], x, anchor, side);
}
}
/* interact the cells (with periodic boundaries) */
for (i = 0; i < 100; ++i) {
cell_i = &cells[i];
for (j = 0; j < 100; ++j) {
if (i != j) {
cell_j = &cells[j];
dx[0] = cell_i->x[0] - cell_j->x[0];
dx[1] = cell_i->x[1] - cell_j->x[1];
/* periodic boundaries */
if (dx[0] >= 0.5f) {
dx[0] -= 1.0f;
}
if (dx[0] < -0.5f) {
dx[0] += 1.0f;
}
if (dx[1] >= 0.5f) {
dx[1] -= 1.0f;
}
if (dx[1] < -0.5f) {
dx[1] += 1.0f;
}
#ifdef VORONOI_VERBOSE
message("Cell %i before:", i);
voronoi_print_cell(&cells[i]);
message("Interacting cell %i with cell %i (%g %g, %g %g", i, j,
cells[i].x[0], cells[i].x[1], cells[j].x[0], cells[j].x[1]);
#endif
voronoi_cell_interact(cell_i, dx, j);
}
}
}
Atot = 0.0f;
/* print the cells to the stdout */
for (i = 0; i < 100; ++i) {
#ifdef VORONOI_VERBOSE
printf("Cell %i:\n", i);
voronoi_print_cell(&cells[i]);
#endif
voronoi_cell_finalize(&cells[i]);
Atot += cells[i].volume;
}
/* Check the total surface area */
assert(fabs(Atot - 1.0f) < 1.e-6);
/* Check the neighbour relations for an arbitrary cell: cell 44 We plotted
the grid and manually found the correct neighbours and their
order. Variation is found when optimizing, so we have two possible
outcomes... */
if (cells[44].nvert == 5) {
assert(cells[44].nvert == 5);
assert(cells[44].ngbs[0] == 43);
assert(cells[44].ngbs[1] == 34);
assert(cells[44].ngbs[2] == 45);
assert(cells[44].ngbs[3] == 55);
} else {
assert(cells[44].nvert == 4);
assert(cells[44].ngbs[0] == 34);
assert(cells[44].ngbs[1] == 45);
assert(cells[44].ngbs[2] == 54);
assert(cells[44].ngbs[3] == 43);
}
message("Done.");
}
return 0;
}