Fixing the test

configure.ac

@@ -134,14 +134,6 @@ fi
 AX_LIB_HILBERT([],[echo "Cannot detect libhilbert, use the --with-libhilbert option if it is not installed in the default location"
                     exit 210])
 
-##########
-
-## Check for PETSC
-
-AX_LIB_PETSC()
-
-## Check for quadmath
-
 
 have_quad_lib=no
 have_quad_head=no
@@ -372,6 +364,12 @@ AX_LAPACK([],[])
 
 AX_SUITESPARSE([],[])
 
+##########
+
+## Check for PETSC
+
+AX_LIB_PETSC()
+
 ###### Checking for EIGEN
 
 AX_EIGEN([],[])

example/Numerics/Vortex_in_cell/main_vic_petsc_opt.cpp

@@ -716,7 +716,7 @@ int main(int argc, char* argv[])
 	// Grid points on x=128 y=64 z=64
 	// if we use Re = 7500
 	//  long int sz[] = {1600,400,400};
-	long int sz[] = {96,96,96};
+	long int sz[] = {128,128,128};
 	size_t szu[] = {(size_t)sz[0],(size_t)sz[1],(size_t)sz[2]};
 
 	periodicity<3> bc = {{PERIODIC,PERIODIC,PERIODIC}};

m4/ax_petsc_lib.m4

@@ -105,7 +105,7 @@ AC_DEFUN([AX_LIB_PETSC], [
 			AX_OPENMP([CFLAGS="$OPENMP_CFLAGS"
 				   LDFLAGS="$OPENMP_LDFLAGS"],[])
                         CFLAGS="$CFLAGS -I$with_petsc/include $HDF5_INCLUDE $METIS_INCLUDE "
-                        LDFLAGS="$LDFLAGS -L$with_petsc/lib $HDF5_LDFLAGS  $HDF5_LIBS $METIS_LIB -lmetis "
+                        LDFLAGS="$LDFLAGS -L$with_petsc/lib $HDF5_LDFLAGS  $HDF5_LIBS $METIS_LIB -lmetis $SUITESPARSE_LIBS"
 			CC=$CXX
 
                         AC_LANG_SAVE

openfpm_io @ a0466487

-Subproject commit 701d310eb6c3347e52327ddcf9a814c22b0dadbd
+Subproject commit a0466487a3c6a1a8c8e29f02afd4db6beb9186ea

script/install_SUITESPARSE.sh

@@ -34,7 +34,7 @@ if [ x"$platform" == x"cygwin" ]; then
 fi
 
 echo "Compiling SuiteSparse without CUDA (old variable $CUDA)"
-make "CUDA=no" "BLAS=-L$1/OPENBLAS/lib -lopenblas" "LAPACK="
+make "CUDA=no" "BLAS=-L$1/OPENBLAS/lib -lopenblas -pthread" "LAPACK="
 if [ $? != 0 ]; then
   echo "Failed to compile SuiteSparse"
   exit 1

src/Makefile.am

-LINKLIBS = $(HDF5_LDFLAGS)  $(HDF5_LIBS) $(OPENMP_LDFLAGS) $(LIBHILBERT_LIB)  $(METIS_LIB) $(PTHREAD_LIBS) $(OPT_LIBS) $(BOOST_LDFLAGS) $(BOOST_IOSTREAMS_LIB) $(CUDA_LIBS) $(PETSC_LIB) $(PARMETIS_LIB) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) $(BOOST_CHRONO_LIB) $(BOOST_TIMER_LIB) $(BOOST_SYSTEM_LIB) $(LIBIFCORE)
+LINKLIBS = $(HDF5_LDFLAGS)  $(HDF5_LIBS) $(OPENMP_LDFLAGS) $(LIBHILBERT_LIB)  $(METIS_LIB) $(PTHREAD_LIBS) $(OPT_LIBS) $(BOOST_LDFLAGS) $(BOOST_IOSTREAMS_LIB) $(CUDA_LIBS) $(PETSC_LIB) $(PARMETIS_LIB) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) $(BOOST_CHRONO_LIB) $(BOOST_TIMER_LIB) $(BOOST_SYSTEM_LIB) $(LIBIFCORE) $(SUITESPARSE_LIBS)
 
 noinst_PROGRAMS = pdata
 pdata_SOURCES = main.cpp Grid/tests/grid_dist_id_HDF5_chckpnt_restart_test.cpp Grid/tests/grid_dist_id_unit_test.cpp Grid/tests/staggered_grid_dist_unit_test.cpp Vector/tests/vector_dist_cell_list_tests.cpp Vector/tests/vector_dist_complex_prp_unit_test.cpp Vector/tests/vector_dist_HDF5_chckpnt_restart_test.cpp Vector/tests/vector_dist_MP_unit_tests.cpp Vector/tests/vector_dist_NN_tests.cpp Vector/tests/vector_dist_unit_test.cpp  pdata_performance.cpp Decomposition/tests/CartDecomposition_unit_test.cpp Decomposition/tests/shift_vect_converter_tests.cpp  lib/pdata.cpp test_multiple_o.cpp ../openfpm_devices/src/memory/HeapMemory.cpp ../openfpm_devices/src/memory/PtrMemory.cpp ../openfpm_vcluster/src/VCluster/VCluster.cpp ../openfpm_devices/src/Memleak_check.cpp

src/Vector/tests/vector_dist_cell_list_tests.cpp

@@ -1854,3 +1854,303 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_list_multi_type )
 	BOOST_REQUIRE_EQUAL(ret,true);
 }
 
+
+
+BOOST_AUTO_TEST_CASE( vector_dist_particle_NN_MP_iteration )
+{
+	Vcluster & v_cl = create_vcluster();
+
+	if (v_cl.getProcessingUnits() > 24)
+	{return;}
+
+	float L = 1000.0;
+
+    // set the seed
+	// create the random generator engine
+    std::default_random_engine eg;
+    eg.seed(v_cl.rank()*4533);
+    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;
+
+	print_test_v("Testing 3D periodic vector symmetric cell-list k=",k);
+	BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list 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);
+
+	// 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);
+	size_t start = vd.init_size_accum(k);
+
+	auto it = vd.getIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		vd.getPosWrite(key)[0] = ud(eg);
+		vd.getPosWrite(key)[1] = ud(eg);
+		vd.getPosWrite(key)[2] = ud(eg);
+
+		// Fill some properties randomly
+
+		vd.getPropWrite<0>(key) = 0;
+		vd.getPropWrite<1>(key) = 0;
+		vd.getPropWrite<2>(key) = key.getKey() + start;
+
+		++it;
+	}
+
+	vd.map();
+
+	// sync the ghost
+	vd.ghost_get<0,2>();
+
+	auto NN = vd.getCellList(r_cut);
+	auto p_it = vd.getDomainIterator();
+
+	while (p_it.isNext())
+	{
+		auto p = p_it.get();
+
+		Point<3,float> xp = vd.getPosRead(p);
+
+		auto Np = NN.getNNIterator(NN.getCell(xp));
+
+		while (Np.isNext())
+		{
+			auto q = Np.get();
+
+			if (p.getKey() == q)
+			{
+				++Np;
+				continue;
+			}
+
+			// repulsive
+
+			Point<3,float> xq = vd.getPosRead(q);
+			Point<3,float> f = (xp - xq);
+
+			float distance = f.norm();
+
+			// Particle should be inside 2 * r_cut range
+
+			if (distance < r_cut )
+			{
+				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);
+			}
+
+			++Np;
+		}
+
+		++p_it;
+	}
+
+	// We now divide the particles on 4 phases
+
+	openfpm::vector<vector_dist<3,float, part_prop >> phases;
+	phases.add( vector_dist<3,float, part_prop >(vd.getDecomposition(),0));
+	phases.add( vector_dist<3,float, part_prop >(phases.get(0).getDecomposition(),0));
+	phases.add( vector_dist<3,float, part_prop >(phases.get(0).getDecomposition(),0));
+	phases.add( vector_dist<3,float, part_prop >(phases.get(0).getDecomposition(),0));
+
+	auto it2 = vd.getDomainIterator();
+
+	while (it2.isNext())
+	{
+		auto p = it2.get();
+
+		if (p.getKey() % 4 == 0)
+		{
+			phases.get(0).add();
+			phases.get(0).getLastPos()[0] = vd.getPos(p)[0];
+			phases.get(0).getLastPos()[1] = vd.getPos(p)[1];
+			phases.get(0).getLastPos()[2] = vd.getPos(p)[2];
+
+			phases.get(0).template getLastProp<2>() = vd.template getProp<2>(p);
+		}
+		else if (p.getKey() % 4 == 1)
+		{
+			phases.get(1).add();
+			phases.get(1).getLastPos()[0] = vd.getPos(p)[0];
+			phases.get(1).getLastPos()[1] = vd.getPos(p)[1];
+			phases.get(1).getLastPos()[2] = vd.getPos(p)[2];
+
+			phases.get(1).template getLastProp<2>() = vd.template getProp<2>(p);
+		}
+		else if (p.getKey() % 4 == 2)
+		{
+			phases.get(2).add();
+			phases.get(2).getLastPos()[0] = vd.getPos(p)[0];
+			phases.get(2).getLastPos()[1] = vd.getPos(p)[1];
+			phases.get(2).getLastPos()[2] = vd.getPos(p)[2];
+
+			phases.get(2).template getLastProp<2>() = vd.template getProp<2>(p);
+		}
+		else
+		{
+			phases.get(3).add();
+			phases.get(3).getLastPos()[0] = vd.getPos(p)[0];
+			phases.get(3).getLastPos()[1] = vd.getPos(p)[1];
+			phases.get(3).getLastPos()[2] = vd.getPos(p)[2];
+
+			phases.get(3).template getLastProp<2>() = vd.template getProp<2>(p);
+		}
+
+		++it2;
+	}
+
+	// now we get all the Cell-lists
+
+	for (size_t i = 0 ; i < phases.size() ; i++)
+	{
+		phases.get(i).ghost_get<0,2>();
+	}
+
+	openfpm::vector<CellList<3, float, Mem_fast<>, shift<3, float> >> NN_ptr;
+
+	for (size_t i = 0 ; i < phases.size() ; i++)
+	{
+		NN_ptr.add(phases.get(i).getCellListSym(r_cut));
+	}
+
+	// We now interact all the phases
+
+	for (size_t i = 0 ; i < phases.size() ; i++)
+	{
+		for (size_t j = 0 ; j < phases.size() ; j++)
+		{
+			auto p_it2 = phases.get(i).getDomainIterator();
+
+			while (p_it2.isNext())
+			{
+				auto p = p_it2.get();
+
+				Point<3,float> xp = phases.get(i).getPosRead(p);
+
+				auto Np = NN_ptr.get(j).getNNIteratorSymMP<NO_CHECK>(NN_ptr.get(j).getCell(xp),p.getKey(),phases.get(i).getPosVector(),phases.get(j).getPosVector());
+
+				while (Np.isNext())
+				{
+					auto q = Np.get();
+
+					if (p.getKey() == q && i == j)
+					{
+						++Np;
+						continue;
+					}
+
+					// repulsive
+
+					Point<3,float> xq = phases.get(j).getPosRead(q);
+					Point<3,float> f = (xp - xq);
+
+					float distance = f.norm();
+
+					// Particle should be inside r_cut range
+
+					if (distance < r_cut )
+					{
+						phases.get(i).getPropWrite<1>(p)++;
+						phases.get(j).getPropWrite<1>(q)++;
+
+						phases.get(i).getPropWrite<4>(p).add();
+						phases.get(j).getPropWrite<4>(q).add();
+
+						phases.get(i).getPropWrite<4>(p).last().xq = xq;
+						phases.get(j).getPropWrite<4>(q).last().xq = xp;
+						phases.get(i).getPropWrite<4>(p).last().id = phases.get(j).getProp<2>(q);
+						phases.get(j).getPropWrite<4>(q).last().id = phases.get(i).getProp<2>(p);
+					}
+
+					++Np;
+				}
+
+				++p_it2;
+			}
+		}
+	}
+
+	for (size_t i = 0 ; i < phases.size() ; i++)
+	{
+		phases.get(i).ghost_put<add_,1>();
+		phases.get(i).ghost_put<merge_,4>();
+	}
+
+	auto p_it3 = vd.getDomainIterator();
+
+	bool ret = true;
+	while (p_it3.isNext())
+	{
+		auto p = p_it3.get();
+
+		int ph;
+
+		if (p.getKey() % 4 == 0)
+		{ph = 0;}
+		else if (p.getKey() % 4 == 1)
+		{ph = 1;}
+		else if (p.getKey() % 4 == 2)
+		{ph = 2;}
+		else
+		{ph = 3;}
+
+		size_t pah = p.getKey()/4;
+		ret &= phases.get(ph).getPropRead<1>(pah) == vd.getPropRead<0>(p);
+
+		vd.getPropRead<3>(p).sort();
+		phases.get(ph).getPropRead<4>(pah).sort();
+
+		ret &= vd.getPropRead<3>(p).size() == phases.get(ph).getPropRead<4>(pah).size();
+
+		for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
+			ret &= vd.getPropRead<3>(p).get(i).id == phases.get(ph).getPropRead<4>(pah).get(i).id;
+
+		if (ret == false)
+		{
+			std::cout << "Error on particle: " << vd.getPropRead<2>(p) << "   " << v_cl.rank() << std::endl;
+
+			std::cout << vd.getPropRead<3>(p).size() << "   " << phases.get(ph).getPropRead<4>(pah).size() << "  " << v_cl.rank() << std::endl;
+
+			for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
+				std::cout << vd.getPropRead<3>(p).get(i).id << "    " << phases.get(ph).getPropRead<4>(pah).get(i).id << "  " << v_cl.rank() << std::endl;
+
+			std::cout << phases.get(ph).getPropRead<1>(pah) << "  A  " << vd.getPropRead<0>(p) << std::endl;
+
+			break;
+		}
+
+		++p_it3;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+}

src/Vector/tests/vector_dist_unit_test.cpp

@@ -1821,6 +1821,70 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost_put )
 	}
 }
 
+BOOST_AUTO_TEST_CASE( vector_fixing_noposition_and_keep_prop )
+{
+	Vcluster & v_cl = create_vcluster();
+
+	if (v_cl.getProcessingUnits() > 48)
+		return;
+
+	// Boundary conditions
+	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+	// Box
+	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+	// ghost
+	Ghost<3,float> ghost(0.1);
+
+	vector_dist<3,float, aggregate<double,double>> vd(4096,box,bc,ghost);
+
+	auto it = vd.getDomainIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		vd.getPos(key)[0] = ((double)rand())/RAND_MAX;
+		vd.getPos(key)[1] = ((double)rand())/RAND_MAX;
+		vd.getPos(key)[2] = ((double)rand())/RAND_MAX;
+
+		++it;
+	}
+
+	vd.map();
+
+	vd.ghost_get<>();
+	size_t local = vd.getPosVector().size();
+
+	vd.ghost_get<>(KEEP_PROPERTIES | NO_POSITION);
+
+	size_t local2 = vd.getPosVector().size();
+
+	BOOST_REQUIRE_EQUAL(local,local2);
+
+	// Check now that map reset
+
+	vd.map();
+
+	local = vd.getPosVector().size();
+	BOOST_REQUIRE_EQUAL(local,vd.size_local());
+	vd.ghost_get<>(KEEP_PROPERTIES  | NO_POSITION);
+
+	local2 = vd.getPosVector().size();
+
+	BOOST_REQUIRE_EQUAL(local,local2);
+
+	vd.ghost_get<>(KEEP_PROPERTIES);
+	BOOST_REQUIRE_EQUAL(local,vd.getPosVector().size());
+	BOOST_REQUIRE_EQUAL(vd.getPropVector().size(),local);
+
+	vd.ghost_get<0>(KEEP_PROPERTIES);
+	BOOST_REQUIRE_EQUAL(local,vd.getPosVector().size());
+	BOOST_REQUIRE_EQUAL(vd.getPropVector().size(),local);
+}
+
+
 BOOST_AUTO_TEST_CASE( vector_of_vector_dist )
 {
 	Vcluster & v_cl = create_vcluster();
@@ -1868,6 +1932,5 @@ BOOST_AUTO_TEST_CASE( vector_of_vector_dist )
 }
 
 
-
 BOOST_AUTO_TEST_SUITE_END()
 

src/Vector/vector_dist.hpp

@@ -123,6 +123,27 @@ struct gcl<dim,St,CellL,Vector,GCL_SYMMETRIC>
 	}
 };
 
+/////////////////// GCL anisotropic ///////////////////////////////////////////
+
+//! General function t get a cell-list
+template<unsigned int dim, typename St, typename CellL, typename Vector, unsigned int impl>
+struct gcl_An
+{
+	/*! \brief Get the Cell list based on the type
+	 *
+	 * \param vd Distributed vector
+	 * \param r_cut Cut-off radius
+	 * \param g Ghost
+	 *
+	 * \return the constructed cell-list
+	 *
+	 */
+	static inline CellL get(Vector & vd, const size_t (& div)[dim], const size_t (& pad)[dim], const Ghost<dim,St> & g)
+	{
+		return vd.template getCellListSym<CellL>(div,pad);
+	}
+};
+
 #define CELL_MEMFAST(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_fast<>, shift<dim, St> >
 #define CELL_MEMBAL(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_bal<>, shift<dim, St> >
 #define CELL_MEMMW(dim,St) CellList_gen<dim, St, Process_keys_lin, Mem_mw<>, shift<dim, St> >
@@ -996,6 +1017,52 @@ public:
 		return cell_list;
 	}
 
+	/*! \brief Construct a cell list symmetric based on a cut of radius
+	 *
+	 * \tparam CellL CellList type to construct
+	 *
+	 * \param r_cut interation radius, or size of each cell
+	 *
+	 * \return the Cell list
+	 *
+	 */
+	template<typename CellL = CellList<dim, St, Mem_fast<>, shift<dim, St> > >
+	CellL getCellListSym(const size_t (& div)[dim],
+						 const size_t (& pad)[dim])
+	{
+#ifdef SE_CLASS1
+		if (!(opt & BIND_DEC_TO_GHOST))
+		{
+			if (getDecomposition().getGhost().getLow(dim-1) == 0.0)
+			{
+				std::cerr << __FILE__ << ":" << __LINE__ << " Error the vector has been constructed without BIND_DEC_TO_GHOST, If you construct a vector without BIND_DEC_TO_GHOST the ghost must be full without reductions " << std::endl;
+				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
+			}
+		}
+#endif
+
+		size_t pad_max = pad[0];
+		for (size_t i = 1 ; i < dim ; i++)
+		{if (pad[i] > pad_max)	{pad_max = pad[i];}}
+
+		// Cell list
+		CellL cell_list;
+
+		CellDecomposer_sm<dim,St,shift<dim,St>> cd_sm;
+		cd_sm.setDimensions(getDecomposition().getDomain(),div,pad_max);
+
+		// Processor bounding box
+		Box<dim, St> pbox = getDecomposition().getProcessorBounds();
+
+		// Ghost padding extension
+		Ghost<dim,size_t> g_ext(0);
+		cell_list.Initialize(cd_sm,pbox,pad_max);
+		cell_list.set_ndec(getDecomposition().get_ndec());
+
+		updateCellListSym(cell_list);
+
+		return cell_list;
+	}
 
 	/*! \brief Construct a cell list starting from the stored particles
 	 *
@@ -1070,7 +1137,7 @@ public:
 		// but in the worst case we take the maximum
 		St r_cut = 0;
 		for (size_t i = 0 ; i < dim ; i++)
-			r_cut = std::max(r_cut,cell_list.getCellBox().getHigh(i));
+		{r_cut = std::max(r_cut,cell_list.getCellBox().getHigh(i));}
 
 		// Here we have to check that the Cell-list has been constructed
 		// from the same decomposition
@@ -1097,18 +1164,13 @@ public:
 	 * \param cell_list Cell list to update
 	 *
 	 */
-	template<typename CellL = CellList<dim, St, Mem_fast<>, shift<dim, St> > > void updateCellListSym(CellL & cell_list)
+	template<typename CellL = CellList<dim, St, Mem_fast<>, shift<dim, St> > >
+	void updateCellListSym(CellL & cell_list)
 	{
 #ifdef SE_CLASS3
 		se3.getNN();
 #endif
 
-		// This function assume equal spacing in all directions
-		// but in the worst case we take the maximum
-		St r_cut = 0;
-		for (size_t i = 0 ; i < dim ; i++)
-		{r_cut = std::max(r_cut,cell_list.getCellBox().getHigh(i));}
-
 		// Here we have to check that the Cell-list has been constructed
 		// from the same decomposition
 		bool to_reconstruct = cell_list.get_ndec() != getDecomposition().get_ndec();
@@ -1121,7 +1183,10 @@ public:
 		}
 		else
 		{
-			CellL cli_tmp = gcl<dim,St,CellL,self,GCL_SYMMETRIC>::get(*this,r_cut,getDecomposition().getGhost());
+			CellL cli_tmp = gcl_An<dim,St,CellL,self,GCL_SYMMETRIC>::get(*this,
+																		 cell_list.getDivWP(),
+																		 cell_list.getPadding(),
+																		 getDecomposition().getGhost());
 
 			cell_list.swap(cli_tmp);
 		}

src/Vector/vector_dist_comm.hpp

@@ -268,26 +268,40 @@ class vector_dist_comm
 	 *
 	 * \param v_pos vector of particle positions
 	 * \param v_prp vector of particles properties
+	 * \param opt options
 	 *
 	 */
 	void local_ghost_from_opart(openfpm::vector<Point<dim, St>> & v_pos,
-			                    openfpm::vector<prop,Memory,typename layout_base<prop>::type,layout_base> & v_prp)
+			                    openfpm::vector<prop,Memory,typename layout_base<prop>::type,layout_base> & v_prp,
+								size_t opt)
 	{
 		// get the shift vectors
 		const openfpm::vector<Point<dim, St>> & shifts = dec.getShiftVectors();
 
-		for (size_t i = 0 ; i < o_part_loc.size() ; i++)
+		if (!(opt & NO_POSITION))
 		{
-			size_t lin_id = o_part_loc.get<1>(i);
-			size_t key = o_part_loc.template get<0>(i);
+			for (size_t i = 0 ; i < o_part_loc.size() ; i++)
+			{
+				size_t lin_id = o_part_loc.get<1>(i);
+				size_t key = o_part_loc.template get<0>(i);
 
-			Point<dim, St> p = v_pos.get(key);
-			// shift
-			p -= shifts.get(lin_id);
+				Point<dim, St> p = v_pos.get(key);
+				// shift
+				p -= shifts.get(lin_id);
+
+				// add this particle shifting its position
+				v_pos.add(p);
+				v_prp.get(lg_m+i) = v_prp.get(key);
+			}
+		}
+		else
+		{