periodic boundary condition is added to pointStructure

src/phasicFlow/CMakeLists.txt

@@ -71,11 +71,14 @@ structuredData/plane/plane.cpp
 structuredData/domain/domain.cpp
 structuredData/domain/simulationDomain.cpp
 structuredData/domain/regularSimulationDomain.cpp
+structuredData/pointStructure/internalPointsKernels.cpp
 structuredData/pointStructure/internalPoints.cpp
 structuredData/pointStructure/pointStructure.cpp
 structuredData/boundaries/boundaryBase/boundaryBase.cpp
+structuredData/boundaries/boundaryBase/boundaryBaseKernels.cpp
 structuredData/boundaries/boundaryExit/boundaryExit.cpp
 structuredData/boundaries/boundaryNone/boundaryNone.cpp
+structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp
 structuredData/boundaries/boundaryList.cpp
 structuredData/pointStructure/selectors/pStructSelector/pStructSelector.cpp
 structuredData/pointStructure/selectors/selectBox/selectBox.cpp

src/phasicFlow/Kokkos/KokkosTypes.hpp

@@ -39,13 +39,6 @@ Licence:
 namespace pFlow
 {
 
-///********DEP
-class DeviceSide{};
-class HostSide{};
-
-template<typename side>
-struct selectSide{};
-/*********///
 
 /// Host memory space 
 using HostSpace 	= Kokkos::HostSpace;
@@ -71,6 +64,20 @@ using DefaultHostExecutionSpace 	= Kokkos::DefaultHostExecutionSpace;
 using DefaultExecutionSpace 		= Kokkos::DefaultExecutionSpace;
 
 
+using deviceRPolicyStatic = 
+	Kokkos::RangePolicy<
+			Kokkos::DefaultExecutionSpace,
+			Kokkos::Schedule<Kokkos::Static>,
+			Kokkos::IndexType<pFlow::uint32> >;
+
+
+using hostRPolicyStatic = 
+	Kokkos::RangePolicy<
+			Kokkos::DefaultExecutionSpace,
+			Kokkos::Schedule<Kokkos::Static>,
+			Kokkos::IndexType<pFlow::uint32> >;
+
+
 /// Pair of two variables
 template<typename T1, typename T2>
 using Pair 	= Kokkos::pair<T1,T2>;

src/phasicFlow/containers/VectorHD/VectorSingle.cpp

@@ -232,7 +232,14 @@ pFlow::VectorSingle<T,MemorySpace>::VectorField()const
 
 template<typename T, typename MemorySpace>
 INLINE_FUNCTION_H 
-auto pFlow::VectorSingle<T,MemorySpace>::deviceViewAll() const
+const auto& pFlow::VectorSingle<T,MemorySpace>::deviceViewAll() const
+{
+    return view_;
+}
+
+template<typename T, typename MemorySpace>
+INLINE_FUNCTION_H 
+auto& pFlow::VectorSingle<T,MemorySpace>::deviceViewAll()
 {
     return view_;
 }

src/phasicFlow/containers/VectorHD/VectorSingle.hpp

@@ -191,7 +191,11 @@ public:
 
 		/// Device view range [0,capcity)
 		INLINE_FUNCTION_H 
-		auto deviceViewAll() const;
+		auto& deviceViewAll();
+
+		/// Device view range [0,capcity)
+		INLINE_FUNCTION_H 
+		const auto& deviceViewAll() const;
 
 		///  Device view range [0, size)
 		INLINE_FUNCTION_H 

src/phasicFlow/containers/pointField/periodicBoundaryField.cpp

@@ -0,0 +1,29 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+template<class T, class MemorySpace>
+    pFlow::periodicBoundaryField<T, MemorySpace>::periodicBoundaryField
+(
+	const boundaryBase& boundary, 
+	InternalFieldType& internal
+)
+:
+    BoundaryFieldType(boundary, internal)
+{}

src/phasicFlow/containers/pointField/periodicBoundaryField.hpp

@@ -0,0 +1,83 @@
+/*------------------------------- 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 __periodicBoundaryField_hpp__
+#define __periodicBoundaryField_hpp__
+
+#include "boundaryField.hpp"
+
+namespace pFlow
+{
+
+template< class T, class MemorySpace = void>
+class periodicBoundaryField
+:
+    public boundaryField<T, MemorySpace> 
+{
+public:
+	
+    using periodicBoundaryFieldType = periodicBoundaryField<T, MemorySpace>;
+
+	using BoundaryFieldType = boundaryField<T, MemorySpace>;
+
+	using InternalFieldType = typename BoundaryFieldType::InternalFieldType;
+
+	using memory_space 		= typename BoundaryFieldType::memory_space;
+
+	using execution_space 	= typename BoundaryFieldType::execution_space;
+
+   
+
+public:
+
+	TypeInfo("boundaryField<periodic>");
+
+	periodicBoundaryField(
+		const boundaryBase& boundary, 
+		InternalFieldType& internal);
+		
+
+	add_vCtor
+	(
+		BoundaryFieldType,
+		periodicBoundaryFieldType,
+		boundaryBase
+	);
+
+
+	bool hearChanges
+	(
+		real t,
+		real dt,
+		uint32 iter,
+		const message& msg, 
+    	const anyList& varList
+	) override
+    {
+		notImplementedFunction;
+		return false;
+	}
+
+};
+
+}
+
+#include "periodicBoundaryField.cpp"
+
+#endif

src/phasicFlow/containers/pointField/pointFields.cpp

@@ -21,69 +21,73 @@ Licence:
 
 #include "pointFields.hpp"
 #include "createBoundaryFields.hpp"
+#include "periodicBoundaryField.hpp"
 
+#define createAllBoundary(DataType, MemorySpaceType)          \
+    template class pFlow::exitBoundaryField<DataType, MemorySpaceType>; \
+    template class pFlow::periodicBoundaryField<DataType, MemorySpaceType>;
 
 template class pFlow::pointField<pFlow::int8, pFlow::HostSpace>; 
 createBaseBoundary(pFlow::int8, pFlow::HostSpace);
-createBoundary(pFlow::int8, pFlow::HostSpace, exit);
+createAllBoundary(pFlow::int8, pFlow::HostSpace);
 
 template class pFlow::pointField<pFlow::int8>; 
 createBaseBoundary(pFlow::int8, void);
-createBoundary(pFlow::int8, void, exit);
+createAllBoundary(pFlow::int8, void);
 
 
 template class pFlow::pointField<pFlow::uint8, pFlow::HostSpace>; 
 createBaseBoundary(pFlow::uint8, pFlow::HostSpace);
-createBoundary(pFlow::uint8, pFlow::HostSpace, exit);
+createAllBoundary(pFlow::uint8, pFlow::HostSpace);
 
 template class pFlow::pointField<pFlow::uint8>; 
 createBaseBoundary(pFlow::uint8, void);
-createBoundary(pFlow::uint8, void, exit);
+createAllBoundary(pFlow::uint8, void);
 
 template class pFlow::pointField<pFlow::int32, pFlow::HostSpace>; 
 createBaseBoundary(pFlow::int32, pFlow::HostSpace);
-createBoundary(pFlow::int32, pFlow::HostSpace, exit);
+createAllBoundary(pFlow::int32, pFlow::HostSpace);
 
 template class pFlow::pointField<pFlow::int32>; 
 createBaseBoundary(pFlow::int32, void);
-createBoundary(pFlow::int32, void, exit);
+createAllBoundary(pFlow::int32, void);
 
 template class pFlow::pointField<pFlow::uint32, pFlow::HostSpace>; 
 createBaseBoundary(pFlow::uint32, pFlow::HostSpace);
-createBoundary(pFlow::uint32, pFlow::HostSpace, exit);
+createAllBoundary(pFlow::uint32, pFlow::HostSpace);
 
 template class pFlow::pointField<pFlow::uint32>; 
 createBaseBoundary(pFlow::uint32, void);
-createBoundary(pFlow::uint32, void, exit);
+createAllBoundary(pFlow::uint32, void);
 
 template class pFlow::pointField<pFlow::real, pFlow::HostSpace>; 
 createBaseBoundary(pFlow::real, pFlow::HostSpace);
-createBoundary(pFlow::real, pFlow::HostSpace, exit);
+createAllBoundary(pFlow::real, pFlow::HostSpace);
 
 
 template class pFlow::pointField<pFlow::real>; 
 createBaseBoundary(pFlow::real, void);
-createBoundary(pFlow::real, void, exit);
+createAllBoundary(pFlow::real, void);
 
 template class pFlow::pointField<pFlow::realx3, pFlow::HostSpace>; 
 createBaseBoundary(pFlow::realx3, pFlow::HostSpace);
-createBoundary(pFlow::realx3, pFlow::HostSpace, exit);
+createAllBoundary(pFlow::realx3, pFlow::HostSpace);
 
 
 template class pFlow::pointField<pFlow::realx3>; 
 createBaseBoundary(pFlow::realx3, void);
-createBoundary(pFlow::realx3, void, exit);
+createAllBoundary(pFlow::realx3, void);
 
 
 template class pFlow::pointField<pFlow::realx4, pFlow::HostSpace>; 
 createBaseBoundary(pFlow::realx4, pFlow::HostSpace);
-createBoundary(pFlow::realx4, pFlow::HostSpace, exit);
+createAllBoundary(pFlow::realx4, pFlow::HostSpace);
 
 
 template class pFlow::pointField<pFlow::realx4>; 
 createBaseBoundary(pFlow::realx4, void);
-createBoundary(pFlow::realx4, void, exit);
+createAllBoundary(pFlow::realx4, void);
 
 template class pFlow::pointField<pFlow::word, pFlow::HostSpace>; 
 createBaseBoundary(pFlow::word, pFlow::HostSpace);
-createBoundary(pFlow::word, pFlow::HostSpace, exit);
+createAllBoundary(pFlow::word, pFlow::HostSpace);

src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp

@@ -22,6 +22,8 @@ Licence:
 #include "dictionary.hpp"
 #include "internalPoints.hpp"
 
+#include "boundaryBaseKernels.hpp"
+
 /*pFlow::boundaryBase::boundaryBase
 (
 	const plane& 		bplane,
@@ -48,21 +50,142 @@ void pFlow::boundaryBase::setNewIndices
 {
 	auto newSize = newIndices.size();
 	setSize(static_cast<uint32>(newSize));
-	copy(indexList_.deviceView(), newIndices);
+	if(newSize>0u)
+	{
+		copy(indexList_.deviceView(), newIndices);
+	}
 }
 
-pFlow::boundaryBase::boundaryBase(
+void pFlow::boundaryBase::appendNewIndices
+(
+	deviceViewType1D<uint32> newIndices
+)
+{
+	auto s = static_cast<uint32>(newIndices.size());
+	if(s == 0) return;
+
+	uint32 oldS = size();
+	uint32 newSize = oldS + s; 
+
+	setSize(newSize);
+	auto appendView = Kokkos::subview(
+		indexList_.deviceViewAll(),
+		Kokkos::make_pair<uint32>(oldS, newSize));
+	copy(appendView, newIndices);
+	
+	// TODO: notify observers about this change 
+
+	// the index list should be sorted  
+	//sort(indexList_.deviceViewAll(), 0, newSize);
+}
+
+
+
+bool pFlow::boundaryBase::removeIndices
+(
+	uint32 numRemove,
+	deviceViewType1D<uint32> removeMask
+)
+{
+    if(removeMask.size() != size()+1 )
+	{
+		fatalErrorInFunction<<"size mismatch"<<endl;
+		return false;
+	}
+
+	deviceViewType1D<uint32> removeIndices("removeIndices", 1);
+	deviceViewType1D<uint32> keepIndices("keepIndices",1);
+
+	pFlow::boundaryBaseKernels::createRemoveKeepIndices
+	(
+		indexList_.deviceView(),
+		numRemove,
+		removeMask,
+		removeIndices,
+		keepIndices
+	);
+	
+	if(!internal_.deletePoints(removeIndices))
+	{
+		fatalErrorInFunction<<
+		"error in deleting points from boundary "<< name()<<endl;
+		return false;
+	}
+
+	setNewIndices(keepIndices);
+
+	//TODO: notify observers about changes 
+
+	return false;
+}
+
+bool pFlow::boundaryBase::transferPoints
+(
+	uint32 numTransfer, 
+	deviceViewType1D<uint32> transferMask, 
+	uint32 transferBoundaryIndex, 
+	realx3 transferVector
+)
+{
+	if(transferMask.size() != size()+1 )
+	{
+		fatalErrorInFunction<<"size mismatch"<<endl;
+		return false;
+	}
+
+	deviceViewType1D<uint32> transferIndices("transferIndices",1);
+	deviceViewType1D<uint32> keepIndices("keepIndices",1);
+
+	pFlow::boundaryBaseKernels::createRemoveKeepIndices
+	(
+		indexList_.deviceView(),
+		numTransfer,
+		transferMask,
+		transferIndices,
+		keepIndices
+	);
+
+	// third, remove the indices from this list 
+	setNewIndices(keepIndices);
+
+	// first, change the flags in the internalPoints
+	if( !internal_.changePointsFlag(
+		transferIndices, 
+		transferBoundaryIndex) ) 
+	{
+		return false;
+	}
+
+	// second, change the position of points 
+	if(!internal_.changePointsPoisition(
+		transferIndices,
+		transferVector))
+	{
+		return false;
+	}
+
+	// fourth, add the indices to the mirror boundary 
+	internal_.boundary(transferBoundaryIndex).
+		appendNewIndices(transferIndices);
+	
+    return true;
+}
+
+pFlow::boundaryBase::boundaryBase
+(
     const dictionary &dict,
     const plane &bplane,
-    internalPoints &internal)
+    internalPoints &internal
-    : subscriber(dict.name()),
+)
-      boundaryPlane_(bplane),
+: 
-      indexList_(groupNames(dict.name(), "indexList")),
+	subscriber(dict.name()),
-      neighborLength_(dict.getVal<real>("neighborLength")),
+	boundaryPlane_(bplane),
-      internal_(internal),
+	indexList_(groupNames(dict.name(), "indexList")),
-      mirrorProcessoNo_(dict.getVal<uint32>("mirrorProcessorNo")),
+	neighborLength_(dict.getVal<real>("neighborLength")),
-      name_(dict.name()),
+	internal_(internal),
-      type_(dict.getVal<word>("type"))
+	mirrorProcessoNo_(dict.getVal<uint32>("mirrorProcessorNo")),
+	name_(dict.name()),
+	type_(dict.getVal<word>("type"))
 {
 }
 

src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp

@@ -69,6 +69,19 @@ protected:
 
 	void setNewIndices(deviceViewType1D<uint32> newIndices);
 
+	void appendNewIndices(deviceViewType1D<uint32> newIndices);
+	
+	bool removeIndices(
+		uint32 numRemove, 
+		deviceViewType1D<uint32> removeMask);
+	
+	bool transferPoints(
+		uint32 numTransfer, 
+		deviceViewType1D<uint32> transferMask,
+		uint32 transferBoundaryIndex,
+		realx3 transferVector);
+	
+	
 public:
 
 	TypeInfo("boundaryBase");
@@ -102,12 +115,17 @@ public:
 		(dict, bplane, internal)
 	);
 
-	inline 
+	virtual 
-	auto neighborLength()const
+	real neighborLength()const
 	{
 		return neighborLength_;
 	}
 
+	virtual 
+	realx3 boundaryExtensionLength()const
+	{
+		return {0,0,0};
+	}
 
 	const word& type()const
 	{
@@ -129,6 +147,7 @@ public:
 		return indexList_.size();
 	}
 
+	virtual 
 	const plane& boundaryPlane()const
 	{
 		return boundaryPlane_;
@@ -168,11 +187,6 @@ public:
 		return indexList_;
 	}
 
-	/*auto& indexList()
-	{
-		return indexList_;
-	}*/
-
     pointFieldAccessType thisPoints();
 
     virtual

src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBaseKernels.cpp

@@ -0,0 +1,90 @@
+/*------------------------------- 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 "boundaryBaseKernels.hpp"
+
+void pFlow::boundaryBaseKernels::createRemoveKeepLists
+(
+    uint32 numTotal, 
+    uint32 numRemove, 
+    deviceViewType1D<uint32> removeMask, 
+    deviceViewType1D<uint32>& removeList, 
+    deviceViewType1D<uint32>& keepList
+)
+{
+    uint32 numKeep = numTotal - numRemove;
+    deviceViewType1D<uint32> rList("rList",numRemove);
+    deviceViewType1D<uint32> kList("kList",numKeep);
+
+    exclusiveScan(removeMask, 0u, numTotal+1, removeMask, 0u);
+
+    Kokkos::parallel_for
+    (
+        "pFlow::boundaryBaseKernels::createRemoveKeepLists",
+        deviceRPolicyStatic(0, numTotal),
+        LAMBDA_HD(uint32 i)
+		{
+			if(removeMask(i)!= removeMask(i+1))
+				rList(removeMask(i)) = i;
+            else
+                kList(i-removeMask(i)) = i;
+		}
+    );
+    Kokkos::fence();
+
+    removeList = rList;
+    keepList = kList;
+
+}
+
+
+void pFlow::boundaryBaseKernels::createRemoveKeepIndices
+(
+    deviceViewType1D<uint32> indices,
+    uint32 numRemove,
+    deviceViewType1D<uint32> removeMask, 
+    deviceViewType1D<uint32>& removeIndices, 
+    deviceViewType1D<uint32>& keepIndices
+)
+{
+    uint32 numTotal = indices.size();
+    uint32 numKeep = numTotal - numRemove;
+    deviceViewType1D<uint32> rIndices("rIndices",numRemove);
+    deviceViewType1D<uint32> kIndices("kIndices",numKeep);
+
+    exclusiveScan(removeMask, 0u, numTotal+1, removeMask, 0u);
+
+    Kokkos::parallel_for
+    (
+        "pFlow::boundaryBaseKernels::createRemoveKeepLists",
+        deviceRPolicyStatic(0, numTotal),
+        LAMBDA_HD(uint32 i)
+		{
+			if(removeMask(i)!= removeMask(i+1))
+				rIndices(removeMask(i)) = indices(i);
+            else
+                kIndices(i-removeMask(i)) = indices(i);
+		}
+    );
+    Kokkos::fence();
+
+    removeIndices = rIndices;
+    keepIndices = kIndices;
+}

src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBaseKernels.hpp

@@ -0,0 +1,46 @@
+/*------------------------------- 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 __boundaryBaseKernels_hpp__
+#define __boundaryBaseKernels_hpp__
+
+#include "phasicFlowKokkos.hpp"
+
+namespace pFlow::boundaryBaseKernels
+{
+
+void createRemoveKeepLists(
+    uint32 numTotal, 
+    uint32 numRemove, 
+    deviceViewType1D<uint32> removeMask, 
+    deviceViewType1D<uint32>& removeList, 
+    deviceViewType1D<uint32>& keepList);
+
+void createRemoveKeepIndices(
+    deviceViewType1D<uint32> indices,
+    uint32 numRemove,
+    deviceViewType1D<uint32> removeMask, 
+    deviceViewType1D<uint32>& removeIndices, 
+    deviceViewType1D<uint32>& keepIndices);
+
+
+}
+
+#endif //__boundaryBaseKernels_hpp__

src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp

@@ -52,106 +52,36 @@ bool pFlow::boundaryExit::beforeIteration
 	}
 
 	uint32 s = size();
-	deviceViewType1D<uint32> delFlags("delFlags",s+1); 
+	deviceViewType1D<uint32> deleteFlags("deleteFlags",s+1); 
-	deviceViewType1D<uint32> keepFlags("keepFlags", s+1);
+	fill(deleteFlags, 0, s+1, 0u);
-	fill(delFlags, 0, s+1, 0u);
-	fill(keepFlags, 0, s+1, 0u);
-
-	using policy = Kokkos::RangePolicy<
-			pFlow::DefaultExecutionSpace,
-			Kokkos::Schedule<Kokkos::Static>,
-			Kokkos::IndexType<pFlow::uint32>>;
 
 	auto points = thisPoints();
-	uint32 numDeleted = 0;
 	auto p = boundaryPlane().infPlane();
-	
+
+	uint32 numDeleted = 0;	
 
 	Kokkos::parallel_reduce
 	(
 		"boundaryExit::beforeIteration",
-		policy(0,s),
+		deviceRPolicyStatic(0,s),
 		LAMBDA_HD(uint32 i, uint32& delToUpdate)
 		{
 			if(p.pointInNegativeSide(points(i)))
 			{
-				delFlags(i)=1;
+				deleteFlags(i)=1;
 				delToUpdate++;
 			}
-			else
-			{
-				keepFlags(i) = 1;
-			}
 		}, 
 		numDeleted
 	);
-	
+		
-	
 	// no point is deleted
 	if(numDeleted == 0 )
 	{
 		return true;
 	}
 
-	exclusiveScan(delFlags, 0u, s+1, delFlags, 0u);
+	return this->removeIndices(numDeleted, deleteFlags);
-	exclusiveScan(keepFlags, 0u, s+1, keepFlags, 0u);
-	
-	deviceViewType1D<uint32> deleteList("deleteList", numDeleted);
-
-	Kokkos::parallel_for
-	(
-		"boundaryExit::parllel_for",
-		policy(0, size()),
-		LAMBDA_HD(uint32 i)
-		{
-			if(delFlags(i)!= delFlags(i+1))
-				deleteList(delFlags(i)) = i;
-		}	
-	);
-	Kokkos::fence();
-
-	deviceScatteredFieldAccess<uint32> deleteIndices(
-			numDeleted, 
-			deleteList,
-			indexList().deviceViewAll());
-	
-	// tell internal to remove these points from its list 
-	if(!internal().deletePoints(deleteIndices))
-	{
-		fatalErrorInFunction<<
-		"error in deleting points from boundary "<< name()<<endl;
-		return false;
-	}
-
-	// delete these indices from your list 
-	if(numDeleted == s )
-	{
-		setSize(0u);
-	}
-	else
-	{
-		uint32 newSize = s-numDeleted;
-		deviceViewType1D<uint32> newIndices("newIndices", newSize);
-		auto oldIndices = indexList().deviceViewAll();
-
-		Kokkos::parallel_for(
-			"fillIndices",
-			policy(0,s),
-			LAMBDA_HD(uint32 i)
-			{
-				if(keepFlags(i)!= keepFlags(i+1))
-				{
-					newIndices(keepFlags(i)) = oldIndices(i);
-				}
-			}
-		);
-		Kokkos::fence();
-		setNewIndices(newIndices);
-	}
-
-	//TODO: notify observers about changes 
-	
-	return true;
 }
 
 bool pFlow::boundaryExit::iterate

src/phasicFlow/structuredData/boundaries/boundaryList.cpp

@@ -40,9 +40,19 @@ bool pFlow::boundaryList::updateLists()
 	dist[4] = boundary(4).neighborLength();
 	dist[5] = boundary(5).neighborLength();
 
+	realx3 lowerExt = 
+		boundary(0).boundaryExtensionLength() + 
+		boundary(2).boundaryExtensionLength() + 
+		boundary(4).boundaryExtensionLength();
 	
+	realx3 upperExt = 
+		boundary(1).boundaryExtensionLength()+
+		boundary(3).boundaryExtensionLength()+
+		boundary(5).boundaryExtensionLength();
+
+	auto extDomain = pStruct_.simDomain().extendThisDomain(lowerExt, upperExt);
 	pStruct_.updateFlag(
-			pStruct_.simDomain().thisDomain(),
+			extDomain,
 			dist);
 	const auto& maskD = pStruct_.activePointsMaskDevice();
 	boundary(0).setSize( maskD.leftSize() );

src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp

@@ -0,0 +1,125 @@
+/*------------------------------- 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 "boundaryPeriodic.hpp"
+#include "internalPoints.hpp"
+#include "dictionary.hpp"
+
+
+pFlow::boundaryPeriodic::boundaryPeriodic
+(
+	const dictionary& dict,
+	const plane&      bplane,
+	internalPoints&   internal
+)
+:
+	boundaryBase(dict, bplane, internal),
+	mirrorBoundaryIndex_(dict.getVal<uint32>("mirrorBoundaryIndex"))
+{
+	extendedPlane_ = boundaryBase::boundaryPlane().parallelPlane(-boundaryBase::neighborLength());
+}
+
+pFlow::real pFlow::boundaryPeriodic::neighborLength() const
+{
+    return 2.0*boundaryBase::neighborLength();
+}
+
+pFlow::realx3 pFlow::boundaryPeriodic::boundaryExtensionLength() const
+{
+    return -neighborLength()*boundaryBase::boundaryPlane().normal();
+}
+
+const pFlow::plane &pFlow::boundaryPeriodic::boundaryPlane() const
+{
+    return extendedPlane_;
+}
+
+bool pFlow::boundaryPeriodic::beforeIteration(
+    uint32 iterNum,
+    real t,
+    real dt)
+{
+	// nothing have to be done
+	if(empty())
+	{
+		return true;
+	}
+
+	uint32 s = size();
+	deviceViewType1D<uint32> transferFlags("transferFlags",s+1); 
+	fill(transferFlags, 0, s+1, 0u);
+	
+	auto points = thisPoints();
+	auto p = boundaryPlane().infPlane();
+	uint32 numTransfered = 0;
+
+	Kokkos::parallel_reduce
+	(
+		"boundaryPeriodic::beforeIteration",
+		deviceRPolicyStatic(0u,s),
+		LAMBDA_HD(uint32 i, uint32& trnasToUpdate)
+		{
+			if(p.pointInNegativeSide(points(i)))
+			{
+				transferFlags(i)=1;
+				trnasToUpdate++;
+			}
+		}, 
+		numTransfered
+	);
+	
+	// no point to be transfered 
+	if(numTransfered == 0 )
+	{
+		return true;
+	}
+	
+	// to obtain the transfer vector 
+	const auto& thisP = boundaryPlane();
+	const auto& mirrorP = internal().boundary(mirrorBoundaryIndex_).boundaryPlane();
+	realx3 transferVec = thisP.normal()*(thisP.d() + mirrorP.d());
+
+	return transferPoints
+	(
+		numTransfered,
+		transferFlags, 
+		mirrorBoundaryIndex_, 
+		transferVec
+	);
+
+}
+
+bool pFlow::boundaryPeriodic::iterate
+(
+	uint32 iterNum, 
+	real t
+)
+{
+	return true;
+}
+
+bool pFlow::boundaryPeriodic::afterIteration
+(
+	uint32 iterNum, 
+	real t
+)
+{
+	return true;
+}

src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp

@@ -0,0 +1,76 @@
+/*------------------------------- 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 __boundaryPeriodic_hpp__
+#define __boundaryPeriodic_hpp__
+
+
+#include "boundaryBase.hpp"
+
+namespace pFlow
+{
+
+class boundaryPeriodic
+:
+ 	public boundaryBase
+{
+private:
+
+	uint32 mirrorBoundaryIndex_;
+
+	plane extendedPlane_;
+
+public:
+
+	TypeInfo("boundary<periodic>");
+
+	boundaryPeriodic(
+		const dictionary& 	dict,
+		const plane&    	bplane,
+		internalPoints& 	internal);
+
+	virtual 
+	~boundaryPeriodic() override= default;
+	
+	add_vCtor
+	(
+		boundaryBase,
+		boundaryPeriodic,
+		dictionary
+	);
+
+	real neighborLength()const override;
+
+	realx3 boundaryExtensionLength()const override;
+
+	const plane& boundaryPlane()const override;
+
+	bool beforeIteration(uint32 iterNum, real t, real dt) override;
+
+	bool iterate(uint32 iterNum, real t) override;
+
+	bool afterIteration(uint32 iterNum, real t) override;
+
+
+};
+
+}
+
+#endif

src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp

@@ -8,8 +8,8 @@ bool pFlow::regularSimulationDomain::createBoundaryDicts()
 {
     auto& boundaries = this->subDict("boundaries");
     
-    this->addDict("regularSimulationDomainBoundaries", boundaries);
+    this->addDict("regularBoundaries", boundaries);
-    auto& rbBoundaries = this->subDict("regularSimulationDomainBoundaries");
+    auto& rbBoundaries = this->subDict("regularBoundaries");
     
 	real neighborLength = boundaries.getVal<real>("neighborLength");
 
@@ -37,7 +37,25 @@ bool pFlow::regularSimulationDomain::createBoundaryDicts()
 			" in dictionary "<< boundaries.globalName()<<endl;
 			return false;
 		}
+		auto bType = bDict.getVal<word>("type");
+		uint32 mirrorIndex = 0;
+		if(bType == "periodic")
+		{
+			if(bName == bundaryName(0)) mirrorIndex = 1;
+			if(bName == bundaryName(1)) mirrorIndex = 0;
+			if(bName == bundaryName(2)) mirrorIndex = 3;
+			if(bName == bundaryName(3)) mirrorIndex = 2;
+			if(bName == bundaryName(4)) mirrorIndex = 5;
+			if(bName == bundaryName(5)) mirrorIndex = 4;
 
+			if(! bDict.addOrReplace("mirrorBoundaryIndex", mirrorIndex))
+			{
+				fatalErrorInFunction<<
+				"canno add entry mirroBoundaryIndex to dictionary "<<
+				bDict.globalName()<<endl;
+				return false;
+			}
+		}
 	}
     
     
@@ -168,7 +186,7 @@ bool pFlow::regularSimulationDomain::initialTransferBlockData
 
 const pFlow::dictionary &pFlow::regularSimulationDomain::thisBoundaryDict() const
 {
-    return this->subDict("regularSimulationDomainBoundaries");
+    return this->subDict("regularBoundaries");
 }
 
 bool pFlow::regularSimulationDomain::requiresDataTransfer()const

src/phasicFlow/structuredData/domain/simulationDomain.cpp

@@ -41,8 +41,19 @@ pFlow::simulationDomain::simulationDomain(systemControl& control)
     
 }
 
-pFlow::uniquePtr<pFlow::simulationDomain> 
+pFlow::domain pFlow::simulationDomain::extendThisDomain
-    pFlow::simulationDomain::create(systemControl& control)
+(
+	const realx3 &lowerPointExtension, 
+	const realx3 &upperPointExtension
+) const
+{
+    realx3 minP = thisDomain().domainBox().minPoint() + lowerPointExtension;
+	realx3 maxP = thisDomain().domainBox().maxPoint() + upperPointExtension;
+	return domain({minP, maxP});
+}
+
+pFlow::uniquePtr<pFlow::simulationDomain>
+pFlow::simulationDomain::create(systemControl &control)
 {
 	word sType = angleBracketsNames(
         "simulationDomain", 

src/phasicFlow/structuredData/domain/simulationDomain.hpp

@@ -162,6 +162,11 @@ public:
 		return thisDomain_;
 	}
 
+	
+	domain extendThisDomain(
+		const realx3& lowerPointExtension, 
+		const realx3& upperPointExtension)const;
+
 	inline 
 	const auto& globalBoundaryDict()const
 	{

src/phasicFlow/structuredData/pointStructure/internalPoints.cpp

@@ -22,6 +22,7 @@ Licence:
 #include "internalPoints.hpp"
 #include "domain.hpp"
 #include "Vectors.hpp"
+#include "internalPointsKernels.hpp"
 
 void pFlow::internalPoints::syncPFlag()const
 {
@@ -32,6 +33,54 @@ void pFlow::internalPoints::syncPFlag()const
 	}
 }
 
+bool pFlow::internalPoints::deletePoints(deviceViewType1D<uint32> delPoints)
+{
+    if(!pFlagsD_.deletePoints(delPoints))
+	{
+		fatalErrorInFunction<<
+		"Error in deleting points from internal points"<<endl;
+		return false;
+	}
+	pFlagSync_ = false;
+	WARNING<<"Notify the observersin in internalPoints"<<END_WARNING;
+	return true;
+}
+
+bool pFlow::internalPoints::changePointsFlag
+(
+	deviceViewType1D<uint32> changePoints, 
+	uint32 boundaryIndex
+)
+{
+	if(boundaryIndex>5)
+	{
+		fatalErrorInFunction<<
+		"Invalid boundary index "<< boundaryIndex<<endl;
+		return false;
+	}
+	
+	pFlagsD_.changeFlags(changePoints, boundaryIndex);
+    pFlagSync_ = false;
+	WARNING<<"NOTIFY About transfering the data "<<END_WARNING;
+	return true;
+}
+
+bool pFlow::internalPoints::changePointsPoisition
+(
+	deviceViewType1D<uint32> changePoints, 
+	realx3 transferVector
+)
+{
+	pFlow::internalPointsKernels::changePosition
+	(
+		pointPosition_.deviceViewAll(),
+		changePoints,
+		transferVector
+	);
+	WARNING<<"Notify about the change in the position of points"<<END_WARNING;
+    return true;
+}
+
 pFlow::internalPoints::internalPoints()
 :
 	subscriber("internalPoints"),

src/phasicFlow/structuredData/pointStructure/internalPoints.hpp

@@ -31,6 +31,8 @@ namespace pFlow
 {
 
 class domain;
+class boundaryBase;
+class Time; 
 
 class internalPoints
 :
@@ -70,6 +72,19 @@ protected:
 
 		void syncPFlag()const;
 
+		friend boundaryBase;
+		
+
+		bool deletePoints(deviceViewType1D<uint32> delPoints);
+
+		bool changePointsFlag(
+			deviceViewType1D<uint32> changePoints, 
+			uint32 boundaryIndex);
+
+		bool changePointsPoisition(
+			deviceViewType1D<uint32> changePoints,
+			realx3 transferVector);
+
 public:
 
 	//friend class dynamicinternalPoints;
@@ -165,6 +180,18 @@ public:
 		{
 			return pFlagsD_.activeRange();
 		}
+		
+		virtual 
+		Time& time() = 0;
+
+		virtual
+		const Time& time() const = 0;
+
+		virtual 
+		boundaryBase& boundary(size_t boundaryIndex ) = 0;
+
+		virtual 
+		const boundaryBase& boundary(size_t boundaryIndex ) const = 0;
 
 		///@brief delete points at indices given in delPoints.
 		/// The default is that delPoints contains sorted indices 

src/phasicFlow/structuredData/pointStructure/internalPointsKernels.cpp

@@ -0,0 +1,44 @@
+/*------------------------------- 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 "internalPointsKernels.hpp"
+
+void pFlow::internalPointsKernels::changePosition
+(
+    deviceViewType1D<realx3> points,
+    deviceViewType1D<uint32> indices, 
+    realx3 transferVector
+)
+{
+    auto s = static_cast<uint32>(indices.size());
+    if(s==0u)return;
+
+    Kokkos::parallel_for
+    (
+        "internalPointsKernels::changePosition",
+        deviceRPolicyStatic(0u, s),
+        LAMBDA_HD(uint32 i)
+        {   
+            points[indices[i]] += transferVector;
+        }
+    );
+    Kokkos::fence();
+    
+}

src/phasicFlow/structuredData/pointStructure/internalPointsKernels.hpp

@@ -0,0 +1,31 @@
+/*------------------------------- 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 "phasicFlowKokkos.hpp"
+
+namespace pFlow::internalPointsKernels
+{
+
+void changePosition(
+    deviceViewType1D<realx3> points,
+    deviceViewType1D<uint32> indices,
+    realx3 transferVector);
+
+}

src/phasicFlow/structuredData/pointStructure/pointFlag.hpp

@@ -51,6 +51,10 @@ class pointFlag
 
 	using device_type 		= typename viewType::device_type;
 
+	using rPolicy = Kokkos::RangePolicy<
+		execution_space,  
+		Kokkos::IndexType<uint32>>;
+
 protected:
 	
 	viewType 	flags_;
@@ -75,7 +79,34 @@ protected:
 
 	//- Protected methods 
 
-	
+	uint8 getBoundaryFlag(uint32 index)const
+	{
+		uint8 flg;
+		switch (index){
+			case 0u:
+				flg = Flag::LEFT;
+				break;
+			case 1u:
+				flg = Flag::RIGHT;
+				break;
+			case 2u:
+				flg = Flag::BOTTOM;
+				break;
+			case 3u:
+				flg = Flag::TOP;
+				break;
+			case 4u:
+				flg = Flag::REAR;
+				break;
+			case 5u:
+				flg = Flag::FRONT;
+				break;
+			default:
+				flg=0;
+		}
+
+		return flg;
+	}
 
 public:
 
@@ -255,9 +286,15 @@ public:
 		const box& validBox,
 		ViewType1D<realx3, memory_space> points);
 
-	
 	bool deletePoints(
 		scatteredFieldAccess<uint32, memory_space> points);
+
+	bool deletePoints(
+		ViewType1D<uint32, memory_space> points);
+	
+	bool changeFlags(
+		ViewType1D<uint32, memory_space> changePoints, 
+		uint32 boundaryIndex);
 	
 	/// @brief mark points based on their position in the domain.
 	/// This should be the first method to be called when updating

src/phasicFlow/structuredData/pointStructure/pointFlagKernels.hpp

@@ -158,12 +158,59 @@ bool pFlow::pointFlag<ExecutionSpace>::deletePoints
 {
 	if(points.empty())return true;
 	
-	uint32 minIndex = points.getFirstCopy();
+	//uint32 minIndex = points.getFirstCopy();
-	uint32 maxIndex = points.getLastCopy();
+	//uint32 maxIndex = points.getLastCopy();
 	
-	if(maxIndex<minIndex) return false;
+	//if(maxIndex<minIndex) return false;
-	if(maxIndex>activeRange_.end())return false;
+	//if(maxIndex>activeRange_.end())return false;
-	if(minIndex<activeRange_.start())return false;
+	//if(minIndex<activeRange_.start())return false;
+
+	using policy = Kokkos::RangePolicy<
+		execution_space,  
+		Kokkos::IndexType<uint32>>;
+
+	uint32 numDeleted = 0;
+	Kokkos::parallel_reduce
+	(
+		"pointFlagKernels::deletePoints",
+		policy(0u, points.size()),
+		CLASS_LAMBDA_HD(uint32 i, uint32& valDelUpdate)
+		{
+			uint32 n = points(i);
+			if(isActive(n))
+			{
+				valDelUpdate++;
+				flags_[n] = Flag::DELETED;
+			}
+		},
+		numDeleted
+	);	
+
+	if(numDeleted >= numActive_)
+	{
+		activeRange_ = {0, 0};
+		numDeleted == numActive_;
+	}
+		
+	numActive_ 	= numActive_ - numDeleted;
+	isAllActive_ = 
+		(activeRange_.numElements() == numActive_)&& numActive_>0;
+	
+	return true;
+}
+template<typename ExecutionSpace>
+bool pFlow::pointFlag<ExecutionSpace>::deletePoints
+(
+	ViewType1D<uint32, memory_space> points
+)
+{
+
+	uint32 s = points.size();
+	if(s==0u)return true;
+	
+	//if(maxIndex<minIndex) return false;
+	//if(maxIndex>activeRange_.end())return false;
+	//if(minIndex<activeRange_.start())return false;
 
 	using policy = Kokkos::RangePolicy<
 		execution_space,  
@@ -200,6 +247,27 @@ bool pFlow::pointFlag<ExecutionSpace>::deletePoints
 }
 
 
+template<typename ExecutionSpace>
+bool pFlow::pointFlag<ExecutionSpace>::changeFlags
+(
+	ViewType1D<uint32, memory_space> changePoints, 
+	uint32 boundaryIndex
+)
+{
+	auto flg = getBoundaryFlag(boundaryIndex);
+	Kokkos::parallel_for
+	(
+		"pointFlag::changeFlags",
+		rPolicy(0, changePoints.size()),
+		LAMBDA_HD(uint32 i)
+		{
+			flags_[changePoints(i)] = flg;
+		}
+	);
+	Kokkos::fence();
+	return true;
+}
+
 template<typename ExecutionSpace>
 pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markPointRegions
 (

src/phasicFlow/structuredData/pointStructure/pointStructure.hpp

@@ -119,12 +119,23 @@ public:
             return boundaries_;
         }
 
-        auto& boundary(size_t i )
+		
+		Time& time() override
+		{
+			return demComponent::time();
+		}
+
+		const Time& time() const override
+		{
+			return demComponent::time();
+		}
+
+        boundaryBase& boundary(size_t i ) override
         {
             return boundaries_[i];
         }
 
-        auto& boundary(size_t i)const
+        const boundaryBase& boundary(size_t i)const override
         {
             return boundaries_[i];
         }
@@ -134,6 +145,7 @@ public:
 			return simulationDomain_();
 		}
 
+		
 	// - IO methods 
 
 		/// @brief read the point structure from the input stream.