diff --git a/src/Particles/particles/particles.cpp b/src/Particles/particles/particles.cpp
index ae6e89c1..d5c3bb45 100644
--- a/src/Particles/particles/particles.cpp
+++ b/src/Particles/particles/particles.cpp
@@ -150,6 +150,35 @@ pFlow::particles::particles
 
 }
 
+bool pFlow::particles::beforeIteration() 
+{
+	auto domain = this->control().domain();
+
+	auto numMarked = dynPointStruct_.markDeleteOutOfBox(domain);
+	
+	if(time_.sortTime())
+	{
+		real min_dx, max_dx;
+		boundingSphereMinMax(min_dx, max_dx);
+		Timer t;
+		t.start();
+		REPORT(0)<<"Performing morton sorting on particles ...."<<endREPORT;
+		if(!pStruct().mortonSortPoints(domain, max_dx))
+		{
+			fatalErrorInFunction<<"Morton sorting was not successful"<<endl;
+			return false;
+		}
+		t.end();
+		REPORT(1)<<"Active range is "<< pStruct().activeRange()<<endREPORT;
+		REPORT(1)<<"It took "<< yellowText(t.totalTime())<<" seconds."<<endREPORT;
+	}
+
+	this->zeroForce();
+	this->zeroTorque();
+
+	return true;
+}
+
 pFlow::uniquePtr<pFlow::List<pFlow::eventObserver*>> 
 pFlow::particles::getFieldObjectList()const
 {
diff --git a/src/Particles/particles/particles.hpp b/src/Particles/particles/particles.hpp
index 0341ef5a..7b207c50 100644
--- a/src/Particles/particles/particles.hpp
+++ b/src/Particles/particles/particles.hpp
@@ -241,18 +241,7 @@ public:
 		return shapeName_;
 	}
 
-	bool beforeIteration() override
-	{
-		auto domain = this->control().domain();
-
-		auto numMarked = dynPointStruct_.markDeleteOutOfBox(domain);
-		
-
-		this->zeroForce();
-		this->zeroTorque();
-
-		return true;
-	}
+	bool beforeIteration() override;
 
 	virtual
 	bool insertParticles
diff --git a/src/phasicFlow/CMakeLists.txt b/src/phasicFlow/CMakeLists.txt
index 34103498..b6b68815 100644
--- a/src/phasicFlow/CMakeLists.txt
+++ b/src/phasicFlow/CMakeLists.txt
@@ -45,12 +45,14 @@ repository/IOobject/IOobject.cpp
 repository/IOobject/IOfileHeader.cpp
 
 structuredData/box/box.cpp
+structuredData/cells/cells.cpp
 structuredData/cylinder/cylinder.cpp
 structuredData/sphere/sphere.cpp
 structuredData/iBox/iBoxs.cpp
 structuredData/line/line.cpp
 structuredData/zAxis/zAxis.cpp
 structuredData/pointStructure/pointStructure.cpp
+structuredData/pointStructure/mortonIndexing.cpp
 structuredData/pointStructure/selectors/pStructSelector/pStructSelector.cpp
 structuredData/pointStructure/selectors/selectBox/selectBox.cpp
 structuredData/pointStructure/selectors/selectRange/selectRange.cpp 
diff --git a/src/phasicFlow/Kokkos/KokkosTypes.hpp b/src/phasicFlow/Kokkos/KokkosTypes.hpp
index b378124d..4f59872f 100644
--- a/src/phasicFlow/Kokkos/KokkosTypes.hpp
+++ b/src/phasicFlow/Kokkos/KokkosTypes.hpp
@@ -26,6 +26,8 @@ Licence:
 #include <Kokkos_DualView.hpp>
 #include <Kokkos_UnorderedMap.hpp>
 
+#include "iOstream.hpp"
+
 
 namespace pFlow
 {
@@ -51,9 +53,12 @@ using DefaultExecutionSpace 		= Kokkos::DefaultExecutionSpace;
 template<typename T1, typename T2>
 	using kPair = Kokkos::pair<T1,T2>;
 
-using range = kPair<int,int>;
+template<typename T>
+	using kRange = kPair<T,T>;
 
-using range64 = kPair<int long,int long>;
+using range = kRange<int>;
+
+using range64 = kRange<int long>;
 
 template<typename T, typename... properties>
 	using ViewTypeScalar = Kokkos::View<T,properties...>;
@@ -132,6 +137,13 @@ using deviceAtomicViewType3D =
 	T***,
 	Kokkos::MemoryTraits<std::is_same<DefaultExecutionSpace,Serial>::value?0:Kokkos::Atomic>>;
 
+template<typename T>
+iOstream& operator <<(iOstream& os, const kRange<T>& rng)
+{
+	os<<"["<<rng.first<<" "<<rng.second<<")";
+	return os;
+}
+
 
 } // pFlow
 
diff --git a/src/phasicFlow/algorithms/stdAlgorithms.hpp b/src/phasicFlow/algorithms/stdAlgorithms.hpp
index 0f8c75a8..94a46804 100644
--- a/src/phasicFlow/algorithms/stdAlgorithms.hpp
+++ b/src/phasicFlow/algorithms/stdAlgorithms.hpp
@@ -195,8 +195,8 @@ void permuteSort(const Type* first, PermuteType* pFirst, int32 numElems)
 	};
 
 	compOperator compare{first};
-	fillSequence<Type, useParallel>(pFirst, numElems, static_cast<PermuteType>(0));
-	sort<Type, compOperator, useParallel>(pFirst, numElems, compare);
+	fillSequence<PermuteType, useParallel>(pFirst, numElems, static_cast<PermuteType>(0));
+	sort<PermuteType, compOperator, useParallel>(pFirst, numElems, compare);
 }
 
 
diff --git a/src/phasicFlow/containers/Vector/Vector.cpp b/src/phasicFlow/containers/Vector/Vector.cpp
index f979f5f3..b9d76479 100644
--- a/src/phasicFlow/containers/Vector/Vector.cpp
+++ b/src/phasicFlow/containers/Vector/Vector.cpp
@@ -197,6 +197,26 @@ bool pFlow::Vector<T, Allocator>::deleteElement
 		return false;
 }
 
+template<typename T, typename Allocator>
+void pFlow::Vector<T, Allocator>::sortItems(
+	const int32IndexContainer& indices)
+{
+	if(indices.size() == 0)
+	{
+		this->resize(0);
+		return;
+	}
+	size_t newSize = indices.size();
+	auto hIndices = indices.hostView();
+	VectorType sortedVec(name(), capacity(), newSize, RESERVE());
+	
+	ForAll(i, hIndices)
+	{
+		sortedVec[i] = vectorType::operator[](i); 
+	}
+	*this = std::move(sortedVec);
+}
+
 template<typename T, typename Allocator>
 bool pFlow::Vector<T, Allocator>::insertSetElement(
 	const  int32IndexContainer& indices,
diff --git a/src/phasicFlow/containers/Vector/Vector.hpp b/src/phasicFlow/containers/Vector/Vector.hpp
index c0b1f17a..2cfccaf0 100644
--- a/src/phasicFlow/containers/Vector/Vector.hpp
+++ b/src/phasicFlow/containers/Vector/Vector.hpp
@@ -323,6 +323,9 @@ public:
 	//   return false if out of range
 	bool deleteElement(label index);
 
+	/// Sort elements based on the indices
+	void sortItems(const int32IndexContainer& indices);
+
 	// - set or insert new elements into the vector 
 	//   return false if it fails
 	bool insertSetElement(const  int32IndexContainer& indices, const T& val);
diff --git a/src/phasicFlow/containers/VectorHD/VectorDual.hpp b/src/phasicFlow/containers/VectorHD/VectorDual.hpp
index 0839601d..f5dd5611 100644
--- a/src/phasicFlow/containers/VectorHD/VectorDual.hpp
+++ b/src/phasicFlow/containers/VectorHD/VectorDual.hpp
@@ -545,6 +545,35 @@ public:
 			return true;
 		}
 
+		INLINE_FUNCTION_H
+		void sortItems(const int32IndexContainer& indices)
+		{
+			if(indices.size() == 0)
+			{
+				setSize(0);
+				return;
+			}
+			size_t newSize = indices.size();
+
+			deviceViewType 	sortedView("sortedView", newSize);
+			auto dVec = deviceVectorAll();
+
+			auto d_indices = indices.deviceView();
+			
+			Kokkos::parallel_for(
+				"sortItems",
+				newSize,
+				LAMBDA_HD(int32 i){
+					sortedView[i] = dVec[d_indices[i]];
+				}
+				);
+			Kokkos::fence();
+			setSize(newSize);
+			copy(deviceVector(), sortedView);
+			modifyOnDevice();
+			syncViews();
+		}
+
 		INLINE_FUNCTION_H
 		bool insertSetElement(const int32IndexContainer& indices, const T& val)
 		{
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
index cbe173aa..7dfd4105 100644
--- a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
+++ b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
@@ -558,6 +558,46 @@ public:
 			return false;
 		}
 
+		INLINE_FUNCTION_H
+		void sortItems(const int32IndexContainer& indices)
+		{
+			if(indices.size() == 0)
+			{
+				setSize(0);
+				return;
+			}
+			
+			size_t 		newSize = indices.size();
+			viewType 	sortedView("sortedView", newSize);
+
+			if constexpr (isHostAccessible_)
+			{
+				auto h_indices = indices.hostView();
+				Kokkos::parallel_for(
+				"sortItems",
+				newSize, 
+				LAMBDA_HD(int32 i){	
+					sortedView[i] = view_[h_indices[i]];
+				});
+			}else
+			{
+				auto d_indices = indices.deviceView();
+				Kokkos::parallel_for(
+				"sortItems",
+				newSize, 
+				LAMBDA_HD(int32 i){	
+					sortedView[i] = view_[d_indices[i]];
+				});
+			}
+
+			Kokkos::fence();
+			setSize(newSize);
+
+			copy(deviceVector(), sortedView);
+			return;
+
+		}
+
 		INLINE_FUNCTION_H
 		bool insertSetElement(const int32IndexContainer& indices, const Vector<T>& vals)
 		{
diff --git a/src/phasicFlow/containers/indexContainer/indexContainer.cpp b/src/phasicFlow/containers/indexContainer/indexContainer.cpp
index 15d717b2..744f45a9 100644
--- a/src/phasicFlow/containers/indexContainer/indexContainer.cpp
+++ b/src/phasicFlow/containers/indexContainer/indexContainer.cpp
@@ -18,4 +18,7 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-#include "indexContainer.hpp"
\ No newline at end of file
+#include "indexContainer.hpp"
+
+
+template class pFlow::indexContainer<pFlow::int32>;
\ No newline at end of file
diff --git a/src/phasicFlow/containers/indexContainer/indexContainer.hpp b/src/phasicFlow/containers/indexContainer/indexContainer.hpp
index 14e88726..7b517a5a 100644
--- a/src/phasicFlow/containers/indexContainer/indexContainer.hpp
+++ b/src/phasicFlow/containers/indexContainer/indexContainer.hpp
@@ -43,6 +43,10 @@ public:
   	// - viewType of data on host
   	using HostViewType 		= typename DualViewType::t_host;
 
+  	using HostType 			= typename HostViewType::device_type;
+
+  	using DeviceType 		= typename DeviceViewType::device_type;
+
   	template<typename ViewType>
   	class IndexAccessor
   	{
@@ -99,7 +103,11 @@ public:
 
 	indexContainer(const indexContainer&) = default;
 
-	indexContainer& operator = (const indexContainer&) = default;	
+	indexContainer& operator = (const indexContainer&) = default;
+
+	indexContainer(indexContainer&&) = default;
+
+	indexContainer& operator = (indexContainer&&) = default;	
 
 	~indexContainer() = default;
 
@@ -150,6 +158,16 @@ public:
 		return views_.d_view;
 	}
 
+	HostViewType& hostView()
+	{
+		return views_.h_view;
+	}
+
+	DeviceViewType& deviceView()
+	{
+		return views_.d_view;
+	}
+
 	auto indicesHost()const
 	{
 		return IndexAccessor<HostViewType>(views_.h_view);
@@ -157,7 +175,45 @@ public:
 
 	auto indicesDevice()const
 	{
-		return IndexAccessor<DeviceViewType>(views_.d_veiw);
+		return IndexAccessor<DeviceViewType>(views_.d_view);
+	}
+
+	void modifyOnHost()
+	{
+		views_.modify_host();
+	}
+
+	void modifyOnDevice()
+	{
+		views_.modify_device();
+	}
+
+	void syncViews()
+	{
+		bool findMinMax = false;
+		if(views_.template need_sync<HostType>())
+		{
+			Kokkos::deep_copy(views_.d_view, views_.h_view);
+			findMinMax = true;
+		}
+		else if(views_.template need_sync<DeviceType>())
+		{
+			Kokkos::deep_copy(views_.h_view, views_.d_view);
+			findMinMax = true;
+		}
+
+		if(findMinMax)
+		{
+			min_ 	= pFlow::min(views_.d_view, 0, size_);
+			max_ 	= pFlow::max(views_.d_view, 0, size_);
+		}
+	}
+
+	size_t setSize(size_t ns)
+	{
+		auto tmp = size_;
+		size_ = ns;
+		return tmp;
 	}
 };
 
diff --git a/src/phasicFlow/containers/pointField/pointField.cpp b/src/phasicFlow/containers/pointField/pointField.cpp
index 41c5230b..bc0a62e0 100644
--- a/src/phasicFlow/containers/pointField/pointField.cpp
+++ b/src/phasicFlow/containers/pointField/pointField.cpp
@@ -123,6 +123,13 @@ bool pFlow::pointField<VectorField, T, MemorySpace>::update(const eventMessage&
 		//Vector<T> vals( newElems.size(), defaultValue_);
 		return this->insertSetElement(newElems, defaultValue_);
 	}
+
+	if(msg.isRearranged())
+	{
+		auto sortedIndex = pStruct().mortonSortedIndex();
+		this->sortItems(sortedIndex);
+		return true;
+	}
 	
 	return true;
 }
diff --git a/src/phasicFlow/repository/Time/timeControl.cpp b/src/phasicFlow/repository/Time/timeControl.cpp
index 6fbad2ad..cb6a162a 100644
--- a/src/phasicFlow/repository/Time/timeControl.cpp
+++ b/src/phasicFlow/repository/Time/timeControl.cpp
@@ -61,8 +61,16 @@ pFlow::timeControl::timeControl
 	(
 		startTime_,
 		dict.getValOrSet("timersReportInterval", 0.04)
+	),
+	performSorting_
+	(
+		dict.getValOrSet("performSorting", Logical("No"))
+	),
+	sortingInterval_
+	(
+		startTime_,
+		dict.getValOrSet("sortingInterval", static_cast<real>(1.0))
 	)
-
 {
 	checkForOutputToFile();
 }
@@ -95,6 +103,15 @@ pFlow::timeControl::timeControl(
 	(
 		startTime_,
 		dict.getValOrSet("timersReportInterval", 0.04)
+	),
+	performSorting_
+	(
+		dict.getValOrSet("performSorting", Logical("No"))
+	),
+	sortingInterval_
+	(
+		startTime_,
+		dict.getValOrSet("sortingInterval", static_cast<real>(1.0))
 	)	
 {
 	checkForOutputToFile();
@@ -160,6 +177,11 @@ bool pFlow::timeControl::timersReportTime()const
 	return timersReportInterval_.isMember(currentTime_, dt_);
 }
 
+bool pFlow::timeControl::sortTime()const
+{
+	return performSorting_()&&sortingInterval_.isMember(currentTime_,dt_);
+}
+
 void pFlow::timeControl::setSaveTimeFolder(
 	bool saveToFile, 
 	const word& timeName)
diff --git a/src/phasicFlow/repository/Time/timeControl.hpp b/src/phasicFlow/repository/Time/timeControl.hpp
index b3ac5c1b..9503d8ea 100644
--- a/src/phasicFlow/repository/Time/timeControl.hpp
+++ b/src/phasicFlow/repository/Time/timeControl.hpp
@@ -74,9 +74,13 @@ protected:
 
 	real 	writeTime_ = 0; // for managedExternamly
 
-	realStridedRange  timersReportInterval_;
+	realStridedRange  	timersReportInterval_;
 
-	int32StridedRagne screenReportInterval_ ={0,100};
+	Logical				performSorting_;
+
+	realStridedRange 	sortingInterval_;
+
+	int32StridedRagne 	screenReportInterval_ ={0,100};
 
 	bool outputToFile_ = false;
 
@@ -164,6 +168,8 @@ public:
 	}
 
 	bool timersReportTime()const;
+
+	bool sortTime()const;
 	
 	bool setOutputToFile(real writeTime, const word& timeName)
 	{
diff --git a/src/phasicFlow/structuredData/cells/cells.cpp b/src/phasicFlow/structuredData/cells/cells.cpp
new file mode 100644
index 00000000..bf5b2d2a
--- /dev/null
+++ b/src/phasicFlow/structuredData/cells/cells.cpp
@@ -0,0 +1,21 @@
+/*------------------------------- 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 "cells.hpp"
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/cells/cells.hpp b/src/phasicFlow/structuredData/cells/cells.hpp
new file mode 100644
index 00000000..c77b1a69
--- /dev/null
+++ b/src/phasicFlow/structuredData/cells/cells.hpp
@@ -0,0 +1,259 @@
+/*------------------------------- 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_hpp__
+#define __cells_hpp__
+
+
+#include "types.hpp"
+#include "box.hpp"
+
+namespace pFlow
+{
+
+template<typename indexType>
+class cells
+{
+public:
+	
+	using CellType = triple<indexType>;
+
+protected:
+
+	// - domain 
+	box 		domain_{realx3(0.0), realx3(1.0)};
+
+	// - cell size 
+	realx3 		cellSize_{1,1,1};
+
+	CellType 	numCells_{1,1,1};
+
+	
+	// - 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_H
+	cells(const box& domain, int32 nx, int32 ny, int32 nz)
+	:
+		domain_(domain),
+		cellSize_(
+			(domain_.maxPoint() - domain_.minPoint())/realx3(nx, ny,  nz)
+			),
+		numCells_(nx, ny, nz)
+	{}
+
+	INLINE_FUNCTION_HD
+	cells(const cells&) = default;
+
+	INLINE_FUNCTION_HD
+	cells& operator = (const cells&) = default;
+
+	INLINE_FUNCTION_HD
+	cells(cells &&) = default;
+
+	INLINE_FUNCTION_HD
+	cells& operator=(cells&&) = default;
+
+	cells getCells()const
+	{
+		return *this;
+	}
+
+	INLINE_FUNCTION_H
+	void setCellSize(real cellSize)
+	{
+		cellSize_ = cellSize;
+		calculate();
+	}
+
+	INLINE_FUNCTION_H
+	void setCellSize(realx3 cellSize)
+	{
+		cellSize_ = cellSize;
+		calculate();
+	}
+
+	INLINE_FUNCTION_HD
+	realx3 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());
+	}
+
+	const auto& domain()const
+	{
+		return domain_;
+	}
+	
+	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
diff --git a/src/phasicFlow/structuredData/pointStructure/mortonIndexing.cpp b/src/phasicFlow/structuredData/pointStructure/mortonIndexing.cpp
new file mode 100644
index 00000000..c3b55e55
--- /dev/null
+++ b/src/phasicFlow/structuredData/pointStructure/mortonIndexing.cpp
@@ -0,0 +1,75 @@
+/*------------------------------- 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 "mortonIndexing.hpp"
+#include "cells.hpp"
+
+
+bool pFlow::getSortedIndex(
+	box boundingBox,
+	real dx,
+	range activeRange, 
+	ViewType1D<realx3> pos, 
+	ViewType1D<int8> flag,
+	int32IndexContainer& sortedIndex)
+{
+
+	// obtain the morton code of the particles
+	cells<size_t> allCells( boundingBox, dx);
+	int32IndexContainer index(activeRange.first, activeRange.second);
+
+	ViewType1D<uint64_t> mortonCode("mortonCode", activeRange.second);
+
+	using rpMorton = 
+		Kokkos::RangePolicy<Kokkos::IndexType<int32>>;
+	int32 numActive = 0;
+	Kokkos::parallel_reduce
+	(
+		"mortonIndexing::getIndex::morton",
+		rpMorton(activeRange.first, activeRange.second),
+		LAMBDA_HD(int32 i, int32& sumToUpdate){
+			if( flag[i] == 1 ) // active point 
+			{
+				auto cellInd = allCells.pointIndex(pos[i]);
+				mortonCode[i] = xyzToMortonCode64(cellInd.x(), cellInd.y(), cellInd.z());
+				sumToUpdate++;
+			}else
+			{
+				mortonCode[i] = xyzToMortonCode64(-1,-1,-1);
+			}
+		},
+		numActive
+	);
+	
+
+	permuteSort(
+		mortonCode, 
+		activeRange.first, 
+		activeRange.second,
+		index.deviceView(),
+		0 );
+	index.modifyOnDevice();
+	index.setSize(numActive);
+	index.syncViews();
+
+	sortedIndex = index;
+
+	return true;
+}
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/pointStructure/mortonIndexing.hpp b/src/phasicFlow/structuredData/pointStructure/mortonIndexing.hpp
new file mode 100644
index 00000000..9ccdc83e
--- /dev/null
+++ b/src/phasicFlow/structuredData/pointStructure/mortonIndexing.hpp
@@ -0,0 +1,88 @@
+/*------------------------------- 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 __mortonIndexing_hpp__
+#define __mortonIndexing_hpp__
+
+#include "types.hpp"
+#include "box.hpp"
+#include "indexContainer.hpp"
+
+namespace pFlow
+{
+
+bool getSortedIndex(
+	box boundingBox,
+	real dx, 
+	range activeRange, 
+	ViewType1D<realx3> pos, 
+	ViewType1D<int8> flag,
+	int32IndexContainer& sortedIndex);
+
+
+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);
+}
+
+struct indexMorton
+{
+	size_t morton;
+	size_t index;
+};
+
+
+}
+
+#endif //__mortonIndexing_hpp__
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp b/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp
index cacb6280..a36c8e9d 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp
@@ -24,7 +24,8 @@ Licence:
 #include "setFieldList.hpp"
 #include "error.hpp"
 #include "iOstream.hpp"
-#include "Time.hpp"
+//#include "Time.hpp"
+#include "mortonIndexing.hpp"
 
 FUNCTION_H
 bool pFlow::pointStructure::evaluatePointStructure()
@@ -231,6 +232,55 @@ bool pFlow::pointStructure::allActive()const
 	return numActivePoints_ == numPoints_;
 }
 
+FUNCTION_H
+bool pFlow::pointStructure::mortonSortPoints(const box& domain, real dx)
+{
+	if( !getSortedIndex(
+		domain, 
+		dx, 
+		activeRange_,
+		pointPosition_.deviceVectorAll(),
+		pointFlag_.deviceVectorAll(),
+		mortonSortedIndex_) )
+	{
+		fatalErrorInFunction<<"failed to perform morton sorting!"<<endl;
+		return false;
+	}
+
+	pointPosition_.sortItems(mortonSortedIndex_);
+	pointFlag_.sortItems(mortonSortedIndex_);
+
+	auto oldSize 			= size();
+	auto oldCapacity 	= capacity();
+	auto oldRange 		= activeRange();
+
+	// update size, range, capacity
+	setNumMaxPoints();
+	activeRange_ = {0, static_cast<int>(mortonSortedIndex_.size())};
+	numActivePoints_ = mortonSortedIndex_.size();
+
+	eventMessage msg(eventMessage::REARRANGE);
+
+	if(oldSize != size() )
+	{
+		msg.add(eventMessage::SIZE_CHANGED);
+	}
+
+	if(oldCapacity != capacity())
+	{
+		msg.add(eventMessage::CAP_CHANGED);
+	}
+
+	if( oldRange != activeRange_)
+	{
+		msg.add(eventMessage::RANGE_CHANGED);
+	}
+
+	// notify all the registered objects except the exclusionList 
+	if( !this->notify(msg) ) return false;
+
+	return true;
+}
 
 FUNCTION_H
 size_t pFlow::pointStructure::markDeleteOutOfBox(const box& domain)
@@ -306,7 +356,6 @@ pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::pointStructure::insertPoints
 )
 {
 
-
 	auto numNew = pos.size();
 	if( numNew==0)
 	{
@@ -359,8 +408,8 @@ pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::pointStructure::insertPoints
 	}
 
 	// changes the active rage based on the new inserted points 
-	activeRange_ = { min(activeRange_.first, minInd ),
-					max(activeRange_.second, maxInd+1)};
+	activeRange_ = { static_cast<int>(min(activeRange_.first, minInd )),
+					 static_cast<int>(max(activeRange_.second, maxInd+1))};
 
 	numActivePoints_ += numNew;
 
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp b/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp
index 4f53f2ce..be70072d 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp
@@ -163,8 +163,11 @@ protected:
 		// index range of active points (half-open range)
 		range 					activeRange_;	
 
-		// - index vector for points to be inserted 
-		int32IndexContainer   		tobeInsertedIndex_;
+		/// Index vector for points to be inserted 
+		int32IndexContainer   	tobeInsertedIndex_;
+
+		/// Sorted index of particles based on morton code 
+		int32IndexContainer 	mortonSortedIndex_;
 	
 
 	//// - protected methods
@@ -298,6 +301,10 @@ public:
 		FUNCTION_H
 		virtual bool updateForDelete();
 
+
+		FUNCTION_H
+		virtual bool mortonSortPoints(const box& domain, real dx);
+
 		///////////////////////////////////////////////////////////////////////////////////////////////////
 		
 		// - const access to points to be newly inserted
@@ -320,6 +327,13 @@ public:
 		}
 
 
+		FUNCTION_H
+		auto mortonSortedIndex()const
+		{
+			return mortonSortedIndex_;
+		}
+
+
 		// - update data structure by inserting/setting new points 
 		//   Notifies all the fields in the registered list of data structure
 		//   and exclude the fields that re in the exclusionList