1/*
2 * CellList_test.hpp
3 *
4 * Created on: Mar 23, 2015
5 * Author: Pietro Incardona
6 */
7
8#include "CellList.hpp"
9#include "CellListM.hpp"
10#include "Grid/grid_sm.hpp"
11
12#ifndef CELLLIST_TEST_HPP_
13#define CELLLIST_TEST_HPP_
14
15
16/*! \brief Test cell structure
17 *
18 * \tparam CellS
19 *
20 */
21template<unsigned int dim, typename T, typename CellS> void Test_cell_s(SpaceBox<dim,T> & box)
22{
23 //! [Declare a cell list]
24 //Space where is living the Cell list
25 //SpaceBox<dim,T> box({0.0f,0.0f,0.0f},{1.0f,1.0f,1.0f});
26
27 // Subdivisions
28 size_t div[dim] = {16,16,16};
29
30 // Origin
31 Point<dim,T> org({0.0,0.0,0.0});
32
33 // id Cell list
34 CellS cl2(box,div);
35 //! [Declare a cell list]
36
37 // grid info
38 grid_sm<dim,void> g_info(div);
39
40 // Test force reallocation in case of Cell list fast
41 for (size_t i = 0 ; i < CELL_REALLOC * 3 ; i++)
42 {
43 cl2.add(org,i);
44 }
45
46 // Check the elements
47 BOOST_REQUIRE_EQUAL(cl2.getNelements(cl2.getCell(org)),CELL_REALLOC * 3ul);
48 for (size_t i = 0 ; i < CELL_REALLOC * 3 ; i++)
49 {
50 BOOST_REQUIRE_EQUAL(cl2.get(cl2.getCell(org),i),i);
51 }
52
53 //! [Usage of cell list]
54
55 // id Cell list
56 CellS cl1(box,div);
57
58 // Create a grid iterator
59 grid_key_dx_iterator<dim> g_it(g_info);
60
61 // Iterate through each element
62 // Add 1 element for each cell
63
64 // Usefull definition of points
65 Point<dim,T> end = box.getP2() - box.getP1();
66 Point<dim,T> middle = end / div / 2.0;
67 Point<dim,T> spacing = end / div;
68
69 Point<dim,T> offset[dim] = {middle,middle,middle};
70
71 // Create offset shift vectors
72 for (size_t i = 0 ; i < dim ; i++)
73 {
74 offset[i].get(i) += (1.0 / div[i]) / 8.0;
75 }
76
77 openfpm::vector<Point<dim,T>> pos;
78 size_t id = 0;
79
80 while (g_it.isNext())
81 {
82 // Add 2 particles on each cell
83
84 Point<dim,T> key = Point<dim,T>(g_it.get().toPoint());
85 key = pmul(key,spacing) + offset[0] + box.getP1();
86 pos.add(key);
87
88 cl1.add(key,id);
89 ++id;
90
91 key = Point<dim,T>(g_it.get().toPoint());
92 key = pmul(key,spacing) + offset[1] + box.getP1();
93 pos.add(key);
94
95 cl1.add(key,id);
96 ++id;
97
98 ++g_it;
99 }
100
101 //! [Usage of cell list]
102
103 // check the cell are correctly filled
104
105 // reset iterator
106 g_it.reset();
107
108 while (g_it.isNext())
109 {
110 // Check that there are 2 particles on each cell
111
112 Point<dim,T> key = Point<dim,T>(g_it.get().toPoint());
113 key = pmul(key,spacing) + offset[2] + box.getP1();
114
115 size_t cell = cl1.getCell(key);
116 size_t n_ele = cl1.getNelements(cell);
117
118 BOOST_REQUIRE_EQUAL(n_ele,2ul);
119 BOOST_REQUIRE_EQUAL((long int)(cl1.get(cell,1) - cl1.get(cell,0)),1);
120
121 ++g_it;
122 }
123
124 // reset itarator
125 g_it.reset();
126
127 //! [remove one particle from each cell]
128
129 while (g_it.isNext())
130 {
131 // remove 1 particle on each cell
132
133 Point<dim,T> key = Point<dim,T>(g_it.get().toPoint());
134 key = pmul(key,spacing) + offset[0] + box.getP1();
135
136 auto cell = cl1.getCell(key);
137
138 // Remove the first particle in the cell
139 cl1.remove(cell,0);
140 ++g_it;
141 }
142
143 //! [remove one particle from each cell]
144
145 // Check we have 1 object per cell
146 g_it.reset();
147
148 while (g_it.isNext())
149 {
150 // remove 1 particle on each cell
151
152 Point<dim,T> key = Point<dim,T>(g_it.get().toPoint());
153 key = pmul(key,spacing) + offset[0] + box.getP1();
154
155 auto cell = cl1.getCell(key);
156 size_t n_ele = cl1.getNelements(cell);
157
158 BOOST_REQUIRE_EQUAL(n_ele,1ul);
159 ++g_it;
160 }
161
162
163 // Check we have 1 object per cell
164
165 // Create a grid iterator
166 grid_key_dx<dim> p1(1,1,1);
167 grid_key_dx<dim> p2(div[0]-2,div[1]-2,div[2]-2);
168 grid_key_dx_iterator_sub<dim> g_it_s(g_info,p1,p2);
169 id = 0;
170
171 while (g_it_s.isNext())
172 {
173 //! [Usage of the neighborhood iterator]
174
175 Point<dim,T> key = Point<dim,T>(g_it_s.get().toPoint());
176 key = pmul(key,spacing) + offset[0] + box.getP1();
177
178 auto NN = cl1.template getNNIterator<NO_CHECK>(cl1.getCell(key));
179 size_t total = 0;
180
181 while(NN.isNext())
182 {
183 // total
184
185 total++;
186
187 ++NN;
188 }
189
190 //! [Usage of the neighborhood iterator]
191
192 BOOST_REQUIRE_EQUAL(total,(size_t)openfpm::math::pow(3,dim));
193
194 // in SE1_CLASS the cell list consider this construction as an hack
195 // disable the test
196
197#ifndef SE_CLASS1
198
199 id = cl1.get(cl1.getCell(key),0);
200 auto NNSym = cl1.template getNNIteratorSym<NO_CHECK>(cl1.getCell(key),id,pos);
201 total = 0;
202
203 while(NNSym.isNext())
204 {
205 // total
206
207 total++;
208
209 ++NNSym;
210 }
211
212 BOOST_REQUIRE_EQUAL(total,(size_t)openfpm::math::pow(3,dim) / 2 + 1);
213
214#endif
215
216 ++g_it_s;
217 }
218
219}
220
221
222/*! \brief Test cell structure
223 *
224 * \tparam CellS
225 *
226 */
227template<unsigned int dim, typename T, typename CellS> void Test_cell_sM(SpaceBox<dim,T> & box)
228{
229 //Space where is living the Cell list
230 //SpaceBox<dim,T> box({0.0f,0.0f,0.0f},{1.0f,1.0f,1.0f});
231
232 // Subdivisions
233 size_t div[dim] = {16,16,16};
234
235 // Origin
236 Point<dim,T> org({0.0,0.0,0.0});
237
238 // grid info
239 grid_sm<dim,void> g_info(div);
240
241 //! [Usage of cell list multi]
242
243 // CellS = CellListM<dim,T,8>
244 CellS cl1(box,div);
245
246 // Create a grid iterator
247 grid_key_dx_iterator<dim> g_it(g_info);
248
249 // Iterate through each element
250 // Add 1 element for each cell
251
252 // Usefull definition of points
253 Point<dim,T> end = box.getP2() - box.getP1();
254 Point<dim,T> middle = end / div / 2.0;
255 Point<dim,T> spacing = end / div;
256
257 Point<dim,T> offset[dim] = {middle,middle,middle};
258
259 // Create offset shift vectors
260 for (size_t i = 0 ; i < dim ; i++)
261 {
262 offset[i].get(i) += (1.0 / div[i]) / 8.0;
263 }
264
265 openfpm::vector<Point<dim,T>> phase1;
266 openfpm::vector<Point<dim,T>> phase2;
267
268 openfpm::vector<pos_v<openfpm::vector<Point<dim,T>>>> phases;
269 phases.add(pos_v<openfpm::vector<Point<dim,T>>>(phase1));
270 phases.add(pos_v<openfpm::vector<Point<dim,T>>>(phase2));
271
272 size_t id = 0;
273
274 while (g_it.isNext())
275 {
276 // Add 2 particles on each cell
277
278 Point<dim,T> key = Point<dim,T>(g_it.get().toPoint());
279 key = pmul(key,spacing) + offset[0] + box.getP1();
280
281 phase1.add(key);
282
283 cl1.add(key,id,0);
284
285 key = Point<dim,T>(g_it.get().toPoint());
286 key = pmul(key,spacing) + offset[1] + box.getP1();
287
288 phase2.add(key);
289 cl1.add(key,id,1);
290 ++id;
291
292 ++g_it;
293 }
294
295 // check the cell are correctly filled
296
297 // reset iterator
298 g_it.reset();
299
300 while (g_it.isNext())
301 {
302 // Add 2 particles on each cell
303
304 Point<dim,T> key = Point<dim,T>(g_it.get().toPoint());
305 key = pmul(key,spacing) + offset[2] + box.getP1();
306
307 size_t cell = cl1.getCell(key);
308 size_t n_ele = cl1.getNelements(cell);
309
310 size_t p1 = cl1.getP(cell,1);
311 size_t p2 = cl1.getP(cell,0);
312
313 size_t v1 = cl1.getV(cell,1);
314 size_t v2 = cl1.getV(cell,0);
315
316 BOOST_REQUIRE_EQUAL(n_ele,2ul);
317 BOOST_REQUIRE_EQUAL((long int)(p1 - p2),0);
318 BOOST_REQUIRE_EQUAL((long int)(v1 - v2),1);
319 ++g_it;
320 }
321}
322
323
324
325template<typename CellList> void Test_CellDecomposer_consistent()
326{
327 Box<2,float> bx({-1.0/3.0,-1.0/3.0},{1.0/3.0,1.0/3.0});
328
329 size_t div[2] = {36,36};
330
331 CellDecomposer_sm<2,float,shift<2,float>> cd(bx,div,1);
332
333 Box<2,float> bx_sub({-1.0/5.0,-1.0/5.0},{1.0/5.0,1.0/5.0});
334
335 size_t bc[2] = {NON_PERIODIC,NON_PERIODIC};
336 CellList cl(cd,bx_sub);
337 Box<2,long int> bx_int = cd.convertDomainSpaceIntoGridUnits(bx_sub,bc);
338
339 BOOST_REQUIRE_EQUAL(bx_int.getLow(0),8);
340 BOOST_REQUIRE_EQUAL(bx_int.getLow(1),8);
341
342 BOOST_REQUIRE_EQUAL(bx_int.getHigh(0),28);
343 BOOST_REQUIRE_EQUAL(bx_int.getHigh(1),28);
344
345 cd.convertCellUnitsIntoDomainSpace(bx_sub);
346
347 BOOST_REQUIRE_EQUAL(bx_sub.getLow(0),-1.0f/5.0f);
348 BOOST_REQUIRE_EQUAL(bx_sub.getLow(1),-1.0f/5.0f);
349
350 BOOST_REQUIRE_EQUAL(bx_sub.getHigh(0),1.0f/5.0f);
351 BOOST_REQUIRE_EQUAL(bx_sub.getHigh(1),1.0f/5.0f);
352}
353
354/*! \brief Test cell structure
355 *
356 * \tparam CellS
357 *
358 */
359template<unsigned int dim, typename T, typename CellS> void Test_NN_iterator_radius(SpaceBox<dim,T> & box)
360{
361 //Space where is living the Cell list
362 //SpaceBox<dim,T> box({0.0f,0.0f,0.0f},{1.0f,1.0f,1.0f});
363
364 // Subdivisions
365 size_t div[dim];
366 size_t div2[dim];
367
368 for (size_t i = 0 ; i < dim ; i++)
369 {
370 div[i] = 17;
371 div2[i] = 34;
372 }
373
374 // grid info
375 grid_sm<dim,void> g_info(div);
376 grid_sm<dim,void> g_info2(div2);
377
378 //! [Usage of cell list multi]
379
380 // CellS = CellListM<dim,T,8>
381 CellS cl1(box,div,2);
382 CellS cl2(box,div2,3);
383
384 T radius = (box.getHigh(0) - box.getLow(0))/div[0];
385
386 cl2.setRadius( radius );
387
388 // create a vector of random point
389
390 openfpm::vector<Point<dim,float>> vrp;
391
392 for (size_t j = 0 ; j < 10000 ; j++)
393 {
394 vrp.add();
395
396 for (size_t i = 0 ; i < dim ; i++)
397 {
398 vrp.template get<0>(j)[i] = ((float)rand() / RAND_MAX)*(box.getHigh(i) - box.getLow(i)) + box.getLow(i);
399 }
400 }
401
402 auto g_it = vrp.getIterator();
403
404 while (g_it.isNext())
405 {
406 Point<dim,T> xp = vrp.get(g_it.get());
407
408 size_t debug = cl1.getCell(xp);
409
410 cl1.add(xp,g_it.get());
411 cl2.add(xp,g_it.get());
412
413 ++g_it;
414 }
415
416 // Get the neighborhood of each particle and compare the numbers of cell
417
418 bool match = true;
419 auto g_it2 = vrp.getIterator();
420
421 size_t number_of_nn = 0;
422 size_t number_of_nn2 = 0;
423
424 while (g_it2.isNext())
425 {
426 Point<dim,T> xp = vrp.get(g_it2.get());
427
428 openfpm::vector<size_t> ids1;
429 openfpm::vector<size_t> ids2;
430
431 size_t local1 = 0;
432 size_t local2 = 0;
433
434 auto NNit = cl1.getNNIterator(cl1.getCell(xp));
435
436 while (NNit.isNext())
437 {
438 auto q = NNit.get();
439
440 // calculate distance
441
442 Point<dim,T> xq = vrp.get(q);
443 Point<dim,T> r = xq - xp;
444
445 if (r.norm() <= radius)
446 {ids1.add(q);}
447
448 number_of_nn++;
449 local1++;
450
451 ++NNit;
452 }
453
454 auto NN2it = cl2.getNNIteratorRadius(cl2.getCell(xp));
455
456 while (NN2it.isNext())
457 {
458 auto q = NN2it.get();
459
460 Point<dim,T> xq = vrp.get(q);
461 Point<dim,T> r = xq - xp;
462
463 if (r.norm() <= radius)
464 {ids2.add(q);}
465
466 number_of_nn2++;
467 local2++;
468
469 ++NN2it;
470 }
471
472 // Sort ids1
473 ids1.sort();
474 ids2.sort();
475
476 match &= ids1.size() == ids2.size();
477
478 for (size_t i = 0 ; i < ids1.size() ; i++)
479 {match &= ids1.get(i) == ids2.get(i);}
480
481 if (match == false)
482 {break;}
483
484 ++g_it2;
485 }
486
487 BOOST_REQUIRE_EQUAL(match,true);
488 BOOST_REQUIRE(number_of_nn2 < number_of_nn);
489}
490
491BOOST_AUTO_TEST_SUITE( CellList_test )
492
493BOOST_AUTO_TEST_CASE ( NN_radius_check )
494{
495 SpaceBox<2,float> box1({0.1,0.1},{0.3,0.5});
496 SpaceBox<3,float> box2({0.0,0.1,0.2},{0.3,0.7,0.5});
497 SpaceBox<3,float> box3({0.0,0.0,0.0},{1.0,1.0,1.0});
498
499 std::cout << "Test cell list radius" << "\n";
500
501 Test_NN_iterator_radius<2,float,CellList<2,float,Mem_fast<>,shift<2,float>>>(box1);
502 Test_NN_iterator_radius<3,float,CellList<3,float,Mem_fast<>,shift<3,float>>>(box2);
503 Test_NN_iterator_radius<3,float,CellList<3,float,Mem_fast<>,shift<3,float>>>(box3);
504
505 std::cout << "End cell list" << "\n";
506}
507
508BOOST_AUTO_TEST_CASE( CellList_use)
509{
510 std::cout << "Test cell list" << "\n";
511
512 SpaceBox<3,double> box({0.0f,0.0f,0.0f},{1.0f,1.0f,1.0f});
513 SpaceBox<3,double> box2({-1.0f,-1.0f,-1.0f},{1.0f,1.0f,1.0f});
514 Test_cell_s<3,double,CellList<3,double,Mem_fast<>>>(box);
515 Test_cell_s<3,double,CellList<3,double,Mem_fast<>,shift<3,double>> >(box2);
516 Test_cell_sM<3,double,CellListM<3,double,8>>(box);
517 Test_cell_sM<3,double,CellListM<3,double,8>>(box2);
518
519
520 Test_cell_s<3,double,CellList<3,double,Mem_bal<>>>(box);
521 Test_cell_s<3,double,CellList<3,double,Mem_mw<>>>(box);
522
523 std::cout << "End cell list" << "\n";
524
525 // Test the cell list
526}
527
528BOOST_AUTO_TEST_CASE( CellList_consistent )
529{
530 Test_CellDecomposer_consistent<CellList<2,float,Mem_fast<>,shift<2,float>>>();
531}
532
533BOOST_AUTO_TEST_CASE( CellList_NNc_csr_calc )
534{
535 openfpm::vector<std::pair<grid_key_dx<3>,grid_key_dx<3>>> cNN;
536
537 NNcalc_csr(cNN);
538
539 BOOST_REQUIRE_EQUAL(cNN.size(),14ul);
540
541 BOOST_REQUIRE(cNN.get(0).first == grid_key_dx<3>(0,0,0));
542 BOOST_REQUIRE(cNN.get(0).second == grid_key_dx<3>(0,0,0));
543
544 BOOST_REQUIRE(cNN.get(1).first == grid_key_dx<3>(0,0,0));
545 BOOST_REQUIRE(cNN.get(1).second == grid_key_dx<3>(0,0,1));
546
547 BOOST_REQUIRE(cNN.get(2).first == grid_key_dx<3>(0,0,1));
548 BOOST_REQUIRE(cNN.get(2).second == grid_key_dx<3>(0,1,0));
549
550 BOOST_REQUIRE(cNN.get(3).first == grid_key_dx<3>(0,0,0));
551 BOOST_REQUIRE(cNN.get(3).second == grid_key_dx<3>(0,1,0));
552
553 BOOST_REQUIRE(cNN.get(4).first == grid_key_dx<3>(0,0,0));
554 BOOST_REQUIRE(cNN.get(4).second == grid_key_dx<3>(0,1,1));
555
556 BOOST_REQUIRE(cNN.get(5).first == grid_key_dx<3>(0,1,1));
557 BOOST_REQUIRE(cNN.get(5).second == grid_key_dx<3>(1,0,0));
558
559 BOOST_REQUIRE(cNN.get(6).first == grid_key_dx<3>(0,1,0));
560 BOOST_REQUIRE(cNN.get(6).second == grid_key_dx<3>(1,0,0));
561
562 BOOST_REQUIRE(cNN.get(7).first == grid_key_dx<3>(0,1,0));
563 BOOST_REQUIRE(cNN.get(7).second == grid_key_dx<3>(1,0,1));
564
565 BOOST_REQUIRE(cNN.get(8).first == grid_key_dx<3>(0,0,1));
566 BOOST_REQUIRE(cNN.get(8).second == grid_key_dx<3>(1,0,0));
567
568 BOOST_REQUIRE(cNN.get(9).first == grid_key_dx<3>(0,0,0));
569 BOOST_REQUIRE(cNN.get(9).second == grid_key_dx<3>(1,0,0));
570
571 BOOST_REQUIRE(cNN.get(10).first == grid_key_dx<3>(0,0,0));
572 BOOST_REQUIRE(cNN.get(10).second == grid_key_dx<3>(1,0,1));
573
574 BOOST_REQUIRE(cNN.get(11).first == grid_key_dx<3>(0,0,1));
575 BOOST_REQUIRE(cNN.get(11).second == grid_key_dx<3>(1,1,0));
576
577 BOOST_REQUIRE(cNN.get(12).first == grid_key_dx<3>(0,0,0));
578 BOOST_REQUIRE(cNN.get(12).second == grid_key_dx<3>(1,1,0));
579
580 BOOST_REQUIRE(cNN.get(13).first == grid_key_dx<3>(0,0,0));
581 BOOST_REQUIRE(cNN.get(13).second == grid_key_dx<3>(1,1,1));
582
583}
584
585
586
587BOOST_AUTO_TEST_SUITE_END()
588
589#endif /* CELLLIST_TEST_HPP_ */
590