src folder

src/Interaction/contactSearch/ContactSearch/ContactSearch.H

@@ -0,0 +1,250 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   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 
+  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 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+#ifndef __ContactSearch_H__
+#define __ContactSearch_H__
+
+
+#include "contactSearch.H"
+#include "box.H"
+
+namespace pFlow
+{
+
+template<
+	template<class, class, class> class BaseMethod,
+	template<class, class, class> class WallMapping	
+>
+class ContactSearch
+:
+	public contactSearch
+{
+public:
+	
+	using IdType 			= typename contactSearch::IdType;
+
+	using IndexType 		= typename contactSearch::IndexType;
+
+	using ExecutionSpace 	= typename contactSearch::ExecutionSpace;
+
+	using PairContainerType = typename contactSearch::PairContainerType;
+
+	using ParticleContactSearchType = 
+		BaseMethod<
+			ExecutionSpace,
+			IdType,
+			IndexType
+			>;
+
+	using WallMappingType =
+		WallMapping<
+			ExecutionSpace,
+			IdType,
+			IndexType
+			>;
+
+protected:
+	
+
+	uniquePtr<ParticleContactSearchType> particleContactSearch_ = nullptr;
+
+	uniquePtr<WallMappingType> 			 wallMapping_ = nullptr;
+
+public:
+
+	TypeNameTemplate2("ContactSearch", ParticleContactSearchType, WallMappingType);
+
+	ContactSearch(
+		const dictionary& csDict,
+		const box& domain,
+	 	const particles& prtcl,
+	 	const geometry& geom,
+	 	Timers& timers)
+	:
+		contactSearch(csDict, domain, prtcl, geom, timers)
+		
+	{
+
+		auto method = dict().getVal<word>("method");
+		auto wmMethod = dict().getVal<word>("wallMapping");
+		
+		auto nbDict = dict().subDictOrCreate(method+"Info");
+
+		real minD, maxD;
+		this->Particles().boundingSphereMinMax(minD, maxD);
+		
+		const auto& position = this->Particles().pointPosition().deviceVectorAll();
+		const auto& diam = this->Particles().boundingSphere().deviceVectorAll();
+
+		particleContactSearch_ = 
+			makeUnique<ParticleContactSearchType>
+			(
+				nbDict,
+				this->domain(),
+				maxD,
+				position,
+				diam
+			);
+		Report(2)<<"Contact search algorithm for particle-particle is "<<
+				 greenText(particleContactSearch_().typeName())<<endReport;
+
+
+		auto wmDict = dict().subDictOrCreate(wmMethod+"Info");
+
+		int32 wnPoints = this->Geometry().numPoints();
+		int32 wnTri    = this->Geometry().size();
+
+		const auto& wPoints = this->Geometry().points().deviceVectorAll();
+		const auto& wVertices = this->Geometry().vertices().deviceVectorAll();
+
+		wallMapping_ = 
+			makeUnique<WallMappingType>(
+				wmDict,
+				particleContactSearch_().getCells(),
+				wnPoints,
+				wnTri, 
+				wPoints,
+				wVertices
+				);
+		Report(2)<<"Wall mapping algorithm for particle-wall is "<<
+				 greenText(wallMapping_().typeName())<< endReport;
+
+	}
+
+
+	add_vCtor(
+		contactSearch,
+		ContactSearch,
+		dictionary);
+	
+	bool broadSearch(
+		PairContainerType& ppPairs,
+		PairContainerType& pwPairs,
+		bool force = false) override
+	{
+
+		//output<<" ContactSearch::broadSearch::PP before.\n";
+		if(particleContactSearch_)
+		{
+			auto activeRange = this->Particles().activeRange();
+
+			sphereSphereTimer_.start();
+
+			if(this->Particles().allActive())
+			{
+				particleContactSearch_().broadSearch(ppPairs, activeRange, force);
+			}
+			else
+			{	
+				particleContactSearch_().broadSearch(ppPairs, activeRange, this->Particles().activePointsMaskD(), force);
+			}
+
+			sphereSphereTimer_.end();
+
+		}
+		else
+			return false;
+		
+		if(wallMapping_)
+		{
+			sphereWallTimer_.start();
+			wallMapping_().broadSearch(pwPairs, particleContactSearch_(), force);
+			sphereWallTimer_.end();
+		}
+		else
+			return false;
+		
+		return true;
+	}
+
+
+	bool ppEnterBroadSearch()const override
+	{
+		if(particleContactSearch_)
+		{
+			return particleContactSearch_().enterBoadSearch(); 
+		}
+		return false;
+	}
+
+	bool pwEnterBroadSearch()const override
+	{
+		if(wallMapping_)
+		{
+			return wallMapping_().enterBoadSearch();
+		}
+		return false;
+	}
+	
+
+	bool ppPerformedBroadSearch()const override
+	{
+		if(particleContactSearch_)
+		{
+			return particleContactSearch_().performedSearch();
+		}
+		return false;
+	}
+
+	
+	bool pwPerformedBroadSearch()const override
+	{
+		if(wallMapping_)
+		{
+			return wallMapping_().performedSearch();
+		}
+		return false;
+	}
+
+	/*bool update(const eventMessage& msg) 
+	{
+		if(msg.isSizeChanged() ) 
+		{
+			auto newSize = this->prtcl().size();
+			if(!particleContactSearch_().objectSizeChanged(newSize))
+			{
+				fatalErrorInFunction<<
+				"erro in changing the size for particleContactSearch_ \n";
+				return false;
+			}
+		}
+
+		if(msg.isCapacityChanged() ) 
+		{
+			auto newSize = this->prtcl().capacity();
+			if(!particleContactSearch_().objectSizeChanged(newSize))
+			{
+				fatalErrorInFunction<<
+				"erro in changing the capacity for particleContactSearch_ \n";
+				return false;
+			}
+		}	
+		
+		return true;
+	}*/
+	
+	
+};
+
+}
+
+
+#endif //__ContactSearch_H__

src/Interaction/contactSearch/ContactSearch/ContactSearchs.C

@@ -0,0 +1,26 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   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 
+  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 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#include "ContactSearch.H"
+
+#include "cellsSimple.H"
+#include "NBS.H"
+
+template class pFlow::ContactSearch<pFlow::NBS, pFlow::cellsSimple>;

src/Interaction/contactSearch/cells.H

@@ -0,0 +1,233 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   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 
+  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 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __cells_H__
+#define __cells_H__
+
+
+#include "types.H"
+#include "box.H"
+
+namespace pFlow
+{
+
+template<typename indexType>
+class cells
+{
+public:
+	
+	using CellType = triple<indexType>;
+
+protected:
+
+	// - domain 
+	box 		domain_;
+
+	// - cell size 
+	real 		cellSize_;
+
+	CellType 	numCells_;
+
+	
+	// - protected methods
+	INLINE_FUNCTION_H
+	void calculate()
+	{
+		numCells_ = (domain_.maxPoint()-domain_.minPoint())/cellSize_ + realx3(1.0);
+		numCells_ = max( numCells_ , CellType(static_cast<indexType>(1)) );
+	}
+
+public:
+
+	INLINE_FUNCTION_HD
+	cells()
+	{}
+
+	INLINE_FUNCTION_H
+	cells(const box& domain, real cellSize)
+	:
+		domain_(domain),
+		cellSize_(cellSize)
+	{
+		calculate();
+	}
+
+	INLINE_FUNCTION_HD
+	cells(const cells&) = default;
+
+	INLINE_FUNCTION_HD
+	cells& operator = (const cells&) = default;
+
+	cells getCells()const
+	{
+		return *this;
+	}
+
+	INLINE_FUNCTION_H
+	void setCellSize(real cellSize)
+	{
+		cellSize_ = cellSize;
+		calculate();
+	}
+
+	INLINE_FUNCTION_HD
+	real cellSize()const
+	{
+		return cellSize_;
+	}
+
+	INLINE_FUNCTION_HD
+	const CellType& numCells()const
+	{
+		return numCells_;
+	}
+
+	INLINE_FUNCTION_HD
+	indexType nx()const
+	{
+		return numCells_.x();
+	}
+
+	INLINE_FUNCTION_HD
+	indexType ny()const
+	{
+		return numCells_.y();
+	}
+
+	INLINE_FUNCTION_HD
+	indexType nz()const
+	{
+		return numCells_.z();
+	}
+
+	INLINE_FUNCTION_HD
+	int64 totalCells()const
+	{
+		return  static_cast<int64>(numCells_.x())*
+				static_cast<int64>(numCells_.y())*
+				static_cast<int64>(numCells_.z());
+	}
+	
+	INLINE_FUNCTION_HD
+	CellType pointIndex(const realx3& p)const
+	{
+		return CellType( (p - domain_.minPoint())/cellSize_ );
+	}
+
+	INLINE_FUNCTION_HD
+	bool pointIndexInDomain(const realx3 p, CellType& index)const
+	{
+		if( !domain_.isInside(p) ) return false;
+		
+		index = this->pointIndex(p);
+		return true;
+	}
+
+	INLINE_FUNCTION_HD
+	bool inDomain(const realx3& p)const
+	{	
+		return domain_.isInside(p);
+	}
+
+	INLINE_FUNCTION_HD
+	bool isInRange(const CellType& cell)const
+	{
+		if(cell.x()<0)return false;
+		if(cell.y()<0)return false;
+		if(cell.z()<0)return false;
+		if(cell.x()>numCells_.x()-1) return false;
+		if(cell.y()>numCells_.y()-1) return false;
+		if(cell.z()>numCells_.z()-1) return false;
+		return true;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isInRange(indexType i, indexType j, indexType k)const
+	{
+		if(i<0)return false;
+		if(j<0)return false;
+		if(k<0)return false;
+		if(i>numCells_.x()-1) return false;
+		if(j>numCells_.y()-1) return false;
+		if(k>numCells_.z()-1) return false;
+		return true;
+	}
+
+	INLINE_FUNCTION_HD
+	void extendBox(
+		const CellType& p1,
+		const CellType& p2,
+		const CellType& p3,
+		indexType extent,
+		CellType& minP,
+		CellType& maxP)const
+	{
+		minP = min( min( p1, p2), p3)-extent;
+		maxP = max( max( p1, p2), p3)+extent;
+
+		minP = bound(minP);
+		maxP = bound(maxP);
+	}
+
+	INLINE_FUNCTION_HD
+	void extendBox(
+		const realx3& p1,
+		const realx3& p2,
+		const realx3& p3, 
+		real  extent,
+		realx3& minP,
+		realx3& maxP)const
+	{
+		minP = min(min(p1,p2),p3) - extent*cellSize_ ;
+		maxP = max(max(p1,p2),p3) + extent*cellSize_ ;
+
+		minP = bound(minP);
+		maxP = bound(maxP);
+	}
+
+	INLINE_FUNCTION_HD
+	CellType bound(CellType p)const
+	{
+		return CellType(
+			min( numCells_.x()-1, max(0,p.x())),
+			min( numCells_.y()-1, max(0,p.y())),
+			min( numCells_.z()-1, max(0,p.z()))
+			);
+	}
+
+	INLINE_FUNCTION_HD
+	realx3 bound(realx3 p)const
+	{
+		return realx3(
+			min( domain_.maxPoint().x(), max(domain_.minPoint().x(),p.x())),
+			min( domain_.maxPoint().y(), max(domain_.minPoint().y(),p.y())),
+			min( domain_.maxPoint().z(), max(domain_.minPoint().z(),p.z()))
+			);
+	}
+
+
+
+};	
+
+
+}
+
+
+#endif

src/Interaction/contactSearch/contactSearch/contactSearch.C

@@ -0,0 +1,78 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   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 
+  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 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#include "contactSearch.H"
+
+
+
+pFlow::contactSearch::contactSearch(
+	const dictionary& dict,
+	const box& domain,
+ 	const particles& prtcl,
+ 	const geometry& geom,
+ 	Timers& timers)
+:
+	interactionBase(prtcl, geom),
+	domain_(domain),
+	dict_(dict),
+	sphereSphereTimer_("particle-particle contact search", &timers),
+	sphereWallTimer_("particle-wall contact search", &timers)
+{
+
+}
+
+
+pFlow::uniquePtr<pFlow::contactSearch> pFlow::contactSearch::create(
+	const dictionary& dict,
+	const box& domain,
+ 	const particles& prtcl,
+ 	const geometry& geom,
+ 	Timers& timers)
+{
+
+	word baseMethName	= dict.getVal<word>("method");
+	word wallMethod 	= dict.getVal<word>("wallMapping");
+
+	auto model = angleBracketsNames2("ContactSearch", baseMethName, wallMethod);
+
+
+	Report(1)<<"Selecting contact search model . . ."<<endReport;
+	
+
+	if( dictionaryvCtorSelector_.search(model))
+	{
+		auto objPtr = dictionaryvCtorSelector_[model] (dict, domain, prtcl, geom, timers);
+		Report(2)<<"Model "<< greenText(model)<<" is created."<<endReport;
+		return objPtr;
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< model << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			dictionaryvCtorSelector_
+		);
+		fatalExit;
+	}
+
+	return nullptr;
+}

src/Interaction/contactSearch/contactSearch/contactSearch.H

@@ -0,0 +1,135 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   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 
+  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 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+#ifndef __contactSearch_H__
+#define __contactSearch_H__
+
+
+#include "interactionBase.H"
+#include "unsortedPairs.H"
+#include "box.H"
+#include "dictionary.H"
+
+namespace pFlow
+{
+
+
+class contactSearch
+:
+	public interactionBase
+{
+public:
+	using IdType 			= typename interactionBase::IdType;
+
+	using IndexType 		= typename interactionBase::IndexType;
+
+	using ExecutionSpace 	= typename interactionBase::ExecutionSpace;
+
+	using PairContainerType   = unsortedPairs<ExecutionSpace, IdType>;
+
+protected:
+
+	const box& 			domain_;
+
+	dictionary 	dict_;
+
+	Timer 		sphereSphereTimer_;
+
+	Timer 		sphereWallTimer_;
+
+	auto& dict()
+	{
+		return dict_;
+	}
+	
+public:
+
+	TypeName("contactSearch");
+
+	contactSearch(
+		const dictionary& dict,
+		const box& domain,
+	 	const particles& prtcl,
+	 	const geometry& geom,
+	 	Timers& timers);
+
+	virtual ~contactSearch()=default;
+
+
+	create_vCtor
+	(
+		contactSearch,
+		dictionary,
+		(
+			const dictionary& dict,
+			const box& domain,
+		 	const particles& prtcl,
+		 	const geometry&  geom,
+		 	Timers& timers
+	 	),
+	 	(dict, domain, prtcl, geom, timers)
+	);
+
+	const auto& domain()const
+	{
+		return domain_;
+	}
+
+	const auto& dict()const
+	{
+		return dict_;
+	}
+
+
+	virtual 
+	bool broadSearch(
+		PairContainerType& ppPairs,
+		PairContainerType& pwPairs,
+		bool force = false) = 0;
+
+	virtual 
+	bool ppEnterBroadSearch()const = 0;
+
+	virtual 
+	bool pwEnterBroadSearch()const = 0;
+
+	virtual 
+	bool ppPerformedBroadSearch()const = 0;
+
+	virtual 
+	bool pwPerformedBroadSearch()const = 0;
+	
+
+	static 
+	uniquePtr<contactSearch> create(
+		const dictionary& dict,
+		const box& domain,
+	 	const particles& prtcl,
+	 	const geometry& geom,
+	 	Timers& timers);
+
+};
+
+
+}
+
+
+#endif //__ContactSearch_H__

src/Interaction/contactSearch/contactSearchFunctions.H

@@ -0,0 +1,112 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   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 
+  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 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __broadSearchFunctions_H__
+#define __broadSearchFunctions_H__
+
+#include "types.H"
+#include "iBox.H"
+
+namespace pFlow
+{
+
+INLINE_FUNCTION_HD
+uint64_t splitBy3(const uint64_t val){
+	uint64_t x = val;
+	x = (x | x << 32) & 0x1f00000000ffff;
+	x = (x | x << 16) & 0x1f0000ff0000ff;
+	x = (x | x << 8) & 0x100f00f00f00f00f;
+	x = (x | x << 4) & 0x10c30c30c30c30c3;
+	x = (x | x << 2) & 0x1249249249249249;
+	return x;
+}
+
+INLINE_FUNCTION_HD
+uint64_t xyzToMortonCode64(uint64_t x, uint64_t y, uint64_t z)
+{
+	return splitBy3(x) | (splitBy3(y) << 1) | (splitBy3(z) << 2);
+}
+
+
+INLINE_FUNCTION_HD
+uint64_t getThirdBits(uint64_t x)
+{
+	x = x & 0x9249249249249249;
+    x = (x | (x >> 2))  & 0x30c30c30c30c30c3;
+    x = (x | (x >> 4))  & 0xf00f00f00f00f00f;
+    x = (x | (x >> 8))  & 0x00ff0000ff0000ff;
+    x = (x | (x >> 16)) & 0xffff00000000ffff;
+    x = (x | (x >> 32)) & 0x00000000ffffffff;
+    return x;
+	
+}
+
+INLINE_FUNCTION_HD
+void mortonCode64Toxyz(uint64_t morton, uint64_t& x, uint64_t& y, uint64_t& z)
+{
+	x = getThirdBits(morton);
+	y = getThirdBits(morton >> 1);
+	z = getThirdBits(morton >> 2);
+}
+
+template<typename indexType, typename cellIndexType>
+INLINE_FUNCTION_HD
+void indexToCell(const indexType idx, const triple<cellIndexType>& extent, triple<cellIndexType>& cell)
+{
+	indexType nxny = extent.x()*extent.y();
+	cell.z() = static_cast<cellIndexType>(idx / nxny);
+	auto rem = idx % nxny;
+	cell.y() = static_cast<cellIndexType>(rem / extent.x());
+	cell.x() = static_cast<cellIndexType>(rem % extent.x());
+}
+
+template<typename cellIndexType>
+INLINE_FUNCTION_HD
+triple<cellIndexType> boxExtent( const iBox<cellIndexType>& box)
+{
+	return triple<cellIndexType>(
+		box.maxPoint().x() - box.minPoint().x() + 1,
+		box.maxPoint().y() - box.minPoint().y() + 1,
+		box.maxPoint().z() - box.minPoint().z() + 1);
+}
+
+template<typename indexType, typename cellIndexType>
+INLINE_FUNCTION_HD
+void indexToCell(const indexType idx, const iBox<cellIndexType>& box, triple<cellIndexType>& cell)
+{
+	auto extent = boxExtent(box);
+	indexToCell(idx, extent, cell);
+	cell+= box.minPoint();
+}
+
+
+
+
+
+
+INLINE_FUNCTION_HD
+bool sphereSphereCheck(const realx3& p1, const realx3 p2, real d1, real d2)
+{
+	return length(p2-p1) < 0.5*(d2+d1);
+}
+
+}
+
+#endif //__broadSearchFunctions_H__

src/Interaction/contactSearch/methods/NBS.H

@@ -0,0 +1,522 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   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 
+  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 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+#ifndef __NBS_H__
+#define __NBS_H__
+
+#include "cells.H"
+#include "contactSearchFunctions.H"
+
+namespace pFlow
+{
+
+
+template<
+	typename executionSpace,
+	typename idType,
+	typename indexType=int32
+	>
+class NBS
+:
+	public cells<indexType>
+{
+public:
+
+	using IdType 		= idType;
+	
+	using IndexType 	= indexType;
+
+	using Cells 		= cells<IndexType>;
+
+	using CellType 		= typename Cells::CellType;
+
+	using ExecutionSpace= executionSpace;
+
+	using memory_space 	= typename ExecutionSpace::memory_space;
+
+	struct TagBuild{};
+
+	struct TagFindPairs{};
+
+protected:
+
+	int32 		capacity_		= 1;
+
+	real 		sizeRatio_ 		= 1.0;
+
+	int32 		updateFrequency_= 1;
+
+	int32 		currentIter_ 	= 0;
+
+	bool 		performedSearch_ = false;
+
+	ViewType1D<realx3, memory_space> 		pointPosition_;
+
+	ViewType1D<real, memory_space>			diameter_;
+
+	ViewType3D<int32, memory_space> 		head_;
+
+	ViewType1D<int32, memory_space> 		next_;
+
+
+	INLINE_FUNCTION_H
+	void nullify()
+	{
+		fill(
+			head_,
+			range(0,this->nx()),
+			range(0,this->ny()),
+			range(0,this->nz()),
+			-1
+			);
+
+		fill(
+			next_,
+			range(0,capacity_),
+			-1
+			);
+	}
+
+	void nullify(range nextRng)
+	{
+		fill(
+			head_,
+			range(0,this->nx()),
+			range(0,this->ny()),
+			range(0,this->nz()),
+			-1
+			);
+
+		fill(
+			next_,
+			nextRng,
+			-1
+			);
+	}
+
+	using mdrPolicyFindPairs = 
+		Kokkos::MDRangePolicy<
+			Kokkos::Rank<3>,
+			Kokkos::Schedule<Kokkos::Dynamic>,
+			ExecutionSpace>;
+
+private:
+
+	void checkAllocateNext(int newCap)
+	{	
+		
+		if( capacity_ < newCap)
+		{
+			capacity_ = newCap;	
+			Kokkos::realloc(next_, capacity_);
+		}
+		
+	}
+
+	void allocateHead()
+	{
+		Kokkos::realloc(head_, this->nx(), this->ny(), this->nz());
+	}
+
+	bool performSearch()
+	{
+		if(currentIter_ % updateFrequency_ == 0)
+		{
+			currentIter_++;
+			return true;
+
+		}else
+		{
+			currentIter_++;
+			return false;
+		}
+	}
+
+	static INLINE_FUNCTION_HD
+	void Swap(int32& x, int32& y)
+	{
+		int32 tmp = x;
+		x = y;
+		y = tmp;
+	}
+
+public:
+
+	TypeNameNV("NBS");
+
+		NBS(
+			const 	box& domain,
+			real 	cellSize,
+			const 	ViewType1D<realx3, memory_space>& position,
+			const 	ViewType1D<real, memory_space>& diam,
+			int32   initialContainerSize = 1)
+		:
+			Cells(domain, cellSize),
+			pointPosition_(position),
+			diameter_(diam),
+			head_("NBS::head_",1,1,1), //, this->nx(), this->ny(), this->nz()),
+			next_("NBS::next_",1) //,position.size()),
+		{
+			checkAllocateNext(pointPosition_.size());
+			allocateHead();
+		}
+
+		NBS(
+			dictionary dict,
+			const 	box& domain,
+			real 	cellSize,
+			const 	ViewType1D<realx3, memory_space>& position,
+			const 	ViewType1D<real, memory_space>& diam,
+			int32   initialContainerSize = 1
+			)
+		:
+			Cells(domain, cellSize),
+			pointPosition_(position),
+			diameter_(diam),
+			head_("NBS::head_",1,1,1), //, this->nx(), this->ny(), this->nz()),
+			next_("NBS::next_",1) //,position.size()),
+		{
+			
+			updateFrequency_ = max( 
+				dict.getVal<int32>("updateFrequency"),1 );
+
+			sizeRatio_ = max(dict.getVal<real>(
+				"sizeRatio"),1.0);
+
+			this->setCellSize(cellSize*sizeRatio_);
+			checkAllocateNext(pointPosition_.size());
+			allocateHead();
+		}
+			
+		INLINE_FUNCTION_HD
+		NBS(const NBS&) = default;
+
+		INLINE_FUNCTION_HD
+		NBS& operator = (const NBS&) = default;
+
+		INLINE_FUNCTION_HD
+		~NBS()=default;
+
+	//// - Methods 
+
+	bool enterBoadSearch()const
+	{
+		return currentIter_%updateFrequency_==0;
+	}
+
+	bool performedSearch()const
+	{
+		return performedSearch_;
+	}
+
+	// - Perform the broad search to find pairs
+	//	 with force = true, perform broad search regardless of 
+	//   updateFrequency_ value
+	//   on all the points in the range of [0,numPoints_)
+	template<typename PairsContainer>
+	bool broadSearch(PairsContainer& pairs, range activeRange, bool force=false)
+	{
+		
+		if(force) currentIter_ = 0;
+		performedSearch_ = false;
+
+		if( !performSearch() ) return true;
+		
+		build(activeRange);
+			
+		findPairs(pairs);
+		
+		performedSearch_ = true;		
+		return true;
+	}
+
+	// - Perform the broad search to find pairs,
+	//   ignore particles with incld(i) = true,
+	//	 with force = true, perform broad search regardless of 
+	//   updateFrequency_ value
+	template<typename PairsContainer, typename IncludeFunction>
+	bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld, bool force = false)
+	{
+		if(force) currentIter_ = 0;
+		performedSearch_ = false;
+
+		if( !performSearch() ) return true;
+
+		build(activeRange, incld);
+	
+		findPairs(pairs);
+		
+		performedSearch_ = true;
+		return true;
+	}
+
+	// - build based on all points in range [0, numPoints_)
+	INLINE_FUNCTION_H
+	void build(range activeRange)
+	{
+		checkAllocateNext(activeRange.second);
+		nullify(activeRange);
+		
+		Cells cellIndex = static_cast<Cells>(*this);
+		auto points = pointPosition_;
+		auto next = next_;
+		auto head = head_;
+
+		Kokkos::RangePolicy<
+			Kokkos::IndexType<int32>,
+			Kokkos::Schedule<Kokkos::Static>,
+			ExecutionSpace> rPolicy(activeRange.first, activeRange.second);
+		
+		Kokkos::parallel_for(
+			"NBS::build",
+			rPolicy,
+			LAMBDA_HD(int32 i){
+				
+				CellType ind = cellIndex.pointIndex(points[i]);	
+				int32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
+				next[i] = old;
+			});
+		Kokkos::fence();
+
+	}
+
+
+	template<typename IncludeFunction>
+	INLINE_FUNCTION_H
+	void build(range activeRange, IncludeFunction incld)
+	{
+		checkAllocateNext(activeRange.second);
+		nullify(activeRange);
+
+		Cells cellIndex = static_cast<Cells>(*this);
+		auto points = pointPosition_;
+		auto next = next_;
+		auto head = head_;
+
+		Kokkos::RangePolicy<
+			Kokkos::IndexType<int32>,
+			Kokkos::Schedule<Kokkos::Dynamic>,
+			ExecutionSpace> rPolicy(activeRange.first, activeRange.second);
+
+		Kokkos::parallel_for(
+			"NBS::build",
+			rPolicy,
+			LAMBDA_HD(int32 i){
+				
+				if( incld(i) )
+				{
+					CellType ind = cellIndex.pointIndex(points[i]);
+					auto old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
+					next[i] = old;
+				}
+			});
+		Kokkos::fence();
+
+	}
+
+	template<typename PairsContainer>	
+	INLINE_FUNCTION_H
+	bool findPairs(PairsContainer& pairs)
+	{
+		mdrPolicyFindPairs
+			mdrPolicy(
+				{0,0,0},
+				{this->nx(),this->ny(),this->nz()} );
+
+		int32 getFull = 1;
+				
+		// loop until the container size fits the numebr of contact pairs
+		while (getFull > 0)
+		{
+
+			getFull = 0;
+			
+			Kokkos::parallel_reduce	(
+				"NBS::broadSearch",
+				mdrPolicy,
+				CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){
+					#include "NBSLoop.H"
+				}, getFull);
+						
+			if(getFull)
+			{
+				// - resize the container
+				//   note that getFull now shows the number of failed insertions.
+				uint32 len = max(getFull,50) ;
+				
+				auto oldCap = pairs.capacity();
+				
+				pairs.increaseCapacityBy(len);
+
+				Info<< "The contact pair container capacity increased from "<<
+				oldCap << " to "<<pairs.capacity()<<" in NBS."<<endInfo;
+				
+			}
+
+			Kokkos::fence();
+		}
+		
+		return true;
+	}
+
+	bool objectSizeChanged(int32 newSize)
+	{
+		checkAllocateNext(newSize);
+		return true;
+	}
+	
+	/*INLINE_FUNCTION_HD
+	void operator()(
+		TagFindPairs,
+		int32 i,
+		int32 j,
+		int32 k,
+		int32& getFullUpdate)const 
+	{
+
+		int32 m = head_(i,j,k);
+		CellType currentCell(i,j,k);
+		int32 n = -1;
+
+		while( m > -1 )
+		{
+			
+			auto p_m = pointPosition_[m];
+			auto d_m = sizeRatio_* diameter_[m];
+			
+			// the same cell
+			n = next_(m);
+			
+			while(n >-1)
+			{
+				
+				auto p_n = pointPosition_[n];
+				auto d_n = sizeRatio_*diameter_[n];
+				
+				if( sphereSphereCheck(p_m, p_n, d_m, d_n) )
+				{
+					auto ln = n;
+					auto lm = m;
+					if(lm>ln) Swap(lm,ln);
+					if( auto res = pairs_.insert(lm,ln); res <0)
+					{
+						getFullUpdate++;
+					}
+				}
+				
+				n = next_(n);
+			}
+
+			// neighbor cells
+			CellType neighborCell;
+			for(int32 ni=0; ni<13; ni++)
+			{
+				if(ni==0)		neighborCell = currentCell + CellType( 0, 0,-1);
+				else if(ni==1)	neighborCell = currentCell + CellType(-1, 0,-1);
+				else if(ni==2)	neighborCell = currentCell + CellType(-1, 0, 0); 
+				else if(ni==3)	neighborCell = currentCell + CellType(-1, 0, 1); 
+				else if(ni==4)	neighborCell = currentCell + CellType( 0,-1,-1); 
+				else if(ni==5)	neighborCell = currentCell + CellType( 0,-1, 0); 
+				else if(ni==6)	neighborCell = currentCell + CellType( 0,-1, 1); 
+				else if(ni==7)	neighborCell = currentCell + CellType(-1,-1,-1); 
+				else if(ni==8)	neighborCell = currentCell + CellType(-1,-1, 0); 
+				else if(ni==9)	neighborCell = currentCell + CellType(-1,-1, 1);
+				else if(ni==10)	neighborCell = currentCell + CellType( 1,-1,-1);  
+				else if(ni==11)	neighborCell = currentCell + CellType( 1,-1, 0);  
+				else if(ni==12)	neighborCell = currentCell + CellType( 1,-1, 1);  
+				 
+				if( this->isInRange(neighborCell) )
+				{
+					
+					n = head_(neighborCell.x(), neighborCell.y(), neighborCell.z());
+					while( n>-1)
+					{
+						
+						auto p_n = pointPosition_[n];
+						auto d_n = sizeRatio_*diameter_[n];
+						
+						if(sphereSphereCheck(p_m, p_n, d_m, d_n))
+						{
+							auto ln = n;
+							auto lm = m;
+							if(lm>ln) Swap(lm,ln);
+							if( auto res = pairs_.insert(lm,ln); res <0)
+							{
+								getFullUpdate++;
+							}
+						}
+						n = next_[n];
+					}
+				}
+
+			}
+			m = next_[m];
+		}
+	}*/
+
+	template <typename PairsContainer, typename teamMemberType>
+	INLINE_FUNCTION_HD
+	int32 addPointsInBoxToList( 
+		const teamMemberType& teamMember,
+		IdType id,
+		const iBox<IndexType>& triBox,
+		const PairsContainer& pairs)const
+	{
+		int32 getFull = 0;
+
+		auto bExtent = boxExtent(triBox);
+		int32 numCellBox = bExtent.x()*bExtent.y()*bExtent.z();
+
+		const auto head = head_;
+		const auto next = next_;
+
+		// perform a loop over all cells in the triBox
+		Kokkos::parallel_reduce(
+			Kokkos::TeamThreadRange(teamMember,numCellBox),
+			[=](int32 linIndex, int32& valToUpdate){
+
+				CellType cell;
+				indexToCell(linIndex, triBox, cell);
+
+				int32 n = head(cell.x(),cell.y(),cell.z());
+				
+				while( n>-1)
+				{
+					
+					// id is wall id the pair is (particle id, wall id)
+					if( pairs.insert(static_cast<IdType>(n), id) < 0 )
+						valToUpdate++;
+					n = next(n);
+				}
+
+			},
+			getFull
+			);
+
+		return getFull;
+
+	}
+
+};
+
+}
+
+#endif

src/Interaction/contactSearch/methods/NBSLoop.H

@@ -0,0 +1,100 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   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 
+  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 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+int32 m = head_(i,j,k);
+CellType currentCell(i,j,k);
+int32 n = -1;
+
+while( m > -1 )
+{
+	
+	auto p_m = pointPosition_[m];
+	auto d_m = sizeRatio_* diameter_[m];
+	
+	// the same cell
+	n = next_(m);
+	
+	while(n >-1)
+	{
+		
+		auto p_n = pointPosition_[n];
+		auto d_n = sizeRatio_*diameter_[n];
+		
+		if( sphereSphereCheck(p_m, p_n, d_m, d_n) )
+		{
+			auto ln = n;
+			auto lm = m;
+			if(lm>ln) Swap(lm,ln);
+			if( auto res = pairs.insert(lm,ln); res <0)
+			{
+				getFullUpdate++;
+			}
+		}
+		
+		n = next_(n);
+	}
+
+	// neighbor cells
+	CellType neighborCell;
+	for(int32 ni=0; ni<13; ni++)
+	{
+		if(ni==0)		neighborCell = currentCell + CellType( 0, 0,-1);
+		else if(ni==1)	neighborCell = currentCell + CellType(-1, 0,-1);
+		else if(ni==2)	neighborCell = currentCell + CellType(-1, 0, 0); 
+		else if(ni==3)	neighborCell = currentCell + CellType(-1, 0, 1); 
+		else if(ni==4)	neighborCell = currentCell + CellType( 0,-1,-1); 
+		else if(ni==5)	neighborCell = currentCell + CellType( 0,-1, 0); 
+		else if(ni==6)	neighborCell = currentCell + CellType( 0,-1, 1); 
+		else if(ni==7)	neighborCell = currentCell + CellType(-1,-1,-1); 
+		else if(ni==8)	neighborCell = currentCell + CellType(-1,-1, 0); 
+		else if(ni==9)	neighborCell = currentCell + CellType(-1,-1, 1);
+		else if(ni==10)	neighborCell = currentCell + CellType( 1,-1,-1);  
+		else if(ni==11)	neighborCell = currentCell + CellType( 1,-1, 0);  
+		else if(ni==12)	neighborCell = currentCell + CellType( 1,-1, 1);  
+		 
+		if( this->isInRange(neighborCell) )
+		{
+			
+			n = head_(neighborCell.x(), neighborCell.y(), neighborCell.z());
+			while( n>-1)
+			{
+				
+				auto p_n = pointPosition_[n];
+				auto d_n = sizeRatio_*diameter_[n];
+				
+				if(sphereSphereCheck(p_m, p_n, d_m, d_n))
+				{
+					auto ln = n;
+					auto lm = m;
+					if(lm>ln) Swap(lm,ln);
+					if( auto res = pairs.insert(lm,ln); res <0)
+					{
+						getFullUpdate++;
+					}
+				}
+				n = next_[n];
+			}
+		}
+
+	}
+	m = next_[m];
+}
+

src/Interaction/contactSearch/wallMappings/cellsSimple.H

@@ -0,0 +1,338 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   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 
+  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 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __cellsSimple_H__
+#define __cellsSimple_H__
+
+#include "types.H"
+#include "KokkosTypes.H"	
+#include "cells.H"
+#include "iBox.H"
+#include "dictionary.H"
+
+
+namespace pFlow
+{
+
+template<
+	typename executionSpace,
+	typename idType,
+	typename indexType = int32
+	>
+class cellsSimple
+:
+	public cells<indexType>
+{
+public:
+
+	using IdType 		= idType;
+
+	using IndexType 	= indexType;
+
+	using Cells 		= cells<IndexType>;
+
+	using CellType 		= typename Cells::CellType;
+
+	using ExecutionSpace = executionSpace;
+
+	using memory_space 	= typename ExecutionSpace::memory_space;
+
+	using iBoxType = iBox<IndexType>;	
+
+	bool constexpr static LOOP_ELEMENT_RANGE = true;
+	
+	class TagFindCellRange2{};
+protected:
+
+	// - box extent
+	real 			cellExtent_ = 0.6;
+
+	// - update frequency 
+	int32 			updateFrequency_=1;
+
+
+	int32 			currentIter_ = 0;
+
+	// - number of triangle elements 
+	int32			numElements_ = 0;
+
+	// - number of points
+	int32 			numPoints_ = 0;
+
+	/// a broad search has been occured during last pass?
+	bool 			performedSearch_ = false;
+
+	// - ref to vectices 
+	ViewType1D<int32x3, memory_space> 	vertices_;
+
+	// - ref to points in the trisurface 
+	ViewType1D<realx3, memory_space> 	points_;
+
+	// cell range of element/triangle bounding box
+	ViewType1D<iBoxType, memory_space> 	elementBox_;
+	
+
+	using tpPWContactSearch = Kokkos::TeamPolicy<
+			ExecutionSpace,
+			Kokkos::Schedule<Kokkos::Dynamic>,
+			Kokkos::IndexType<int32>
+			>;
+
+	using rpFindCellRange2Type =
+		Kokkos::RangePolicy<TagFindCellRange2, ExecutionSpace, Kokkos::IndexType<int32>>;
+
+	
+	FUNCTION_H
+	void allocateArrays()
+	{
+		Kokkos::realloc( elementBox_, numElements_);
+	}
+
+private:
+
+	bool performSearch()
+	{
+		if(currentIter_ % updateFrequency_ == 0)
+		{
+			currentIter_++;
+			return true;
+
+		}else
+		{
+			currentIter_++;
+			return false;
+		}
+	}
+
+public:
+
+	TypeNameNV("cellsSimple");
+
+	FUNCTION_H
+	cellsSimple(
+		Cells 		ppCells,
+		real 		cellExtent,
+		int32 		numPoints,
+		int32 		numElements,
+		const ViewType1D<realx3,memory_space>& 		points, 
+		const ViewType1D<int32x3,memory_space>& 	vertices
+		)
+	:
+		Cells(ppCells),
+		cellExtent_( max(cellExtent, 0.6 ) ),
+		numElements_(numElements),
+		numPoints_(numPoints),
+		vertices_(vertices),
+		points_(points)
+	{
+
+		allocateArrays();
+	}
+
+	cellsSimple(
+		dictionary  dict,
+		Cells 		ppCells,
+		int32 		numPoints,
+		int32 		numElements,
+		const ViewType1D<realx3,memory_space>& 		points, 
+		const ViewType1D<int32x3,memory_space>& 	vertices
+		)
+	:
+		Cells(ppCells),
+		numElements_(numElements),
+		numPoints_(numPoints),
+		vertices_(vertices),
+		points_(points)
+	{
+
+		updateFrequency_ = dict.getVal<int32>(
+			"updateFrequency" );
+		updateFrequency_ = max(updateFrequency_,1);
+
+		cellExtent_ = dict.getVal<real>(
+			"cellExtent");
+
+		cellExtent_ = max(cellExtent_,0.6);
+
+		allocateArrays();
+		
+	}
+
+	constexpr bool loopElementRange()const
+	{
+		return LOOP_ELEMENT_RANGE;
+	}
+
+	// - host call
+	//   reset triangle elements if they have changed
+	FUNCTION_H
+	bool resetElements( 
+		int32 		numElements,
+		int32       numPoints,
+		ViewType1D<realx3, memory_space>& points, 
+		ViewType1D<int32x3, memory_space>& vertices	)
+	{
+
+		numElements_ = numElements;
+		numPoints_ = numPoints;
+		points_ 	= points;
+		vertices_ 	= vertices;
+		
+		allocateArrays();
+
+		return true;
+	}
+
+	INLINE_FUNCTION_HD
+	iBoxType elementBox(int32 i)const
+	{
+		return elementBox_[i];
+	}
+
+	INLINE_FUNCTION_HD
+	int32 numElements()const
+	{
+		return numElements_;
+	}
+
+	bool enterBoadSearch()const
+	{
+		return currentIter_%updateFrequency_==0;
+	}
+
+	bool performedSearch()const
+	{
+		return performedSearch_;
+	}
+		
+	template<typename PairsContainer, typename particleMapType>
+	bool broadSearch(PairsContainer& pairs, particleMapType& particleMap, bool force=false)
+	{
+		if(force) currentIter_ = 0;
+		performedSearch_= false;
+		if(!performSearch())return true;
+	
+		// map walls onto the cells
+		this->build();
+		
+		this->particleWallFindPairs(pairs, particleMap);
+		
+		performedSearch_ = true;
+		return true;
+	}
+
+	bool build()
+	{
+		Kokkos::parallel_for(
+			"cellsSimple::findcellrange2",
+			rpFindCellRange2Type(0,numElements_),
+			*this);
+		Kokkos::fence();
+		return true;
+	}
+
+	template<typename PairsContainer, typename particleMapType>
+	bool particleWallFindPairs(PairsContainer& pairs, particleMapType& particleMap)
+	{
+
+		int32 getFull = 1;
+
+		while (getFull)
+		{
+			
+			getFull = findPairsElementRange(pairs, particleMap);
+			
+			if(getFull)
+			{
+				// - resize the container
+				//   note that getFull now shows the number of failed insertions.
+				uint32 len = max(getFull, 50);
+				auto oldCap = pairs.capacity();
+				pairs.increaseCapacityBy(len);
+
+				Info<<"Contact pair container capacity increased from "<<
+						oldCap << " to " 
+						<< pairs.capacity() <<" in cellsSimple."<<endInfo;
+				
+				
+				
+				Kokkos::fence();
+			}
+		}
+		
+		return true;
+	}
+
+	template<typename PairsContainer, typename particleMapType>
+	int32 findPairsElementRange(PairsContainer& pairs, particleMapType& particleMap) 
+	{
+		int32 getFull =0;
+		
+		const auto pwPairs = pairs;
+		const auto elementBox = elementBox_;
+
+		
+		Kokkos::parallel_reduce(
+			"cellsSimple::findPairsElementRange",
+			tpPWContactSearch(numElements_, Kokkos::AUTO),
+			LAMBDA_HD(
+				const typename tpPWContactSearch::member_type & teamMember,
+				int32& valueToUpdate){
+				int32 i  = teamMember.league_rank();
+				IdType id = i;
+				const auto triBox = elementBox[i];
+				
+				valueToUpdate +=
+					particleMap.addPointsInBoxToList(
+						teamMember,
+						id,
+						triBox,
+						pwPairs
+						);
+			},
+			getFull
+			);
+
+		return getFull;
+	}
+
+	
+	INLINE_FUNCTION_HD
+	void operator()(TagFindCellRange2, int32 i) const
+	{
+		auto v = vertices_[i];
+		auto p1 = points_[v.x()];
+		auto p2 = points_[v.y()];
+		auto p3 = points_[v.z()];
+		
+		realx3 minP, maxP;
+
+		this->extendBox(p1, p2, p3, cellExtent_, minP, maxP);
+		//output<< minP << " maxP "<< maxP<<endl;
+		elementBox_[i] = iBoxType(this->pointIndex(minP), this->pointIndex(maxP));
+
+	}
+
+};
+
+}
+
+
+#endif