diff --git a/cmake/globals.cmake b/cmake/globals.cmake
index b2c084a3..34c075b4 100644
--- a/cmake/globals.cmake
+++ b/cmake/globals.cmake
@@ -4,9 +4,9 @@ list(APPEND validFiles *.C *.cpp *.cxx *.c *.cu *.H *.hpp *.hxx *.h *.cuh)
 
 macro(Kokkos_cmake_settings)
 
-mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
+#mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
 mark_as_advanced(FORCE var Kokkos_CXX_STANDARD)
-mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_CONSTEXPR)
+#mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_CONSTEXPR)
 mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMP)
 mark_as_advanced(FORCE var Kokkos_ENABLE_SERIAL)
 mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA)
diff --git a/src/phasicFlow/CMakeLists.txt b/src/phasicFlow/CMakeLists.txt
index 5114f190..d9d701eb 100644
--- a/src/phasicFlow/CMakeLists.txt
+++ b/src/phasicFlow/CMakeLists.txt
@@ -21,7 +21,8 @@ streams/TStream/iTstream.cpp
 streams/TStream/oTstream.cpp
 streams/Fstream/iFstream.cpp
 streams/Fstream/oFstream.cpp
-streams/Fstream/fileStream.cpp 
+streams/Fstream/fileStream.cpp
+streams/dataIO/dataIO.cpp 
 streams/streams.cpp 
 
 fileSystem/fileSystem.cpp
@@ -33,6 +34,12 @@ dictionary/entry/iEntry.cpp
 dictionary/entry/dataEntry.cpp 
 dictionary/twoPartEntry/twoPartEntry.cpp
 
+containers/Vector/Vectors.cpp
+
+repository/IOobject/objectFile.cpp
+repository/IOobject/IOfileHeader.cpp
+repository/IOobject/IOobject.cpp
+
 )
 
 set(link_libs)
diff --git a/src/phasicFlow/Kokkos/KokkosUtilities.hpp b/src/phasicFlow/Kokkos/KokkosUtilities.hpp
index 1b42bba7..451920ab 100644
--- a/src/phasicFlow/Kokkos/KokkosUtilities.hpp
+++ b/src/phasicFlow/Kokkos/KokkosUtilities.hpp
@@ -25,6 +25,7 @@ Licence:
 #include "KokkosTypes.hpp"
 #include "pFlowMacros.hpp"
 #include "types.hpp"
+#include "span.hpp"
 #include "iOstream.hpp"
 
 namespace pFlow
@@ -41,7 +42,7 @@ template<typename ExecutionSpace>
 INLINE_FUNCTION_H
 bool constexpr isDeviceAccessible()
 {
-	return Kokkos::SpaceAccessibility<ExecutionSpace,DefaultExecutionSpace>::accessible;
+	return Kokkos::SpaceAccessibility<ExecutionSpace,DefaultExecutionSpace::memory_space>::accessible;
 }
 
 /// Is MemoerySpace accessible from ExecutionSpace
@@ -56,7 +57,7 @@ template <
 	typename Type,
 	typename... Properties>
 INLINE_FUNCTION_H
-void reallocInit( ViewType1D<Type,Properties...>& view, int32 len)
+void reallocInit( ViewType1D<Type,Properties...>& view, uint32 len)
 {
 	Kokkos::realloc(Kokkos::WithoutInitializing, view, len);	
 }
@@ -65,7 +66,7 @@ template <
 	typename Type,
 	typename... Properties>
 INLINE_FUNCTION_H
-void reallocNoInit(ViewType1D<Type,Properties...>& view, int32 len)
+void reallocNoInit(ViewType1D<Type,Properties...>& view, uint32 len)
 {
 	Kokkos::realloc(Kokkos::WithoutInitializing, view, len);
 }
@@ -74,7 +75,7 @@ template <
 	typename Type,
 	typename... Properties>
 INLINE_FUNCTION_H
-void reallocFill( ViewType1D<Type,Properties...>& view, int32 len, Type val)
+void reallocFill( ViewType1D<Type,Properties...>& view, uint32 len, Type val)
 {
 	reallocNoInit(view, len);
 	Kokkos::deep_copy(view, val);
@@ -84,7 +85,7 @@ template <
 	typename Type,
 	typename... Properties>
 INLINE_FUNCTION_H
-void reallocInit( ViewType2D<Type,Properties...>& view, int32 len1, int32 len2)
+void reallocInit( ViewType2D<Type,Properties...>& view, uint32 len1, uint32 len2)
 {
 	Kokkos::realloc(view, len1, len2);
 }
@@ -93,7 +94,7 @@ template <
 	typename Type,
 	typename... Properties>
 INLINE_FUNCTION_H
-void reallocNoInit(ViewType2D<Type,Properties...>& view, int32 len1, int32 len2)
+void reallocNoInit(ViewType2D<Type,Properties...>& view, uint32 len1, uint32 len2)
 {	
   	Kokkos::realloc(Kokkos::WithoutInitializing, view, len1, len2);
 }
@@ -102,7 +103,7 @@ template <
 	typename Type,
 	typename... Properties>
 INLINE_FUNCTION_H
-void reallocFill( ViewType2D<Type,Properties...>& view, int32 len1, int32 len2, Type val)
+void reallocFill( ViewType2D<Type,Properties...>& view, uint32 len1, uint32 len2, Type val)
 {
 	reallocNoInit(view, len1, len2);
 	Kokkos::deep_copy(view, val);
@@ -112,7 +113,7 @@ template <
 	typename Type,
 	typename... Properties>
 INLINE_FUNCTION_H
-void reallocInit( ViewType3D<Type,Properties...>& view, int32 len1, int32 len2, int32 len3)
+void reallocInit( ViewType3D<Type,Properties...>& view, uint32 len1, uint32 len2, uint32 len3)
 {
 	Kokkos::realloc(view, len1, len2, len3);
 }
@@ -121,7 +122,7 @@ template <
 	typename Type,
 	typename... Properties>
 INLINE_FUNCTION_H
-void reallocNoInit(ViewType3D<Type,Properties...>& view, int32 len1, int32 len2, int32 len3)
+void reallocNoInit(ViewType3D<Type,Properties...>& view, uint32 len1, uint32 len2, uint32 len3)
 {
 	
   	Kokkos::realloc(Kokkos::WithoutInitializing, view, len1, len2, len3);
@@ -131,7 +132,7 @@ template <
 	typename Type,
 	typename... Properties>
 INLINE_FUNCTION_H
-void reallocFill( ViewType3D<Type,Properties...>& view, int32 len1, int32 len2, int32 len3, Type val)
+void reallocFill( ViewType3D<Type,Properties...>& view, uint32 len1, uint32 len2, uint32 len3, Type val)
 {
 	reallocNoInit(view, len1, len2, len3);
 	Kokkos::deep_copy(view, val);
@@ -141,7 +142,7 @@ template <
 	typename Type,
 	typename... Properties>
 INLINE_FUNCTION_H
-void resizeInit(ViewType1D<Type,Properties...>& view, int32 newLen)
+void resizeInit(ViewType1D<Type,Properties...>& view, uint32 newLen)
 {
 	Kokkos::resize(view, newLen);
 }
@@ -150,7 +151,7 @@ template <
 	typename Type,
 	typename... Properties>
 INLINE_FUNCTION_H
-void resizeNoInit(ViewType1D<Type,Properties...>& view, int32 newLen)
+void resizeNoInit(ViewType1D<Type,Properties...>& view, uint32 newLen)
 {
 	Kokkos::resize(Kokkos::WithoutInitializing, view, newLen);
 }
@@ -177,6 +178,22 @@ iOstream& operator <<(iOstream& os, const Pair<T1,T2>& p)
 	return os;
 }
 
+template<typename T, typename... properties>
+INLINE_FUNCTION_H
+iOstream& operator <<(iOstream& os, const ViewType1D<T, properties...> & v)
+{
+	
+	using ExSpace = typename ViewType1D<T, properties...>::execution_space;
+
+	static_assert(isHostAccessible<ExSpace>(), "View memory is not accessible from Host");
+
+	span<T> spn(v.data(), v.size());
+	os<<spn;
+	
+	return os;
+}
+
+
 } // pFlow
 
 #endif  //__KokkosUtilities_hpp__
diff --git a/src/phasicFlow/Kokkos/Range.hpp b/src/phasicFlow/Kokkos/Range.hpp
index 35b1babe..72ecd20b 100644
--- a/src/phasicFlow/Kokkos/Range.hpp
+++ b/src/phasicFlow/Kokkos/Range.hpp
@@ -129,6 +129,12 @@ public Kokkos::pair<T,T>
 			return end()-start();
 		}
 
+		INLINE_FUNCTION_HD
+		auto getPair()const
+		{
+			return Pair(this->first, this->second);
+		}
+
 };
 
 template<typename T>
@@ -143,6 +149,10 @@ using range32 	= Range<int32>;
 
 using range64 	= Range<int64>;
 
+using rangeU32 	= Range<uint32>;
+
+using rangeU64 	= Range<uint64>;
+
 
 } // pFlow
 
diff --git a/src/phasicFlow/Kokkos/ViewAlgorithms.hpp b/src/phasicFlow/Kokkos/ViewAlgorithms.hpp
index 112569c6..7e8134a9 100644
--- a/src/phasicFlow/Kokkos/ViewAlgorithms.hpp
+++ b/src/phasicFlow/Kokkos/ViewAlgorithms.hpp
@@ -35,16 +35,16 @@ namespace pFlow
 { 	
 template<typename T, typename... properties>
 INLINE_FUNCTION_H
-int32 count(
+uint32 count(
 	const ViewType1D<T, properties...>& view,
-	int32 start,
-	int32 end,
+	uint32 start,
+	uint32 end,
 	const T& val)
 {
 	
 	using ExecutionSpace = typename ViewType1D<T, properties...>::execution_space;
 	
-	int32 numElems = end-start;
+	uint32 numElems = end-start;
 
 	return  pFlow::algorithms::KOKKOS::count<T, ExecutionSpace>
 	(
@@ -59,11 +59,11 @@ INLINE_FUNCTION_H
 void fill
 (
 	ViewType1D<T, properties...>& view,
-	range32 span,
+	rangeU32 span,
 	T val
 )
 {
-	auto subV = Kokkos::subview(view, span);
+	auto subV = Kokkos::subview(view, span.getPair() );
 	Kokkos::deep_copy(subV, val);
 }
 
@@ -71,12 +71,12 @@ template<typename T, typename... properties>
 void fill
 (
 	ViewType1D<T, properties...>& view,
-	int32 start,
-	int32 end,
+	uint32 start,
+	uint32 end,
 	T val
 )
 {
-	fill(view, range32(start,end),val);
+	fill(view, rangeU32(start, end),val);
 }
 
 template<
@@ -84,14 +84,14 @@ template<
 	typename... properties>
 void fillSequence(
 	ViewType1D<Type, properties...>& view,
-	int32 start,
-	int32 end,
+	uint32 start,
+	uint32 end,
 	const Type startVal
 	)
 {
 
 	using ExecutionSpace = typename ViewType1D<Type, properties...>::execution_space;
-	int32 numElems = end-start;
+	uint32 numElems = end-start;
 
 	pFlow::algorithms::KOKKOS::fillSequence<Type, ExecutionSpace>(
 			view.data()+start,
@@ -110,24 +110,27 @@ template<
 bool fillSelected
 (
 	ViewType1D<Type, properties...> view,
-	const ViewType1D<indexType, indexProperties...> indices,
-	const int32 numElems,
-	const Type val
+	ViewType1D<indexType, indexProperties...> indices,
+	uint32 numElems,
+	Type val
 )
 {
 	static_assert(
       areAccessible<
 			typename ViewType1D<Type, properties...>::execution_space,
 			typename ViewType1D<indexType, indexProperties...>::memory_space>(),
-      "In fillSelected arguments view and indices must have similar spaces");
+      "In fillSelected, arguments view and indices must have similar spaces");
 
-	using ExecutionSpace = typename ViewType1D<Type, properties...>::execution_space;
-	
-	pFlow::algorithms::KOKKOS::fillSelected<Type, indexType, ExecutionSpace>(
-			view.data(),
-			indices.data(),
-			numElems,
-			val);
+	using ExSpace = typename ViewType1D<Type, properties...>::execution_space;
+	using policy = Kokkos::RangePolicy<ExSpace,Kokkos::IndexType<uint32> >;
+
+	Kokkos::parallel_for(
+		"ViewAlgorithms::fillSelected",
+		policy(0,numElems),
+		LAMBDA_HD(uint32 i){
+			//view[indices[i]]= val;
+		});
+	Kokkos::fence();
 
 	return true;
 }
@@ -141,7 +144,7 @@ bool fillSelected(
 	ViewType1D<Type, properties...> view,
 	const ViewType1D<indexType, indexProperties...> indices,
 	const ViewType1D<Type, properties...> vals,
-	const int32 numElems )
+	const uint32 numElems )
 {
 
 	static_assert(
@@ -167,14 +170,15 @@ template<typename T, typename... properties>
 INLINE_FUNCTION_H
 T min( 
 	const ViewType1D<T, properties...>& view,
-	int32 start,
-	int32 end)
+	uint32 start,
+	uint32 end)
 {
 
 	using ExecutionSpace = typename ViewType1D<T, properties...>::execution_space;
 	
-	int32 numElems = end-start;
- 
+	uint32 numElems = end-start;
+ 	
+ 	return 
 	pFlow::algorithms::KOKKOS::min<T, ExecutionSpace>(
 		view.data()+start,
 		numElems);
@@ -184,13 +188,13 @@ template<typename T, typename... properties>
 INLINE_FUNCTION_H
 T max( 
 	const ViewType1D<T, properties...>& view,
-	int32 start,
-	int32 end)
+	uint32 start,
+	uint32 end)
 {
 
 	using ExecutionSpace = typename ViewType1D<T, properties...>::execution_space;
 	
-	int32 numElems = end-start;
+	uint32 numElems = end-start;
 
 	return  
 	pFlow::algorithms::KOKKOS::max<T, ExecutionSpace>(
@@ -220,10 +224,10 @@ template <
 INLINE_FUNCTION_H
 void copy(
 	const ViewType1D<dType, dProperties...>& dst,
-	int32 dStart,
+	uint32 dStart,
 	const ViewType1D<sType, sProperties...>& src,
-	int32 sStart,
-	int32 sEnd
+	uint32 sStart,
+	uint32 sEnd
 	)
 {
 
@@ -244,7 +248,7 @@ INLINE_FUNCTION_H
 void getNth(
 	dType& dst,
 	const ViewType1D<sType, sProperties...>& src,
-	const int32 n
+	const uint32 n
 	)
 {
 	range32 span(n,n+1);
@@ -259,12 +263,12 @@ template<typename T, typename... properties>
 INLINE_FUNCTION_H
 void sort( 
 	ViewType1D<T, properties...>& view,
-	int32 start,
-	int32 end)
+	uint32 start,
+	uint32 end)
 {
 	using ExecutionSpace = typename ViewType1D<T, properties...>::execution_space;
 	
-	int32 numElems = end-start;
+	uint32 numElems = end-start;
 
 	if constexpr( isHostAccessible<ExecutionSpace>())
 	{		
@@ -291,14 +295,14 @@ template<typename T, typename... properties, typename CompareFunc>
 INLINE_FUNCTION_H
 void sort( 
 	ViewType1D<T, properties...>& view,
-	int32 start,
-	int32 end,
+	uint32 start,
+	uint32 end,
 	CompareFunc compare)
 {
 
 	using ExecutionSpace = typename ViewType1D<T, properties...>::execution_space;
 	
-	int32 numElems = end-start;
+	uint32 numElems = end-start;
 
 	if constexpr( isHostAccessible<ExecutionSpace>())
 	{
@@ -330,10 +334,10 @@ template<
 	typename... permProperties>
 void permuteSort(
 	const ViewType1D<Type, properties...>& view,
-	int32 start,
-	int32 end,
+	uint32 start,
+	uint32 end,
 	ViewType1D<permType, permProperties...>& permuteView,
-	int32 permStart )
+	uint32 permStart )
 {
 	static_assert(
 		areAccessible<
@@ -343,7 +347,7 @@ void permuteSort(
 
 	using ExecutionSpace = typename ViewType1D<Type, properties...>::execution_space;
 	
-	int32 numElems = end-start;
+	uint32 numElems = end-start;
 		
 	pFlow::algorithms::STD::permuteSort<Type,permType,true>(
 		view.data()+start,
@@ -367,8 +371,9 @@ void permuteSort(
 
 template<typename T>
 INLINE_FUNCTION_HD
-int binarySearch_(const T* array, int length, const T& val)
+int32 binarySearch_(const T* array, int32 length, const T& val)
 {
+    if(length <= 0) return -1;
     
     int low = 0;            
     int high = length - 1;  
@@ -401,8 +406,8 @@ template<
 INLINE_FUNCTION_HD
 int32 binarySearch(
 	const ViewType1D<Type, properties...>& view,
-	int32 start,
-	int32 end,
+	uint32 start,
+	uint32 end,
 	const Type& val)
 {
 	
@@ -424,10 +429,10 @@ template<
 	typename... dProperties>
 void exclusiveScan(
 	const ViewType1D<Type, properties...>& view,
-	int32 start,
-	int32 end,
+	uint32 start,
+	uint32 end,
 	ViewType1D<dType, dProperties...>& dView,
-	int32 dStart )
+	uint32 dStart )
 {
 
 	static_assert
@@ -441,7 +446,7 @@ void exclusiveScan(
 
 	using ExecutionSpace = typename ViewType1D<Type, properties...>::execution_space;
 	
-	int32 numElems = end-start;
+	uint32 numElems = end-start;
 	
 	pFlow::algorithms::KOKKOS::exclusiveScan<Type,dType,ExecutionSpace>(
 		view.data()+start,
@@ -457,10 +462,10 @@ template<
 	typename... dProperties>
 void inclusiveScan(
 	const ViewType1D<Type, properties...>& view,
-	int32 start,
-	int32 end,
+	uint32 start,
+	uint32 end,
 	ViewType1D<dType, dProperties...>& dView,
-	int32 dStart)
+	uint32 dStart)
 {
 	using ExecutionSpace = typename ViewType1D<Type, properties...>::execution_space;
 	
@@ -473,7 +478,7 @@ void inclusiveScan(
 	);
 
 
-	int32 numElems = end-start;
+	uint32 numElems = end-start;
 
 	pFlow::algorithms::KOKKOS::inclusiveScan<Type,dType,ExecutionSpace>(
 		view.data()+start,
diff --git a/src/phasicFlow/Kokkos/phsicFlowKokkos.hpp b/src/phasicFlow/Kokkos/phasicFlowKokkos.hpp
similarity index 97%
rename from src/phasicFlow/Kokkos/phsicFlowKokkos.hpp
rename to src/phasicFlow/Kokkos/phasicFlowKokkos.hpp
index 22666e49..3372ba8a 100644
--- a/src/phasicFlow/Kokkos/phsicFlowKokkos.hpp
+++ b/src/phasicFlow/Kokkos/phasicFlowKokkos.hpp
@@ -21,7 +21,7 @@ Licence:
 #ifndef __phsicFlowKokkos_hpp__
 #define __phsicFlowKokkos_hpp__
 
-#include "kokkosTypes.hpp"
+#include "KokkosTypes.hpp"
 #include "KokkosUtilities.hpp"
 #include "ViewAlgorithms.hpp"
 #include "Range.hpp"
diff --git a/src/phasicFlow/containers/List/ListPtr/ListPtr.hpp b/src/phasicFlow/containers/List/ListPtr/ListPtr.hpp
index 8d12d1a6..ed150eed 100644
--- a/src/phasicFlow/containers/List/ListPtr/ListPtr.hpp
+++ b/src/phasicFlow/containers/List/ListPtr/ListPtr.hpp
@@ -64,16 +64,16 @@ protected:
 		bool copy(const ListPtrType& src);
 
 		// - return ith pointer
-		T* ptr(label i);
+		T* ptr(size_t i);
 
 		// - return ith const poiter
-		const T* ptr(label i)const;
+		const T* ptr(size_t i)const;
 
 		// - iterator position of ith element
-		auto pos(label i);
+		auto pos(size_t i);
 		
 		// - const iterator position of ith element
-		auto pos(label i) const;
+		auto pos(size_t i) const;
 
 public:
 
@@ -141,15 +141,15 @@ public:
 	//// - Methods 
 
 		// - set the ith element 
-		T* set(label i, T* ptr);
+		T* set(size_t i, T* ptr);
 		
 		// - set the ith element and take the ownership from uniquePtr
-		uniquePtr<T> set(label i, uniquePtr<T>& ptr );
+		uniquePtr<T> set(size_t i, uniquePtr<T>& ptr );
 		
 		// - create the object in-place and set the pointer in ith position
 		//   if oject creation fails, uniquePtr deletes the memeory
 		template<typename... Args>
-	    uniquePtr<T> setSafe(label i, Args&&... args);
+	    uniquePtr<T> setSafe(size_t i, Args&&... args);
 	    
 	    // - put the pointer at the end
 	    void push_back(T* ptr);
@@ -163,11 +163,11 @@ public:
 	    
 		// - access to ith element
 		//   fatalexit if out of range or nullptr    
-		T& operator[](label i);
+		T& operator[](size_t i);
 		
 		// - const access to ith element
 		//   fatalexit if out of range or nullptr    
-		const T& operator[](label i) const;
+		const T& operator[](size_t i) const;
 
 		// size of container
 	    size_t size()const;
@@ -176,13 +176,13 @@ public:
 		auto empty() const;
 
 		// release the ownership of ith pointer 
-		uniquePtr<T> release(label i);
+		uniquePtr<T> release(size_t i);
 		
 		// - clear the content of list and delete objects
 		void clear();
 		
 		// - clear the ith element
-		void clear(label i);	
+		void clear(size_t i);	
 		
 		// - clone the object 
 		
diff --git a/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp b/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp
index 2d476154..57f42c2b 100644
--- a/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp
+++ b/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp
@@ -39,7 +39,7 @@ inline bool pFlow::ListPtr<T>::copy(const ListPtrType& src)
 }
 
 template<typename T>
-T* pFlow::ListPtr<T>::ptr(label i) 
+T* pFlow::ListPtr<T>::ptr(size_t i) 
 {
 
 	if(i >= size() ) return nullptr;	
@@ -53,7 +53,7 @@ T* pFlow::ListPtr<T>::ptr(label i)
 template<typename T>
 const T* pFlow::ListPtr<T>::ptr
 (
-	label i
+	size_t i
 ) const
 {
 
@@ -68,7 +68,7 @@ const T* pFlow::ListPtr<T>::ptr
 template<typename T>
 auto pFlow::ListPtr<T>::pos
 (
-	label i
+	size_t i
 )
 {
 	if(i >= size() ) 
@@ -86,7 +86,7 @@ auto pFlow::ListPtr<T>::pos
 template<typename T>
 auto pFlow::ListPtr<T>::pos
 (
-	label i
+	size_t i
 )const 
 {
 	if(i >= size() ) 
@@ -146,7 +146,7 @@ pFlow::ListPtr<T>& pFlow::ListPtr<T>::operator=
 template<typename T>
 T* pFlow::ListPtr<T>::set
 (
-	label i, T* ptr
+	size_t i, T* ptr
 )
 {
 	uniquePtr<T> uptr(ptr);
@@ -157,7 +157,7 @@ T* pFlow::ListPtr<T>::set
 template<typename T>
 pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
 (
-	label i,
+	size_t i,
 	uniquePtr<T>& ptr
 )
 {
@@ -182,7 +182,7 @@ template<typename T>
 template<typename... Args>
 pFlow::uniquePtr<T> pFlow::ListPtr<T>::setSafe
 (
-	label i,
+	size_t i,
 	Args&&... args
 )
 {
@@ -216,7 +216,7 @@ void pFlow::ListPtr<T>::push_backSafe(Args&&... args)
 template<typename T>
 T& pFlow::ListPtr<T>::operator[]
 (
-	label i
+	size_t i
 )
 {
 	T* p = ptr(i);
@@ -234,7 +234,7 @@ T& pFlow::ListPtr<T>::operator[]
 template<typename T>
 const T& pFlow::ListPtr<T>::operator[]
 (
-	label i
+	size_t i
 ) const
 {
 	const T* p = ptr(i);
@@ -263,7 +263,7 @@ auto pFlow::ListPtr<T>::empty() const
 template<typename T>
 pFlow::uniquePtr<T> pFlow::ListPtr<T>::release
 (
-	label i
+	size_t i
 )
 {
 	auto p = ptr(i);
@@ -288,7 +288,7 @@ void pFlow::ListPtr<T>::clear()
 template<typename T>
 void pFlow::ListPtr<T>::clear
 (
-	label i
+	size_t i
 )
 {
 	T* p = ptr(i);
diff --git a/src/phasicFlow/containers/Map/Map/Map.hpp b/src/phasicFlow/containers/Map/Map/Map.hpp
index 16c0b9d3..ffeef6b0 100644
--- a/src/phasicFlow/containers/Map/Map/Map.hpp
+++ b/src/phasicFlow/containers/Map/Map/Map.hpp
@@ -138,7 +138,7 @@ template<typename T>
 using wordMap = Map<word,T>;
 
 template<typename T>
-using labelMap = Map<label,T>;
+using uint64Map = Map<uint64,T>;
 
 template<typename T>
 using uint32Map = Map<uint32,T>;
@@ -155,7 +155,7 @@ template<typename T>
 inline iOstream& printKeys(iOstream& os, const wordMap<T> & m);
 
 template<typename T>
-inline iOstream& printKeys(iOstream& os, const labelMap<T> & m);
+inline iOstream& printKeys(iOstream& os, const uint64Map<T> & m);
 
 template<typename T>
 inline iOstream& printKeys(iOstream& os, const uint32Map<T> & m);
diff --git a/src/phasicFlow/containers/Map/Map/MapI.hpp b/src/phasicFlow/containers/Map/Map/MapI.hpp
index 7d52ea7e..7b28efcd 100644
--- a/src/phasicFlow/containers/Map/Map/MapI.hpp
+++ b/src/phasicFlow/containers/Map/Map/MapI.hpp
@@ -92,7 +92,7 @@ inline iOstream& printKeys(iOstream& os, const uint32Map<T> & m)
 }
 
 template<typename T>
-inline iOstream& printKeys(iOstream& os, const labelMap<T> & m)
+inline iOstream& printKeys(iOstream& os, const uint64Map<T> & m)
 {
 	if (m.empty())
 		return os<<"labelMap is empty"<<endl;
diff --git a/src/phasicFlow/containers/Map/hashMap/hashMap.hpp b/src/phasicFlow/containers/Map/hashMap/hashMap.hpp
index b6eabc69..676c7022 100644
--- a/src/phasicFlow/containers/Map/hashMap/hashMap.hpp
+++ b/src/phasicFlow/containers/Map/hashMap/hashMap.hpp
@@ -140,7 +140,7 @@ template<typename T>
 using wordHashMap = hashMap<word,T>;
 
 template<typename T>
-using labelHashMap = hashMap<label,T>;
+using uint64HashMap = hashMap<uint64,T>;
 
 template<typename T>
 using uint32HashMap = hashMap<uint32,T>;
@@ -155,7 +155,7 @@ template<typename T>
 inline iOstream& printKeys(iOstream& os, const wordHashMap<T> & m);
 
 template<typename T>
-inline iOstream& printKeys(iOstream& os, const labelHashMap<T> & m);
+inline iOstream& printKeys(iOstream& os, const uint64HashMap<T> & m);
 
 template<typename T>
 inline iOstream& printKeys(iOstream& os, const uint32HashMap<T> & m);
diff --git a/src/phasicFlow/containers/Map/hashMap/hashMapI.hpp b/src/phasicFlow/containers/Map/hashMap/hashMapI.hpp
index 74e85e6d..4c45687a 100644
--- a/src/phasicFlow/containers/Map/hashMap/hashMapI.hpp
+++ b/src/phasicFlow/containers/Map/hashMap/hashMapI.hpp
@@ -80,7 +80,7 @@ inline iOstream& printKeys(iOstream& os, const wordHashMap<T> & m)
 
 
 template<typename T>
-inline iOstream& printKeys(iOstream& os, const labelHashMap<T> & m)
+inline iOstream& printKeys(iOstream& os, const uint64HashMap<T> & m)
 {
 	if (m.empty())
 		return os<<"labelHashMap is empty"<<endl;
diff --git a/src/phasicFlow/containers/Vector/Vector.cpp b/src/phasicFlow/containers/Vector/Vector.cpp
index 7c0c9db3..71dc312a 100644
--- a/src/phasicFlow/containers/Vector/Vector.cpp
+++ b/src/phasicFlow/containers/Vector/Vector.cpp
@@ -19,12 +19,6 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 
-template<typename T, typename Allocator>
-pFlow::Vector<T, Allocator>::Vector(iIstream& is)
-{
-	readVector(is);
-}
-
 template<typename T, typename Allocator>
 bool pFlow::Vector<T, Allocator>::readVector
 (
@@ -33,7 +27,6 @@ bool pFlow::Vector<T, Allocator>::readVector
 )
 {
 	
-
 	if(is.isBinary() && !std::is_same_v<T,word>)
 	{
 		this->resize(len);
@@ -90,8 +83,6 @@ bool pFlow::Vector<T, Allocator>::readVector
 		}
 	}
 
-	
-
 	return true;
 }
 
@@ -102,8 +93,12 @@ bool pFlow::Vector<T, Allocator>::writeVector
 	iOstream& os
 ) const
 {
+	
+	span<T> s( const_cast<T*>(this->data()), this->size());
+	os<<s;
+
 	// start of 
-	if( os.isBinary() && !std::is_same_v<T,word>)
+	/*if( os.isBinary() && !std::is_same_v<T,word>)
 	{
 		os.write(reinterpret_cast<const char*>(this->data()), this->size()*sizeof(T));
 	}
@@ -113,12 +108,12 @@ bool pFlow::Vector<T, Allocator>::writeVector
 		auto len = size();
 		auto stride  = getVectorStride(len);
 		os << token::BEGIN_LIST;
-		label i = 0;
+		size_t i = 0;
 		while( i<len )
 		{
 
 			os << this->operator[](i++);
-			for(label j=0; j<stride-1 && i<len; j++ )
+			for(size_t j=0; j<stride-1 && i<len; j++ )
 			{
 				os << token::SPACE << this->operator[](i++);	
 			}
@@ -130,12 +125,12 @@ bool pFlow::Vector<T, Allocator>::writeVector
 	    os << token::END_LIST;
 
 	    os.check(FUNCTION_NAME);
-	}
+	}*/
 
     return true;
 }
 
-template<typename T, typename Allocator>
+/*template<typename T, typename Allocator>
 bool pFlow::Vector<T, Allocator>::deleteElement_sorted
 (
 	const Vector<label>& indices
@@ -384,4 +379,4 @@ inline bool pFlow::Vector<T, Allocator>::insertSetElement
 	}
 	
 	return true;
-}
\ No newline at end of file
+}*/
\ No newline at end of file
diff --git a/src/phasicFlow/containers/Vector/Vector.hpp b/src/phasicFlow/containers/Vector/Vector.hpp
index c8d354ee..44e3beb3 100644
--- a/src/phasicFlow/containers/Vector/Vector.hpp
+++ b/src/phasicFlow/containers/Vector/Vector.hpp
@@ -29,7 +29,7 @@ Licence:
 #include "error.hpp"
 #include "uniquePtr.hpp"
 #include "stdAlgorithms.hpp"
-#include "indexContainer.hpp"
+#include "span.hpp"
 #include "iOstream.hpp"
 #include "iIstream.hpp"
 
@@ -93,16 +93,6 @@ protected:
 	// - name of the vector  
 	word name_;
 
-	static inline size_t getVectorStride(const size_t& len)
-	{
-		size_t stride = 1;
-		if( len < 6 ) 		stride = len;
-		else if( len <16 )  stride = 3;
-		else if( len < 31)  stride = 2;
-		else				stride = 1;
-
-		return stride;
-	}
 
 	static constexpr bool  isHostAccessible_ = true;
 
@@ -119,37 +109,28 @@ public:
 
 	//// - Constructors
 
-		// - empty Vector
+		/// Empty Vector
 		inline Vector()
 		:
 			Vector("Vector")
 		{}
 
+
+		/// Empty Vector with a name 
 		inline Vector(const word& name)
 		:
 			name_(name)
 		{}
-		// - with sepcified length
-		inline Vector(const size_t len)
-		:
-			Vector("Vector",len)
-		{}
 
-		// - with specified length and name
+		
+		/// Vector with specified length and name
 		inline Vector(const word& name, size_t len)
 		:
 			vectorType(len),
 			name_(name)
-		{
-			
-		}
-
-		// - with length and value
-		inline Vector(size_t len, const T& val)
-		:
-			Vector("Vector", len, val)
 		{}
 
+		/// Vector with name, length and value
 		inline Vector(const word& name, size_t len, const T& val)
 		:
 			Vector(name, len)
@@ -157,20 +138,7 @@ public:
 			this->assign(len, val);
 		}
 
-		// - zero length with specified capacity, use Logical
-		//   to make it different from previous constructor.
-		inline Vector(const size_t cap, RESERVE ):
-			Vector("Vector", cap, 0, RESERVE())
-		{
-		}
-
-		inline Vector(const size_t cap, const size_t len, RESERVE )
-		:
-			Vector("Vector", cap, len, RESERVE())
-		{
-			
-		}
-
+		/// Vector with name, size and reserved capacity 
 		Vector(const word& name, size_t cap, size_t len, RESERVE ):
 			name_(name)
 		{
@@ -178,142 +146,157 @@ public:
 			this->resize(len);
 		}
 
-		inline Vector(const size_t cap, const size_t len, const T& val, RESERVE )
-		{
-			name_ = "Vector";
-			reserve(cap);
-			this->assign(len, val);
-		}
 
-
-		// from initializer list 
-		inline Vector(const initList &l)
+		/// Vector from name and initializer list 
+		inline Vector(const word& name, const initList &l)
 		:
-			vectorType(l)
+			vectorType(l),
+			name_(name)
 		{}
 
-		// - from src and a new name
+		/// Construct with a name and form std::vector (host memory)
+		inline Vector(const word& name, const vectorType& src)
+		:
+			vectorType(src),
+			name_(name)
+		{}
+
+		/// Construct with a name and form std::vector (host memory)
+		/// and with a desired capacity. 
+		inline Vector(const word& name, const vectorType& src, size_t cap)
+		:
+			Vector(name, cap, src.size(), RESERVE())
+		{
+			this->assign(src.begin(), src.end());
+		}
+
+		/// Copy construct 
+		inline Vector(const VectorType& src) = default;
+
+		/// Copy from src with a new name 
 		inline Vector(const word name, const Vector<T>& src):
 			vectorType(src),
 			name_(name)
 		{}
 
-		// copy construct 
-		inline Vector(const VectorType& src) = default;
-		
-		// move construct 
-		inline Vector( VectorType && mv) = default;
-
-		inline Vector(const vectorType& src)
-		:
-			vectorType(src),
-			name_("Vector")
-		{
-			
-		}
-		
-		// copy assignment
+		/// Copy assignment
 		inline VectorType& operator=( const VectorType& rhs ) = default;
 
+		/// Copy assignment from std::vector 
 		inline VectorType& operator=(const vectorType& rhs)
 		{
 			Vector::assign(rhs.begin(), rhs.end());
 			return *this;
 		}
-		
-		// move assignment 
+
+		/// Move construct 
+		inline Vector( VectorType && mv) = default;
+	
+		/// Move assignment 
 		inline VectorType& operator=( VectorType && mv) = default;
 	
-		// scalar assignment 
+		/// Scalar assignment 
 		inline void operator=(const T& val)
 		{
 			fill(val);
 		}
 
+		/// Destructor 
 		inline ~Vector()
 		{
 			vectorType::clear();
 		}
 
+		/// Clone as a uniquePtr
 		inline uniquePtr<VectorType> clone() const
 		{
 			return makeUnique<VectorType>(*this);
 		}
 
+		/// Clone as a pointer 
 		inline VectorType* clonePtr()const
 		{
 			return new VectorType(*this);
 		}
 
-	inline auto clear()
-	{
-		return vectorType::clear();
-	}
+	
+	//// - Methods 
 
-	// access to this, mostly used by derived classes 
-	const VectorType& VectorField() const
-	{
-		return *this;
-	}
- 
-	VectorType& VectorField()
-	{
-		return *this;
-	}
+		/// Access to this, mostly used by derived classes 
+		const VectorType& VectorField() const
+		{
+			return *this;
+		}
+	 
+		VectorType& VectorField()
+		{
+			return *this;
+		}
 
-	const vectorType& vectorField()const
-	{
-		return *this;
-	}
+		const vectorType& vectorField()const
+		{
+			return *this;
+		}
 
-	vectorType& vectorField()
-	{
-		return *this;
-	}
+		vectorType& vectorField()
+		{
+			return *this;
+		}
 
-	auto& deviceVectorAll()
-	{
-		return *this;
-	}
+		auto& deviceVectorAll()
+		{
+			return *this;
+		}
 
-	const auto& deviceVectorAll()const
-	{
-		return *this;
-	}
+		const auto& deviceVectorAll()const
+		{
+			return *this;
+		}
 
-	auto& deviceVector()
-	{
-		return *this;
-	}
+		auto& deviceVector()
+		{
+			return *this;
+		}
 
-	const auto& deviceVector()const
-	{
-		return *this;
-	}
+		const auto& deviceVector()const
+		{
+			return *this;
+		}
 
 
+		/// Name of the vector 
+		const word& name()const
+		{
+			return name_;
+		}
 
-	const word& name()const
-	{
-		return name_;
-	}
+		/// Size of the vector 
+		inline auto size()const
+		{
+			return vectorType::size();
+		}
 
-	inline auto size()const
-	{
-		return vectorType::size();
-	}
+		/// Capacity of the vector 
+		inline auto capacity()const
+		{
+			return vectorType::capacity();
+		}
 
-	inline auto capacity()const
-	{
-		return vectorType::capacity();
-	}
+		/// If vector is empty 
+		inline bool empty()const
+		{
+			return vectorType::empty();
+		}
 
-	inline auto reserve(label len)
-	{
-		return vectorType::reserve(len);
-	}
+		/// Reserve capacity for vector 
+		/// Preserve the content.
+		inline void reserve(size_t cap)
+		{
+			vectorType::reserve(cap);
+		}
 
-	// - delete elemens of vector based on sorted indices 
+
+	/*// - delete elemens of vector based on sorted indices 
 	//   return false if out of range
 	bool deleteElement_sorted(const Vector<label>& indices );
 
@@ -346,7 +329,7 @@ public:
 
 	// - set or insert a new element into the vecor 
 	//   return false if it fails 
-	inline bool insertSetElement(int32 idx, const T& val);
+	inline bool insertSetElement(int32 idx, const T& val);*/
 
 	// - fill the whole content of vector, [begin, end), with val 
 	inline void fill( const T& val);
@@ -368,12 +351,7 @@ public:
 
 	inline VectorType operator -()const;
 
-	// from iIstream and specified size 
-	//Vector(iIstream & is, size_t len);
-
-	// from iIstream and free size
-	Vector(iIstream& is);
-
+	
 	bool readVector(iIstream& is, size_t len=0);
 
 	bool writeVector(iOstream& os) const;
diff --git a/src/phasicFlow/containers/Vector/VectorI.hpp b/src/phasicFlow/containers/Vector/VectorI.hpp
index 6bacd5dd..c4426169 100644
--- a/src/phasicFlow/containers/Vector/VectorI.hpp
+++ b/src/phasicFlow/containers/Vector/VectorI.hpp
@@ -77,7 +77,7 @@ inline void pFlow::Vector<T, Allocator>::operator +=( const Vector<T, Allocator>
 		fatalExit;
 	}
 
-	for(label i=0; i<v.size(); i++)
+	for(size_t i=0; i<v.size(); i++)
 	{
 		this->operator[](i) += v[i];
 	}
@@ -96,7 +96,7 @@ inline void pFlow::Vector<T, Allocator>::operator -=( const Vector<T, Allocator>
 		fatalExit;
 	}
 
-	for(label i=0; i<v.size(); i++)
+	for(size_t i=0; i<v.size(); i++)
 	{
 		this->operator[](i) -= v[i];
 	}
@@ -115,7 +115,7 @@ inline void pFlow::Vector<T, Allocator>::operator *=( const Vector<T, Allocator>
 		fatalExit;
 	}
 
-	for(label i=0; i<v.size(); i++)
+	for(size_t i=0; i<v.size(); i++)
 	{
 		this->operator[](i) *= v[i];
 	}
@@ -134,7 +134,7 @@ inline void pFlow::Vector<T, Allocator>::operator /=( const Vector<T, Allocator>
 		fatalExit;
 	}
 
-	for(label i=0; i<v.size(); i++)
+	for(size_t i=0; i<v.size(); i++)
 	{
 		this->operator[](i) /= v[i];
 	}
diff --git a/src/phasicFlow/containers/Vector/VectorMath.hpp b/src/phasicFlow/containers/Vector/VectorMath.hpp
index 22c7087d..153fdfff 100644
--- a/src/phasicFlow/containers/Vector/VectorMath.hpp
+++ b/src/phasicFlow/containers/Vector/VectorMath.hpp
@@ -18,12 +18,14 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+namespace pFlow
+{
 
 #define VecFunc(fnName) 									\
 template<typename T, typename Allocator>										\
-inline pFlow::Vector<T, Allocator> pFlow::fnName(const Vector<T,Allocator>& v)	\
+inline Vector<T, Allocator> fnName(const Vector<T,Allocator>& v)	\
 {															\
-	Vector<T, Allocator> res(v.capacity(), Logical());		\
+	Vector<T, Allocator> res(v.name(), v.capacity(), 0 ,RESERVE());		\
 	for(auto& e:v)											\
 	{														\
 		res.push_back( fnName(e) );							\
@@ -31,10 +33,10 @@ inline pFlow::Vector<T, Allocator> pFlow::fnName(const Vector<T,Allocator>& v)	\
 	return std::move(res);									\
 }															\
 template<typename T, typename Allocator, typename indexFunc>					\
-inline pFlow::Vector<T, Allocator> pFlow::fnName(const Vector<T, Allocator>& v, indexFunc iFn) 	\
+inline Vector<T, Allocator> fnName(const Vector<T, Allocator>& v, indexFunc iFn) 	\
 {																		\
-	Vector<T, Allocator> res(v.capacity(), Logical());								\
-	for(label i=0; i<v.size(); i++)										\
+	Vector<T, Allocator> res(v.name(), v.capacity(), 0, RESERVE());								\
+	for(size_t i=0; i<v.size(); i++)										\
 	{																	\
 		if( iFn(i) )													\
 			res.push_back(fnName(v[i]));								\
@@ -46,20 +48,20 @@ inline pFlow::Vector<T, Allocator> pFlow::fnName(const Vector<T, Allocator>& v,
 
 #define VecFunc2(fnName) 									\
 template<typename T, typename Allocator>										\
-inline pFlow::Vector<T, Allocator> pFlow::fnName(const Vector<T, Allocator>& v1, const Vector<T, Allocator>& v2)	\
+inline Vector<T, Allocator> fnName(const Vector<T, Allocator>& v1, const Vector<T, Allocator>& v2)	\
 {															\
-	Vector<T, Allocator> res(v1.capacity(), Logical());					\
-	for(label i=0; i<v1.size(); i++)							\
+	Vector<T, Allocator> res(v1.name(), v1.capacity(), 0 ,RESERVE());					\
+	for(size_t i=0; i<v1.size(); i++)							\
 	{														\
 		res.push_back( fnName(v1[i], v2[i]));				\
 	}														\
 	return std::move(res);									\
 }															\
 template<typename T, typename Allocator, typename indexFunc>					\
-inline pFlow::Vector<T, Allocator> pFlow::fnName(const Vector<T, Allocator>& v1, const Vector<T, Allocator>& v2, indexFunc iFn) 	\
+inline Vector<T, Allocator> fnName(const Vector<T, Allocator>& v1, const Vector<T, Allocator>& v2, indexFunc iFn) 	\
 {																		\
-	Vector<T, Allocator> res(v1.capacity(), Logical());							\
-	for(label i=0; i<v1.size(); i++)									\
+	Vector<T, Allocator> res(v1.name(), v1.capacity(), 0 ,RESERVE());							\
+	for(size_t i=0; i<v1.size(); i++)									\
 	{																	\
 		if( iFn(i) )													\
 			res.push_back(fnName(v1[i], v2[i]));						\
@@ -77,6 +79,8 @@ inline pFlow::Vector<T, Allocator> pFlow::fnName(const Vector<T, Allocator>& v1,
 // min, max
 //* * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+
+
 VecFunc(abs);
 VecFunc2(mod);
 VecFunc(exp);
@@ -100,6 +104,9 @@ VecFunc(acosh);
 VecFunc(atanh);
 
 #undef VecFunc
+}
+
+
 
 //// special implementation of math functions 
 namespace pFlow
@@ -108,7 +115,7 @@ namespace pFlow
 template<typename T, typename Allocator>
 inline Vector<T, Allocator> pow(const Vector<T, Allocator>& v, T e)
 {
-	Vector<T, Allocator> res(v.capacity(), Logical());
+	Vector<T, Allocator> res(v.name(), v.capacity(), 0 ,RESERVE());
 	for(auto& elm:v)
 	{
 		res.push_back(pow(elm,e));
@@ -119,8 +126,8 @@ inline Vector<T, Allocator> pow(const Vector<T, Allocator>& v, T e)
 template<typename T, typename Allocator, typename indexFunc>
 inline Vector<T,Allocator> pow(const Vector<T,Allocator>& v, T e, indexFunc iFn)
 {
-	Vector<T, Allocator> res(v.capacity(), Logical());
-	for(label i=0; i<v.size(); i++)
+	Vector<T, Allocator> res(v.name(), v.capacity(), 0 ,RESERVE());
+	for(size_t i=0; i<v.size(); i++)
 	{
 		if(iFn(i))
 		{
@@ -149,7 +156,7 @@ template<typename T, typename Allocator, typename indexFunc>
 inline T min(const Vector<T, Allocator>& v, indexFunc iFn)
 {
 	T minVal(largestPositive<T>());
-	for(label i=0; i<v.size(); i++)
+	for(size_t i=0; i<v.size(); i++)
 	{
 		if(iFn(i))
 		{
@@ -175,7 +182,7 @@ template<typename T, typename Allocator ,typename indexFunc>
 inline T max(const Vector<T, Allocator>& v, indexFunc iFn)
 {
 	T maxVal(largestNegative<T>());
-	for(label i=0; i<v.size(); i++)
+	for(size_t i=0; i<v.size(); i++)
 	{
 		if(iFn(i))
 		{
@@ -201,7 +208,7 @@ template<typename T, typename Allocator, typename indexFunc>
 inline T sum(const Vector<T, Allocator>& v, indexFunc iFn)
 {
 	T s = static_cast<T>(0);
-	for(label i=0; i<v.size(); ++i)
+	for(size_t i=0; i<v.size(); ++i)
 	{
 		if(iFn(i))
 			s += v[i];
diff --git a/src/phasicFlow/containers/Vector/Vectors.cpp b/src/phasicFlow/containers/Vector/Vectors.cpp
index dd958b08..43d2952f 100644
--- a/src/phasicFlow/containers/Vector/Vectors.cpp
+++ b/src/phasicFlow/containers/Vector/Vectors.cpp
@@ -21,25 +21,20 @@ Licence:
 #include "Vectors.hpp"
 
 // instantiation just for numeral types 
-template class pFlow::Vector<pFlow::int8>;
 
-template class pFlow::Vector<pFlow::int16>;
+template class pFlow::Vector<pFlow::int8>;
 
 template class pFlow::Vector<pFlow::int32>;
 
-template class pFlow::Vector<pFlow::int64>;
-
 template class pFlow::Vector<pFlow::uint32>;
 
-template class pFlow::Vector<pFlow::label>;
-
 template class pFlow::Vector<pFlow::real>;
 
 template class pFlow::Vector<pFlow::realx3>;
 
 template class pFlow::Vector<pFlow::realx3x3>;
 
-//template class pFlow::Vector<pFlow::word>;
+
 
  
 
diff --git a/src/phasicFlow/containers/Vector/Vectors.hpp b/src/phasicFlow/containers/Vector/Vectors.hpp
index 8eba5863..523444e4 100644
--- a/src/phasicFlow/containers/Vector/Vectors.hpp
+++ b/src/phasicFlow/containers/Vector/Vectors.hpp
@@ -31,69 +31,25 @@ namespace pFlow
 
 using int8Vector 	= Vector<int8>;
 
-using int16Vector 	= Vector<int16>;
-
 using int32Vector 	= Vector<int32>;
 
 using int64Vector 	= Vector<int64>;
 
-using uint32Vector 	= Vector<uint32>;
+using uint8Vector  = Vector<uint8>;
 
-using labelVector 	= Vector<label>;
+using uint32Vector   = Vector<uint32>;
+
+using uint64Vector   = Vector<uint64>;
 
 using realVector 	= Vector<real> 	;
 
 using realx3Vector 	= Vector<realx3>;
 
-using uint16x3Vector= Vector<uint16x3>;
-
-using uint32x3Vector= Vector<uint32x3>;
-
-using int32x3Vector = Vector<int32x3>;
-
-using int64x3Vector = Vector<int64x3>;
-
-using uint16x3x3Vector 	= Vector<uint16x3x3>;
-
-using uint32x3x3Vector 	= Vector<uint32x3x3>;
-
-using int32x3x3Vector 	= Vector<int32x3x3>;
-
 using realx3x3Vector 	= Vector<realx3x3>;
 
 using wordVector 		= Vector<word>;
 
-/*template<>
-inline word Vector<label>::TYPENAME() { return "labelVector"; }
-typedef Vector<label> 	labelVector;
 
-template<>
-inline word Vector<real>::TYPENAME() { return "realVector"; }
-typedef Vector<real> 	realVector;
-
-template<>
-inline word Vector<unit3>::TYPENAME() { return "unit3Vector"; }
-typedef Vector<unit3> 	unit3Vector;
-
-template<>
-inline word Vector<real3>::TYPENAME() { return "real3Vector"; }
-typedef Vector<real3>	real3Vector;
-
-template<>
-inline word Vector<real33>::TYPENAME() { return "real33Vector"; }
-typedef Vector<real33>	real33Vector;
-
-template<>
-inline word Vector<bool>::TYPENAME() { return "boolVector"; }
-typedef Vector<bool> 	boolVector;
-
-template<>
-inline word Vector<word>::TYPENAME() { return "wordVector"; }
-typedef Vector<word> 	wordVector;
-
-template<>
-inline word Vector<sint>::TYPENAME() { return "sintVector"; }
-typedef Vector<sint> 	sintVector;*/
 
 }
 
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
index 8dafeeb7..cf1069fd 100644
--- a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
+++ b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
@@ -22,16 +22,20 @@ Licence:
 #ifndef __VectorSingle_hpp__
 #define __VectorSingle_hpp__
 
+#include <vector>
+
+
 #include "globalSettings.hpp"
-#include "types.hpp"
 #include "typeInfo.hpp"
-#include "Vector.hpp"
+#include "types.hpp"
+#include "error.hpp"
 #include "indexContainer.hpp"
+#include "iOstream.hpp"
+#include "iIstream.hpp"
+#include "span.hpp"
+#include "Vector.hpp"
+#include "phasicFlowKokkos.hpp"
 
-#include "streams.hpp"
-
-#include "KokkosTypes.hpp"
-#include "ViewAlgorithms.hpp"
 
 
 #ifndef __RESERVE__
@@ -42,103 +46,131 @@ Licence:
 namespace pFlow
 {
 
-
+//- Forward 
 template<typename T, typename MemorySpace>
 class VectorSingle;
 
 
+
 template<typename T, typename MemorySpace=void>
 class VectorSingle
 {
 public:
 	
-	// viewType (view of data host and device)
-  	using VectorType		= VectorSingle<T, MemorySpace>;
+	//- typedefs for accessing data
 
-  	using iterator        = T*;
+		using VectorType		= VectorSingle<T, MemorySpace>;
 
-  	using constIterator   = const T*;
+		using iterator        = T*;
+
+		using constIterator   = const T*;
+	  	
+		using reference       = T&;
+	  	
+	 	using constReference  = const T&;
+
+		using valueType       = T;
+	  	
+		using pointer         = T*;
+	  	
+		using constPointer    = const T*;
   	
-	using reference       = T&;
-  	
-  	using constReference  = const T&;
+	//- typedefs related to memory management 
 
-	using valueType       = T;
-  	
-  	using pointer         = T*;
-  	
-  	using constPointer    = const T*;
-  	
-  	// type defs related to Kokkos 
-  	using viewType 			= ViewType1D<T, MemorySpace>;
+		using viewType 			= ViewType1D<T, MemorySpace>;
 
-  	using deviceType 		= typename viewType::device_type;
+		using deviceType 		= typename viewType::device_type;
 
-  	using memory_space 		= typename viewType::memory_space;
+		using memorySpace 		= typename viewType::memory_space;
 
-  	using execution_space 	= typename viewType::execution_space;
+		using executionSpace 	= typename viewType::execution_space;
 
 protected:
 	
-	size_t 			size_ = 0;
+	// - Data members
 
-	size_t 			capacity_ = 0;
+		/// Size of the vector 
+		uint32 			size_ = 0;
 
-	viewType   		view_;
+		/// view of the vector 
+		viewType   		view_;
 
-	mutable viewType	subView_;
+	// - protected members and methods 
 
-	mutable bool		subViewUpdated_ = false;
+		/// growth factor for vector 
+		static const inline 
+		real 	growthFactor_ = vectorGrowthFactor__;
 
-	static const inline real 	growthFactor_ = vectorGrowthFactor__;
+		/// Is the memory of this vector accessible from Host
+		static constexpr 
+		bool  isHostAccessible_ = isHostAccessible<executionSpace>();
+			//Kokkos::SpaceAccessibility<executionSpace,HostSpace>::accessible;
 
-	static constexpr bool  isHostAccessible_ =
-  	Kokkos::SpaceAccessibility<execution_space,Kokkos::HostSpace>::accessible;
+		/// Is the memory of this vector accessiple from Divce
+		static constexpr
+		bool isDeviceAccessible_ = isDeviceAccessible<executionSpace>();
 
-  	constexpr static inline const char* memoerySpaceName()
-  	{
-  		return memory_space::name();
-  	}
 
-	static INLINE_FUNCTION_H size_t evalCapacity(size_t n)
+		/// Name of the memory space
+	  	constexpr static inline 
+	  	const char* memoerySpaceName()
+	  	{
+	  		return memorySpace::name();
+	  	}
+
+	  	/// Evaluate capacity based on the input size 
+		static INLINE_FUNCTION_H uint32 evalCapacity(uint32 n)
+		{
+			return static_cast<uint32>(n*growthFactor_+1);
+		}
+
+	/// Change size to n and preserve the conetent if realloc 
+	/// occurs 
+	INLINE_FUNCTION_H 
+	uint32 changeSize(uint32 n, bool withInit= false)
 	{
-		return static_cast<size_t>(n*growthFactor_+1);
+		if(n > this->capacity() )
+		{
+			uint32 newCap = evalCapacity(n);
+			changeCapacity(newCap, withInit);
+		}
+		return setSize(n);
 	}
 
-	// use actualCap = true only for reserve
-	INLINE_FUNCTION_H void changeSize(size_t n, bool actualCap=false)
+	INLINE_FUNCTION_H
+	uint32 changeCapacitySize(uint32 actualCap, uint32 n, bool withInit= false)
 	{
-		if(n > capacity_ )
-		{
-			if(actualCap)
-				capacity_ = n;
-			else
-				capacity_ = evalCapacity(n);
+		changeCapacity(actualCap, withInit);
+		return setSize(n);
+	}
 
-			Kokkos::resize(view_, capacity_);
-			subViewUpdated_ = false;
-		}
-		if(!actualCap) 
+	INLINE_FUNCTION_H
+	void changeCapacity(uint32 actualCap, bool withInit= false)
+	{
+		if(withInit)
 		{
-			setSize(n);
+			resizeInit(view_, actualCap);
+		}
+		else
+		{
+			resizeNoInit(view_, actualCap);
 		}
 	}
 
-	INLINE_FUNCTION_H void setSize(size_t n)
+	INLINE_FUNCTION_H
+	uint32 reallocateCapacitySize(uint32 cap, uint32 s)
 	{
+		reallocNoInit(view_, cap);
+		return setSize(s);
+	}
+
+	INLINE_FUNCTION_H 
+	uint32 setSize(uint32 n)
+	{
+		auto os = size_;
 		size_ = n;
-		subViewUpdated_ = false;
+		return os;
 	}
-
-	// - update subview 
-	INLINE_FUNCTION_H void updateSubView()const
-	{
-		if(subViewUpdated_) return;
-
-		subView_ 	= Kokkos::subview(view_, Kokkos::make_pair(0,int(size_)));
-		subViewUpdated_ = true;
-	}
-	
 	
 public:
 
@@ -147,567 +179,301 @@ public:
 
 	//// - Constructors
 
-		// - empty constructor
+		/// Empty vector
 		VectorSingle()
 		:
 			VectorSingle("VectorSingle")
 		{}
 
-		// empty vector with a name 
+		/// Empty vector with a name 
 		VectorSingle(const word& name)
 		:
 			size_(0),
-			capacity_(2),
-			view_(name,capacity_)
-		{
-			changeSize(size_);
-		}
-
-		// - a vector with size n
-		VectorSingle(size_t n)
-		:
-			VectorSingle("VectorSingle",n)
+			view_(name,2)
 		{}
 
-		// - a vector with name and size n
-		VectorSingle(const word& name, size_t n)
+		
+		/// Vector with name and size n
+		VectorSingle(const word& name, uint32 n)
 		:
 			size_(n),
-			capacity_(evalCapacity(n)),
-			view_(name, capacity_)
-		{
-			changeSize(size_);
-		}
-
-		// a vector with size and value 
-		VectorSingle(size_t n, const T& val)
-		:
-			VectorSingle("VectorSingle", n , val)
+			view_(name, evalCapacity(n))
 		{}
 
-		// a vector with name, size and value
-		VectorSingle(const word& name, size_t n, const T& val)
+		
+
+		/// Vector with name, size and value
+		VectorSingle(const word& name, uint32 n, const T& val)
 		:
 			VectorSingle(name, n)
 		{
-			assign(n, val);
+			assign(n, val); 
 		}
 
-		// a vector with name and reserved capacity 
-		VectorSingle(size_t cap, size_t n, RESERVE )
+
+		/// Vector with name, size (n) and reserved capacity 
+		VectorSingle(const word& name, uint32 cap, uint32 n, RESERVE )
 		:
-			VectorSingle("VectorSingle", cap, n, RESERVE())
+			size_(n),
+			view_(name, cap)
 		{}
 
-
-		// a vector with name and reserved capacity 
-		VectorSingle(const word& name, size_t cap, size_t n, RESERVE )
-		:
-			VectorSingle(name)
-		{
-			reallocate(cap);
-			size_ = n;
-		}
-
-		// - construct from pFlow::Vector (host memory)
-		VectorSingle(const Vector<T> & src)
-		:
-			VectorSingle("VectorSingle", src)
-		{}
-
-		// - construct from pFlow::Vector and name 
-		VectorSingle(const word& name, const Vector<T> & src)
+		/// Construct with a name and form std::vector (host memory)
+		VectorSingle(const word& name, const std::vector<T> & src) 
 		:
 			VectorSingle(name)
 		{
 			assign(src);
 		}
 
-		// - copy construct (perform deep copy)
-		VectorSingle(const VectorSingle& src)
+		/// Construct with a name and form std::vector (host memory)
+		/// and with a desired capacity. 
+		VectorSingle(const word& name, const std::vector<T> & src, uint32 cap) 
 		:
-			VectorSingle(src.name(), src.capacity(), src.size(), RESERVE())
-		{	
-			copy(deviceVectorAll(), src.deviceVectorAll());
-			////Kokkos::deep_copy(deviceVectorAll(), src.deviceVectorAll());
+			VectorSingle(name)
+		{
+			assign(src, cap);
 		}
 
-		// - copy construct with a new name 
+		/// Copy construct (performs deep copy)
+		VectorSingle(const VectorSingle& src)
+		:
+			VectorSingle(src.name(), src)
+		{	
+			copy(deviceVector(), src.deviceVector());
+		}
+
+		/// Copy construct with a new name (perform deep copy)
 		VectorSingle(const word& name, const VectorSingle& src)
 		:
 			VectorSingle(name, src.capacity(), src.size(), RESERVE())
 		{
-			copy(deviceVectorAll(), src.deviceVectorAll());	
-			////Kokkos::deep_copy(deviceVectorAll(), src.deviceVectorAll());
+			copy(deviceVector(), src.deviceVector());	
 		}
 
-		// - copy assignment 
-		VectorSingle& operator = (const VectorSingle& rhs)
+		/// Copy assignment (perform deep copy from rhs to *this)
+		VectorSingle& operator = (const VectorSingle& rhs) 
 		{
 			if(&rhs == this) return *this;
 			VectorSingle temp(rhs);
-			capacity_   = temp.capacity();
-			size_ 		= temp.size();	
+
+			size_ 	= temp.size();	
 			view_   = temp.view_;
-			subViewUpdated_ = false;
 
 			return *this;
 		}
 
-		// no move construct
-		VectorSingle(VectorSingle&&) = delete;
+		/// Move construct
+		VectorSingle(VectorSingle&&) = default;   
 
-		// no move assignment
-		VectorSingle& operator= (VectorSingle&&) = delete;
+		/// Move assignment
+		VectorSingle& operator= (VectorSingle&&) = default; 
 
-
-		// - clone as a uniquePtr
-		INLINE_FUNCTION_H 
-		uniquePtr<VectorSingle> clone() const
+		/// This may not destroy the underlying memory, sice view is 
+		/// shared_ptr and maybe referenced by another object too 
+		~VectorSingle()
 		{
-			return makeUnique<VectorSingle>(*this);
+			view_ = viewType();
+			size_ = 0;
 		}
 
-		// - clone as a pointer
-		INLINE_FUNCTION_H
-		VectorSingle* clonePtr()const
+		/// Clone as a uniquePtr (perform deep copy)
+		INLINE_FUNCTION_H 
+		uniquePtr<VectorSingle> clone() const  
 		{
-			return new VectorSingle(*this);
+			return makeUnique<VectorSingle>( this->name(), *this);
+		}
+
+		/// Clone as a pointer (perform deep copy)
+		INLINE_FUNCTION_H
+		VectorSingle* clonePtr()const   
+		{
+			return new VectorSingle(this->name(), *this);
 		}
 
 	//// - Methods
 
-		// - return *this
+		/// Return *this
 		INLINE_FUNCTION_H
 		VectorType& VectorField()
 		{
 			return *this;
 		}
 
-		// - return *this
+		/// Return *this
 		INLINE_FUNCTION_H
 		const VectorType& VectorField()const
 		{
 			return *this;
 		}
 
-
-		// - Device vector range [0,capcity)
+		/// Device vector range [0,capcity)
 		INLINE_FUNCTION_H 
-		viewType& deviceVectorAll(){
+		viewType deviceVectorAll() const{
 			return view_;
 		}
 
-		// - Device vector range [0,capacity)
+		///  Device vector range [0, size)
 		INLINE_FUNCTION_H 
-		const viewType& deviceVectorAll() const {
-			return view_;
+		viewType deviceVector()const
+		{
+			return Kokkos::subview(view_, Kokkos::make_pair(0,int(size_)));
 		}
 
-		//  - Device vector range [0, size)
+		/// Return a vector accessible on Host in range [0,capacity)
 		INLINE_FUNCTION_H 
-		viewType& deviceVector(){
-			updateSubView();
-			return subView_;
-		}
-
-		//  - Device vector range [0, size)
-		INLINE_FUNCTION_H 
-		const viewType& deviceVector()const{
-			updateSubView();
-			return subView_;
-		}
-
-		INLINE_FUNCTION_H 
-		const auto hostVectorAll()const
+		auto hostVectorAll()const
 		{
 			auto hView = Kokkos::create_mirror_view(view_);
 			copy(hView, view_);
 			return hView;
 		}
 
+		/// Return a vector accessible on Host in range [0,size)
 		INLINE_FUNCTION_H 
-		auto hostVectorAll()
-		{
-			auto hView = Kokkos::create_mirror_view(view_);
-			copy(hView, view_);
-			return hView;
-		}
-
-		INLINE_FUNCTION_H 
-		const auto hostVector()const
+		auto hostVector()const
 		{
 			auto hView = Kokkos::create_mirror_view(deviceVector());
 			copy(hView, deviceVector());
 			return hView;
 		}
 
-		INLINE_FUNCTION_H 
-		auto hostVector()
-		{
-			auto hView = Kokkos::create_mirror_view(deviceVector());
-			copy(hView, deviceVector());
-			return hView;
-		}
-
-		// - name of vector 
+		/// Name of the vector 
 		INLINE_FUNCTION_H 
 		const word name()const
 		{
 			return view_.label();
 		}
 		
-		// - size of vector 
+		/// Size of the vector 
 		INLINE_FUNCTION_H 
-		size_t size()const
+		uint32 size()const
 		{
 			return size_;
 		}
 
-		// - capcity of vector 
+		// Capcity of the vector 
 		INLINE_FUNCTION_H 
-		size_t capacity()const
+		uint32 capacity()const
 		{
-			return capacity_;
+			return view_.extent(0);
 		}
 
-		// - if vector is empty 
+		/// If vector is empty 
 		INLINE_FUNCTION_H 
 		bool empty()const
 		{
 			return size_==0;
 		}
 
-		// - reserve capacity for vector 
-		//   preserve the content
+		/// Reserve capacity for vector 
+		/// Preserve the content.
 		INLINE_FUNCTION_H 
-		void reserve(size_t cap)
+		void reserve(uint32 cap) 
 		{
-			changeSize(cap, true);
+			changeCapacity(cap);
 		}
 
-		// - reallocate memory
-		INLINE_FUNCTION_H void reallocate(size_t cap)
-		{
-			capacity_ = cap;
-			size_ = 0;
-			reallocNoInit(view_, capacity_);
-			subViewUpdated_ = false;
-		}
-
-		INLINE_FUNCTION_H void reallocate(size_t cap, size_t size)
-		{
-			capacity_ = cap;
-			size_ = size;
-			reallocNoInit(view_, capacity_);
-			subViewUpdated_ = false;
-		}
-
-		// resize the vector  
+		/// Reallocate memory to new cap and set size to 0. 
 		INLINE_FUNCTION_H 
-		void resize(size_t n){
+		void reallocate(uint32 cap) 
+		{
+			reallocateCapacitySize(cap, 0);			
+		}
+
+		/// Reallocate memory to new cap and set size to newSize. 
+		/// Do not preserve the content
+		INLINE_FUNCTION_H 
+		void reallocate(uint32 cap, uint32 newSize) 
+		{
+			reallocateCapacitySize(cap, newSize);			
+		}
+
+		/// Resize the vector and preserve the content
+		INLINE_FUNCTION_H   
+		void resize(uint32 n){
 			changeSize(n);
 		}
 
-		// resize the view and assign value to the most recent view (most updated)
-		INLINE_FUNCTION_H 
-		void resize(size_t n, const T& val) {
+		/// Resize the vector and assign the value to it. 
+		INLINE_FUNCTION_H  
+		void resize(uint32 n, const T& val) {
 			assign(n, val);
 		}
 
-		// - clear the vector 
+		/// Clear the vector, but keep the allocated memory unchanged 
 		INLINE_FUNCTION_H 
-		void clear() {
+		void clear() 
+		{
 			 size_ = 0;
-			 subViewUpdated_ = false;
-
 		}
 
-		// - fill the range [0,size) with val
+		/// Fill the range [0,size) with val
 		INLINE_FUNCTION_H 
 		void fill(const T& val)
 		{
 			if(empty())return;
-			pFlow::fill(deviceVectorAll(),0 ,size_ ,val);
+			pFlow::fill(view_, rangeU32(0 ,size_) ,val);
 		}
 		
-		// - host calls only
-		// - assign n first elements to val
-		//   resize view
-		//   assign value to either side (device/host)
+		/// Change size of the vector and assign val to vector and 
 		INLINE_FUNCTION_H
 		void assign(size_t n, const T& val)
 		{
-			if(capacity()<n)
+			if( n > capacity() )
 			{
-				this->reallocate(evalCapacity(n));
+				reallocateCapacitySize(evalCapacity(n), n);
+			}
+			else
+			{
+				setSize(n);	
 			}
-			size_ = n;
 			this->fill(val);
 		}
 
-		// - host calls only
-		// - assign source vector
-		//   resize views
-		//   assign to both sides (device&host)
-		INLINE_FUNCTION_H void assign(const Vector<T>& src)
+		
+		/// Assign source vector with specified capacity.
+		/// The size of *this becomes the size of src. 
+		INLINE_FUNCTION_H 
+		void assign(const std::vector<T>& src, uint32 cap) 
 		{
-			auto srcSize = src.size();
-			if( capacity() < srcSize ) 
+			uint32 srcSize = src.size();
+			
+			if(cap != capacity())
 			{
-				this->reallocate( src.capacity() );
+				reallocateCapacitySize(cap, srcSize);
+			}
+			else
+			{
+				setSize(srcSize);
 			}
-			size_ = srcSize;
 			
 			// - unmanaged view in the host
 			hostViewType1D<const T> temp(src.data(), srcSize );
 			copy(deviceVector(), temp);
 		}
 
-
-		//TODO: change it to parallel version 
-		// - delete elements from vector 
-		//   similar memory spaces 
-		/*template<class indT, class MSpace>
+		/// Assign source vector.
+		/// The size of *this becomes the size of src. 
+		/// The capacity of *this becomes the capacity of src.
 		INLINE_FUNCTION_H 
-			typename std::enable_if<
-				Kokkos::SpaceAccessibility<
-					execution_space, typename VectorSingle<indT,MSpace>::memory_space>::accessible,
-					bool>::type
-		deleteElement
-		(
-			const VectorSingle<indT,MSpace>& sortedIndices
-		)
+		void assign(const std::vector<T>& src) 
 		{
-			
-			auto& indices = sortedIndices.deviceVectorAll();
-			auto& dVec = deviceVectorAll();
-			indT  numInd = sortedIndices.size();
-			indT  oldSize = this->size();
-
-			if( numInd == 0 )return true;
-			
-			// an scalalr
-			
-			Kokkos::parallel_for(1, LAMBDA_HD(int nn)
-			 	{
-			 		(void)nn;	
-			 		indT n = 0;
-					indT nextInd = indices[0];
-					indT j = indices[0];
-					for(label i=indices[0]; i < oldSize; ++i)
-					{
-						if( n < numInd && i == nextInd )
-						{
-							++n;
-							nextInd = indices[n];
-						}
-						else
-						{
-							dVec[j] = dVec[i];
-							++j;	
-						}				
-					}
-
-			 	});
-			typename viewType::execution_space().fence();
-			size_ = oldSize - indices.size();
-			subViewUpdated_ = false;
-
-			return true;
+			assign(src, src.capacity());
 		}
 
-		// different memory spaces 
-		template<class indT, class MSpace>
-		INLINE_FUNCTION_H 
-		typename std::enable_if<
-				! Kokkos::SpaceAccessibility<
-					execution_space, typename VectorSingle<indT,MSpace>::memory_space>::accessible,
-					bool>::type
-		deleteElement
-		(
-			const VectorSingle<indT,MSpace>& sortedIndices
-		)
-		{
-			
-			notImplementedFunction;
-		}*/
 
 		INLINE_FUNCTION_H
-		bool insertSetElement(const int32IndexContainer& indices, const T& val)
-		{
-			if(indices.empty()) return true;
-			auto maxInd = indices.max();
-
-			if(this->empty() || maxInd > size()-1 )
-			{
-				resize(maxInd+1);
-			}
-
-			
-			using policy = Kokkos::RangePolicy<
-			execution_space,
-			Kokkos::IndexType<int32> >;
-			auto dVec = deviceVectorAll();
-			auto dIndex = indices.deviceView();
-			
-			Kokkos::parallel_for(
-				"insertSetElement",
-				policy(0,indices.size()), LAMBDA_HD(int32 i){	
-				dVec(dIndex(i))= val;
-				});
-			Kokkos::fence();
-
-			return true;
-		}
-
-		INLINE_FUNCTION_H
-		void sortItems(const int32IndexContainer& indices)
-		{
-			if(indices.size() == 0)
-			{
-				setSize(0);
-				return;
-			}
-			
-			size_t 		newSize = indices.size();
-			viewType 	sortedView("sortedView", newSize);
-
-			using policy = Kokkos::RangePolicy<
-			execution_space,
-			Kokkos::IndexType<int32> >;
-
-			auto d_indices = indices.deviceView();
-			auto d_view = view_;
-			Kokkos::parallel_for(
-				"sortItems",
-				newSize, 
-				LAMBDA_HD(int32 i){	
-					sortedView(i) = d_view(d_indices(i));
-				});
-
-			Kokkos::fence();
-
-			setSize(newSize);
-
-			copy(deviceVector(), sortedView);
-			
-			return;
-
-		}
-
-		INLINE_FUNCTION_H
-		bool insertSetElement(const int32IndexContainer& indices, const Vector<T>& vals)
-		{
-
-			if(indices.size() == 0)return true;
-			if(indices.size() != vals.size())return false;
-
-			auto maxInd = indices.max(); 
-			
-			if(this->empty() || maxInd > size()-1 )
-			{
-				resize(maxInd+1);
-			} 
-
-				
-			hostViewType1D<const T> hVecVals( vals.data(), vals.size());
-			deviceViewType1D<T> dVecVals("dVecVals", indices.size());
-
-			copy(dVecVals, hVecVals);
-
-			using policy = Kokkos::RangePolicy<
-			execution_space,
-			Kokkos::IndexType<int32> >;
-			auto dVec = deviceVectorAll();
-			auto dIndex = indices.deviceView();
+		bool insertSetElement(uint32IndexContainer indices, const T& val);
 		
-			Kokkos::parallel_for(
-				"insertSetElement",
-				policy(0,indices.size()), LAMBDA_HD(int32 i){	
-				dVec(dIndex(i))= dVecVals(i);
-				});
-			Kokkos::fence();
-
-			return true;
-			
-		}
-
 		INLINE_FUNCTION_H
-		bool insertSetElement(const Vector<int32>& indices, const T& val)
-		{
-			if(indices.empty()) return true;
-
-			auto maxInd = max(indices);
-			
-			if(this->empty() || maxInd > size()-1 )
-			{
-				resize(maxInd+1);
-			}
-
-			if constexpr (isHostAccessible_)
-			{
-				hostViewType1D<int32> hostView(const_cast<int32*>(indices.data()), indices.size());
-				fillSelected(deviceVectorAll(), hostView, indices.size(), val);
-				return true;
-			
-			}else
-			{
-
-				// TODO: remove the const_cast 
-				hostViewType1D<int32> hostView(const_cast<int32*>(indices.data()), indices.size());
-				deviceViewType1D<int32> dView("dView", indices.size());
-				copy(dView, hostView);
-				fillSelected(deviceVectorAll(), dView, indices.size(), val);
-				return true;
-			}
-
-			return false;
-		}
-
+		bool insertSetElement(uint32IndexContainer indices, const std::vector<T>& vals);
+		
 		INLINE_FUNCTION_H
-		bool insertSetElement(const Vector<int32>& indices, const Vector<T>& vals)
-		{
-			if(indices.size() == 0)return true;
-			if(indices.size() != vals.size())return false;
+		bool reorderItems(uint32IndexContainer indices);
 
-			auto maxInd = max(indices); 
-			
-			if(this->empty() || maxInd > size()-1 )
-			{
-				resize(maxInd+1);
-			} 
 
-			if constexpr (isHostAccessible_)
-			{
-				// TODO: remove const_cast
-				hostViewType1D<int32> hVecInd( const_cast<int32*>(indices.data()), indices.size());
-				hostViewType1D<T> hVecVals( const_cast<T*>(vals.data()), vals.size());
-
-				fillSelected(deviceVectorAll(), hVecInd, hVecVals, indices.size());
-				return true;
-			
-			}else
-			{
-				
-				// TODO: remove const_cast
-				hostViewType1D<int32> hVecInd( const_cast<int32*>(indices.data()), indices.size());
-				deviceViewType1D<int32> dVecInd("dVecInd", indices.size());
-
-				hostViewType1D<T> hVecVals( const_cast<T*>(vals.data()), vals.size());
-				deviceViewType1D<T> dVecVals("dVecVals", indices.size());
-				
-				copy(dVecVals, hVecVals);
-				copy(dVecInd, hVecInd);
-
-				fillSelected(deviceVectorAll(), dVecInd, dVecVals, indices.size());
-				return true;
-			}
-
-			return false;
-		}
-
-		INLINE_FUNCTION_H
+		/*INLINE_FUNCTION_H
 		bool append(const deviceViewType1D<T>& dVec, size_t numElems)
 		{
 
@@ -734,106 +500,87 @@ public:
 		bool append(const VectorSingle& Vec)
 		{
 			return append(Vec.deviceVector(), Vec.size());
-		}
+		}*/
 
 		// - host calls only
 		//   push a new element at the end
 		//   resize if necessary
 		//   works on host accessible vector 
 		template<bool Enable = true>
-		typename std::enable_if<
-			isHostAccessible_ && Enable,
-			void>::type
+		INLINE_FUNCTION_H
+		typename std::enable_if<isHostAccessible_ && Enable, void>::type
 		push_back(const T& val)
 		{
-			if(size_ == capacity_) 
-			{
-				changeSize(evalCapacity(capacity_), true);
-			}
-			data()[size_++] = val;
-			subViewUpdated_ = false;
+			auto n = changeSize(size_+1);
+			data()[n] = val;	
 		}
 
-		INLINE_FUNCTION_H pointer data(){
+		INLINE_FUNCTION_H 
+		pointer data(){
 			return view_.data();
 		}
 
-		INLINE_FUNCTION_H constPointer data()const{
+		INLINE_FUNCTION_H 
+		constPointer data()const{
 			return view_.data();
 		}
 
-		// - host calls only
-		//   works on host accessible vector 
-		// 	 returns begin iterator 
+		/// Return begin iterator. it works when host is accessible.
 		template<bool Enable = true>
 		INLINE_FUNCTION_H 
-		typename std::enable_if_t<
-			isHostAccessible_ && Enable,
-			iterator>
+		typename std::enable_if_t<isHostAccessible_ && Enable,iterator>
 		begin(){
 			return data();
 		}
 
-		// - host calls only
-		//   works on host accessible vector 
-		//   returns begin iterator 
+		/// Return begin iterator. it works when host is accessible.
 		template<bool Enable = true>
 		INLINE_FUNCTION_H 
-		typename std::enable_if<
-			isHostAccessible_ && Enable,
-			constIterator>::type
+		typename std::enable_if_t<isHostAccessible_ && Enable,constIterator>
 		begin()const {
 			return data();
 		}
 
-		// - host calls only
-		//   works on host accessible vector 
-		//   returns end iterator
+		
+		/// Return end iterator. it works when host is accessible.
 		template<bool Enable = true> 
 		INLINE_FUNCTION_H
-		typename std::enable_if<
-			isHostAccessible_ && Enable,
-			iterator>::type
+		typename std::enable_if_t<isHostAccessible_ && Enable,iterator>
 		end(){
 			return size_ > 0 ? data() + size_: data();
 		}
 
-		// host call
-		// returns end iterator 
+		
+		/// Return end iterator. it works when host is accessible.
 		template<bool Enable = true>
 		INLINE_FUNCTION_H
-		typename std::enable_if<
-			isHostAccessible_ && Enable,
-			constIterator>::type
+		typename std::enable_if_t<isHostAccessible_ && Enable,constIterator>
 		end()const{
 			return size_ > 0 ? data() + size_: data();
 		}
 				
-		// operator to be used on host side vectors
+		/// Return reference to element i. it works when host is accessible.
 		template<bool Enable = true>
 		INLINE_FUNCTION_H 
-		typename std::enable_if<
-			isHostAccessible_ && Enable,
-			reference>::type
-		operator[](label i){
+		typename std::enable_if_t<isHostAccessible_ && Enable,reference>
+		operator[](size_t i){
 			return view_[i];
 		}
 
+		/// Return reference to element i. it works when host is accessible.
 		template<bool Enable = true>
 		INLINE_FUNCTION_H 
-		typename std::enable_if<
-			isHostAccessible_ && Enable,
-			constReference>::type
-		 operator[](label i)const{
+		typename std::enable_if_t<isHostAccessible_ && Enable,constReference>
+		 operator[](size_t i)const{
 			return view_[i];
 		}
 
 	//// - IO operations
 
+		/// Read vector from is stream (ASCII or Binary)
+		/// For binary read, len should be provided 
 		FUNCTION_H
-		bool readVector(
-			iIstream& is,
-			size_t len=0)
+		bool readVector(iIstream& is,	size_t len=0)
 		{
 			Vector<T> vecFromFile;
 			if( !vecFromFile.readVector(is,len) ) return false;
@@ -843,23 +590,21 @@ public:
 			return true;
 		}
 
+		/// Read vector from stream (ASCII)
 		FUNCTION_H
 		bool read(iIstream& is)
 		{
 			return readVector(is);
 		}
 
+		/// Write the vector to os (ASCII or Binary)
 		FUNCTION_H
 		bool write(iOstream& os)const
 		{
-			
-			Vector<T, noConstructAllocator<T>> vecToFile(this->size());
-			
-			const auto dVec = Kokkos::subview(view_, Kokkos::make_pair(0,int(size_)));
-			hostViewType1D<T> mirror(vecToFile.data(), vecToFile.size());
-			copy(mirror,dVec);
-
-			return vecToFile.write(os);
+			auto hVec = hostVector();
+			auto sp = span<T>( const_cast<T*>(hVec.data()), hVec.size());
+			os<<sp;
+			return true;
 		}
 
 }; // class VectorSingle
@@ -894,6 +639,7 @@ inline iOstream& operator << (iOstream& os, const VectorSingle<T, MemorySpace>&
 } // - pFlow
 
 #include "VectorSingleAlgorithms.hpp"
+#include "VectorSingleI.hpp"
 
 
 #endif //__VectorSingle_hpp__
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingleI.hpp b/src/phasicFlow/containers/VectorHD/VectorSingleI.hpp
new file mode 100644
index 00000000..61217a96
--- /dev/null
+++ b/src/phasicFlow/containers/VectorHD/VectorSingleI.hpp
@@ -0,0 +1,197 @@
+/*------------------------------- 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.
+
+-----------------------------------------------------------------------------*/
+
+
+template<typename T, typename MemorySpace>
+INLINE_FUNCTION_H
+bool pFlow::VectorSingle<T,MemorySpace>::insertSetElement(uint32IndexContainer indices, const T& val)
+{
+	if(indices.empty()) return true;
+	
+	auto maxInd = indices.max();
+
+	if(this->empty() || maxInd > size()-1 )
+	{
+		resize(maxInd+1);
+	}
+
+	using policy = Kokkos::RangePolicy<executionSpace,Kokkos::IndexType<uint32>>;
+
+	if constexpr( isDeviceAccessible_ )
+	{
+		auto v = view_;
+		auto ind = indices.deviceView();
+
+		Kokkos::parallel_for(
+			"VectorSingle::insertSetElement",
+			policy(0,indices.size()),
+			LAMBDA_HD(uint32 i){
+				v[ind[i]]= val;
+			});
+
+		Kokkos::fence();
+	}
+	else
+	{
+
+		auto v = view_;
+		auto ind = indices.hostView();
+		
+		Kokkos::parallel_for(
+			"VectorSingle::insertSetElement",
+			policy(0,indices.size()),
+			LAMBDA_HD(uint32 i){
+				v[ind[i]]= val;
+			});
+
+		Kokkos::fence();
+
+	}
+
+	return true;
+}
+
+
+template<typename T, typename MemorySpace>
+INLINE_FUNCTION_H
+bool pFlow::VectorSingle<T,MemorySpace>::insertSetElement
+(
+	const uint32IndexContainer indices, 
+	const std::vector<T>& vals
+)
+{
+
+	if(indices.size() == 0)return true;
+	if(indices.size() != vals.size())return false;
+
+	auto maxInd = indices.max(); 
+	
+	if(this->empty() || maxInd > size()-1 )
+	{
+		resize(maxInd+1);
+	} 
+	
+	using policy = Kokkos::RangePolicy<executionSpace,Kokkos::IndexType<uint32>>;
+
+	hostViewType1D<const T> hVals( vals.data(), vals.size());
+
+	if constexpr( isDeviceAccessible_ )
+	{
+		deviceViewType1D<T> dVals("dVals", indices.size());
+		copy(dVals, hVals);
+		auto dVec = view_;
+		auto ind = indices.deviceView();
+
+		Kokkos::parallel_for(
+			"VectorSingle::insertSetElement",
+			policy(0,indices.size()), LAMBDA_HD(int32 i){	
+				dVec(ind(i))= dVals(i);}
+			);
+
+		Kokkos::fence();
+
+	}
+	else
+	{
+		auto dVec = view_;
+		auto ind = indices.hostView();
+
+		Kokkos::parallel_for(
+			"VectorSingle::insertSetElement",
+			policy(0,indices.size()), LAMBDA_HD(int32 i){	
+				dVec(ind(i))= hVals(i);}
+			);
+
+		Kokkos::fence();
+
+	}
+		
+	return true;
+}
+
+
+template<typename T, typename MemorySpace>
+INLINE_FUNCTION_H
+bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(uint32IndexContainer indices)
+{
+	if(indices.size() == 0)
+	{
+		setSize(0);
+		return true;
+	}
+	
+	auto maxInd = indices.max();
+
+	if(maxInd >= this->size())
+	{
+		fatalErrorInFunction<<"In reordering the VectorSingle ("
+		<< this->name()<< ") maximum index ("<< maxInd << 
+		") exceeds the size of the vector ("	<< this->size()<<")"<<endl;
+		return false;
+	}
+
+	uint32 		newSize = indices.size();
+	
+	setSize(newSize);
+
+	viewType 	sortedView(this->name(), newSize);
+
+	using policy = Kokkos::RangePolicy< executionSpace,Kokkos::IndexType<uint32>>;
+
+	if constexpr( isDeviceAccessible_)
+	{
+		auto d_indices = indices.deviceView();
+		auto d_view = view_;
+		
+		Kokkos::parallel_for
+		(
+			"VectorSingle::sortItems",
+			policy(0,newSize), 
+			LAMBDA_HD(uint32 i)
+			{	
+				sortedView(i) = d_view(d_indices(i));
+			}
+		);
+
+		Kokkos::fence();
+
+	}
+	else
+	{
+		auto h_indices = indices.hostView();
+		auto d_view = view_;
+		
+		Kokkos::parallel_for
+		(
+			"VectorSingle::sortItems",
+			policy(0,newSize), 
+			LAMBDA_HD(uint32 i)
+			{	
+				sortedView(i) = d_view(h_indices(i));
+			}
+		);
+		
+		Kokkos::fence();
+	}
+
+	copy(deviceVector(), sortedView);
+	
+	return true;
+}
\ No newline at end of file
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingles.hpp b/src/phasicFlow/containers/VectorHD/VectorSingles.hpp
index ecae1398..d0c3acdf 100644
--- a/src/phasicFlow/containers/VectorHD/VectorSingles.hpp
+++ b/src/phasicFlow/containers/VectorHD/VectorSingles.hpp
@@ -33,10 +33,6 @@ typedef VectorSingle<int8> 				int8Vector_D;
 
 typedef VectorSingle<int8, HostSpace> 	int8Vector_H;
 
-typedef VectorSingle<int16> 			int16Vector_D;
-
-typedef VectorSingle<int16, HostSpace> 	int16Vector_H;
-
 typedef VectorSingle<int32> 			int32Vector_D;
 
 typedef VectorSingle<int32, HostSpace> 	int32Vector_H;
@@ -45,13 +41,17 @@ typedef VectorSingle<int64> 			int64Vector_D;
 
 typedef VectorSingle<int64, HostSpace> 	int64Vector_H;
 
+typedef VectorSingle<uint8> 			uint8Vector_D;
+
+typedef VectorSingle<uint8, HostSpace> 	uint8Vector_H;
+
 typedef VectorSingle<uint32> 			uint32Vector_D;
 
 typedef VectorSingle<uint32, HostSpace> uint32Vector_H;
 
-typedef VectorSingle<label> 			labelVector_D;
+typedef VectorSingle<uint64> 			uint64Vector_D;
 
-typedef VectorSingle<label, HostSpace> 	labelVector_H;
+typedef VectorSingle<uint64, HostSpace> uint64Vector_H;
 
 typedef VectorSingle<real> 				realVector_D;
 
@@ -61,14 +61,6 @@ typedef VectorSingle<realx3> 				realx3Vector_D;
 
 typedef VectorSingle<realx3, HostSpace> 	realx3Vector_H;
 
-typedef VectorSingle<uint16x3> 				uint16x3Vector_D;
-
-typedef VectorSingle<uint16x3, HostSpace> 	uint16x3Vector_H;
-
-typedef VectorSingle<uint32x3> 				uint32x3Vector_D;
-
-typedef VectorSingle<uint32x3, HostSpace> 	uint32x3Vector_H;
-
 typedef VectorSingle<int32x3> 				int32x3Vector_D;
 
 typedef VectorSingle<int32x3, HostSpace> 	int32x3Vector_H;
@@ -77,6 +69,10 @@ typedef VectorSingle<int64x3> 				int64x3Vector_D;
 
 typedef VectorSingle<int64x3, HostSpace> 	int64x3Vector_H;
 
+typedef VectorSingle<uint32x3> 				uint32x3Vector_D;
+
+typedef VectorSingle<uint32x3, HostSpace> 	uint32x3Vector_H;
+
 typedef VectorSingle<realx3x3> 				realx3x3Vector_D;
 
 typedef VectorSingle<realx3x3, HostSpace> 	realx3x3Vector_H;
diff --git a/src/phasicFlow/containers/indexContainer/indexContainer.hpp b/src/phasicFlow/containers/indexContainer/indexContainer.hpp
index 7b517a5a..10fce100 100644
--- a/src/phasicFlow/containers/indexContainer/indexContainer.hpp
+++ b/src/phasicFlow/containers/indexContainer/indexContainer.hpp
@@ -21,32 +21,41 @@ Licence:
 #ifndef __indexContainer_hpp__ 
 #define __indexContainer_hpp__
 
-#include "span.hpp"
-#include "KokkosTypes.hpp"
-#include "KokkosUtilities.hpp"
-#include "ViewAlgorithms.hpp"
+#include <vector>
+
+#include "phasicFlowKokkos.hpp"
 
 
 namespace pFlow
 {
 
+/**
+ * It holds two vectors of indecis on Host and Device.
+ * Host vector should be used for threads running on
+ * Host and Device vector should be used for threads
+ * running on Device.
+ */
 template<typename IndexType>
 class indexContainer
 {
 public:
 	
-	using DualViewType 		= Kokkos::DualView<IndexType*>;
+	using DVType 					= DualViewType1D<IndexType>;
+	///  Device type on device 
+	using DeviceViewType 	= typename DVType::t_dev;
 
-	// - viewType of data on device 
-  	using DeviceViewType 	= typename DualViewType::t_dev;
+	/// Host type on device 
+	using HostViewType 		= typename DVType::t_host;
 
-  	// - viewType of data on host
-  	using HostViewType 		= typename DualViewType::t_host;
+	/// Host memory type
+	using HostType 			= typename HostViewType::device_type;
 
-  	using HostType 			= typename HostViewType::device_type;
-
-  	using DeviceType 		= typename DeviceViewType::device_type;
+	/// Device memory ype
+	using DeviceType 		= typename DeviceViewType::device_type;
 
+	/**
+	 * Helper class for accessing index on host or device 
+	 */
   	template<typename ViewType>
   	class IndexAccessor
   	{
@@ -58,169 +67,207 @@ public:
   		view_(v){}
 
   		INLINE_FUNCTION_HD
-  		IndexType operator()(int32 i)const
+  		IndexType operator()(uint32 i)const
   		{
   			return view_[i];
   		}
+
+  		uint32 size()const
+  		{
+  			return view_.extent(0);
+  		}
   	};
 
 protected:
 
-	int32 		min_  	= 0;
-	int32 		max_ 	= 0;
-	size_t 		size_	= 0; 	
+	/// min value in indices
+	IndexType 		min_  	= 0;
 
-	DualViewType 	views_;
+	/// max value in the indices
+	IndexType 		max_ 	= 0;
+
+	/// number/size of index vector 
+	uint32 			size_	= 0; 	
+
+	/// views to hold indices on Host and Device
+	DVType 			views_;
 
 public:
 
-	indexContainer(){}
+	//// - Constructors 
 
-	// half open [begin,end)
-	indexContainer(IndexType begin, IndexType end)
-	:
-		min_(begin),
-		max_(end-1),
-		size_(end-begin),
-		views_("indexContainer", size_)
-	{
-		pFlow::fillSequence(views_.h_view, 0, size_, min_);
-		copy(views_.d_view, views_.h_view);
-	}
+		/// Default
+		indexContainer(){}
 
+		/// Construct from a Range (half open)
+		template<typename T>
+		indexContainer(const Range<T>& rng)
+		:
+			indexContainer
+			(
+				static_cast<IndexType>(rng.begin()), 
+				static_cast<IndexType>(rng.end())
+			)
+		{}
 
-	indexContainer(IndexType* data, int32 numElems)
-	:
-		size_(numElems),
-		views_("indexContainer", numElems)
-	{
-		HostViewType hData(data, numElems);
-		copy(views_.h_view, hData);
-		copy(views_.d_view, views_.h_view);
-		min_ 	= pFlow::min(views_.d_view, 0, numElems);
-		max_ 	= pFlow::max(views_.d_view, 0, numElems);
-	}
-
-	indexContainer(const indexContainer&) = default;
-
-	indexContainer& operator = (const indexContainer&) = default;
-
-	indexContainer(indexContainer&&) = default;
-
-	indexContainer& operator = (indexContainer&&) = default;	
-
-	~indexContainer() = default;
-
-	INLINE_FUNCTION_HD
-	size_t size()const
-	{
-		return size_;
-	}
-
-	INLINE_FUNCTION_HD
-	size_t empty()const
-	{
-		return size_==0;
-	}
-	
-	INLINE_FUNCTION_HD
-	IndexType min()const
-	{
-		return min_;
-	}
-
-	INLINE_FUNCTION_HD
-	IndexType max()const
-	{
-		return max_;
-	}
-
-	template<typename executionSpace>
-	INLINE_FUNCTION_HD
-	IndexType operator()(selectSide<executionSpace>,int32 i)const
-	{	
-		if constexpr (isHostAccessible<executionSpace>())
+		/// Construct half open [begin,end)
+		indexContainer(IndexType begin, IndexType end)
+		:
+			min_(begin),
+			max_(end-1),
+			size_(end-begin),
+			views_("indexContainer", size_)
 		{
-			return views_.h_view(i);
-		}else
-		{
-			return views_.d_view(i);
-		}
-	}
-
-	const HostViewType& hostView()const
-	{
-		return views_.h_view;
-	}
-
-	const DeviceViewType& deviceView()const
-	{
-		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);
-	}
-
-	auto indicesDevice()const
-	{
-		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;
+			pFlow::fillSequence(views_.h_view, 0, size_, min_);
+			copy(views_.d_view, views_.h_view);
 		}
 
-		if(findMinMax)
+		/// From data and number of elements in data. 
+		/// data is a pointer in the Host memory 
+		indexContainer(IndexType* data, int32 numElems)
+		:
+			size_(numElems),
+			views_("indexContainer", numElems)
 		{
-			min_ 	= pFlow::min(views_.d_view, 0, size_);
-			max_ 	= pFlow::max(views_.d_view, 0, size_);
+			HostViewType hData(data, numElems);
+			copy(views_.h_view, hData);
+			copy(views_.d_view, views_.h_view);
+			min_ 	= pFlow::min(views_.d_view, 0, numElems);
+			max_ 	= pFlow::max(views_.d_view, 0, numElems);
+		}
+
+		indexContainer(std::vector<IndexType> &ind)
+		:
+			indexContainer(ind.data(), ind.size())
+		{}
+
+		/// Copy
+		indexContainer(const indexContainer&) = default;
+
+		/// Copy assignment
+		indexContainer& operator = (const indexContainer&) = default;
+
+		/// Move
+		indexContainer(indexContainer&&) = default;
+
+		/// Mover assignement 
+		indexContainer& operator = (indexContainer&&) = default;	
+
+		/// Destructor
+		~indexContainer() = default;
+
+	//// - Methods 
+
+		/// Size 
+		INLINE_FUNCTION_HD
+		auto size()const
+		{
+			return size_;
+		}
+
+		/// If the container empty
+		INLINE_FUNCTION_HD
+		bool empty()const
+		{
+			return size_==0;
+		}
+		
+		/// Min value of indices 
+		INLINE_FUNCTION_HD
+		IndexType min()const
+		{
+			return min_;
+		}
+
+		/// Max value of indices 
+		INLINE_FUNCTION_HD
+		IndexType max()const
+		{
+			return max_;
+		}
+
+		/// Return Host veiw
+		const HostViewType& hostView()const
+		{
+			return views_.h_view;
+		}
+
+		/// Return Device view
+		const DeviceViewType& deviceView()const
+		{
+			return views_.d_view;
+		}
+
+		/// Return Host veiw
+		HostViewType& hostView()
+		{
+			return views_.h_view;
+		}
+
+		/// Return Device veiw
+		DeviceViewType& deviceView()
+		{
+			return views_.d_view;
+		}
+
+		/// Return index accessor that works on Host
+		auto indicesHost()const
+		{
+			return IndexAccessor<HostViewType>(views_.h_view);
+		}
+
+		/// Return index accessor that works on Device 
+		auto indicesDevice()const
+		{
+			return IndexAccessor<DeviceViewType>(views_.d_view);
+		}
+
+		/// Mark host is modified 
+		void modifyOnHost()
+		{
+			views_.modify_host();
+		}
+
+		/// Mark device is modified
+		void modifyOnDevice()
+		{
+			views_.modify_device();
+		}
+
+		/// synchronize views 
+		void syncViews()
+		{
+			bool findMinMax = false;
+			if(views_.template need_sync<HostType>())
+			{
+				Kokkos::deep_copy(views_.d_view, views_.h_view);
+				views_.clear_sync_state();
+				findMinMax = true;
+			}
+			else if(views_.template need_sync<DeviceType>())
+			{
+				Kokkos::deep_copy(views_.h_view, views_.d_view);
+				views_.clear_sync_state();
+				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;
-	}
 };
 
 
+
 using int32IndexContainer = indexContainer<int32>;
 using int64IndexContainer = indexContainer<int64>;
 
+using uint32IndexContainer = indexContainer<uint32>;
+using uint64IndexContainer = indexContainer<uint64>;
+
 
 }
 
diff --git a/src/phasicFlow/containers/span/span.hpp b/src/phasicFlow/containers/span/span.hpp
index 63c89190..0ee41e71 100644
--- a/src/phasicFlow/containers/span/span.hpp
+++ b/src/phasicFlow/containers/span/span.hpp
@@ -51,7 +51,7 @@ protected:
 
     T* data_ 		= nullptr;
 
-    label 	size_ 	= 0;
+    uint32 	size_ 	= 0;
 
 public:
     
@@ -63,7 +63,7 @@ public:
 
 
     INLINE_FUNCTION_HD
-    span(T* data, label size)
+    span(T* data, uint32 size)
         : data_(data), size_(size)
     {}
 
@@ -98,7 +98,7 @@ public:
 
     /// Returns the number of elements in the span
     INLINE_FUNCTION_HD
-    label size() const
+    uint32 size() const
     {
         return size_;
     }
@@ -131,6 +131,18 @@ public:
         return data_ + size_;
     }
 
+    INLINE_FUNCTION_HD
+    T& operator[](uint32 i)
+    {
+    	return data_[i];
+    }
+
+    INLINE_FUNCTION_HD
+    const T& operator[](uint32 i)const
+    {
+    	return data_[i];
+    }
+
     INLINE_FUNCTION_HD
     T& operator[](int32 i)
     {
@@ -143,16 +155,35 @@ public:
     	return data_[i];
     }
 
-    INLINE_FUNCTION_HD
-    T& operator[](label i)
+    /// Write in ASCII format, can be sued for all variable types
+    INLINE_FUNCTION_H
+    bool writeASCII(iOstream& os) const
     {
-    	return data_[i];
+        os << token::BEGIN_LIST;
+        if(size()>0)
+        {
+            for(uint32 i=0; i<size()-1; i++)
+            {
+                os << data_[i]<<token::NL;
+            }
+            os << data_[size()-1] << token::END_LIST;
+        }
+        else
+        {
+            os<< token::END_LIST;
+        }
+        
+        os.check(FUNCTION_NAME);
+        return true;
     }
 
-    INLINE_FUNCTION_HD
-    const T& operator[](label i)const
+    /// Write in Binray format, can be used for numeral types (Not word or similar ones)
+    INLINE_FUNCTION_H
+    bool writeBinary(iOstream& os) const
     {
-    	return data_[i];
+        os.write(reinterpret_cast<const char*>(data_), this->size()*sizeof(T));
+        os.check(FUNCTION_NAME);
+        return true;
     }
 
 };
@@ -161,15 +192,15 @@ template<typename T>
 inline 
 iOstream& operator<<(iOstream& os, const span<T>& s)
 {
-    os << token::BEGIN_LIST;
-    for(size_t i=0; i<s.size(); i++)
-    {
-        os << s[i]<<token::NL;
-    }
-    
-    os << token::END_LIST;
 
-    os.check(FUNCTION_NAME);
+    if( os.isBinary() && !std::is_same_v<T,word>)
+    {
+        s.writeBinary(os);
+    }
+    else
+    {
+        s.writeASCII(os);
+    }
 
     return os;
 }
diff --git a/src/phasicFlow/repository/IOobject/IOfileHeader.cpp b/src/phasicFlow/repository/IOobject/IOfileHeader.cpp
index 2aae5642..8c92c944 100644
--- a/src/phasicFlow/repository/IOobject/IOfileHeader.cpp
+++ b/src/phasicFlow/repository/IOobject/IOfileHeader.cpp
@@ -12,12 +12,12 @@ Licence:
   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 
+  phasicFlow is distribute+d 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 "processors.hpp"
 #include "IOfileHeader.hpp"
 #include "repository.hpp"
 
@@ -133,9 +133,24 @@ bool pFlow::IOfileHeader::readIfPresent()const
 	return fileExist() && isReadIfPresent();
 }
 
-
-bool pFlow::IOfileHeader::writeHeader(iOstream& os, const word& typeName) const
+bool pFlow::IOfileHeader::writeHeader()const
 {
+	if( !this->readWriteHeader() ) return false;
+	if( !processors::isMaster() ) return false;
+	if( !implyWrite() ) return false;
+
+	return true;
+}
+
+bool pFlow::IOfileHeader::writeHeader
+(
+	iOstream& os, 
+	const word& typeName, 
+	bool forceWrite
+) const
+{
+
+	if(!forceWrite && !writeHeader()) return true;
 
 	writeBanner(os);
 
@@ -158,14 +173,22 @@ bool pFlow::IOfileHeader::writeHeader(iOstream& os, const word& typeName) const
 	return true;
 }
 
-bool pFlow::IOfileHeader::writeHeader(iOstream& os) const
+bool pFlow::IOfileHeader::writeHeader(iOstream& os, bool forceWrite) const
 {
-	return writeHeader(os, objectType_);
+	return writeHeader(os, objectType_, forceWrite);
+}
+
+bool pFlow::IOfileHeader::readHeader()const
+{
+	if( !this->readWriteHeader() ) return false;
+	return true;
 }
 
 bool pFlow::IOfileHeader::readHeader(iIstream& is, bool silent)
 {
 	
+	if(!readHeader()) return true;
+
 	if( !is.findTokenAndNextSilent("objectName", objectName_) )
 	{
 		if(!silent)
@@ -200,6 +223,19 @@ bool pFlow::IOfileHeader::readHeader(iIstream& is, bool silent)
 		return false;
 	}
 
+	return true;
+} 
+
+bool pFlow::IOfileHeader::writeData()const
+{
+	if( processors::isMaster() || this->differentDataOnProcessors())
+		return true;
+	else
+		return false;
+}
+
+bool pFlow::IOfileHeader::readData()const
+{
 	return true;
 }
 
diff --git a/src/phasicFlow/repository/IOobject/IOfileHeader.hpp b/src/phasicFlow/repository/IOobject/IOfileHeader.hpp
index fea5cb66..e019a79e 100644
--- a/src/phasicFlow/repository/IOobject/IOfileHeader.hpp
+++ b/src/phasicFlow/repository/IOobject/IOfileHeader.hpp
@@ -95,26 +95,45 @@ public:
 	//   read the header of the file to check if it is ok
 	bool headerOk(bool silent = false);
 	
-	// - imply read 
+	/// Imply read 
 	bool implyRead() const;
 
-	// - imply write 
+	/// Imply write 
 	bool implyWrite() const;
 	
-	// - check if file exists
+	/// Check if file exists
 	bool fileExist() const;
 	
-	// - check read if present 
+	/// Check read if present 
 	bool readIfPresent()const;
 
-	// - write the header in the  file , typeName comes from caller
-	bool writeHeader(iOstream& os, const word& typeName) const;
+	/// Check if the header should be written to file
+	/// True: on master + implyWrite + readWriteHeader = true
+	/// False: slave or NOT implyWrite
+	bool writeHeader()const;
 
-	// - write the header in the  file, typeName comes from the one read from file 
-	bool writeHeader(iOstream& os) const;
-	
-	// - read the header in the file 
+	/// Write the header to the file , typeName comes from caller
+	bool writeHeader(iOstream& os, const word& typeName, bool forceWrite = false ) const;
+
+	/// Write the header to the file, typeName comes from the one read from file 
+	bool writeHeader(iOstream& os, bool forceWrite = false) const;
+
+	/// Check if the data should be written to file 
+	/// True: on master or differentDataOnProcessor is true
+	/// False: otherwise 
+	bool writeData()const;
+
+	/// Check if header should be read from file
+	/// True: All processors, read the file header 
+	/// False: readWriteHeader = false
+	bool readHeader()const;
+
+	/// Read the header in the file
 	bool readHeader(iIstream& is, bool silent=false);
+
+	/// Check if data should be read from file 
+	/// Always return true
+	bool readData()const;
 	
 	// - write the banner 
 	bool writeBanner(iOstream& os)const;
diff --git a/src/phasicFlow/repository/IOobject/IOobject.hpp b/src/phasicFlow/repository/IOobject/IOobject.hpp
index 7b2ab200..2b947946 100644
--- a/src/phasicFlow/repository/IOobject/IOobject.hpp
+++ b/src/phasicFlow/repository/IOobject/IOobject.hpp
@@ -61,23 +61,25 @@ private:
 	:
 		public iObject
 	{
-	public:
+	protected:
 		dataType 	data_;
 
 	public:
 
 		template<typename... Args,
-          		typename = std::enable_if_t<!std::is_constructible<object_t, Args&&...>::value>>
+          		typename = std::enable_if_t<std::is_constructible<dataType, Args&&...>::value>>
 		object_t(Args&&... args)
 		:
 			data_(std::forward<Args>(args)...)
-		{}
+		{
+			/*constexpr word msg(dataType::TYPENAME()+"input is not a member function.")
+			static_assert(std::is_member_function_pointer<decltype(&dataType::write)>::value,
+                  msg.c_str());*/
+		}
 
 		// cunstruct by copying data
-		object_t(const dataType& data): data_(data){}
-
-		// construct by moving data
-		//object_t(dataType&& data): data_(std::move(data)){}
+		object_t(const dataType& data): data_(data)
+		{}
 
 		
 		virtual uniquePtr<iObject> clone() const
@@ -100,6 +102,16 @@ private:
 			return data_.write(os);
 		}
 
+		auto& data()
+		{
+			return data_;
+		}
+
+		const auto& data()const
+		{
+			return data_;
+		}
+
 	};
 
 protected:
@@ -107,8 +119,9 @@ protected:
 	//// - data members
 
 		// underlaying data object 
-		uniquePtr<iObject> object_;
+		uniquePtr<iObject> object_ = nullptr;
 
+		
 
 public:
 
@@ -127,6 +140,42 @@ public:
 		// - construct from components, transfer the ownership of IOobject to the owner (no read happens)
 		IOobject(const objectFile& objf, const repository* owner, uniquePtr<IOobject>&& obj);
 		
+		
+		template<typename T, 
+			typename = std::enable_if_t<
+				!std::is_same<T, uniquePtr<IOobject::iObject>>::value && 
+				!std::is_same<T, uniquePtr<IOobject>>::value>>
+		IOobject(const objectFile& objf, const repository* owner, const T& data)
+		:
+			IOfileHeader(objf, owner),
+			object_(makeUnique<object_t<T>>(data))
+		{
+			if(!read(this->readWriteHeader()))
+			{
+				fatalErrorInFunction<<
+				"error in reading " << name() << " from path " << path()<<endl;
+				fatalExit;
+			}
+		}
+
+		
+		template<typename T, 
+			typename = std::enable_if_t<
+				!std::is_same<T, uniquePtr<IOobject::iObject>>::value && 
+				!std::is_same<T, uniquePtr<IOobject>>::value>>
+		IOobject(const objectFile& objf, const repository* owner,  T&& data)
+		:
+			IOfileHeader(objf, owner),
+			object_(makeUnique<object_t<T>>(data))
+		{
+			if(!read(this->readWriteHeader()))
+			{
+				fatalErrorInFunction<<
+				"error in reading " << name() << " from path " << path()<<endl;
+				fatalExit;
+			}
+		}		
+
 		// - copy construct 
 		IOobject(const IOobject& src)=delete;	
 
@@ -142,7 +191,7 @@ public:
 		
 		// - construct object_t with the Args as the arguments of object constructor 
 		template<typename T, typename... Args>
-		static auto make_object_t(Args&&... args);
+		static uniquePtr<iObject> make_object_t(Args&&... args);
 		
 
 	//// - Access to data object 
@@ -158,7 +207,6 @@ public:
 		template<typename T>
 		const auto& getObject()const;
 	
-		
 
 	//// - IO operations 
 
diff --git a/src/phasicFlow/repository/IOobject/IOobjectTemplates.cpp b/src/phasicFlow/repository/IOobject/IOobjectTemplates.cpp
index 4c8f20e1..706f5ced 100644
--- a/src/phasicFlow/repository/IOobject/IOobjectTemplates.cpp
+++ b/src/phasicFlow/repository/IOobject/IOobjectTemplates.cpp
@@ -18,6 +18,7 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+
 template<typename T, typename... Args>
 auto pFlow::IOobject::make
 (
@@ -32,7 +33,7 @@ auto pFlow::IOobject::make
 }
 
 template<typename T, typename... Args>
-auto pFlow::IOobject::make_object_t(Args&&... args)
+pFlow::uniquePtr<pFlow::IOobject::iObject> pFlow::IOobject::make_object_t(Args&&... args)
 {
 	auto ptr = makeUnique<object_t<T>>(std::forward<Args>(args)...);
 	return ptr;
@@ -47,7 +48,7 @@ auto& pFlow::IOobject::getObject()
 		"accessing an invalid objecct "<< name() <<endl;
 		fatalExit;
 	}
-	return  dynamic_cast<object_t<T>&>(*object_).data_;
+	return  static_cast<object_t<T>&>(*object_).data();
 }
 
 
@@ -60,5 +61,5 @@ const auto& pFlow::IOobject::getObject()const
 		"accessing an invalid objecct "<< name() <<endl;
 		fatalExit;
 	}
-	return  dynamic_cast<const object_t<T>&>(*object_).data_;	
+	return  static_cast<const object_t<T>&>(*object_).data();	
 }
\ No newline at end of file
diff --git a/src/phasicFlow/repository/IOobject/objectFile.cpp b/src/phasicFlow/repository/IOobject/objectFile.cpp
index 55776975..c9957daf 100644
--- a/src/phasicFlow/repository/IOobject/objectFile.cpp
+++ b/src/phasicFlow/repository/IOobject/objectFile.cpp
@@ -25,10 +25,7 @@ pFlow::objectFile::objectFile
 	const word& name
 )
 :
-	name_(name),
-	rFlag_(READ_NEVER),
-	wFlag_(WRITE_NEVER),
-	localPath_("")
+	name_(name)
 {}
 
 pFlow::objectFile::objectFile
@@ -37,13 +34,15 @@ pFlow::objectFile::objectFile
 	const fileSystem& 		localPath,
 	const readFlag&   		rf, 
 	const writeFlag&  		wf,
-	bool  rwHeader
+	bool  diffDataOnProcessors,
+	bool  rwHdr
 )
 :
 	name_(name),
 	rFlag_(rf),
 	wFlag_(wf),
 	localPath_(localPath),
-	readWriteHeader_(rwHeader)
+	differentDataOnProcessors_(diffDataOnProcessors),
+	readWriteHeader_(rwHdr)	
 {
 }
\ No newline at end of file
diff --git a/src/phasicFlow/repository/IOobject/objectFile.hpp b/src/phasicFlow/repository/IOobject/objectFile.hpp
index 7ab6f153..6dd70ca2 100644
--- a/src/phasicFlow/repository/IOobject/objectFile.hpp
+++ b/src/phasicFlow/repository/IOobject/objectFile.hpp
@@ -49,18 +49,30 @@ public:
 
 protected:
 
-	// name of the entity
+	/// Name of the entity
 	word 	name_;
 
-	// read flag
-	readFlag  rFlag_;
+	/// Read flag
+	readFlag  rFlag_ = READ_NEVER;
 
-	// write flag
-	writeFlag wFlag_;
+	/// Write flag
+	writeFlag wFlag_ = WRITE_NEVER;
 
-	// local path of entity
-	fileSystem localPath_;
+	/// Local path of entity
+	fileSystem localPath_ = "";
 
+	/// Number of bytes used for writing/reading real variable (mostly used for binray)
+	int 		numBytesForReal_ = numBytesForReal__;
+	
+	/// All processors have the different set of data or not?
+	/// True: Each processor should read its own part of data from file and should write 
+	/// its own part of data to the file.  
+	/// Flase: All processors read the same data from file and in writing, only master data
+	/// write the data to the file. 
+
+	bool differentDataOnProcessors_ = true;
+
+	/// Does the objectFile write the header or not
 	bool readWriteHeader_ = true;
 
 public:
@@ -80,7 +92,8 @@ public:
 		const fileSystem& localPath,
 		const readFlag&   rf = READ_NEVER,
 		const writeFlag&  wf = WRITE_NEVER,
-		bool  rwHeader = true
+		bool  diffDataOnProcessors = true,
+		bool  rwHdr = true
 	);
 
 	// copy construct 
@@ -92,6 +105,7 @@ public:
 
 	objectFile& operator = (objectFile && rhs) = default;
 
+
 	virtual ~objectFile()=default;
 
 	virtual word name() const
@@ -104,41 +118,55 @@ public:
 		return localPath_;
 	}
 
+	inline
 	readFlag rFlag()const
 	{
 		return rFlag_;
 	}
 
+	inline
 	writeFlag wFlag()const
 	{
 		return wFlag_;
 	}
 
+	inline
 	bool isReadAlways()const
 	{
 		return rFlag_ == READ_ALWAYS;
 	}
 
+	inline
 	bool isReadNever()const
 	{
 		return rFlag_ == READ_NEVER;
 	}
 
+	inline
 	bool isReadIfPresent()const
 	{
 		return rFlag_ == READ_IF_PRESENT;
 	}
 
+	inline
 	bool isWriteAlways()const
 	{
 		return wFlag_ == WRITE_ALWAYS;
 	}
 
+	inline
 	bool isWriteNever()const
 	{
 		return wFlag_ == WRITE_NEVER;
 	}
 
+	inline
+	bool differentDataOnProcessors()const
+	{
+		return differentDataOnProcessors_;
+	}
+
+	inline 
 	bool readWriteHeader()const
 	{
 		return readWriteHeader_;
diff --git a/src/phasicFlow/streams/Stream/Istream.cpp b/src/phasicFlow/streams/Stream/Istream.cpp
index cda6117e..be2ea9cb 100755
--- a/src/phasicFlow/streams/Stream/Istream.cpp
+++ b/src/phasicFlow/streams/Stream/Istream.cpp
@@ -858,6 +858,36 @@ pFlow::iIstream& pFlow::Istream::read
     return *this;
 }
 
+size_t pFlow::Istream::findBinaryBlockStart()
+{
+    size_t pos = 0;
+    char getChar = 'a';
+    unsigned char bFlag = 255;
+    int numFound = 0;
+
+    while( is_.good() && !is_.eof() )
+    {
+        getChar = is_.get();
+        pos++;
+
+        if( numFound <3 &&
+         static_cast<unsigned char>(getChar) == bFlag ) 
+        {
+            numFound++;
+        }
+        else if(numFound == 3 && static_cast<unsigned char>(getChar) == 0 )
+        {
+            return pos;
+        }
+        else
+        {
+            numFound = 0;
+        }
+    }
+
+    return static_cast<size_t>(-1);
+}
+
 void pFlow::Istream::rewind()
 {
     lineNumber_ = 1;      // Reset line number
diff --git a/src/phasicFlow/streams/Stream/Istream.hpp b/src/phasicFlow/streams/Stream/Istream.hpp
index 2b3b4879..2c324214 100755
--- a/src/phasicFlow/streams/Stream/Istream.hpp
+++ b/src/phasicFlow/streams/Stream/Istream.hpp
@@ -153,6 +153,8 @@ public:
         iIstream& read(char* buffer, std::streamsize count) override;
 
     
+    size_t findBinaryBlockStart()override;
+
     /// Rewind the stream so that it may be read again
     virtual void rewind();
 
diff --git a/src/phasicFlow/streams/Stream/Ostream.cpp b/src/phasicFlow/streams/Stream/Ostream.cpp
index aa89bd54..d3c74dbc 100755
--- a/src/phasicFlow/streams/Stream/Ostream.cpp
+++ b/src/phasicFlow/streams/Stream/Ostream.cpp
@@ -222,6 +222,13 @@ pFlow::iOstream& pFlow::Ostream::write(const double val)
     return *this;
 }
 
+pFlow::iOstream& pFlow::Ostream::write(const size_t val)
+{
+    os_ << val;
+    setState(os_.rdstate());
+    return *this;
+}
+
 pFlow::iOstream& pFlow::Ostream::write
 (
     const char* binaryData, 
@@ -239,11 +246,19 @@ pFlow::iOstream& pFlow::Ostream::write
     os_.write(binaryData, count);
     os_ << token::END_LIST;
 
-    setState(os_.rdstate());
-
     return *this;
 }
 
+pFlow::iOstream& pFlow::Ostream::writeBinaryBlockFlag()
+{
+    write( static_cast<char>(255));
+    write( static_cast<char>(255));
+    write( static_cast<char>(255));
+    write( static_cast<char>(0));
+    return *this;
+}
+
+
 void pFlow::Ostream::indent()
 {
     for (unsigned short i = 0; i < indentLevel_*indentSize_; ++i)
diff --git a/src/phasicFlow/streams/Stream/Ostream.hpp b/src/phasicFlow/streams/Stream/Ostream.hpp
index 25b2b09c..4d5f3205 100755
--- a/src/phasicFlow/streams/Stream/Ostream.hpp
+++ b/src/phasicFlow/streams/Stream/Ostream.hpp
@@ -124,9 +124,14 @@ public:
         /// Write double
         iOstream& write(const double val) override;
 
+        /// Write size_t
+        iOstream& write(const size_t val) override;
+
         /// Write a block of binray data 
         iOstream& write(const char* binaryData, std::streamsize count) override;
 
+        iOstream& writeBinaryBlockFlag() override;
+
         /// Add indentation characters
         void indent() override;
 
diff --git a/src/phasicFlow/streams/TStream/oTstream.cpp b/src/phasicFlow/streams/TStream/oTstream.cpp
index 2ef897ca..381dfaa1 100755
--- a/src/phasicFlow/streams/TStream/oTstream.cpp
+++ b/src/phasicFlow/streams/TStream/oTstream.cpp
@@ -155,6 +155,14 @@ pFlow::iOstream& pFlow::oTstream::write(const double val)
     return *this;
 }
 
+pFlow::iOstream& pFlow::oTstream::write(const size_t val)
+{
+    append(token(static_cast<int64>(val))); // tokenType::INT64
+
+    return *this;
+}
+
+
 pFlow::iOstream& pFlow::oTstream::write
 (
     const char* binaryData, 
diff --git a/src/phasicFlow/streams/TStream/oTstream.hpp b/src/phasicFlow/streams/TStream/oTstream.hpp
index c86b0090..23752a5d 100755
--- a/src/phasicFlow/streams/TStream/oTstream.hpp
+++ b/src/phasicFlow/streams/TStream/oTstream.hpp
@@ -129,6 +129,9 @@ public:
         /// Write double
         virtual iOstream& write(const double val) override;
 
+        /// Write double
+        virtual iOstream& write(const size_t val) override;
+
         /// Write a block of binray data
         iOstream& write(
             const char* binaryData, 
diff --git a/src/phasicFlow/streams/dataIO/dataIO.cpp b/src/phasicFlow/streams/dataIO/dataIO.cpp
new file mode 100644
index 00000000..9789c56e
--- /dev/null
+++ b/src/phasicFlow/streams/dataIO/dataIO.cpp
@@ -0,0 +1,612 @@
+
+#include "phasicFlowConfig.H"
+#ifdef pFlow_Build_MPI
+	#include <mpi.h>
+#endif
+
+#include <cstdio>
+#include <numeric>
+
+
+#include "dataIO.hpp"
+#include "fileSystem.hpp"
+
+
+#include "streams.hpp"
+
+static const size_t numcharFlag = 8;
+
+static inline unsigned char binaryFlag__[numcharFlag] = {255, 255, 255, 255, 255, 255, 255 ,0};
+
+inline 
+pFlow::uint64 chunkSizeOffeset(pFlow::uint64 procNo)
+{
+	// 1 is added for number of processors 
+	return (1+procNo)*sizeof(pFlow::uint64);
+}
+
+bool writeBinaryBlockFlagSTD(std::FILE* fh)
+{
+
+	std::size_t nw = std::fwrite( binaryFlag__ , sizeof(unsigned char), numcharFlag, fh );
+
+	if(nw < numcharFlag) 
+		return false;
+	else
+		return true;
+}
+
+pFlow::uint64 findBindaryBlockFlagSTD(std::FILE* fh)
+{
+	// get the current postion 
+	long fpos;
+	if( fpos = std::ftell( fh) ; fpos == -1L )
+	{
+		fatalErrorInFunction;
+		return -1;
+	}
+
+	pFlow::uint64 filePos = static_cast<pFlow::uint64>(fpos);
+
+	// start reading char by char 
+	int c; 
+	unsigned char ch;
+	int currPos = 0;
+    while ( std::fread(&ch, sizeof(ch), 1, fh) == 1 ) 
+    {
+   		if(std::ferror(fh)) return -1;	
+   		if(std::feof(fh))return -1;
+
+   		filePos++;
+
+   		if(ch == binaryFlag__[currPos] ) 
+   		{
+   			currPos++;
+
+   			if(currPos == numcharFlag) return filePos;
+   		}
+   		else
+   		{
+   			currPos = 0;
+   		}   
+    }
+
+    return -1;
+}
+
+
+#ifdef pFlow_Build_MPI
+pFlow::uint64 findBindaryBlockFlagMPI(MPI_File fh)
+{
+	// get the current postion 
+	MPI_Offset fpos;
+	if( MPI_File_get_position( fh, &fpos) != MPI_SUCCESS )
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+	pFlow::uint64 filePos = static_cast<pFlow::uint64>(fpos);
+
+	// start reading char by char 
+	unsigned char ch;
+	int currPos = 0;
+
+    while( MPI_File_read(
+    	fh, 
+    	&ch, 
+    	1, 
+    	MPI_UNSIGNED_CHAR, 
+    	MPI_STATUS_IGNORE) == MPI_SUCCESS ) 
+    {
+   		
+   		filePos++;
+   		
+   		if(ch == binaryFlag__[currPos] ) 
+   		{
+   			currPos++;
+
+   			if(currPos == numcharFlag) return filePos;
+   		}
+   		else
+   		{
+   			currPos = 0;
+   		}   
+    }
+
+    return -1;
+}
+
+#endif
+
+
+bool pFlow::dataIO::writeDataToFileEndSTD
+(
+	const fileSystem& filePath, 
+	const span<unsigned char>& data
+)
+{
+
+	if(!processors::isMaster()) return true;
+
+	// openfile 
+	word wFile = filePath.wordPath();
+	auto fh = std::fopen(wFile.c_str(), "ab");
+
+	if(!fh)
+	{
+		fatalErrorInFunction<<
+		"Error in Opening file "<< filePath <<endl;
+		std::fclose(fh);
+		return false;
+	}
+
+	// write the start of binary block flag
+	if(std::fseek(fh, 0 , SEEK_END)!=0)
+	{
+		fatalErrorInFunction<<
+		"error at reaching end of file "<<filePath<<endl;
+		std::fclose(fh);
+		return false;
+	}
+
+	if(!writeBinaryBlockFlagSTD(fh) )
+	{
+		fatalErrorInFunction<<
+		"Error in writing to file "<< filePath<<endl;
+		std::fclose(fh);
+		return false;
+	}
+		
+	uint64 numChunks = 1;	
+
+	// write number of data chunks
+	auto wc = std::fwrite(&numChunks, sizeof(numChunks), 1, fh);
+
+	if(wc < 1)
+	{
+		fatalErrorInFunction<<
+		"Error in writing numChunks to file "<< filePath <<endl;
+		std::fclose(fh);
+		return false;
+	}
+
+
+	// write size of each data chunk
+	uint64 sizeOfData = data.size();
+
+	wc = std::fwrite(&sizeOfData, sizeof(sizeOfData), 1, fh);
+	if(wc <1)
+	{
+		fatalErrorInFunction<<
+		"Error in writing size of data chunk to file "<< filePath <<endl;
+		std::fclose(fh);
+		return false;	
+	}
+
+	if(sizeOfData>0 )
+	{
+		// write data chunk to file
+		wc = std::fwrite(data.data(), sizeof(unsigned char), sizeOfData, fh);
+		if(wc < sizeOfData )
+		{
+			fatalErrorInFunction<<
+			"Error in writing size of data to file "<< filePath <<endl;
+			std::fclose(fh);
+			return false;	
+		}	
+	}
+
+	// close the file 
+	std::fclose(fh);
+
+	return true;
+}
+
+bool pFlow::dataIO::writeDataToFileEndMPI
+(
+	const fileSystem& filePath, 
+	const span<unsigned char>& data
+)
+{	
+
+#ifdef pFlow_Build_MPI
+	
+	// collect information from other processes
+	uint64 numProc 	= processors::globalSize();
+	uint64 thisSize = data.size();
+	uint64 offset;
+
+	CheckMPI
+	(
+		MPI_Scan(&thisSize, &offset, 1, MPI_UINT64_T,  MPI_SUM, MPI_COMM_WORLD),
+		true
+	);
+
+	word wFile = filePath.wordPath();
+
+	MPI_File fh;
+	
+	if( MPI_File_open(
+			MPI_COMM_WORLD,
+			wFile.c_str(),
+			MPI_MODE_WRONLY+MPI_MODE_APPEND,
+			MPI_INFO_NULL,
+			&fh) != MPI_SUCCESS)
+	{
+		fatalErrorInFunction<<
+		"Cannot open file "<< filePath<<endl;
+		return false;
+	}	
+
+	uint64 startPos;
+
+	if( processors::isMaster())
+	{
+		// set position to the end of file 
+		MPI_File_seek(fh, 0 , MPI_SEEK_END);		
+
+		if( MPI_File_write(
+			fh, 
+			binaryFlag__,
+			numcharFlag,
+			MPI_UNSIGNED_CHAR,
+			MPI_STATUS_IGNORE
+			) != MPI_SUCCESS )
+		{
+			fatalErrorInFunction<<
+			"Cannot write binary block flag into "<< filePath<<endl;
+			return false;
+		}
+		
+		MPI_Offset posOfBlock;
+		if( MPI_File_get_position(fh, &posOfBlock) != MPI_SUCCESS )
+		{
+			fatalErrorInFunction<<
+			"Cannot get the end pos of file "<< filePath<<endl;
+			return false;
+		}
+
+		
+		startPos = static_cast<uint64>(posOfBlock);
+	}
+
+		
+	if( MPI_Bcast(
+		&startPos, 
+		1, 
+		MPI_UINT64_T, 
+		processors::masterNo(), 
+		MPI_COMM_WORLD) != MPI_SUCCESS)
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+ 
+
+	if( processors::isMaster())
+	{
+		uint64 numProc = processors::globalSize();
+		if(MPI_File_write(
+			fh,
+			&numProc,
+			1,
+			MPI_UINT64_T,
+			MPI_STATUS_IGNORE) != MPI_SUCCESS)
+		{
+			fatalErrorInFunction<<
+			"Cannot write number of chunks into "<<filePath<<endl;
+			return false;
+		}
+	}
+
+	uint64 sizeOffset = startPos  + chunkSizeOffeset(processors::globalRank());
+
+	
+	if(MPI_File_write_at_all(
+		fh,
+		sizeOffset, 
+		&thisSize,  
+		1,
+		MPI_UINT64_T,
+		MPI_STATUS_IGNORE)!= MPI_SUCCESS)
+	{
+		fatalErrorInFunction<<
+		"Cannot write size of chunk into "<<filePath<<endl;
+		return false;
+	}
+
+	offset -= thisSize;
+
+	uint64 chunkOffset = startPos + chunkSizeOffeset(processors::globalSize()) + offset;
+
+	
+	if(MPI_File_write_at_all(
+		fh,
+		chunkOffset,
+		data.data(),
+		thisSize,
+		MPI_UNSIGNED_CHAR,
+		MPI_STATUS_IGNORE) != MPI_SUCCESS)
+	{
+		fatalErrorInFunction<<
+		"Cannot write data into "<<filePath<<endl;
+		return false;;
+	}
+
+	MPI_File_close(&fh);
+	MPI_Barrier(MPI_COMM_WORLD);
+	return true;
+
+#else
+	return writeDataToFileEndSTD(filePath, data);
+#endif
+
+	
+}
+
+
+bool pFlow::dataIO::readDataSTD
+(
+	const fileSystem& filePath,
+	const std::vector<uint64> chunkSizes,
+	span<unsigned char>& data,
+	uint64 binaryBlockStart
+)
+{
+	// sum of all chuncks
+	uint64 toRecv = std::accumulate(
+		chunkSizes.begin(), 
+		chunkSizes.end(), 
+		static_cast<uint64>(0));
+
+	if( data.size() != toRecv )
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+	word wFile = filePath.wordPath();
+	auto fh = std::fopen(wFile.c_str(), "rb");
+
+	if(!fh)
+	{
+		fatalErrorInFunction<<
+		"Error in Opening file "<< filePath<<endl;
+		return false;
+	}
+
+	// start of data chunks 
+	uint64 offset = binaryBlockStart + chunkSizeOffeset(chunkSizes.size());
+	
+	if(auto res = std::fseek(fh, offset, SEEK_SET); res!= 0 )
+	{
+		fatalErrorInFunction<<
+		"Error in file seek "<< filePath<<endl;
+		std::fclose(fh);
+		return false;
+	}
+
+	if(auto res = std::fread(
+			data.data(),
+			sizeof(unsigned char),
+			data.size(),
+			fh);
+		res!= data.size() )
+	{
+		fatalErrorInFunction<<
+		"Error in reading file "<< filePath<<endl;
+		std::fclose(fh);
+		return false;
+	}
+
+	std::fclose(fh);
+	return true;
+}
+
+
+bool pFlow::dataIO::readDataMPI
+(
+	const fileSystem& filePath, 
+	const std::vector<uint64> chunkSizes,
+	span<unsigned char>& data,
+	uint64 binaryBlockStart
+)
+{
+
+#ifdef pFlow_Build_MPI
+
+	if(chunkSizes.size() != processors::globalSize() )
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+	
+	word wFile = filePath.wordPath();
+	
+	MPI_File fh;
+
+	if( MPI_File_open(
+		MPI_COMM_WORLD,
+		wFile.c_str(),
+		MPI_MODE_RDONLY,
+		MPI_INFO_NULL,
+		&fh))
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+		
+	
+	auto procNo = processors::globalRank();
+
+	uint64 toRecv = chunkSizes[procNo];
+
+	// start of data chunks 
+	uint64 offset = binaryBlockStart + chunkSizeOffeset(processors::globalSize());
+
+	for(auto i=0; i<procNo; i++)
+	{
+		offset += chunkSizes[i];
+	}
+
+	if( data.size() != toRecv )
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+	
+	MPI_Status status;
+	
+	if( MPI_File_read_at_all(
+		fh,
+		offset,
+		data.data(),
+		data.size(),
+		MPI_UNSIGNED_CHAR,
+		&status) != MPI_SUCCESS )
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+#else
+
+#endif
+
+	return true;
+}
+
+
+bool pFlow::dataIO::readMetaMPI
+(
+	const fileSystem& filePath, 
+	std::vector<uint64>& chunkSizes,
+	uint64 &startPosBinaryBlock
+)
+{
+	word wFile = filePath.wordPath();
+
+#ifdef pFlow_Build_MPI	
+	MPI_File fh;
+
+	if(MPI_File_open(
+		MPI_COMM_WORLD, 
+		wFile.c_str(),
+		MPI_MODE_RDONLY,
+		MPI_INFO_NULL,
+		&fh) != MPI_SUCCESS)
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+	
+
+	uint64 startPos = findBindaryBlockFlagMPI(fh);
+	if( startPos == -1 )
+	{
+		fatalErrorInFunction;
+		return false;	
+	} 
+
+	startPosBinaryBlock  = startPos; 
+
+	uint64 numProcInFile;
+
+	if( MPI_File_read_at_all(
+		fh,
+		startPos,
+		&numProcInFile,
+		1,
+		MPI_UINT64_T,
+		MPI_STATUS_IGNORE) != MPI_SUCCESS)
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+	
+	chunkSizes.resize(numProcInFile);
+	
+	if(MPI_File_read_at_all(
+		fh,
+		startPos + sizeof(numProcInFile),
+		chunkSizes.data(), 
+		chunkSizes.size(), 
+		MPI_UINT64_T,
+		MPI_STATUS_IGNORE) != MPI_SUCCESS)
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+	MPI_File_close(&fh);
+
+#endif
+
+	return true;
+}
+
+bool pFlow::dataIO::readMetaSTD
+(
+	const fileSystem& filePath, 
+	std::vector<uint64>& chunkSizes,
+	uint64 &startPosBinaryBlock
+)
+{
+	// only on master
+	if( !processors::isMaster()) return true;
+
+	word wFile = filePath.wordPath();
+	std::FILE *fh = std::fopen(wFile.c_str(), "rb");
+
+	if(!fh)
+	{
+		fatalErrorInFunction<<
+		"Error in Opening file "<< filePath<<endl;
+		return false;
+	}
+
+	uint64 startPos = findBindaryBlockFlagSTD(fh);
+	if(startPos == -1 )
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+	startPosBinaryBlock = startPos;
+
+
+	std::fseek(fh, startPos, SEEK_SET); 
+
+	uint64 numProcInFile;
+	
+	auto res = std::fread(&numProcInFile, sizeof(numProcInFile), 1, fh);
+	
+
+	if(res != 1 )
+	{
+		fatalErrorInFunction<<
+		"Error in reading file "<< filePath<<endl;
+		std::fclose(fh);
+		return false;
+	}
+
+	
+	chunkSizes.resize(numProcInFile);
+
+	res = std::fread(chunkSizes.data(), sizeof(numProcInFile), numProcInFile, fh);
+
+	if(res!= numProcInFile)
+	{
+		fatalErrorInFunction<<
+		"Error in reading chunkSizes from file "<< filePath<<endl;
+		std::fclose(fh);
+		return false;
+	}
+
+	std::fclose(fh);
+	return true;
+}
+
+
diff --git a/src/phasicFlow/streams/dataIO/dataIO.hpp b/src/phasicFlow/streams/dataIO/dataIO.hpp
new file mode 100644
index 00000000..9a9d8453
--- /dev/null
+++ b/src/phasicFlow/streams/dataIO/dataIO.hpp
@@ -0,0 +1,103 @@
+
+#ifndef __dataIO_hpp__
+#define __dataIO_hpp__
+
+#include <vector>
+
+#include "types.hpp"
+#include "span.hpp"
+#include "processors.hpp"
+
+namespace pFlow
+{
+
+class fileSystem;
+
+
+class dataIO
+{
+
+public:
+
+	/**
+	 * Type of input/output
+	 * MasterProcessor: Read or write is done on master processor
+	 *    and the data on master processor is affected.
+	 * AllProcessorsDifferent: Read or write is done on all processors, 
+	 *    and each processor munipulates its own data.
+	 * AllProcessorsSimilar: Read is done on all processors but the read data
+	 *    is similar on processors. Write is done on master processor, since 
+	 *    all processors have the same copy of data. 
+	 */
+	enum IOType : int
+	{
+		MasterProcessor,
+		AllProcessorsDifferent,
+		AllProcessorsSimilar
+	};
+protected:
+
+	IOType IOType_;
+
+	
+	bool writeDataToFileEndSTD(
+		const fileSystem& filePath, 
+		const span<unsigned char>& data);
+
+	bool writeDataToFileEndMPI(
+		const fileSystem& filePath, 
+		const span<unsigned char>& data);
+	
+	bool readDataSTD(
+		const fileSystem& filePath,
+		const std::vector<uint64> chunkSizes,
+		span<unsigned char>& data,
+		uint64 binaryBlockStart);
+	
+	bool readDataMPI(
+		const fileSystem& filePath, 
+		const std::vector<uint64> chunkSizes,
+		span<unsigned char>& data,
+		uint64 binaryBlockStart);
+
+	bool readMetaMPI(
+		const fileSystem& filePath, 
+		std::vector<uint64>& chunkSizes,
+		uint64 &startPosBinaryBlock);
+
+	bool readMetaSTD(
+		const fileSystem& filePath, 
+		std::vector<uint64>& chunkSizes,
+		uint64 &startPosBinaryBlock);
+
+public:
+
+	dataIO(IOType ioT)
+	:
+		IOType_(ioT)
+	{}
+
+	~dataIO()=default;
+
+	template<typename T> 
+	bool writeDataEnd(
+		const fileSystem& filePath,
+		const span<T>& data);
+
+	template<typename T>
+	bool readData(
+		const fileSystem& filePath, 
+		std::vector<T>& data);
+
+	
+
+};
+
+}
+
+#include "dataIOTemplate.cpp"
+
+#endif
+
+
+
diff --git a/src/phasicFlow/streams/dataIO/dataIOTemplate.cpp b/src/phasicFlow/streams/dataIO/dataIOTemplate.cpp
new file mode 100644
index 00000000..b8450f97
--- /dev/null
+++ b/src/phasicFlow/streams/dataIO/dataIOTemplate.cpp
@@ -0,0 +1,112 @@
+
+template<typename T>
+bool pFlow::dataIO::writeDataEnd(
+	const fileSystem& filePath,
+	const span<T>& data)
+{
+
+	span<unsigned char> charSpan(
+		reinterpret_cast<unsigned char*> (const_cast<T*>(data.data())),
+		data.size()*sizeof(T));
+
+	switch (IOType_)
+	{
+		case MasterProcessor:
+		case AllProcessorsSimilar:
+		{
+			// This	means that only master processor should write
+			// in this case we perform write on master processor only
+			// this means that the data 
+			if(processors::isMaster())
+			{
+				return writeDataToFileEndSTD(filePath, charSpan);
+			}
+			else
+			{
+				return true;
+			}
+			break;
+		}
+		case AllProcessorsDifferent:
+		{
+			// This means that all processors should write their own
+			// copy of data 
+			return writeDataToFileEndMPI(filePath, charSpan);
+			break;
+		}
+	}
+
+	return false;
+}
+
+template<typename T>
+bool pFlow::dataIO::readData
+(
+	const fileSystem& filePath, 
+	std::vector<T>& data
+)
+{
+	std::vector<uint64> chunkSizes;
+	uint64 startPosBinaryBlock;
+
+	// read meta data
+	switch (IOType_)
+	{
+		case MasterProcessor:
+		{
+			if(!readMetaSTD(
+				filePath, 
+				chunkSizes, 
+				startPosBinaryBlock))
+			{
+				return false;
+			}
+			break;
+		}
+		case AllProcessorsDifferent:
+		case AllProcessorsSimilar:
+		{
+			if(!readMetaMPI(
+				filePath,
+				chunkSizes,
+				startPosBinaryBlock))
+			{
+				return false;
+			}
+			break;
+		}
+	}
+
+	data.clear();
+	if(IOType_ == MasterProcessor)
+	{
+		auto sizeOfData = std::accumulate( 
+			chunkSizes.begin(), 
+			chunkSizes.end(),
+			static_cast<uint64>(0));
+
+		data.resize(sizeOfData/sizeof(T));
+
+		span<unsigned char> charSpan(
+			reinterpret_cast<unsigned char*>(data.data()),
+			data.size()*sizeof(T));
+
+		readDataSTD(filePath, chunkSizes, charSpan, startPosBinaryBlock);
+	}
+	if( IOType_ == AllProcessorsDifferent )
+	{
+		auto thisProc = processors::globalRank();
+
+		data.resize(chunkSizes[thisProc]/sizeof(T));
+
+		span<unsigned char> charSpan(
+			reinterpret_cast<unsigned char*>(data.data()),
+			data.size()*sizeof(T));
+
+		std::cout<<"MPI part"<<std::endl;
+		readDataMPI(filePath, chunkSizes, charSpan, startPosBinaryBlock);
+
+	}
+	
+	return true;
+}
\ No newline at end of file
diff --git a/src/phasicFlow/streams/iStream/iIstream.cpp b/src/phasicFlow/streams/iStream/iIstream.cpp
index f91798b2..b909b999 100755
--- a/src/phasicFlow/streams/iStream/iIstream.cpp
+++ b/src/phasicFlow/streams/iStream/iIstream.cpp
@@ -75,6 +75,12 @@ bool pFlow::iIstream::peekBack(token& tok)
     return putBack_;
 }
 
+size_t pFlow::iIstream::findBinaryBlockStart()
+{
+    notImplementedFunction;
+    return 0;
+}
+
 bool pFlow::iIstream::findToken( const word & w )
 {
     rewind();
diff --git a/src/phasicFlow/streams/iStream/iIstream.hpp b/src/phasicFlow/streams/iStream/iIstream.hpp
index 90a2f81d..7f0ddb43 100755
--- a/src/phasicFlow/streams/iStream/iIstream.hpp
+++ b/src/phasicFlow/streams/iStream/iIstream.hpp
@@ -128,8 +128,8 @@ public:
 	    /// Read a doubleScalar
 	    virtual iIstream& read(double&) = 0;
 
+	    /// Write a block of binray data 
 	    virtual iIstream& read(char* buffer, std::streamsize count) =0;
-
 	    
 	    /// Rewind the stream so that it may be read again
 	    virtual void rewind() = 0;
@@ -137,6 +137,11 @@ public:
 
     ////- find and lookups
 
+	    /// It seek for a character sequence that indicates the start of 
+	    /// a binary block
+	    /// char sequence is 255 255 255 0 
+	    virtual size_t findBinaryBlockStart();
+
         /// search for all tokesn and find the first word token tbat matchs w
         virtual bool findToken( const word & w );
 
diff --git a/src/phasicFlow/streams/iStream/iOstream.cpp b/src/phasicFlow/streams/iStream/iOstream.cpp
index 622ff494..dfe82685 100644
--- a/src/phasicFlow/streams/iStream/iOstream.cpp
+++ b/src/phasicFlow/streams/iStream/iOstream.cpp
@@ -21,9 +21,16 @@ Licence:
 
 #include "iOstream.hpp"
 #include "token.hpp"
+#include "error.hpp"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+pFlow::iOstream& pFlow::iOstream::writeBinaryBlockFlag()
+{
+    notImplementedFunction;
+    return *this;
+}
+
 void pFlow::iOstream::decrIndent()
 {
     if (!indentLevel_)
diff --git a/src/phasicFlow/streams/iStream/iOstream.hpp b/src/phasicFlow/streams/iStream/iOstream.hpp
index 76141a01..1a82f8b1 100644
--- a/src/phasicFlow/streams/iStream/iOstream.hpp
+++ b/src/phasicFlow/streams/iStream/iOstream.hpp
@@ -135,8 +135,14 @@ public:
         /// Write double
         virtual iOstream& write(const double val) = 0;
        
+        /// Write size_t
+        virtual iOstream& write(const size_t val) = 0;
+         
         /// Write a block of binray data 
         virtual iOstream& write(const char* binaryData, std::streamsize count) = 0;
+
+        /// Write the flag to indicate the start of a binary block
+        virtual iOstream& writeBinaryBlockFlag();
         
     
     //// - Indent 
@@ -411,6 +417,11 @@ inline iOstream& operator<<( iOstream& os, const double& val)
     return os.write(val);
 }
 
+inline iOstream& operator<<( iOstream& os, const size_t& val)
+{
+    return os.write(val);
+}
+
 // Useful aliases for tab and newline characters
 constexpr char tab = '\t';
 constexpr char nl = '\n';
diff --git a/src/phasicFlow/streams/masterOstream/masterOstream.cpp b/src/phasicFlow/streams/masterOstream/masterOstream.cpp
index 6c02c7a2..24d0931a 100755
--- a/src/phasicFlow/streams/masterOstream/masterOstream.cpp
+++ b/src/phasicFlow/streams/masterOstream/masterOstream.cpp
@@ -153,6 +153,15 @@ pFlow::iOstream& pFlow::masterOstream::write(const double val)
     return *this; 
 }
 
+pFlow::iOstream& pFlow::masterOstream::write(const size_t val)
+{
+    if(showOutput())
+    {
+        Ostream::write(val);
+    }
+    return *this; 
+}
+
 
 pFlow::iOstream& pFlow::masterOstream::write
 (
diff --git a/src/phasicFlow/streams/masterOstream/masterOstream.hpp b/src/phasicFlow/streams/masterOstream/masterOstream.hpp
index 6568bb0d..60931ea4 100755
--- a/src/phasicFlow/streams/masterOstream/masterOstream.hpp
+++ b/src/phasicFlow/streams/masterOstream/masterOstream.hpp
@@ -64,11 +64,20 @@ public:
 
 
     //// - Methods
+        /// Set if this processor is slave or master 
         void setMasterSlave(bool master)
         {
             isThisMaster_ = master;
         }
 
+        /// Write token to stream or otherwise handle it.
+        /// return false if the token type was not handled by this method
+        bool write(const token& tok)
+        // to prevent compiler warning, this method is overriden agian 
+        {
+            return Ostream::write(tok);
+        }
+
         /// Write character
         iOstream& write(const char c)override;
 
@@ -105,6 +114,9 @@ public:
 
         /// Write double
         iOstream& write(const double val) override;
+
+        /// Write double
+        iOstream& write(const size_t val) override;
         
         /// Write a block of binray data 
         iOstream& write(const char* binaryData, std::streamsize count) override;
diff --git a/src/phasicFlow/streams/processorOstream/processorOstream.cpp b/src/phasicFlow/streams/processorOstream/processorOstream.cpp
index 5f39d158..b91950cd 100755
--- a/src/phasicFlow/streams/processorOstream/processorOstream.cpp
+++ b/src/phasicFlow/streams/processorOstream/processorOstream.cpp
@@ -134,6 +134,11 @@ pFlow::iOstream& pFlow::processorOstream::write(const double val)
     return Ostream::write(val);
 }
 
+pFlow::iOstream& pFlow::processorOstream::write(const size_t val)
+{
+    checkForPrefix();
+    return Ostream::write(val);
+}
 
 pFlow::iOstream& pFlow::processorOstream::write
 (
diff --git a/src/phasicFlow/streams/processorOstream/processorOstream.hpp b/src/phasicFlow/streams/processorOstream/processorOstream.hpp
index 6f6dd7e6..d448f082 100755
--- a/src/phasicFlow/streams/processorOstream/processorOstream.hpp
+++ b/src/phasicFlow/streams/processorOstream/processorOstream.hpp
@@ -85,6 +85,14 @@ public:
             printPrefix_ = true;
         }
         
+        /// Write token to stream or otherwise handle it.
+        /// return false if the token type was not handled by this method
+        bool write(const token& tok)
+        // to prevent compiler warning, this method is overriden agian 
+        {
+            return Ostream::write(tok);
+        }
+
         /// Write character
         iOstream& write(const char c)override;
 
@@ -121,6 +129,9 @@ public:
 
         /// Write double
         iOstream& write(const double val) override;
+
+        /// Write double
+        iOstream& write(const size_t val) override;
         
         /// Write a block of binray data 
         iOstream& write(const char* binaryData, std::streamsize count) override;
diff --git a/src/phasicFlow/types/basicTypes/builtinTypes.hpp b/src/phasicFlow/types/basicTypes/builtinTypes.hpp
index 340179ec..d87a0a73 100755
--- a/src/phasicFlow/types/basicTypes/builtinTypes.hpp
+++ b/src/phasicFlow/types/basicTypes/builtinTypes.hpp
@@ -64,6 +64,8 @@ using size_t  = std::size_t;
 
 using word    = std::string;
 
+inline const int  numBytesForReal__ = sizeof(real);
+
 inline 
 word floatingPointDescription()
 {
diff --git a/src/phasicFlow/types/basicTypes/math.hpp b/src/phasicFlow/types/basicTypes/math.hpp
index 0ac62f6c..cd61e2d3 100755
--- a/src/phasicFlow/types/basicTypes/math.hpp
+++ b/src/phasicFlow/types/basicTypes/math.hpp
@@ -23,11 +23,15 @@ Licence:
 
 #include <cmath>
 
+#include "pFlowMacros.hpp"
+#include "builtinTypes.hpp"
+
+
 #ifdef __CUDACC__
 #include "math.h"
 #endif
 
-#include "builtinTypes.hpp"
+
 
 //* * * * * * * * * * *  List of functinos * * * * * * * * //
 // abs, mod, exp, log, log10, pow, sqrt, cbrt