diff --git a/src/Particles/CMakeLists.txt b/src/Particles/CMakeLists.txt
index 44fd8ff3..36649ec6 100644
--- a/src/Particles/CMakeLists.txt
+++ b/src/Particles/CMakeLists.txt
@@ -9,6 +9,7 @@ particles/regularParticleIdHandler/regularParticleIdHandler.cpp
 SphereParticles/sphereShape/sphereShape.cpp
 SphereParticles/sphereParticles/sphereParticles.cpp
 SphereParticles/sphereParticles/sphereParticlesKernels.cpp
+Insertion/collisionCheck/collisionCheck.cpp
 Insertion/insertionRegion/insertionRegion.cpp
 Insertion/insertion/insertion.cpp
 Insertion/shapeMixture/shapeMixture.cpp
diff --git a/src/Particles/Insertion/collisionCheck/collisionCheck.cpp b/src/Particles/Insertion/collisionCheck/collisionCheck.cpp
new file mode 100644
index 00000000..b5abe370
--- /dev/null
+++ b/src/Particles/Insertion/collisionCheck/collisionCheck.cpp
@@ -0,0 +1,91 @@
+#include "collisionCheck.hpp"
+
+bool
+pFlow::collisionCheck::build()
+{
+	fill(head_, static_cast<uint32>(-1));
+	fill(next_, static_cast<uint32>(-1));
+
+	for (auto i = 0uL; i < position_.size(); i++)
+	{
+		if(!searchBox_.isInside(position_[i]))
+        {
+            fatalErrorInFunction<<
+            "Point "<< position_[i]<< "is not in search box"<<endl;
+        }
+
+        const auto ind                   = pointIndex(position_[i]);
+		next_[i]                         = head_(ind.x(), ind.y(), ind.z());
+		head_(ind.x(), ind.y(), ind.z()) = static_cast<uint32>(i);
+	}
+
+    return true;
+}
+
+pFlow::collisionCheck::collisionCheck(
+  box                 sBox,
+  real                dx,
+  const realx3Vector& pos,
+  const realVector&   diam
+)
+  : searchBox_(sBox),
+    dx_(dx),
+    nCells_((sBox.maxPoint() - sBox.minPoint()) / dx + realx3(1.0)),
+    position_(pos),
+    diameters_(diam),
+    next_("next", pos.size()),
+    head_("head", nCells_.x(), nCells_.y(), nCells_.z())
+{
+	if(!build())
+    {
+        fatalExit;
+    }
+}
+
+bool
+pFlow::collisionCheck::checkPoint(const realx3& p, const real d) const
+{
+    
+    if(!searchBox_.isInside(p)) return false;
+
+    const auto ind = pointIndex(p);
+    
+    uint32 n = head_(ind.x(), ind.y(), ind.z());
+    while( n != -1)
+    {
+        if( (position_[n]-p).length() - 0.5*(diameters_[n]+d )<0.0 )
+        {
+            return false;
+        }
+        n = next_[n];
+    }
+	return true;
+}
+
+bool
+pFlow::collisionCheck::mapLastAddedParticle()
+{
+    size_t n = position_.size()-1;
+    if( next_.size() != n )
+    {
+        fatalErrorInFunction<<
+        "size mismatch of next and position"<<endl;
+        return false;
+    }
+    next_.push_back(-1);
+    const auto& p = position_[n];
+
+    if(!searchBox_.isInside(p))
+    {
+        fatalErrorInFunction<< 
+        "Point "<< p <<" is not inside the search box"<<endl;
+        return false;
+    }
+
+    const auto ind = pointIndex(p);
+
+    next_[n]       = head_(ind.x(), ind.y(), ind.z());
+    head_(ind.x(), ind.y(), ind.z()) = static_cast<uint32>(n);
+
+	return true;
+}
diff --git a/src/Particles/Insertion/collisionCheck/collisionCheck.hpp b/src/Particles/Insertion/collisionCheck/collisionCheck.hpp
new file mode 100644
index 00000000..f0a6b267
--- /dev/null
+++ b/src/Particles/Insertion/collisionCheck/collisionCheck.hpp
@@ -0,0 +1,54 @@
+
+#ifndef __collisionCheck_hpp__
+#define __collisionCheck_hpp__
+
+#include "Vectors.hpp"
+#include "VectorSingles.hpp"
+#include "box.hpp"
+
+
+namespace pFlow
+{
+
+class collisionCheck
+{
+private:
+
+	box                           searchBox_;
+
+	real                          dx_;
+
+	int32x3                       nCells_;
+
+    const realx3Vector&          position_;
+
+	const realVector&            diameters_;
+    
+	uint32Vector                  next_;
+
+	ViewType3D<uint32, HostSpace> head_;
+
+    int32x3 pointIndex(const realx3& p)const
+	{
+		return int32x3( (p - searchBox_.minPoint())/dx_ );
+	}
+
+    bool build();
+
+public:
+
+	collisionCheck(
+	  box                 sBox,
+	  real                dx,
+	  const realx3Vector& pos,
+	  const realVector&   diam
+	);
+
+    bool checkPoint(const realx3& p, const real d)const;
+
+    bool mapLastAddedParticle();
+};
+
+}
+
+#endif //__collisionCheck_hpp__
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/peakableRegion/PeakableRegion/PeakableRegion.hpp b/src/phasicFlow/structuredData/peakableRegion/PeakableRegion/PeakableRegion.hpp
index 90251e6c..16a12e1b 100644
--- a/src/phasicFlow/structuredData/peakableRegion/PeakableRegion/PeakableRegion.hpp
+++ b/src/phasicFlow/structuredData/peakableRegion/PeakableRegion/PeakableRegion.hpp
@@ -65,9 +65,19 @@ public:
 
 	//// - Methods
 
-	virtual bool   isInside(const realx3& point) const override;
+	bool   isInside(const realx3& point) const override;
 
-	virtual realx3 peek() const override;
+	realx3 peek() const override;
+
+	const realx3& minPoint()const override
+	{
+		return region_.minPoint();
+	}
+
+	const realx3& maxPoint()const override
+	{
+		return region_.maxPoint();
+	}
 
 	//// - IO operatoins
 
diff --git a/src/phasicFlow/structuredData/peakableRegion/geometricRegion/geometricRegion.hpp b/src/phasicFlow/structuredData/peakableRegion/geometricRegion/geometricRegion.hpp
index f6dbc7f7..b598a36b 100644
--- a/src/phasicFlow/structuredData/peakableRegion/geometricRegion/geometricRegion.hpp
+++ b/src/phasicFlow/structuredData/peakableRegion/geometricRegion/geometricRegion.hpp
@@ -62,6 +62,16 @@ public:
 
 	realx3 peek() const;
 
+	const auto& minPoint()const
+	{
+		return minPoint_;
+	}
+
+	const auto& maxPoint()const
+	{
+		return maxPoint_;
+	}
+
 	//// IO operation
 	bool   read(const dictionary& dict);
 
diff --git a/src/phasicFlow/structuredData/peakableRegion/peakableRegion/peakableRegion.hpp b/src/phasicFlow/structuredData/peakableRegion/peakableRegion/peakableRegion.hpp
index e1449902..14378494 100644
--- a/src/phasicFlow/structuredData/peakableRegion/peakableRegion/peakableRegion.hpp
+++ b/src/phasicFlow/structuredData/peakableRegion/peakableRegion/peakableRegion.hpp
@@ -57,6 +57,10 @@ public:
 
 	virtual realx3 peek() const = 0;
 
+	virtual const realx3& minPoint()const = 0;
+
+	virtual const realx3& maxPoint()const = 0;
+
 	//// - IO operatoins
 
 	virtual bool   read(const dictionary& dict) = 0;
diff --git a/utilities/particlesPhasicFlow/CMakeLists.txt b/utilities/particlesPhasicFlow/CMakeLists.txt
index 649d9ba0..1ee6f3f8 100644
--- a/utilities/particlesPhasicFlow/CMakeLists.txt
+++ b/utilities/particlesPhasicFlow/CMakeLists.txt
@@ -3,7 +3,7 @@ set(source_files
 particlesPhasicFlow.cpp
 positionParticles/positionParticles.cpp 
 positionOrdered/positionOrdered.cpp 
-#positionRandom/positionRandom.cpp
+positionRandom/positionRandom.cpp
 empty/empty.cpp
 )
 #set(link_lib phasicFlow Kokkos::kokkos Interaction Utilities)
diff --git a/utilities/particlesPhasicFlow/empty/empty.cpp b/utilities/particlesPhasicFlow/empty/empty.cpp
index 5ed799a2..d0da9b3b 100755
--- a/utilities/particlesPhasicFlow/empty/empty.cpp
+++ b/utilities/particlesPhasicFlow/empty/empty.cpp
@@ -30,7 +30,7 @@ pFlow::empty::empty(
 	positionParticles(control, dict),
 	position_
 	(
-		"empty",maxNumberOfParticles_, 0, RESERVE()
+		"empty",maxNumberOfParticles(), 0, RESERVE()
 	)
 {
 }
\ No newline at end of file
diff --git a/utilities/particlesPhasicFlow/empty/empty.hpp b/utilities/particlesPhasicFlow/empty/empty.hpp
index a4583b7c..9b296893 100755
--- a/utilities/particlesPhasicFlow/empty/empty.hpp
+++ b/utilities/particlesPhasicFlow/empty/empty.hpp
@@ -22,7 +22,6 @@ Licence:
 #define __empty_hpp__
 
 #include "positionParticles.hpp"
-#include "box.hpp"
 
 namespace pFlow
 {
@@ -32,10 +31,7 @@ class empty
 :
 	public positionParticles
 {
-protected:
-
-	dictionary 	emptyDict_;
-
+private:
 
 	realx3Vector 	position_;
 
@@ -54,33 +50,33 @@ public:
 		empty,
 		dictionary);
 
-	virtual ~empty() = default;
+	~empty() final = default;
 
 	//// - Methods 
 
-	virtual uint64 numPoints()const
+	uint32 numPoints()const final
 	{
 		return 0;
 	}
 
-	virtual uint64 size()const
+	uint32 size()const final
 	{
 		return 0;
 	}
 
-	real maxDiameter() const override
+	real maxDiameter() const final
 	{
 		return 1.0;
 	}
 
 	// - const access to position
-	virtual const realx3Vector& position()const 
+	const realx3Vector& position()const final 
 	{
 		return position_;
 	}
 
 	// - access to position 
-	virtual realx3Vector& position()
+	realx3Vector& position() final
 	{
 		return position_;
 	}
diff --git a/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.cpp b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.cpp
index 77d624b4..1fb60b87 100755
--- a/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.cpp
+++ b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.cpp
@@ -18,11 +18,9 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-#include "positionOrdered.hpp"
 #include "error.hpp"
-
-
-#include "streams.hpp"
+#include "dictionary.hpp"
+#include "positionOrdered.hpp"
 
 
 bool pFlow::positionOrdered::findAxisIndex()
@@ -45,9 +43,9 @@ bool pFlow::positionOrdered::findAxisIndex()
 		return false;
 	}
 
-	realx3 uV[3];
+	std::array<realx3,3> uV;
 	size_t i=0;
-	for(auto& ca: axisOrder_)
+	for(const auto& ca: axisOrder_)
 	{
 		if(ca == "x")
 		{
@@ -82,15 +80,16 @@ bool pFlow::positionOrdered::positionPointsOrdered()
 	position_.clear();
 
 	realx3 dl(diameter_);
-	auto minP =  region_->minPoint();
-	auto maxP =  region_->maxPoint();
+	const auto& region = pRegion();
+	auto minP =  region.minPoint();
+	auto maxP =  region.maxPoint();
 
 	auto cntr = minP;
 
 	size_t n = 0;
 	while( n < numPoints_ )
 	{
-		if(region_->isInside(cntr))
+		if(region.isInside(cntr))
 		{
 			position_.push_back(cntr);
 			n++;
@@ -130,7 +129,7 @@ pFlow::positionOrdered::positionOrdered
 	positionParticles(control, dict),
 	poDict_
 	(  
-		dict.subDict("positionOrderedInfo")
+		dict.subDict("orderedInfo")
 	),
 	diameter_
 	(
@@ -146,7 +145,10 @@ pFlow::positionOrdered::positionOrdered
 	),
 	position_
 	(
-		"positionOrdered", maxNumberOfParticles_, numPoints_ ,RESERVE()
+		"positionOrdered", 
+		max(maxNumberOfParticles(), numPoints_), 
+		numPoints_ ,
+		RESERVE()
 	)
 {
 	
@@ -155,13 +157,6 @@ pFlow::positionOrdered::positionOrdered
 		fatalExit;
 	}
 
-	if(!region_)
-	{
-		fatalErrorInFunction<<"You must provided a region (box, cylinder, ...) for positioning particles in dictionary "<<
-		dict.globalName()<<endl;
-		fatalExit;
-	}
-
 	if(!positionPointsOrdered())
 	{
 		fatalExit;
diff --git a/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.hpp b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.hpp
index 5d18ae57..caf711dd 100755
--- a/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.hpp
+++ b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.hpp
@@ -31,13 +31,13 @@ class positionOrdered
 :
 	public positionParticles
 {
-protected:
+private:
 
 	dictionary 	poDict_;
 
 	real 		diameter_;
 
-	size_t 		numPoints_;
+	uint32 		numPoints_;
 	
 	wordList 	axisOrder_;
 
@@ -60,7 +60,7 @@ protected:
 public:
 
 	// - type Info
-	TypeInfo("positionOrdered");
+	TypeInfo("ordered");
 
 	positionOrdered(
 		systemControl& control,
@@ -72,33 +72,33 @@ public:
 		positionOrdered,
 		dictionary);
 
-	virtual ~positionOrdered() = default;
+	~positionOrdered() final = default;
 
 	//// - Methods 
 
-	virtual uint64 numPoints()const
+	uint32 numPoints()const final
 	{
-		return position_.size();
+		return static_cast<uint32>(position_.size());
 	}
 
-	virtual uint64 size()const
+	uint32 size()const final
 	{
-		return position_.size();
+		return static_cast<uint32>(position_.size());
 	}
 
-	real maxDiameter() const override
+	real maxDiameter() const final
 	{
 		return diameter_;
 	}
 	
 	// - const access to position
-	virtual const realx3Vector& position()const 
+	const realx3Vector& position()const final 
 	{
 		return position_;
 	}
 
 	// - access to position 
-	virtual realx3Vector& position()
+	realx3Vector& position() final
 	{
 		return position_;
 	}
diff --git a/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp b/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp
index a7d3879a..d99cb217 100755
--- a/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp
+++ b/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp
@@ -20,16 +20,11 @@ Licence:
 
 #include "positionParticles.hpp"
 #include "vocabs.hpp"
-#include "box.hpp"
-#include "cylinder.hpp"
-#include "sphere.hpp"
-#include "cells.hpp"
-//#include "contactSearchFunctions.hpp"
+#include "dictionary.hpp"
+#include "systemControl.hpp"
 
-#include "streams.hpp"
-
-
-pFlow::realx3Vector pFlow::positionParticles::sortByMortonCode(realx3Vector& position)const
+pFlow::realx3Vector pFlow::positionParticles::sortByMortonCode(
+	const realx3Vector& position)const
 {
 	struct indexMorton
 	{
@@ -81,25 +76,20 @@ pFlow::positionParticles::positionParticles
 	systemControl& control,
 	const dictionary& dict
 )
+:
+	regionType_(dict.getValOrSet<word>("regionType", "domain")),
+	maxNumberOfParticles_(dict.getValOrSet(
+		"maxNumberOfParticles", 
+		static_cast<uint32>(10000))),
+	mortonSorting_(dict.getValOrSet("mortonSorting", Logical("Yes")))
 {
-	maxNumberOfParticles_ = dict.getValOrSet("maxNumberOfParticles", static_cast<uint64>(10000));
-	
-	mortonSorting_ = dict.getValOrSet("mortonSorting", Logical("Yes"));
 
-	if( dict.containsDictionay("box") )	
+	if( regionType_ != "domain" )	
 	{
-		region_ = makeUnique<region<box>>(dict.subDict("box"));
-	}
-	else if(dict.containsDictionay("cylinder"))
-	{
-        WARNING<<"cylinder region is not set!"<<END_WARNING;
-		//region_ = makeUnique<region<cylinder>>(dict.subDict("cylinder"));
-	}
-	else if(dict.containsDictionay("sphere"))
-	{
-        WARNING<<"sphere region is not set!"<<END_WARNING;
-		//region_ = makeUnique<region<sphere>>(dict.subDict("sphere"));
-	}
+		pRegion_ = peakableRegion::create(
+			regionType_,
+			dict.subDict(regionType_+"Info"));
+	}	
 	else
 	{
 		fileDictionary domainDict
@@ -113,7 +103,7 @@ pFlow::positionParticles::positionParticles
 			},
 			&control.settings()
 		);
-		region_ = makeUnique<region<box>>( domainDict.subDict("globalBox"));
+		pRegion_ = peakableRegion::create(regionType_,domainDict.subDict("globalBox"));
 	}
 	
 }
@@ -128,9 +118,8 @@ pFlow::realx3Vector pFlow::positionParticles::getFinalPosition()
 	else
 	{
 		realx3Vector vec("final position",position().capacity(), 0 , RESERVE());
-		vec.assign( position().begin(), position().end());
-		
-		return std::move(vec);
+		vec.assign( position().begin(), position().end());	
+		return vec;
 	}
 }
 
diff --git a/utilities/particlesPhasicFlow/positionParticles/positionParticles.hpp b/utilities/particlesPhasicFlow/positionParticles/positionParticles.hpp
index ddd12455..85c34e54 100755
--- a/utilities/particlesPhasicFlow/positionParticles/positionParticles.hpp
+++ b/utilities/particlesPhasicFlow/positionParticles/positionParticles.hpp
@@ -23,96 +23,39 @@ Licence:
 
 #include "virtualConstructor.hpp"
 #include "Vectors.hpp"
-#include "dictionary.hpp"
-#include "systemControl.hpp"
+#include "peakableRegion.hpp"
 
 namespace pFlow
 {
 
-class regionBase
-{
-public:
-	
-	regionBase() = default;
-	
-	regionBase(const regionBase&) = default;
+class dictionary;
+class systemControl;
 
-	regionBase& operator =(const regionBase&) = default;
-	
-	virtual ~regionBase() = default;
-
-	virtual bool isInside(const realx3 point)const = 0;
-
-	virtual realx3 minPoint()const =0;
-
-	virtual realx3 maxPoint()const =0;
-
-	virtual word name()const =0;
-
-};
-
-template<typename T>
-class region
-:
-	public regionBase
-{
-	protected:
-
-		T 		region_;
-
-	public:
-		
-		region(const T& rgn)
-		:
-			region_(rgn)
-		{}
-
-		region(const dictionary& dict)
-		:
-			region_(dict)
-		{}
-
-		region(const region&) = default;
-		
-		region& operator =(const region&) = default;
-
-		virtual ~region()=default;
-
-		bool isInside(const realx3 point) const override
-		{
-			return region_.isInside(point);
-		}
-
-		realx3 minPoint()const override
-		{
-			return region_.minPoint();
-		}
-
-		realx3 maxPoint()const override
-		{
-			return region_.maxPoint();
-		}
-
-		word name()const override
-		{
-			return region_.typeName();
-		}
-
-};
 
 class positionParticles
 {
-protected:
+private:
 
-	uniquePtr<regionBase> 		region_ = nullptr;
+	uniquePtr<peakableRegion> 		pRegion_ = nullptr;
 
-	size_t  		maxNumberOfParticles_ = 10000;
+	word 			regionType_;
+
+	uint32  		maxNumberOfParticles_ = 10000;
 
 	Logical 		mortonSorting_;
 
-	static const size_t numReports_ = 40;
+	
 
-	realx3Vector sortByMortonCode(realx3Vector& position)const;
+	realx3Vector sortByMortonCode(const realx3Vector& position)const;
+
+protected:
+
+	static const uint32 numReports_ = 40;
+
+	const auto& pRegion()const
+	{
+		return pRegion_();
+	}
 	
 public:
 
@@ -133,11 +76,22 @@ public:
 
 	virtual ~positionParticles() = default;
 
-	//// - Methods 
+	//// - Methods
 
-	virtual uint64 numPoints()const = 0;
+	bool mortonSorting()const
+	{
+		return mortonSorting_();
+	} 
 
-	virtual uint64 size()const = 0;
+	inline
+	auto maxNumberOfParticles()const
+	{
+		return maxNumberOfParticles_;
+	}
+
+	virtual uint32 numPoints()const = 0;
+
+	virtual uint32 size()const = 0;
 
 	virtual real maxDiameter() const = 0;
 
diff --git a/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp b/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp
index 750509f3..adb387fe 100755
--- a/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp
+++ b/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp
@@ -20,141 +20,81 @@ Licence:
 
 
 #include "positionRandom.hpp"
-#include "uniformRandomReal.hpp"
-#include "NBSLevel0.hpp"
-#include "unsortedPairs.hpp"
-#include "box.hpp"
+#include "collisionCheck.hpp"
 
 
 
-namespace pFlow
+bool pFlow::positionRandom::positionOnePass(collisionCheck& collCheck)
 {
 
-using SearchType = NBSLevel0<DefaultExecutionSpace> ;
-using ContainerType = unsortedPairs<DefaultExecutionSpace, int32>;
+	auto const& region = pRegion();
+	
+	auto n = static_cast<uint32>(position_.size());
 
-
-
-int32 findCollisions(
-		ContainerType& pairs,
-		int32Vector_HD& flags);
-
-int32 findCollisions(int32 num, realx3Vector_HD& points, real diam)
-{
-	int32 res =0;
-	for(auto i=0; i<num;i++)
+	for(uint32 iter = 0; iter<numPoints_; iter++)
 	{
-		for(auto j=i+1; j<num; j++)
+		realx3 p = region.peek();
+
+		if( collCheck.checkPoint(p, diameter_) )
 		{
-			if(sphereSphereCheck(points[i],points[j],diam,diam))res++;
-		}
-	}
-
-	return res;
-}
-
-}
-
-
-bool pFlow::positionRandom::positionOnePass(int32 pass, int32 startNum)
-{
-
-	realVector_D diameter(startNum , diameter_);
-	int32Vector_HD  flagHD(startNum, 0);
-	realx3Vector_HD positionHD(startNum);
-	
-	auto minP = region_->minPoint();
-	auto maxP = region_->maxPoint();
-
-	SearchType search(
-					box(minP, maxP),
-					diameter_,
-					positionHD.deviceViewAll(),
-					diameter.deviceViewAll());
-
-	ContainerType pairs(3*startNum);
-		
-	REPORT(1)<< "Positioning "<<
-		greenText("(Pass #"<< pass+1<<")")<< 
-		": started with "<< startNum <<" points."<<endREPORT;
-	
-	fillPoints(startNum, positionHD, flagHD);
-		
-	search.broadSearch(pairs, range(0, startNum));
-
-	
-	int32 numCollisions = findCollisions(pairs, flagHD);
-	
-	
-	REPORT(2)<< "Positioned " << cyanText(startNum - numCollisions) <<
-	" without collision \n"<<endREPORT; 
-
-	if(startNum-numCollisions >= numPoints_ )
-	{
-		
-		REPORT(1)<<"Selected "<< cyanText(numPoints_)<< " for the final field.\n"<<endREPORT;
-
-		positionHD.syncViews();
-		position_.clear();
-		int32 n=0;
-		for(int32 i=0; i<startNum; i++)
-		{
-			if(flagHD[i] == 0 )
-			{
-				position_.push_back( positionHD[i]);
-				n++;
-				if(n==numPoints_)break;	
-			}
+			position_.push_back(p);
+			diameters_.push_back(diameter_);
 			
+			if(!collCheck.mapLastAddedParticle())
+			{
+				fatalErrorInFunction;
+				return false;
+			}
+			n++;
+
+			if(n == numPoints_) break;
 		}
-		
-		return true;
 	}
 	
-	
-	return false;
-	
-	
+	return true;	
 }
 
 bool pFlow::positionRandom::positionPointsRandom()
 {
 	
 	position_.clear();
-
+	diameters_.clear();
 	if(numPoints_ == 0)return true;
 	
-	size_t pass = 0;
-	int32 startNum = numPoints_;
+	uint32 pass = 0;
+	collisionCheck collCheck(
+		box(pRegion().minPoint(), pRegion().maxPoint()),
+		diameter_,
+		position_,
+		diameters_);
+	
 
 	while ( pass <maxIterations_)
 	{
-		if( positionOnePass(pass, startNum) )return true;
-		startNum = 1.1*startNum+1;
+		if( !positionOnePass(collCheck) )
+		{
+			return false;
+		}
 		pass++;
+		REPORT(1)<<"Positioning "<< Green_Text("(Pass #"<< pass<<")")<<
+		": number of non-colliding spheres is "<<
+		Yellow_Text(position_.size())<<END_REPORT;
+
+		if( position_.size() == numPoints_ ) 
+		{
+			REPORT(0)<<END_REPORT;
+			return true;
+		}
 	}	
 	
-	
 	fatalErrorInFunction<< 
-	"  cannot position "<< numPoints_ << " in the domain in " << maxIterations_ << " iterations.\n" <<
+	"  Cannot position "<< numPoints_ << " in the domain in " <<
+	maxIterations_ << " iterations.\n" <<
 	"  you may increase maxIterations for positioning points.\n";
 
 	return false;
 }
 
-bool pFlow::positionRandom::inCollision
-(
-	const realx3 &cntr,
-	real diam
-)
-{
-	for(const auto& cp: position_)
-	{
-		if( length(cp-cntr) <= diam ) return true;
-	}
-
-	return false;
-}
 
 pFlow::positionRandom::positionRandom
 (
@@ -165,7 +105,7 @@ pFlow::positionRandom::positionRandom
 	positionParticles(control, dict),
 	prDict_
 	(  
-		dict.subDict("positionRandomInfo")
+		dict.subDict("randomInfo")
 	),
 	diameter_
 	(
@@ -173,107 +113,35 @@ pFlow::positionRandom::positionRandom
 	),
 	numPoints_
 	(
-		prDict_.getVal<size_t>("numPoints")
+		prDict_.getVal<uint32>("numPoints")
 	),
 	maxIterations_
 	(
-		prDict_.getValOrSet("maxIterations", 10)
+		prDict_.getValOrSet("maxIterations", 10u)
 	),
 	position_
 	(
-		maxNumberOfParticles_, RESERVE()
+		"position",
+		maxNumberOfParticles(),
+		0,
+		RESERVE()
+	),
+	diameters_
+	(
+		"diameters",
+		maxNumberOfParticles(),
+		0,
+		RESERVE()
 	)
 {
 
-	reportInterval_ = max(numPoints_/numReports_, static_cast<size_t>(2));
+	reportInterval_ = max(numPoints_/numReports_, static_cast<uint32>(2));
 
 	if( !positionPointsRandom() )
 	{
 		fatalExit;
 	}
 
-	if(!region_)
-	{
-		fatalErrorInFunction<<"You must provided a region (box, cylinder, ...) for positioning particles in dictionary "<<
-		dict.globalName()<<endl;
-		fatalExit;
-	}
-
-}
-
-void pFlow::positionRandom::fillPoints(
-	uint numPoints,
-	realx3Vector_HD& points,
-	int32Vector_HD& flags )
-{
-
-	uniformRandomReal rand;
-	
-	auto minP = region_().minPoint();
-	auto maxP = region_().maxPoint();
-
-	for(size_t i=0; i<numPoints; i++)
-	{
-		if(flags[i] == 0)
-		{
-			bool loop=true;
-			size_t n=0;
-			while (loop)
-			{
-
-				auto pos = rand(minP, maxP);
-				if( region_().isInside(pos))
-				{
-					points[i] =pos;
-					loop = false;		
-				}
-				n++;
-
-				if(n>100)
-				{
-					fatalErrorInFunction<< 
-					"could not find a point inside region"<<region_->name()<<endl;
-					fatalExit;
-				}
-				
-			}
-			
-		}
-	}
-	points.modifyOnHost();
-	points.syncViews();
-}
-
-pFlow::int32 pFlow::findCollisions(
-		ContainerType& pairs,
-		int32Vector_HD& flags)
-{
-	auto allPairs = pairs.getPairs();
-	auto num = pairs.capacity();
-	auto dFlags = flags.deviceView();
-	
-
-	int32 numCollisions = 0;
-
-	Kokkos::parallel_reduce(
-		"positionRandom::findCollisions",
-		num,
-		LAMBDA_HD(int32 i, int32& valueToUpdate){
-			if(allPairs.isValid(i))
-			{
-				auto pair = allPairs.getPair(i);
-				if( dFlags[pair.first] ==0 )
-				{
-					dFlags[pair.first] = 1;
-					valueToUpdate++;	
-				}
-			}
-		}, numCollisions);
-
-	flags.modifyOnDevice();
-	flags.syncViews();	
-
-	return numCollisions;
 }
 
 
diff --git a/utilities/particlesPhasicFlow/positionRandom/positionRandom.hpp b/utilities/particlesPhasicFlow/positionRandom/positionRandom.hpp
index ff301c54..6faa81d9 100755
--- a/utilities/particlesPhasicFlow/positionRandom/positionRandom.hpp
+++ b/utilities/particlesPhasicFlow/positionRandom/positionRandom.hpp
@@ -2,90 +2,79 @@
       O        C enter of
      O O       E ngineering and
     O   O      M ultiscale modeling of
-   OOOOOOO     F luid flow       
+   OOOOOOO     F luid flow
 ------------------------------------------------------------------------------
   Copyright (C): www.cemf.ir
   email: hamid.r.norouzi AT gmail.com
-------------------------------------------------------------------------------  
+------------------------------------------------------------------------------
 Licence:
-  This file is part of phasicFlow code. It is a free software for simulating 
+  This file is part of phasicFlow code. It is a free software for simulating
   granular and multiphase flows. You can redistribute it and/or modify it under
-  the terms of GNU General Public License v3 or any other later versions. 
- 
-  phasicFlow is distributed to help others in their research in the field of 
+  the terms of GNU General Public License v3 or any other later versions.
+
+  phasicFlow is distributed to help others in their research in the field of
   granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
-
-#ifndef __positionOrdered_hpp__
-#define __positionOrdered_hpp__
+#ifndef __positionRandom_hpp__
+#define __positionRandom_hpp__
 
 #include "positionParticles.hpp"
-#include "VectorSingles.hpp"
-#include "VectorDuals.hpp"
-
+#include "dictionary.hpp"
 
 namespace pFlow
 {
 
+class collisionCheck;
 
-
-class positionRandom
-:
-	public positionParticles
+class positionRandom : public positionParticles
 {
-protected:
+private:
 
-	dictionary 	prDict_;
+	dictionary   prDict_;
 
-	real 		diameter_;
+	real         diameter_;
 
-	size_t 		numPoints_;
+	uint32       numPoints_;
 
-	size_t 		maxIterations_;
-	
-	realx3Vector 	position_;
+	uint32       maxIterations_;
 
-	size_t 		reportInterval_;
+	realx3Vector position_;
 
-	bool positionOnePass(int32 pass, int32 startNum);
-	
-	bool positionPointsRandom();
+	realVector   diameters_;
 
-	bool inCollision(const realx3 &cntr, real diam);
+	uint32       reportInterval_;
 
-	void fillPoints(
-		uint numPoints,
-		realx3Vector_HD& points,
-		int32Vector_HD& flags );
+	bool         positionOnePass(collisionCheck& collCheck);
+
+	bool         positionPointsRandom();
 
 public:
 
 	// - type Info
-	TypeInfo("positionRandom");
+	TypeInfo("random");
 
-	positionRandom(
-		systemControl& control,
-		const dictionary& dict);
+	positionRandom(systemControl& control, const dictionary& dict);
 
-	// - add this class to vCtor selection table 
-	add_vCtor(
-		positionParticles,
-		positionRandom,
-		dictionary);
+	// - add this class to vCtor selection table
+	add_vCtor
+	(
+		positionParticles, 
+		positionRandom, 
+		dictionary
+	);
 
-	virtual ~positionRandom() = default;
+	~positionRandom() final = default;
 
-	//// - Methods 
+	//// - Methods
 
-	virtual uint64 numPoints()const
+	uint32 numPoints() const final
 	{
 		return position_.size();
 	}
 
-	virtual uint64 size()const
+	uint32 size() const final
 	{
 		return position_.size();
 	}
@@ -96,25 +85,18 @@ public:
 	}
 
 	// - const access to position
-	virtual const realx3Vector& position()const 
+	const realx3Vector& position() const final
 	{
 		return position_;
 	}
 
-	// - access to position 
-	virtual realx3Vector& position()
+	// - access to position
+	realx3Vector& position() final
 	{
 		return position_;
 	}
-
-
-
-	
-};	
-
+};
 
 }
 
-
-
 #endif // __positionOrdered_hpp__