testVoronoi2D.c 6.57 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/*******************************************************************************
 * 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 <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
19
20
21
#include "../config.h"

/* Local headers. */
22
#include "hydro/Shadowswift/voronoi2d_algorithm.h"
23
#include "tools.h"
24

25
26
27
/* Number of cells used to test the 2D interaction algorithm */
#define TESTVORONOI2D_NUMCELL 100

28
int main(int argc, char *argv[]) {
29

30
  /* initialize simulation box */
31
32
  double anchor[3] = {-0.5f, -0.5f, -0.5f};
  double side[3] = {2.0f, 2.0f, 2.0f};
33

34
35
  /* test initialization and finalization algorithms */
  {
36
37
    message("Testing initialization and finalization algorithm...");

38
39
40
    struct voronoi_cell cell;
    double x[3] = {0.5, 0.5, 0.5};

41
    voronoi_cell_init(&cell, x, anchor, side);
42
43

    float maxradius = voronoi_cell_finalize(&cell);
44

45
    assert(maxradius == 2.0f * sqrtf(2.0f));
46

47
    assert(cell.volume == 4.0f);
48

49
50
    assert(cell.centroid[0] == 0.5f);
    assert(cell.centroid[1] == 0.5f);
51
52

    message("Done.");
53
  }
54

55
  /* test interaction algorithm: normal case */
56
  {
57
58
59
    message("Testing %i cell grid with random positions...",
            TESTVORONOI2D_NUMCELL);

60
61
62
63
64
65
66
    /* 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;
67

68
    for (i = 0; i < TESTVORONOI2D_NUMCELL; ++i) {
69
70
71
      x[0] = random_uniform(0., 1.);
      x[1] = random_uniform(0., 1.);
      voronoi_cell_init(&cells[i], x, anchor, side);
72
73
74
75
#ifdef VORONOI_VERBOSE
      message("cell[%i]: %g %g", i, x[0], x[1]);
#endif
    }
76

77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
    /* 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);
        }
      }
    }
108

109
    Atot = 0.0f;
110
111
    /* print the cells to the stdout */
    for (i = 0; i < TESTVORONOI2D_NUMCELL; ++i) {
112
#ifdef VORONOI_VERBOSE
113
      printf("Cell %i:\n", i);
114
      voronoi_print_cell(&cells[i]);
115
#endif
116
117
118
      voronoi_cell_finalize(&cells[i]);
      Atot += cells[i].volume;
    }
119

120
    /* Check the total surface area */
121
    assert(fabs(Atot - 1.0f) < 1.e-6);
122
123
124
125
126
127
128
129
130
131
132
133
134
135

    /* 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.");
136
  }
137

138
139
  /* test interaction algorithm */
  {
140
141
    message("Testing 100 cell grid with Cartesian mesh positions...");

142
143
144
145
146
147
148
149
150
151
152
    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;
153
        voronoi_cell_init(&cells[10 * i + j], x, anchor, side);
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
      }
    }

    /* 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);
        }
      }
    }

189
    Atot = 0.0f;
190
191
    /* print the cells to the stdout */
    for (i = 0; i < 100; ++i) {
192
#ifdef VORONOI_VERBOSE
193
      printf("Cell %i:\n", i);
194
      voronoi_print_cell(&cells[i]);
195
#endif
196
197
198
199
      voronoi_cell_finalize(&cells[i]);
      Atot += cells[i].volume;
    }

200
    /* Check the total surface area */
201
    assert(fabs(Atot - 1.0f) < 1.e-6);
202

203
204
205
206
207
    /* 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) {
Peter W. Draper's avatar
Peter W. Draper committed
208
209
210
211
212
      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);
213
214

    } else {
Peter W. Draper's avatar
Peter W. Draper committed
215
216
217
218
219
      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);
220
    }
221
    message("Done.");
222
223
  }

224
225
  return 0;
}