Merge pull request #109 from PhasicFlow/develop

develop branch update for changes in MPI branch for data transfer stage

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()

src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp

@@ -160,6 +160,7 @@ public:
 		csPairContainerType& pwPairs,
 		bool force = false)override
 	{
+		if(i==0u)
 			Particles().boundingSphere().updateBoundaries(DataDirection::SlaveToMaster);
 		return csBoundaries_[i].broadSearch(
 			iter, 
@@ -184,6 +185,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();
+	}
+	
 };
 
 }

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,

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,

src/Interaction/contactSearch/methods/cellBased/NBS/NBS.hpp

@@ -135,6 +135,16 @@ public:
 			return 1;
 		}
 		
+		real sizeRatio()const
+		{
+			return sizeRatio_;
+		}
+
+		real cellExtent()const
+		{
+			return cellExtent_;
+		}
+
 		auto getCellIterator([[maybe_unused]] uint32 lvl)const
 		{
 			return NBSLevel0_.getCellIterator();

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();
+	}
+
 };
 
 

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) )
 	{

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>

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,

src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp

@@ -31,6 +31,9 @@ pFlow::periodicBoundarySphereInteraction<cFM, gMM>::periodicBoundarySphereIntera
 	if(boundary.thisBoundaryIndex()%2==1)
     {
     	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;
 }

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;
 	
 };
 

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++)
+	{
+		auto& BI =boundaryInteraction_[i];
+		if(BI.ppPairsAllocated())
 		{
 			if( !contactSearch_().boundaryBroadSearch(
 				i,
 				iter,
 				t,
 				dt,
-			boundaryInteraction_[i].ppPairs(),
-			boundaryInteraction_[i].pwPairs()))
+				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();
+	{
+	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++)
 		{
-		boundaryInteraction_[i].sphereSphereInteraction(
+			if(requireStep[i])
+			{
+				requireStep[i] = boundaryInteraction_[i].sphereSphereInteraction(
 					dt, 
-			this->forceModel_());
+					this->forceModel_(),
+					step
+				);
+			}		
+		}
+		step++;
 	}
 	boundaryInteractionTimer_.end();
+	}
+	
 	return true;
 }
 

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();

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;
 }

src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp

@@ -79,7 +79,7 @@ private:
 	/// timer for integration computations (correction step)
 	Timer                  intCorrectTimer_;
 
-	
+	Timer  				   fieldUpdateTimer_;
 
 private:
 	bool initializeParticles();

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()

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_;
 

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;
 }
 

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 };
 

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

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)
 	{
@@ -833,8 +894,6 @@ bool pFlow::VectorSingle<T,MemorySpace>::reorderItems(uint32IndexContainer indic
 
 	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();
 	}
 
+    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;
 }

src/phasicFlow/containers/VectorHD/VectorSingle.hpp

@@ -159,6 +159,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.

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();
+		}
+
 };
 
 

src/phasicFlow/containers/pointField/boundaryField/boundaryField/boundaryField.hpp

@@ -121,6 +121,17 @@ public:
 		return true;
 	}
 
+	inline
+	InternalFieldType& internal()
+	{
+		return internal_;
+	}
+
+	inline
+	const InternalFieldType& internal()const
+	{
+		return internal_;
+	}
 
 	FieldAccessType thisField()const
 	{
@@ -144,11 +155,6 @@ public:
 	virtual
 	const ProcVectorType& neighborProcField()const;
 
-	void fill(const std::any& val)override
-	{
-		return;
-	}
-
 	virtual
 	void fill(const T& val)
 	{

src/phasicFlow/containers/pointField/boundaryField/generalBoundary/generalBoundary.hpp

@@ -148,10 +148,6 @@ public:
 
     const Time& time()const;
 	
-	virtual
-	void fill(const std::any& val)=0;
-
-	
 };
 
 } // pFlow

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;
 }
 

src/phasicFlow/containers/pointField/internalField/internalField.hpp

@@ -73,6 +73,8 @@ protected:
 
 	bool insert(const anyList& varList);
 
+	bool rearrange(const anyList& varList);
+
 public:
 
 	internalField(

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:

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_)

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_;

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 			

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),

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

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];

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())
 	{

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;
 }

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);
 

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;
 
 
 };

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();
 }
 
 
@@ -58,6 +61,8 @@ bool pFlow::boundaryPeriodic::beforeIteration(
 		return true;
 	}
 
+	if( !boundaryListUpdate(iterNum))return true;
+	
 	uint32 s = size();
 	uint32Vector_D transferFlags("transferFlags",s+1, s+1, RESERVE()); 
 	transferFlags.fill(0u);
@@ -92,7 +97,7 @@ bool pFlow::boundaryPeriodic::beforeIteration(
 	// to obtain the transfer vector 
 	realx3 transferVec = displacementVectroToMirror();
 	
-	return transferPoints
+	return transferPointsToMirror
 	(
 		numTransfered,
 		transferFlags, 

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:
 

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,21 +173,6 @@ 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(
@@ -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()))
 			);
 	}
 };	

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 <<

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"),
+    : subscriber("internalPoints"),
       pointPosition_("internalPoints", "internalPoints", initialCapacity_, 0, RESERVE()),
-	pFlagsD_(initialCapacity_, 0 , 0)
+      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

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;
@@ -236,6 +238,8 @@ public:
 		
 		bool insertPoints(const realx3Vector& points, anyList& varList);
 
+		bool insertPointsOnly(const realx3Vector_D& points, message& msg, anyList& varList);
+		
 	
 	//// - IO operations 
 		

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;
 
@@ -285,6 +307,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
 	/// and return number of deleted points 
 	/// @param validBox the box whose inside is valid 

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

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)
+	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;
 }

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);
+	const ViewType1D<realx3>& pos, 
+	const pFlagTypeDevice& flag);
 
 
 INLINE_FUNCTION_HD

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);
+}

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

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))
+    {
+        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;
 }
 

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);