vector_dist_cell_list_tests.cpp 48.5 KB
Newer Older
1 2 3 4 5 6 7
/*
 * vector_dist_cell_list_tests.hpp
 *
 *  Created on: Aug 16, 2016
 *      Author: i-bird
 */

8
#include "config.h"
incardon's avatar
incardon committed
9 10 11 12 13
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include "Point_test.hpp"
#include "Vector/performance/vector_dist_performance_common.hpp"
#include "Vector/vector_dist.hpp"
14

incardon's avatar
incardon committed
15 16
extern void print_test_v(std::string test, size_t sz);
extern long int decrement(long int k, long int step);
17 18 19

///////////////////////// test hilb ///////////////////////////////

20
void test_reorder_sfc(reorder_opt opt)
21
{
22
	Vcluster<> & v_cl = create_vcluster();
23 24 25 26 27 28 29 30 31 32

	if (v_cl.getProcessingUnits() > 48)
		return;

    // set the seed
	// create the random generator engine
	std::srand(v_cl.getProcessUnitID());
    std::default_random_engine eg;
    std::uniform_real_distribution<float> ud(0.0f, 1.0f);

33 34 35
#ifdef TEST_COVERAGE_MODE
    long int k = 24288 * v_cl.getProcessingUnits();
#else
36
    long int k = 524288 * v_cl.getProcessingUnits();
37
#endif
38 39 40 41

	long int big_step = k / 4;
	big_step = (big_step == 0)?1:big_step;

42
	print_test_v( "Testing 2D vector with sfc curve reordering k<=",k);
43 44 45 46

	// 2D test
	for ( ; k >= 2 ; k-= decrement(k,big_step) )
	{
47
		BOOST_TEST_CHECKPOINT( "Testing 2D vector with sfc curve reordering k=" << k );
48 49 50 51 52 53

		Box<2,float> box({0.0,0.0},{1.0,1.0});

		// Boundary conditions
		size_t bc[2]={NON_PERIODIC,NON_PERIODIC};

incardon's avatar
incardon committed
54
		vector_dist<2,float, Point_test<float> > vd(k,box,bc,Ghost<2,float>(0.01));
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69

		auto it = vd.getIterator();

		while (it.isNext())
		{
			auto key = it.get();

			vd.getPos(key)[0] = ud(eg);
			vd.getPos(key)[1] = ud(eg);

			++it;
		}

		vd.map();

70 71
		// in case of SE_CLASS3 get a cell-list without ghost get is an error

72 73
		// Create first cell list

74
		auto NN1 = vd.getCellList(0.01,true);
75 76 77 78 79

		//An order of a curve
		int32_t m = 6;

		//Reorder a vector
80
		vd.reorder(m,opt);
81 82

		// Create second cell list
83
		auto NN2 = vd.getCellList(0.01,true);
84 85 86 87 88 89 90 91 92 93 94 95

		//Check equality of cell sizes
		for (size_t i = 0 ; i < NN1.getGrid().size() ; i++)
		{
			size_t n1 = NN1.getNelements(i);
			size_t n2 = NN2.getNelements(i);

			BOOST_REQUIRE_EQUAL(n1,n2);
		}
	}
}

incardon's avatar
incardon committed
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
///////////////////////// test hilb ///////////////////////////////

void test_reorder_cl()
{
	Vcluster<> & v_cl = create_vcluster();

	if (v_cl.getProcessingUnits() > 48)
		return;

    // set the seed
	// create the random generator engine
	std::srand(v_cl.getProcessUnitID());
    std::default_random_engine eg;
    std::uniform_real_distribution<float> ud(0.0f, 1.0f);

#ifdef TEST_COVERAGE_MODE
    long int k = 24288 * v_cl.getProcessingUnits();
#else
    long int k = 524288 * v_cl.getProcessingUnits();
#endif

	long int big_step = k / 4;
	big_step = (big_step == 0)?1:big_step;

	print_test_v( "Testing 2D vector with sfc curve reordering k<=",k);

	// 2D test
	for ( ; k >= 2 ; k-= decrement(k,big_step) )
	{
		BOOST_TEST_CHECKPOINT( "Testing 2D vector with sfc curve reordering k=" << k );

		Box<2,float> box({0.0,0.0},{1.0,1.0});

		// Boundary conditions
		size_t bc[2]={NON_PERIODIC,NON_PERIODIC};

		vector_dist<2,float, Point_test<float> > vd(k,box,bc,Ghost<2,float>(0.01));

		auto it = vd.getIterator();

		while (it.isNext())
		{
			auto key = it.get();

			vd.getPos(key)[0] = ud(eg);
			vd.getPos(key)[1] = ud(eg);

			++it;
		}

		vd.map();

		// in case of SE_CLASS3 get a cell-list without ghost get is an error

		// Create first cell list

		auto NN1 = vd.getCellList(0.01,true);

		//Reorder a vector
		vd.reorder_rcut(0.01);

		// Create second cell list
		auto NN2 = vd.getCellList(0.01,true);

		//Check equality of cell sizes
		for (size_t i = 0 ; i < NN1.getGrid().size() ; i++)
		{
			size_t n1 = NN1.getNelements(i);
			size_t n2 = NN2.getNelements(i);

			BOOST_REQUIRE_EQUAL(n1,n2);
		}
	}
}

incardon's avatar
incardon committed
171 172
BOOST_AUTO_TEST_SUITE( vector_dist_cell_list_test_suite )

173 174 175 176 177 178
BOOST_AUTO_TEST_CASE( vector_dist_reorder_2d_test )
{
	test_reorder_sfc(reorder_opt::HILBERT);
	test_reorder_sfc(reorder_opt::LINEAR);
}

incardon's avatar
incardon committed
179 180 181 182 183
BOOST_AUTO_TEST_CASE( vector_dist_reorder_cl_test )
{
	test_reorder_cl();
}

184 185
BOOST_AUTO_TEST_CASE( vector_dist_cl_random_vs_hilb_forces_test )
{
186
	Vcluster<> & v_cl = create_vcluster();
187 188

	if (v_cl.getProcessingUnits() > 48)
incardon's avatar
incardon committed
189
	{return;}
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216

	///////////////////// INPUT DATA //////////////////////

	// Dimensionality of the space
	const size_t dim = 3;
	// Cut-off radiuses. Can be put different number of values
	openfpm::vector<float> cl_r_cutoff {0.05};
	// The starting amount of particles (remember that this number is multiplied by number of processors you use for testing)
	size_t cl_k_start = 10000;
	// The lower threshold for number of particles
	size_t cl_k_min = 1000;
	// Ghost part of distributed vector
	double ghost_part = 0.05;

	///////////////////////////////////////////////////////

	//For different r_cut
	for (size_t r = 0; r < cl_r_cutoff.size(); r++ )
	{
		//Cut-off radius
		float r_cut = cl_r_cutoff.get(r);

		//Number of particles
		size_t k = cl_k_start * v_cl.getProcessingUnits();

		std::string str("Testing " + std::to_string(dim) + "D vector's forces (random vs hilb celllist) k<=");

incardon's avatar
incardon committed
217
		print_test_v(str,k);
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237

		//For different number of particles
		for (size_t k_int = k ; k_int >= cl_k_min ; k_int/=2 )
		{
			BOOST_TEST_CHECKPOINT( "Testing " << dim << "D vector's forces (random vs hilb celllist) k<=" << k_int );

			Box<dim,float> box;

			for (size_t i = 0; i < dim; i++)
			{
				box.setLow(i,0.0);
				box.setHigh(i,1.0);
			}

			// Boundary conditions
			size_t bc[dim];

			for (size_t i = 0; i < dim; i++)
				bc[i] = PERIODIC;

incardon's avatar
incardon committed
238
			vector_dist<dim,float, aggregate<float[dim]> > vd(k_int,box,bc,Ghost<dim,float>(ghost_part));
239

incardon's avatar
incardon committed
240
			vector_dist<dim,float, aggregate<float[dim]> > vd2(k_int,box,bc,Ghost<dim,float>(ghost_part));
241 242 243 244

			// Initialize dist vectors
			vd_initialize_double<dim>(vd, vd2, v_cl, k_int);

245 246
			vd.ghost_get<0>();
			vd2.ghost_get<0>();
247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274

			//Get a cell list

			auto NN = vd.getCellList(r_cut);

			//Calculate forces

			calc_forces<dim>(NN,vd,r_cut);

			//Get a cell list hilb

			auto NN_hilb = vd2.getCellList_hilb(r_cut);

			//Calculate forces
			calc_forces_hilb<dim>(NN_hilb,vd2,r_cut);

			// Calculate average
			size_t count = 1;
			Point<dim,float> avg;
			for (size_t i = 0 ; i < dim ; i++)	{avg.get(i) = 0.0;}

			auto it_v2 = vd.getIterator();
			while (it_v2.isNext())
			{
				//key
				vect_dist_key_dx key = it_v2.get();

				for (size_t i = 0; i < dim; i++)
incardon's avatar
incardon committed
275
				{avg.get(i) += fabs(vd.getProp<0>(key)[i]);}
276 277 278 279 280 281 282 283 284 285 286 287 288 289 290

				++count;
				++it_v2;
			}

			for (size_t i = 0 ; i < dim ; i++)	{avg.get(i) /= count;}

			auto it_v = vd.getIterator();
			while (it_v.isNext())
			{
				//key
				vect_dist_key_dx key = it_v.get();

				for (size_t i = 0; i < dim; i++)
				{
291 292
					auto a1 = vd.getProp<0>(key)[i];
					auto a2 = vd2.getProp<0>(key)[i];
293 294 295 296 297 298

					//Check that the forces are (almost) equal
					float per = 0.1;
					if (a1 != 0.0)
						per = fabs(0.1*avg.get(i)/a1);

incardon's avatar
incardon committed
299
					BOOST_REQUIRE_CLOSE((float)a1,(float)a2,per);
300 301 302 303 304 305 306 307 308 309
				}

				++it_v;
			}
		}
	}
}

BOOST_AUTO_TEST_CASE( vector_dist_cl_random_vs_reorder_forces_test )
{
310
	Vcluster<> & v_cl = create_vcluster();
311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340

	if (v_cl.getProcessingUnits() > 48)
		return;

	///////////////////// INPUT DATA //////////////////////

	// Dimensionality of the space
	const size_t dim = 3;
	// Cut-off radiuses. Can be put different number of values
	openfpm::vector<float> cl_r_cutoff {0.01};
	// The starting amount of particles (remember that this number is multiplied by number of processors you use for testing)
	size_t cl_k_start = 10000;
	// The lower threshold for number of particles
	size_t cl_k_min = 1000;
	// Ghost part of distributed vector
	double ghost_part = 0.01;

	///////////////////////////////////////////////////////

	//For different r_cut
	for (size_t r = 0; r < cl_r_cutoff.size(); r++ )
	{
		//Cut-off radius
		float r_cut = cl_r_cutoff.get(r);

		//Number of particles
		size_t k = cl_k_start * v_cl.getProcessingUnits();

		std::string str("Testing " + std::to_string(dim) + "D vector's forces (random vs reorder) k<=");

incardon's avatar
incardon committed
341
		print_test_v(str,k);
342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361

		//For different number of particles
		for (size_t k_int = k ; k_int >= cl_k_min ; k_int/=2 )
		{
			BOOST_TEST_CHECKPOINT( "Testing " << dim << "D vector's forces (random vs reorder) k<=" << k_int );

			Box<dim,float> box;

			for (size_t i = 0; i < dim; i++)
			{
				box.setLow(i,0.0);
				box.setHigh(i,1.0);
			}

			// Boundary conditions
			size_t bc[dim];

			for (size_t i = 0; i < dim; i++)
				bc[i] = PERIODIC;

incardon's avatar
incardon committed
362
			vector_dist<dim,float, aggregate<float[dim], float[dim]> > vd(k_int,box,bc,Ghost<dim,float>(ghost_part));
363 364

			// Initialize vd
incardon's avatar
incardon committed
365
			vd_initialize<dim,decltype(vd)>(vd, v_cl);
366

367
			vd.ghost_get<0>();
368 369 370 371 372 373 374 375 376 377 378 379 380

			//Get a cell list

			auto NN1 = vd.getCellList(r_cut);

			//Calculate forces '0'

			calc_forces<dim>(NN1,vd,r_cut);

			//Reorder and get a cell list again

			vd.reorder(4);

381
			vd.ghost_get<0>();
382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399

			auto NN2 = vd.getCellList(r_cut);

			//Calculate forces '1'
			calc_forces<dim,1>(NN2,vd,r_cut);

			// Calculate average (For Coverty scan we start from 1)
			size_t count = 1;
			Point<dim,float> avg;
			for (size_t i = 0 ; i < dim ; i++)	{avg.get(i) = 0.0;}

			auto it_v2 = vd.getIterator();
			while (it_v2.isNext())
			{
				//key
				vect_dist_key_dx key = it_v2.get();

				for (size_t i = 0; i < dim; i++)
400
					avg.get(i) += fabs(vd.getProp<0>(key)[i]);
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417

				++count;
				++it_v2;
			}

			for (size_t i = 0 ; i < dim ; i++)	{avg.get(i) /= count;}

			//Test for equality of forces
			auto it_v = vd.getDomainIterator();

			while (it_v.isNext())
			{
				//key
				vect_dist_key_dx key = it_v.get();

				for (size_t i = 0; i < dim; i++)
				{
incardon's avatar
incardon committed
418 419
					float a1 = vd.getProp<0>(key)[i];
					float a2 = vd.getProp<1>(key)[i];
420 421 422 423 424 425 426 427 428 429 430 431 432 433 434

					//Check that the forces are (almost) equal
					float per = 0.1;
					if (a1 != 0.0)
						per = fabs(0.1*avg.get(i)/a1);

					BOOST_REQUIRE_CLOSE(a1,a2,per);
				}

				++it_v;
			}
		}
	}
}

incardon's avatar
incardon committed
435 436
BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
{
437
	Vcluster<> & v_cl = create_vcluster();
incardon's avatar
incardon committed
438 439 440 441

	if (v_cl.getProcessingUnits() > 24)
		return;

442 443
	float L = 1000.0;

incardon's avatar
incardon committed
444 445
    // set the seed
	// create the random generator engine
446
	std::srand(0);
incardon's avatar
incardon committed
447
    std::default_random_engine eg;
448
    std::uniform_real_distribution<float> ud(-L,L);
incardon's avatar
incardon committed
449 450 451 452 453 454

    long int k = 4096 * v_cl.getProcessingUnits();

	long int big_step = k / 4;
	big_step = (big_step == 0)?1:big_step;

incardon's avatar
incardon committed
455
	print_test_v("Testing 3D periodic vector symmetric cell-list k=",k);
incardon's avatar
incardon committed
456 457
	BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );

458
	Box<3,float> box({-L,-L,-L},{L,L,L});
incardon's avatar
incardon committed
459 460 461 462

	// Boundary conditions
	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};

463
	float r_cut = 100.0;
incardon's avatar
incardon committed
464 465 466 467

	// ghost
	Ghost<3,float> ghost(r_cut);

468 469 470 471 472 473
	// Point and global id
	struct point_and_gid
	{
		size_t id;
		Point<3,float> xq;

474
		bool operator<(const struct point_and_gid & pag) const
475 476 477 478 479 480
		{
			return (id < pag.id);
		}
	};

	typedef  aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
incardon's avatar
incardon committed
481 482

	// Distributed vector
incardon's avatar
incardon committed
483
	vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
484
	size_t start = vd.init_size_accum(k);
incardon's avatar
incardon committed
485 486 487 488 489 490 491

	auto it = vd.getIterator();

	while (it.isNext())
	{
		auto key = it.get();

492 493 494
		vd.getPosWrite(key)[0] = ud(eg);
		vd.getPosWrite(key)[1] = ud(eg);
		vd.getPosWrite(key)[2] = ud(eg);
incardon's avatar
incardon committed
495 496 497

		// Fill some properties randomly

498 499 500
		vd.getPropWrite<0>(key) = 0;
		vd.getPropWrite<1>(key) = 0;
		vd.getPropWrite<2>(key) = key.getKey() + start;
incardon's avatar
incardon committed
501 502 503 504 505 506 507

		++it;
	}

	vd.map();

	// sync the ghost
508
	vd.ghost_get<0,2>();
incardon's avatar
incardon committed
509

510
	auto NN = vd.getCellList(r_cut);
incardon's avatar
incardon committed
511
	auto p_it = vd.getDomainIterator();
512

incardon's avatar
incardon committed
513 514 515 516
	while (p_it.isNext())
	{
		auto p = p_it.get();

517
		Point<3,float> xp = vd.getPosRead(p);
incardon's avatar
incardon committed
518

incardon's avatar
incardon committed
519
		auto Np = NN.getNNIterator(NN.getCell(xp));
incardon's avatar
incardon committed
520 521 522 523 524

		while (Np.isNext())
		{
			auto q = Np.get();

525 526 527 528 529 530
			if (p.getKey() == q)
			{
				++Np;
				continue;
			}

incardon's avatar
incardon committed
531 532
			// repulsive

533
			Point<3,float> xq = vd.getPosRead(q);
incardon's avatar
incardon committed
534 535 536 537 538 539 540
			Point<3,float> f = (xp - xq);

			float distance = f.norm();

			// Particle should be inside 2 * r_cut range

			if (distance < r_cut )
541
			{
542 543 544 545
				vd.getPropWrite<0>(p)++;
				vd.getPropWrite<3>(p).add();
				vd.getPropWrite<3>(p).last().xq = xq;
				vd.getPropWrite<3>(p).last().id = vd.getPropWrite<2>(q);
546
			}
incardon's avatar
incardon committed
547 548 549 550 551 552 553 554 555

			++Np;
		}

		++p_it;
	}

	// We now try symmetric  Cell-list

556
	auto NN2 = vd.getCellListSym(r_cut);
incardon's avatar
incardon committed
557 558 559 560 561 562 563

	auto p_it2 = vd.getDomainIterator();

	while (p_it2.isNext())
	{
		auto p = p_it2.get();

564
		Point<3,float> xp = vd.getPosRead(p);
incardon's avatar
incardon committed
565

incardon's avatar
incardon committed
566
		auto Np = NN2.getNNIteratorSym<NO_CHECK>(NN2.getCell(xp),p.getKey(),vd.getPosVector());
incardon's avatar
incardon committed
567 568 569 570 571

		while (Np.isNext())
		{
			auto q = Np.get();

572 573 574 575 576 577
			if (p.getKey() == q)
			{
				++Np;
				continue;
			}

incardon's avatar
incardon committed
578 579
			// repulsive

580
			Point<3,float> xq = vd.getPosRead(q);
incardon's avatar
incardon committed
581 582 583 584 585 586 587 588
			Point<3,float> f = (xp - xq);

			float distance = f.norm();

			// Particle should be inside r_cut range

			if (distance < r_cut )
			{
589 590
				vd.getPropWrite<1>(p)++;
				vd.getPropWrite<1>(q)++;
591

592 593
				vd.getPropWrite<4>(p).add();
				vd.getPropWrite<4>(q).add();
594

595 596 597 598
				vd.getPropWrite<4>(p).last().xq = xq;
				vd.getPropWrite<4>(q).last().xq = xp;
				vd.getPropWrite<4>(p).last().id = vd.getProp<2>(q);
				vd.getPropWrite<4>(q).last().id = vd.getProp<2>(p);
incardon's avatar
incardon committed
599 600 601 602 603 604 605 606 607
			}

			++Np;
		}

		++p_it2;
	}

	vd.ghost_put<add_,1>();
608 609 610 611 612 613 614 615 616
	vd.ghost_put<merge_,4>();

	auto p_it3 = vd.getDomainIterator();

	bool ret = true;
	while (p_it3.isNext())
	{
		auto p = p_it3.get();

617
		ret &= vd.getPropRead<1>(p) == vd.getPropRead<0>(p);
618

incardon's avatar
incardon committed
619 620
		vd.getPropWrite<3>(p).sort();
		vd.getPropWrite<4>(p).sort();
621

622
		ret &= vd.getPropRead<3>(p).size() == vd.getPropRead<4>(p).size();
623

624 625
		for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
			ret &= vd.getPropRead<3>(p).get(i).id == vd.getPropRead<4>(p).get(i).id;
626

incardon's avatar
incardon committed
627
		if (ret == false)
incardon's avatar
incardon committed
628
		{
629
			std::cout << vd.getPropRead<3>(p).size() << "   " << vd.getPropRead<4>(p).size() << std::endl;
incardon's avatar
incardon committed
630

631 632
			for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
				std::cout << vd.getPropRead<3>(p).get(i).id << "    " << vd.getPropRead<4>(p).get(i).id << std::endl;
incardon's avatar
incardon committed
633

634
			std::cout << vd.getPropRead<1>(p) << "  A  " << vd.getPropRead<0>(p) << std::endl;
incardon's avatar
incardon committed
635

incardon's avatar
incardon committed
636
			break;
incardon's avatar
incardon committed
637
		}
incardon's avatar
incardon committed
638

639 640 641 642
		++p_it3;
	}

	BOOST_REQUIRE_EQUAL(ret,true);
incardon's avatar
incardon committed
643
}
644

645 646
BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
{
647
	Vcluster<> & v_cl = create_vcluster();
648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663

	if (v_cl.getProcessingUnits() > 24)
		return;

	float L = 1000.0;

    // set the seed
	// create the random generator engine
    std::default_random_engine eg(1132312*v_cl.getProcessUnitID());
    std::uniform_real_distribution<float> ud(-L,L);

    long int k = 4096 * v_cl.getProcessingUnits();

	long int big_step = k / 4;
	big_step = (big_step == 0)?1:big_step;

incardon's avatar
incardon committed
664
	print_test_v("Testing 3D periodic vector symmetric crs cell-list k=",k);
incardon's avatar
incardon committed
665
	BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric crs cell-list k=" << k );
666 667 668 669 670 671 672 673 674 675

	Box<3,float> box({-L,-L,-L},{L,L,L});

	// Boundary conditions
	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};

	float r_cut = 100.0;

	// ghost
	Ghost<3,float> ghost(r_cut);
incardon's avatar
incardon committed
676 677 678 679
	Ghost<3,float> ghost2(r_cut);
	ghost2.setLow(0,0.0);
	ghost2.setLow(1,0.0);
	ghost2.setLow(2,0.0);
680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696

	// Point and global id
	struct point_and_gid
	{
		size_t id;
		Point<3,float> xq;

		bool operator<(const struct point_and_gid & pag) const
		{
			return (id < pag.id);
		}
	};

	typedef  aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;

	// Distributed vector
	vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
incardon's avatar
incardon committed
697
	vector_dist<3,float, part_prop > vd2(k,box,bc,ghost2,BIND_DEC_TO_GHOST);
698 699 700 701 702 703 704 705
	size_t start = vd.init_size_accum(k);

	auto it = vd.getIterator();

	while (it.isNext())
	{
		auto key = it.get();

706 707 708
		vd.getPosWrite(key)[0] = ud(eg);
		vd.getPosWrite(key)[1] = ud(eg);
		vd.getPosWrite(key)[2] = ud(eg);
709

710 711 712
		vd2.getPosWrite(key)[0] = vd.getPosRead(key)[0];
		vd2.getPosWrite(key)[1] = vd.getPosRead(key)[1];
		vd2.getPosWrite(key)[2] = vd.getPosRead(key)[2];
incardon's avatar
incardon committed
713

714 715
		// Fill some properties randomly

716 717 718
		vd.getPropWrite<0>(key) = 0;
		vd.getPropWrite<1>(key) = 0;
		vd.getPropWrite<2>(key) = key.getKey() + start;
719

720 721 722
		vd2.getPropWrite<0>(key) = 0;
		vd2.getPropWrite<1>(key) = 0;
		vd2.getPropWrite<2>(key) = key.getKey() + start;
incardon's avatar
incardon committed
723

724 725 726 727
		++it;
	}

	vd.map();
incardon's avatar
incardon committed
728
	vd2.map();
729 730 731

	// sync the ghost
	vd.ghost_get<0,2>();
incardon's avatar
incardon committed
732 733
	vd2.ghost_get<0,2>();

734 735 736 737 738 739 740
	auto NN = vd.getCellList(r_cut);
	auto p_it = vd.getDomainIterator();

	while (p_it.isNext())
	{
		auto p = p_it.get();

741
		Point<3,float> xp = vd.getPosRead(p);
742

incardon's avatar
incardon committed
743
		auto Np = NN.getNNIterator(NN.getCell(xp));
744 745 746 747 748 749 750 751 752 753 754 755 756

		while (Np.isNext())
		{
			auto q = Np.get();

			if (p.getKey() == q)
			{
				++Np;
				continue;
			}

			// repulsive

757
			Point<3,float> xq = vd.getPosRead(q);
758 759 760 761 762 763 764 765
			Point<3,float> f = (xp - xq);

			float distance = f.norm();

			// Particle should be inside 2 * r_cut range

			if (distance < r_cut )
			{
766 767 768 769
				vd.getPropWrite<0>(p)++;
				vd.getPropWrite<3>(p).add();
				vd.getPropWrite<3>(p).last().xq = xq;
				vd.getPropWrite<3>(p).last().id = vd.getPropRead<2>(q);
770 771 772 773 774 775 776 777 778 779
			}

			++Np;
		}

		++p_it;
	}

	// We now try symmetric  Cell-list

incardon's avatar
incardon committed
780
	auto NN2 = vd2.getCellListSym(r_cut);
781 782 783

	// In case of CRS we have to iterate particles within some cells
	// here we define whichone
incardon's avatar
incardon committed
784
	auto p_it2 = vd2.getParticleIteratorCRS_Cell(NN2);
785 786 787 788 789 790

	// For each particle
	while (p_it2.isNext())
	{
		auto p = p_it2.get();

791
		Point<3,float> xp = vd2.getPosRead(p);
792

incardon's avatar
incardon committed
793
		auto Np = p_it2.getNNIteratorCSR(vd2.getPosVector());
794 795 796 797 798 799 800 801 802 803 804 805 806

		while (Np.isNext())
		{
			auto q = Np.get();

			if (p == q)
			{
				++Np;
				continue;
			}

			// repulsive

807
			Point<3,float> xq = vd2.getPosRead(q);
808 809 810 811 812 813 814 815
			Point<3,float> f = (xp - xq);

			float distance = f.norm();

			// Particle should be inside r_cut range

			if (distance < r_cut )
			{
816 817
				vd2.getPropWrite<1>(p)++;
				vd2.getPropWrite<1>(q)++;
818

819 820
				vd2.getPropWrite<4>(p).add();
				vd2.getPropWrite<4>(q).add();
821

822 823 824 825
				vd2.getPropWrite<4>(p).last().xq = xq;
				vd2.getPropWrite<4>(q).last().xq = xp;
				vd2.getPropWrite<4>(p).last().id = vd2.getPropRead<2>(q);
				vd2.getPropWrite<4>(q).last().id = vd2.getPropRead<2>(p);
826 827 828 829 830 831 832 833
			}

			++Np;
		}

		++p_it2;
	}

incardon's avatar
incardon committed
834
	vd2.ghost_put<add_,1>(NO_CHANGE_ELEMENTS);
incardon's avatar
incardon committed
835
	vd2.ghost_put<merge_,4>();
836

837 838 839 840
#ifdef SE_CLASS3
	vd2.getDomainIterator();
#endif

841 842 843 844 845 846 847
	auto p_it3 = vd.getDomainIterator();

	bool ret = true;
	while (p_it3.isNext())
	{
		auto p = p_it3.get();

848
		ret &= vd2.getPropRead<1>(p) == vd.getPropRead<0>(p);
849 850


851 852
		vd.getPropWrite<3>(p).sort();
		vd2.getPropWrite<4>(p).sort();
853

854
		ret &= vd.getPropRead<3>(p).size() == vd2.getPropRead<4>(p).size();
855

856 857
		for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
			ret &= vd.getPropRead<3>(p).get(i).id == vd2.getPropRead<4>(p).get(i).id;
858 859 860 861 862 863 864 865 866

		if (ret == false)
			break;

		++p_it3;
	}

	BOOST_REQUIRE_EQUAL(ret,true);
}
incardon's avatar
incardon committed
867

868 869
template<typename VerletList>
void test_vd_symmetric_verlet_list()
870
{
871
	Vcluster<> & v_cl = create_vcluster();
incardon's avatar
incardon committed
872 873

	if (v_cl.getProcessingUnits() > 24)
874 875
		return;

incardon's avatar
incardon committed
876
	float L = 1000.0;
877

incardon's avatar
incardon committed
878 879 880 881 882 883 884 885 886
    // set the seed
	// create the random generator engine
	std::srand(0);
    std::default_random_engine eg;
    std::uniform_real_distribution<float> ud(-L,L);

    long int k = 4096 * v_cl.getProcessingUnits();

	long int big_step = k / 4;
887 888
	big_step = (big_step == 0)?1:big_step;

incardon's avatar
incardon committed
889
	print_test_v("Testing 3D periodic vector symmetric cell-list k=",k);
incardon's avatar
incardon committed
890
	BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric verlet-list k=" << k );
891

incardon's avatar
incardon committed
892
	Box<3,float> box({-L,-L,-L},{L,L,L});
893

incardon's avatar
incardon committed
894 895
	// Boundary conditions
	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
896

incardon's avatar
incardon committed
897
	float r_cut = 100.0;
898

incardon's avatar
incardon committed
899 900
	// ghost
	Ghost<3,float> ghost(r_cut);
901

incardon's avatar
incardon committed
902 903 904 905 906
	// Point and global id
	struct point_and_gid
	{
		size_t id;
		Point<3,float> xq;
907

incardon's avatar
incardon committed
908
		bool operator<(const struct point_and_gid & pag) const
909
		{
incardon's avatar
incardon committed
910 911 912
			return (id < pag.id);
		}
	};
913

incardon's avatar
incardon committed
914
	typedef  aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
915

incardon's avatar
incardon committed
916
	// Distributed vector
incardon's avatar
incardon committed
917
	vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
incardon's avatar
incardon committed
918
	size_t start = vd.init_size_accum(k);
919

incardon's avatar
incardon committed
920
	auto it = vd.getIterator();
921

incardon's avatar
incardon committed
922 923 924
	while (it.isNext())
	{
		auto key = it.get();
925

926 927 928
		vd.getPosWrite(key)[0] = ud(eg);
		vd.getPosWrite(key)[1] = ud(eg);
		vd.getPosWrite(key)[2] = ud(eg);
929

incardon's avatar
incardon committed
930
		// Fill some properties randomly
931

932 933 934
		vd.template getPropWrite<0>(key) = 0;
		vd.template getPropWrite<1>(key) = 0;
		vd.template getPropWrite<2>(key) = key.getKey() + start;
935

incardon's avatar
incardon committed
936 937
		++it;
	}
938

incardon's avatar
incardon committed
939
	vd.map();
940

incardon's avatar
incardon committed
941
	// sync the ghost
942
	vd.template ghost_get<0,2>();
943

944
	auto NN = vd.template getVerlet<VerletList>(r_cut);
incardon's avatar
incardon committed
945
	auto p_it = vd.getDomainIterator();
946

incardon's avatar
incardon committed
947 948 949
	while (p_it.isNext())
	{
		auto p = p_it.get();
950

951
		Point<3,float> xp = vd.getPosRead(p);
952

incardon's avatar
incardon committed
953
		auto Np = NN.getNNIterator(p.getKey());
954

incardon's avatar
incardon committed
955 956 957
		while (Np.isNext())
		{
			auto q = Np.get();
958

incardon's avatar
incardon committed
959 960 961 962 963
			if (p.getKey() == q)
			{
				++Np;
				continue;
			}
964

incardon's avatar
incardon committed
965
			// repulsive
966

967
			Point<3,float> xq = vd.getPosRead(q);
incardon's avatar
incardon committed
968 969 970 971 972
			Point<3,float> f = (xp - xq);

			float distance = f.norm();

			// Particle should be inside 2 * r_cut range
973

incardon's avatar
incardon committed
974 975
			if (distance < r_cut )
			{
976 977 978 979
				vd.template getPropWrite<0>(p)++;
				vd.template getPropWrite<3>(p).add();
				vd.template getPropWrite<3>(p).last().xq = xq;
				vd.template getPropWrite<3>(p).last().id = vd.template getPropRead<2>(q);
980 981
			}

incardon's avatar
incardon committed
982
			++Np;
983 984
		}

incardon's avatar
incardon committed
985 986
		++p_it;
	}
987

incardon's avatar
incardon committed
988
	// We now try symmetric  Cell-list
989

990
	auto NN2 = vd.template getVerletSym<VerletList>(r_cut);
991

incardon's avatar
incardon committed
992 993 994 995 996
	auto p_it2 = vd.getDomainIterator();

	while (p_it2.isNext())
	{
		auto p = p_it2.get();
997

998
		Point<3,float> xp = vd.getPosRead(p);
999

1000
		auto Np = NN2.template getNNIterator<NO_CHECK>(p.getKey());
1001

incardon's avatar
incardon committed
1002 1003 1004
		while (Np.isNext())
		{
			auto q = Np.get();
1005

incardon's avatar
incardon committed
1006
			if (p.getKey() == q)
1007
			{
incardon's avatar
incardon committed
1008 1009 1010
				++Np;
				continue;
			}
1011

incardon's avatar
incardon committed
1012
			// repulsive
1013

1014
			Point<3,float> xq = vd.getPosRead(q);
incardon's avatar
incardon committed
1015
			Point<3,float> f = (xp - xq);
1016

incardon's avatar
incardon committed
1017
			float distance = f.norm();
1018

incardon's avatar
incardon committed
1019
			// Particle should be inside r_cut range
1020

incardon's avatar
incardon committed
1021 1022
			if (distance < r_cut )
			{
1023 1024
				vd.template getPropWrite<1>(p)++;
				vd.template getPropWrite<1>(q)++;
1025

1026 1027
				vd.template getPropWrite<4>(p).add();
				vd.template getPropWrite<4>(q).add();
incardon's avatar
incardon committed
1028

1029 1030 1031 1032
				vd.template getPropWrite<4>(p).last().xq = xq;
				vd.template getPropWrite<4>(q).last().xq = xp;
				vd.template getPropWrite<4>(p).last().id = vd.template getPropRead<2>(q);
				vd.template getPropWrite<4>(q).last().id = vd.template getPropRead<2>(p);
1033 1034
			}

incardon's avatar
incardon committed
1035 1036
			++Np;
		}
1037

incardon's avatar
incardon committed
1038 1039
		++p_it2;
	}
1040

1041 1042
	vd.template ghost_put<add_,1>();
	vd.template ghost_put<merge_,4>();
1043

incardon's avatar
incardon committed
1044
	auto p_it3 = vd.getDomainIterator();
1045

incardon's avatar
incardon committed
1046 1047 1048 1049
	bool ret = true;
	while (p_it3.isNext())
	{
		auto p = p_it3.get();
1050

1051
		ret &= vd.template getPropRead<1>(p) == vd.template getPropRead<0>(p);
1052

1053 1054
		vd.template getPropWrite<3>(p).sort();
		vd.template getPropWrite<4>(p).sort();
incardon's avatar
incardon committed
1055

1056
		ret &= vd.template getPropRead<3>(p).size() == vd.template getPropRead<4>(p).size();
incardon's avatar
incardon committed
1057

1058 1059
		for (size_t i = 0 ; i < vd.template getPropRead<3>(p).size() ; i++)
			ret &= vd.template getPropRead<3>(p).get(i).id == vd.template getPropRead<4>(p).get(i).id;
incardon's avatar
incardon committed
1060 1061 1062 1063 1064

		if (ret == false)
			break;

		++p_it3;
1065
	}
incardon's avatar
incardon committed
1066 1067

	BOOST_REQUIRE_EQUAL(ret,true);
1068 1069
}

1070 1071 1072 1073 1074 1075 1076 1077 1078
BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
{
	test_vd_symmetric_verlet_list<VERLET_MEMFAST(3,float)>();
	test_vd_symmetric_verlet_list<VERLET_MEMBAL(3,float)>();
	test_vd_symmetric_verlet_list<VERLET_MEMMW(3,float)>();
}

template<typename VerletList>
void vector_sym_verlet_list_nb()
1079
{
1080
	Vcluster<> & v_cl = create_vcluster();
1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097

	if (v_cl.getProcessingUnits() > 24)
		return;

	float L = 1000.0;

    // set the seed
	// create the random generator engine
	std::srand(0);
    std::default_random_engine eg;
    std::uniform_real_distribution<float> ud(-L,L);

    long int k = 4096 * v_cl.getProcessingUnits();

	long int big_step = k / 4;
	big_step = (big_step == 0)?1:big_step;

incardon's avatar
incardon committed
1098
	print_test_v("Testing 3D periodic vector symmetric cell-list no bottom k=",k);
1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126
	BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list no bottom k=" << k );

	Box<3,float> box({-L,-L,-L},{L,L,L});

	// Boundary conditions
	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};

	float r_cut = 100.0;

	// ghost
	Ghost<3,float> ghost(r_cut);
	Ghost<3,float> ghost2(r_cut);
	ghost2.setLow(2,0.0);

	// Point and global id
	struct point_and_gid
	{
		size_t id;
		Point<3,float> xq;

		bool operator<(const struct point_and_gid & pag) const
		{
			return (id < pag.id);
		}
	};

	typedef  aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;

1127 1128
	// 3D test
	for (size_t s = 0 ; s < 8 ; s++)
1129 1130
	{

1131 1132 1133 1134
		// Distributed vector
		vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
		vector_dist<3,float, part_prop > vd2(k,box,bc,ghost2,BIND_DEC_TO_GHOST);
		size_t start = vd.init_size_accum(k);
1135

1136
		auto it = vd.getIterator();
1137

1138 1139 1140
		while (it.isNext())
		{
			auto key = it.get();
1141

1142 1143 1144
			vd.getPosWrite(key)[0] = ud(eg);
			vd.getPosWrite(key)[1] = ud(eg);
			vd.getPosWrite(key)[2] = ud(eg);
1145

1146 1147 1148
			vd2.getPosWrite(key)[0] = vd.getPosRead(key)[0];
			vd2.getPosWrite(key)[1] = vd.getPosRead(key)[1];
			vd2.getPosWrite(key)[2] = vd.getPosRead(key)[2];
1149

1150
			// Fill some properties randomly
1151

1152 1153 1154
			vd.template getPropWrite<0>(key) = 0;
			vd.template getPropWrite<1>(key) = 0;
			vd.template getPropWrite<2>(key) = key.getKey() + start;
1155

1156 1157 1158
			vd2.template getPropWrite<0>(key) = 0;
			vd2.template getPropWrite<1>(key) = 0;
			vd2.template getPropWrite<2>(key) = key.getKey() + start;
1159

1160 1161
			++it;
		}
1162

1163 1164
		vd.map();
		vd2.map();
1165

1166
		// sync the ghost
1167 1168
		vd.template ghost_get<0,2>();
		vd2.template ghost_get<0,2>();
1169

1170
		auto NN = vd.template getVerlet<VerletList>(r_cut);
1171
		auto p_it = vd.getDomainIterator();
1172

1173
		while (p_it.isNext())
1174
		{
1175
			auto p = p_it.get();
1176

1177
			Point<3,float> xp = vd.getPosRead(p);
1178 1179 1180 1181

			auto Np = NN.getNNIterator(p.getKey());

			while (Np.isNext())
1182
			{
1183
				auto q = Np.get();
1184

1185 1186 1187 1188 1189
				if (p.getKey() == q)
				{
					++Np;
					continue;
				}
1190

1191
				// repulsive
1192

1193
				Point<3,float> xq = vd.getPosRead(q);
1194
				Point<3,float> f = (xp - xq);
1195

1196
				float distance = f.norm();
1197

1198 1199 1200 1201
				// Particle should be inside 2 * r_cut range

				if (distance < r_cut )
				{
1202 1203 1204 1205
					vd.template getPropWrite<0>(p)++;
					vd.template getPropWrite<3>(p).add();
					vd.template getPropWrite<3>(p).last().xq = xq;
					vd.template getPropWrite<3>(p).last().id = vd.template getPropRead<2>(q);
1206 1207 1208
				}

				++Np;
1209 1210
			}

1211
			++p_it;
1212 1213
		}

1214
		// We now try symmetric  Cell-list
1215

1216
		auto NN2 = vd2.template getVerletSym<VerletList>(r_cut);
1217

1218
		auto p_it2 = vd2.getDomainIterator();
1219

1220 1221 1222
		while (p_it2.isNext())
		{
			auto p = p_it2.get();
1223

1224
			Point<3,float> xp = vd2.getPosRead(p);
1225

1226
			auto Np = NN2.template getNNIterator<NO_CHECK>(p.getKey());
1227

1228 1229 1230
			while (Np.isNext())
			{
				auto q = Np.get();
1231

1232 1233 1234 1235 1236
				if (p.getKey() == q)
				{
					++Np;
					continue;
				}
1237

1238
				// repulsive
1239

1240
				Point<3,float> xq = vd2.getPosRead(q);
1241
				Point<3,float> f = (xp - xq);
1242

1243
				float distance = f.norm();
1244

1245
				// Particle should be inside r_cut range
1246

1247 1248
				if (distance < r_cut )
				{
1249 1250
					vd2.template getPropWrite<1>(p)++;
					vd2.template getPropWrite<1>(q)++;
1251

1252 1253
					vd2.template getPropWrite<4>(p).add();
					vd2.template getPropWrite<4>(q).add();
1254

1255 1256 1257 1258
					vd2.template getPropWrite<4>(p).last().xq = xq;
					vd2.template getPropWrite<4>(q).last().xq = xp;
					vd2.template getPropWrite<4>(p).last().id = vd2.template getPropRead<2>(q);
					vd2.template getPropWrite<4>(q).last().id = vd2.template getPropRead<2>(p);
1259
				}
1260

1261
				++Np;
1262 1263
			}

1264 1265

			++p_it2;
1266 1267
		}

1268 1269
		vd2.template ghost_put<add_,1>();
		vd2.template ghost_put<merge_,4>();
1270

1271 1272 1273 1274
#ifdef SE_CLASS3
		vd2.getDomainIterator();
#endif

1275
		auto p_it3 = vd.getDomainIterator();
1276

1277 1278 1279 1280
		bool ret = true;
		while (p_it3.isNext())
		{
			auto p = p_it3.get();
1281

1282
			ret &= vd2.template getPropRead<1>(p) == vd.template getPropRead<0>(p);
1283

1284 1285
			vd.template getPropWrite<3>(p).sort();
			vd2.template getPropWrite<4>(p).sort();
1286

1287
			ret &= vd.template getPropRead<3>(p).size() == vd2.template getPropRead<4>(p).size();
1288

1289 1290
			for (size_t i = 0 ; i < vd.template getPropRead<3>(p).size() ; i++)
				ret &= vd.template getPropRead<3>(p).get(i).id == vd2.template getPropRead<4>(p).get(i).id;
1291

1292 1293
			if (ret == false)
				break;
1294

1295 1296
			++p_it3;
		}
1297

1298
		BOOST_REQUIRE_EQUAL(ret,true);
1299 1300
	}
}
1301

1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313
BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
{
	vector_sym_verlet_list_nb<VERLET_MEMFAST(3,float)>();
	vector_sym_verlet_list_nb<VERLET_MEMBAL(3,float)>();
	vector_sym_verlet_list_nb<VERLET_MEMMW(3,float)>();

	vector_sym_verlet_list_nb<VERLET_MEMFAST_INT(3,float)>();
	vector_sym_verlet_list_nb<VERLET_MEMBAL_INT(3,float)>();
	vector_sym_verlet_list_nb<VERLET_MEMMW_INT(3,float)>();
}

template<typename VerletList, typename part_prop> void test_crs_full(vector_dist<3,float, part_prop > & vd,
1314 1315 1316 1317 1318
		                                        vector_dist<3,float, part_prop > & vd2,
												std::default_random_engine & eg,
												std::uniform_real_distribution<float> & ud,
												size_t start,
												float r_cut)
incardon's avatar
incardon committed
1319 1320 1321 1322 1323 1324 1325
{
	auto it = vd.getIterator();

	while (it.isNext())
	{
		auto key = it.get();

1326 1327 1328
		vd.getPosWrite(key)[0] = ud(eg);
		vd.getPosWrite(key)[1] = ud(eg);
		vd.getPosWrite(key)[2] = ud(eg);
incardon's avatar
incardon committed
1329

1330 1331 1332
		vd2.getPosWrite(key)[0] = vd.getPosRead(key)[0];
		vd2.getPosWrite(key)[1] = vd.getPosRead(key)[1];
		vd2.getPosWrite(key)[2] = vd.getPosRead(key)[2];
incardon's avatar
incardon committed
1333 1334 1335

		// Fill some properties randomly

1336 1337 1338
		vd.template getPropWrite<0>(key) = 0;
		vd.template getPropWrite<1>(key) = 0;
		vd.template getPropWrite<2>(key) = key.getKey() + start;
incardon's avatar
incardon committed
1339

1340 1341 1342
		vd2.template getPropWrite<0>(key) = 0;
		vd2.template getPropWrite<1>(key) = 0;
		vd2.template getPropWrite<2>(key) = key.getKey() + start;
incardon's avatar
incardon committed
1343 1344 1345 1346 1347 1348 1349 1350

		++it;
	}

	vd.map();
	vd2.map();

	// sync the ghost
1351 1352
	vd.template ghost_get<0,2>();
	vd2.template ghost_get<0,2>();
incardon's avatar
incardon committed
1353

1354
	auto NN = vd.template getVerlet<VerletList>(r_cut);
incardon's avatar
incardon committed
1355 1356 1357 1358 1359 1360
	auto p_it = vd.getDomainIterator();

	while (p_it.isNext())
	{
		auto p = p_it.get();

1361
		Point<3,float> xp = vd.getPosRead(p);
1362

incardon's avatar
incardon committed
1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376
		auto Np = NN.getNNIterator(p.getKey());

		while (Np.isNext())
		{
			auto q = Np.get();

			if (p.getKey() == q)
			{
				++Np;
				continue;
			}

			// repulsive

1377
			Point<3,float> xq = vd.getPosRead(q);
incardon's avatar
incardon committed
1378 1379 1380 1381 1382 1383 1384 1385
			Point<3,float> f = (xp - xq);

			float distance = f.norm();

			// Particle should be inside 2 * r_cut range

			if (distance < r_cut )
			{
1386 1387 1388 1389
				vd.template getPropWrite<0>(p)++;
				vd.template getPropWrite<3>(p).add();
				vd.template getPropWrite<3>(p).last().xq = xq;
				vd.template getPropWrite<3>(p).last().id = vd.template getPropRead<2>(q);
incardon's avatar
incardon committed
1390 1391 1392 1393 1394 1395 1396 1397 1398 1399
			}

			++Np;
		}

		++p_it;
	}

	// We now try symmetric Verlet-list Crs scheme

1400
	auto NN2 = vd2.template getVerletCrs<VerletList>(r_cut);
incardon's avatar
incardon committed
1401 1402

	// Because iterating across particles in the CSR scheme require a Cell-list
incardon's avatar
incardon committed
1403
	auto p_it2 = vd2.getParticleIteratorCRS_Cell(NN2.getInternalCellList());
incardon's avatar
incardon committed
1404 1405 1406 1407 1408

	while (p_it2.isNext())
	{
		auto p = p_it2.get();

1409
		Point<3,float> xp = vd2.getPosRead(p);