From 614b2f732e18fb1e446e79393f185cd659c87d1a Mon Sep 17 00:00:00 2001
From: HRN <hamid.r.norouzi@gmail.com>
Date: Sun, 12 May 2024 19:10:04 +0330
Subject: [PATCH 1/2] develop branch update for changes in MPI branch for data
 transfer stage mortong indexing is added (messaging is not complete)

---
 .../contactLists/unsortedContactList.hpp      |   5 +-
 .../ContactSearch/ContactSearch.hpp           |  20 ++-
 .../boundaryContactSearch.hpp                 |   5 -
 .../contactSearch/contactSearch.hpp           |   9 ++
 .../methods/cellBased/NBS/NBS.hpp             |  10 ++
 .../methods/particleWallContactSearchs.hpp    |  15 +++
 src/Interaction/interaction/interaction.cpp   |   2 +-
 .../boundarySphereInteraction.cpp             |  17 ++-
 .../boundarySphereInteraction.hpp             |  38 ++++--
 .../periodicBoundarySphereInteraction.cpp     |  12 +-
 .../periodicBoundarySphereInteraction.hpp     |   3 +-
 .../sphereInteraction/sphereInteraction.cpp   | 107 +++++++++++++---
 .../sphereInteraction/sphereInteraction.hpp   |   3 +
 .../sphereParticles/sphereParticles.cpp       |   7 +-
 .../sphereParticles/sphereParticles.hpp       |   2 +-
 .../dynamicPointStructure.cpp                 |  11 +-
 .../dynamicPointStructure.hpp                 |   2 +
 src/Particles/particles/particles.cpp         |  11 +-
 src/Particles/particles/particles.hpp         |   2 +
 src/phasicFlow/CMakeLists.txt                 |   2 +
 .../containers/VectorHD/VectorSingle.cpp      |  88 ++++++++++++-
 .../containers/VectorHD/VectorSingle.hpp      |   8 +-
 .../indexContainer/indexContainer.hpp         |  20 ++-
 .../boundaryField/boundaryField.hpp           |  18 ++-
 .../generalBoundary/generalBoundary.hpp       |   4 -
 .../internalField/internalField.cpp           |  32 ++++-
 .../internalField/internalField.hpp           |   2 +
 src/phasicFlow/eventManagement/message.hpp    |  12 +-
 .../repository/Time/baseTimeControl.cpp       |  10 ++
 .../repository/Time/baseTimeControl.hpp       |   7 +
 src/phasicFlow/streams/dataIO/dataIO.cpp      |   4 +-
 .../boundaries/boundaryBase/boundaryBase.cpp  |   5 +-
 .../boundaries/boundaryBase/boundaryBase.hpp  |  39 +++++-
 .../boundaryBase/scatteredFieldAccess.hpp     |  10 ++
 .../boundaries/boundaryExit/boundaryExit.cpp  |   3 +
 .../boundaries/boundaryList.cpp               |  36 +++---
 .../boundaries/boundaryList.hpp               |  14 +-
 .../boundaries/boundaryNone/boundaryNone.hpp  |   8 +-
 .../boundaryPeriodic/boundaryPeriodic.cpp     |  15 ++-
 .../boundaryPeriodic/boundaryPeriodic.hpp     |   3 +-
 src/phasicFlow/structuredData/cells/cells.hpp | 107 +++++-----------
 .../domain/regularSimulationDomain.cpp        |  13 +-
 .../internalPoints/internalPoints.cpp         | 121 +++++++++++++++++-
 .../internalPoints/internalPoints.hpp         |   4 +
 .../internalPoints/pointFlag.hpp              |  27 +++-
 .../internalPoints/pointFlagKernels.hpp       |  45 +++++++
 .../{ => pointSorting}/mortonIndexing.cpp     |  60 ++++-----
 .../{ => pointSorting}/mortonIndexing.hpp     |  12 +-
 .../pointSorting/pointSorting.cpp             |  52 ++++++++
 .../pointSorting/pointSorting.hpp             |  68 ++++++++++
 .../pointStructure/pointStructure.cpp         |  64 ++++++++-
 .../pointStructure/pointStructure.hpp         |  11 +-
 52 files changed, 937 insertions(+), 268 deletions(-)
 rename src/phasicFlow/structuredData/pointStructure/pointStructure/{ => pointSorting}/mortonIndexing.cpp (59%)
 rename src/phasicFlow/structuredData/pointStructure/pointStructure/{ => pointSorting}/mortonIndexing.hpp (94%)
 create mode 100644 src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp
 create mode 100644 src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.hpp

diff --git a/src/Interaction/contactLists/unsortedContactList.hpp b/src/Interaction/contactLists/unsortedContactList.hpp
index 5922c4dc..8a61b2ce 100644
--- a/src/Interaction/contactLists/unsortedContactList.hpp
+++ b/src/Interaction/contactLists/unsortedContactList.hpp
@@ -95,7 +95,8 @@ public:
 		// swap conainer and values 
 		swapViews(values0_, values_);
 		swapViews(container0_, this->container_);
-		return UnsortedPairs::beforeBroadSearch();
+		UnsortedPairs::beforeBroadSearch();
+		return true;
 	}
 
 	bool afterBroadSearch()
@@ -110,7 +111,7 @@ public:
 			rpFillPairs(0,this->capacity()),
 			*this);
 		Kokkos::fence();
-
+		
 		return true; 
 	}
 
diff --git a/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp b/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp
index d80813ca..f2dda414 100644
--- a/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp
+++ b/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp
@@ -160,7 +160,8 @@ public:
 		csPairContainerType& pwPairs,
 		bool force = false)override
 	{
-		Particles().boundingSphere().updateBoundaries(DataDirection::SlaveToMaster);
+		if(i==0u)
+			Particles().boundingSphere().updateBoundaries(DataDirection::SlaveToMaster);
 		return csBoundaries_[i].broadSearch(
 			iter, 
 			t, 
@@ -183,6 +184,23 @@ public:
 	{
 		return ppwContactSearch_().performedSearch();
 	}
+
+	
+	uint32 updateInterval()const override
+	{
+		return ppwContactSearch_().updateInterval();
+	}
+
+	real sizeRatio()const override
+	{
+		return ppwContactSearch_().sizeRatio();
+	}
+
+	 
+	real cellExtent()const override
+	{
+		return ppwContactSearch_().cellExtent();
+	}
 	
 };
 
diff --git a/src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.hpp b/src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.hpp
index 78c61bb8..55ecb47e 100644
--- a/src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.hpp
+++ b/src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.hpp
@@ -72,11 +72,6 @@ public:
         return contactSearch_;
     }
 
-    void fill(const std::any &val) override
-    {
-        return;
-    }
-
     bool hearChanges(
         real t,
         real dt,
diff --git a/src/Interaction/contactSearch/contactSearch/contactSearch.hpp b/src/Interaction/contactSearch/contactSearch/contactSearch.hpp
index 235403c4..865ec291 100644
--- a/src/Interaction/contactSearch/contactSearch/contactSearch.hpp
+++ b/src/Interaction/contactSearch/contactSearch/contactSearch.hpp
@@ -140,6 +140,15 @@ public:
 	virtual 
 	bool performedBroadSearch()const = 0;
 
+	virtual
+	uint32 updateInterval()const = 0;
+
+	virtual 
+	real sizeRatio()const = 0;
+
+	virtual 
+	real cellExtent()const = 0;
+	
 	static 
 	uniquePtr<contactSearch> create(
 		const dictionary& dict,
diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp b/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp
index 11d2e108..9cd5ab6f 100644
--- a/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp
+++ b/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp
@@ -134,6 +134,16 @@ public:
 		{
 			return 1;
 		}
+		
+		real sizeRatio()const
+		{
+			return sizeRatio_;
+		}
+
+		real cellExtent()const
+		{
+			return cellExtent_;
+		}
 
 		auto getCellIterator([[maybe_unused]] uint32 lvl)const
 		{
diff --git a/src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp b/src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp
index d72f9404..1bf5dab3 100644
--- a/src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp
+++ b/src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp
@@ -113,6 +113,21 @@ public:
 		return false;		
 	}
 
+	uint32 updateInterval()const
+	{
+		return updateInterval_;
+	}
+
+	real sizeRatio()const
+	{
+		return getMethod().sizeRatio();
+	}
+
+	real cellExtent()const
+	{
+		return getMethod().cellExtent();
+	}
+
 };
 
 
diff --git a/src/Interaction/interaction/interaction.cpp b/src/Interaction/interaction/interaction.cpp
index c6099405..03abbb4b 100644
--- a/src/Interaction/interaction/interaction.cpp
+++ b/src/Interaction/interaction/interaction.cpp
@@ -72,7 +72,7 @@ pFlow::uniquePtr<pFlow::interaction> pFlow::interaction::create
 
 
 	gSettings::sleepMiliSeconds(
-			1000*(pFlowProcessors().localSize()-pFlowProcessors().localRank()-1));
+			100*(pFlowProcessors().localSize()-pFlowProcessors().localRank()-1));
 	pOutput.space(2)<<"Creating interaction "<<Green_Text(interactionModel)<<" . . ."<<END_REPORT;
 	if( systemControlvCtorSelector_.search(interactionModel) )
 	{
diff --git a/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.cpp b/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.cpp
index f9f6374a..74a50bd3 100644
--- a/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.cpp
+++ b/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.cpp
@@ -1,3 +1,4 @@
+#include "boundarySphereInteraction.hpp"
 /*------------------------------- phasicFlow ---------------------------------
       O        C enter of
      O O       E ngineering and
@@ -18,6 +19,20 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+template <typename cFM, typename gMM>
+void pFlow::boundarySphereInteraction<cFM, gMM>::allocatePPPairs()
+{
+	ppPairs_.reset(nullptr);
+	ppPairs_ = makeUnique<ContactListType>(1);
+}
+
+template <typename cFM, typename gMM>
+void pFlow::boundarySphereInteraction<cFM, gMM>::allocatePWPairs()
+{
+	pwPairs_.reset(nullptr);
+	pwPairs_ = makeUnique<ContactListType>(1);
+}
+
 
 template <typename cFM, typename gMM>
 pFlow::boundarySphereInteraction<cFM, gMM>::boundarySphereInteraction(
@@ -28,8 +43,6 @@ pFlow::boundarySphereInteraction<cFM, gMM>::boundarySphereInteraction(
 	  geometryMotion_(geomMotion),
 	  sphParticles_(sphPrtcls)
 {
-	ppPairs_ = makeUnique<ContactListType>(1);
-	pwPairs_ = makeUnique<ContactListType>(1);
 }
 
 template <typename cFM, typename gMM>
diff --git a/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.hpp b/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.hpp
index 56058b6d..708fb61e 100644
--- a/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.hpp
+++ b/src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.hpp
@@ -22,7 +22,7 @@ Licence:
 
 #include "virtualConstructor.hpp"
 #include "generalBoundary.hpp"
-#include "unsortedContactList.hpp"
+#include "sortedContactList.hpp"
 #include "sphereParticles.hpp"
 
 namespace pFlow
@@ -51,7 +51,7 @@ public:
 	using IndexType    		    = uint32;
     
     using ContactListType = 
-        unsortedContactList<ModelStorage, DefaultExecutionSpace, IdType>;
+        sortedContactList<ModelStorage, DefaultExecutionSpace, IdType>;
 
 private:
 
@@ -60,9 +60,15 @@ private:
 	/// const reference to sphere particles 
 	const sphereParticles& 				sphParticles_;
     
-	uniquePtr<ContactListType> 			ppPairs_;
+	uniquePtr<ContactListType> 			ppPairs_ = nullptr;
 
-	uniquePtr<ContactListType>			pwPairs_;
+	uniquePtr<ContactListType>			pwPairs_ = nullptr;
+
+protected:
+
+	void allocatePPPairs();
+
+	void allocatePWPairs();
 
 public:
 
@@ -124,15 +130,30 @@ public:
 		return pwPairs_();
 	}
 
+	bool ppPairsAllocated()const
+	{
+		if( ppPairs_)return true;
+		return false;
+	}
+
+	bool pwPairsAllocated()const
+	{
+		if( pwPairs_)return true;
+		return false;
+	}
+
 	virtual
 	bool sphereSphereInteraction(
 		real dt,
-		const ContactForceModel& cfModel)
+		const ContactForceModel& cfModel,
+		uint32 step)
 	{
 		// for default boundary, no thing to be done 
-		return true;
+		return false;
 	}
 
+
+
 	bool hearChanges
 	(
 		real t,
@@ -149,11 +170,6 @@ public:
 		return true;
 	}
 
-	void fill(const std::any& val)override
-	{
-		notImplementedFunction;
-	}
-
 	static
 	uniquePtr<BoundarySphereInteractionType> create(
 		const boundaryBase& boundary,
diff --git a/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp b/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp
index 38a4336b..3200c4e7 100644
--- a/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp
+++ b/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp
@@ -30,7 +30,10 @@ pFlow::periodicBoundarySphereInteraction<cFM, gMM>::periodicBoundarySphereIntera
 {
 	if(boundary.thisBoundaryIndex()%2==1)
     {
-       masterInteraction_ = true;
+    	masterInteraction_ = true;
+		this->allocatePPPairs();
+		this->allocatePWPairs();
+	   
     }
     else
     {
@@ -42,10 +45,11 @@ template <typename cFM, typename gMM>
 bool pFlow::periodicBoundarySphereInteraction<cFM, gMM>::sphereSphereInteraction
 (
 	real dt,
-	const ContactForceModel &cfModel
+	const ContactForceModel &cfModel,
+	uint32 step
 )
 {
-	if(!masterInteraction_) return true;
+	if(!masterInteraction_) return false;
 	
 	pFlow::periodicBoundarySIKernels::sphereSphereInteraction(
 		dt,
@@ -61,5 +65,5 @@ bool pFlow::periodicBoundarySphereInteraction<cFM, gMM>::sphereSphereInteraction
 		this->sphParticles().contactForce().deviceViewAll(),
 		this->sphParticles().contactTorque().deviceViewAll());
 
-	return true;
+	return false;
 }
\ No newline at end of file
diff --git a/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.hpp b/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.hpp
index e2887787..5f1e7a45 100644
--- a/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.hpp
+++ b/src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.hpp
@@ -83,7 +83,8 @@ public:
 
     bool sphereSphereInteraction(
         real dt,
-		const ContactForceModel& cfModel)override;
+		const ContactForceModel& cfModel,
+        uint32 step)override;
 	
 };
 
diff --git a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp
index b4546521..2aa10ef8 100644
--- a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp
+++ b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp
@@ -137,7 +137,9 @@ pFlow::sphereInteraction<cFM,gMM, cLT>::sphereInteraction
 	ppInteractionTimer_("sphere-sphere interaction", &this->timers()),
 	pwInteractionTimer_("sphere-wall interaction", &this->timers()),
 	contactListMangementTimer_("contact-list management", &this->timers()),
-	boundaryInteractionTimer_("interaction for boundary", &this->timers())
+	boundaryContactSearchTimer_("contact search for boundary", &this->timers()),
+	boundaryInteractionTimer_("interaction for boundary", &this->timers()),
+	contactListBoundaryTimer_("contact-list management for boundary", &this->timers())
 {
 	
 	if(!createSphereInteraction())
@@ -179,11 +181,14 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 		contactListMangementTimer_.pause();
 	} 
 
+	contactListBoundaryTimer_.start();
 	for(uint32 i=0; i<6u; i++)
 	{
-		boundaryInteraction_[i].ppPairs().beforeBroadSearch();
-		boundaryInteraction_[i].pwPairs().beforeBroadSearch();
+		auto& BI = boundaryInteraction_[i];
+		if(BI.ppPairsAllocated()) BI.ppPairs().beforeBroadSearch();
+		if(BI.pwPairsAllocated()) BI.pwPairs().beforeBroadSearch();
 	}
+	contactListBoundaryTimer_.pause();
 	
 	if( sphParticles_.numActive()<=0)return true;
 	
@@ -200,21 +205,27 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 		fatalExit;
 	}
 	
+	boundaryContactSearchTimer_.start();
 	for(uint32 i=0; i<6u; i++)
 	{
-		if( !contactSearch_().boundaryBroadSearch(
-			i,
-			iter,
-			t,
-			dt,
-			boundaryInteraction_[i].ppPairs(),
-			boundaryInteraction_[i].pwPairs()))
+		auto& BI =boundaryInteraction_[i];
+		if(BI.ppPairsAllocated())
 		{
-			fatalErrorInFunction<<
-			"failed to perform broadSearch for boundary index "<<i<<endl;
-			return false;
-		}
+			if( !contactSearch_().boundaryBroadSearch(
+				i,
+				iter,
+				t,
+				dt,
+				BI.ppPairs(),
+				BI.pwPairs()))
+			{
+				fatalErrorInFunction<<
+				"failed to perform broadSearch for boundary index "<<i<<endl;
+				return false;
+			}
+		} 
 	}
+	boundaryContactSearchTimer_.end();
 	
 	if(broadSearch && contactSearch_().performedBroadSearch())
 	{
@@ -224,12 +235,44 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 		contactListMangementTimer_.end();
 	}
 
+	contactListBoundaryTimer_.resume();
 	for(uint32 i=0; i<6u; i++)
 	{
-		boundaryInteraction_[i].ppPairs().afterBroadSearch();
-		boundaryInteraction_[i].pwPairs().afterBroadSearch();
+		auto& BI = boundaryInteraction_[i];
+		if(BI.ppPairsAllocated()) BI.ppPairs().afterBroadSearch();
+		if(BI.pwPairsAllocated()) BI.pwPairs().afterBroadSearch();
 	}
+	contactListBoundaryTimer_.end();
 	
+	
+	{
+	boundaryInteractionTimer_.start();
+	std::array<bool,6> requireStep{
+		boundaryInteraction_[0].isBoundaryMaster(),
+		boundaryInteraction_[1].isBoundaryMaster(),
+		boundaryInteraction_[2].isBoundaryMaster(),
+		boundaryInteraction_[3].isBoundaryMaster(),
+		boundaryInteraction_[4].isBoundaryMaster(),
+		boundaryInteraction_[5].isBoundaryMaster()};
+	int step = 1;
+	const auto& cfModel = this->forceModel_();
+	while( std::any_of(requireStep.begin(), requireStep.end(), [](bool v) { return v==true; }))
+	{
+		for(uint32 i=0; i<6u; i++)
+		{
+			if(step==1u || requireStep[i] )
+			{
+				requireStep[i] = boundaryInteraction_[i].sphereSphereInteraction(
+					dt, 
+					this->forceModel_(),
+					step
+				);
+			}		
+		}
+		step++;
+	}
+	boundaryInteractionTimer_.pause();
+	}
 	ppInteractionTimer_.start();
 		sphereSphereInteraction();
 	ppInteractionTimer_.end();
@@ -239,14 +282,36 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 		sphereWallInteraction();
 	pwInteractionTimer_.end();
 
-	boundaryInteractionTimer_.start();
-	for(uint32 i=0; i<6u; i++)
 	{
-		boundaryInteraction_[i].sphereSphereInteraction(
-			dt, 
-			this->forceModel_());
+	boundaryInteractionTimer_.resume();
+	std::array<bool,6> requireStep{
+		!boundaryInteraction_[0].isBoundaryMaster(),
+		!boundaryInteraction_[1].isBoundaryMaster(),
+		!boundaryInteraction_[2].isBoundaryMaster(),
+		!boundaryInteraction_[3].isBoundaryMaster(),
+		!boundaryInteraction_[4].isBoundaryMaster(),
+		!boundaryInteraction_[5].isBoundaryMaster()};
+
+	int step = 2;
+	const auto& cfModel = this->forceModel_();
+	while(std::any_of(requireStep.begin(), requireStep.end(), [](bool v) { return v==true; }))
+	{
+		for(uint32 i=0; i<6u; i++)
+		{
+			if(requireStep[i])
+			{
+				requireStep[i] = boundaryInteraction_[i].sphereSphereInteraction(
+					dt, 
+					this->forceModel_(),
+					step
+				);
+			}		
+		}
+		step++;
 	}
 	boundaryInteractionTimer_.end();
+	}
+	
 	return true;
 }
 
diff --git a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.hpp b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.hpp
index 4ccb9836..74cb2442 100644
--- a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.hpp
+++ b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.hpp
@@ -96,9 +96,12 @@ private:
 	/// timer for managing contact lists (only inernal points)
 	Timer 		contactListMangementTimer_;
 
+	Timer 		boundaryContactSearchTimer_;
 	/// timer for boundary interaction time 
 	Timer 		boundaryInteractionTimer_;
 
+	Timer 		contactListBoundaryTimer_;
+
 	
 
 	bool createSphereInteraction();
diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp b/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp
index ea29aee8..df7b98b8 100644
--- a/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp
+++ b/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp
@@ -389,7 +389,9 @@ pFlow::sphereParticles::sphereParticles(
 	intPredictTimer_(
 		"Integration-predict", &this->timers() ),
 	intCorrectTimer_(
-		"Integration-correct", &this->timers() )
+		"Integration-correct", &this->timers() ),
+	fieldUpdateTimer_(
+		"fieldUpdate", &this->timers() )
 {
 	
 	auto intMethod = control.settingsDict().getVal<word>("integrationMethod");
@@ -514,13 +516,14 @@ bool pFlow::sphereParticles::beforeIteration()
 		rVelIntegration_().predict(dt,rVelocity_, rAcceleration_);
 	intPredictTimer_.end();
 
+	fieldUpdateTimer_.start();
 	propertyId_.updateBoundariesSlaveToMasterIfRequested();
 	diameter_.updateBoundariesSlaveToMasterIfRequested();
 	mass_.updateBoundariesSlaveToMasterIfRequested();
 	I_.updateBoundariesSlaveToMasterIfRequested();
 	rVelocity_.updateBoundariesSlaveToMasterIfRequested();
 	rAcceleration_.updateBoundariesSlaveToMasterIfRequested();
-	
+	fieldUpdateTimer_.end();
 	
 	return true;
 }
diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp b/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp
index fc0aa82c..0f78e498 100644
--- a/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp
+++ b/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp
@@ -79,7 +79,7 @@ private:
 	/// timer for integration computations (correction step)
 	Timer                  intCorrectTimer_;
 
-	
+	Timer  				   fieldUpdateTimer_;
 
 private:
 	bool initializeParticles();
diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp b/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp
index c57475ec..6ad16556 100644
--- a/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp
+++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp
@@ -38,6 +38,7 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 		*this,
 		zero3
 	),
+	velocityUpdateTimer_("velocity boundary update", &timers()),
 	integrationMethod_
 	(
 		control.settingsDict().getVal<word>("integrationMethod")
@@ -81,11 +82,11 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 
 bool pFlow::dynamicPointStructure::beforeIteration()
 {
-	return pointStructure::beforeIteration();
-	/*real dt = this->dt();
-
-    auto& acc = time().lookupObject<realx3PointField_D>("acceleration");
-	return predict(dt, acc);*/
+	if(!pointStructure::beforeIteration())return false;
+	velocityUpdateTimer_.start();
+	velocity_.updateBoundariesSlaveToMasterIfRequested();
+	velocityUpdateTimer_.end();
+	return true;
 }
 
 bool pFlow::dynamicPointStructure::iterate()
diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp b/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp
index 0e6a9fcd..f82f1dca 100644
--- a/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp
+++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp
@@ -44,6 +44,8 @@ private:
 
 	uniquePtr<integration>  integrationVel_ = nullptr;
 
+	Timer 					velocityUpdateTimer_;
+
 	/// @brief integration method for velocity and position
 	word 					integrationMethod_;
 
diff --git a/src/Particles/particles/particles.cpp b/src/Particles/particles/particles.cpp
index 187f45da..a8fc970a 100644
--- a/src/Particles/particles/particles.cpp
+++ b/src/Particles/particles/particles.cpp
@@ -74,7 +74,8 @@ pFlow::particles::particles(systemControl& control)
       dynPointStruct_,
       zero3
     ),
-    idHandler_(particleIdHandler::create(dynPointStruct_))
+    idHandler_(particleIdHandler::create(dynPointStruct_)),
+    baseFieldBoundaryUpdateTimer_("baseFieldBoundaryUpdate",&timers())
 {
 	this->addToSubscriber(dynPointStruct_);
 
@@ -84,18 +85,18 @@ pFlow::particles::particles(systemControl& control)
 bool
 pFlow::particles::beforeIteration()
 {
-	zeroForce();
-	zeroTorque();
-
 	if( !dynPointStruct_.beforeIteration())
   {
     return false;
   }
 
+  zeroForce();
+	zeroTorque();
+  baseFieldBoundaryUpdateTimer_.start();
   shapeIndex_.updateBoundariesSlaveToMasterIfRequested();
   accelertion_.updateBoundariesSlaveToMasterIfRequested();
   idHandler_().updateBoundariesSlaveToMasterIfRequested();
-
+  baseFieldBoundaryUpdateTimer_.end();
   return true;
 }
 
diff --git a/src/Particles/particles/particles.hpp b/src/Particles/particles/particles.hpp
index 1f3c7ab1..12e3da8c 100644
--- a/src/Particles/particles/particles.hpp
+++ b/src/Particles/particles/particles.hpp
@@ -54,6 +54,8 @@ private:
 	/// handling new ids for new particles 
 	uniquePtr<particleIdHandler> idHandler_ = nullptr;
 
+	Timer 						 baseFieldBoundaryUpdateTimer_;
+
 	/// messages for this objects
 	static inline const message  defaultMessage_{ message::DEFAULT };
 
diff --git a/src/phasicFlow/CMakeLists.txt b/src/phasicFlow/CMakeLists.txt
index de526e20..408fc042 100644
--- a/src/phasicFlow/CMakeLists.txt
+++ b/src/phasicFlow/CMakeLists.txt
@@ -94,6 +94,8 @@ structuredData/pointStructure/selectors/selectorStridedRange/selectorStridedRang
 structuredData/pointStructure/selectors/selectorRandomPoints/selectorRandomPoints.cpp
 structuredData/pointStructure/selectors/selectBox/selectBox.cpp
 structuredData/pointStructure/selectors/selectorGeometric/selectorGeometrics.cpp
+structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.cpp
+structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp
 
 triSurface/subSurface.cpp
 triSurface/triSurface.cpp
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingle.cpp b/src/phasicFlow/containers/VectorHD/VectorSingle.cpp
index 24248594..dd093270 100644
--- a/src/phasicFlow/containers/VectorHD/VectorSingle.cpp
+++ b/src/phasicFlow/containers/VectorHD/VectorSingle.cpp
@@ -226,6 +226,32 @@ pFlow::VectorSingle<T,MemorySpace>::VectorSingle
     }
 }
 
+template<typename T, typename MemorySpace>
+pFlow::VectorSingle<T,MemorySpace>::VectorSingle
+(
+    const word& name, 
+    const ViewType1D<T, MemorySpace>& src
+)
+:
+    VectorSingle(name, src.size(), src.size(), RESERVE())
+{
+    if constexpr(isTriviallyCopyable_)
+    {
+        copy(deviceView(), src);
+    }
+    else if constexpr( isHostAccessible_)
+    {
+        for(auto i=0u; i<size(); i++)
+        {
+            view_[i] = src[i];
+        }
+    }
+    else
+    {
+        static_assert("This constructor is not valid for non-trivially copyable data type on device memory");
+    }
+}
+
 template<typename T, typename MemorySpace>
 pFlow::VectorSingle<T,MemorySpace>& 
     pFlow::VectorSingle<T,MemorySpace>::operator = (const VectorSingle& rhs) 
@@ -367,6 +393,7 @@ template<typename T, typename MemorySpace>
 INLINE_FUNCTION_H 
 void pFlow::VectorSingle<T,MemorySpace>::reserve(uint32 cap) 
 {
+    if(cap == capacity() ) return;
     changeCapacity(cap);
 }
 
@@ -571,6 +598,40 @@ void pFlow::VectorSingle<T,MemorySpace>::assign
     }
 }
 
+
+template <typename T, typename MemorySpace>
+INLINE_FUNCTION_H
+void pFlow::VectorSingle<T, MemorySpace>::append(const ViewType1D<T,MemorySpace>& appVec)
+{
+    uint32 appSize  = appVec.size();
+	if(appSize == 0) return;
+
+	uint32 oldS = size();
+	uint32 newSize = oldS + appSize; 
+
+	changeSize(newSize);
+
+	auto appendView = Kokkos::subview(
+		view_,
+		Kokkos::make_pair<uint32>(oldS, newSize));
+	
+    if constexpr( isTriviallyCopyable_)
+    {
+        copy(appendView, appVec);
+    }
+    else if constexpr( isHostAccessible_)
+    {
+        for(auto i=0u; i<appVec.size(); i++)
+        {
+            appendView[i] = appVec[i];
+        }
+    }
+    else
+    {
+        static_assert("not a valid operation for this data type on device memory");
+    }
+}
+
 template <typename T, typename MemorySpace>
 INLINE_FUNCTION_H
 void pFlow::VectorSingle<T, MemorySpace>::append
@@ -813,7 +874,7 @@ bool pFlow::VectorSingle<T,MemorySpace>::insertSetElement
 
 template<typename T, typename MemorySpace>
 INLINE_FUNCTION_H
-bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(uint32IndexContainer indices)
+bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(const uint32IndexContainer& indices)
 {
 	if(indices.size() == 0)
 	{
@@ -831,10 +892,8 @@ bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(uint32IndexContainer indic
 		return false;
 	}
 
-	uint32 		newSize = indices.size();
+	uint32 	newSize = indices.size();
 	
-	setSize(newSize);
-
 	viewType 	sortedView(this->name(), newSize);
 
 	using policy = Kokkos::RangePolicy< execution_space,Kokkos::IndexType<uint32>>;
@@ -875,7 +934,24 @@ bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(uint32IndexContainer indic
 		Kokkos::fence();
 	}
 
-	copy(deviceView(), sortedView);
-	
+    setSize(newSize);
+
+    
+    if constexpr( isTriviallyCopyable_ )
+    {
+        copy(deviceView(), sortedView);
+    }
+    else if constexpr( isHostAccessible_)
+    {
+        for(auto i=0u; i<newSize; i++)
+        {
+            view_[i] = sortedView[i];
+        }   
+    }
+    else
+    {
+        static_assert("Not a valid operation for this data type on memory device");
+    }
+    
 	return true;
 }
\ No newline at end of file
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
index 256ccb6c..d30d1589 100644
--- a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
+++ b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
@@ -158,6 +158,9 @@ public:
 		
     	/// Copy construct with a new name (perform deep copy)
 		VectorSingle(const word& name, const VectorSingle& src);
+
+		/// Copy construct with a new name (perform deep copy)
+		VectorSingle(const word& name, const ViewType1D<T, MemorySpace>& src);
 		
 		/// Copy assignment (perform deep copy from rhs to *this)
 		VectorSingle& operator = (const VectorSingle& rhs) ;
@@ -307,6 +310,9 @@ public:
 			}
 		}
 
+		INLINE_FUNCTION_H
+		void append(const ViewType1D<T,MemorySpace>& appVec);
+
 		INLINE_FUNCTION_H
 		void append(const std::vector<T>& appVec);
 
@@ -331,7 +337,7 @@ public:
 			const ViewType1D<T, memory_space> vals);
 
 		INLINE_FUNCTION_H
-		bool reorderItems(uint32IndexContainer indices);
+		bool reorderItems(const uint32IndexContainer& indices);
 
 		/// @brief push a new element at the end (host call only)
 		///  resize if necessary and works on host accessible vector.
diff --git a/src/phasicFlow/containers/indexContainer/indexContainer.hpp b/src/phasicFlow/containers/indexContainer/indexContainer.hpp
index 72118609..4486cd2c 100644
--- a/src/phasicFlow/containers/indexContainer/indexContainer.hpp
+++ b/src/phasicFlow/containers/indexContainer/indexContainer.hpp
@@ -143,6 +143,16 @@ public:
 			indexContainer(ind.data(), ind.size())
 		{}
 
+		indexContainer(const DeviceViewType& ind):
+			size_(ind.size()),
+			views_("indexContainer", size_)
+		{
+			copy(views_.h_view, ind);
+			copy(views_.d_view, ind);
+			min_ = pFlow::min(views_.d_view, 0, size_);
+			max_ = pFlow::max(views_.d_view, 0, size_);
+		}
+
 		/// Copy
 		indexContainer(const indexContainer&) = default;
 
@@ -240,13 +250,13 @@ public:
 		void syncViews()
 		{
 			bool findMinMax = false;
-			if(views_.template need_sync<HostType>())
+			if(views_.template need_sync<DeviceType>())
 			{
 				Kokkos::deep_copy(views_.d_view, views_.h_view);
 				views_.clear_sync_state();
 				findMinMax = true;
 			}
-			else if(views_.template need_sync<DeviceType>())
+			else if(views_.template need_sync<HostType>())
 			{
 				Kokkos::deep_copy(views_.h_view, views_.d_view);
 				views_.clear_sync_state();
@@ -260,6 +270,12 @@ public:
 			}
 		}
 
+		void syncViews(uint32 newSize)
+		{
+			size_ = newSize;
+			syncViews();
+		}
+
 };
 
 
diff --git a/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp b/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp
index d30b3f18..9570f86b 100644
--- a/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp
+++ b/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp
@@ -121,7 +121,18 @@ public:
 		return true;
 	}
 
-	
+	inline
+	InternalFieldType& internal()
+	{
+		return internal_;
+	}
+
+	inline
+	const InternalFieldType& internal()const
+	{
+		return internal_;
+	}
+
 	FieldAccessType thisField()const
 	{
 		return FieldAccessType(
@@ -144,11 +155,6 @@ public:
 	virtual
 	const ProcVectorType& neighborProcField()const;
 
-	void fill(const std::any& val)override
-	{
-		return;
-	}
-
 	virtual
 	void fill(const T& val)
 	{
diff --git a/src/phasicFlow/containers/pointField/boundaryField/generalBoundary/generalBoundary.hpp b/src/phasicFlow/containers/pointField/boundaryField/generalBoundary/generalBoundary.hpp
index c55a7f71..6a1f6808 100644
--- a/src/phasicFlow/containers/pointField/boundaryField/generalBoundary/generalBoundary.hpp
+++ b/src/phasicFlow/containers/pointField/boundaryField/generalBoundary/generalBoundary.hpp
@@ -147,10 +147,6 @@ public:
     }
 
     const Time& time()const;
-
-	virtual
-	void fill(const std::any& val)=0;
-
 	
 };
 
diff --git a/src/phasicFlow/containers/pointField/internalField/internalField.cpp b/src/phasicFlow/containers/pointField/internalField/internalField.cpp
index c59faf3b..4b7c01e9 100644
--- a/src/phasicFlow/containers/pointField/internalField/internalField.cpp
+++ b/src/phasicFlow/containers/pointField/internalField/internalField.cpp
@@ -52,9 +52,29 @@ bool pFlow::internalField<T, MemorySpace>::insert(const anyList& varList)
 	}
 
 	return true;
-
 }
 
+template<class T, class MemorySpace>
+bool pFlow::internalField<T, MemorySpace>::rearrange(const anyList& varList)
+{
+	const word eventName = message::eventName(message::ITEM_REARRANGE);
+
+	const auto& indices = varList.getObject<uint32IndexContainer>(
+		eventName);
+	
+	field_.reserve( internalPoints_.capacity());
+	field_.resize(internalPoints_.size());
+	if(!field_.reorderItems(indices))
+	{
+		fatalErrorInFunction<<
+		"cannot reorder items in field "<< name()<<endl;
+		return false;
+	}
+
+	return true;
+}
+
+
 template<class T, class MemorySpace>
 pFlow::internalField<T, MemorySpace>::internalField
 (
@@ -166,15 +186,15 @@ bool pFlow::internalField<T, MemorySpace>:: hearChanges
 	{
 		// do nothing
 	}
-	if(msg.equivalentTo(message::ITEM_REARRANGE))
-	{
-		notImplementedFunction;
-		return false;
-	}
 	if(msg.equivalentTo(message::ITEM_INSERT))
 	{
 		return insert(varList);
 	}
+	if(msg.equivalentTo(message::ITEM_REARRANGE))
+	{
+		return rearrange(varList);
+
+	}
 	return true;
 }
 
diff --git a/src/phasicFlow/containers/pointField/internalField/internalField.hpp b/src/phasicFlow/containers/pointField/internalField/internalField.hpp
index 423f1ef7..ab26c663 100644
--- a/src/phasicFlow/containers/pointField/internalField/internalField.hpp
+++ b/src/phasicFlow/containers/pointField/internalField/internalField.hpp
@@ -73,6 +73,8 @@ protected:
 
 	bool insert(const anyList& varList);
 
+	bool rearrange(const anyList& varList);
+
 public:
 
 	internalField(
diff --git a/src/phasicFlow/eventManagement/message.hpp b/src/phasicFlow/eventManagement/message.hpp
index 3ece7411..0593221a 100644
--- a/src/phasicFlow/eventManagement/message.hpp
+++ b/src/phasicFlow/eventManagement/message.hpp
@@ -48,14 +48,16 @@ public:
 		BNDR_RESET		= 10,  // boundary indices reset entirely  
 		BNDR_DELETE 	= 11,  // boundary indices deleted
 		BNDR_APPEND		= 12,  // 
-		BNDR_PROCTRANS1 = 13,  // transfer of data between processors step 1 
-		BNDR_PROCTRANS2 = 14   // transfer of data between processors step 2
+		BNDR_PROCTRANSFER_SEND = 13,    // transfer of data between processors step 1 
+		BNDR_PROCTRANSFER_RECIEVE = 14,  // transfer of data between processors step 2
+		BNDR_PROCTRANSFER_WAITFILL = 15,  // wait for data and fill the field 
+		BNDR_PROC_SIZE_CHANGED = 16
 	};
 
 	
 private:
 
-	static constexpr size_t numberOfEvents_ = 15;
+	static constexpr size_t numberOfEvents_ = 17;
 
 	std::bitset<numberOfEvents_> events_{0x0000};
 	
@@ -76,7 +78,9 @@ private:
 		"deletedIndices",
 		"appendedIndices",
 		"transferredIndices",
-		"numberRecieved"
+		"numToRecieve",
+		"insertedIndices",
+		"size"
 	};
 
 public:
diff --git a/src/phasicFlow/repository/Time/baseTimeControl.cpp b/src/phasicFlow/repository/Time/baseTimeControl.cpp
index 8f4f2ef2..78d04225 100644
--- a/src/phasicFlow/repository/Time/baseTimeControl.cpp
+++ b/src/phasicFlow/repository/Time/baseTimeControl.cpp
@@ -65,6 +65,16 @@ pFlow::baseTimeControl::baseTimeControl
     
 }
 
+pFlow::baseTimeControl::baseTimeControl(int32 start, int32 end, int32 stride, const word &intervalPrefix)
+:
+	isTimeStep_(true),
+	iRange_(start, end, max(stride,1)),
+	intervalPrefix_(
+		intervalPrefix.size()==0uL? word("interval"): intervalPrefix+"Interval"
+	)
+{
+}
+
 bool pFlow::baseTimeControl::timeEvent(uint32 iter, real t, real dt) const
 {
 	if(isTimeStep_)
diff --git a/src/phasicFlow/repository/Time/baseTimeControl.hpp b/src/phasicFlow/repository/Time/baseTimeControl.hpp
index 192973fc..5a7634d5 100644
--- a/src/phasicFlow/repository/Time/baseTimeControl.hpp
+++ b/src/phasicFlow/repository/Time/baseTimeControl.hpp
@@ -46,6 +46,13 @@ public:
 	  real              defStartTime   = 0.0
 	);
 
+	baseTimeControl(
+		int32 start,
+		int32 end,
+		int32 stride,
+		const word& intervalPrefix = ""
+	);
+
 	inline bool isTimeStep() const
 	{
 		return isTimeStep_;
diff --git a/src/phasicFlow/streams/dataIO/dataIO.cpp b/src/phasicFlow/streams/dataIO/dataIO.cpp
index 39fbd9e0..851ec4a5 100644
--- a/src/phasicFlow/streams/dataIO/dataIO.cpp
+++ b/src/phasicFlow/streams/dataIO/dataIO.cpp
@@ -22,7 +22,7 @@ template<typename T>
 bool pFlow::writeDataAsciiBinary(iOstream& os, span<T> data)
 {
 
-	if( os.isBinary() )
+	if( std::is_trivially_copyable_v<T> && os.isBinary() )
 	{
 		// first write the number of data
 		uint64 numData = data.size();
@@ -83,7 +83,7 @@ bool pFlow::readDataAsciiBinary
 )
 {
 	
-	if(is.isBinary())
+	if(std::is_trivially_copyable_v<T> && is.isBinary())
 	{
 		data.clear();
 		// read length of data 			
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp
index f8307478..906d5fd9 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp
@@ -42,6 +42,7 @@ void pFlow::boundaryBase::setNewIndices
 	unSyncLists();
 }
 
+
 bool pFlow::boundaryBase::appendNewIndices
 (
 	const uint32Vector_D& newIndices
@@ -144,7 +145,7 @@ bool pFlow::boundaryBase::setRemoveKeepIndices
     return true;
 }
 
-bool pFlow::boundaryBase::transferPoints
+bool pFlow::boundaryBase::transferPointsToMirror
 (
 	uint32 numTransfer, 
 	const uint32Vector_D& transferMask, 
@@ -203,6 +204,8 @@ pFlow::boundaryBase::boundaryBase
 	indexList_(groupNames("indexList", dict.name())),
 	indexListHost_(groupNames("hostIndexList",dict.name())),
 	neighborLength_(dict.getVal<real>("neighborLength")),
+	updateInetrval_(dict.getVal<uint32>("updateInterval")),
+	boundaryExtntionLengthRatio_(dict.getVal<real>("boundaryExtntionLengthRatio")),
 	internal_(internal),
 	boundaries_(bndrs),
 	thisBoundaryIndex_(thisIndex),
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp
index 9b63c9df..978f674b 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp
@@ -68,6 +68,12 @@ private:
 	/// The length defined for creating neighbor list
 	real                   neighborLength_;
 
+	/// time steps between successive update of boundary lists  
+	uint32 				   updateInetrval_;
+
+	/// the extra boundary extension beyound actual limits of boundary 
+	real 				   boundaryExtntionLengthRatio_;
+
 	/// a reference to internal points
 	internalPoints&        internal_;
 
@@ -106,7 +112,7 @@ protected:
 		const uint32Vector_D& removeIndices,
 		const uint32Vector_D& keepIndices);
 	
-	bool transferPoints(
+	bool transferPointsToMirror(
 		uint32 numTransfer, 
 		const uint32Vector_D& transferMask,
 		uint32 transferBoundaryIndex,
@@ -126,6 +132,12 @@ protected:
 		}
 	}
 
+	/// Is this iter the right time for updating bounday list
+	bool boundaryListUpdate(uint32 iter)const
+	{
+		return iter%updateInetrval_ == 0u || iter == 0u;
+	}
+
 	/// Update this boundary data in two steps (1 and 2).
 	/// This is called after calling beforeIteration for 
 	/// all boundaries, so any particle addition, deletion,
@@ -133,7 +145,7 @@ protected:
 	/// This two-step update help to have a flexible mechanism
 	/// for data transfer, mostly for MPI related jobs. 
 	virtual 
-	bool updataBoundary(int step)
+	bool updataBoundaryData(int step)
 	{
 		return true;
 	}
@@ -145,11 +157,13 @@ protected:
 	/// to complete. false: if the operation is complete and no need for
 	/// additional step in operation. 
 	virtual 
-	bool transferData(int step)
+	bool transferData(uint32 iter, int step)
 	{
 		return false;
 	}
 	
+	 
+
 public:
 
 	TypeInfo("boundaryBase");
@@ -189,20 +203,21 @@ public:
 	/// The length from boundary plane into the domain 
 	/// where beyond that distance internal points exist.
 	/// By conventions is it always equal to neighborLength_  
+	inline
 	real neighborLengthIntoInternal()const
 	{
 		return neighborLength_;
 	}
 
 	/// The distance length from boundary plane 
-	/// where neighbor particles exist in that distance. 
+	/// where neighbor particles still exist in that distance. 
 	/// This length may be modified in each boundary type 
 	/// as required. In this case the boundaryExtensionLength
 	/// method should also be modified accordingly. 
 	virtual 
 	real neighborLength()const
 	{
-		return neighborLength_;
+		return (1+boundaryExtntionLengthRatio_)*neighborLength_;
 	}
 
 	/// The extention length (in vector form) for the boundary
@@ -213,7 +228,7 @@ public:
 	virtual 
 	realx3 boundaryExtensionLength()const
 	{
-		return {0,0,0};
+		return -boundaryExtntionLengthRatio_*neighborLength_ * boundaryPlane_.normal();
 	}
 
 	inline
@@ -351,6 +366,18 @@ public:
 	virtual
 	const realx3Vector_D& neighborProcPoints()const;
     
+	virtual 
+	uint32 numToTransfer()const
+	{
+		return 0u;
+	}
+
+	virtual 
+	uint32 numToRecieve()const
+	{
+		return 0u;
+	}
+
 	/// - static create 
 	static
 	uniquePtr<boundaryBase> create
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp
index f26438db..e051c803 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp
@@ -114,6 +114,16 @@ public:
 		return fieldVals_;
 	}
 
+	auto& indices()
+	{
+		return indices_;
+	}
+
+	const auto& indices()const
+	{
+		return indices_;
+	}
+
 	uint32 index(uint32 i)const
 	{
 		return indices_[i];
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp b/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp
index 3b753877..9937433c 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp
@@ -47,6 +47,9 @@ bool pFlow::boundaryExit::beforeIteration
 	real dt
 )
 {
+
+	if( !boundaryListUpdate(iterNum))return true;
+
 	// nothing have to be done
 	if(empty())
 	{
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryList.cpp b/src/phasicFlow/structuredData/boundaries/boundaryList.cpp
index 6cccb194..473e523d 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryList.cpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryList.cpp
@@ -91,18 +91,18 @@ pFlow::boundaryList::updateNeighborLists()
 pFlow::boundaryList::boundaryList(pointStructure& pStruct)
   : ListPtr<boundaryBase>(pStruct.simDomain().sizeOfBoundaries()),
     pStruct_(pStruct),
-    timeControl_(
-      pStruct.simDomain().subDict("boundaries"),
-      "update",
-      pStruct.currentTime()
-    )
+    neighborListUpdateInterval_(
+		pStruct.simDomain().subDict("boundaries").getVal<uint32>(
+			"neighborListUpdateInterval"
+		)
+	)
 {
 }
 
 bool
-pFlow::boundaryList::updateNeighborLists(uint32 iter, real t, real dt)
+pFlow::boundaryList::updateNeighborLists(uint32 iter, bool force)
 {
-	if (timeControl_.timeEvent(iter, t, dt))
+	if(iter%neighborListUpdateInterval_==0u || iter == 0u || force)
 	{
 		return updateNeighborLists();
 	}
@@ -155,15 +155,18 @@ pFlow::boundaryList::internalDomainBox() const
 }
 
 bool
-pFlow::boundaryList::beforeIteration(uint32 iter, real t, real dt)
+pFlow::boundaryList::beforeIteration(uint32 iter, real t, real dt, bool force)
 {
 	// it is time to update lists
-	if (timeControl_.timeEvent(iter, t, dt) && !updateNeighborLists())
+	if(!updateNeighborLists(iter , force) )
 	{
 		fatalErrorInFunction;
 		return false;
 	}
 
+	// this forces performing updating the list on each boundary
+	if(force) iter = 0; 
+	
 	for (auto bdry : *this)
 	{
 		if (!bdry->beforeIteration(iter, t, dt))
@@ -176,12 +179,12 @@ pFlow::boundaryList::beforeIteration(uint32 iter, real t, real dt)
 
 	for (auto bdry : *this)
 	{
-		bdry->updataBoundary(1);
+		bdry->updataBoundaryData(1);
 	}
 
 	for (auto bdry : *this)
 	{
-		bdry->updataBoundary(2);
+		bdry->updataBoundaryData(2);
 	}
 
 	return true;
@@ -215,20 +218,11 @@ pFlow::boundaryList::afterIteration(uint32 iter, real t, real dt)
 		{
 			if(requireStep[i])
 			{
-				requireStep[i] = this->operator[](i).transferData(step);
+				requireStep[i] = this->operator[](i).transferData(iter, step);
 			}
 		}
 		step++;
 	}
 
-	/*for (auto& bdry : *this)
-	{
-		if (!bdry->afterIteration(iter, t, dt))
-		{
-			fatalErrorInFunction << "Error in afterIteration in boundary "
-			                     << bdry->name() << endl;
-			return false;
-		}
-	}*/
 	return true;
 }
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryList.hpp b/src/phasicFlow/structuredData/boundaries/boundaryList.hpp
index 5eba6399..8d4c8710 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryList.hpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryList.hpp
@@ -17,15 +17,12 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
-
 #ifndef __boundaryList_hpp__
 #define __boundaryList_hpp__
 
 #include "domain.hpp"
 #include "boundaryBase.hpp"
 #include "ListPtr.hpp"
-#include "baseTimeControl.hpp"
 
 
 namespace pFlow
@@ -42,7 +39,7 @@ private:
 	//// - data members
 	pointStructure& pStruct_;
 
-	baseTimeControl timeControl_;
+	uint32 			neighborListUpdateInterval_;
 
 	domain          extendedDomain_;
 
@@ -72,7 +69,7 @@ public:
 
 	/// @brief update neighbor list of boundaries based on 
 	/// the time intervals
-	bool updateNeighborLists(uint32 iter, real t, real dt);
+	bool updateNeighborLists(uint32 iter, bool force = false);
 
 	bool createBoundaries();
 
@@ -94,11 +91,6 @@ public:
 		return ListPtr<boundaryBase>::operator[](i);
 	}
 
-	inline
-	const baseTimeControl& timeControl()const
-	{
-		return timeControl_;
-	}
 
 	inline
 	const auto& extendedDomain()const
@@ -114,7 +106,7 @@ public:
 
 	box internalDomainBox()const;
 	
-	bool beforeIteration(uint32 iter, real t, real dt);
+	bool beforeIteration(uint32 iter, real t, real dt, bool force = false);
 
 	bool iterate(uint32 iter, real t, real dt);
 
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryNone/boundaryNone.hpp b/src/phasicFlow/structuredData/boundaries/boundaryNone/boundaryNone.hpp
index 372a1c28..8e103d39 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryNone/boundaryNone.hpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryNone/boundaryNone.hpp
@@ -43,7 +43,7 @@ public:
 		boundaryList	&bndrs,
 		uint32 			thisIndex);
 
-	~boundaryNone() override= default;
+	~boundaryNone() final= default;
 	
 	add_vCtor
 	(
@@ -52,11 +52,11 @@ public:
 		dictionary
 	);
 
-	bool beforeIteration(uint32 iterNum, real t, real dt) override;
+	bool beforeIteration(uint32 iterNum, real t, real dt) final;
 
-	bool iterate(uint32 iterNum, real t, real dt) override;
+	bool iterate(uint32 iterNum, real t, real dt) final;
 
-	bool afterIteration(uint32 iterNum, real t, real dt) override;
+	bool afterIteration(uint32 iterNum, real t, real dt) final;
 
 
 };
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp
index 30fdc039..cbdecf71 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp
@@ -33,17 +33,20 @@ pFlow::boundaryPeriodic::boundaryPeriodic
 )
 :
 	boundaryBase(dict, bplane, internal, bndrs, thisIndex),
-	mirrorBoundaryIndex_(dict.getVal<uint32>("mirrorBoundaryIndex"))
-{}
+	mirrorBoundaryIndex_(dict.getVal<uint32>("mirrorBoundaryIndex")),
+	extensionLength_(dict.getVal<real>("boundaryExtntionLengthRatio"))
+{
+	extensionLength_ = max(extensionLength_, static_cast<real>(0.1));
+}
 
 pFlow::real pFlow::boundaryPeriodic::neighborLength() const
 {
-    return (1+extensionLength_)*boundaryBase::neighborLength();
+    return (1+extensionLength_)*neighborLengthIntoInternal();
 }
 
 pFlow::realx3 pFlow::boundaryPeriodic::boundaryExtensionLength() const
 {
-    return -extensionLength_*neighborLength()*boundaryBase::boundaryPlane().normal();
+    return -extensionLength_*neighborLengthIntoInternal()*boundaryBase::boundaryPlane().normal();
 }
 
 
@@ -57,6 +60,8 @@ bool pFlow::boundaryPeriodic::beforeIteration(
 	{
 		return true;
 	}
+
+	if( !boundaryListUpdate(iterNum))return true;
 	
 	uint32 s = size();
 	uint32Vector_D transferFlags("transferFlags",s+1, s+1, RESERVE()); 
@@ -92,7 +97,7 @@ bool pFlow::boundaryPeriodic::beforeIteration(
 	// to obtain the transfer vector 
 	realx3 transferVec = displacementVectroToMirror();
 	
-	return transferPoints
+	return transferPointsToMirror
 	(
 		numTransfered,
 		transferFlags, 
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp
index b4746a5e..88983fa3 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp
@@ -35,8 +35,7 @@ private:
 
 	uint32 mirrorBoundaryIndex_;
 
-	static
-	inline const real extensionLength_ = 0.1;
+	real extensionLength_ = 0.1;
 
 public:
 
diff --git a/src/phasicFlow/structuredData/cells/cells.hpp b/src/phasicFlow/structuredData/cells/cells.hpp
index c77b1a69..c4b50813 100644
--- a/src/phasicFlow/structuredData/cells/cells.hpp
+++ b/src/phasicFlow/structuredData/cells/cells.hpp
@@ -28,58 +28,40 @@ Licence:
 namespace pFlow
 {
 
-template<typename indexType>
 class cells
 {
-public:
-	
-	using CellType = triple<indexType>;
-
-protected:
+private:
 
 	// - domain 
-	box 		domain_{realx3(0.0), realx3(1.0)};
+	box 		domainBox_{realx3(0.0), realx3(1.0)};
 
 	// - cell size 
-	realx3 		cellSize_{1,1,1};
-
-	CellType 	numCells_{1,1,1};
+	real 		celldx_{1};
 
+	int32x3 	numCells_{1,1,1};
 	
 	// - protected methods
 	INLINE_FUNCTION_H
 	void calculate()
 	{
-		numCells_ = (domain_.maxPoint()-domain_.minPoint())/cellSize_ + realx3(1.0);
-		numCells_ = max( numCells_ , CellType(static_cast<indexType>(1)) );
+		numCells_ = (domainBox_.maxPoint()-domainBox_.minPoint())/celldx_ + realx3(1.0);
+		numCells_ = max( numCells_ , int32x3(1) );
 	}
 
 public:
 
 	INLINE_FUNCTION_HD
-	cells()
-	{}
+	cells() = default;
 
 	INLINE_FUNCTION_H
 	cells(const box& domain, real cellSize)
 	:
-		domain_(domain),
-		cellSize_(cellSize)
+		domainBox_(domain),
+		celldx_(cellSize)
 	{
 		calculate();
 	}
 
- 
-	INLINE_FUNCTION_H
-	cells(const box& domain, int32 nx, int32 ny, int32 nz)
-	:
-		domain_(domain),
-		cellSize_(
-			(domain_.maxPoint() - domain_.minPoint())/realx3(nx, ny,  nz)
-			),
-		numCells_(nx, ny, nz)
-	{}
-
 	INLINE_FUNCTION_HD
 	cells(const cells&) = default;
 
@@ -100,43 +82,36 @@ public:
 	INLINE_FUNCTION_H
 	void setCellSize(real cellSize)
 	{
-		cellSize_ = cellSize;
-		calculate();
-	}
-
-	INLINE_FUNCTION_H
-	void setCellSize(realx3 cellSize)
-	{
-		cellSize_ = cellSize;
+		celldx_ = cellSize;
 		calculate();
 	}
 
 	INLINE_FUNCTION_HD
-	realx3 cellSize()const
+	real cellSize()const
 	{
-		return cellSize_;
+		return celldx_;
 	}
 
 	INLINE_FUNCTION_HD
-	const CellType& numCells()const
+	const int32x3& numCells()const
 	{
 		return numCells_;
 	}
 
 	INLINE_FUNCTION_HD
-	indexType nx()const
+	int32 nx()const
 	{
 		return numCells_.x();
 	}
 
 	INLINE_FUNCTION_HD
-	indexType ny()const
+	int32 ny()const
 	{
 		return numCells_.y();
 	}
 
 	INLINE_FUNCTION_HD
-	indexType nz()const
+	int32 nz()const
 	{
 		return numCells_.z();
 	}
@@ -149,22 +124,21 @@ public:
 				static_cast<int64>(numCells_.z());
 	}
 
-	const auto& domain()const
+	const auto& domainBox()const
 	{
-		return domain_;
+		return domainBox_;
 	}
 	
 	INLINE_FUNCTION_HD
-	CellType pointIndex(const realx3& p)const
+	int32x3 pointIndex(const realx3& p)const
 	{
-		return CellType( (p - domain_.minPoint())/cellSize_ );
+		return int32x3( (p - domainBox_.minPoint())/celldx_ );
 	}
 
 	INLINE_FUNCTION_HD
-	bool pointIndexInDomain(const realx3 p, CellType& index)const
+	bool pointIndexInDomain(const realx3 p, int32x3& index)const
 	{
-		if( !domain_.isInside(p) ) return false;
-		
+		if(!inDomain(p))return false;
 		index = this->pointIndex(p);
 		return true;
 	}
@@ -172,11 +146,11 @@ public:
 	INLINE_FUNCTION_HD
 	bool inDomain(const realx3& p)const
 	{	
-		return domain_.isInside(p);
+		return domainBox_.isInside(p);
 	}
 
 	INLINE_FUNCTION_HD
-	bool isInRange(const CellType& cell)const
+	bool inCellRange(const int32x3& cell)const
 	{
 		if(cell.x()<0)return false;
 		if(cell.y()<0)return false;
@@ -188,7 +162,7 @@ public:
 	}
 
 	INLINE_FUNCTION_HD
-	bool isInRange(indexType i, indexType j, indexType k)const
+	bool inCellRange(int32 i, int32 j, int32 k)const
 	{
 		if(i<0)return false;
 		if(j<0)return false;
@@ -199,22 +173,7 @@ public:
 		return true;
 	}
 
-	INLINE_FUNCTION_HD
-	void extendBox(
-		const CellType& p1,
-		const CellType& p2,
-		const CellType& p3,
-		indexType extent,
-		CellType& minP,
-		CellType& maxP)const
-	{
-		minP = min( min( p1, p2), p3)-extent;
-		maxP = max( max( p1, p2), p3)+extent;
-
-		minP = bound(minP);
-		maxP = bound(maxP);
-	}
-
+	
 	INLINE_FUNCTION_HD
 	void extendBox(
 		const realx3& p1,
@@ -224,17 +183,17 @@ public:
 		realx3& minP,
 		realx3& maxP)const
 	{
-		minP = min(min(p1,p2),p3) - extent*cellSize_ ;
-		maxP = max(max(p1,p2),p3) + extent*cellSize_ ;
+		minP = min(min(p1,p2),p3) - extent*celldx_ ;
+		maxP = max(max(p1,p2),p3) + extent*celldx_ ;
 
 		minP = bound(minP);
 		maxP = bound(maxP);
 	}
 
 	INLINE_FUNCTION_HD
-	CellType bound(CellType p)const
+	int32x3 bound(int32x3 p)const
 	{
-		return CellType(
+		return int32x3(
 			min( numCells_.x()-1, max(0,p.x())),
 			min( numCells_.y()-1, max(0,p.y())),
 			min( numCells_.z()-1, max(0,p.z()))
@@ -245,9 +204,9 @@ public:
 	realx3 bound(realx3 p)const
 	{
 		return realx3(
-			min( domain_.maxPoint().x(), max(domain_.minPoint().x(),p.x())),
-			min( domain_.maxPoint().y(), max(domain_.minPoint().y(),p.y())),
-			min( domain_.maxPoint().z(), max(domain_.minPoint().z(),p.z()))
+			min( domainBox_.maxPoint().x(), max(domainBox_.minPoint().x(),p.x())),
+			min( domainBox_.maxPoint().y(), max(domainBox_.minPoint().y(),p.y())),
+			min( domainBox_.maxPoint().z(), max(domainBox_.minPoint().z(),p.z()))
 			);
 	}
 };	
diff --git a/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp b/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp
index e94d111e..a84325a9 100644
--- a/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp
+++ b/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp
@@ -31,7 +31,10 @@ bool pFlow::regularSimulationDomain::createBoundaryDicts()
     this->addDict("regularBoundaries", boundaries);
     auto& rbBoundaries = this->subDict("regularBoundaries");
     
-	real neighborLength = boundaries.getVal<real>("neighborLength");
+	auto neighborLength = boundaries.getVal<real>("neighborLength");
+	auto boundaryExtntionLengthRatio = 
+		boundaries.getValOrSet<real>("boundaryExtntionLengthRatio", 0.1);
+	auto updateIntercal = boundaries.getValOrSet<uint32>("updateInterval", 1u); 
 
 	for(uint32 i=0; i<sizeOfBoundaries(); i++)
 	{
@@ -51,6 +54,14 @@ bool pFlow::regularSimulationDomain::createBoundaryDicts()
 			return false;
 		}
 
+		if(!bDict.addOrReplace("updateInterval", updateIntercal))
+		{
+			fatalErrorInFunction<<"error in adding updateIntercal to "<< bName <<
+			"in dictionary "<< boundaries.globalName()<<endl;
+		}
+
+		bDict.addOrReplace("boundaryExtntionLengthRatio", boundaryExtntionLengthRatio);
+
 		if(!bDict.add("neighborProcessorNo",(uint32) processors::globalRank()))
 		{
 			fatalErrorInFunction<<"error in adding neighborProcessorNo to "<< bName <<
diff --git a/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.cpp b/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.cpp
index 6155a33d..4dc8ff0d 100644
--- a/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.cpp
+++ b/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.cpp
@@ -142,13 +142,24 @@ bool pFlow::internalPoints::changePointsFlagPosition
 	return true;
 }
 
+bool pFlow::internalPoints::sortPoints(const uint32IndexContainer &sortedIndices)
+{
+	if(!pointPosition_.reorderItems(sortedIndices))
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+    
+	pFlagsD_.resetFlags(pFlagsD_.capacity(), 0u, sortedIndices.size());
+	unSyncFlag();
 
+	return true;
+}
 
 pFlow::internalPoints::internalPoints()
-:
-	subscriber("internalPoints"),
-	pointPosition_("internalPoints", "internalPoints", initialCapacity_, 0, RESERVE()),
-	pFlagsD_(initialCapacity_, 0 , 0)
+    : subscriber("internalPoints"),
+      pointPosition_("internalPoints", "internalPoints", initialCapacity_, 0, RESERVE()),
+      pFlagsD_(initialCapacity_, 0, 0)
 {}
 
 
@@ -398,6 +409,108 @@ pFlow::internalPoints::insertPoints(
 	return true;
 }
 
+bool pFlow::internalPoints::insertPointsOnly(
+	const realx3Vector_D &points,
+	message& msg, 
+	anyList &varList
+)
+{
+
+	uint32 numNew = static_cast<uint32>(points.size());
+
+	auto aRange = pFlagsD_.activeRange();
+	uint32 emptySpots = pFlagsD_.capacity() - pFlagsD_.numActive();
+
+	if(emptySpots!= 0) emptySpots--;
+	
+	if( numNew > emptySpots )
+	{
+		// increase the capacity to hold new points 
+		aRange = pFlagsD_.activeRange();
+		uint32 newCap = pFlagsD_.changeCapacity(numNew);
+		unSyncFlag();
+		varList.emplaceBack(
+			msg.addAndName(message::CAP_CHANGED),
+			newCap);
+	}
+	
+	
+	// first check if it is possible to add to the beggining of list
+	if(numNew <= aRange.start())
+	{
+		varList.emplaceBack<uint32IndexContainer>(
+			msg.addAndName(message::ITEM_INSERT),
+			0u, numNew);
+	}// check if it is possible to add to the end of the list  
+	else if( numNew <= pFlagsD_.capacity() - aRange.end() )
+	{	 	
+		varList.emplaceBack<uint32IndexContainer>(
+			msg.addAndName(message::ITEM_INSERT),
+			aRange.end(), aRange.end()+numNew);
+	}// we should fill the scattered empty spots
+	else
+	{
+		pOutput<<"numNew to be inserted "<< numNew <<endl;
+		auto newIndices = pFlagsD_.getEmptyPoints(numNew);
+		if(numNew != newIndices.size())
+		{
+			fatalErrorInFunction<<"not enough empty points in pointFlag"<<
+			numNew<< " "<<newIndices.size() <<endl;
+			pOutput<< pFlagsD_.capacity()<<endl;
+			pOutput<< pFlagsD_.numActive()<<endl;
+			return false;
+		}
+		pOutput<<newIndices<<endl;
+		varList.emplaceBack<uint32IndexContainer>(
+			msg.addAndName(message::ITEM_INSERT),
+			newIndices
+		);
+	}
+
+	const auto& indices = varList.getObject<uint32IndexContainer>(
+		message::eventName(message::ITEM_INSERT)
+	);
+
+	auto nAdded = pFlagsD_.addInternalPoints(indices.deviceView());
+	unSyncFlag();
+
+	if(nAdded != numNew )
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+	pointPosition_.reserve( pFlagsD_.capacity() );
+
+	if(!pointPosition_.insertSetElement(indices, points.deviceView()))
+	{
+		fatalErrorInFunction<<
+		"Error in inserting new positions into pointPosition"<<
+		" internalPoints field"<<endl;
+		return false;
+	}
+
+	auto newARange = pFlagsD_.activeRange();
+
+	if( aRange.end() != newARange.end() )
+	{
+		varList.emplaceBack(
+			msg.addAndName(message::RANGE_CHANGED),
+			newARange);
+		varList.emplaceBack(
+			msg.addAndName(message::SIZE_CHANGED),
+			newARange.end());
+	}
+	else if(aRange.start() != newARange.start())
+	{
+		varList.emplaceBack(
+			msg.addAndName(message::RANGE_CHANGED),
+			newARange);
+	}
+	
+    return true;
+}
+
 bool pFlow::internalPoints::read
 (
 	iIstream& is
diff --git a/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.hpp b/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.hpp
index 2a8260af..c03a027d 100644
--- a/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.hpp
+++ b/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.hpp
@@ -101,6 +101,8 @@ protected:
 		pFlagsD_ = pFlagTypeDevice(cap, start, end);
 	}
 
+	bool sortPoints(const uint32IndexContainer& sortedIndices);
+
 public:
 
 	//friend class dynamicinternalPoints;
@@ -235,6 +237,8 @@ public:
 		
 		
 		bool insertPoints(const realx3Vector& points, anyList& varList);
+
+		bool insertPointsOnly(const realx3Vector_D& points, message& msg, anyList& varList);
 		
 	
 	//// - IO operations 
diff --git a/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlag.hpp b/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlag.hpp
index 28509010..609bd0f6 100644
--- a/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlag.hpp
+++ b/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlag.hpp
@@ -57,6 +57,8 @@ class pointFlag
 
 protected:
 	
+	friend class internalPoints;
+
 	viewType 	flags_;
 
 	uint32 		numActive_ 		= 0;
@@ -108,9 +110,29 @@ protected:
 		return flg;
 	}
 
+	void resetFlags(uint32 cap, uint32 start, uint32 end)
+	{
+		if(cap != capacity() )
+		{
+			reallocFill(flags_, cap, static_cast<uint8>(Flag::DELETED));
+		}
+		else
+		{
+			fill(flags_, end, cap, static_cast<uint8>(Flag::DELETED));
+		}
+
+		activeRange_ = {start, end};
+		numActive_ = activeRange_.numElements();
+		isAllActive_ = true;
+		fill(flags_, activeRange_, static_cast<uint8>(Flag::INTERNAL));
+
+		nLeft_ = nRight_ = 0;
+		nBottom_ = nTop_ = 0;
+		nRear_ = nFront_ = 0;
+	}
+
 public:
 
-	friend class internalPoints;
 
 	pointFlag() = default;
 
@@ -283,6 +305,9 @@ public:
 	}
 
 	ViewType1D<uint32, memory_space> getActivePoints();
+
+
+	ViewType1D<uint32, memory_space> getEmptyPoints(uint32 numToGet)const;
 	
 
 	/// @brief Loop over the active points and mark those out of the box
diff --git a/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp b/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp
index 91764cf6..de8940ba 100644
--- a/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp
+++ b/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp
@@ -50,6 +50,51 @@ pFlow::ViewType1D<pFlow::uint32, typename pFlow::pointFlag<ExecutionSpace>::memo
 	return aPoints;
 }
 
+template<typename ExecutionSpace>
+pFlow::ViewType1D<pFlow::uint32, typename pFlow::pointFlag<ExecutionSpace>::memory_space>
+	pFlow::pointFlag<ExecutionSpace>::getEmptyPoints(uint32 numToGet)const
+{
+
+	uint32 cap = capacity();
+	using rpAPoints = Kokkos::RangePolicy<execution_space,  
+			Kokkos::IndexType<uint32>>;
+
+	ViewType1D<uint32,memory_space> indices("indices", cap+1);
+
+	ViewType1D<uint32,memory_space> emptyPoints("emptyPoints", numToGet);
+
+	
+	Kokkos::parallel_for(
+		"getEmptyPoints",
+		rPolicy(0, cap),
+		LAMBDA_HD(uint32 i){
+			indices(i) = flags_[i] == DELETED;
+	});
+	Kokkos::fence();
+
+	exclusiveScan(indices, 0u, cap, indices, 0u);
+	uint32 numFound = 0;
+	Kokkos::parallel_reduce(
+		"fillEmptyPoints",
+		rpAPoints(0, cap-1),
+		LAMBDA_HD(uint32 i, uint32& numFoundUpdate){
+			if(indices(i)!= indices(i+1) && indices(i)< numToGet)
+			{
+				emptyPoints(indices(i)) = i;
+				numFoundUpdate++;
+			}
+		},
+		numFound
+	);
+	
+	if(numFound < numToGet)
+	{
+		return Kokkos::subview(emptyPoints, Kokkos::make_pair(0u, numFound));
+	}
+	return emptyPoints;
+
+}
+
 
 template<typename ExecutionSpace>
 pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markOutOfBoxDelete
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.cpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.cpp
similarity index 59%
rename from src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.cpp
rename to src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.cpp
index 63969747..2405bca4 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.cpp
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.cpp
@@ -17,66 +17,60 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
 #include "mortonIndexing.hpp"
 #include "cells.hpp"
 
-#include "streams.hpp"
 
-bool pFlow::getSortedIndex(
+pFlow::uint32IndexContainer pFlow::getSortedIndices(
 	box boundingBox,
-	real dx,
-	range activeRange, 
-	ViewType1D<realx3> pos, 
-	ViewType1D<int8> flag,
-	int32IndexContainer& sortedIndex)
+	real dx,  
+	const ViewType1D<realx3>& pos, 
+	const pFlagTypeDevice& flag
+)
 {
+	if(flag.numActive() == 0u)return uint32IndexContainer();
 
 	// obtain the morton code of the particles
-	cells<size_t> allCells( boundingBox, dx);
-	int32IndexContainer index(activeRange.first, activeRange.second);
+	cells allCells( boundingBox, dx);
+	auto aRange = flag.activeRange();
 
-	ViewType1D<uint64_t> mortonCode("mortonCode", activeRange.second);
+	uint32IndexContainer sortedIndex(aRange.start(), aRange.end());
 
-	output<<"before first kernel"<<endl;;
+	ViewType1D<uint64_t> mortonCode("mortonCode", aRange.end());
 
-	using rpMorton = 
-		Kokkos::RangePolicy<Kokkos::IndexType<int32>>;
-	int32 numActive = 0;
-	Kokkos::parallel_reduce
+
+	Kokkos::parallel_for
 	(
 		"mortonIndexing::getIndex::morton",
-		rpMorton(activeRange.first, activeRange.second),
-		LAMBDA_HD(int32 i, int32& sumToUpdate){
-			if( flag[i] == 1 ) // active point 
+		deviceRPolicyStatic(aRange.start(), aRange.end()),
+		LAMBDA_HD(uint32 i){
+			if( flag.isActive(i)) // active point 
 			{
 				auto cellInd = allCells.pointIndex(pos[i]);
 				mortonCode[i] = xyzToMortonCode64(cellInd.x(), cellInd.y(), cellInd.z());
-				sumToUpdate++;
 			}else
 			{
 				mortonCode[i] = xyzToMortonCode64
 				(
-					static_cast<uint64_t>(-1),
-					static_cast<uint64_t>(-1),
-					static_cast<uint64_t>(-1)
+					largestPosInt32,
+					largestPosInt32,
+					largestPosInt32
 				);
 			}
-		},
-		numActive
+		}
 	);
 	
+	Kokkos::fence();
+
 	permuteSort(
 		mortonCode, 
-		activeRange.first, 
-		activeRange.second,
-		index.deviceView(),
+		aRange.start(), 
+		aRange.end(),
+		sortedIndex.deviceView(),
 		0 );
-	index.modifyOnDevice();
-	index.setSize(numActive);
-	index.syncViews();
 
-	sortedIndex = index;
+	sortedIndex.modifyOnDevice();
+	sortedIndex.syncViews(flag.numActive());
 
-	return true;
+	return sortedIndex;
 }
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.hpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.hpp
similarity index 94%
rename from src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.hpp
rename to src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.hpp
index 9ccdc83e..0d8c58d7 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointStructure/mortonIndexing.hpp
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/mortonIndexing.hpp
@@ -17,24 +17,22 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
 #ifndef __mortonIndexing_hpp__
 #define __mortonIndexing_hpp__
 
 #include "types.hpp"
 #include "box.hpp"
 #include "indexContainer.hpp"
+#include "pointFlag.hpp"
 
 namespace pFlow
 {
 
-bool getSortedIndex(
+uint32IndexContainer getSortedIndices(
 	box boundingBox,
-	real dx, 
-	range activeRange, 
-	ViewType1D<realx3> pos, 
-	ViewType1D<int8> flag,
-	int32IndexContainer& sortedIndex);
+	real dx,  
+	const ViewType1D<realx3>& pos, 
+	const pFlagTypeDevice& flag);
 
 
 INLINE_FUNCTION_HD
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp
new file mode 100644
index 00000000..c8dda774
--- /dev/null
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.cpp
@@ -0,0 +1,52 @@
+/*------------------------------- 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 "pointSorting.hpp"
+#include "mortonIndexing.hpp"
+
+pFlow::pointSorting::pointSorting(const dictionary & dict)
+:
+    performSorting_(dict.getValOrSet("active", Logical(false))),
+    timeControl_(
+        performSorting_()? 
+            baseTimeControl(dict, "sorting"): 
+            baseTimeControl(0,1,1, "sorting")
+    ),
+    dx_(
+        performSorting_()?
+            dict.getVal<real>("dx"):
+            1.0
+    )
+{
+    if( performSorting_() )
+        REPORT(1)<<"Point sorting is "<<Yellow_Text("active")<<" in simulation"<<END_REPORT;
+    else
+       REPORT(1)<<"Point sorting is "<<Yellow_Text("inactive")<<" in simulation"<<END_REPORT;
+
+}
+pFlow::uint32IndexContainer 
+pFlow::pointSorting::getSortedIndices(
+    const box& boundingBox, 
+    const ViewType1D<realx3> &pos, 
+    const pFlagTypeDevice &flag
+) const
+{
+    return pFlow::getSortedIndices(boundingBox, dx_, pos, flag);
+}
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.hpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.hpp
new file mode 100644
index 00000000..d8d9f195
--- /dev/null
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointSorting/pointSorting.hpp
@@ -0,0 +1,68 @@
+/*------------------------------- 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 __pointSorting_hpp__
+#define __pointSorting_hpp__
+
+#include "baseTimeControl.hpp"
+#include "indexContainer.hpp"
+#include "pointFlag.hpp"
+
+
+namespace pFlow
+{
+
+class box;
+
+
+class pointSorting
+{
+private:
+
+	Logical                 performSorting_;
+
+	baseTimeControl         timeControl_;
+
+	real 					dx_;
+
+public:
+
+	explicit pointSorting(const dictionary& dict);
+
+	bool performSorting()const
+	{
+		return performSorting_();
+	}
+
+	bool sortTime(uint32 iter, real t, real dt)const
+	{
+		return performSorting_() && timeControl_.timeEvent(iter, t, dt);
+	}
+
+	uint32IndexContainer getSortedIndices(
+		const box& boundingBox, 
+		const ViewType1D<realx3>& pos,
+		const pFlagTypeDevice& flag
+	)const;
+	
+};
+
+}
+
+#endif
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.cpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.cpp
index 9a20395c..bc0e7976 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.cpp
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.cpp
@@ -22,6 +22,7 @@ Licence:
 #include "pointStructure.hpp"
 #include "systemControl.hpp"
 #include "vocabs.hpp"
+#include "anyList.hpp"
 
 bool pFlow::pointStructure::setupPointStructure(const realx3Vector& points)
 {
@@ -73,6 +74,7 @@ bool pFlow::pointStructure::setupPointStructure(const PointsTypeHost &points)
     }
 
     boundaries_.createBoundaries();
+    boundaries_.updateNeighborLists(0u, true);
 
     return true;
 }
@@ -111,10 +113,13 @@ pFlow::pointStructure::pointStructure
 	(
 		simulationDomain::create(control)
 	),
+    pointSorting_(simulationDomain_->subDictOrCreate("pointSorting")),
 	boundaries_
 	(
         *this
-	)
+	),
+    boundaryUpdateTimer_("boundaryUpdate", &timers()),
+    boundaryDataTransferTimer_("boundaryDataTransferTimer", &timers())
 {
     REPORT(1)<< "Reading "<< Green_Text("pointStructure")<<
     " from "<<IOobject::path()<<END_REPORT;
@@ -150,6 +155,7 @@ pFlow::pointStructure::pointStructure(
 	(
 		simulationDomain::create(control)
 	),
+    pointSorting_(simulationDomain_->subDictOrCreate("pointSorting")),
 	boundaries_
 	(
 		*this
@@ -165,12 +171,58 @@ pFlow::pointStructure::pointStructure(
 
 bool pFlow::pointStructure::beforeIteration()
 {
-    if( !boundaries_.beforeIteration(currentIter(), currentTime(), dt()) )
+    uint32 iter = currentIter();
+    real t = currentTime();
+    real deltat = dt();
+
+    if(pointSorting_.sortTime(iter, t, deltat))
     {
-        fatalErrorInFunction<<
-        "Unable to perform beforeIteration for boundaries"<<endl;
-        return false;
+        auto sortedIndices = pointSorting_.getSortedIndices(
+            simulationDomain_().globalBox(),
+            pointPositionDevice(),
+            activePointsMaskDevice()
+        );
+        
+        if( !sortPoints(sortedIndices) )
+        {
+            fatalErrorInFunction;
+            return false;
+        }
+        
+        boundaryUpdateTimer_.start();
+        boundaries_.beforeIteration(iter, t, deltat, true); 
+        boundaryUpdateTimer_.end();
+
+        INFORMATION<<"Reordering of particles has been done. New active range for particles is "<<
+        activeRange()<<END_INFO; 
+        message msg;
+        anyList varList;
+
+        varList.emplaceBack(
+            msg.addAndName(message::ITEM_REARRANGE),
+            sortedIndices);
+
+        if(!notify(iter, t, deltat, msg, varList))
+        {
+            fatalErrorInFunction<<
+            "cannot notify for reordering items."<<endl;
+            return false;
+        }
+        return true;
+        // notify others about this change 
     }
+    else
+    {
+        boundaryUpdateTimer_.start();
+        if( !boundaries_.beforeIteration(iter, t, deltat) )
+        {
+            fatalErrorInFunction<<
+            "Unable to perform beforeIteration for boundaries"<<endl;
+            return false;
+        }
+        boundaryUpdateTimer_.end();
+    }
+    
     return true;
 }
 
@@ -188,12 +240,14 @@ bool pFlow::pointStructure::iterate()
 
 bool pFlow::pointStructure::afterIteration()
 {
+    boundaryDataTransferTimer_.start();
     if( !boundaries_.afterIteration(currentIter(), currentTime(), dt()) )
     {
         fatalErrorInFunction<<
         "Unable to perform afterIteration for boundaries"<<endl;
         return false;
     }
+    boundaryDataTransferTimer_.end();
     return true;
 }
 
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.hpp b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.hpp
index 065e58cd..e1a60ff2 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.hpp
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure/pointStructure.hpp
@@ -24,6 +24,7 @@ Licence:
 #include "demComponent.hpp"
 #include "internalPoints.hpp"
 #include "simulationDomain.hpp"
+#include "pointSorting.hpp"
 #include "boundaryList.hpp"
 #include "streams.hpp"
 
@@ -53,8 +54,16 @@ private:
 	//// - data members 
 	uniquePtr<simulationDomain> simulationDomain_ = nullptr;
 
+	pointSorting 				pointSorting_;
+
 	boundaryList 				boundaries_;
-    
+	
+	Timer 						boundaryUpdateTimer_;
+
+	Timer  						boundaryDataTransferTimer_;
+
+	
+
 	bool setupPointStructure(const realx3Vector& points);
 
     bool setupPointStructure(const PointsTypeHost& points);

From 4e8b9215143b19a254bc4e1f62ea4a589535dc32 Mon Sep 17 00:00:00 2001
From: Hamidreza Norouzi <hamid.r.norouzi@gmail.com>
Date: Sat, 18 May 2024 18:40:25 +0330
Subject: [PATCH 2/2] bug remove for GPU run after CPU MPI parallelization

- specialization of VectorSingle for word
- dummyFile creation to solve write to file in MPI mode
---
 solvers/iterateGeometry/iterateGeometry.cpp   |  10 +-
 .../createDEMComponents.hpp                   |   2 +-
 .../iterateSphereParticles.cpp                |  10 +-
 .../sphereGranFlow/createDEMComponents.hpp    |   6 +-
 solvers/sphereGranFlow/sphereGranFlow.cpp     |  13 +-
 src/Geometry/geometry/geometry.cpp            |   4 +-
 .../geometryMotion/geometryMotion.cpp         |  34 +-
 src/Interaction/CMakeLists.txt                |   4 +-
 .../contactLists/sortedContactList.hpp        |   4 +-
 .../contactLists/unsortedContactList.hpp      |   6 +-
 .../contactLists/unsortedPairs.hpp            |   6 +-
 .../ppwBndryContactSearchKernels.cpp          |   4 +-
 .../wallBoundaryContactSearch.cpp             |  13 +-
 .../methods/cellBased/NBS/NBSLoop.hpp         |   4 +-
 .../methods/cellBased/NBS/cellsWallLevel0.cpp |  37 +-
 .../sphereInteractionKernels.hpp              |   2 +-
 .../entities/stationary/stationary.hpp        |   2 +-
 .../rotatingAxisMotion/rotatingAxisMotion.hpp |   6 +-
 .../stationaryWall/stationaryWall.hpp         |   2 +
 .../vibratingMotion/vibratingMotion.hpp       |   6 +-
 .../collisionCheck/collisionCheck.cpp         |   4 +-
 .../Insertion/insertion/insertion.cpp         |   3 -
 .../Insertion/insertion/insertion.hpp         |   6 +-
 .../sphereParticles/sphereParticles.hpp       |   7 +-
 .../regularParticleIdHandler.cpp              |   7 +-
 .../regularParticleIdHandler.hpp              |   2 +-
 .../particles/shape/baseShapeNames.hpp        |   4 +-
 src/phasicFlow/CMakeLists.txt                 |   1 +
 src/phasicFlow/Kokkos/ViewAlgorithms.hpp      |   4 +-
 src/phasicFlow/commandLine/CLI/TypeTools.hpp  |   4 +-
 src/phasicFlow/containers/Field/Field.hpp     |   1 +
 src/phasicFlow/containers/Field/Fields.cpp    |   3 +-
 .../containers/List/ListPtr/ListPtrI.hpp      |  25 +-
 .../containers/VectorHD/VectorSingle.cpp      | 239 ++------
 .../containers/VectorHD/VectorSingle.hpp      |  25 +-
 .../containers/VectorHD/VectorSingles.hpp     |   3 +
 .../containers/VectorHD/wordVectorHost.cpp    |  23 +
 .../containers/VectorHD/wordVectorHost.hpp    | 558 ++++++++++++++++++
 .../boundaryField/boundaryField.hpp           |  29 +-
 .../boundaryField/boundaryFieldList.hpp       |   4 +-
 src/phasicFlow/dictionary/dictionary.hpp      |   4 +-
 src/phasicFlow/dictionary/fileDictionary.hpp  |   3 +
 src/phasicFlow/globals/pFlowMacros.hpp        |   4 +
 src/phasicFlow/processors/localProcessors.hpp |   9 +
 src/phasicFlow/processors/processors.cpp      |   1 -
 .../repository/IOobject/IOfileHeader.cpp      |  22 +-
 .../repository/IOobject/IOfileHeader.hpp      |   2 +
 .../repository/IOobject/IOobject.cpp          |  29 +-
 .../repository/Time/baseTimeControl.cpp       |   2 +-
 .../repository/Time/timeControl.cpp           |  12 +-
 src/phasicFlow/smartPointers/uniquePtr.hpp    |   7 +-
 src/phasicFlow/streams/TStream/iTstream.cpp   |   2 +-
 .../boundaryBase/scatteredFieldAccess.hpp     |   1 +
 .../boundaryReflective/boundaryReflective.cpp |   4 +-
 .../structuredData/cylinder/cylinder.cpp      |   2 +-
 .../internalPoints/internalPoints.cpp         |  23 +-
 .../internalPoints/pointFlagKernels.hpp       |  21 +-
 .../structuredData/zAxis/array2D.hpp          |   4 +-
 .../triSurface/triangleFunctions.hpp          |   2 +-
 src/phasicFlow/types/basicTypes/math.hpp      | 191 ++++--
 src/setHelpers/setProperty.hpp                |   2 +-
 utilities/CMakeLists.txt                      |   2 +-
 utilities/Utilities/CMakeLists.txt            |   3 +-
 .../geometryPhasicFlow/stlWall}/stlFile.cpp   |  10 +-
 .../geometryPhasicFlow/stlWall}/stlFile.hpp   |   4 +-
 utilities/checkPhasicFlow/checkPhasicFlow.cpp |   9 +-
 .../geometryPhasicFlow/geometryPhasicFlow.cpp |  36 +-
 utilities/pFlowToVTK/pFlowToVTK.cpp           |  26 +-
 .../particlesPhasicFlow.cpp                   |  50 +-
 69 files changed, 1124 insertions(+), 490 deletions(-)
 create mode 100644 src/phasicFlow/containers/VectorHD/wordVectorHost.cpp
 create mode 100644 src/phasicFlow/containers/VectorHD/wordVectorHost.hpp
 rename {src/phasicFlow/triSurface => utilities/Utilities/geometryPhasicFlow/stlWall}/stlFile.cpp (97%)
 rename {src/phasicFlow/triSurface => utilities/Utilities/geometryPhasicFlow/stlWall}/stlFile.hpp (97%)

diff --git a/solvers/iterateGeometry/iterateGeometry.cpp b/solvers/iterateGeometry/iterateGeometry.cpp
index d541a5bf..1279d1c6 100755
--- a/solvers/iterateGeometry/iterateGeometry.cpp
+++ b/solvers/iterateGeometry/iterateGeometry.cpp
@@ -35,12 +35,10 @@ Licence:
 #include "commandLine.hpp"
 //#include "readControlDict.hpp"
 
-using namespace pFlow;
-
 int main( int argc, char* argv[] )
 {
 
-commandLine cmds(
+pFlow::commandLine cmds(
     "iterateGeometry",
     "Performs simulation without particles, only geometry is solved");
 
@@ -55,8 +53,8 @@ commandLine cmds(
    
 
 // this should be palced in each main 
-processors::initProcessors(argc, argv);
-initialize_pFlowProcessors();
+pFlow::processors::initProcessors(argc, argv);
+pFlow::initialize_pFlowProcessors();
 #include "initialize_Control.hpp"
 	
 	#include "setProperty.hpp"
@@ -71,7 +69,7 @@ initialize_pFlowProcessors();
 
 // this should be palced in each main 
 #include "finalize.hpp"
-processors::finalizeProcessors();
+pFlow::processors::finalizeProcessors();
 
 }
 
diff --git a/solvers/iterateSphereParticles/createDEMComponents.hpp b/solvers/iterateSphereParticles/createDEMComponents.hpp
index 85cb798a..e9975f44 100755
--- a/solvers/iterateSphereParticles/createDEMComponents.hpp
+++ b/solvers/iterateSphereParticles/createDEMComponents.hpp
@@ -20,6 +20,6 @@ Licence:
 
 //
 REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
-sphereParticles sphParticles(Control, proprties);
+pFlow::sphereParticles sphParticles(Control, proprties);
 
 WARNING<<"Particle insertion has not been set yet!"<<END_WARNING;
diff --git a/solvers/iterateSphereParticles/iterateSphereParticles.cpp b/solvers/iterateSphereParticles/iterateSphereParticles.cpp
index 30c181ba..1af6c7eb 100755
--- a/solvers/iterateSphereParticles/iterateSphereParticles.cpp
+++ b/solvers/iterateSphereParticles/iterateSphereParticles.cpp
@@ -36,7 +36,7 @@ Licence:
 #include "systemControl.hpp"
 #include "commandLine.hpp"
 
-using namespace pFlow;
+
 
 /**
  * DEM solver for simulating granular flow of cohesion-less particles.
@@ -47,7 +47,7 @@ using namespace pFlow;
 int main( int argc, char* argv[])
 {
 
-commandLine cmds(
+pFlow::commandLine cmds(
 		"sphereGranFlow",
 		"DEM solver for non-cohesive spherical particles with particle insertion "
 		"mechanism and moving geometry");
@@ -57,8 +57,8 @@ bool isCoupling = false;
 if(!cmds.parse(argc, argv)) return 0;
 
 // this should be palced in each main 
-processors::initProcessors(argc, argv); 
-initialize_pFlowProcessors();
+pFlow::processors::initProcessors(argc, argv); 
+pFlow::initialize_pFlowProcessors();
 #include "initialize_Control.hpp"
 	
 	#include "setProperty.hpp"
@@ -82,7 +82,7 @@ initialize_pFlowProcessors();
 
 // this should be palced in each main 
 #include "finalize.hpp"
-processors::finalizeProcessors();
+pFlow::processors::finalizeProcessors();
 
 
 }	
diff --git a/solvers/sphereGranFlow/createDEMComponents.hpp b/solvers/sphereGranFlow/createDEMComponents.hpp
index 57c5866d..ddbceafe 100755
--- a/solvers/sphereGranFlow/createDEMComponents.hpp
+++ b/solvers/sphereGranFlow/createDEMComponents.hpp
@@ -20,7 +20,7 @@ Licence:
 
 //
 REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
-sphereParticles sphParticles(Control, proprties);
+pFlow::sphereParticles sphParticles(Control, proprties);
 
 //
 REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
@@ -36,12 +36,12 @@ REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
 		sphParticles.shapes()
 	);*/
 
-auto sphInsertion =  sphereInsertion(
+auto sphInsertion =  pFlow::sphereInsertion(
 	sphParticles, 
 	sphParticles.spheres());
 
 REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<END_REPORT;
-auto interactionPtr = interaction::create(
+auto interactionPtr = pFlow::interaction::create(
 	Control,
 	sphParticles,
 	surfGeometry
diff --git a/solvers/sphereGranFlow/sphereGranFlow.cpp b/solvers/sphereGranFlow/sphereGranFlow.cpp
index 9def5903..d854b417 100755
--- a/solvers/sphereGranFlow/sphereGranFlow.cpp
+++ b/solvers/sphereGranFlow/sphereGranFlow.cpp
@@ -40,12 +40,8 @@ Licence:
 #include "interaction.hpp"
 #include "Insertions.hpp"
 
-
-
 //#include "readControlDict.hpp"
 
-using namespace pFlow;
-
 /**
  * DEM solver for simulating granular flow of cohesion-less particles.
  *
@@ -55,7 +51,7 @@ using namespace pFlow;
 int main( int argc, char* argv[])
 {
 
-commandLine cmds(
+pFlow::commandLine cmds(
 		"sphereGranFlow",
 		"DEM solver for non-cohesive spherical particles with particle insertion "
 		"mechanism and moving geometry");
@@ -65,13 +61,12 @@ bool isCoupling = false;
 if(!cmds.parse(argc, argv)) return 0;
 
 // this should be palced in each main 
-processors::initProcessors(argc, argv);
-initialize_pFlowProcessors();
+pFlow::processors::initProcessors(argc, argv);
+pFlow::initialize_pFlowProcessors();
 #include "initialize_Control.hpp"
 	 
 	#include "setProperty.hpp"
 	#include "setSurfaceGeometry.hpp"
-	
 
 	#include "createDEMComponents.hpp"
 	
@@ -117,7 +112,7 @@ initialize_pFlowProcessors();
 
 // this should be palced in each main 
 #include "finalize.hpp"
-processors::finalizeProcessors();
+pFlow::processors::finalizeProcessors();
 
 
 }	
diff --git a/src/Geometry/geometry/geometry.cpp b/src/Geometry/geometry/geometry.cpp
index 00855156..e9cbab45 100644
--- a/src/Geometry/geometry/geometry.cpp
+++ b/src/Geometry/geometry/geometry.cpp
@@ -249,7 +249,9 @@ pFlow::geometry::geometry
 	if( this->numSurfaces() != motionComponentName_.size() )
 	{
 		fatalErrorInFunction<<
-		"Number of surfaces is not equal to number of motion component names"<<endl;
+		"Number of surfaces ("<< this->numSurfaces() <<
+        ") is not equal to number of motion component names("<<
+        motionComponentName_.size()<<")"<<endl;
 		fatalExit;
 	}
 
diff --git a/src/Geometry/geometryMotion/geometryMotion.cpp b/src/Geometry/geometryMotion/geometryMotion.cpp
index 6420dadb..22124c47 100644
--- a/src/Geometry/geometryMotion/geometryMotion.cpp
+++ b/src/Geometry/geometryMotion/geometryMotion.cpp
@@ -1,4 +1,3 @@
-#include "geometryMotion.hpp"
 /*------------------------------- phasicFlow ---------------------------------
       O        C enter of
      O O       E ngineering and
@@ -70,6 +69,28 @@ bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
 	return true;
 }
 
+namespace pFlow::GMotion
+{
+    template<typename ModelInterface>
+    void moveGeometry(
+        real dt,
+        uint32 nPoints,
+        const ModelInterface& mModel,
+        const deviceViewType1D<uint32>& pointMIndexD,
+        const deviceViewType1D<realx3>& pointsD
+    )
+    {
+        Kokkos::parallel_for(
+         "geometryMotion<MotionModel>::movePoints",
+         deviceRPolicyStatic(0, nPoints),
+         LAMBDA_HD(uint32 i){
+             auto newPos = mModel.transferPoint(pointMIndexD[i], pointsD[i], dt);
+             pointsD[i] = newPos;
+         });
+
+        Kokkos::fence();
+    }
+}
 
 template<typename MotionModel>
  bool pFlow::geometryMotion<MotionModel>::moveGeometry()
@@ -84,8 +105,15 @@ template<typename MotionModel>
     auto& pointMIndexD= pointMotionIndex_.deviceViewAll();
     auto& pointsD = points().deviceViewAll();
     
+    pFlow::GMotion::moveGeometry(
+        dt, 
+        numPoints(), 
+        motionModel_.getModelInterface(iter, t, dt),
+        pointMotionIndex_.deviceViewAll(),
+        points().deviceViewAll()
+    );
 
-    Kokkos::parallel_for(
+    /*Kokkos::parallel_for(
          "geometryMotion<MotionModel>::movePoints",
          deviceRPolicyStatic(0, numPoints()),
          LAMBDA_HD(uint32 i){
@@ -93,7 +121,7 @@ template<typename MotionModel>
              pointsD[i] = newPos;
          });
 
-    Kokkos::fence();
+    Kokkos::fence();*/
 
     // move the motion components
     motionModel_.move(iter, t,dt);
diff --git a/src/Interaction/CMakeLists.txt b/src/Interaction/CMakeLists.txt
index c0265808..96197f6a 100644
--- a/src/Interaction/CMakeLists.txt
+++ b/src/Interaction/CMakeLists.txt
@@ -7,8 +7,8 @@ contactSearch/methods/cellBased/NBS/NBS.cpp
 contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
 
 contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.cpp
-contactSearch/boundaries/twoPartContactSearch/twoPartContactSearchKernels.cpp
-contactSearch/boundaries/twoPartContactSearch/twoPartContactSearch.cpp
+#contactSearch/boundaries/twoPartContactSearch/twoPartContactSearchKernels.cpp
+#contactSearch/boundaries/twoPartContactSearch/twoPartContactSearch.cpp
 contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp
 contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearch.cpp
 contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp
diff --git a/src/Interaction/contactLists/sortedContactList.hpp b/src/Interaction/contactLists/sortedContactList.hpp
index 73f0ce4b..4f41caac 100644
--- a/src/Interaction/contactLists/sortedContactList.hpp
+++ b/src/Interaction/contactLists/sortedContactList.hpp
@@ -138,7 +138,7 @@ public:
 				start,
 				end,
 				newPair);
-				idx0!=-1)
+				idx0!=static_cast<uint32>(-1))
 		{
 			values_[idx] = values0_[idx0];
 		}
@@ -147,7 +147,7 @@ public:
 			start, 
 			end,
 			newPair);
-			idx0!=-1)
+			idx0!= static_cast<uint32>(-1) )
 		{
 			values_[idx] = values0_[idx0];
 			
diff --git a/src/Interaction/contactLists/unsortedContactList.hpp b/src/Interaction/contactLists/unsortedContactList.hpp
index 8a61b2ce..7aa98237 100644
--- a/src/Interaction/contactLists/unsortedContactList.hpp
+++ b/src/Interaction/contactLists/unsortedContactList.hpp
@@ -124,7 +124,7 @@ public:
 	INLINE_FUNCTION_HD
 	bool getValue(const PairType& p, ValueType& val)const
 	{
-		if(auto idx = this->find(p); idx!=-1)
+		if(auto idx = this->find(p); idx!=static_cast<uint32>(-1))
 		{
 			val = getValue(idx); 
 			return true;
@@ -141,7 +141,7 @@ public:
 	INLINE_FUNCTION_HD
 	bool setValue(const PairType& p, const ValueType& val)const
 	{
-		if(uint32 idx = this->find(p); idx!=-1)
+		if(uint32 idx = this->find(p); idx!=static_cast<uint32>(-1))
 		{
 			setValue(idx, val);
 			return true;;
@@ -156,7 +156,7 @@ public:
 		{
 			if( uint32 idx0 = 
 					container0_.find(this->getPair(idx));
-					idx0!=-1 )
+					idx0!= static_cast<uint32>(-1) )
 			{
 				values_[idx] = values0_[idx0];
 			}
diff --git a/src/Interaction/contactLists/unsortedPairs.hpp b/src/Interaction/contactLists/unsortedPairs.hpp
index fd4f117d..fb622bdb 100644
--- a/src/Interaction/contactLists/unsortedPairs.hpp
+++ b/src/Interaction/contactLists/unsortedPairs.hpp
@@ -105,7 +105,7 @@ public:
 	uint32 insert(idType i, idType j)const
 	{
 		if(auto insertResult = container_.insert(PairType(i,j)); insertResult.failed())
-			return -1;
+			return static_cast<uint32>(-1);
 		else
 			return insertResult.index();
 			
@@ -115,7 +115,7 @@ public:
 	uint32 insert(const PairType& p)const
 	{
 		if(auto insertResult = container_.insert(p); insertResult.failed())
-			return -1;
+			return static_cast<uint32>(-1);
 		else
 			return insertResult.index();
 
@@ -154,7 +154,7 @@ public:
 			idx != Kokkos::UnorderedMapInvalidIndex )
 			return idx;
 		else
-			return -1;
+			return static_cast<uint32>(-1);
 	}
 
 	INLINE_FUNCTION_HD
diff --git a/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp b/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp
index ea266e78..a0d47e78 100644
--- a/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp
+++ b/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp
@@ -80,14 +80,14 @@ pFlow::uint32 pFlow::pweBndryContactSearchKernels::broadSearchPP
 				if(!searchCells.inCellRange(ind))continue;
 
 				uint32 thisI = head(ind.x(),ind.y(),ind.z());
-				while (thisI!=-1)
+				while (thisI!=static_cast<uint32>(-1))
 				{				
 					
 					auto d_n = sizeRatio*diams[thisI];
 					
                     // first item is for this boundary and second itme, for mirror 
 					if(sphereSphereCheckB(p_m, points[thisI], d_m, d_n)&&
-                       ppPairs.insert(thisI,mrrI) == -1)
+                       ppPairs.insert(thisI,mrrI) == static_cast<uint32>(-1))
 					{
 						getFullUpdate++;	
 					}
diff --git a/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp b/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp
index 99bbb8ac..acfad226 100644
--- a/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp
+++ b/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp
@@ -113,7 +113,10 @@ pFlow::uint32 pFlow::wallBoundaryContactSearch::findPairsElementRangeCount
 
     uint32 nNotInserted = 0;
     uint32 nThis = pPoints.size();
-    
+    const auto& numElements = numElements_;
+    const auto& elementBox = elementBox_;
+    const auto& validBox = validBox_;
+
     Kokkos::parallel_reduce(
         "pFlow::wallBoundaryContactSearch::findPairsElementRangeCount",
         deviceRPolicyDynamic(0,nThis),
@@ -123,11 +126,11 @@ pFlow::uint32 pFlow::wallBoundaryContactSearch::findPairsElementRangeCount
             int32x3 ind;
             if( searchCells.pointIndexInDomain(p, ind) )
             {
-                for(uint32 nTri=0; nTri<numElements_; nTri++)
+                for(uint32 nTri=0; nTri<numElements; nTri++)
                 {
-                    if( validBox_[nTri]== 0)continue;
-                    if( elementBox_[nTri].isInside(ind)&& 
-                        pairs.insert(i,nTri+baseTriIndex) == -1)
+                    if( validBox[nTri]== 0)continue;
+                    if( elementBox[nTri].isInside(ind)&& 
+                        pairs.insert(i,nTri+baseTriIndex) == static_cast<uint32>(-1))
                     {
                         notInsertedUpdate++;    
                     }
diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLoop.hpp b/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLoop.hpp
index ab3eeea7..ecbc63e5 100644
--- a/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLoop.hpp
+++ b/src/Interaction/contactSearch/methods/cellBased/NBS/NBSLoop.hpp
@@ -43,7 +43,7 @@ while( m !=  mapperNBS::NoPos)
 			auto lm = m;
 
 			if(lm>ln) Swap(lm,ln);
-			if( pairs.insert(lm,ln) == -1)
+			if( pairs.insert(lm,ln) == static_cast<uint32>(-1))
 			{
 				getFullUpdate++;
 			}
@@ -86,7 +86,7 @@ while( m !=  mapperNBS::NoPos)
 					auto ln = n;
 					auto lm = m;
 					if(lm>ln) Swap(lm,ln);
-					if( pairs.insert(lm,ln) == -1)
+					if( pairs.insert(lm,ln) == static_cast<uint32>(-1))
 					{
 						getFullUpdate++;
 					}
diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp b/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
index 00932281..e1b6e3a7 100644
--- a/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
+++ b/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
@@ -85,21 +85,26 @@ bool pFlow::cellsWallLevel0::broadSearch
 bool pFlow::cellsWallLevel0::build(const cells & searchBox)
 {
 	
+    const auto& points = points_;
+    const auto& vertices = vertices_;
+    const auto& elementBox = elementBox_;
+    const auto cellExtent = cellExtent_;
+
 	Kokkos::parallel_for(
 		"pFlow::cellsWallLevel0::build",
 		deviceRPolicyStatic(0,numElements_),
-		CLASS_LAMBDA_HD(uint32 i)
+		LAMBDA_HD(uint32 i)
 		{
-			auto v = vertices_[i];
-			auto p1 = points_[v.x()];
-			auto p2 = points_[v.y()];
-			auto p3 = points_[v.z()];
+			auto v = vertices[i];
+			auto p1 = points[v.x()];
+			auto p2 = points[v.y()];
+			auto p3 = points[v.z()];
 			
 			realx3 minP;
 			realx3 maxP;
 
-			searchBox.extendBox(p1, p2, p3, cellExtent_, minP, maxP);
-			elementBox_[i] = iBoxType(searchBox.pointIndex(minP), searchBox.pointIndex(maxP));
+			searchBox.extendBox(p1, p2, p3, cellExtent, minP, maxP);
+			elementBox[i] = iBoxType(searchBox.pointIndex(minP), searchBox.pointIndex(maxP));
 		});
 	Kokkos::fence();
 
@@ -153,7 +158,12 @@ pFlow::int32 pFlow::cellsWallLevel0::findPairsElementRangeCount
 {
 	uint32 getFull =0;
 			
-	
+	const auto& elementBox = elementBox_;
+    const auto& normals = normals_;
+    const auto& points = points_;
+    const auto& vertices = vertices_;
+    const auto cellExtent = cellExtent_;
+
 	Kokkos::parallel_reduce(
 		"pFlow::cellsWallLevel0::findPairsElementRangeCount",
 		tpPWContactSearch(numElements_, Kokkos::AUTO),
@@ -163,10 +173,10 @@ pFlow::int32 pFlow::cellsWallLevel0::findPairsElementRangeCount
 			
 			const uint32 iTri = teamMember.league_rank();
 
-			const auto triBox = elementBox_[iTri];
+			const auto triBox = elementBox[iTri];
 			const auto triPlane = infinitePlane(
-				normals_[iTri], 
-				points_[vertices_[iTri].x()]);
+				normals[iTri], 
+				points[vertices[iTri].x()]);
 
 			uint32 getFull2 = 0;
 
@@ -186,11 +196,12 @@ pFlow::int32 pFlow::cellsWallLevel0::findPairsElementRangeCount
 						while( n != particleMap.NoPos)
 						{
 							// id is wall id the pair is (particle id, wall id)
-							if( abs(triPlane.pointFromPlane(pPoints[n]))< pDiams[n]*sizeRatio*cellExtent_)
+							if( abs(triPlane.pointFromPlane(pPoints[n]))< pDiams[n]*sizeRatio*cellExtent)
 							{
 								if( pairs.insert(
 									static_cast<csIdType>(n),
-									static_cast<csIdType>(iTri) ) == -1 )
+									static_cast<csIdType>(iTri) ) == static_cast<csIdType>(-1) 
+                                )
 									innerUpdate++;
 							}
 							n = particleMap.next(n);
diff --git a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteractionKernels.hpp b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteractionKernels.hpp
index c581d598..751c0fec 100644
--- a/src/Interaction/sphereInteraction/sphereInteraction/sphereInteractionKernels.hpp
+++ b/src/Interaction/sphereInteraction/sphereInteraction/sphereInteractionKernels.hpp
@@ -261,7 +261,7 @@ struct pwInteractionFunctor
 			int32 propId_i = propId_[i];
 			int32 wPropId_j = wPropId_[tj];
 
-			realx3 FCn, FCt, Mri, Mrj, Mij, Mji;
+			realx3 FCn, FCt, Mri, Mrj, Mij;
 			//output<< "before "<<history.overlap_t_<<endl;
 			// calculates contact force 
 			forceModel_.contactForce(
diff --git a/src/MotionModel/entities/stationary/stationary.hpp b/src/MotionModel/entities/stationary/stationary.hpp
index 723ae31f..2970e197 100644
--- a/src/MotionModel/entities/stationary/stationary.hpp
+++ b/src/MotionModel/entities/stationary/stationary.hpp
@@ -62,7 +62,7 @@ public:
 	INLINE_FUNCTION_HD
 	realx3 linVelocityPoint(const realx3 &)const
 	{
-		return zero3;
+		return realx3(0);
 	}
 
 	INLINE_FUNCTION_HD
diff --git a/src/MotionModel/rotatingAxisMotion/rotatingAxisMotion.hpp b/src/MotionModel/rotatingAxisMotion/rotatingAxisMotion.hpp
index a46350bd..5a191f11 100644
--- a/src/MotionModel/rotatingAxisMotion/rotatingAxisMotion.hpp
+++ b/src/MotionModel/rotatingAxisMotion/rotatingAxisMotion.hpp
@@ -76,8 +76,6 @@ protected:
 	{
 		return true;
 	}
-
-	void impl_setTime(uint32 iter, real t, real dt)const;
 	
 public:
 
@@ -90,6 +88,7 @@ public:
 		const dictionary& dict, 
 		repository* owner);
 
+    using fileDictionary::write;
 
 	bool write(iOstream& os, const IOPattern& iop)const override;
 
@@ -98,6 +97,9 @@ public:
     {
       return rotatingAxis({0,0,0}, {1,0,0}, 0.0);
     }
+
+    // TODO: make this method protected
+    void impl_setTime(uint32 iter, real t, real dt)const;
 };
 
 } // pFlow
diff --git a/src/MotionModel/stationaryWall/stationaryWall.hpp b/src/MotionModel/stationaryWall/stationaryWall.hpp
index 514cba1d..5e0167f0 100644
--- a/src/MotionModel/stationaryWall/stationaryWall.hpp
+++ b/src/MotionModel/stationaryWall/stationaryWall.hpp
@@ -84,6 +84,8 @@ public:
 		repository* owner);
 
 
+    using fileDictionary::write;
+
 	bool write(iOstream& os, const IOPattern& iop)const override;
 
 	static
diff --git a/src/MotionModel/vibratingMotion/vibratingMotion.hpp b/src/MotionModel/vibratingMotion/vibratingMotion.hpp
index 8c64e3a5..79dbc6f4 100644
--- a/src/MotionModel/vibratingMotion/vibratingMotion.hpp
+++ b/src/MotionModel/vibratingMotion/vibratingMotion.hpp
@@ -82,7 +82,6 @@ protected:
 		return true;
 	}
 
-	void impl_setTime(uint32 iter, real t, real dt)const;
 
 public:
 	
@@ -99,15 +98,20 @@ public:
 	/// Destructor 
 	~vibratingMotion()override = default;
 
+    using fileDictionary::write;
     
 	bool write(iOstream& os, const IOPattern& iop)const override;
 
+    
+
 	static
     auto noneComponent()
     {
       return vibrating();
     }
 
+    // TODO: make this protected
+    void impl_setTime(uint32 iter, real t, real dt)const;
 };
 
 } // pFlow
diff --git a/src/Particles/Insertion/collisionCheck/collisionCheck.cpp b/src/Particles/Insertion/collisionCheck/collisionCheck.cpp
index 5cb7b470..ba8d4c9d 100644
--- a/src/Particles/Insertion/collisionCheck/collisionCheck.cpp
+++ b/src/Particles/Insertion/collisionCheck/collisionCheck.cpp
@@ -60,7 +60,7 @@ pFlow::collisionCheck::checkPoint(const realx3& p, const real d) const
             {
                 uint32 n = head_(i, j, k);
     
-                while( n != -1)
+                while( n != static_cast<uint32>(-1))
                 {
                     if( ((position_[n]-p).length() - 0.5*(diameters_[n]+d )) <= 0.0 )
                     {
@@ -85,7 +85,7 @@ pFlow::collisionCheck::mapLastAddedParticle()
         "size mismatch of next and position"<<endl;
         return false;
     }
-    next_.push_back(-1);
+    next_.push_back(static_cast<uint32>(-1));
     const auto& p = position_[n];
 
     if(!searchBox_.isInside(p))
diff --git a/src/Particles/Insertion/insertion/insertion.cpp b/src/Particles/Insertion/insertion/insertion.cpp
index d8d0881d..bb551a47 100644
--- a/src/Particles/Insertion/insertion/insertion.cpp
+++ b/src/Particles/Insertion/insertion/insertion.cpp
@@ -74,9 +74,6 @@ pFlow::insertion::pStruct() const
 	return particles_.pStruct();
 }
 
-
-
-
 bool
 pFlow::insertion::readInsertionDict()
 {
diff --git a/src/Particles/Insertion/insertion/insertion.hpp b/src/Particles/Insertion/insertion/insertion.hpp
index f8419aa3..f39f5c35 100644
--- a/src/Particles/Insertion/insertion/insertion.hpp
+++ b/src/Particles/Insertion/insertion/insertion.hpp
@@ -32,7 +32,9 @@ class pointStructure;
 /**
  * Base class for particle insertion
  */
-class insertion : public fileDictionary
+class insertion 
+: 
+    public fileDictionary
 {
 private:
 
@@ -118,6 +120,8 @@ public:
 	/*/// read from iIstream
 	virtual bool read(iIstream& is) = 0;*/
 
+    using fileDictionary::write;
+
 	/// write to iOstream
 	bool write(iOstream& os, const IOPattern& iop)const override ;
 };
diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp b/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp
index 0f78e498..d71e751a 100644
--- a/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp
+++ b/src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp
@@ -82,7 +82,7 @@ private:
 	Timer  				   fieldUpdateTimer_;
 
 private:
-	bool initializeParticles();
+	
 
 	bool getParticlesInfoFromShape(
 	  const wordVector& shapeNames,
@@ -119,11 +119,14 @@ public:
 	 */
 	/*bool insertParticles
 	(
-	    const realx3Vector& position,
+    const realx3Vector& position,
 	    const wordVector&  shapes,
 	    const setFieldList& setField
 	) override ;*/
 
+    // TODO: make this method private later 
+    bool initializeParticles();
+
 	/// const reference to shapes object
 	const auto& spheres() const
 	{
diff --git a/src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.cpp b/src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.cpp
index 1916f430..b8a55aba 100644
--- a/src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.cpp
+++ b/src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.cpp
@@ -15,8 +15,11 @@ pFlow::regularParticleIdHandler::regularParticleIdHandler
 pFlow::Pair<pFlow::uint32, pFlow::uint32> 
     pFlow::regularParticleIdHandler::getIdRange(uint32 nNewParticles)
 {
+    
+    if(nNewParticles==0) return {0,0};
+
     uint32 startId;
-    if(maxId_==-1)
+    if(maxId_== static_cast<uint32>(-1))
     {
         startId = 0;
     }
@@ -37,7 +40,7 @@ bool pFlow::regularParticleIdHandler::initialIdCheck()
     uint32 maxId = max( *this );
     
     /// particles should get ids from 0 to size-1
-    if(maxId == -1)
+    if(maxId == static_cast<uint32>(-1))
     {
         fillSequence(*this,0u);
         maxId_ = size()-1;
diff --git a/src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.hpp b/src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.hpp
index a237e0da..3e4f8b41 100644
--- a/src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.hpp
+++ b/src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.hpp
@@ -31,7 +31,7 @@ class regularParticleIdHandler
 {
 private:
     
-    uint32 maxId_ = -1;
+    uint32 maxId_ = static_cast<uint32>(-1);
 
   bool initialIdCheck()override; 
 public:
diff --git a/src/Particles/particles/shape/baseShapeNames.hpp b/src/Particles/particles/shape/baseShapeNames.hpp
index bc513b06..7e0a00c3 100644
--- a/src/Particles/particles/shape/baseShapeNames.hpp
+++ b/src/Particles/particles/shape/baseShapeNames.hpp
@@ -132,7 +132,7 @@ public:
 				return true;
 			}
 		}
-		idx = -1;
+		idx = static_cast<uint32>(-1);
 		return false;
 	}
 
@@ -144,6 +144,8 @@ public:
 
 	// - IO
 
+    using fileDictionary::write;
+
 	bool write(iOstream& os)const override;
 
 };
diff --git a/src/phasicFlow/CMakeLists.txt b/src/phasicFlow/CMakeLists.txt
index 408fc042..7de9ceda 100644
--- a/src/phasicFlow/CMakeLists.txt
+++ b/src/phasicFlow/CMakeLists.txt
@@ -59,6 +59,7 @@ Timer/Timers.cpp
 
 
 containers/Vector/Vectors.cpp
+containers/VectorHD/wordVectorHost.cpp
 containers/VectorHD/VectorSingles.cpp
 containers/Field/Fields.cpp
 containers/symArrayHD/symArrays.cpp
diff --git a/src/phasicFlow/Kokkos/ViewAlgorithms.hpp b/src/phasicFlow/Kokkos/ViewAlgorithms.hpp
index db071f0b..28e32abd 100644
--- a/src/phasicFlow/Kokkos/ViewAlgorithms.hpp
+++ b/src/phasicFlow/Kokkos/ViewAlgorithms.hpp
@@ -421,10 +421,10 @@ binarySearch(
 )
 {
 	if (end <= start)
-		return -1;
+		return static_cast<uint32>(-1);
 
 	if (auto res = binarySearch_(view.data() + start, end - start, val);
-	    res != -1)
+	    res != static_cast<uint32>(-1))
 	{
 		return res + start;
 	}
diff --git a/src/phasicFlow/commandLine/CLI/TypeTools.hpp b/src/phasicFlow/commandLine/CLI/TypeTools.hpp
index 2b87ec60..667d3eac 100644
--- a/src/phasicFlow/commandLine/CLI/TypeTools.hpp
+++ b/src/phasicFlow/commandLine/CLI/TypeTools.hpp
@@ -146,11 +146,11 @@ template <typename T, typename C> class is_direct_constructible {
     static auto test(int, std::true_type) -> decltype(
 // NVCC warns about narrowing conversions here
 #ifdef __CUDACC__
-#pragma diag_suppress 2361
+#pragma nv_diag_suppress 2361
 #endif
         TT { std::declval<CC>() }
 #ifdef __CUDACC__
-#pragma diag_default 2361
+#pragma nv_diag_default 2361
 #endif
         ,
         std::is_move_assignable<TT>());
diff --git a/src/phasicFlow/containers/Field/Field.hpp b/src/phasicFlow/containers/Field/Field.hpp
index 887b994d..738811bc 100644
--- a/src/phasicFlow/containers/Field/Field.hpp
+++ b/src/phasicFlow/containers/Field/Field.hpp
@@ -24,6 +24,7 @@ Licence:
 
 #include "types.hpp"
 #include "VectorSingle.hpp"
+#include "wordVectorHost.hpp"
 #include "Vector.hpp"
 #include "streams.hpp"
 
diff --git a/src/phasicFlow/containers/Field/Fields.cpp b/src/phasicFlow/containers/Field/Fields.cpp
index b0cbad75..e76daa11 100644
--- a/src/phasicFlow/containers/Field/Fields.cpp
+++ b/src/phasicFlow/containers/Field/Fields.cpp
@@ -31,5 +31,4 @@ template class pFlow::Field<pFlow::real>;
 
 template class pFlow::Field<pFlow::realx3>;
 
-
-
+template class pFlow::Field<pFlow::word, pFlow::HostSpace>;
diff --git a/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp b/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp
index 77ef5a99..861239de 100644
--- a/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp
+++ b/src/phasicFlow/containers/List/ListPtr/ListPtrI.hpp
@@ -39,6 +39,7 @@ inline bool pFlow::ListPtr<T>::copy(const ListPtrType& src)
 }
 
 template<typename T>
+inline
 T* pFlow::ListPtr<T>::ptr(size_t i) 
 {
 
@@ -51,6 +52,7 @@ T* pFlow::ListPtr<T>::ptr(size_t i)
 }
 
 template<typename T>
+inline
 const T* pFlow::ListPtr<T>::ptr
 (
 	size_t i
@@ -66,6 +68,7 @@ const T* pFlow::ListPtr<T>::ptr
 }
 
 template<typename T>
+inline
 auto pFlow::ListPtr<T>::pos
 (
 	size_t i
@@ -84,6 +87,7 @@ auto pFlow::ListPtr<T>::pos
 }
 
 template<typename T>
+inline
 auto pFlow::ListPtr<T>::pos
 (
 	size_t i
@@ -102,6 +106,7 @@ auto pFlow::ListPtr<T>::pos
 }
 
 template<typename T>
+inline
 pFlow::ListPtr<T>::ListPtr
 (
 	const ListPtrType& src
@@ -119,6 +124,7 @@ pFlow::ListPtr<T>::ListPtr
 
 
 template<typename T>
+inline
 pFlow::ListPtr<T>& pFlow::ListPtr<T>::operator=
 (
 	const ListPtrType& rhs
@@ -144,6 +150,7 @@ pFlow::ListPtr<T>& pFlow::ListPtr<T>::operator=
 }
 
 template<typename T>
+inline
 pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
 (
 	size_t i, T* ptr
@@ -155,6 +162,7 @@ pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
 
 
 template<typename T>
+inline
 pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
 (
 	size_t i,
@@ -179,6 +187,7 @@ pFlow::uniquePtr<T> pFlow::ListPtr<T>::set
 
 template<typename T>
 template<typename... Args>
+inline
 pFlow::uniquePtr<T> pFlow::ListPtr<T>::setSafe
 (
 	size_t i,
@@ -190,6 +199,7 @@ pFlow::uniquePtr<T> pFlow::ListPtr<T>::setSafe
 }
 
 template<typename T>
+inline
 void pFlow::ListPtr<T>::push_back
 (
 	T* ptr
@@ -199,6 +209,7 @@ void pFlow::ListPtr<T>::push_back
 }
 
 template<typename T>
+inline
 void pFlow::ListPtr<T>::push_back(uniquePtr<T>&& ptr)
 {
 	list_.push_back( ptr.release() );
@@ -206,13 +217,15 @@ void pFlow::ListPtr<T>::push_back(uniquePtr<T>&& ptr)
 
 template<typename T>
 template<typename... Args>
+inline
 void pFlow::ListPtr<T>::push_backSafe(Args&&... args)
 {
-	auto ptr=makeUnique<T>(std::forward<Args>(args)...) ;
+	uniquePtr<T> ptr = makeUnique<T>(std::forward<Args>(args)...) ;
 	push_back(ptr);
 }
 
 template<typename T>
+inline
 T& pFlow::ListPtr<T>::operator[]
 (
 	size_t i
@@ -231,6 +244,7 @@ T& pFlow::ListPtr<T>::operator[]
 }
 
 template<typename T>
+inline
 const T& pFlow::ListPtr<T>::operator[]
 (
 	size_t i
@@ -248,18 +262,21 @@ const T& pFlow::ListPtr<T>::operator[]
 }
 
 template<typename T>
+inline
 size_t pFlow::ListPtr<T>::size()const
 {
 	return list_.size();
 }
 
 template<typename T>
+inline
 auto pFlow::ListPtr<T>::empty() const
 {
 	return list_.emtpy();
 }
 
 template<typename T>
+inline
 pFlow::uniquePtr<T> pFlow::ListPtr<T>::release
 (
 	size_t i
@@ -273,15 +290,14 @@ pFlow::uniquePtr<T> pFlow::ListPtr<T>::release
 
 
 template<typename T>
+inline
 void pFlow::ListPtr<T>::clear()
 {
 	
-	int i =0;
 	for( auto iter = list_.begin(); iter != list_.end(); ++iter )
 	{
 		if(*iter != nullptr)
-		{
-	
+		{	
 			delete *iter;
 			*iter = nullptr;
 		}	
@@ -291,6 +307,7 @@ void pFlow::ListPtr<T>::clear()
 }
 
 template<typename T>
+inline
 void pFlow::ListPtr<T>::clear
 (
 	size_t i
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingle.cpp b/src/phasicFlow/containers/VectorHD/VectorSingle.cpp
index dd093270..0e6f02b2 100644
--- a/src/phasicFlow/containers/VectorHD/VectorSingle.cpp
+++ b/src/phasicFlow/containers/VectorHD/VectorSingle.cpp
@@ -55,30 +55,14 @@ void pFlow::VectorSingle<T,MemorySpace>::changeCapacity
     bool withInit
 )
 {
-    if constexpr( isTriviallyCopyable_ )
+    
+    if(withInit)
     {
-        if(withInit)
-        {
-            resizeInit(view_, actualCap);
-        }
-        else
-        {
-            resizeNoInit(view_, actualCap);
-        }
-    }
-    else if constexpr( isHostAccessible_ )
-    {
-        viewType newView(view_.label(), actualCap);
-
-        for(auto i=0u; i<min(size_,actualCap); i++)
-        {
-            newView(i) = view_(i);
-        }
-        view_ = newView;
+        resizeInit(view_, actualCap);
     }
     else
     {
-        static_assert("changeCapacity is not a valid operation for non-trivially-copyable data type on device memory");
+        resizeNoInit(view_, actualCap);
     }
     
 }
@@ -91,17 +75,8 @@ pFlow::uint32 pFlow::VectorSingle<T,MemorySpace>::reallocateCapacitySize
     uint32 s
 )
 {
-    if constexpr (isTriviallyCopyable_)
-    {
-        reallocNoInit(view_, cap);
-    }
-    else
-    {
-        viewType newView(view_.label(), cap);
-        view_ = newView;
-    }
-    
-    return setSize(s);
+    reallocNoInit(view_, cap);
+    return setSize(s);   
 }
 
 template<typename T, typename MemorySpace>
@@ -209,21 +184,7 @@ pFlow::VectorSingle<T,MemorySpace>::VectorSingle
 :
     VectorSingle(name, src.capacity(), src.size(), RESERVE())
 {
-    if constexpr(isTriviallyCopyable_)
-    {
-        copy(deviceView(), src.deviceView());
-    }
-    else if constexpr( isHostAccessible_)
-    {
-        for(auto i=0u; i<src.size(); i++)
-        {
-            view_[i] = src.view_[i];
-        }
-    }
-    else
-    {
-        static_assert("This constructor is not valid for non-trivially copyable data type on device memory");
-    }
+    copy(deviceView(), src.deviceView());
 }
 
 template<typename T, typename MemorySpace>
@@ -235,21 +196,7 @@ pFlow::VectorSingle<T,MemorySpace>::VectorSingle
 :
     VectorSingle(name, src.size(), src.size(), RESERVE())
 {
-    if constexpr(isTriviallyCopyable_)
-    {
-        copy(deviceView(), src);
-    }
-    else if constexpr( isHostAccessible_)
-    {
-        for(auto i=0u; i<size(); i++)
-        {
-            view_[i] = src[i];
-        }
-    }
-    else
-    {
-        static_assert("This constructor is not valid for non-trivially copyable data type on device memory");
-    }
+    copy(deviceView(), src);
 }
 
 template<typename T, typename MemorySpace>
@@ -325,18 +272,7 @@ INLINE_FUNCTION_H
 auto pFlow::VectorSingle<T,MemorySpace>::hostViewAll()const
 {
     auto hView = Kokkos::create_mirror_view(view_);
-    if constexpr(isTriviallyCopyable_)
-    {
-        copy(hView, view_);
-    }
-    else if constexpr( isHostAccessible_ )
-    {
-        // nothing to be done, since it is already a host memory
-    }
-    else
-    {
-        static_assert("hostViewAll is not valid for non-trivially copyable data type on device memory");
-    }
+    copy(hView, view_);
     return hView;
 }
 
@@ -345,19 +281,7 @@ INLINE_FUNCTION_H
 auto pFlow::VectorSingle<T,MemorySpace>::hostView()const
 {
     auto hView = Kokkos::create_mirror_view(deviceView());
-    if constexpr(isTriviallyCopyable_)
-    {
-        copy(hView, deviceView());
-    }
-    else if constexpr( isHostAccessible_ )
-    {
-        // nothing to be done, since it is already a host memory
-    }
-    else
-    {
-        static_assert("hostView is not valid for non-trivially copyable data type on device memory");
-    }
-    
+    copy(hView, deviceView());
     return hView;
 }
 
@@ -497,23 +421,9 @@ void pFlow::VectorSingle<T,MemorySpace>::assign
         changeSize(srcSize);
     }
     
-    if constexpr( isTriviallyCopyable_ )
-    {
-        // - unmanaged view in the host
-        hostViewType1D<const T> temp(src.data(), srcSize );
-        copy(deviceView(), temp);
-    }
-    else if constexpr( isHostAccessible_)
-    {
-        for(auto i=0u; i<srcSize; i++)
-        {
-            view_[i] = src[i];
-        }   
-    }
-    else
-    {
-        static_assert("Not a valid operation for this data type on memory device");
-    }
+    // - unmanaged view in the host
+    hostViewType1D<const T> temp(src.data(), srcSize );
+    copy(deviceView(), temp);
     
 }
 
@@ -543,21 +453,8 @@ void pFlow::VectorSingle<T,MemorySpace>::assignFromHost(const VectorTypeHost& sr
         changeSize(srcSize);
     }
 	
-    if constexpr(isTriviallyCopyable_)
-    {
-        copy(deviceView(), src.hostView());
-    }
-    else if constexpr( isHostAccessible_)
-    {
-        for(auto i=0u; i<src.size(); i++)
-        {
-            view_[i] = src.view_[i];
-        }
-    }
-    else
-    {
-        static_assert("Not a valid operation for this data type on device memory");
-    }
+    copy(deviceView(), src.hostView());
+    
 }
 
 template<typename T, typename MemorySpace>
@@ -580,25 +477,30 @@ void pFlow::VectorSingle<T,MemorySpace>::assign
         changeSize(srcSize);
     }
     
-
-    if constexpr(isTriviallyCopyable_)
-    {
-        copy(deviceView(), src.deviceView());
-    }
-    else if constexpr( isHostAccessible_)
-    {
-        for(auto i=0u; i<src.size(); i++)
-        {
-            view_[i] = src.view_[i];
-        }
-    }
-    else
-    {
-        static_assert("Not a valid operation for this data type on device memory");
-    }
+    copy(deviceView(), src.deviceView());
 }
 
 
+template<typename T, typename MemorySpace>
+template<typename MSpace>
+INLINE_FUNCTION_H
+void pFlow::VectorSingle<T,MemorySpace>::assignFromDevice(
+    const VectorSingle<T, MSpace>& src, 
+    bool srcCapacity
+)
+{
+    uint32 srcSize = src.size();
+    uint32 srcCap = src.capacity();
+
+    if(srcCapacity && srcCap != capacity()){
+        reallocateCapacitySize(srcCap, srcSize);
+    }
+    else {
+        changeSize(srcSize);
+    }
+    copy(deviceView(), src.deviceView());
+}
+
 template <typename T, typename MemorySpace>
 INLINE_FUNCTION_H
 void pFlow::VectorSingle<T, MemorySpace>::append(const ViewType1D<T,MemorySpace>& appVec)
@@ -615,21 +517,8 @@ void pFlow::VectorSingle<T, MemorySpace>::append(const ViewType1D<T,MemorySpace>
 		view_,
 		Kokkos::make_pair<uint32>(oldS, newSize));
 	
-    if constexpr( isTriviallyCopyable_)
-    {
-        copy(appendView, appVec);
-    }
-    else if constexpr( isHostAccessible_)
-    {
-        for(auto i=0u; i<appVec.size(); i++)
-        {
-            appendView[i] = appVec[i];
-        }
-    }
-    else
-    {
-        static_assert("not a valid operation for this data type on device memory");
-    }
+    copy(appendView, appVec);
+    
 }
 
 template <typename T, typename MemorySpace>
@@ -650,22 +539,7 @@ void pFlow::VectorSingle<T, MemorySpace>::append
     hostViewType1D<const T> temp(appVec.data(), srcSize );
     auto dest = Kokkos::subview(view_, Kokkos::make_pair<uint32>(oldSize,newSize));
 
-    if constexpr( isTriviallyCopyable_)
-    {
-        copy(dest, temp);
-    }
-    else if constexpr( isHostAccessible_)
-    {
-        for(auto i=0u; i<appVec.size(); i++)
-        {
-            dest[i] = appVec[i];
-        }
-    }
-    else
-    {
-        static_assert("not a valid operation for this data type on device memory");
-    }
-    
+    copy(dest, temp);    
 }
 
 template <typename T, typename MemorySpace>
@@ -687,21 +561,8 @@ void pFlow::VectorSingle<T, MemorySpace>::append
 		view_,
 		Kokkos::make_pair<uint32>(oldS, newSize));
 	
-    if constexpr( isTriviallyCopyable_)
-    {
-        copy(appendView, appVec.deviceView());
-    }
-    else if constexpr( isHostAccessible_)
-    {
-        for(auto i=0u; i<appVec.size(); i++)
-        {
-            appendView[i] = appVec.view_[i];
-        }
-    }
-    else
-    {
-        static_assert("not a valid operation for this data type on device memory");
-    }
+    copy(appendView, appVec.deviceView());
+    
 }
 
 template <typename T, typename MemorySpace>
@@ -936,22 +797,8 @@ bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(const uint32IndexContainer
 
     setSize(newSize);
 
-    
-    if constexpr( isTriviallyCopyable_ )
-    {
-        copy(deviceView(), sortedView);
-    }
-    else if constexpr( isHostAccessible_)
-    {
-        for(auto i=0u; i<newSize; i++)
-        {
-            view_[i] = sortedView[i];
-        }   
-    }
-    else
-    {
-        static_assert("Not a valid operation for this data type on memory device");
-    }
+    copy(deviceView(), sortedView);
+
     
 	return true;
 }
\ No newline at end of file
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
index d30d1589..f64f78b1 100644
--- a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
+++ b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp
@@ -39,9 +39,6 @@ Licence:
 namespace pFlow
 {
 
-//- Forward 
-template<typename T, typename MemorySpace>
-class VectorSingle;
 
 template<typename T, typename MemorySpace=void>
 class VectorSingle
@@ -101,6 +98,8 @@ private:
 		static constexpr
 		bool isTriviallyCopyable_ = std::is_trivially_copyable_v<T>;
 
+        static_assert(isTriviallyCopyable_, "This type is not trivially copyable");
+
 	  	/// Evaluate capacity based on the input size 
 		static INLINE_FUNCTION_H uint32 evalCapacity(uint32 n)
 		{
@@ -290,25 +289,7 @@ public:
 
 		template<typename MSpace>
 		INLINE_FUNCTION_H
-		void assignFromDevice(const VectorSingle<T, MSpace>& src, bool srcCapacity = true)
-		{
-			uint32 srcSize = src.size();
-			uint32 srcCap = src.capacity();
-
-			if(srcCapacity && srcCap != capacity()){
-				reallocateCapacitySize(srcCap, srcSize);
-			}
-			else {
-				changeSize(srcSize);
-			}
-
-			if constexpr(isTriviallyCopyable_){
-				copy(deviceView(), src.deviceView());
-			}
-			else{
-				static_assert("Not a valid operation for this data type ");
-			}
-		}
+		void assignFromDevice(const VectorSingle<T, MSpace>& src, bool srcCapacity = true);
 
 		INLINE_FUNCTION_H
 		void append(const ViewType1D<T,MemorySpace>& appVec);
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingles.hpp b/src/phasicFlow/containers/VectorHD/VectorSingles.hpp
index d0c3acdf..2479715a 100644
--- a/src/phasicFlow/containers/VectorHD/VectorSingles.hpp
+++ b/src/phasicFlow/containers/VectorHD/VectorSingles.hpp
@@ -25,6 +25,7 @@ Licence:
 
 #include "types.hpp"
 #include "VectorSingle.hpp"
+#include "wordVectorHost.hpp"
 
 namespace pFlow
 {
@@ -77,6 +78,8 @@ typedef VectorSingle<realx3x3> 				realx3x3Vector_D;
 
 typedef VectorSingle<realx3x3, HostSpace> 	realx3x3Vector_H;
 
+typedef VectorSingle<word, HostSpace>        wordVector_H;
+
 }
 
 #endif
\ No newline at end of file
diff --git a/src/phasicFlow/containers/VectorHD/wordVectorHost.cpp b/src/phasicFlow/containers/VectorHD/wordVectorHost.cpp
new file mode 100644
index 00000000..1941a239
--- /dev/null
+++ b/src/phasicFlow/containers/VectorHD/wordVectorHost.cpp
@@ -0,0 +1,23 @@
+/*------------------------------- 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 "wordVectorHost.hpp"
+
diff --git a/src/phasicFlow/containers/VectorHD/wordVectorHost.hpp b/src/phasicFlow/containers/VectorHD/wordVectorHost.hpp
new file mode 100644
index 00000000..e638bd27
--- /dev/null
+++ b/src/phasicFlow/containers/VectorHD/wordVectorHost.hpp
@@ -0,0 +1,558 @@
+/*------------------------------- 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 __wordVectorHost_hpp__
+#define __wordVectorHost_hpp__
+
+#include "VectorSingle.hpp"
+#include "Vector.hpp"
+
+namespace pFlow
+{
+
+
+template<>
+class VectorSingle<word, HostSpace>
+{
+public:
+	
+	//- typedefs for accessing data
+
+		using VectorType		= VectorSingle<word, HostSpace>;
+
+		using VectorTypeHost 	= VectorSingle<word, HostSpace>; 		
+
+		using iterator        = word*;
+
+		using const_iterator   = const word*;
+	  	
+		using reference       = word&;
+	  	
+	 	using const_reference  = const word&;
+
+		using value_type       = word;
+	  	
+		using pointer         = word*;
+	  	
+		using const_pointer    = const word*;
+  	
+	//- typedefs related to memory management 
+
+		using viewType 			= ViewType1D<word, HostSpace>;
+
+		using device_type 		= typename viewType::device_type;
+
+		using memory_space 		= typename viewType::memory_space;
+
+		using execution_space 	= typename viewType::execution_space;
+
+private:
+	
+	// - Data members
+
+		Vector<word>              container_;
+
+        mutable viewType          unManagedView_;
+
+	// - protected members and methods 		
+
+		/// Is the memory of this vector accessible from Host
+		static constexpr 
+		bool  isHostAccessible_ = true; 
+
+		/// Is the memory of this vector accessiple from Divce
+		static constexpr
+		bool isDeviceAccessible_ = false; 
+
+		static constexpr
+		bool isTriviallyCopyable_ = false; 
+
+	
+public:
+
+	/// Type info
+	TypeInfoTemplateNV111("VectorSingle", word , memoerySpaceName());
+
+	//// - Constructors
+
+		/// Empty vector
+		VectorSingle() = default;
+		
+		/// Empty vector with a name (capacity = 2) 
+		explicit VectorSingle(const word& name)
+        :
+            container_(name)
+        {}
+		
+		/// Vector with name and size n
+		VectorSingle(const word& name, uint32 n)
+        :
+            container_(name, n)
+        {}
+		
+		/// Vector with name, size and value
+		VectorSingle(const word& name, uint32 n, const word& val)
+        :
+            container_(name, n, val)
+        {}
+
+		/// Vector with name, size (n) and reserved capacity 
+		VectorSingle(const word& name, uint32 cap, uint32 n, const RESERVE& r )
+        :
+            container_(name, cap, n, r)
+        {}
+		
+		/// Construct with a name and form std::vector (host memory)
+		VectorSingle(const word& name, const std::vector<word> & src)
+        :
+            container_(name, src)
+        {} 
+
+		/// Construct with a name and form std::vector (host memory) and with a desired capacity. 
+		VectorSingle(const word& name, const std::vector<word> & src, uint32 cap)
+        :
+            container_(name, src, cap)
+        {}
+		
+		/// Copy construct (performs deep copy)
+		VectorSingle(const VectorSingle& src)=default;
+		
+    	/// Copy construct with a new name (perform deep copy)
+		VectorSingle(const word& name, const VectorSingle& src)
+        :
+            container_(name, src.container_)
+        {}
+
+				
+		/// Copy assignment (perform deep copy from rhs to *this)
+		VectorSingle& operator = (const VectorSingle& rhs)=default;
+		
+		/// Move construct
+		VectorSingle(VectorSingle&&) = default;   
+
+		/// Move assignment
+		VectorSingle& operator= (VectorSingle&&) = default; 
+
+		/// @brief Descructor
+        /// This may not destroy the underlying memory, sice view is 
+		/// shared_ptr and maybe referenced by another object too 
+		~VectorSingle() = default;
+
+		/// Clone as a uniquePtr (perform deep copy)
+		INLINE_FUNCTION_H 
+		uniquePtr<VectorSingle> clone() const
+        {
+            return makeUnique<VectorSingle>(*this);
+        }
+
+			
+	//// - Methods
+
+		/// Return *this
+		INLINE_FUNCTION_H
+		VectorType& VectorField()
+        {
+            return *this;
+        }
+
+		/// Return *this
+		INLINE_FUNCTION_H
+		const VectorType& VectorField()const
+        {
+            return *this;
+        }
+
+		/// Device view range [0,capcity)
+		INLINE_FUNCTION_H 
+		auto& deviceViewAll()
+        {
+            // return un-managed view 
+            unManagedView_ = viewType(container_.data(), container_.capacity());
+            return unManagedView_;
+        }
+
+		/// Device view range [0,capcity)
+		INLINE_FUNCTION_H 
+		const auto& deviceViewAll() const
+        {
+            // return un-managed view 
+            unManagedView_ = viewType(const_cast<word*>(container_.data()), container_.capacity());
+            return unManagedView_;
+        }
+
+		///  Device view range [0, size)
+		INLINE_FUNCTION_H 
+		auto deviceView()const
+        {
+            return viewType(const_cast<word*>(container_.data()), container_.size());
+        }
+
+		/// Return a view accessible on Host in range [0,capacity)
+		INLINE_FUNCTION_H 
+		auto hostViewAll()const
+        {
+            return viewType(const_cast<word*>(container_.data()), container_.capacity());
+        }
+		
+
+		/// Return a view accessible on Host in range [0,size)
+		INLINE_FUNCTION_H 
+		auto hostView()const
+        {
+            return viewType(const_cast<word*>(container_.data()), container_.size());
+        }
+
+		/// Name of the vector 
+		INLINE_FUNCTION_H 
+		word name()const
+        {
+            return container_.name();
+        }
+		
+		
+		/// Size of the vector 
+		INLINE_FUNCTION_H 
+		uint32 size()const
+        {
+            return container_.size();
+        }
+		
+
+		// Capcity of the vector 
+		INLINE_FUNCTION_H 
+		uint32 capacity()const
+        {
+            return container_.capacity();
+        }
+		
+		/// If vector is empty 
+		INLINE_FUNCTION_H 
+		bool empty()const
+        {
+            return container_.size()==0uL;
+        }
+
+		
+		/// Reserve capacity for vector 
+		/// Preserve the content.
+		INLINE_FUNCTION_H 
+		void reserve(uint32 cap)
+        {
+            container_.reserve(cap);
+        } 
+
+		/// Reallocate memory to new cap and set size to 0. 
+		/*INLINE_FUNCTION_H 
+		void reallocate(uint32 cap); 
+		
+		/// Reallocate memory to new cap and set size to newSize. 
+		/// Do not preserve the content
+		INLINE_FUNCTION_H 
+		void reallocate(uint32 cap, uint32 newSize);*/
+		
+		/// Resize the vector and preserve the content
+		INLINE_FUNCTION_H   
+		void resize(uint32 n)
+        {
+            container_.resize(n);
+        }
+
+		/// Resize the vector and assign the value to it. 
+		INLINE_FUNCTION_H  
+		void resize(uint32 n, const word& val)
+        {
+            container_.resize(n, val);
+        }
+
+		/// Clear the vector, but keep the allocated memory unchanged 
+		INLINE_FUNCTION_H 
+		void clear()
+        {
+            container_.clear();
+        }
+		
+		/// Fill the range [0,size) with val
+		INLINE_FUNCTION_H 
+		void fill(const word& val)
+        {
+            container_.fill(val);
+        }
+
+		INLINE_FUNCTION_H
+		void fill(rangeU32 r,  const word& val)
+        {
+            container_.fill(r.start(), r.end(), val);
+        }
+		
+		/// Change size of the vector and assign val to vector and 
+		INLINE_FUNCTION_H
+		void assign(size_t n, const word& val)
+        {
+            container_.assign(n, val);
+        }
+		
+		/// Assign source vector with specified capacity.
+		/// The size of *this becomes the size of src. 
+		INLINE_FUNCTION_H 
+		void assign(const std::vector<word>& src, uint32 cap)
+        {
+            container_.reserve(cap);
+            this->assign(src);
+        }
+		
+
+		/// Assign source vector.
+		/// The size of *this becomes the size of src. 
+		/// The capacity of *this becomes the capacity of src.
+		INLINE_FUNCTION_H 
+		void assign(const std::vector<word>& src)
+        {
+            container_.assign(src.begin(), src.end());
+        }
+
+		/// Assign source vector from host side.
+		/// The size of *this becomes the size of src. 
+		/// The capacity of *this becomes the capacity of src.
+		INLINE_FUNCTION_H
+		void assignFromHost(const VectorTypeHost& src)
+        {
+            notImplementedFunction;
+            fatalExit;
+        }
+
+		INLINE_FUNCTION_H
+		void assign(const VectorType& src, bool srcCapacity = true)
+        {
+            notImplementedFunction;
+            fatalExit;
+        }
+
+		template<typename MSpace>
+		INLINE_FUNCTION_H
+		void assignFromDevice(const VectorSingle<word, MSpace>& src, bool srcCapacity = true)
+		{
+			notImplementedFunction;
+            fatalExit;
+		}
+
+		/*INLINE_FUNCTION_H
+		void append(const ViewType1D<T,MemorySpace>& appVec);
+
+		INLINE_FUNCTION_H
+		void append(const std::vector<T>& appVec);
+
+		INLINE_FUNCTION_H
+		void append(const VectorType& appVec);*/
+		
+		INLINE_FUNCTION_H
+		auto getSpan()
+        {
+            return span<word>(container_.data(), container_.size());
+        }
+				
+		INLINE_FUNCTION_H 
+		auto getSpan()const
+        {
+            return span<word>(const_cast<word*>(container_.data()), container_.size());
+        }
+		
+		INLINE_FUNCTION_H
+		bool insertSetElement(const uint32IndexContainer& indices, const word& val)
+        {
+            notImplementedFunction;
+            return false;
+        }
+		
+		INLINE_FUNCTION_H
+		bool insertSetElement(const uint32IndexContainer& indices, const std::vector<word>& vals)
+        {
+            notImplementedFunction;
+            return false;
+        }
+		
+		INLINE_FUNCTION_H
+		bool insertSetElement(
+			const uint32IndexContainer& indices, 
+			const ViewType1D<word, memory_space> vals)
+        {
+            notImplementedFunction;
+            return false;
+        }
+
+		INLINE_FUNCTION_H
+		bool reorderItems(const uint32IndexContainer& indices)
+        {
+            notImplementedFunction;
+            return false;
+        }
+
+		/// @brief push a new element at the end (host call only)
+		///  resize if necessary and works on host accessible vector.
+		
+		void push_back(const word& val)
+		{
+			container_.push_back(val);
+		}
+
+		INLINE_FUNCTION_H 
+		pointer data(){
+			return container_.data();
+		}
+
+		INLINE_FUNCTION_H 
+		const_pointer data()const{
+			return container_.data();
+		}
+
+		/// Return begin iterator. It works when devices is host accessible.
+		auto
+		begin(){
+			return container_.begin();
+		}
+
+		/// Return begin iterator. it works when host is accessible.
+		const auto
+		begin()const {
+			return container_.begin();
+		}
+		
+		auto end(){
+			return container_.end();
+		}
+
+		/// Return end iterator. it works when host is accessible.
+		const auto end()const{
+			return container_.end();
+		}
+				
+		/// Return reference to element i. it works when host is accessible.		
+		word& operator[](size_t i){
+			return container_[i];
+		}
+
+		const word& operator[](size_t i)const{
+			return container_[i];
+		}
+
+	//// - IO operations
+
+		/// Read vector from stream
+		FUNCTION_H
+		bool read(iIstream& is)
+		{
+			return container_.read(is);
+		}
+
+		/// Read vector from stream
+		FUNCTION_H
+		bool read(iIstream& is, const IOPattern& iop)
+		{
+			return container_.read(is, iop);
+		}
+
+		/// Write the vector to os
+		FUNCTION_H
+		bool write(iOstream& os, const IOPattern& iop)const
+		{
+			return container_.write(os, iop);
+		}
+
+		FUNCTION_H
+		bool write(iOstream& os)const
+		{
+			return container_.write(os);	
+		}
+
+		template<typename HostMask>
+		FUNCTION_H
+		bool write(iOstream& os, const IOPattern& iop, const HostMask& mask)const
+		{
+			notImplementedFunction;
+			return false;
+		}
+
+
+		/// Name of the memory space
+		static  
+		constexpr const char* memoerySpaceName()
+		{
+			return memory_space::name();
+		}
+
+}; // class wordVectorHost
+
+inline iIstream& operator >> (iIstream & is, VectorSingle<word, HostSpace> & ivec )
+{
+	if( !ivec.read(is) )
+	{
+		ioErrorInFile (is.name(), is.lineNumber());
+		fatalExit;
+	}
+	return is;
+}
+
+inline iOstream& operator << (iOstream& os, const VectorSingle<word, HostSpace>& ovec )
+{
+	
+	if( !ovec.write(os) )
+	{
+		ioErrorInFile(os.name(), os.lineNumber());
+		fatalExit;
+	}
+
+	return os; 
+}
+
+
+} // - pFlow
+
+
+
+#endif
+
+
+/*
+INLINE_FUNCTION_H
+		bool append(const deviceViewType1D<T>& dVec, size_t numElems)
+		{
+
+			if(numElems == 0 )return true;
+			auto oldSize = size_;
+			auto newSize = size_ + numElems;
+
+			if(this->empty() || newSize > capacity() )
+			{
+				resize(newSize);
+			}
+			else
+			{
+				size_ = size_+numElems;
+			}
+			
+			auto dSubView = Kokkos::subview(view_, Kokkos::make_pair(oldSize, newSize)); 
+			copy(dSubView, dVec);
+
+			return true;
+		}
+
+		INLINE_FUNCTION_H
+		bool append(const VectorSingle& Vec)
+		{
+			return append(Vec.deviceView(), Vec.size());
+		}*/
\ No newline at end of file
diff --git a/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp b/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp
index 9570f86b..20267f33 100644
--- a/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp
+++ b/src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp
@@ -135,18 +135,31 @@ public:
 
 	FieldAccessType thisField()const
 	{
-		return FieldAccessType(
-			this->size(),
-			this->indexList().deviceViewAll(),
-			internal_.deviceViewAll());
+        if constexpr(isDeviceAccessible<execution_space>())
+            return FieldAccessType(
+                this->size(),
+                this->indexList().deviceViewAll(),
+                internal_.deviceViewAll());
+        else
+            return FieldAccessType(
+                this->size(),
+                this->boundary().indexListHost().deviceViewAll(),
+                internal_.deviceViewAll());
 	}
 
 	FieldAccessType mirrorField()const
 	{
-		return FieldAccessType(
-			this->mirrorBoundary().size(),
-			this->mirrorBoundary().indexList().deviceViewAll(),
-			internal_.deviceViewAll());
+        if constexpr(isDeviceAccessible<execution_space>())
+            return FieldAccessType(
+                this->mirrorBoundary().size(),
+                this->mirrorBoundary().indexList().deviceViewAll(),
+                internal_.deviceViewAll());
+        else
+            return FieldAccessType(
+                this->mirrorBoundary().size(),
+                this->mirrorBoundary().indexListHost().deviceViewAll(),
+                internal_.deviceViewAll());
+
 	}
 	
 	virtual 
diff --git a/src/phasicFlow/containers/pointField/boundaryField/boundaryFieldList.hpp b/src/phasicFlow/containers/pointField/boundaryField/boundaryFieldList.hpp
index 31a13661..9ef6fcb1 100644
--- a/src/phasicFlow/containers/pointField/boundaryField/boundaryFieldList.hpp
+++ b/src/phasicFlow/containers/pointField/boundaryField/boundaryFieldList.hpp
@@ -43,7 +43,7 @@ protected:
 
     const boundaryList&         boundaries_;
 
-    uint32                      slaveToMasterUpdateIter_ = -1;
+    uint32                      slaveToMasterUpdateIter_ = static_cast<uint32>(-1);
 
 public:
 
@@ -95,7 +95,7 @@ public:
 
     bool slaveToMasterUpdateRequested()const
     {
-        return slaveToMasterUpdateIter_ != -1;
+        return slaveToMasterUpdateIter_ != static_cast<uint32>(-1);
     }
 
     
diff --git a/src/phasicFlow/dictionary/dictionary.hpp b/src/phasicFlow/dictionary/dictionary.hpp
index 9baba109..1684987c 100644
--- a/src/phasicFlow/dictionary/dictionary.hpp
+++ b/src/phasicFlow/dictionary/dictionary.hpp
@@ -302,10 +302,10 @@ public:
 	//// - IO operations 
 
 		/// read from stream
-		virtual bool read(iIstream& is);
+		bool read(iIstream& is) override;
 
 		/// write to stream
-		virtual bool write(iOstream& os) const;
+		bool write(iOstream& os) const override;
 
 };
 
diff --git a/src/phasicFlow/dictionary/fileDictionary.hpp b/src/phasicFlow/dictionary/fileDictionary.hpp
index b8685dfd..a74c989e 100644
--- a/src/phasicFlow/dictionary/fileDictionary.hpp
+++ b/src/phasicFlow/dictionary/fileDictionary.hpp
@@ -46,6 +46,9 @@ public:
     const dictionary& dict, 
     repository* owner=nullptr);
 
+    using dictionary::read;
+    using dictionary::write;
+
     /// read from stream 
     bool read(iIstream& is, const IOPattern& iop) override;
 
diff --git a/src/phasicFlow/globals/pFlowMacros.hpp b/src/phasicFlow/globals/pFlowMacros.hpp
index 544fc150..8973b6ba 100755
--- a/src/phasicFlow/globals/pFlowMacros.hpp
+++ b/src/phasicFlow/globals/pFlowMacros.hpp
@@ -37,6 +37,10 @@ Licence:
 
 #define CONSUME_PARAM(x) (void)(x);
 
+#if defined(pFlow_Build_Cuda) && !defined(__CUDACC__)
+#define __CUDACC__
+#endif
+
 #ifdef __CUDACC__
     #define INLINE_FUNCTION_HD inline __host__ __device__
     #define INLINE_FUNCTION_D  inline __device__
diff --git a/src/phasicFlow/processors/localProcessors.hpp b/src/phasicFlow/processors/localProcessors.hpp
index eab18c08..1e1d6d3e 100644
--- a/src/phasicFlow/processors/localProcessors.hpp
+++ b/src/phasicFlow/processors/localProcessors.hpp
@@ -145,6 +145,15 @@ public:
         }
     #endif
 
+    static
+    bool builtForMPI()
+    {
+        #ifdef pFlow_Build_MPI
+            return true;
+        #else
+            return false;
+        #endif
+    }
 };
 
 }
diff --git a/src/phasicFlow/processors/processors.cpp b/src/phasicFlow/processors/processors.cpp
index f5eb066a..078c9a5d 100644
--- a/src/phasicFlow/processors/processors.cpp
+++ b/src/phasicFlow/processors/processors.cpp
@@ -29,7 +29,6 @@ Licence:
 #include "processors.hpp"
 #include "streams.hpp"
 
-static int numVarsInitialized__ = 0;
 
 #ifdef pFlow_Build_MPI
 
diff --git a/src/phasicFlow/repository/IOobject/IOfileHeader.cpp b/src/phasicFlow/repository/IOobject/IOfileHeader.cpp
index 96a49f67..f67eec28 100644
--- a/src/phasicFlow/repository/IOobject/IOfileHeader.cpp
+++ b/src/phasicFlow/repository/IOobject/IOfileHeader.cpp
@@ -42,12 +42,22 @@ pFlow::uniquePtr<pFlow::oFstream> pFlow::IOfileHeader::outStream()const
 	return osPtr;
 }
 
-pFlow::IOfileHeader::IOfileHeader
-(
-	const objectFile& objf
-)
-:
-	objectFile(objf)
+pFlow::uniquePtr<pFlow::oFstream> pFlow::IOfileHeader::dummyOutStream() const
+{
+    auto osPtr = makeUnique<oFstream>( CWD()+word("dummyFile") , outFileBinary());
+    
+    if(osPtr && owner())
+	{
+	 	auto outPrecision = owner()->outFilePrecision();
+	 	osPtr->precision(static_cast<int>(outPrecision));
+	}
+
+	return osPtr;
+}
+
+pFlow::IOfileHeader::IOfileHeader(
+    const objectFile &objf)
+    : objectFile(objf)
 {}
 
 pFlow::fileSystem pFlow::IOfileHeader::path() const
diff --git a/src/phasicFlow/repository/IOobject/IOfileHeader.hpp b/src/phasicFlow/repository/IOobject/IOfileHeader.hpp
index 1d2a1e4c..d4549c01 100644
--- a/src/phasicFlow/repository/IOobject/IOfileHeader.hpp
+++ b/src/phasicFlow/repository/IOobject/IOfileHeader.hpp
@@ -58,6 +58,8 @@ protected:
 		// - ouput file stream
 		uniquePtr<oFstream> outStream()const;
 
+        uniquePtr<oFstream> dummyOutStream()const;
+
 public:
 
 	// with owner 
diff --git a/src/phasicFlow/repository/IOobject/IOobject.cpp b/src/phasicFlow/repository/IOobject/IOobject.cpp
index b5857ea1..6fae2f0a 100644
--- a/src/phasicFlow/repository/IOobject/IOobject.cpp
+++ b/src/phasicFlow/repository/IOobject/IOobject.cpp
@@ -122,18 +122,35 @@ bool pFlow::IOobject::writeObject() const
 
 	if(ioPattern().thisCallWrite())
 	{
-		
-        if(auto ptrOS = outStream(); ptrOS )
+		if( ioPattern().thisProcWriteData())
         {
-            return writeObject(ptrOS());
+            if(auto ptrOS = outStream(); ptrOS )
+            {
+                return writeObject(ptrOS());
+            }
+            else
+            {
+                warningInFunction<< 
+                "error in opening file "<< path() <<endl;
+                return false;
+            }
         }
         else
         {
-            warningInFunction<< 
-            "error in opening file "<< path() <<endl;
-            return false;
+            
+            if(auto ptrOS = dummyOutStream(); ptrOS )
+            {
+                return writeObject(ptrOS());
+            }
+            else
+            {
+                warningInFunction<< 
+                "error in opening file "<< path() <<endl;
+                return false;
+            }
         }
         
+        
 	}
 
 	return true;
diff --git a/src/phasicFlow/repository/Time/baseTimeControl.cpp b/src/phasicFlow/repository/Time/baseTimeControl.cpp
index 78d04225..d44b9aac 100644
--- a/src/phasicFlow/repository/Time/baseTimeControl.cpp
+++ b/src/phasicFlow/repository/Time/baseTimeControl.cpp
@@ -68,7 +68,7 @@ pFlow::baseTimeControl::baseTimeControl
 pFlow::baseTimeControl::baseTimeControl(int32 start, int32 end, int32 stride, const word &intervalPrefix)
 :
 	isTimeStep_(true),
-	iRange_(start, end, max(stride,1)),
+	iRange_(start, end, std::max(stride,1)),
 	intervalPrefix_(
 		intervalPrefix.size()==0uL? word("interval"): intervalPrefix+"Interval"
 	)
diff --git a/src/phasicFlow/repository/Time/timeControl.cpp b/src/phasicFlow/repository/Time/timeControl.cpp
index 206da94e..b6305545 100644
--- a/src/phasicFlow/repository/Time/timeControl.cpp
+++ b/src/phasicFlow/repository/Time/timeControl.cpp
@@ -18,7 +18,7 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-
+#include "math.hpp"
 #include "timeControl.hpp"
 #include "dictionary.hpp"
 
@@ -137,14 +137,14 @@ pFlow::word pFlow::timeControl::timeName()const
 bool pFlow::timeControl::finalTime()const
 {
 	if( currentTime_ >= endTime_ ) return true;
-	if( abs(currentTime_-endTime_) < 0.5*dt_ )return true;
+	if( std::abs(currentTime_-endTime_) < 0.5*dt_ )return true;
 	return false;	
 }
 
 bool pFlow::timeControl::reachedStopAt()const
 {
 	if( currentTime_ >= stopAt_ ) return true;
-	if( abs(currentTime_-stopAt_) < 0.5*dt_ )return true;
+	if( std::abs(currentTime_-stopAt_) < 0.5*dt_ )return true;
 	return false;
 }
 
@@ -154,7 +154,7 @@ void pFlow::timeControl::checkForOutputToFile()
 	bool save = false;
 	if(managedExternaly_)
 	{
-		if( abs(currentTime_-writeTime_) < 0.5*dt_)
+		if( std::abs(currentTime_-writeTime_) < 0.5*dt_)
 		{
 			save = true;
 			lastSaved_ = currentTime_;
@@ -162,12 +162,12 @@ void pFlow::timeControl::checkForOutputToFile()
 	}
 	else
 	{
-		if ( abs(currentTime_ - lastSaved_ - saveInterval_) < 0.5 * dt_ )
+		if ( std::abs(currentTime_ - lastSaved_ - saveInterval_) < 0.5 * dt_ )
 		{	
 			lastSaved_ = currentTime_;
 			save = true;
 		}
-		else if( abs(currentTime_ - lastSaved_) < min( pow(10.0,-1.0*timePrecision_), 0.5 *dt_) )
+		else if( std::abs(currentTime_ - lastSaved_) < std::min( pow(10.0,-1.0*timePrecision_), 0.5 *dt_) )
 		{
 			lastSaved_ = currentTime_;
 			save = true;
diff --git a/src/phasicFlow/smartPointers/uniquePtr.hpp b/src/phasicFlow/smartPointers/uniquePtr.hpp
index c1363cde..7bb6a6b9 100644
--- a/src/phasicFlow/smartPointers/uniquePtr.hpp
+++ b/src/phasicFlow/smartPointers/uniquePtr.hpp
@@ -37,16 +37,15 @@ namespace pFlow
 
 
 template<
-	typename T, 
-	typename Deleter = std::default_delete<T>
+	typename T
 >
 class uniquePtr
 :
-	public std::unique_ptr<T, Deleter>
+	public std::unique_ptr<T>
 {
 public:
 
-	using uniquePtrType = std::unique_ptr<T, Deleter>;
+	using uniquePtrType = std::unique_ptr<T>;
 
 	// using base constructors 
 	using uniquePtrType::unique_ptr;
diff --git a/src/phasicFlow/streams/TStream/iTstream.cpp b/src/phasicFlow/streams/TStream/iTstream.cpp
index 49b80cee..ac8e503e 100755
--- a/src/phasicFlow/streams/TStream/iTstream.cpp
+++ b/src/phasicFlow/streams/TStream/iTstream.cpp
@@ -325,7 +325,7 @@ void pFlow::iTstream::reset()
 size_t pFlow::iTstream::tell() 
 {
     notImplementedFunction;
-    return -1;
+    return static_cast<size_t>(-1);
 }
 
 const pFlow::tokenList& pFlow::iTstream::tokens()const
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp
index e051c803..c1967bc8 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/scatteredFieldAccess.hpp
@@ -124,6 +124,7 @@ public:
 		return indices_;
 	}
 
+    INLINE_FUNCTION_HD
 	uint32 index(uint32 i)const
 	{
 		return indices_[i];
diff --git a/src/phasicFlow/structuredData/boundaries/boundaryReflective/boundaryReflective.cpp b/src/phasicFlow/structuredData/boundaries/boundaryReflective/boundaryReflective.cpp
index 6d818a23..74ee4d75 100644
--- a/src/phasicFlow/structuredData/boundaries/boundaryReflective/boundaryReflective.cpp
+++ b/src/phasicFlow/structuredData/boundaries/boundaryReflective/boundaryReflective.cpp
@@ -125,6 +125,8 @@ bool pFlow::boundaryReflective::afterIteration
 
     const auto& velocity = time().lookupObject<realx3PointField_D>(velocityName_);
     const auto& velocityD = velocity.deviceViewAll();
+    
+    const auto restitution = restitution_;
 
     Kokkos::parallel_for(
         "pFlow::boundaryReflective::velocityChange",
@@ -137,7 +139,7 @@ bool pFlow::boundaryReflective::afterIteration
             if(vn < 0)
             {
                 realx3 vt = vel - vn*p.normal();
-                vel = restitution_*(vt - vn*p.normal());
+                vel = restitution*(vt - vn*p.normal());
             }
         }
     );
diff --git a/src/phasicFlow/structuredData/cylinder/cylinder.cpp b/src/phasicFlow/structuredData/cylinder/cylinder.cpp
index 8c706d3e..014462c2 100644
--- a/src/phasicFlow/structuredData/cylinder/cylinder.cpp
+++ b/src/phasicFlow/structuredData/cylinder/cylinder.cpp
@@ -18,7 +18,7 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-
+#include "math.hpp"
 #include "cylinder.hpp"
 #include "zAxis.hpp"
 #include "streams.hpp"
diff --git a/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.cpp b/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.cpp
index 4dc8ff0d..0ea527f6 100644
--- a/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.cpp
+++ b/src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.cpp
@@ -317,6 +317,9 @@ pFlow::internalPoints::insertPoints(
 
 	auto aRange = pFlagsD_.activeRange();
 	uint32 emptySpots = pFlagsD_.capacity() - pFlagsD_.numActive();
+
+    if(emptySpots!= 0) emptySpots--;
+
 	message msg;
 
 	if( numNew > emptySpots )
@@ -348,8 +351,20 @@ pFlow::internalPoints::insertPoints(
 	// we should fill the scattered empty spots
 	else
 	{
-		notImplementedFunction;
-		return false;
+		auto newIndices = pFlagsD_.getEmptyPoints(numNew);
+		if(numNew != newIndices.size())
+		{
+			fatalErrorInFunction<<"not enough empty points in pointFlag"<<
+			numNew<< " "<<newIndices.size() <<endl;
+			pOutput<< pFlagsD_.capacity()<<endl;
+			pOutput<< pFlagsD_.numActive()<<endl;
+			return false;
+		}
+		
+		varList.emplaceBack<uint32IndexContainer>(
+			msg.addAndName(message::ITEM_INSERT),
+			newIndices
+		);
 	}
 
 	const auto& indices = varList.getObject<uint32IndexContainer>(
@@ -450,7 +465,7 @@ bool pFlow::internalPoints::insertPointsOnly(
 	}// we should fill the scattered empty spots
 	else
 	{
-		pOutput<<"numNew to be inserted "<< numNew <<endl;
+		
 		auto newIndices = pFlagsD_.getEmptyPoints(numNew);
 		if(numNew != newIndices.size())
 		{
@@ -460,7 +475,7 @@ bool pFlow::internalPoints::insertPointsOnly(
 			pOutput<< pFlagsD_.numActive()<<endl;
 			return false;
 		}
-		pOutput<<newIndices<<endl;
+		
 		varList.emplaceBack<uint32IndexContainer>(
 			msg.addAndName(message::ITEM_INSERT),
 			newIndices
diff --git a/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp b/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp
index de8940ba..e76533b9 100644
--- a/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp
+++ b/src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp
@@ -63,12 +63,13 @@ pFlow::ViewType1D<pFlow::uint32, typename pFlow::pointFlag<ExecutionSpace>::memo
 
 	ViewType1D<uint32,memory_space> emptyPoints("emptyPoints", numToGet);
 
-	
+	auto flags = flags_;
+
 	Kokkos::parallel_for(
 		"getEmptyPoints",
 		rPolicy(0, cap),
 		LAMBDA_HD(uint32 i){
-			indices(i) = flags_[i] == DELETED;
+			indices(i) = flags[i] == DELETED;
 	});
 	Kokkos::fence();
 
@@ -151,7 +152,7 @@ pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markOutOfBoxDelete
 	{
 		minRange = 0;
 		maxRange = 0;
-		numDeleted == numActive_;
+		numDeleted = numActive_;
 	}
 	else
 	{
@@ -180,7 +181,8 @@ pFlow::pointFlag<ExecutionSpace>::addInternalPoints
 	uint32 maxRange;
 	uint32 numAdded = 0;
 	
-	
+	const auto& flag = flags_;
+
 	Kokkos::parallel_reduce(
 		"pointFlagKernels::addInternalPoints",
 		deviceRPolicyStatic(0, points.extent(0)),
@@ -191,10 +193,10 @@ pFlow::pointFlag<ExecutionSpace>::addInternalPoints
 			uint32& addToUpdate)
 		{
 			uint32 idx = points(i);
-			if( flags_[idx] <= DELETED) addToUpdate ++;
+			if( flag[idx] <= DELETED) addToUpdate ++;
 			minUpdate = min(minUpdate,idx);
 			maxUpdate = max(maxUpdate,idx);
-			flags_[idx] = INTERNAL;	
+			flag[idx] = INTERNAL;	
 		},
 		Kokkos::Min<uint32>(minRange), 
 		Kokkos::Max<uint32>(maxRange),
@@ -247,7 +249,7 @@ bool pFlow::pointFlag<ExecutionSpace>::deletePoints
 	if(numDeleted >= numActive_)
 	{
 		activeRange_ = {0, 0};
-		numDeleted == numActive_;
+		numDeleted = numActive_;
 	}
 		
 	numActive_ 	= numActive_ - numDeleted;
@@ -290,7 +292,7 @@ bool pFlow::pointFlag<ExecutionSpace>::deletePoints
 	if(numDeleted >= numActive_)
 	{
 		activeRange_ = {0, 0};
-		numDeleted == numActive_;
+		numDeleted = numActive_;
 	}
 		
 	numActive_ 	= numActive_ - numDeleted;
@@ -309,13 +311,14 @@ bool pFlow::pointFlag<ExecutionSpace>::changeFlags
 )
 {
 	auto flg = getBoundaryFlag(boundaryIndex);
+    const auto& flags = flags_;
 	Kokkos::parallel_for
 	(
 		"pointFlag::changeFlags",
 		rPolicy(0, changePoints.size()),
 		LAMBDA_HD(uint32 i)
 		{
-			flags_[changePoints(i)] = flg;
+			flags[changePoints(i)] = flg;
 		}
 	);
 	Kokkos::fence();
diff --git a/src/phasicFlow/structuredData/zAxis/array2D.hpp b/src/phasicFlow/structuredData/zAxis/array2D.hpp
index f76e6080..3a1bb059 100644
--- a/src/phasicFlow/structuredData/zAxis/array2D.hpp
+++ b/src/phasicFlow/structuredData/zAxis/array2D.hpp
@@ -33,12 +33,12 @@ struct array2D
 
     constexpr size_t nCols()const noexcept
     {
-        return nCols;
+        return nCol;
     }
 
     constexpr size_t nRows()const noexcept
     {
-        return nRows;
+        return nRow;
     }
     
 };
diff --git a/src/phasicFlow/triSurface/triangleFunctions.hpp b/src/phasicFlow/triSurface/triangleFunctions.hpp
index 7345c07f..c6812d2f 100644
--- a/src/phasicFlow/triSurface/triangleFunctions.hpp
+++ b/src/phasicFlow/triSurface/triangleFunctions.hpp
@@ -39,7 +39,7 @@ realx3 normal(const realx3& p1, const realx3& p2, const realx3& p3)
 {
 	auto n = cross(p2-p1, p3-p1);
 	if( equal(n.length(), 0.0) )
-    	return zero3;
+    	return realx3(0);
 	else
 		return normalize(n);
 }
diff --git a/src/phasicFlow/types/basicTypes/math.hpp b/src/phasicFlow/types/basicTypes/math.hpp
index ee140058..fbe3a6f1 100755
--- a/src/phasicFlow/types/basicTypes/math.hpp
+++ b/src/phasicFlow/types/basicTypes/math.hpp
@@ -21,15 +21,16 @@ Licence:
 #ifndef __math_hpp__
 #define __math_hpp__
 
+
+#include "builtinTypes.hpp"
+#include "pFlowMacros.hpp"
+
 #ifdef __CUDACC__
 #include "math.h"
 #else
 #include <cmath>
 #endif
 
-#include "builtinTypes.hpp"
-#include "pFlowMacros.hpp"
-
 //* * * * * * * * * * *  List of functinos * * * * * * * * //
 // abs, mod, exp, log, log10, pow, sqrt, cbrt
 // sin, cos, tan, asin, acos, atan, atan2
@@ -51,22 +52,26 @@ abs(real x)
 #endif
 }
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 int64
 abs(int64 x)
 {
+#ifdef __CUDACC__
+    return ::abs(x);
+#else
 	return std::abs(x);
-}
 #endif
+}
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD int32
 abs(int32 x)
 {
+#ifdef __CUDACC__
+    return ::abs(x);
+#else
 	return std::abs(x);
-}
 #endif
+}
 
 INLINE_FUNCTION_HD real
 mod(real x, real y)
@@ -78,6 +83,7 @@ mod(real x, real y)
 #endif
 }
 
+
 INLINE_FUNCTION_HD int64
 mod(int64 x, int64 y)
 {
@@ -102,147 +108,188 @@ mod(uint32 x, uint32 y)
 	return x % y;
 }
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD real
 remainder(real x, real y)
 {
+#ifdef __CUDACC__
+    return ::remainder(x,y);
+#else
 	return std::remainder(x, y);
-}
 #endif
+}
 
-#ifndef __CUDACC__
-INLINE_FUNCTION_H
+
+INLINE_FUNCTION_HD
 real
 exp(real x)
 {
+#ifdef __CUDACC__
+    return ::exp(x);
+#else
 	return std::exp(x);
-}
 #endif
+}
+
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 log(real x)
 {
+#ifdef __CUDACC__
+    return ::log(x);
+#else
 	return std::log(x);
-}
 #endif
+}
+
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 log10(real x)
 {
+#ifdef __CUDACC__
+    return ::log10(x);
+#else
 	return std::log10(x);
-}
 #endif
+}
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 pow(real x, real y)
 {
+#ifdef __CUDACC__
+	return ::pow(x,y);
+#else
 	return std::pow(x, y);
-}
 #endif
+}
+
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 sqrt(real x)
 {
+#ifdef __CUDACC__
+    return ::sqrt(x);
+#else
 	return std::sqrt(x);
-}
 #endif
+}
+
+
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 cbrt(real x)
 {
+#ifdef __CUDACC__
+    return ::cbrt(x);
+#else
 	return std::cbrt(x);
+#endif 
 }
-#endif
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 sin(real x)
 {
+#ifdef __CUDACC__
+    return ::sin(x);
+#else
 	return std::sin(x);
-}
 #endif
+}
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 cos(real x)
 {
+#ifdef __CUDACC__
+    return ::cos(x);
+#else
 	return std::cos(x);
-}
 #endif
+}
+
+
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 tan(real x)
 {
+#ifdef __CUDACC__
+    return ::tan(x);
+#else
 	return std::tan(x);
-}
 #endif
+}
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 asin(real x)
 {
+#ifdef __CUDACC__
+    return ::asin(x);
+#else
 	return std::asin(x);
-}
 #endif
+}
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 acos(real x)
 {
+#ifdef __CUDACC__
+    return ::acos(x);
+#else
 	return std::acos(x);
-}
 #endif
+}
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 atan(real x)
 {
+#ifdef __CUDACC__
+    return ::atan(x);
+#else
 	return std::atan(x);
-}
 #endif
+}
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD real
 atan2(real y, real x)
 {
+#ifdef __CUDACC__
+    return ::atan2(y,x);
+#else
 	return std::atan2(y, x);
-}
 #endif
+}
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 sinh(real x)
 {
+#ifdef __CUDACC__
+    return ::sinh(x);
+#else
 	return std::sinh(x);
-}
 #endif
+}
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 cosh(real x)
 {
+#ifdef __CUDACC__
+    return ::cosh(x);
+#else
 	return std::cosh(x);
-}
 #endif
+}
+
 
 INLINE_FUNCTION_HD real
 tanh(real x)
@@ -274,14 +321,17 @@ acosh(real x)
 #endif
 }
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 real
 atanh(real x)
 {
+#ifdef __CUDACC__
+    return ::atanh(x);
+#else
 	return std::atanh(x);
-}
 #endif
+}
+
 
 INLINE_FUNCTION_HD real
 min(real x, real y)
@@ -293,40 +343,54 @@ min(real x, real y)
 #endif
 }
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 int64
 min(int32 x, int32 y)
 {
+#ifdef __CUDACC__
+    return ::min(x,y);
+#else
 	return std::min(x, y);
-}
 #endif
+}
+
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 int64
 min(int64 x, int64 y)
 {
+#ifdef __CUDACC__
+    return ::min(x,y);
+#else
 	return std::min(x, y);
-}
 #endif
+}
+
+
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD
 uint64
 min(uint64 x, uint64 y)
 {
+#ifdef __CUDACC__
+    return ::min(x,y);
+#else
 	return std::min(x, y);
-}
 #endif
+}
+
+
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD uint32
 min(uint32 x, uint32 y)
 {
+#ifdef __CUDACC__
+    return ::min(x,y);
+#else
 	return std::min(x, y);
-}
 #endif
+}
+
 
 INLINE_FUNCTION_HD real
 max(real x, real y)
@@ -338,37 +402,48 @@ max(real x, real y)
 #endif
 }
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD int64
 max(int64 x, int64 y)
 {
+#ifdef __CUDACC__
+	return ::max(x, y);
+#else
 	return std::max(x, y);
-}
 #endif
+}
+
+
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD int32
 max(int32 x, int32 y)
 {
+#ifdef __CUDACC__
+	return ::max(x, y);
+#else
 	return std::max(x, y);
-}
 #endif
+}
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD uint64
 max(uint64 x, uint64 y)
 {
+#ifdef __CUDACC__
+	return ::max(x, y);
+#else
 	return std::max(x, y);
-}
 #endif
+}
 
-#ifndef __CUDACC__
 INLINE_FUNCTION_HD uint32
 max(uint32 x, uint32 y)
 {
+#ifdef __CUDACC__
+	return ::max(x, y);
+#else
 	return std::max(x, y);
-}
 #endif
+}
+
 
 } // pFlow
 
diff --git a/src/setHelpers/setProperty.hpp b/src/setHelpers/setProperty.hpp
index 43ec5ef2..978662d4 100644
--- a/src/setHelpers/setProperty.hpp
+++ b/src/setHelpers/setProperty.hpp
@@ -23,6 +23,6 @@ Licence:
 
 REPORT(0)<<"\nReading proprties . . . "<<END_REPORT;
 
-auto proprties = property(propertyFile__, Control.caseSetup().path());
+auto proprties = pFlow::property(pFlow::propertyFile__, Control.caseSetup().path());
 
 #endif // __setProperty_hpp__
diff --git a/utilities/CMakeLists.txt b/utilities/CMakeLists.txt
index 37a70ad1..18f870cf 100644
--- a/utilities/CMakeLists.txt
+++ b/utilities/CMakeLists.txt
@@ -1,5 +1,5 @@
 
-#add_subdirectory(checkPhasicFlow)
+add_subdirectory(checkPhasicFlow)
 
 add_subdirectory(particlesPhasicFlow)
 
diff --git a/utilities/Utilities/CMakeLists.txt b/utilities/Utilities/CMakeLists.txt
index 729958ed..598a542a 100644
--- a/utilities/Utilities/CMakeLists.txt
+++ b/utilities/Utilities/CMakeLists.txt
@@ -5,7 +5,8 @@ readFromTimeFolder.cpp
 vtkFile/vtkFile.cpp
 geometryPhasicFlow/Wall/Wall.cpp 
 geometryPhasicFlow/planeWall/planeWall.cpp 
-#geometryPhasicFlow/stlWall/stlWall.cpp
+geometryPhasicFlow/stlWall/stlFile.cpp
+geometryPhasicFlow/stlWall/stlWall.cpp
 geometryPhasicFlow/cylinderWall/cylinderWall.cpp
 geometryPhasicFlow/cuboidWall/cuboidWall.cpp
 )
diff --git a/src/phasicFlow/triSurface/stlFile.cpp b/utilities/Utilities/geometryPhasicFlow/stlWall/stlFile.cpp
similarity index 97%
rename from src/phasicFlow/triSurface/stlFile.cpp
rename to utilities/Utilities/geometryPhasicFlow/stlWall/stlFile.cpp
index 0dc3cde6..c0366851 100644
--- a/src/phasicFlow/triSurface/stlFile.cpp
+++ b/utilities/Utilities/geometryPhasicFlow/stlWall/stlFile.cpp
@@ -272,7 +272,7 @@ void pFlow::stlFile::addSolid
 	const realx3x3Vector& vertecies
 )
 {
-	solids_.push_backSafe(vertecies);
+	solids_.push_back(makeUnique<realx3x3Vector>(vertecies));
 	solidNames_.push_back(name);
 }
 
@@ -283,7 +283,7 @@ void pFlow::stlFile::addSolid
 	realx3x3Vector&& vertecies
 )
 {
-	solids_.push_backSafe(std::move(vertecies));
+	solids_.push_back(makeUnique<realx3x3Vector>(vertecies));
 	solidNames_.push_back(name);
 }
 
@@ -325,7 +325,7 @@ bool pFlow::stlFile::write()const
 {
 	oFstream os(file_);
 	os.precision(8);
-	for(label i=0; i<size(); i++)
+	for(size_t i=0; i<size(); i++)
 	{
 		writeSolid(os, solids_[i], solidNames_[i]);
 	}
@@ -353,7 +353,7 @@ size_t pFlow::stlFile::size()const
 
 const pFlow::realx3x3Vector& pFlow::stlFile::solid
 (
-	label i
+	size_t i
 )const
 {
 	if(i >= size() )
@@ -369,7 +369,7 @@ const pFlow::realx3x3Vector& pFlow::stlFile::solid
 
 const pFlow::word& pFlow::stlFile::name
 (
-	label i
+	size_t i
 )const
 {
 	if(i >= size() )
diff --git a/src/phasicFlow/triSurface/stlFile.hpp b/utilities/Utilities/geometryPhasicFlow/stlWall/stlFile.hpp
similarity index 97%
rename from src/phasicFlow/triSurface/stlFile.hpp
rename to utilities/Utilities/geometryPhasicFlow/stlWall/stlFile.hpp
index 6f47fbd5..eea2f928 100644
--- a/src/phasicFlow/triSurface/stlFile.hpp
+++ b/utilities/Utilities/geometryPhasicFlow/stlWall/stlFile.hpp
@@ -106,10 +106,10 @@ public:
 		size_t size()const;
 		
 		// - vertecies of ith solid
-		const realx3x3Vector& solid(uint64 i)const;
+		const realx3x3Vector& solid(size_t i)const;
 
 		// - name of ith solid 
-		const word& name(uint64 i)const;
+		const word& name(size_t i)const;
 
 };
 
diff --git a/utilities/checkPhasicFlow/checkPhasicFlow.cpp b/utilities/checkPhasicFlow/checkPhasicFlow.cpp
index 3ffe4168..cbd36d03 100755
--- a/utilities/checkPhasicFlow/checkPhasicFlow.cpp
+++ b/utilities/checkPhasicFlow/checkPhasicFlow.cpp
@@ -20,6 +20,7 @@ Licence:
 
 #include "KokkosTypes.hpp"
 #include "systemControl.hpp"
+#include "localProcessors.hpp"
 #include "commandLine.hpp"
 
 
@@ -40,9 +41,11 @@ if(!cmds.parse(argc, argv)) return 0;
 #include "initialize.hpp"
 
 	output<<endl;
-	REPORT(1)<< "You are using "<<yellowText(cmds.productNameCopyright())<<endREPORT;
-  REPORT(1)<< yellowText(pFlow::floatingPointDescription())<<endREPORT;
-	
+	REPORT(1)<< "You are using "<<Yellow_Text(cmds.productNameCopyright())<<END_REPORT;
+	REPORT(1)<< Yellow_Text(pFlow::floatingPointDescription())<<END_REPORT;
+    REPORT(1)<< (pFlow::localProcessors::builtForMPI()?
+        "This is a built for MPI execution":
+        "This is not a build for MPI execution")<<END_REPORT;
 
 // this should be palced in each main 
 #include "finalize.hpp"
diff --git a/utilities/geometryPhasicFlow/geometryPhasicFlow.cpp b/utilities/geometryPhasicFlow/geometryPhasicFlow.cpp
index 4fa51770..8db0b791 100755
--- a/utilities/geometryPhasicFlow/geometryPhasicFlow.cpp
+++ b/utilities/geometryPhasicFlow/geometryPhasicFlow.cpp
@@ -26,12 +26,10 @@ Licence:
 //#include "readControlDict.hpp"
 
 
-using namespace pFlow;
-
 int main( int argc, char* argv[] )
 {
 
-	commandLine cmds(
+	pFlow::commandLine cmds(
     	"geometryPhasicFlow",
     	"Converts the supplied informaiton for sufraces in"
     	" geometryDict into PhasicFlow geometry data structure");
@@ -51,13 +49,13 @@ int main( int argc, char* argv[] )
 	#include "setProperty.hpp"
 
 	REPORT(0)<<"\nReading "<<"geometryDict"<<" . . ."<<END_REPORT;
-	auto geometryDict = fileDictionary(
-		objectFile
+	auto geometryDict = pFlow::fileDictionary(
+		pFlow::objectFile
 		(
 			"geometryDict",
 			Control.settings().path(),
-			objectFile::READ_ALWAYS,
-			objectFile::WRITE_NEVER
+			pFlow::objectFile::READ_ALWAYS,
+			pFlow::objectFile::WRITE_NEVER
 		),
 		nullptr
 	);
@@ -67,33 +65,33 @@ int main( int argc, char* argv[] )
 	auto wallsDictName = surfacesDict.dictionaryKeywords();
 
 	
-	word mSurfaceName = word("geometryPhasicFlow_")+word(triSurfaceFile__);
-	multiTriSurface surface
+    auto mSurfaceName = pFlow::word("geometryPhasicFlow_")+pFlow::word(pFlow::triSurfaceFile__);
+	pFlow::multiTriSurface surface
 	(
-		objectFile
+		pFlow::objectFile
 		(
 			mSurfaceName,
 			"",
-			objectFile::READ_NEVER,
-			objectFile::WRITE_ALWAYS
+			pFlow::objectFile::READ_NEVER,
+			pFlow::objectFile::WRITE_ALWAYS
 		),
 		nullptr
 	);
 	
-	wordVector materials;
-	wordList materialsList;
+	pFlow::wordVector materials;
+	pFlow::wordList materialsList;
 
-	wordVector motion;
-	wordList motionList;
+	pFlow::wordVector motion;
+	pFlow::wordList motionList;
 
 	for(auto& name:wallsDictName)
 	{
 		REPORT(1)<<"Creating wall "<<Green_Text(name)<<" from dictionary "<<surfacesDict.globalName() <<END_REPORT;
-		auto wallPtr = Wall::create( surfacesDict.subDict(name));
+		auto wallPtr = pFlow::Wall::create( surfacesDict.subDict(name));
 		auto& wall = wallPtr();
 		REPORT(1)<<"wall type is "<<Green_Text(wall.typeName())<<'\n'<<END_REPORT;		
 
-		realx3x3Vector trinalges(wall.name(), wall.triangles());
+		pFlow::realx3x3Vector trinalges(wall.name(), wall.triangles());
 		
 		surface.appendSurface(wall.name(), trinalges);
 		materials.push_back(wall.materialName());
@@ -108,7 +106,7 @@ int main( int argc, char* argv[] )
 	REPORT(1)<<"Selected wall motion components are "<<Cyan_Text(motionList)<<'\n'<<END_REPORT;
 		
 	REPORT(0)<< "\nCreating geometry . . ."<<END_REPORT;
-	auto geomPtr = geometry::create(
+	auto geomPtr = pFlow::geometry::create(
 		Control, 
 		proprties, 
 		surface, 
diff --git a/utilities/pFlowToVTK/pFlowToVTK.cpp b/utilities/pFlowToVTK/pFlowToVTK.cpp
index 6dfcae3c..8edbbc3c 100755
--- a/utilities/pFlowToVTK/pFlowToVTK.cpp
+++ b/utilities/pFlowToVTK/pFlowToVTK.cpp
@@ -30,18 +30,16 @@ Licence:
 //#include "readControlDict.hpp"
 
 
-using namespace pFlow;
-
 int main(int argc, char** argv )
 {
-	word outFolder = (pFlow::CWD()/word("VTK")).wordPath();
+	pFlow::word outFolder = (pFlow::CWD()/pFlow::word("VTK")).wordPath();
 
-	commandLine cmds(
+	pFlow::commandLine cmds(
 		"pFlowToVTK",
 		"Convrtes the saved pointField and geometry"
 		" date in time folders into vtk file format.");
 
-	wordVector times;
+	pFlow::wordVector times;
 		 
 	bool noGoem = false;
 	cmds.add_flag(
@@ -65,7 +63,7 @@ int main(int argc, char** argv )
 		separateSurfaces,
 		"use this when you want to have sub-surfaces in separate files");
 
-	wordVector fields;
+	pFlow::wordVector fields;
 	bool 			 allFields = true;
 	cmds.addOption("-f,--fields",
 		fields.vectorField(),
@@ -87,10 +85,10 @@ int main(int argc, char** argv )
 #include "initialize_Control.hpp"
 
 
-	timeFolder folders(Control);
-	auto destFolder = fileSystem(outFolder)/word(geometryFolder__);
-	auto destFolderField = fileSystem(outFolder);
-	wordList geomfiles{"triSurface"};
+	pFlow::timeFolder folders(Control);
+	auto destFolder = pFlow::fileSystem(outFolder)/pFlow::word(pFlow::geometryFolder__);
+	auto destFolderField = pFlow::fileSystem(outFolder);
+	pFlow::wordList geomfiles{"triSurface"};
 
 
 	if(cmds.count("--fields"))
@@ -98,7 +96,7 @@ int main(int argc, char** argv )
 		allFields = false;
 	}
 
-	realCombinedRange validRange;
+	pFlow::realCombinedRange validRange;
 	if( cmds.count("--time") )
 	{
 		if(!validRange.addRanges(times))
@@ -116,7 +114,7 @@ int main(int argc, char** argv )
 		Control.time().setTime(folders.time());
 		if( !validRange.isMember( folders.time() ) )continue;
 		
-		output<< "time: " << Cyan_Text( folders.time() )<<" s" <<endl;
+		pFlow::output<< "time: " << Cyan_Text( folders.time() )<<" s" <<pFlow::endl;
 		if(!noGoem)
 		{	
 			
@@ -155,13 +153,13 @@ int main(int argc, char** argv )
 			}	
 		}
 		
-		output<<endl;
+		pFlow::output<<pFlow::endl;
 
 	}
 	while( folders++ );
 
 
-	output<< "\nFinished successfully.\n";
+	pFlow::output<< "\nFinished successfully.\n";
 
 	
 // this should be palced in each main 
diff --git a/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp b/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp
index 5cdced7c..637fcf65 100755
--- a/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp
+++ b/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp
@@ -29,12 +29,10 @@ Licence:
 
 //#include "readControlDict.hpp"
 
-using namespace pFlow;
-
 int main( int argc, char* argv[] )
 {
 
-	commandLine cmds(
+	pFlow::commandLine cmds(
 		"createParticles",
 		"Read the dictionary createParticles and create particles"
 		" based on the two sub-dictionaries positionParticles and setFields.\n"
@@ -73,10 +71,10 @@ int main( int argc, char* argv[] )
 // this should be palced in each main 
 #include "initialize_Control.hpp"
 
-	fileDictionary cpDict("particlesDict", Control.settings().path());
+	pFlow::fileDictionary cpDict("particlesDict", Control.settings().path());
 	
 
-	uniquePtr<pointStructure> pStructPtr = nullptr;
+	pFlow::uniquePtr<pFlow::pointStructure> pStructPtr = nullptr;
 
 
 	if(!setOnly)
@@ -84,13 +82,13 @@ int main( int argc, char* argv[] )
 
 		// position particles based on the dict content 
 		REPORT(0)<< "Positioning points . . . \n"<<END_REPORT;
-		auto pointPosition = positionParticles::create(Control, cpDict.subDict("positionParticles"));
+		auto pointPosition = pFlow::positionParticles::create(Control, cpDict.subDict("positionParticles"));
 
-		fileSystem pStructPath = Control.time().path()+pointStructureFile__;
+		pFlow::fileSystem pStructPath = Control.time().path()+pFlow::pointStructureFile__;
 
 		auto finalPos = pointPosition().getFinalPosition();
 
-		pStructPtr = makeUnique<pointStructure>(Control, finalPos);
+		pStructPtr = pFlow::makeUnique<pFlow::pointStructure>(Control, finalPos);
 
 
 		REPORT(1)<< "Created pStruct with "<< pStructPtr().size() << " points and capacity "<<
@@ -100,11 +98,11 @@ int main( int argc, char* argv[] )
     else
 	{
 		// read the content of pStruct from 0/pStructure
-		pStructPtr = makeUnique<pointStructure>(Control);
+		pStructPtr = pFlow::makeUnique<pFlow::pointStructure>(Control);
 
 	}
 
-	List<uniquePtr<IOobject>> allObjects;
+	pFlow::List<pFlow::uniquePtr<pFlow::IOobject>> allObjects;
 
 	if(!positionOnly)
 	{
@@ -113,7 +111,7 @@ int main( int argc, char* argv[] )
 
 		auto& sfDict = cpDict.subDict("setFields");
 
-		setFieldList defValueList(sfDict.subDict("defaultValue"));
+		pFlow::setFieldList defValueList(sfDict.subDict("defaultValue"));
 
 		for(auto& sfEntry: defValueList)
 		{
@@ -128,7 +126,7 @@ int main( int argc, char* argv[] )
 			}
 		}
 		
-		output<<endl;
+		pFlow::output<<pFlow::endl;
 
 		auto& selectorsDict = sfDict.subDict("selectors");
 
@@ -145,45 +143,45 @@ int main( int argc, char* argv[] )
 				ERR<<"\n error occured in setting selector. \n"<<END_ERR;
 				return 1;
 			}
-			output<<endl;
+			pFlow::output<<pFlow::endl;
 		}
 	}
 
 	Control.clearIncludeExclude();
 	Control.addExclude("shapeName");
 
-	uint64PointField_H shapeHash
+	pFlow::uint64PointField_H shapeHash
 	(
-		objectFile
+		pFlow::objectFile
 		(
 			"shapeHash",
 			"",
-			objectFile::READ_NEVER,
-			objectFile::WRITE_ALWAYS
+			pFlow::objectFile::READ_NEVER,
+			pFlow::objectFile::WRITE_ALWAYS
 		),
 		pStructPtr(),
 		0u
 	);
 
-	uint32PointField_H shapeIndex
+	pFlow::uint32PointField_H shapeIndex
 	(
-		objectFile
+		pFlow::objectFile
 		(
 			"shapeIndex",
 			"",
-			objectFile::READ_NEVER,
-			objectFile::WRITE_ALWAYS
+			pFlow::objectFile::READ_NEVER,
+			pFlow::objectFile::WRITE_ALWAYS
 		),
 		pStructPtr(),
 		0u
 	);
 
-	baseShapeNames shapes(
-		shapeFile__,
+	pFlow::baseShapeNames shapes(
+		pFlow::shapeFile__,
 		&Control.caseSetup()
 	);
 	
-	auto& shapeName = Control.time().template lookupObject<wordPointField_H>("shapeName");
+	auto& shapeName = Control.time().template lookupObject<pFlow::wordPointField_H>("shapeName");
 
 	REPORT(0)<< "Converting shapeName field to shapeIndex field"<<END_REPORT;
 
@@ -196,7 +194,7 @@ int main( int argc, char* argv[] )
 
 	ForAll(i, shapeHash)
 	{
-		if(uint32 index; shapes.shapeNameToIndex(shapeName_D[i], index))
+		if(pFlow::uint32 index; shapes.shapeNameToIndex(shapeName_D[i], index))
 		{
 			shapeHash_D[i] = shapes.hashes()[index];
 			shapeIndex_D[i] = index;
@@ -206,7 +204,7 @@ int main( int argc, char* argv[] )
 			fatalErrorInFunction<<"Found shape name "<< Yellow_Text(shapeName_D[i])<<
 			"in shapeName field. But the list of shape names in file "<< 
 			shapes.globalName()<<" is : \n"<<
-			shapes.shapeNames()<<endl;
+			shapes.shapeNames()<<pFlow::endl;
 		}
 	}