diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3ce683a4..e411fe79 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -74,7 +74,7 @@ add_subdirectory(src)
 
 #add_subdirectory(solvers)
 
-#add_subdirectory(utilities)
+add_subdirectory(utilities)
 
 #add_subdirectory(DEMSystems)
 add_subdirectory(testIO)
diff --git a/src/phasicFlow/CMakeLists.txt b/src/phasicFlow/CMakeLists.txt
index eb52ef39..2129f985 100644
--- a/src/phasicFlow/CMakeLists.txt
+++ b/src/phasicFlow/CMakeLists.txt
@@ -37,20 +37,35 @@ containers/Vector/Vectors.cpp
 containers/Field/Fields.cpp
 containers/List/anyList/anyList.cpp
 
+Timer/Timer.cpp
+Timer/Timers.cpp
+
 repository/IOobject/objectFile.cpp
 repository/IOobject/IOfileHeader.cpp
 repository/IOobject/IOobject.cpp
 repository/IOobject/IOPattern.cpp
+repository/repository/repository.cpp
+repository/Time/Time.cpp
+repository/Time/timeControl.cpp
+repository/systemControl/systemControl.cpp
+repository/systemControl/dynamicLinkLibs.cpp
 
 eventManagement/subscriber.cpp
 eventManagement/observer.cpp
 
+structuredData/box/box.cpp
 structuredData/line/line.cpp
 structuredData/infinitePlane/infinitePlane.cpp
 structuredData/plane/plane.cpp
 structuredData/domain/domain.cpp
+structuredData/domain/simulationDomain.cpp
+structuredData/domain/regularSimulationDomain.cpp
+structuredData/pointStructure/boundaryList.cpp
 structuredData/pointStructure/internalPoints.cpp
+structuredData/pointStructure/pointStructure.cpp
 structuredData/boundaries/boundaryBase/boundaryBase.cpp
+structuredData/boundaries/boundaryExit/boundaryExit.cpp
+
 
 )
 
diff --git a/src/phasicFlow/containers/Field/Field.hpp b/src/phasicFlow/containers/Field/Field.hpp
index d183b5c8..aafbfbd2 100644
--- a/src/phasicFlow/containers/Field/Field.hpp
+++ b/src/phasicFlow/containers/Field/Field.hpp
@@ -48,17 +48,17 @@ public:
 
 	using iterator        = typename VectorType::iterator;
 
-  	using constIterator   = typename VectorType::constIterator;
+	using const_iterator   = typename VectorType::const_iterator;
 
 	using reference       = typename VectorType::reference;
   	
-  	using constReference  = typename VectorType::constReference;
+	using const_reference  = typename VectorType::const_reference;
 
-	using valueType       = typename VectorType::valueType;
+	using value_type       = typename VectorType::value_type;
   	
-  	using pointer         = typename VectorType::pointer;
+	using pointer         = typename VectorType::pointer;
   	
-  	using constPointer    = typename VectorType::constPointer;
+	using const_pointer    = typename VectorType::const_pointer;
   	
 
 protected:
diff --git a/src/phasicFlow/containers/List/ListPtr/ListPtr.hpp b/src/phasicFlow/containers/List/ListPtr/ListPtr.hpp
index ed150eed..fe1f1850 100644
--- a/src/phasicFlow/containers/List/ListPtr/ListPtr.hpp
+++ b/src/phasicFlow/containers/List/ListPtr/ListPtr.hpp
@@ -30,7 +30,6 @@ Licence:
 #include "error.hpp"
 #include "iOstream.hpp"
 
-
 namespace pFlow
 {
 
@@ -91,7 +90,7 @@ public:
 		// - a list with initial length of len
 		ListPtr(size_t len)
 		:
-			list_(len)
+			list_(len, nullptr)
 		{}
 
 		// - copy construct, create new objects out of the pointers in the src 
@@ -141,10 +140,10 @@ public:
 	//// - Methods 
 
 		// - set the ith element 
-		T* set(size_t i, T* ptr);
+		uniquePtr<T> set(size_t i, T* ptr);
 		
 		// - set the ith element and take the ownership from uniquePtr
-		uniquePtr<T> set(size_t 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
@@ -155,7 +154,7 @@ public:
 	    void push_back(T* ptr);
 	    
 	    // - put the pointer at the end
-	    void push_back(uniquePtr<T>& ptr);
+	    void push_back(uniquePtr<T>&& ptr);
 	    
 	    // - safely create (in-place) and put the pointer at the end
 	    template<typename... Args>
diff --git a/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp b/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp
index 57f42c2b..77ef5a99 100644
--- a/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp
+++ b/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp
@@ -144,7 +144,7 @@ pFlow::ListPtr<T>& pFlow::ListPtr<T>::operator=
 }
 
 template<typename T>
-T* pFlow::ListPtr<T>::set
+pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
 (
 	size_t i, T* ptr
 )
@@ -158,7 +158,7 @@ template<typename T>
 pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
 (
 	size_t i,
-	uniquePtr<T>& ptr
+	uniquePtr<T>&& ptr
 )
 {
 	if( i > size() )
@@ -171,10 +171,9 @@ pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
 
 	auto iter = list_.begin();
 	std::advance(iter, i);
-	auto oldIter = iter;
+	T* oldPtr = *iter;
 	*iter = ptr.release();
-	return *oldIter;
-
+	return uniquePtr<T>(oldPtr);
 }
 
 
@@ -200,7 +199,7 @@ void pFlow::ListPtr<T>::push_back
 }
 
 template<typename T>
-void pFlow::ListPtr<T>::push_back(uniquePtr<T>& ptr)
+void pFlow::ListPtr<T>::push_back(uniquePtr<T>&& ptr)
 {
 	list_.push_back( ptr.release() );
 }
@@ -271,17 +270,23 @@ pFlow::uniquePtr<T> pFlow::ListPtr<T>::release
 	return p;
 }
 
+
+
 template<typename T>
 void pFlow::ListPtr<T>::clear()
 {
+	
+	int i =0;
 	for( auto iter = list_.begin(); iter != list_.end(); ++iter )
 	{
 		if(*iter != nullptr)
 		{
+	
 			delete *iter;
 			*iter = nullptr;
-		}
+		}	
 	}
+
 	list_.clear();
 }
 
diff --git a/src/phasicFlow/containers/Vector/Vector.hpp b/src/phasicFlow/containers/Vector/Vector.hpp
index b025a148..05275ce5 100644
--- a/src/phasicFlow/containers/Vector/Vector.hpp
+++ b/src/phasicFlow/containers/Vector/Vector.hpp
@@ -65,19 +65,19 @@ public:
 
 	typedef typename vectorType::iterator 			iterator;
 	
-	typedef typename vectorType::const_iterator 	constIterator;
+	typedef typename vectorType::const_iterator 	const_iterator;
 	
 	typedef typename vectorType::reference 			reference;
 	
-	typedef typename vectorType::const_reference 	constReference;
+	typedef typename vectorType::const_reference 	const_reference;
 	
-	typedef  T 										valueType; 		
+	typedef  T 										value_type; 		
 
 	typedef  T* 									pointer;
 
-	typedef const T*  								constPointer;
+	typedef const T*  								const_pointer;
 
-	typedef typename std::initializer_list<T> 		initList;
+	typedef typename std::initializer_list<T> 		init_list;
 
 protected:
 
@@ -139,7 +139,7 @@ public:
 
 
 		/// Vector from name and initializer list 
-		inline Vector(const word& name, const initList &l)
+		inline Vector(const word& name, const init_list &l)
 		:
 			vectorType(l),
 			name_(name)
@@ -325,11 +325,25 @@ public:
 	// - fill the whole content of vector, [begin, end), with val 
 	inline void fill( const T& val);
 
+	inline
+	auto getSpan()
+	{
+		return span<T>(this->data(), this->size());
+	}
+	
+	inline 
+	auto getSpan()const
+	{
+		return span<const T>(this->data(), this->size());
+	}
+
 	static constexpr bool isHostAccessible()
 	{
 		return isHostAccessible_;
 	}
 
+	
+
 	inline void operator +=( const T& val);
 	inline void operator -=( const T& val);
 	inline void operator *=( const T& val);
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
index 46b461b3..4e3b3d2d 100644
--- a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
+++ b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
@@ -62,27 +62,27 @@ public:
 
 		using iterator        = T*;
 
-		using constIterator   = const T*;
+		using const_iterator   = const T*;
 	  	
 		using reference       = T&;
 	  	
-	 	using constReference  = const T&;
+	 	using const_reference  = const T&;
 
-		using valueType       = T;
+		using value_type       = T;
 	  	
 		using pointer         = T*;
 	  	
-		using constPointer    = const T*;
+		using const_pointer    = const T*;
   	
 	//- typedefs related to memory management 
 
 		using viewType 			= ViewType1D<T, MemorySpace>;
 
-		using deviceType 		= typename viewType::device_type;
+		using device_type 		= typename viewType::device_type;
 
-		using memorySpace 		= typename viewType::memory_space;
+		using memory_space 		= typename viewType::memory_space;
 
-		using executionSpace 	= typename viewType::execution_space;
+		using execution_space 	= typename viewType::execution_space;
 
 protected:
 	
@@ -102,19 +102,19 @@ protected:
 
 		/// Is the memory of this vector accessible from Host
 		static constexpr 
-		bool  isHostAccessible_ = isHostAccessible<executionSpace>();
-			//Kokkos::SpaceAccessibility<executionSpace,HostSpace>::accessible;
+		bool  isHostAccessible_ = isHostAccessible<execution_space>();
+			//Kokkos::SpaceAccessibility<execution_space,HostSpace>::accessible;
 
 		/// Is the memory of this vector accessiple from Divce
 		static constexpr
-		bool isDeviceAccessible_ = isDeviceAccessible<executionSpace>();
+		bool isDeviceAccessible_ = isDeviceAccessible<execution_space>();
 
 
 		/// Name of the memory space
 	  	constexpr static inline 
 	  	const char* memoerySpaceName()
 	  	{
-	  		return memorySpace::name();
+	  		return memory_space::name();
 	  	}
 
 	  	/// Evaluate capacity based on the input size 
@@ -461,6 +461,17 @@ public:
 			assign(src, src.capacity());
 		}
 
+		INLINE_FUNCTION_H
+		auto getSpan()
+		{
+			return span<T>(view_.data(), size());
+		}
+		
+		INLINE_FUNCTION_H 
+		auto getSpan()const
+		{
+			return span<const T>(view_.data(), this->size());
+		}
 
 		INLINE_FUNCTION_H
 		bool insertSetElement(uint32IndexContainer indices, const T& val);
@@ -520,7 +531,7 @@ public:
 		}
 
 		INLINE_FUNCTION_H 
-		constPointer data()const{
+		const_pointer data()const{
 			return view_.data();
 		}
 
@@ -535,7 +546,7 @@ public:
 		/// Return begin iterator. it works when host is accessible.
 		template<bool Enable = true>
 		INLINE_FUNCTION_H 
-		typename std::enable_if_t<isHostAccessible_ && Enable,constIterator>
+		typename std::enable_if_t<isHostAccessible_ && Enable,const_iterator>
 		begin()const {
 			return data();
 		}
@@ -553,7 +564,7 @@ public:
 		/// Return end iterator. it works when host is accessible.
 		template<bool Enable = true>
 		INLINE_FUNCTION_H
-		typename std::enable_if_t<isHostAccessible_ && Enable,constIterator>
+		typename std::enable_if_t<isHostAccessible_ && Enable,const_iterator>
 		end()const{
 			return size_ > 0 ? data() + size_: data();
 		}
@@ -569,7 +580,7 @@ public:
 		/// Return reference to element i. it works when host is accessible.
 		template<bool Enable = true>
 		INLINE_FUNCTION_H 
-		typename std::enable_if_t<isHostAccessible_ && Enable,constReference>
+		typename std::enable_if_t<isHostAccessible_ && Enable,const_reference>
 		 operator[](size_t i)const{
 			return view_[i];
 		}
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingleI.hpp b/src/phasicFlow/containers/VectorHD/VectorSingleI.hpp
index 61217a96..2dabf14a 100644
--- a/src/phasicFlow/containers/VectorHD/VectorSingleI.hpp
+++ b/src/phasicFlow/containers/VectorHD/VectorSingleI.hpp
@@ -32,7 +32,7 @@ bool pFlow::VectorSingle<T,MemorySpace>::insertSetElement(uint32IndexContainer i
 		resize(maxInd+1);
 	}
 
-	using policy = Kokkos::RangePolicy<executionSpace,Kokkos::IndexType<uint32>>;
+	using policy = Kokkos::RangePolicy<execution_space,Kokkos::IndexType<uint32>>;
 
 	if constexpr( isDeviceAccessible_ )
 	{
@@ -88,7 +88,7 @@ bool pFlow::VectorSingle<T,MemorySpace>::insertSetElement
 		resize(maxInd+1);
 	} 
 	
-	using policy = Kokkos::RangePolicy<executionSpace,Kokkos::IndexType<uint32>>;
+	using policy = Kokkos::RangePolicy<execution_space,Kokkos::IndexType<uint32>>;
 
 	hostViewType1D<const T> hVals( vals.data(), vals.size());
 
@@ -153,7 +153,7 @@ bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(uint32IndexContainer indic
 
 	viewType 	sortedView(this->name(), newSize);
 
-	using policy = Kokkos::RangePolicy< executionSpace,Kokkos::IndexType<uint32>>;
+	using policy = Kokkos::RangePolicy< execution_space,Kokkos::IndexType<uint32>>;
 
 	if constexpr( isDeviceAccessible_)
 	{
diff --git a/src/phasicFlow/dictionary/dictionary.cpp b/src/phasicFlow/dictionary/dictionary.cpp
index a52d1d76..156f99a2 100644
--- a/src/phasicFlow/dictionary/dictionary.cpp
+++ b/src/phasicFlow/dictionary/dictionary.cpp
@@ -378,7 +378,8 @@ bool pFlow::dictionary::isFileDict()const
 bool  pFlow::dictionary::addPtr
 (
 	const word& keyword,
-	uniquePtr<iEntry>& entry
+	uniquePtr<iEntry>& entry,
+	bool warning
 )
 {
 
@@ -390,9 +391,13 @@ bool  pFlow::dictionary::addPtr
 	// search all entries for repeated keyword
 	if(auto [ptr, exist] = entries_.find(keyword); exist )
 	{
-		warningInFunction <<
-		"keyword " << keyword << " already exists in the dicrionary " <<
-		this->globalName() << ". The old entry will be replaced by the new one. \n";
+		if(warning)
+		{
+			warningInFunction <<
+			"keyword " << keyword << " already exists in the dicrionary " <<
+			this->globalName() << ". The old entry will be replaced by the new one. \n";
+		}
+		
 		// store the old pointer to entry
 		oldEntryPtr = ptr;
 	}
@@ -400,6 +405,7 @@ bool  pFlow::dictionary::addPtr
 	
 	if( entries_.insertReplace(keyword, entry) )
 	{
+	
 		if(oldEntryPtr)
 		{
 			// this should be replaced 
diff --git a/src/phasicFlow/dictionary/dictionary.hpp b/src/phasicFlow/dictionary/dictionary.hpp
index 661882d5..fcd55ba4 100644
--- a/src/phasicFlow/dictionary/dictionary.hpp
+++ b/src/phasicFlow/dictionary/dictionary.hpp
@@ -175,9 +175,10 @@ public:
 		/// if dictionary is file dictionary, return false
 		virtual bool isFileDict()const;
 
-		/// add a pointer entry (dictionary/dataEntry)
-		//   replaces this entry with existing one
-		bool addPtr(const word& keyword, uniquePtr<iEntry>& etry );
+		/// @brief add a pointer entry (dictionary/dataEntry)
+		///   replaces this entry with existing one and issue a warning
+		bool addPtr(const word& keyword, uniquePtr<iEntry>& etry, bool warning = true );
+
 
 		/// add a float dataEntry 
 		bool add(const word& keyword, const float& v);
@@ -206,13 +207,16 @@ public:
 		/// add a uint8 dataEntry
 		bool add(const word& keyword, const uint8& v);
 
-		// add a dictionary with the specifiedd keyword
+		/// @brief add a dictionary with the specifiedd keyword, if 
+		/// it exists, replace it.
 		bool addDict(const word& keyword, const dictionary& dict);
 
 		/// add a dataEntry of type T
 		template<typename T>
 		bool add(const word& keyword, const T& v );
 
+		template<typename T>
+		bool addOrKeep(const word& keyword, const T& v);
 
 		void clear();
 
@@ -311,6 +315,19 @@ bool dictionary::add(const word& keyword, const T& v )
 }
 
 
+template<typename T>
+bool dictionary::addOrKeep(const word& keyword, const T& v )
+{
+	if(!containsDataEntry(keyword))
+	{
+		uniquePtr<iEntry> ptr = makeUnique<dataEntry>(keyword, *this ,v);
+		return addPtr(keyword, ptr);
+	}
+	
+	return true;
+}
+
+
 template<typename T>
 bool dictionary::readDataEntry
 (
diff --git a/src/phasicFlow/eventManagement/eventManagement.hpp b/src/phasicFlow/eventManagement/eventManagement.hpp
new file mode 100644
index 00000000..3b07b2be
--- /dev/null
+++ b/src/phasicFlow/eventManagement/eventManagement.hpp
@@ -0,0 +1,16 @@
+#ifndef __eventManagement_hpp__
+#define __eventManagement_hpp__
+
+#include "subscriber.hpp"
+#include "observer.hpp"
+#include "pointMessage.hpp"
+
+namespace pFlow
+{
+
+using pointSubcriber = subscriber<pointMessage>
+
+}
+
+
+#endif
\ No newline at end of file
diff --git a/src/phasicFlow/eventManagement/subscriber.hpp b/src/phasicFlow/eventManagement/subscriber.hpp
index 3c424a3a..48090597 100644
--- a/src/phasicFlow/eventManagement/subscriber.hpp
+++ b/src/phasicFlow/eventManagement/subscriber.hpp
@@ -40,11 +40,23 @@ protected:
 	// - list of subsribed objectd that recieve updage messages 
 	mutable std::array<List<observer*>,16> observerList_;
 
+	word subName_;
+
 public:
 
-	subscriber()
+	subscriber(const word& name)
+	:
+		subName_(name)
 	{}
 
+	subscriber(const subscriber&) = delete;
+
+	subscriber(subscriber&&) = default;
+
+	subscriber& operator = (const subscriber&) = delete;
+
+	subscriber& operator = (subscriber&&) = default;
+
 	virtual ~subscriber();
 
 	virtual bool subscribe(message msg, observer* obsevr)const;
@@ -55,6 +67,7 @@ public:
 
 	//bool notify(const eventMessage& msg, const List<eventObserver*>& exclutionList );
 
+
 	
 	bool notify(const message msg, const anyList& varList);
 	
diff --git a/src/phasicFlow/globals/vocabs.hpp b/src/phasicFlow/globals/vocabs.hpp
index 8ac68538..5afac001 100755
--- a/src/phasicFlow/globals/vocabs.hpp
+++ b/src/phasicFlow/globals/vocabs.hpp
@@ -37,6 +37,7 @@ const inline char* integrationFolder__  	= "integration";
 
 // file names
 const inline char* settingsFile__			= "settingsDict";
+const inline char* domainFile__     = "domainDict";
 const inline char* insertionFile__ 			= "particleInsertion";
 const inline char* sphereShapeFile__  		= "sphereShape";
 const inline char* pointStructureFile__ 	= "pStructure";
diff --git a/src/phasicFlow/processors/processors.hpp b/src/phasicFlow/processors/processors.hpp
index 494787a6..98d9314c 100644
--- a/src/phasicFlow/processors/processors.hpp
+++ b/src/phasicFlow/processors/processors.hpp
@@ -95,6 +95,18 @@ public:
 		return processors::globalSize()>1;
 	}
 	
+	static inline
+	const char* runTypeName()
+	{
+		if(isParallel())
+		{
+			return "MPI";
+		}
+		else
+		{
+			return "regular";
+		}
+	}
 	/// Is MPI initialized?
 	static
 	bool isInitialized();
diff --git a/src/phasicFlow/ranges/stridedRange.hpp b/src/phasicFlow/ranges/stridedRange.hpp
index 5ecf7367..ca1cf7ba 100644
--- a/src/phasicFlow/ranges/stridedRange.hpp
+++ b/src/phasicFlow/ranges/stridedRange.hpp
@@ -76,9 +76,9 @@ public:
 
 	stridedRange(const dictionary& dict)
 	:
-		begin_(dict.getVal<label>("begin")),
-		end_(dict.getVal<label>("end")),
-		stride_(dict.getVal<label>("stride"))
+		begin_(dict.getVal<T>("begin")),
+		end_(dict.getVal<T>("end")),
+		stride_(dict.getVal<T>("stride"))
 	{}
 		
 	inline 
diff --git a/src/phasicFlow/repository/Time/Time.cpp b/src/phasicFlow/repository/Time/Time.cpp
index e99679ba..e4367cee 100644
--- a/src/phasicFlow/repository/Time/Time.cpp
+++ b/src/phasicFlow/repository/Time/Time.cpp
@@ -111,7 +111,7 @@ bool pFlow::Time::write
 {
 	if(outputToFile())
 	{
-		REPORT(0)<<"\nWriting to file at time: "<< cyanText(timeName())<<endREPORT;
+		REPORT(0)<<"\nWriting to file at time: "<< Cyan_Text(timeName())<<END_REPORT;
 		return repository::write(verbose);
 	}
 	return true;
diff --git a/src/phasicFlow/repository/Time/Time.hpp b/src/phasicFlow/repository/Time/Time.hpp
index 915c9478..f1887b38 100644
--- a/src/phasicFlow/repository/Time/Time.hpp
+++ b/src/phasicFlow/repository/Time/Time.hpp
@@ -72,7 +72,7 @@ public:
 	//// -  Methods
 		virtual fileSystem localPath()const
 		{
-			return timeName();
+			return fileSystem(timeName());
 		}
 	
 		// - geometry repository 
diff --git a/src/phasicFlow/repository/Time/timeControl.cpp b/src/phasicFlow/repository/Time/timeControl.cpp
index cb6a162a..d9182701 100644
--- a/src/phasicFlow/repository/Time/timeControl.cpp
+++ b/src/phasicFlow/repository/Time/timeControl.cpp
@@ -203,7 +203,7 @@ bool pFlow::timeControl::operator ++(int)
 	currentTime_ += dt_;
 	if(screenReport() && !managedExternaly_)
 	{
-		REPORT(0)<<"Time (s): "<<cyanText( currentTimeWord() )<<endREPORT;
+		REPORT(0)<<"Time (s): "<<Cyan_Text( currentTimeWord() )<<END_REPORT;
 	}
 	// switch outputToFile_ on/off
 	checkForOutputToFile();
diff --git a/src/phasicFlow/repository/repository/repository.cpp b/src/phasicFlow/repository/repository/repository.cpp
index 520afa38..edbf45ba 100644
--- a/src/phasicFlow/repository/repository/repository.cpp
+++ b/src/phasicFlow/repository/repository/repository.cpp
@@ -245,7 +245,7 @@ bool pFlow::repository::write
 	{
 		if(verbose)
 		{
-			REPORT(1)<< "Writing to " << obj.second.path()<<endREPORT;
+			REPORT(1)<< "Writing to " << obj.second.path()<<END_REPORT;
 		}
 		
 		if(!obj.second.write())
diff --git a/src/phasicFlow/repository/systemControl/dynamicLinkLibs.cpp b/src/phasicFlow/repository/systemControl/dynamicLinkLibs.cpp
index 6004dd77..f74b4303 100644
--- a/src/phasicFlow/repository/systemControl/dynamicLinkLibs.cpp
+++ b/src/phasicFlow/repository/systemControl/dynamicLinkLibs.cpp
@@ -37,7 +37,7 @@ pFlow::dynamicLinkLibs::dynamicLinkLibs(
 	}
 	
 
-	REPORT(1)<< "libs are "<< greenText(libNames)<<endREPORT;
+	REPORT(1)<< "libs are "<< Green_Text(libNames)<<END_REPORT;
 
 	for(auto libName:libNames)
 	{
@@ -68,7 +68,7 @@ pFlow::dynamicLinkLibs::~dynamicLinkLibs()
 
 void* pFlow::dynamicLinkLibs::open(word libName)
 {
-	REPORT(2)<<"Loading "<< greenText(libName)<<endREPORT;
+	REPORT(2)<<"Loading "<< Green_Text(libName)<<END_REPORT;
 
 	void* handle = dlopen(libName.c_str(), RTLD_LAZY|RTLD_GLOBAL);
 	
diff --git a/src/phasicFlow/repository/systemControl/systemControl.cpp b/src/phasicFlow/repository/systemControl/systemControl.cpp
index 4cd60260..974567fb 100644
--- a/src/phasicFlow/repository/systemControl/systemControl.cpp
+++ b/src/phasicFlow/repository/systemControl/systemControl.cpp
@@ -137,7 +137,7 @@ pFlow::systemControl::systemControl
 	),
 	libs_(settingsDict_),
 	outFilePrecision_(
-		settingsDict_.getValOrSet("outFilePrecision", static_cast<size_t>(6))
+		settingsDict_.getValOrSet("outFilePrecision", static_cast<uint64>(6))
 		),
 	Time_
 	(
@@ -231,6 +231,22 @@ pFlow::systemControl::systemControl(
 	writeToFileTimer_("Write to file", &timers_)	
 {}
 
+pFlow::dictionary& pFlow::systemControl::domainDict()
+{
+	return 
+	settings().emplaceObjectOrGet<dictionary>
+	(
+		objectFile
+		(
+			domainFile__,
+			"",
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_NEVER
+		),
+		domainFile__,
+		true
+	);
+}
 
 bool pFlow::systemControl::operator ++(int)
 {
diff --git a/src/phasicFlow/repository/systemControl/systemControl.hpp b/src/phasicFlow/repository/systemControl/systemControl.hpp
index 4556ff40..8d6dfc3b 100644
--- a/src/phasicFlow/repository/systemControl/systemControl.hpp
+++ b/src/phasicFlow/repository/systemControl/systemControl.hpp
@@ -52,25 +52,25 @@ protected:
 
 
 
-	// - settings folder repository
+	/// settings folder repository
 	repository 			settings_;
 
-	// - caseSetup folder repository
+	/// caseSetup folder repository
 	repository 			caseSetup_;
 
-	// - settingsDict fileDictionary
+	/// settingsDict fileDictionary
 	dictionary& 		settingsDict_;
 
-	// - extra libs to be loaded
+	/// extra libs to be loaded
 	dynamicLinkLibs 	libs_;
 	
-	// - precision for writing to file 
+	/// precision for writing to file 
 	size_t 			outFilePrecision_ = 6;
 
-	// - time repository
+	/// time repository
 	Time 			Time_;
 
-	// - if time control is managed externaly
+	/// if time control is managed externaly
 
 	bool 				externalTimeControl_ = false;
 
@@ -163,6 +163,9 @@ public:
 		return settingsDict_;
 	}
 
+	dictionary& domainDict();
+	
+
 	virtual word runName() const
 	{
 		return runName_;
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp
index f2f1eeb9..353912d5 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp
@@ -1,10 +1,27 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
-
+-----------------------------------------------------------------------------*/
 
 #include "boundaryBase.hpp"
+#include "dictionary.hpp"
 
-
-pFlow::boundaryBase::boundaryBase
+/*pFlow::boundaryBase::boundaryBase
 (
 	const plane& 		bplane,
 	uint32 				mirrorProc, 
@@ -13,11 +30,60 @@ pFlow::boundaryBase::boundaryBase
 	internalPoints& 	internal
 )
 :
-boundaryPlane_(bplane),
-name_(name),
-type_(type),
-mirrorProcessoNo_(mirrorProc),
-internal_(internal)
+	subscriber(groupNames(name,type)),
+	boundaryPlane_(bplane),
+	name_(name),
+	type_(type),
+	mirrorProcessoNo_(mirrorProc),
+	internal_(internal)
 {
 	
+}*/
+
+pFlow::boundaryBase::boundaryBase
+(
+	const dictionary& 	dict,
+	const plane& 		bplane,
+	internalPoints& 	internal
+)
+:
+	subscriber(dict.name()),
+	boundaryPlane_(bplane),
+	neighborLength_(dict.getVal<real>("neighborLength")),
+	internal_(internal),
+	indexList_(groupNames(dict.name(),"indexList")),
+	mirrorProcessoNo_(dict.getVal<uint32>("mirrorProcessorNo")),
+	name_(dict.name()),
+	type_(dict.getVal<word>("type"))
+{
+
+}
+
+
+pFlow::uniquePtr<pFlow::boundaryBase> pFlow::boundaryBase::create
+(
+	const dictionary& 	dict,
+	const plane& 		bplane,
+	internalPoints& 	internal
+)
+{
+	word type = dict.getVal<word>("type");
+	word bType = angleBracketsNames("boundary", type);
+
+	if( dictionaryvCtorSelector_.search(bType) )
+	{
+		return dictionaryvCtorSelector_[bType] (dict, bplane, internal);
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< bType << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			dictionaryvCtorSelector_
+		);
+		fatalExit;
+	}
+	return nullptr;
 }
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp
index 709eb413..b66ebc47 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp
@@ -1,12 +1,32 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
 
 #ifndef __boundaryBase_hpp__
 #define __boundaryBase_hpp__
 
-
+#include "virtualConstructor.hpp"
+#include "subscriber.hpp"
 #include "VectorSingles.hpp"
 #include "plane.hpp"
 
-
+#include "streams.hpp"
 
 namespace pFlow
 {
@@ -14,8 +34,11 @@ namespace pFlow
 // forward 
 
 class internalPoints;
+class dictionary;
 
 class boundaryBase
+:
+	public subscriber
 {
 public:
 
@@ -26,45 +49,77 @@ public:
 
 protected:
 
-	plane 		boundaryPlane_;
+	const plane& 	boundaryPlane_;
+
+	/// The length defined for creating neighbor list 
+	real  			neighborLength_;	
+
+	/// a reference to 
+	internalPoints& internal_;
+
+	/// list of particles indices on device 
+	uint32Vector_D 	indexList_;
+
+	uint32 			mirrorProcessoNo_;
 
 	word 		name_;
 
 	word 		type_;
 
-	uint32 		mirrorProcessoNo_;
-
-	/// list of particles indices on device 
-	uint32Vector_D 	indexList_;
-
-	/// a reference to 
-	internalPoints& internal_;
 
 public:
 
 	TypeInfo("boundaryBase");
 
-	boundaryBase(
+	/*boundaryBase(
 		const plane& 		bplane,
 		uint32 				mirrorProc, 
 		const word& 		name,
 		const word& 		type,
+		internalPoints& 	internal);*/
+
+	boundaryBase(
+		const dictionary& 	dict,
+		const plane& 		bplane,
 		internalPoints& 	internal);
 
-	boundaryBase(const boundaryBase&) = default;
+	boundaryBase(const boundaryBase&) = delete;
 
-	boundaryBase& operator=(const boundaryBase&) = default;
+	boundaryBase& operator=(const boundaryBase&) = delete;
 
-	boundaryBase(boundaryBase&&) = delete;
+	boundaryBase(boundaryBase&&) = default;
 
-	boundaryBase& operator=(boundaryBase&&) = delete;
+	boundaryBase& operator=(boundaryBase&&) = default;
 
 	virtual ~boundaryBase() = default;
+	
 
-	virtual bool update() = 0 ;
+	create_vCtor
+	(
+		boundaryBase,
+		dictionary,
+		(
+			const dictionary& 	dict,
+			const plane& 		bplane,
+			internalPoints& 	internal
+		),
+		(dict, bplane, internal)
+	);
+
+	virtual bool beforeIteratoin(uint32 iterNum, real t) = 0 ;
 
 	virtual bool iterate(uint32 iterNum, real t) = 0;
 
+	virtual bool afterIteration(uint32 iterNum, real t) = 0;
+
+	static
+	uniquePtr<boundaryBase> create
+	(
+		const dictionary& 	dict,
+		const plane& 		bplane,
+		internalPoints& 	internal
+	);
+
 	
 };
 
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp b/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp
new file mode 100644
index 00000000..db5d4d5c
--- /dev/null
+++ b/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp
@@ -0,0 +1,64 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#include "boundaryExit.hpp"
+#include "dictionary.hpp"
+
+pFlow::boundaryExit::boundaryExit
+(
+	const dictionary& dict,
+	const plane&      bplane,
+	internalPoints&   internal
+)
+:
+	boundaryBase(dict, bplane, internal)
+{
+	exitMarginLength_ = max( 
+		dict.getValOrSet("exitMarginLength",0.0), 0.0);
+	checkForExitInterval_ = max(
+		dict.getValOrSet("checkForExitInterval", 1), 1);
+}
+
+bool pFlow::boundaryExit::beforeIteratoin
+(
+	uint32 iterNum, 
+	real t
+)
+{
+	return true;
+}
+
+bool pFlow::boundaryExit::iterate
+(
+	uint32 iterNum, 
+	real t
+)
+{
+	return true;
+}
+
+bool pFlow::boundaryExit::afterIteration
+(
+	uint32 iterNum, 
+	real t
+)
+{
+	return true;
+}
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.hpp b/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.hpp
new file mode 100644
index 00000000..0d07f48b
--- /dev/null
+++ b/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.hpp
@@ -0,0 +1,75 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __boundaryExit_hpp__
+#define __boundaryExit_hpp__
+
+
+#include "boundaryBase.hpp"
+
+namespace pFlow
+{
+
+class boundaryExit
+:
+ 	public boundaryBase
+{
+protected:
+
+	/// @brief Length of margin in the negative side of 
+	/// boundary to still allow points stay in the 
+	/// simulation. 
+	real 	exitMarginLength_ = 0;
+
+	/// @brief between successive check for particle exit
+	/// (iteration-based).
+	uint32 	checkForExitInterval_ = 1;
+
+public:
+
+	TypeInfo("boundary<exit>");
+
+	boundaryExit(
+		const dictionary& 	dict,
+		const plane&    	bplane,
+		internalPoints& 	internal);
+
+	virtual 
+	~boundaryExit() = default;
+	
+	add_vCtor
+	(
+		boundaryBase,
+		boundaryExit,
+		dictionary
+	);
+
+	bool beforeIteratoin(uint32 iterNum, real t) override;
+
+	bool iterate(uint32 iterNum, real t) override;
+
+	bool afterIteration(uint32 iterNum, real t) override;
+
+
+};
+
+}
+
+#endif
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/domain/domain.hpp b/src/phasicFlow/structuredData/domain/domain.hpp
index 220eda2f..15f57090 100644
--- a/src/phasicFlow/structuredData/domain/domain.hpp
+++ b/src/phasicFlow/structuredData/domain/domain.hpp
@@ -124,6 +124,17 @@ public:
 			return front_;
 		}
 
+		INLINE_FUNCTION_H
+		const auto& boundaryPlane(uint32 i)const 
+		{
+			if(i==0) return left_;
+			if(i==1) return right_;
+			if(i==2) return bottom_;
+			if(i==3) return top_;
+			if(i==4) return rear_;
+			return front_;
+		}
+
 	
 }; // domain
 
diff --git a/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp b/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp
new file mode 100644
index 00000000..74fd61da
--- /dev/null
+++ b/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp
@@ -0,0 +1,114 @@
+#include <cstring>
+
+#include "regularSimulationDomain.hpp"
+#include "processors.hpp"
+#include "streams.hpp"
+
+bool pFlow::regularSimulationDomain::createBoundaryDicts()
+{
+    auto& boundaries = globalDomainDict_.subDict("boundaries");
+    
+    globalDomainDict_.addDict("regularSimulationDomainBoundaries", boundaries);
+    auto& rbBoundaries = globalDomainDict_.subDict("regularSimulationDomainBoundaries");
+    
+	real neighborLength = boundaries.getVal<real>("neighborLength");
+
+	for(uint32 i=0; i<sizeOfBoundaries(); i++)
+	{
+		word bName = bundaryName(i);
+		if( !boundaries.containsDictionay(bName) )
+		{
+			fatalErrorInFunction<<"dictionary "<< bName<<
+			"does not exist in "<< boundaries.globalName()<<endl;
+			return false;
+		}
+		auto& bDict = rbBoundaries.subDict(bName);
+
+		if(!bDict.addOrKeep("neighborLength", neighborLength))
+		{
+			fatalErrorInFunction<<"error in adding neighborLength to "<< bName <<
+			"in dictionary "<< boundaries.globalName()<<endl;
+			return false;
+		}
+
+		if(!bDict.add("mirrorProcessorNo",(uint32) processors::globalRank()))
+		{
+			fatalErrorInFunction<<"error in adding mirrorProcessorNo to "<< bName <<
+			" in dictionary "<< boundaries.globalName()<<endl;
+			return false;
+		}
+
+	}
+    
+    
+    return true;
+}
+
+bool pFlow::regularSimulationDomain::setThisDomain()
+{
+	thisDomain_ = domain(globalBox_);
+    return true;
+}
+
+pFlow::regularSimulationDomain::regularSimulationDomain(
+    const dictionary &dict)
+    : simulationDomain(dict)
+{
+}
+
+bool pFlow::regularSimulationDomain::updateDomains
+(
+	const realx3Vector& pointPos
+)
+{
+	thisNumPoints_ = pointPos.size();
+	if(!createBoundaryDicts()) return false;
+	if(!setThisDomain()) return false;
+	
+    return true;
+}
+
+pFlow::uint32 pFlow::regularSimulationDomain::thisNumPoints() const
+{
+	return thisNumPoints_;	
+}
+
+bool pFlow::regularSimulationDomain::transferBlockData
+(
+	span<char> src, 
+	span<char> dst,
+	uint32 sizeOfBlock
+)
+{
+	size_t requiredSize = sizeOfBlock*thisNumPoints();
+	if(dst.size() < requiredSize)
+	{
+		fatalErrorInFunction<< 
+		"size of destination data block [in byte]"<< dst.size()<<
+		"is smaller than the required size [in byte]"<< requiredSize<<endl;
+		return false;
+	}
+	
+	if(src.size() < requiredSize )
+	{	
+		fatalErrorInFunction<< 
+		"size of source data block [in byte]"<< src.size()<<
+		"is smaller than the required size [in byte]"<< requiredSize<<endl;
+		return false;
+
+	}
+
+	std::memcpy(dst.data(), src.data(), requiredSize);
+
+    return true;
+}
+
+const pFlow::dictionary &pFlow::regularSimulationDomain::thisBoundaryDict() const
+{
+    return globalDomainDict_.subDict("regularSimulationDomainBoundaries");
+}
+
+bool pFlow::regularSimulationDomain::requiresDataTransfer()const
+{
+    return false;
+}
diff --git a/src/phasicFlow/structuredData/domain/regularSimulationDomain.hpp b/src/phasicFlow/structuredData/domain/regularSimulationDomain.hpp
new file mode 100644
index 00000000..011e8d2c
--- /dev/null
+++ b/src/phasicFlow/structuredData/domain/regularSimulationDomain.hpp
@@ -0,0 +1,56 @@
+
+#ifndef __regularSimulationDomain_hpp__
+#define __regularSimulationDomain_hpp__
+
+#include "simulationDomain.hpp"
+
+namespace pFlow
+{
+
+class regularSimulationDomain
+:
+	public simulationDomain
+{
+protected:
+
+	uint32 		thisNumPoints_=0;
+
+	bool createBoundaryDicts() override;
+
+	bool setThisDomain()override;
+
+public:
+
+	TypeInfo("simulationDomain<regular>");
+
+	regularSimulationDomain(const dictionary& dict);
+
+	virtual
+	~regularSimulationDomain()=default;
+	
+	add_vCtor
+	(
+		simulationDomain,
+		regularSimulationDomain,
+		dictionary
+	);
+
+	bool updateDomains( 
+		const realx3Vector& pointPos)override;
+
+	uint32 thisNumPoints()const override;
+
+	 
+	bool transferBlockData(
+		span<char> src, 
+		span<char> dst,
+		uint32 sizeOfBlock) override;
+	
+	const dictionary& thisBoundaryDict()const override;
+
+	bool requiresDataTransfer() const override;
+};
+
+}
+
+#endif //__regularSimulationDomain_hpp__
diff --git a/src/phasicFlow/structuredData/domain/simulationDomain.cpp b/src/phasicFlow/structuredData/domain/simulationDomain.cpp
new file mode 100644
index 00000000..1fb0867d
--- /dev/null
+++ b/src/phasicFlow/structuredData/domain/simulationDomain.cpp
@@ -0,0 +1,56 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#include "simulationDomain.hpp"
+#include "processors.hpp"
+
+pFlow::simulationDomain::simulationDomain(const dictionary& dict)
+:
+    globalBox_(dict.subDict("globalBox")),
+    globalDomainDict_(dict)
+{
+
+}
+
+pFlow::uniquePtr<pFlow::simulationDomain> 
+    pFlow::simulationDomain::create(const dictionary &dict)
+{
+    
+	word sType = angleBracketsNames(
+        "simulationDomain", 
+        processors::runTypeName());
+
+	if( dictionaryvCtorSelector_.search(sType) )
+	{
+		return dictionaryvCtorSelector_[sType] (dict);
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< sType << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			dictionaryvCtorSelector_
+		);
+		fatalExit;
+	}
+    return nullptr;
+}
diff --git a/src/phasicFlow/structuredData/domain/simulationDomain.hpp b/src/phasicFlow/structuredData/domain/simulationDomain.hpp
new file mode 100644
index 00000000..5844f437
--- /dev/null
+++ b/src/phasicFlow/structuredData/domain/simulationDomain.hpp
@@ -0,0 +1,142 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __simulationDomain_hpp__
+#define __simulationDomain_hpp__
+
+#include "domain.hpp"
+#include "virtualConstructor.hpp"
+#include "Vectors.hpp"
+
+#include "streams.hpp"
+
+namespace pFlow
+{
+
+class simulationDomain
+{
+private:
+
+static inline constexpr uint32 sizeOfBoundaries_ = 6;
+
+static 
+inline const std::array<word,6>  boundaryNames_ = 
+{
+	"left", "right", 
+	"top", "bottom",
+	"rear", "front"
+};
+
+protected:	
+
+	/// @brief acutal limits of the simulation domain
+	box 			globalBox_;
+
+	/// @brief the actual limits of this processor domain 
+	domain 			thisDomain_;
+
+	dictionary 		globalDomainDict_ ;
+
+	bool 			boundariesSet_ = false;
+	
+	virtual bool createBoundaryDicts() = 0;
+
+	virtual bool setThisDomain() = 0;
+
+public:
+
+	// - type info
+	TypeInfo("simulationDomain");
+
+	/// Constrcut from components 
+	simulationDomain(const dictionary& dict);
+	
+	/// Destructor 
+	virtual 
+	~simulationDomain() = default;
+	
+
+	create_vCtor
+	(
+		simulationDomain,
+		dictionary,
+		(const dictionary& dict),
+		(dict)	
+	);
+
+	virtual 
+	const dictionary& thisBoundaryDict()const = 0;
+	
+
+	virtual 
+	bool updateDomains(
+		const realx3Vector& pointPos) = 0;
+
+	//virtual getImportExportList() = 0;
+
+	//virtual getNumberInThisList() = 0;
+
+	virtual
+	uint32 thisNumPoints()const = 0;
+
+	virtual 
+	bool transferBlockData(
+		span<char> src, 
+		span<char> dst,
+		uint32 sizeOfBlock) = 0;
+
+	virtual 
+	bool requiresDataTransfer() const = 0;
+
+	inline
+	const dictionary& boundaryDict(uint32 i)const
+	{
+		return thisBoundaryDict().subDict(bundaryName(i));
+	}
+	
+	inline
+	const auto& boundaryPlane(uint32 i)const
+	{
+		return thisDomain_.boundaryPlane(i);
+	}
+
+	static inline
+	const word& bundaryName(uint32 i)
+	{
+		return boundaryNames_[i];
+	}
+
+	static inline 
+	uint32 sizeOfBoundaries() 
+	{
+		return boundaryNames_.size();
+	}
+	
+	static
+	uniquePtr<simulationDomain> create(const dictionary& dict);
+
+}; // simulationDomain
+
+
+
+}
+
+
+#endif //__simulationDomain_hpp__
diff --git a/src/phasicFlow/structuredData/domain/subDomain.hpp b/src/phasicFlow/structuredData/domain/subDomain.hpp
new file mode 100644
index 00000000..b90ccb61
--- /dev/null
+++ b/src/phasicFlow/structuredData/domain/subDomain.hpp
@@ -0,0 +1,104 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __subDomain_hpp__
+#define __subDomain_hpp__
+
+#include "domain.hpp"
+
+
+namespace pFlow
+{
+
+class subDomain
+{
+protected:	
+
+	/// the actual limits of global domain 
+	const domain& 		globalDomain_;
+
+	/// @brief dictionary of boundaries
+	const dictionary& 	boundaryDict_;
+
+
+	/// the actual limits of this processor domain 
+	domain 		thisDomain_;
+
+
+
+public:
+
+	// - type info
+	TypeInfo("subDomain");
+
+
+	subDomain(
+		const domain& gDomain,
+		const domain& thisD, 
+		const dictionary& bDicts)
+	:
+		globalDomain_(gDomain),
+		thisDomain_(thisD),
+		boundaryDict_(bDicts)
+	{}
+	
+	const auto& thisDomain()const
+	{
+		return thisDomain_;
+	}
+
+	const dictionary& leftDict()const
+	{
+		return boundaryDict_.subDict("left");
+	}
+
+	const dictionary& rightDict()const
+	{
+		return boundaryDict_.subDict("right");
+	}
+
+	const dictionary& bottomDict()const
+	{
+		return boundaryDict_.subDict("bottom");;
+	}
+
+	const dictionary& topDict()const
+	{
+		return boundaryDict_.subDict("top");;
+	}
+
+	const dictionary& rearDict()const
+	{
+		return boundaryDict_.subDict("rear");;
+	}
+
+	const dictionary& frontDict()const
+	{
+		return boundaryDict_.subDict("front");;
+	}
+	
+}; // subDomain
+
+
+
+}
+
+
+#endif //__subDomain_hpp__
diff --git a/src/phasicFlow/structuredData/pointStructure/boundaryList.cpp b/src/phasicFlow/structuredData/pointStructure/boundaryList.cpp
new file mode 100644
index 00000000..cd25a6b2
--- /dev/null
+++ b/src/phasicFlow/structuredData/pointStructure/boundaryList.cpp
@@ -0,0 +1,53 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#include "boundaryList.hpp"
+#include "simulationDomain.hpp"
+
+pFlow::boundaryList::boundaryList
+(
+	const simulationDomain& simD,
+	internalPoints& internal
+)
+:
+	simDomain_(simD),
+	internal_(internal),
+	boundaryList_(simD.sizeOfBoundaries())
+{}
+
+bool pFlow::boundaryList::updateLists()
+{
+
+	for(auto i=0; i<simDomain_.sizeOfBoundaries();i++)
+	{
+		boundaryList_.set
+		(
+			i,
+			boundaryBase::create
+			(
+				simDomain_.boundaryDict(i),
+				simDomain_.boundaryPlane(i),
+				internal_
+			)
+		);
+	}
+	listSet_ = true;
+	return true;
+}
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/pointStructure/boundaryList.hpp b/src/phasicFlow/structuredData/pointStructure/boundaryList.hpp
new file mode 100644
index 00000000..cb52807f
--- /dev/null
+++ b/src/phasicFlow/structuredData/pointStructure/boundaryList.hpp
@@ -0,0 +1,77 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+#ifndef __boundaryList_hpp__
+#define __boundaryList_hpp__
+
+#include "boundaryBase.hpp"
+#include "ListPtr.hpp"
+
+#include "streams.hpp"
+
+namespace pFlow
+{
+
+class simulationDomain;
+
+class boundaryList
+{
+
+protected:
+
+	//// - data members
+        internalPoints& 		internal_;
+		
+        const simulationDomain& simDomain_;
+
+        ListPtr<boundaryBase> 	boundaryList_;
+
+        bool                    listSet_ = false;
+
+    
+
+public:
+	
+	/// type info
+	TypeInfo("boundaryList");
+	
+
+	//// - Constructors
+
+	boundaryList(
+		const simulationDomain& simD, 
+		internalPoints& internal);
+
+	//~boundaryList() = default;
+
+  ~boundaryList() = default;
+
+  bool updateLists();
+};
+
+} // pFlow
+
+
+
+
+#endif //__boundaryList_hpp__
+
+
diff --git a/src/phasicFlow/structuredData/pointStructure/internalPoints.cpp b/src/phasicFlow/structuredData/pointStructure/internalPoints.cpp
index d551a7ea..d797993c 100644
--- a/src/phasicFlow/structuredData/pointStructure/internalPoints.cpp
+++ b/src/phasicFlow/structuredData/pointStructure/internalPoints.cpp
@@ -20,7 +20,7 @@ Licence:
 
 
 #include "internalPoints.hpp"
-
+#include "Vector.hpp"
 
 void pFlow::internalPoints::syncPFlag()const
 {
@@ -31,11 +31,6 @@ void pFlow::internalPoints::syncPFlag()const
 	}
 }
 
-/*#include "setFieldList.hpp"
-#include "error.hpp"
-#include "iOstream.hpp"
-
-#include "mortonIndexing.hpp"*/
 
 /*FUNCTION_H
 bool pFlow::internalPoints::evaluateinternalPoints()
@@ -153,6 +148,7 @@ pFlow::uniquePtr<pFlow::int32IndexContainer>
 
 pFlow::internalPoints::internalPoints()
 :
+	subscriber("internalPoints"),
 	pointPosition_("internalPoints", "internalPoints", initialCapacity_, 0, RESERVE()),
 	pFlagsD_(initialCapacity_, 0 , 0)
 {
@@ -162,9 +158,10 @@ pFlow::internalPoints::internalPoints()
 
 pFlow::internalPoints::internalPoints
 (
-	const std::vector<realx3>& posVec
+	const realx3Vector& posVec
 )
 :
+	subscriber("internalPoints"),
 	pointPosition_("internalPoints", "internalPoints", posVec.capacity(), 0, RESERVE()),
 	pFlagsD_(posVec.capacity(), 0, posVec.size())
 {
@@ -194,14 +191,19 @@ FUNCTION_H
 void pFlow::internalPoints::updateFlag
 (
 	const domain& dm, 
-	real dist
+	const std::array<real,6>& dist
 )
 {
-	pFlagsD_.markDeleteInDomain
+	pFlagsD_.markPointRegions
 	(
 		dm,
 		pointPosition_.deviceVectorAll(),
-		dist
+		dist[0],
+		dist[1],
+		dist[2],
+		dist[3],
+		dist[4],
+		dist[5]
 	);
 }
 
@@ -214,13 +216,20 @@ bool pFlow::internalPoints::read
 )
 {
 	
-	if( !pointPosition_.read(is, iotype))
+	Field<Vector, realx3 , vecAllocator<realx3>> fRead("internalPoints", "internalPoints");
+
+	if( !fRead.read(is, iotype))
 	{
 		fatalErrorInFunction<<
 		"Error in reading pointPosition from stream "<< is.name()<<endl;
 		return false;
 	}
 
+	/// here there should be some mechanism for field distribution between procesors
+	
+
+	pointPosition_.assign(fRead.vectorField());
+
 	pFlagsD_ = pFlagTypeDevice(pointPosition_.capacity(), 0, pointPosition_.size());
 	pFlagSync_ = false;
 	syncPFlag();
diff --git a/src/phasicFlow/structuredData/pointStructure/internalPoints.hpp b/src/phasicFlow/structuredData/pointStructure/internalPoints.hpp
index cec42337..095bdff1 100644
--- a/src/phasicFlow/structuredData/pointStructure/internalPoints.hpp
+++ b/src/phasicFlow/structuredData/pointStructure/internalPoints.hpp
@@ -17,29 +17,23 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
-
 #ifndef __internalPoints_hpp__
 #define __internalPoints_hpp__
 
+#include "subscriber.hpp"
+#include "Vectors.hpp"
 #include "Fields.hpp"
 #include "pointFlag.hpp"
 
 
 //#include "indexContainer.hpp"
 
-
-
 namespace pFlow
 {
 
-//forward
-/*class box;
-class setFieldList;
-class repository;*/
-
 class internalPoints
-
+:
+	public subscriber
 {
 public:
 	
@@ -47,11 +41,11 @@ public:
 
 	using pointsType 		= realx3Field_D;
 
-	using device_type 		= typename pointsType::deviceType;
+	using device_type 		= typename pointsType::device_type;
 
-	using memory_space 		= typename pointsType::memorySpace;
+	using memory_space 		= typename pointsType::memory_space;
 
-	using execution_space 	= typename pointsType::executionSpace;
+	using execution_space 	= typename pointsType::execution_space;
 
 	using pFlagTypeDevice 	= pointFlag<execution_space>;
 
@@ -62,22 +56,25 @@ protected:
 
 	//// - data members
 
-		// position of points on device
-		realx3Field_D 			pointPosition_;
+		/// Position of points on device
+		realx3Field_D 				pointPosition_;
 
-		// flag of points on device
-		mutable pFlagTypeDevice 		pFlagsD_;
+		/// flag of points on device
+		mutable pFlagTypeDevice 	pFlagsD_;
 
-		mutable pFlagTypeHost 			pFlagsH_;
+		/// flag of points on host
+		mutable pFlagTypeHost 		pFlagsH_;
 
-		mutable bool 					pFlagSync_ = false;
+		/// if both host and device flags sync
+		mutable bool 				pFlagSync_ = false;
 	
+	//// - protected members
 
-	void syncPFlag()const;
+		void syncPFlag()const;
 
 public:
 
-	friend class dynamicinternalPoints;
+	//friend class dynamicinternalPoints;
 	
 	// - type info
 	TypeInfo("internalPoints");
@@ -88,23 +85,21 @@ public:
 		// - an empty internalPoints, good for reading from file 
 		internalPoints();
 
-		/// Construct from components
-		//internalPoints(const int8Vector& flgVec, const realx3Vector& posVec);
-
+		
 		/// Construct from point positions, assume all points are active
-		internalPoints(const std::vector<realx3>& posVec);
+		internalPoints(const realx3Vector& posVec);
 
-		/// Copy construct
-		internalPoints(const internalPoints&) = default;
+		/// No Copy construct
+		internalPoints(const internalPoints&) = delete;
 
 		/// Move construct 
 		internalPoints(internalPoints&&) = default;
 
-		/// Copy assignment 
-		internalPoints& operator=(const internalPoints&) = default;
+		/// No Copy assignment 
+		internalPoints& operator=(const internalPoints&) = delete;
 
 		/// Move assignment 
-		internalPoints& operator=(internalPoints&&) = delete;
+		internalPoints& operator=(internalPoints&&) = default;
 
 		/// Destructor 
 		virtual ~internalPoints() = default;
@@ -116,7 +111,7 @@ public:
 
 
 		FUNCTION_H
-		const pFlagTypeHost& 	activePointsMaskH()const;
+		const pFlagTypeHost& activePointsMaskH()const;
 
 		// - Const access pointPosition
 		FUNCTION_H
@@ -166,7 +161,9 @@ public:
 
 		
 		FUNCTION_H
-		void updateFlag(const domain& dm, real dist);
+		void updateFlag(
+			const domain& dm, 
+			const std::array<real,6>& dist);
 
 		/*FUNCTION_H
 		size_t markDeleteOutOfBox(const box& domain);*/
@@ -227,6 +224,7 @@ public:
 
 };
 
+inline
 iOstream& operator<<(iOstream& os, const internalPoints& ip)
 {
 	if( !ip.write(os, IOPattern::AllProcessorsDifferent) )
diff --git a/src/phasicFlow/structuredData/pointStructure/pStcutold.txt b/src/phasicFlow/structuredData/pointStructure/pStcutold.txt
new file mode 100644
index 00000000..7d50ca4c
--- /dev/null
+++ b/src/phasicFlow/structuredData/pointStructure/pStcutold.txt
@@ -0,0 +1,540 @@
+#include "pointStructureKernels.hpp"
+#include "setFieldList.hpp"
+#include "error.hpp"
+#include "iOstream.hpp"
+//#include "Time.hpp"
+#include "mortonIndexing.hpp"
+
+FUNCTION_H
+bool pFlow::pointStructure::evaluatePointStructure()
+{
+	if(pointFlag_.size() != pointPosition_.size())
+	{
+		fatalErrorInFunction<<
+		"number of elements in pointFlag and pointPosition is not equal \n";
+		return false;
+	}
+
+	setNumMaxPoints();
+
+	int32 minActive, maxActive;
+	numActivePoints_ = pFlow::pointStructureKernels::scanPointFlag(
+		0,
+		numPoints_,
+		static_cast<int8>(pointStructure::ACTIVE),
+		pointFlag_.deviceVectorAll(),
+		minActive,
+		maxActive
+		);
+	
+	activeRange_ = {minActive, maxActive};	
+	
+	return true;
+}
+
+FUNCTION_H
+void pFlow::pointStructure::setNumMaxPoints()
+{
+	maxPoints_ = pointFlag_.capacity();
+	numPoints_ = pointFlag_.size();
+}
+
+FUNCTION_H
+pFlow::realx3Field_D& pFlow::pointStructure::pointPosition()
+{
+	return pointPosition_;
+}
+
+FUNCTION_H
+pFlow::int8Field_HD& pFlow::pointStructure::pointFlag()
+{
+	return pointFlag_;
+}
+		
+pFlow::uniquePtr<pFlow::int32IndexContainer> 
+	pFlow::pointStructure::getNewPointsIndices(int32 numNewPoints)const
+{
+
+
+	if( capacity() - activeRange_.second >=  numNewPoints )
+	{
+		// fill the sequence starting from activeRange_.second-1
+		return makeUnique<int32IndexContainer>(
+			activeRange_.second,
+			activeRange_.second+numNewPoints);
+
+	}
+
+	// second, check if there is space at the beginning
+	if( activeRange_.first >= numNewPoints)
+	{
+		return makeUnique<int32IndexContainer>(
+			0,
+			numNewPoints);	
+	}
+
+	// otherwise scan the points from first to the end to find empty spaces 
+	int32Vector newPoints(
+		numNewPoints,
+		RESERVE());
+
+	newPoints.clear();
+	int32 numAdded = 0;
+	ForAll(i, pointFlag_)
+	{
+		if(!isActive(i))
+		{
+			newPoints.push_back(static_cast<int32>(i));
+			numAdded++;
+		}
+
+		if(numAdded == numNewPoints)
+		{
+			return makeUnique<int32IndexContainer>(
+				newPoints.data(),
+				numNewPoints);
+		}
+	}
+
+	// check if there is space at the end for the remaining of points
+	if( numAdded <numNewPoints && capacity() - size() >= numNewPoints - numAdded )
+	{
+		int32 ind = size();
+		for(int32 i=numAdded; i<numNewPoints; i++)
+		{
+			newPoints.push_back(ind);
+			ind++;
+		}
+		
+		return makeUnique<int32IndexContainer>(
+			newPoints.data(),
+			numNewPoints);
+	}
+	else
+	{
+		fatalErrorInFunction<<"not enough capacity for inserting particles into the point structure\n";
+		return nullptr;
+	}
+
+	return nullptr;
+}
+
+pFlow::pointStructure::pointStructure()
+:
+	eventSubscriber(),	
+	pointFlag_("pointFlag", "pointFlag", maxPoints_, 0 , RESERVE()),
+	pointPosition_("pointPosition", "pointPosition", maxPoints_, 0, RESERVE()),
+	activeRange_(0,0)
+{}
+
+pFlow::pointStructure::pointStructure
+(
+	const int8Vector& flgVec,
+	const realx3Vector& posVec
+)
+:
+	eventSubscriber(),	
+	maxPoints_(posVec.capacity()),
+	pointFlag_("pointFlag", "pointFlag", maxPoints_, 0 , RESERVE()),
+	pointPosition_("pointPosition", "pointPosition", maxPoints_, 0, RESERVE()),
+	activeRange_(0,0)
+{
+	
+	pointFlag_.assign(flgVec);
+
+	pointPosition_.assign(posVec);
+
+	if( !evaluatePointStructure() )
+	{
+		fatalExit;
+	}
+}
+
+pFlow::pointStructure::pointStructure
+(
+	const realx3Vector& posVec
+)
+:
+	eventSubscriber(),	
+	maxPoints_(posVec.capacity()),
+	pointFlag_("pointFlag", "pointFlag", maxPoints_, 0 , RESERVE()),
+	pointPosition_("pointPosition", "pointPosition", maxPoints_, 0, RESERVE()),
+	activeRange_(0,0)
+{
+
+	pointPosition_.assign(posVec);
+
+	pointFlag_.resize(pointPosition_.size(), static_cast<int8>(pointStructure::ACTIVE) );
+	//pointFlag_.syncViews();
+
+	if( !evaluatePointStructure() )
+	{
+		fatalExit;
+	}
+}
+
+FUNCTION_H
+const pFlow::realx3Field_D& pFlow::pointStructure::pointPosition()const
+{
+	return pointPosition_;
+}
+
+FUNCTION_H
+const pFlow::int8Field_HD& pFlow::pointStructure::pointFlag()const
+{
+	return pointFlag_;
+}
+
+// - size of data structure
+FUNCTION_H
+pFlow::label pFlow::pointStructure::size()const
+{
+	return numPoints_;
+}
+
+FUNCTION_H
+pFlow::label pFlow::pointStructure::capacity()const
+{
+	return maxPoints_;
+}
+
+FUNCTION_H
+pFlow::label pFlow::pointStructure::numActive() const
+{
+	return numActivePoints_;
+}
+
+FUNCTION_H
+bool pFlow::pointStructure::allActive()const
+{
+	return numActivePoints_ == numPoints_;
+}
+
+FUNCTION_H
+bool pFlow::pointStructure::mortonSortPoints(const box& domain, real dx)
+{
+	if( !getSortedIndex(
+		domain, 
+		dx, 
+		activeRange_,
+		pointPosition_.deviceVectorAll(),
+		pointFlag_.deviceVectorAll(),
+		mortonSortedIndex_) )
+	{
+		fatalErrorInFunction<<"failed to perform morton sorting!"<<endl;
+		return false;
+	}
+
+	pointPosition_.sortItems(mortonSortedIndex_);
+	pointFlag_.sortItems(mortonSortedIndex_);
+
+	auto oldSize 			= size();
+	auto oldCapacity 	= capacity();
+	auto oldRange 		= activeRange();
+
+	// update size, range, capacity
+	setNumMaxPoints();
+	activeRange_ = {0, static_cast<int>(mortonSortedIndex_.size())};
+	numActivePoints_ = mortonSortedIndex_.size();
+
+	eventMessage msg(eventMessage::REARRANGE);
+
+	if(oldSize != size() )
+	{
+		msg.add(eventMessage::SIZE_CHANGED);
+	}
+
+	if(oldCapacity != capacity())
+	{
+		msg.add(eventMessage::CAP_CHANGED);
+	}
+
+	if( oldRange != activeRange_)
+	{
+		msg.add(eventMessage::RANGE_CHANGED);
+	}
+
+	// notify all the registered objects except the exclusionList 
+	if( !this->notify(msg) ) return false;
+
+	return true;
+}
+
+FUNCTION_H
+size_t pFlow::pointStructure::markDeleteOutOfBox(const box& domain)
+{
+	if(numPoints_==0) return 0;
+
+	int32 minRange, maxRange;
+	int32 numMarked = 
+	pFlow::pointStructureKernels::markDeleteOutOfBox(
+		domain,
+		activeRange_.first,
+		activeRange_.second,
+		pointStructure::DELETED,
+		pointPosition_.deviceVectorAll(),
+		pointFlag_.deviceVectorAll(),
+		activePointsMaskD(),
+		minRange, maxRange );
+
+	if(numMarked)
+	{
+		pointFlag_.modifyOnDevice();
+		pointFlag_.syncViews();	
+	}
+
+	if( numMarked<=numActivePoints_)
+	{
+		numActivePoints_ -= numMarked;
+	}
+	else
+	{
+		fatalErrorInFunction<<
+		"number of points marked as delete ("<<numMarked<<
+		") is greater than the activePoints ("<<numActivePoints_<<
+		").\n";
+		fatalExit;
+	}
+
+	range newRange = {minRange, maxRange};
+	
+	if( activeRange_ != newRange )
+	{
+		activeRange_ = newRange;
+		eventMessage msg(eventMessage::RANGE_CHANGED);
+		
+		// notify all the registered objects about active range change
+		if( !this->notify(msg) ) 
+		{
+			fatalExit<<
+			"something went wrong in notifying registered object about range change. \n";
+			fatalExit;
+		}
+	}	
+
+	return numMarked;
+	
+}
+
+FUNCTION_H
+bool pFlow::pointStructure::updateForDelete()
+{
+	notImplementedFunction;
+
+	return true;
+}
+
+/*FUNCTION_H
+pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::pointStructure::insertPoints
+(
+	const realx3Vector& pos,
+	const setFieldList& setField,
+	repository& owner,
+	const List<eventObserver*>& exclusionList
+)
+{
+
+	auto numNew = pos.size();
+	if( numNew==0)
+	{
+		return makeUnique<int32IndexContainer>();
+	}
+		
+	auto newPointsPtr = getNewPointsIndices( numNew );
+
+	if(!newPointsPtr)return nullptr;
+
+	auto oldSize 			= size();
+	auto oldCapacity 	= capacity();
+	auto oldRange 		= activeRange();
+
+	tobeInsertedIndex_ = newPointsPtr();
+
+	// set the position of new points
+	
+	if(!pointPosition_.insertSetElement(
+		newPointsPtr(),
+		pos)
+		)return nullptr;
+
+	
+	if(!pointFlag_.insertSetElement(
+		newPointsPtr(),
+		static_cast<int8>(PointFlag::ACTIVE))
+		)return nullptr;
+
+	
+	setNumMaxPoints();
+	auto minInd = newPointsPtr().min();
+	auto maxInd = newPointsPtr().max();
+
+
+	List<eventObserver*> localExlusion(exclusionList);
+
+	for(auto sfEntry:setField)
+	{
+
+		if(void* fieldPtr = 
+			sfEntry.setPointFieldSelectedAll( 
+				owner, 
+				newPointsPtr(), 
+				false );
+			fieldPtr)
+
+			localExlusion.push_back( 
+				static_cast<eventObserver*>(fieldPtr)
+				);
+		else
+			return nullptr;
+	}
+
+	// changes the active rage based on the new inserted points 
+	activeRange_ = { static_cast<int>(min(activeRange_.first, minInd )),
+					 static_cast<int>(max(activeRange_.second, maxInd+1))};
+
+	numActivePoints_ += numNew;
+
+	eventMessage msg(eventMessage::INSERT);
+
+	if( oldRange != activeRange_ )
+		msg.add(eventMessage::RANGE_CHANGED);
+
+	if( oldSize != size() )
+		msg.add(eventMessage::SIZE_CHANGED);
+
+	if( oldCapacity != capacity() )
+		msg.add(eventMessage::CAP_CHANGED);
+	
+	// notify all the registered objects except the exclusionList 
+	if( !this->notify(msg, localExlusion) ) return nullptr;
+
+	return newPointsPtr;
+}
+
+
+
+/*FUNCTION_H
+bool pFlow::pointStructure::readPointStructure
+(
+	iIstream& is
+)
+{
+	auto psCapacity = is.lookupDataOrSet("pointStructureCapacity", maxSizeDefault_);
+	pointPosition_.reallocate(psCapacity);
+	pointFlag_.reallocate(psCapacity);
+
+	if( !pointPosition_.read(is))
+	{
+		ioErrorInFile(is.name(), is.lineNumber())<<
+		"Error in reading pointPosition in pointStructure \n";
+		return false;	
+	}
+
+	if(! pointFlag_.read(is, true))
+	{
+		ioErrorInFile(is.name(), is.lineNumber())<<
+		"Error in reading pointFlag in pointStructure \n";
+		return false;
+	}
+
+	
+	return evaluatePointStructure();
+}*/
+
+/*FUNCTION_H
+bool pFlow::pointStructure::writePointStructure
+(
+	iOstream& os
+)const
+{
+	os.writeWordEntry("pointStructureCapacity", maxPoints_);
+
+	if(!pointPosition_.write(os))
+	{
+		ioErrorInFile(os.name(), os.lineNumber())<<
+		"error in writing pointPosition to file \n";
+		return false;
+	}
+
+	if(!pointFlag_.write(os))
+	{
+		ioErrorInFile(os.name(), os.lineNumber())<<
+		"error in writing pointFlag to file \n";
+		return false;
+	}
+	return true;
+}*/
+
+
+/*pFlow::uniquePtr<pFlow::int32Vector> 
+pFlow::pointStructure::newPointsIndices(
+	int32 numNewPoints
+)const
+{
+	
+	auto newPointsPtr = makeUnique<int32Vector>(
+		"pointStructure::newPointsPtr",
+		numNewPoints);
+
+	auto& newPoints = newPointsPtr();
+	
+
+	// first, check if there is space at the end
+	if( capacity() - activeRange_.second >=  numNewPoints )
+	{
+		// fill the sequence starting from activeRange_.second-1
+		fillSequence(newPoints, activeRange_.second-1);
+		return newPointsPtr;
+	}
+
+	// second, check if there is space at the beggining
+	if( activeRange_.first >= numNewPoints)
+	{
+		// fill the sequence starting from 0
+		fillSequence(newPoints, 0);
+		return newPointsPtr;
+	}
+
+	// otherwise scan the points from first to the end to find empty spaces 
+	newPoints.clear();
+	int32 numAdded = 0;
+	ForAll(i, pointFlag_)
+	{
+		if(!isActive(i))
+		{
+			newPoints.push_back(static_cast<int32>(i));
+			numAdded++;
+		}
+
+		if(numAdded == numNewPoints)
+		{
+			
+			return newPointsPtr;
+		}
+	}
+
+	// check if there is space at the end for the remaining of points
+	if( capacity() - size() >= numNewPoints - numAdded )
+	{
+		int32 ind = size();
+		for(int32 i=numAdded; i<numNewPoints; i++)
+		{
+			newPoints.push_back(ind);
+			ind++;
+		}
+		
+		return newPointsPtr;
+	}
+	else
+	{
+		fatalErrorInFunction<<"not enough capacity for inserting particles into the point structure\n";
+		return nullptr;
+	}
+
+	return nullptr;
+}*/
+
+pFlow::pointStructure::pointStructure(const dictionary &domainDict)
+{
+
+}
diff --git a/src/phasicFlow/structuredData/pointStructure/pointFlag.hpp b/src/phasicFlow/structuredData/pointStructure/pointFlag.hpp
index 974a2866..62861957 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointFlag.hpp
+++ b/src/phasicFlow/structuredData/pointStructure/pointFlag.hpp
@@ -24,12 +24,9 @@ Licence:
 #include "phasicFlowKokkos.hpp"
 #include "domain.hpp"
 
-
-
 namespace pFlow
 {
 
-
 template<typename ExecutionSpace>
 class pointFlag
 {
@@ -61,7 +58,7 @@ protected:
 
 	rangeU32 	activeRange_ 	= {0,0};
 
-	bool 			isAllActive_	= false;
+	bool 		isAllActive_	= false;
 
 	uint32 		nLeft_ 	= 0;
 
@@ -161,35 +158,19 @@ public:
 		return flg > INTERNAL;
 	}
 
-	INLINE_FUNCTION_HD
-	bool isLeft(uint32 i)const
-	{
-		return (flags_[i]&LEFT) == LEFT;
-	}
-
+	
 	INLINE_FUNCTION_HD
 	bool isLeft(uint8 flg)const
 	{
 		return (flg&LEFT) == LEFT;
 	}
 
-	INLINE_FUNCTION_HD
-	bool isRight(uint32 i)const
-	{
-		return (flags_[i]&&RIGHT) == RIGHT;
-	}
-
 	INLINE_FUNCTION_HD
 	bool isRight(uint8 flg)const
 	{
 		return (flg&&RIGHT) == RIGHT;
 	}
 
-	INLINE_FUNCTION_HD
-	bool isBottom(uint32 i)const
-	{
-		return (flags_[i]&&BOTTOM) == BOTTOM;
-	}
 
 	INLINE_FUNCTION_HD
 	bool isBottom(uint8 flg)const
@@ -197,43 +178,24 @@ public:
 		return (flg&&BOTTOM) == BOTTOM;
 	}
 
-	INLINE_FUNCTION_HD
-	bool isTop(uint32 i)const
-	{
-		return (flags_[i]&&TOP) == TOP;
-	}
-
 	INLINE_FUNCTION_HD
 	bool isTop(uint8 flg)const
 	{
 		return (flg&&TOP) == TOP;
 	}
 
-	INLINE_FUNCTION_HD
-	bool isRear(uint32 i)const
-	{
-		return (flags_[i]&&REAR) == REAR;
-	}
-
 	INLINE_FUNCTION_HD
 	bool isRear(uint8 flg)const
 	{
 		return (flg&&REAR) == REAR;
 	}
 
-	INLINE_FUNCTION_HD
-	bool isFront(uint32 i)const
-	{
-		return (flags_[i]&&FRONT) == FRONT;
-	}
-
 	INLINE_FUNCTION_HD
 	bool isFront(uint8 flg)const
 	{
 		return (flg&&FRONT) == FRONT;
 	}
 
-
 	template<typename ExeSpace>
 	pointFlag<ExeSpace> clone()const
 	{
@@ -254,10 +216,15 @@ public:
 	
 	uint32 scanPointFlag();
 
-	uint32 markDeleteInDomain(
+	uint32 markPointRegions(
 		domain 								dm,
 		ViewType1D<realx3, memory_space> 	points,
-		real dist);
+		real leftLength,
+		real rightLength,
+		real bottomLength,
+		real topLength,
+		real rearLength,
+		real frontLength);
 
 	void fillNeighborsLists(
 		ViewType1D<uint32, memory_space> leftList,
@@ -267,6 +234,12 @@ public:
 		ViewType1D<uint32, memory_space> rearList,
 		ViewType1D<uint32, memory_space> frontList);
 
+	uint32 markDelete(
+		const plane& 		pln,
+		ViewType1D<realx3, memory_space> 	points,
+		ViewType1D<uint32, memory_space> 	indices,
+		ViewType1D<uint32, memory_space> 	onOff);
+
 
 
 };
diff --git a/src/phasicFlow/structuredData/pointStructure/pointFlagKernels.hpp b/src/phasicFlow/structuredData/pointStructure/pointFlagKernels.hpp
index 107ac2e9..95fcd517 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointFlagKernels.hpp
+++ b/src/phasicFlow/structuredData/pointStructure/pointFlagKernels.hpp
@@ -81,11 +81,16 @@ pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::scanPointFlag()
 
 
 template<typename ExecutionSpace>
-pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markDeleteInDomain
+pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markPointRegions
 (
 	domain 								dm,
 	ViewType1D<realx3, memory_space> 	points,
-	real dist)
+	real leftLength,
+	real rightLength,
+	real bottomLength,
+	real topLength,
+	real rearLength,
+	real frontLength)
 {
 	
 	using rpMark = Kokkos::RangePolicy<execution_space,  
@@ -134,37 +139,37 @@ pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markDeleteInDomain
 						minUpdate = min(minUpdate,i);
 						maxUpdate = max(maxUpdate,i);
 
-						if(dm.left().inPositiveDistance(p, dist))
+						if(dm.left().inPositiveDistance(p, leftLength))
 						{
 							flg += LEFT;
 							valLeft++;	
 						} 
 
-						if(dm.right().inPositiveDistance(p, dist))
+						if(dm.right().inPositiveDistance(p, rightLength))
 						{
 							flg += RIGHT;
 							valRight++;
 						}
 						
-						if(dm.bottom().inPositiveDistance(p, dist))
+						if(dm.bottom().inPositiveDistance(p, bottomLength))
 						{
 							flg += BOTTOM;
 							valBottom++;
 						} 
 
-						if(dm.top().inPositiveDistance(p, dist)) 
+						if(dm.top().inPositiveDistance(p, topLength)) 
 						{
 							flg += TOP;
 							valTop++;
 						}
 						
-						if(dm.rear().inPositiveDistance(p, dist)) 
+						if(dm.rear().inPositiveDistance(p, rearLength)) 
 						{
 							flg += REAR;
 							valRear++;
 						}
 
-						if(dm.front().inPositiveDistance(p, dist))
+						if(dm.front().inPositiveDistance(p, frontLength))
 						{
 							flg += FRONT;
 							valFront++;
@@ -257,4 +262,76 @@ void pFlow::pointFlag<ExecutionSpace>::fillNeighborsLists
 
 }
 
+template<typename ExecutionSpace>
+pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markDelete
+(
+	const plane& 		pln,
+	ViewType1D<realx3, memory_space> 	points,
+	ViewType1D<uint32, memory_space> 	indices,
+	ViewType1D<uint32, memory_space> 	onOff
+)
+{
+	// if it is empty, do nothing
+	if(indices.size() == 0 )return 0;
+
+	using rpMark = Kokkos::RangePolicy<execution_space,  
+			Kokkos::IndexType<uint32>>;
+
+	uint32 start = activeRange().start();
+	uint32 end   = activeRange().end();
+
+	uint32 minRange = end;
+	uint32 maxRange = start;
+
+	uint32 numMarked = 0;
+
+	if(start<end)
+	{
+		Kokkos::parallel_reduce(
+			"pointFlagKernels::markDelete",
+			rpScanFlag(0, indices.size()),
+			CLASS_LAMBDA_HD(
+				uint32 i, 
+				uint32& minUpdate, 
+				uint32& maxUpdate, 
+				uint32& sumToUpdate)
+			{
+				auto indx = indices(i);
+				
+				if( pln.pointInNegativeSide(points(indx)))
+				{
+					flags_[indx] = DELETED;
+					sumToUpdate++;
+					onOff[i] = 0;
+				}
+				else
+				{
+					minUpdate = min(minUpdate,i);
+					maxUpdate = max(maxUpdate,i);
+					onOff[i] = 1;
+				}
+			},
+			Kokkos::Min<uint32>(minRange), 
+			Kokkos::Max<uint32>(maxRange),
+			numMarked);	
+	}
+
+	// means either range was empty or all points have been deleted.
+	if(minRange<start || maxRange>end)
+	{
+		minRange = 0;
+		maxRange = 0;
+	}
+	else
+	{
+		maxRange++; // add one to make it half
+	}
+
+	activeRange_ = {minRange, maxRange};
+	isAllActive_ = isAllActive_ && numMarked == 0; 
+	numActive_ 	-= numMarked;
+
+	return numMarked;
+}
+
 #endif // __pointFlagKernels_hpp__
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp b/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp
index 8dd2494e..a2d4b690 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp
@@ -18,540 +18,128 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-
 #include "pointStructure.hpp"
-#include "pointStructureKernels.hpp"
-#include "setFieldList.hpp"
-#include "error.hpp"
-#include "iOstream.hpp"
-//#include "Time.hpp"
-#include "mortonIndexing.hpp"
+#include "simulationDomain.hpp"
 
-FUNCTION_H
-bool pFlow::pointStructure::evaluatePointStructure()
-{
-	if(pointFlag_.size() != pointPosition_.size())
-	{
-		fatalErrorInFunction<<
-		"number of elements in pointFlag and pointPosition is not equal \n";
-		return false;
-	}
+#include "streams.hpp"
 
-	setNumMaxPoints();
-
-	int32 minActive, maxActive;
-	numActivePoints_ = pFlow::pointStructureKernels::scanPointFlag(
-		0,
-		numPoints_,
-		static_cast<int8>(pointStructure::ACTIVE),
-		pointFlag_.deviceVectorAll(),
-		minActive,
-		maxActive
-		);
-	
-	activeRange_ = {minActive, maxActive};	
-	
-	return true;
-}
-
-FUNCTION_H
-void pFlow::pointStructure::setNumMaxPoints()
-{
-	maxPoints_ = pointFlag_.capacity();
-	numPoints_ = pointFlag_.size();
-}
-
-FUNCTION_H
-pFlow::realx3Field_D& pFlow::pointStructure::pointPosition()
-{
-	return pointPosition_;
-}
-
-FUNCTION_H
-pFlow::int8Field_HD& pFlow::pointStructure::pointFlag()
-{
-	return pointFlag_;
-}
-		
-pFlow::uniquePtr<pFlow::int32IndexContainer> 
-	pFlow::pointStructure::getNewPointsIndices(int32 numNewPoints)const
-{
-
-
-	if( capacity() - activeRange_.second >=  numNewPoints )
-	{
-		// fill the sequence starting from activeRange_.second-1
-		return makeUnique<int32IndexContainer>(
-			activeRange_.second,
-			activeRange_.second+numNewPoints);
-
-	}
-
-	// second, check if there is space at the beginning
-	if( activeRange_.first >= numNewPoints)
-	{
-		return makeUnique<int32IndexContainer>(
-			0,
-			numNewPoints);	
-	}
-
-	// otherwise scan the points from first to the end to find empty spaces 
-	int32Vector newPoints(
-		numNewPoints,
-		RESERVE());
-
-	newPoints.clear();
-	int32 numAdded = 0;
-	ForAll(i, pointFlag_)
-	{
-		if(!isActive(i))
-		{
-			newPoints.push_back(static_cast<int32>(i));
-			numAdded++;
-		}
-
-		if(numAdded == numNewPoints)
-		{
-			return makeUnique<int32IndexContainer>(
-				newPoints.data(),
-				numNewPoints);
-		}
-	}
-
-	// check if there is space at the end for the remaining of points
-	if( numAdded <numNewPoints && capacity() - size() >= numNewPoints - numAdded )
-	{
-		int32 ind = size();
-		for(int32 i=numAdded; i<numNewPoints; i++)
-		{
-			newPoints.push_back(ind);
-			ind++;
-		}
-		
-		return makeUnique<int32IndexContainer>(
-			newPoints.data(),
-			numNewPoints);
-	}
-	else
-	{
-		fatalErrorInFunction<<"not enough capacity for inserting particles into the point structure\n";
-		return nullptr;
-	}
-
-	return nullptr;
-}
-
-pFlow::pointStructure::pointStructure()
+pFlow::pointStructure::pointStructure
+(
+  const dictionary& domainDict
+)
 :
-	eventSubscriber(),	
-	pointFlag_("pointFlag", "pointFlag", maxPoints_, 0 , RESERVE()),
-	pointPosition_("pointPosition", "pointPosition", maxPoints_, 0, RESERVE()),
-	activeRange_(0,0)
+  internalPoints(),
+  simulationDomain_
+  (
+    simulationDomain::create(domainDict)
+  ),
+  boundaries_
+  (
+    simulationDomain_(),
+    *this
+  )
 {}
 
-pFlow::pointStructure::pointStructure
+bool pFlow::pointStructure::distributePoints
 (
-	const int8Vector& flgVec,
-	const realx3Vector& posVec
+	const realx3Vector &points
 )
-:
-	eventSubscriber(),	
-	maxPoints_(posVec.capacity()),
-	pointFlag_("pointFlag", "pointFlag", maxPoints_, 0 , RESERVE()),
-	pointPosition_("pointPosition", "pointPosition", maxPoints_, 0, RESERVE()),
-	activeRange_(0,0)
 {
-	
-	pointFlag_.assign(flgVec);
+	uint32 thisN = simulationDomain_->thisNumPoints();
 
-	pointPosition_.assign(posVec);
+    Field<Vector, realx3 , vecAllocator<realx3>> internal(
+        "internalPoints", 
+        "internalPoints", 
+        thisN, 
+        thisN, 
+        RESERVE());
+    
+    
+    auto l = sizeof(realx3);
+    auto pointsSpan = span<char>(
+        reinterpret_cast<char*>
+		(
+			const_cast<realx3*>(points.data())
+		), 
+        points.size()*l );
+
+    auto internalSpan = span<char>(reinterpret_cast<char*>(internal.data()), internal.size()*l);
+
+    if(!simulationDomain_->transferBlockData(pointsSpan, internalSpan, l))
+    {
+        fatalErrorInFunction<<
+        "Error in transfering the block data "<<endl;
+        return false;
+    }
+    
+    pointPosition_.assign(internal.vectorField());
+
+    
+    pFlagsD_ = pFlagTypeDevice(pointPosition_.capacity(), 0, pointPosition_.size());
+    pFlagSync_ = false;
+    syncPFlag();
+    
+    
+
+    return true;
 
-	if( !evaluatePointStructure() )
-	{
-		fatalExit;
-	}
 }
 
-pFlow::pointStructure::pointStructure
+pFlow::pointStructure::pointStructure(
+    const dictionary &domainDict,
+    const realx3Vector &posVec)
+    : internalPoints(),
+      simulationDomain_(
+          simulationDomain::create(domainDict)),
+      boundaries_(
+          simulationDomain_(),
+          *this)
+{
+    if(!simulationDomain_->updateDomains(posVec) )
+    {
+        fatalErrorInFunction<<
+        "Failed to update domains."<<endl;
+        fatalExit;
+    }
+
+    boundaries_.updateLists();
+
+    if( !distributePoints(posVec) )
+    {
+        fatalErrorInFunction<<"Faild to distributes poinst"<<endl;
+        fatalExit;
+    }
+
+    
+}
+
+FUNCTION_H
+bool pFlow::pointStructure::read
 (
-	const realx3Vector& posVec
+    iIstream& is, 
+    IOPattern::IOType iotype
 )
-:
-	eventSubscriber(),	
-	maxPoints_(posVec.capacity()),
-	pointFlag_("pointFlag", "pointFlag", maxPoints_, 0 , RESERVE()),
-	pointPosition_("pointPosition", "pointPosition", maxPoints_, 0, RESERVE()),
-	activeRange_(0,0)
 {
+    Field<Vector, realx3 , vecAllocator<realx3>> fRead("file_internalPoints", "internalPoints");
+    
+    
 
-	pointPosition_.assign(posVec);
-
-	pointFlag_.resize(pointPosition_.size(), static_cast<int8>(pointStructure::ACTIVE) );
-	//pointFlag_.syncViews();
-
-	if( !evaluatePointStructure() )
-	{
-		fatalExit;
-	}
-}
-
-FUNCTION_H
-const pFlow::realx3Field_D& pFlow::pointStructure::pointPosition()const
-{
-	return pointPosition_;
-}
-
-FUNCTION_H
-const pFlow::int8Field_HD& pFlow::pointStructure::pointFlag()const
-{
-	return pointFlag_;
-}
-
-// - size of data structure
-FUNCTION_H
-pFlow::label pFlow::pointStructure::size()const
-{
-	return numPoints_;
-}
-
-FUNCTION_H
-pFlow::label pFlow::pointStructure::capacity()const
-{
-	return maxPoints_;
-}
-
-FUNCTION_H
-pFlow::label pFlow::pointStructure::numActive() const
-{
-	return numActivePoints_;
-}
-
-FUNCTION_H
-bool pFlow::pointStructure::allActive()const
-{
-	return numActivePoints_ == numPoints_;
-}
-
-FUNCTION_H
-bool pFlow::pointStructure::mortonSortPoints(const box& domain, real dx)
-{
-	if( !getSortedIndex(
-		domain, 
-		dx, 
-		activeRange_,
-		pointPosition_.deviceVectorAll(),
-		pointFlag_.deviceVectorAll(),
-		mortonSortedIndex_) )
-	{
-		fatalErrorInFunction<<"failed to perform morton sorting!"<<endl;
-		return false;
-	}
-
-	pointPosition_.sortItems(mortonSortedIndex_);
-	pointFlag_.sortItems(mortonSortedIndex_);
-
-	auto oldSize 			= size();
-	auto oldCapacity 	= capacity();
-	auto oldRange 		= activeRange();
-
-	// update size, range, capacity
-	setNumMaxPoints();
-	activeRange_ = {0, static_cast<int>(mortonSortedIndex_.size())};
-	numActivePoints_ = mortonSortedIndex_.size();
-
-	eventMessage msg(eventMessage::REARRANGE);
-
-	if(oldSize != size() )
-	{
-		msg.add(eventMessage::SIZE_CHANGED);
-	}
-
-	if(oldCapacity != capacity())
-	{
-		msg.add(eventMessage::CAP_CHANGED);
-	}
-
-	if( oldRange != activeRange_)
-	{
-		msg.add(eventMessage::RANGE_CHANGED);
-	}
-
-	// notify all the registered objects except the exclusionList 
-	if( !this->notify(msg) ) return false;
-
-	return true;
-}
-
-FUNCTION_H
-size_t pFlow::pointStructure::markDeleteOutOfBox(const box& domain)
-{
-	if(numPoints_==0) return 0;
-
-	int32 minRange, maxRange;
-	int32 numMarked = 
-	pFlow::pointStructureKernels::markDeleteOutOfBox(
-		domain,
-		activeRange_.first,
-		activeRange_.second,
-		pointStructure::DELETED,
-		pointPosition_.deviceVectorAll(),
-		pointFlag_.deviceVectorAll(),
-		activePointsMaskD(),
-		minRange, maxRange );
-
-	if(numMarked)
-	{
-		pointFlag_.modifyOnDevice();
-		pointFlag_.syncViews();	
-	}
-
-	if( numMarked<=numActivePoints_)
-	{
-		numActivePoints_ -= numMarked;
-	}
-	else
+	if( !fRead.read(is, iotype))
 	{
 		fatalErrorInFunction<<
-		"number of points marked as delete ("<<numMarked<<
-		") is greater than the activePoints ("<<numActivePoints_<<
-		").\n";
-		fatalExit;
-	}
-
-	range newRange = {minRange, maxRange};
-	
-	if( activeRange_ != newRange )
-	{
-		activeRange_ = newRange;
-		eventMessage msg(eventMessage::RANGE_CHANGED);
-		
-		// notify all the registered objects about active range change
-		if( !this->notify(msg) ) 
-		{
-			fatalExit<<
-			"something went wrong in notifying registered object about range change. \n";
-			fatalExit;
-		}
-	}	
-
-	return numMarked;
-	
-}
-
-FUNCTION_H
-bool pFlow::pointStructure::updateForDelete()
-{
-	notImplementedFunction;
-
-	return true;
-}
-
-FUNCTION_H
-pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::pointStructure::insertPoints
-(
-	const realx3Vector& pos,
-	const setFieldList& setField,
-	repository& owner,
-	const List<eventObserver*>& exclusionList
-)
-{
-
-	auto numNew = pos.size();
-	if( numNew==0)
-	{
-		return makeUnique<int32IndexContainer>();
-	}
-		
-	auto newPointsPtr = getNewPointsIndices( numNew );
-
-	if(!newPointsPtr)return nullptr;
-
-	auto oldSize 			= size();
-	auto oldCapacity 	= capacity();
-	auto oldRange 		= activeRange();
-
-	tobeInsertedIndex_ = newPointsPtr();
-
-	// set the position of new points
-	
-	if(!pointPosition_.insertSetElement(
-		newPointsPtr(),
-		pos)
-		)return nullptr;
-
-	
-	if(!pointFlag_.insertSetElement(
-		newPointsPtr(),
-		static_cast<int8>(PointFlag::ACTIVE))
-		)return nullptr;
-
-	
-	setNumMaxPoints();
-	auto minInd = newPointsPtr().min();
-	auto maxInd = newPointsPtr().max();
-
-
-	List<eventObserver*> localExlusion(exclusionList);
-
-	for(auto sfEntry:setField)
-	{
-
-		if(void* fieldPtr = 
-			sfEntry.setPointFieldSelectedAll( 
-				owner, 
-				newPointsPtr(), 
-				false );
-			fieldPtr)
-
-			localExlusion.push_back( 
-				static_cast<eventObserver*>(fieldPtr)
-				);
-		else
-			return nullptr;
-	}
-
-	// changes the active rage based on the new inserted points 
-	activeRange_ = { static_cast<int>(min(activeRange_.first, minInd )),
-					 static_cast<int>(max(activeRange_.second, maxInd+1))};
-
-	numActivePoints_ += numNew;
-
-	eventMessage msg(eventMessage::INSERT);
-
-	if( oldRange != activeRange_ )
-		msg.add(eventMessage::RANGE_CHANGED);
-
-	if( oldSize != size() )
-		msg.add(eventMessage::SIZE_CHANGED);
-
-	if( oldCapacity != capacity() )
-		msg.add(eventMessage::CAP_CHANGED);
-	
-	// notify all the registered objects except the exclusionList 
-	if( !this->notify(msg, localExlusion) ) return nullptr;
-
-	return newPointsPtr;
-}
-
-
-
-FUNCTION_H
-bool pFlow::pointStructure::readPointStructure
-(
-	iIstream& is
-)
-{
-	auto psCapacity = is.lookupDataOrSet("pointStructureCapacity", maxSizeDefault_);
-	pointPosition_.reallocate(psCapacity);
-	pointFlag_.reallocate(psCapacity);
-
-	if( !pointPosition_.read(is))
-	{
-		ioErrorInFile(is.name(), is.lineNumber())<<
-		"Error in reading pointPosition in pointStructure \n";
-		return false;	
-	}
-
-	if(! pointFlag_.read(is, true))
-	{
-		ioErrorInFile(is.name(), is.lineNumber())<<
-		"Error in reading pointFlag in pointStructure \n";
+		"Error in reading pointPosition from stream "<< is.name()<<endl;
 		return false;
 	}
 
-	
-	return evaluatePointStructure();
-}
+    if(!simulationDomain_->updateDomains(fRead))
+    {
+        fatalErrorInFunction<<
+        "error in updating domains"<<endl;
+        return false;
+    }
 
-FUNCTION_H
-bool pFlow::pointStructure::writePointStructure
-(
-	iOstream& os
-)const
-{
-	os.writeWordEntry("pointStructureCapacity", maxPoints_);
+    boundaries_.updateLists();
 
-	if(!pointPosition_.write(os))
-	{
-		ioErrorInFile(os.name(), os.lineNumber())<<
-		"error in writing pointPosition to file \n";
-		return false;
-	}
-
-	if(!pointFlag_.write(os))
-	{
-		ioErrorInFile(os.name(), os.lineNumber())<<
-		"error in writing pointFlag to file \n";
-		return false;
-	}
-	return true;
-}
-
-
-/*pFlow::uniquePtr<pFlow::int32Vector> 
-pFlow::pointStructure::newPointsIndices(
-	int32 numNewPoints
-)const
-{
-	
-	auto newPointsPtr = makeUnique<int32Vector>(
-		"pointStructure::newPointsPtr",
-		numNewPoints);
-
-	auto& newPoints = newPointsPtr();
-	
-
-	// first, check if there is space at the end
-	if( capacity() - activeRange_.second >=  numNewPoints )
-	{
-		// fill the sequence starting from activeRange_.second-1
-		fillSequence(newPoints, activeRange_.second-1);
-		return newPointsPtr;
-	}
-
-	// second, check if there is space at the beggining
-	if( activeRange_.first >= numNewPoints)
-	{
-		// fill the sequence starting from 0
-		fillSequence(newPoints, 0);
-		return newPointsPtr;
-	}
-
-	// otherwise scan the points from first to the end to find empty spaces 
-	newPoints.clear();
-	int32 numAdded = 0;
-	ForAll(i, pointFlag_)
-	{
-		if(!isActive(i))
-		{
-			newPoints.push_back(static_cast<int32>(i));
-			numAdded++;
-		}
-
-		if(numAdded == numNewPoints)
-		{
-			
-			return newPointsPtr;
-		}
-	}
-
-	// check if there is space at the end for the remaining of points
-	if( capacity() - size() >= numNewPoints - numAdded )
-	{
-		int32 ind = size();
-		for(int32 i=numAdded; i<numNewPoints; i++)
-		{
-			newPoints.push_back(ind);
-			ind++;
-		}
-		
-		return newPointsPtr;
-	}
-	else
-	{
-		fatalErrorInFunction<<"not enough capacity for inserting particles into the point structure\n";
-		return nullptr;
-	}
-
-	return nullptr;
-}*/
\ No newline at end of file
+	return distributePoints(fRead);
+    
+}
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp b/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp
index 09afa2af..e66a3b61 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp
@@ -22,46 +22,31 @@ Licence:
 #ifndef __pointStructure_hpp__
 #define __pointStructure_hpp__
 
-
-#include "Lists.hpp"
 #include "internalPoints.hpp"
-
-#include "Vectors.hpp"
-#include "VectorSingles.hpp"
-#include "VectorDuals.hpp"
-#include "Fields.hpp"
-#include "eventSubscriber.hpp"
-#include "indexContainer.hpp"
-
+#include "simulationDomain.hpp"
+#include "boundaryList.hpp"
+#include "uniquePtr.hpp"
 
 
 namespace pFlow
 {
 
-//forward
-class boundaryBase;
+class simulationDomain;
 
 
 class pointStructure
 :
-	public eventSubscriber
+	public internalPoints
 {
-public:
-	
-	using boundaryListType = ListPtr<boundaryBase>;
-
 protected:
 
 	//// - data members 
-		domain 				globalDomain_;
+	uniquePtr<simulationDomain> simulationDomain_ = nullptr;
 
-		domain 				thisDomain_;
-
-		internalPoints 		internal_;
-
-		boundaryListType 	boundaries_;
-
-		wordList 			boundaryTypes_;	
+	boundaryList 				boundaries_;
+		
+	
+	bool distributePoints(const realx3Vector& points);
 
 public:
 
@@ -73,21 +58,22 @@ public:
 	//// - Constructors 
 
 		// - an empty pointStructure, good for reading from file 
-		pointStructure();
+		pointStructure(const dictionary& domainDict);
 
 		// - construct from components
-		pointStructure(const int8Vector& flgVec, const realx3Vector& posVec);
+		//pointStructure(const int8Vector& flgVec, const realx3Vector& posVec);
 
-		// - construct from point positions, assume all points are active
-		pointStructure(const realx3Vector& posVec);
+		/// construct from point positions, assume all points are active
+		pointStructure(
+			const dictionary& domainDict, 
+			const realx3Vector& posVec);
 
 		// - copy construct
 		//
 		//   should be changed, may causs undefined behavior 
 		//////////////////////////////////////////////////////////////
-		pointStructure(const pointStructure&) = default;
+		pointStructure(const pointStructure&) = delete;
 
-		
 		// - no move construct 
 		pointStructure(pointStructure&&) = delete;
 
@@ -95,49 +81,25 @@ public:
 		//
 		//   should be changed, may causs undefined behavior 
 		//////////////////////////////////////////////////////////////
-		pointStructure& operator=(const pointStructure&) = default;
+		pointStructure& operator=(const pointStructure&) = delete;
 
 		// - no move assignment 
 		pointStructure& operator=(pointStructure&&) = delete;
 
 		// - destructor 
-		virtual ~pointStructure() = default
-
-		// - size of data structure
-		FUNCTION_H
-		label size()const;
+		virtual ~pointStructure() = default;
 		
-		// - maximum capacity of data structure 
+
+	/// Read 
 		FUNCTION_H
-		label capacity()const;
+		bool read(iIstream& is, IOPattern::IOType iotype);
+		
+
+		/// Write 
+		/*FUNCTION_H
+		bool write(iOstream& os, IOPattern::IOType iotype)const;	*/
 		
 		
-		// - update data structure by inserting/setting new points 
-		//   Notifies all the fields in the registered list of data structure
-		//   and exclude the fields that re in the exclusionList
-		//   retrun nullptr if it fails 
-		FUNCTION_H
-		/*virtual uniquePtr<int32IndexContainer> insertPoints(
-			const realx3Vector& pos,
-			const setFieldList& setField,
-			repository& owner,
-			const List<eventObserver*>& exclusionList={nullptr}
-		);*/
-
-		boundaryBase& boundary(size_t i)
-		{
-			return boundaries_[i];
-		}
-
-
-		const boundaryBase& boundary(size_t i)const
-		{
-			return boundaries_[i];
-		}
-
-
-	
-	//// - IO operations 
 		
 		
 
diff --git a/src/setHelpers/initialize.hpp b/src/setHelpers/initialize.hpp
index 8ee05c99..2bdd041b 100644
--- a/src/setHelpers/initialize.hpp
+++ b/src/setHelpers/initialize.hpp
@@ -23,8 +23,8 @@ Licence:
 
 // initilized and finalize should be placed in onc scope 
 REPORT(0)<<"Initializing host/device execution spaces . . . \n";
-REPORT(1)<<"Host execution space is "<< greenText(pFlow::DefaultHostExecutionSpace::name())<<endREPORT;
-REPORT(1)<<"Device execution space is "<<greenText(pFlow::DefaultExecutionSpace::name())<<endREPORT;
+REPORT(1)<<"Host execution space is "<< Green_Text(pFlow::DefaultHostExecutionSpace::name())<<END_REPORT;
+REPORT(1)<<"Device execution space is "<<Green_Text(pFlow::DefaultExecutionSpace::name())<<END_REPORT;
 
 Kokkos::initialize( argc, argv );
 {
diff --git a/src/setHelpers/initialize_Control.hpp b/src/setHelpers/initialize_Control.hpp
index b488a6b7..dfaa5ce7 100644
--- a/src/setHelpers/initialize_Control.hpp
+++ b/src/setHelpers/initialize_Control.hpp
@@ -24,10 +24,10 @@ Licence:
 // initilized and finalize should be placed in onc scope 
 #include "initialize.hpp"
 
-REPORT(0)<<"\nCreating Control repository . . ."<<endREPORT;
+REPORT(0)<<"\nCreating Control repository . . ."<<END_REPORT;
 pFlow::uniquePtr<pFlow::systemControl> ControlPtr = nullptr;
 
-if(isCoupling)
+/*if(isCoupling)
 {
 	pFlow::readControlDict controlDict;
 	
@@ -36,13 +36,13 @@ if(isCoupling)
 		controlDict.startTime(),
 		controlDict.endTime(),
 		controlDict.saveInterval(),
-		controlDict.startTimeName()
+		controlDict.startTimeNreseame()
 	);	
-}
-else
-{
+}*/
+/*else
+{*/
   	ControlPtr = pFlow::makeUnique<pFlow::systemControl>();
-}
+//}
 
 auto& Control = ControlPtr();
 
diff --git a/utilities/CMakeLists.txt b/utilities/CMakeLists.txt
index 3a7e766f..25d955d0 100644
--- a/utilities/CMakeLists.txt
+++ b/utilities/CMakeLists.txt
@@ -1,13 +1,13 @@
 
-add_subdirectory(checkPhasicFlow)
+#add_subdirectory(checkPhasicFlow)
 
-add_subdirectory(particlesPhasicFlow)
+#add_subdirectory(particlesPhasicFlow)
 
-add_subdirectory(geometryPhasicFlow)
+#add_subdirectory(geometryPhasicFlow)
 
-add_subdirectory(pFlowToVTK)
+#add_subdirectory(pFlowToVTK)
 
-add_subdirectory(Utilities)
+#add_subdirectory(Utilities)
 
-add_subdirectory(postprocessPhasicFlow)
+#add_subdirectory(postprocessPhasicFlow)
 
diff --git a/utilities/particlesPhasicFlow/empty/empty.cpp b/utilities/particlesPhasicFlow/empty/empty.cpp
index 1b938c00..71bd478a 100755
--- a/utilities/particlesPhasicFlow/empty/empty.cpp
+++ b/utilities/particlesPhasicFlow/empty/empty.cpp
@@ -23,10 +23,11 @@ Licence:
 
 
 pFlow::empty::empty(
+  systemControl& control,
 	const dictionary& dict
 )
 :
-	positionParticles(dict),
+	positionParticles(control, dict),
 	position_
 	(
 		maxNumberOfParticles_, 0, RESERVE()
diff --git a/utilities/particlesPhasicFlow/empty/empty.hpp b/utilities/particlesPhasicFlow/empty/empty.hpp
index fc3bada1..a4583b7c 100755
--- a/utilities/particlesPhasicFlow/empty/empty.hpp
+++ b/utilities/particlesPhasicFlow/empty/empty.hpp
@@ -44,7 +44,9 @@ public:
 	// - type Info
 	TypeInfo("empty");
 
-	empty(const dictionary& dict);
+	empty(
+		systemControl& control,
+		const dictionary& dict);
 
 	// - add this class to vCtor selection table 
 	add_vCtor(
@@ -56,12 +58,12 @@ public:
 
 	//// - Methods 
 
-	virtual label numPoints()const
+	virtual uint64 numPoints()const
 	{
 		return 0;
 	}
 
-	virtual label size()const
+	virtual uint64 size()const
 	{
 		return 0;
 	}
diff --git a/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp b/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp
index 683d33e4..173f09ce 100755
--- a/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp
+++ b/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp
@@ -24,7 +24,7 @@ Licence:
 #include "setFields.hpp"
 #include "systemControl.hpp"
 #include "commandLine.hpp"
-#include "readControlDict.hpp"
+//#include "readControlDict.hpp"
 
 using pFlow::output;
 using pFlow::endl;
diff --git a/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.cpp b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.cpp
index dc4d69e7..45a1635f 100755
--- a/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.cpp
+++ b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.cpp
@@ -123,6 +123,7 @@ bool pFlow::positionOrdered::positionPointsOrdered()
 
 pFlow::positionOrdered::positionOrdered
 (
+	systemControl& control,
 	const dictionary& dict
 )
 :
diff --git a/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.hpp b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.hpp
index 66452ab0..aa5a3881 100755
--- a/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.hpp
+++ b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.hpp
@@ -62,7 +62,9 @@ public:
 	// - type Info
 	TypeInfo("positionOrdered");
 
-	positionOrdered(const dictionary& dict);
+	positionOrdered(
+		systemControl& control,
+		const dictionary& dict);
 
 	// - add this class to vCtor selection table 
 	add_vCtor(
diff --git a/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp b/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp
index fcfb4616..632fceb7 100755
--- a/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp
+++ b/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp
@@ -76,6 +76,7 @@ pFlow::realx3Vector pFlow::positionParticles::sortByMortonCode(realx3Vector& pos
 
 pFlow::positionParticles::positionParticles
 (
+	systemControl& control,
 	const dictionary& dict
 )
 {
@@ -95,6 +96,11 @@ pFlow::positionParticles::positionParticles
 	{
 		region_ = makeUnique<region<sphere>>(dict.subDict("sphere"));
 	}
+	else
+	{
+		region_ = makeUnique<region<box>>( control.domainDict().subDict("globalBox"));
+	}
+	
 }
 
 
diff --git a/utilities/particlesPhasicFlow/positionParticles/positionParticles.hpp b/utilities/particlesPhasicFlow/positionParticles/positionParticles.hpp
index 48fc222b..37053a52 100755
--- a/utilities/particlesPhasicFlow/positionParticles/positionParticles.hpp
+++ b/utilities/particlesPhasicFlow/positionParticles/positionParticles.hpp
@@ -24,6 +24,7 @@ Licence:
 #include "virtualConstructor.hpp"
 #include "Vectors.hpp"
 #include "dictionary.hpp"
+#include "systemControl.hpp"
 
 namespace pFlow
 {
@@ -118,7 +119,9 @@ public:
 	// - type Info
 	TypeInfo("positionParticles");
 
-	positionParticles(const dictionary& dict);
+	positionParticles(
+		systemControl& control, 
+		const dictionary& dict);
 
 	create_vCtor
 	(
@@ -132,9 +135,9 @@ public:
 
 	//// - Methods 
 
-	virtual label numPoints()const = 0;
+	virtual uint64 numPoints()const = 0;
 
-	virtual label size()const = 0;
+	virtual uint64 size()const = 0;
 
 	virtual real maxDiameter() const = 0;
 
diff --git a/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp b/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp
index 616dbbca..d0053791 100755
--- a/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp
+++ b/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp
@@ -158,10 +158,11 @@ bool pFlow::positionRandom::inCollision
 
 pFlow::positionRandom::positionRandom
 (
+	systemControl& control,
 	const dictionary& dict
 )
 :
-	positionParticles(dict),
+	positionParticles(control, dict),
 	prDict_
 	(  
 		dict.subDict("positionRandomInfo")
diff --git a/utilities/particlesPhasicFlow/positionRandom/positionRandom.hpp b/utilities/particlesPhasicFlow/positionRandom/positionRandom.hpp
index 33bc9899..ff301c54 100755
--- a/utilities/particlesPhasicFlow/positionRandom/positionRandom.hpp
+++ b/utilities/particlesPhasicFlow/positionRandom/positionRandom.hpp
@@ -66,7 +66,9 @@ public:
 	// - type Info
 	TypeInfo("positionRandom");
 
-	positionRandom(const dictionary& dict);
+	positionRandom(
+		systemControl& control,
+		const dictionary& dict);
 
 	// - add this class to vCtor selection table 
 	add_vCtor(
@@ -78,12 +80,12 @@ public:
 
 	//// - Methods 
 
-	virtual label numPoints()const
+	virtual uint64 numPoints()const
 	{
 		return position_.size();
 	}
 
-	virtual label size()const
+	virtual uint64 size()const
 	{
 		return position_.size();
 	}