This is the first merge from main into MPI branch

Merge branch 'main' into local-MPI

src/CMakeLists.txt

@@ -13,3 +13,5 @@ add_subdirectory(Interaction)
 
 add_subdirectory(MotionModel)
 
+add_subdirectory(PostprocessData)
+

src/Geometry/geometryMotion/geometryMotion.cpp

@@ -21,7 +21,6 @@ Licence:
 template<typename MotionModel>
 bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
 {
-	
 	if(motionComponentName().size() != numSurfaces() )
 	{
 		fatalErrorInFunction<<

src/Geometry/geometryMotion/geometryMotions.cpp

@@ -26,4 +26,6 @@ template class pFlow::geometryMotion<pFlow::rotatingAxisMotion>;
 
 template class pFlow::geometryMotion<pFlow::stationaryWall>;
 
-//template class pFlow::geometryMotion<pFlow::multiRotatingAxisMotion>;
+template class pFlow::geometryMotion<pFlow::conveyorBeltMotion>;
+
+template class pFlow::geometryMotion<pFlow::multiRotatingAxisMotion>;

src/Geometry/geometryMotion/geometryMotions.hpp

@@ -24,7 +24,8 @@ Licence:
 #include "geometryMotion.hpp"
 #include "stationaryWall.hpp"
 #include "rotatingAxisMotion.hpp"
-//#include "multiRotatingAxisMotion.hpp"
+#include "conveyorBeltMotion.hpp"
+#include "multiRotatingAxisMotion.hpp"
 #include "vibratingMotion.hpp"
 
 
@@ -37,10 +38,9 @@ using rotationAxisMotionGeometry = geometryMotion<rotatingAxisMotion>;
 
 using stationaryGeometry = geometryMotion<stationaryWall>;
 
-//typedef geometryMotion<multiRotatingAxisMotion> multiRotationAxisMotionGeometry;
-
-
+using conveyorBeltMotionGeometry = geometryMotion<conveyorBeltMotion>;
 
+using multiRotationAxisMotionGeometry = geometryMotion<multiRotatingAxisMotion>;
 
 }
 

src/Integration/AdamsBashforth2/AdamsBashforth2.cpp

@@ -37,7 +37,8 @@ bool intAllActive(
 	real dt, 
 	realx3Field_D& y, 
 	realx3PointField_D& dy,
-	realx3PointField_D& dy1)
+	realx3PointField_D& dy1,
+	real damping = 1.0)
 {
 
 	auto d_dy = dy.deviceView();
@@ -49,7 +50,7 @@ bool intAllActive(
 		"AdamsBashforth2::correct",
 		rpIntegration (activeRng.start(), activeRng.end()),
 		LAMBDA_HD(uint32 i){
-			d_y[i] +=  dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
+			d_y[i] +=  damping * dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
 			d_dy1[i] = d_dy[i];
 		});
 	Kokkos::fence();
@@ -62,7 +63,8 @@ bool intScattered
 	real dt, 
 	realx3Field_D& y,
 	realx3PointField_D& dy,
-	realx3PointField_D& dy1
+	realx3PointField_D& dy1,
+	real damping = 1.0
 )
 {
 
@@ -78,7 +80,7 @@ bool intScattered
 		LAMBDA_HD(uint32 i){
 			if( activeP(i))
 			{
-				d_y[i] +=  dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
+				d_y[i] +=  damping * dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
 				d_dy1[i] = d_dy[i];
 			}
 		});
@@ -90,17 +92,16 @@ bool intScattered
 
 }
 
-
-
 pFlow::AdamsBashforth2::AdamsBashforth2
 (
 	const word& baseName,
 	pointStructure& pStruct,
 	const word& method,
-	const realx3Field_D& initialValField
+	const realx3Field_D& initialValField,
+	bool  keepHistory
 )
 :
-	integration(baseName, pStruct, method, initialValField),
+	integration(baseName, pStruct, method, initialValField, keepHistory),
 	realx3PointField_D
 	(
 		objectFile
@@ -108,20 +109,27 @@ pFlow::AdamsBashforth2::AdamsBashforth2
 			groupNames(baseName,"dy1"),
 			pStruct.time().integrationFolder(),
 			objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS
+			keepHistory?objectFile::WRITE_ALWAYS:objectFile::WRITE_NEVER
 		),
 		pStruct,
 		zero3,
 		zero3
-	)
-{}
+	),
+	initialValField_(initialValField),
+	boundaryList_(pStruct, method, *this)
+{
+	realx3PointField_D::addEvent(message::ITEMS_INSERT);
+}
 
-bool pFlow::AdamsBashforth2::predict
-(
-	real UNUSED(dt),
-	realx3PointField_D& UNUSED(y),
-	realx3PointField_D& UNUSED(dy)
-)
+void pFlow::AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested()
+{
+	realx3PointField_D::updateBoundariesSlaveToMasterIfRequested();
+}
+
+bool pFlow::AdamsBashforth2::predict(
+    real UNUSED(dt),
+    realx3PointField_D &UNUSED(y),
+    realx3PointField_D &UNUSED(dy))
 {
 	return true;
 }
@@ -140,31 +148,63 @@ bool pFlow::AdamsBashforth2::correct
 (
 	real dt,
 	realx3PointField_D& y,
-	realx3PointField_D& dy
+	realx3PointField_D& dy,
+	real damping
 )
-{
-	return correct(dt, y.field(), dy);
-}
-
-bool pFlow::AdamsBashforth2::correct(real dt, realx3Field_D &y, realx3PointField_D &dy)
 {
 	auto& dy1l = dy1();
-
+	bool success = false;
 	if(dy1l.isAllActive())
 	{
-		return intAllActive(dt, y, dy, dy1l);
+		success = intAllActive(dt, y.field(), dy, dy1(), damping);
 	}
 	else
 	{
-		return intScattered(dt, y, dy, dy1l);
+		success = intScattered(dt, y.field(), dy, dy1(), damping);
 	}
-    return false;
+
+	success = success && boundaryList_.correct(dt, y, dy);
+
+	return success;
+
 }
 
-bool pFlow::AdamsBashforth2::setInitialVals(
-	const int32IndexContainer& newIndices,
-	const realx3Vector& y)
+bool pFlow::AdamsBashforth2::correctPStruct(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel)
 {
-	return true;
+	auto& dy1l = dy1();
+	bool success = false;
+	if(dy1l.isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel,  dy1());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel,  dy1());
+	}
+
+	success = success && boundaryList_.correctPStruct(dt, pStruct, vel);
+
+	return success;
 }
 
+/*bool pFlow::AdamsBashforth2::hearChanges
+(
+	const timeInfo &ti, 
+	const message &msg, 
+	const anyList &varList
+)
+{
+	if(msg.equivalentTo(message::ITEMS_INSERT))
+	{
+		
+		return insertValues(varList, initialValField_.deviceViewAll(), dy1());
+	}
+	else
+	{
+		return realx3PointField_D::hearChanges(ti, msg, varList);
+	}
+    
+}*/

src/Integration/AdamsBashforth2/AdamsBashforth2.hpp

@@ -17,13 +17,13 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
 #ifndef __AdamsBashforth2_hpp__
 #define __AdamsBashforth2_hpp__
 
 
 #include "integration.hpp"
 #include "pointFields.hpp"
+#include "boundaryIntegrationList.hpp"
 
 namespace pFlow
 {
@@ -41,10 +41,33 @@ class AdamsBashforth2
 {
 private:
 
+	const realx3Field_D& initialValField_;
+
+	boundaryIntegrationList boundaryList_;
+
+	friend class processorAB2BoundaryIntegration;
+
+protected:
+	
+	const auto& dy1()const
+	{
+		return static_cast<const realx3PointField_D&>(*this);
+	}
+
 	auto& dy1()
 	{
 		return static_cast<realx3PointField_D&>(*this);
 	}
+
+	auto& initialValField()
+	{
+		return initialValField_;
+	}
+
+	boundaryIntegrationList& boundaryList()
+	{
+		return boundaryList_;
+	}
 	
 public:
 
@@ -58,10 +81,11 @@ public:
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 
 		/// Destructor 
-		~AdamsBashforth2()final = default;
+		~AdamsBashforth2()override = default;
 
 		/// Add this to the virtual constructor table 
 		add_vCtor(
@@ -71,6 +95,9 @@ public:
 
 
 	// - Methods
+
+		void updateBoundariesSlaveToMasterIfRequested()override;
+
 		/// return integration method 
 		word method()const override
 		{
@@ -90,31 +117,21 @@ public:
 		bool correct(
 			real dt, 
 			realx3PointField_D& y, 
-			realx3PointField_D& dy) final;
+			realx3PointField_D& dy,
+			real damping = 1.0) override;
 
-		bool correct(
+		bool correctPStruct(
 			real dt, 
-			realx3Field_D& y, 
-			realx3PointField_D& dy) final;
+			pointStructure& pStruct, 
+			realx3PointField_D& vel) override;
 		
 		/*bool hearChanges
 		(
-			real t,
-			real dt,
-			uint32 iter,
+			const timeInfo& ti,
 			const message& msg, 
 			const anyList& varList
 		) override;*/
 
-		bool setInitialVals(
-			const int32IndexContainer& newIndices,
-			const realx3Vector& y) final;
-
-		bool needSetInitialVals()const final
-		{
-			return false;
-		}
-
 };
 
 

src/Integration/AdamsBashforth3/AdamsBashforth3.cpp

@@ -20,90 +20,185 @@ Licence:
 
 #include "AdamsBashforth3.hpp"
 
+namespace pFlow
+{
+
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
+
+bool intAllActive(
+	real dt, 
+	realx3Field_D& y, 
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	real damping = 1.0)
+{
+
+	auto d_dy = dy.deviceView();
+	auto d_y  = y.deviceView();
+	auto d_dy1= dy1.deviceView();
+	auto d_dy2= dy2.deviceView();
+	auto activeRng = dy1.activeRange();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth3::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			d_y[i] += damping* dt*( static_cast<real>(23.0 / 12.0) * d_dy[i] 
+						- static_cast<real>(16.0 / 12.0) * d_dy1[i] 
+						+ static_cast<real>(5.0 / 12.0) * d_dy2[i]) ;
+			
+			d_dy2[i] = d_dy1[i];
+			d_dy1[i] = d_dy[i];
+
+		});
+	Kokkos::fence();
+
+	return true;	
+}
+
+bool intScattered
+(
+	real dt, 
+	realx3Field_D& y,
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	real damping = 1.0
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto d_dy2		= dy2.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth3::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] += damping * dt*( static_cast<real>(23.0 / 12.0) * d_dy[i] 
+						- static_cast<real>(16.0 / 12.0) * d_dy1[i] 
+						+ static_cast<real>(5.0 / 12.0) * d_dy2[i]);
+			
+				d_dy2[i] = d_dy1[i];
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+}
+
 //const real AB3_coef[] = { 23.0 / 12.0, 16.0 / 12.0, 5.0 / 12.0 };
 
 pFlow::AdamsBashforth3::AdamsBashforth3
 (
 	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
+	pointStructure& pStruct,
+	const word& method,
+	const realx3Field_D& initialValField,
+	bool  keepHistory
 )
 :
-	integration(baseName, owner, pStruct, method),
-	history_(
-		owner.emplaceObject<pointField<VectorSingle,AB3History>>(
-			objectFile(
-				groupNames(baseName,"AB3History"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			AB3History({zero3,zero3})))
+	AdamsBashforth2(baseName, pStruct, method, initialValField, keepHistory),
+	dy2_
+	(
+		objectFile
+		(
+			groupNames(baseName,"dy2"),
+			pStruct.time().integrationFolder(),
+			objectFile::READ_IF_PRESENT,
+			keepHistory ? objectFile::WRITE_ALWAYS : objectFile::WRITE_NEVER
+		),
+		pStruct,
+		zero3,
+		zero3
+	)
 
 {
 	
 }
 
-bool pFlow::AdamsBashforth3::predict
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
+void pFlow::AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested()
 {
-
-	return true;
+	AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested();
+	dy2_.updateBoundariesSlaveToMasterIfRequested();
 }
 
+
+
 bool pFlow::AdamsBashforth3::correct
 (
-	real dt,
-	realx3Vector_D& y,
-	realx3Vector_D& dy
+	real dt, 
+	realx3PointField_D& y, 
+	realx3PointField_D& dy,
+	real damping
 )
 {
 	
-	if(this->pStruct().allActive())
+	bool success = false;
+	if(y.isAllActive())
 	{
-		return intAll(dt, y, dy, this->pStruct().activeRange());
+		success = intAllActive(dt, y.field(), dy, dy1(), dy2(), damping);
 	}
 	else
 	{
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
+		success = intScattered(dt, y.field(), dy, dy1(), dy2(), damping);
 	}
 
-	return true;
+	success = success && boundaryList().correct(dt, y, dy);
+
+	return success;
 }
 
-bool pFlow::AdamsBashforth3::setInitialVals(
-	const int32IndexContainer& newIndices,
-	const realx3Vector& y)
-{
-	return true;
-}
-
-bool pFlow::AdamsBashforth3::intAll(
+bool pFlow::AdamsBashforth3::correctPStruct(
 	real dt, 
-	realx3Vector_D& y, 
-	realx3Vector_D& dy, 
-	range activeRng)
+	pointStructure &pStruct, 
+	realx3PointField_D &vel)
 {
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
+	
+	bool success = false;
+	if(dy2().isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2());
+	}
 
-	Kokkos::parallel_for(
-		"AdamsBashforth3::correct",
-		rpIntegration (activeRng.first, activeRng.second),
-		LAMBDA_HD(int32 i){
-			auto ldy = d_dy[i];
-			d_y[i] += dt*( static_cast<real>(23.0 / 12.0) * ldy 
-				- static_cast<real>(16.0 / 12.0) * d_history[i].dy1_ 
-				+ static_cast<real>(5.0 / 12.0) * d_history[i].dy2_);
-			d_history[i] = {ldy ,d_history[i].dy1_};
-		});
-	Kokkos::fence();
+	success = success && boundaryList().correctPStruct(dt, pStruct, vel);
 
-	return true;	
-}
+	return success;
+}
+
+/*bool pFlow::AdamsBashforth3::hearChanges
+(
+	const timeInfo &ti, 
+	const message &msg, 
+	const anyList &varList
+)
+{
+	if(msg.equivalentTo(message::ITEMS_INSERT))
+	{
+		return insertValues(varList, initialValField().deviceViewAll(), dy1()) &&
+		insertValues(varList, initialValField().deviceViewAll(), dy2());
+	}
+	else
+	{
+		return realx3PointField_D::hearChanges(ti, msg, varList);
+	}
+}*/

src/Integration/AdamsBashforth3/AdamsBashforth3.hpp

@@ -22,48 +22,14 @@ Licence:
 #define __AdamsBashforth3_hpp__
 
 
-#include "integration.hpp"
+#include "AdamsBashforth2.hpp"
 #include "pointFields.hpp"
+#include "boundaryIntegrationList.hpp"
+
 
 namespace pFlow
 {
 
-struct AB3History
-{
-	TypeInfoNV("AB3History");
-
-	realx3 dy1_={0,0,0};
-	realx3 dy2_={0,0,0};
-};
-
-
-INLINE_FUNCTION
-iIstream& operator>>(iIstream& str, AB3History& ab3)
-{
-	str.readBegin("AB3History");
-
-    str >> ab3.dy1_;
-    str >> ab3.dy2_;
-    
-    str.readEnd("AB3History");
-
-    str.check(FUNCTION_NAME);
-
-	return str;
-	
-}
-
-INLINE_FUNCTION
-iOstream& operator<<(iOstream& str, const AB3History& ab3)
-{
-	str << token::BEGIN_LIST << ab3.dy1_
-	    << token::SPACE << ab3.dy2_
-		<< token::END_LIST;	    
-
-	str.check(FUNCTION_NAME);
-
-	return str;
-}
 
 /**
  * Third order Adams-Bashforth integration method for solving ODE
@@ -72,19 +38,26 @@ iOstream& operator<<(iOstream& str, const AB3History& ab3)
  */
 class AdamsBashforth3
 :
-	public integration
+	public AdamsBashforth2
 {
-protected:
+private:
 	
-	/// Integration history
-	pointField<VectorSingle,AB3History>& history_;
+	friend class processorAB3BoundaryIntegration;
 
-	/// Range policy for integration kernel
-	using rpIntegration = Kokkos::RangePolicy<
-			DefaultExecutionSpace,
-			Kokkos::Schedule<Kokkos::Static>,
-			Kokkos::IndexType<int32>
-			>;
+	realx3PointField_D 		dy2_;
+
+protected:
+
+	const auto& dy2()const
+	{
+		return dy2_;
+	}
+
+	auto& dy2()
+	{
+		return dy2_;
+	}
+	
 
 public:
 
@@ -96,17 +69,14 @@ public:
 		/// Construct from components
 		AdamsBashforth3(
 			const word& baseName,
-			repository& owner,
-			const pointStructure& pStruct,
-			const word& method);
+			pointStructure& pStruct,
+			const word& method,
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 		
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth3>(*this);
-		}
 
 		/// Destructor 
-		virtual ~AdamsBashforth3()=default;
+		~AdamsBashforth3() override =default;
 
 		/// Add this to the virtual constructor table 
 		add_vCtor(
@@ -117,43 +87,39 @@ public:
 
 	// - Methods
 
-		bool predict(
-			real UNUSED(dt),
-			realx3Vector_D & UNUSED(y),
-			realx3Vector_D& UNUSED(dy)) override;
+		void updateBoundariesSlaveToMasterIfRequested()override;
 
-		bool correct(real dt, 
-			realx3Vector_D & y,
-			realx3Vector_D& dy) override;
-
-		bool setInitialVals(
-			const int32IndexContainer& newIndices,
-			const realx3Vector& y) override;
-
-		bool needSetInitialVals()const override
+		/// return integration method 
+		word method()const override
 		{
-			return false;
+			return "AdamsBashforth3";
 		}
 
-		/// Integrate on all points in the active range 
-		bool intAll(
-			real dt, 
-			realx3Vector_D& y, 
-			realx3Vector_D& dy, 
-			range activeRng);
+		
 
-		/// Integrate on active points in the active range 
-		template<typename activeFunctor>
-		bool intRange(
-			real dt,
-			realx3Vector_D& y,
-			realx3Vector_D& dy,
-			activeFunctor activeP );
+		bool correct(
+			real dt, 
+			realx3PointField_D& y, 
+			realx3PointField_D& dy,
+			real damping = 1.0) override;
+
+		bool correctPStruct(
+			real dt, 
+			pointStructure& pStruct, 
+			realx3PointField_D& vel) override;
+			
+		
+		/*bool hearChanges
+		(
+			const timeInfo& ti,
+			const message& msg, 
+			const anyList& varList
+		) override;*/
 
 };
 
 
-template<typename activeFunctor>
+/*template<typename activeFunctor>
 bool pFlow::AdamsBashforth3::intRange(
 	real dt, 
 	realx3Vector_D& y,
@@ -181,7 +147,7 @@ bool pFlow::AdamsBashforth3::intRange(
 	Kokkos::fence();
 
 	return true;
-}
+}*/
 
 } // pFlow
 

src/Integration/AdamsBashforth4/AdamsBashforth4.cpp

@@ -20,96 +20,187 @@ Licence:
 
 #include "AdamsBashforth4.hpp"
 
-
-
-pFlow::AdamsBashforth4::AdamsBashforth4
-(
-	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
-)
-:
-	integration(baseName, owner, pStruct, method),
-	history_(
-		owner.emplaceObject<pointField<VectorSingle,AB4History>>(
-			objectFile(
-				groupNames(baseName,"AB4History"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			AB4History({zero3,zero3, zero3})))
-
+namespace pFlow
 {
 
-}
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
 
-bool pFlow::AdamsBashforth4::predict
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
-{
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth4::correct
-(
-	real dt,
-	realx3Vector_D& y,
-	realx3Vector_D& dy
-)
-{
-	
-	if(this->pStruct().allActive())
-	{
-		return intAll(dt, y, dy, this->pStruct().activeRange());
-	}
-	else
-	{
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
-	}
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth4::setInitialVals(
-	const int32IndexContainer& newIndices,
-	const realx3Vector& y)
-{
-	return true;
-}
-
-bool pFlow::AdamsBashforth4::intAll(
+bool intAllActive(
 	real dt, 
-	realx3Vector_D& y, 
-	realx3Vector_D& dy, 
-	range activeRng)
+	realx3Field_D& y, 
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	real damping = 1.0)
 {
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
+
+	auto d_dy = dy.deviceView();
+	auto d_y  = y.deviceView();
+	auto d_dy1= dy1.deviceView();
+	auto d_dy2= dy2.deviceView();
+	auto d_dy3= dy3.deviceView();
+	auto activeRng = dy1.activeRange();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth4::correct",
-		rpIntegration (activeRng.first, activeRng.second),
-		LAMBDA_HD(int32 i){				
-			d_y[i] += dt*( 
-						static_cast<real>(55.0 / 24.0) * d_dy[i]
-					- static_cast<real>(59.0 / 24.0) * d_history[i].dy1_ 
-					+ static_cast<real>(37.0 / 24.0) * d_history[i].dy2_
-					- static_cast<real>( 9.0 / 24.0) * d_history[i].dy3_
-					);
-			d_history[i].dy3_ = d_history[i].dy2_;
-			d_history[i].dy2_ = d_history[i].dy1_;
-			d_history[i].dy1_ = d_dy[i];
-			
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			d_y[i] += damping * dt*( 
+				  static_cast<real>(55.0 / 24.0) * d_dy[i]
+				- static_cast<real>(59.0 / 24.0) * d_dy1[i] 
+				+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
+				- static_cast<real>( 9.0 / 24.0) * d_dy3[i]);
 
+			d_dy3[i] = d_dy2[i];
+			d_dy2[i] = d_dy1[i];
+			d_dy1[i] = d_dy[i];
 		});
 	Kokkos::fence();
 
 	return true;	
-}
+}
+
+bool intScattered
+(
+	real dt, 
+	realx3Field_D& y,
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	real damping = 1.0
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto d_dy2		= dy2.deviceView();
+	auto d_dy3		= dy3.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth4::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] += damping*  dt*( 
+				  static_cast<real>(55.0 / 24.0) * d_dy[i]
+				- static_cast<real>(59.0 / 24.0) * d_dy1[i] 
+				+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
+				- static_cast<real>( 9.0 / 24.0) * d_dy3[i] );
+				d_dy3[i] = d_dy2[i];
+				d_dy2[i] = d_dy1[i];
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+}
+
+pFlow::AdamsBashforth4::AdamsBashforth4
+(
+	const word& baseName,
+	pointStructure& pStruct,
+	const word& method,
+	const realx3Field_D& initialValField,
+	bool keepHistory
+)
+:
+	AdamsBashforth3(baseName, pStruct, method, initialValField, keepHistory),
+	dy3_
+	(
+		objectFile
+		(
+			groupNames(baseName,"dy3"),
+			pStruct.time().integrationFolder(),
+			objectFile::READ_IF_PRESENT,
+			keepHistory?objectFile::WRITE_ALWAYS:objectFile::WRITE_NEVER
+		),
+		pStruct,
+		zero3,
+		zero3
+	)
+{}
+
+void pFlow::AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested()
+{
+	AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested();
+	dy3_.updateBoundariesSlaveToMasterIfRequested();
+}
+
+bool pFlow::AdamsBashforth4::correct
+(
+	real dt, 
+	realx3PointField_D& y, 
+	realx3PointField_D& dy,
+	real damping
+)
+{
+	
+	bool success = false;
+	if(y.isAllActive())
+	{
+		success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3(), damping);
+	}
+	else
+	{
+		success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3(), damping);
+	}
+
+	success = success && boundaryList().correct(dt, y, dy);
+
+	return success;
+}
+
+bool pFlow::AdamsBashforth4::correctPStruct(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel)
+{
+	
+	bool success = false;
+	if(dy2().isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
+	}
+
+	success = success && boundaryList().correctPStruct(dt, pStruct, vel);
+
+	return success;
+}
+
+/*bool pFlow::AdamsBashforth4::hearChanges
+(
+	const timeInfo &ti, 
+	const message &msg, 
+	const anyList &varList
+)
+{
+	if(msg.equivalentTo(message::ITEMS_INSERT))
+	{
+		return insertValues(varList, initialValField().deviceViewAll(), dy1()) &&
+		insertValues(varList, initialValField().deviceViewAll(), dy2()) && 
+		insertValues(varList, initialValField().deviceViewAll(), dy3());
+	}
+	else
+	{
+		return realx3PointField_D::hearChanges(ti, msg, varList);
+	}
+}*/

src/Integration/AdamsBashforth4/AdamsBashforth4.hpp

@@ -22,53 +22,14 @@ Licence:
 #define __AdamsBashforth4_hpp__
 
 
-#include "integration.hpp"
+#include "AdamsBashforth3.hpp"
 #include "pointFields.hpp"
+#include "boundaryIntegrationList.hpp"
+
 
 namespace pFlow
 {
 
-struct AB4History
-{
-	TypeInfoNV("AB4History");
-
-	realx3 dy1_={0,0,0};
-	realx3 dy2_={0,0,0};
-	realx3 dy3_={0,0,0};
-	
-};
-
-
-INLINE_FUNCTION
-iIstream& operator>>(iIstream& str, AB4History& ab4)
-{
-	str.readBegin("AB4History");
-
-    str >> ab4.dy1_;
-    str >> ab4.dy2_;
-    str >> ab4.dy3_;
-    
-    str.readEnd("AB4History");
-
-    str.check(FUNCTION_NAME);
-
-	return str;
-	
-}
-
-INLINE_FUNCTION
-iOstream& operator<<(iOstream& str, const AB4History& ab4)
-{
-	str << token::BEGIN_LIST << ab4.dy1_
-	    << token::SPACE << ab4.dy2_
-	    << token::SPACE << ab4.dy3_
-		<< token::END_LIST;	    
-
-	str.check(FUNCTION_NAME);
-
-	return str;
-}
-
 /**
  * Fourth order Adams-Bashforth integration method for solving ODE
  * 
@@ -76,19 +37,25 @@ iOstream& operator<<(iOstream& str, const AB4History& ab4)
  */
 class AdamsBashforth4
 :
-	public integration
+	public AdamsBashforth3
 {
+private:
+
+	friend class processorAB4BoundaryIntegration;
+
+	realx3PointField_D 		dy3_;
+
 protected:
 
-	/// Integration history
-	pointField<VectorSingle,AB4History>& history_;
+	const auto& dy3()const
+	{
+		return dy3_;
+	}
 
-	/// Range policy for integration kernel
-	using rpIntegration = Kokkos::RangePolicy<
-			DefaultExecutionSpace,
-			Kokkos::Schedule<Kokkos::Static>,
-			Kokkos::IndexType<int32>
-			>;
+	auto& dy3()
+	{
+		return dy3_;
+	}
 
 public:
 
@@ -100,17 +67,15 @@ public:
 		/// Construct from components
 		AdamsBashforth4(
 			const word& baseName,
-			repository& owner,
-			const pointStructure& pStruct,
-			const word& method);
+			pointStructure& pStruct,
+			const word& method,
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
+		
 		
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth4>(*this);
-		}
 
 		/// Destructor
-		virtual ~AdamsBashforth4()=default;
+		~AdamsBashforth4() override =default;
 
 		/// Add a this to the virtual constructor table 
 		add_vCtor(
@@ -121,77 +86,36 @@ public:
 
 	// - Methods
 		
-		bool predict(
-			real UNUSED(dt),
-			realx3Vector_D & UNUSED(y),
-			realx3Vector_D& UNUSED(dy)) override;
+		void updateBoundariesSlaveToMasterIfRequested()override;
+
+		/// return integration method 
+		word method()const override
+		{
+			return "AdamsBashforth4";
+		}
 
 		bool correct(
 			real dt, 
-			realx3Vector_D & y,
-			realx3Vector_D& dy) override;
+			realx3PointField_D& y, 
+			realx3PointField_D& dy,
+			real damping = 1.0) override;
 
-		bool setInitialVals(
-			const int32IndexContainer& newIndices,
-			const realx3Vector& y) override;
-
-		bool needSetInitialVals()const override
-		{
-			return false;
-		}
-
-		/// Integrate on all points in the active range 
-		bool intAll(
+		bool correctPStruct(
 			real dt, 
-			realx3Vector_D& y, 
-			realx3Vector_D& dy, 
-			range activeRng);
-
-		/// Integrate on active points in the active range 
-		template<typename activeFunctor>
-		bool intRange(
-			real dt,
-			realx3Vector_D& y,
-			realx3Vector_D& dy,
-			activeFunctor activeP );
+			pointStructure& pStruct, 
+			realx3PointField_D& vel) override;
+		
+		/*bool hearChanges
+		(
+			const timeInfo& ti,
+			const message& msg, 
+			const anyList& varList
+		) override;*/
 
 };
 
 
-template<typename activeFunctor>
-bool pFlow::AdamsBashforth4::intRange(
-	real dt, 
-	realx3Vector_D& y,
-	realx3Vector_D& dy,
-	activeFunctor activeP )
-{
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
-	auto activeRng = activeP.activeRange();
 
-	Kokkos::parallel_for(
-		"AdamsBashforth4::correct",
-		rpIntegration (activeRng.first, activeRng.second),
-		LAMBDA_HD(int32 i){
-			if( activeP(i))
-			{
-				
-				d_y[i] += dt*( 
-							static_cast<real>(55.0 / 24.0) * d_dy[i]
-						- static_cast<real>(59.0 / 24.0) * d_history[i].dy1_ 
-						+ static_cast<real>(37.0 / 24.0) * d_history[i].dy2_
-						- static_cast<real>( 9.0 / 24.0) * d_history[i].dy3_
-						);
-				d_history[i].dy3_ = d_history[i].dy2_;
-				d_history[i].dy2_ = d_history[i].dy1_;
-				d_history[i].dy1_ = d_dy[i];
-			}
-		});
-	Kokkos::fence();
-
-	return true;
-}
 
 } // pFlow
 

src/Integration/AdamsBashforth5/AdamsBashforth5.cpp

@@ -20,93 +20,179 @@ Licence:
 
 #include "AdamsBashforth5.hpp"
 
-
-
-pFlow::AdamsBashforth5::AdamsBashforth5
-(
-	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
-)
-:
-	integration(baseName, owner, pStruct, method),
-	history_(
-		owner.emplaceObject<pointField<VectorSingle,AB5History>>(
-			objectFile(
-				groupNames(baseName,"AB5History"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			AB5History({zero3,zero3, zero3})))
-
+namespace pFlow
 {
 
-}
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
 
-bool pFlow::AdamsBashforth5::predict
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
-{
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth5::correct
-(
-	real dt,
-	realx3Vector_D& y,
-	realx3Vector_D& dy
-)
-{
-	
-	if(this->pStruct().allActive())
-	{
-		return intAll(dt, y, dy, this->pStruct().activeRange());
-	}
-	else
-	{
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
-	}
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth5::setInitialVals(
-	const int32IndexContainer& newIndices,
-	const realx3Vector& y)
-{
-	return true;
-}
-
-bool pFlow::AdamsBashforth5::intAll(
+bool intAllActive(
 	real dt, 
-	realx3Vector_D& y, 
-	realx3Vector_D& dy, 
-	range activeRng)
+	realx3Field_D& y, 
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	realx3PointField_D& dy4,
+	real damping = 1.0)
 {
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
+
+	auto d_dy = dy.deviceView();
+	auto d_y  = y.deviceView();
+	auto d_dy1= dy1.deviceView();
+	auto d_dy2= dy2.deviceView();
+	auto d_dy3= dy3.deviceView();
+	auto d_dy4= dy4.deviceView();
+	auto activeRng = dy1.activeRange();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth5::correct",
-		rpIntegration (activeRng.first, activeRng.second),
-		LAMBDA_HD(int32 i){				
-			d_y[i] += dt*( 
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			d_y[i] += damping*  dt*( 
 				  static_cast<real>(1901.0 / 720.0) * d_dy[i]
-				- static_cast<real>(2774.0 / 720.0) * d_history[i].dy1_ 
-				+ static_cast<real>(2616.0 / 720.0) * d_history[i].dy2_
-				- static_cast<real>(1274.0 / 720.0) * d_history[i].dy3_
-				+ static_cast<real>( 251.0 / 720.0) * d_history[i].dy4_
-				);
-			d_history[i] = {d_dy[i] ,d_history[i].dy1_, d_history[i].dy2_, d_history[i].dy3_};
+				- static_cast<real>(2774.0 / 720.0) * d_dy1[i]
+				+ static_cast<real>(2616.0 / 720.0) * d_dy2[i]
+				- static_cast<real>(1274.0 / 720.0) * d_dy3[i]
+				+ static_cast<real>( 251.0 / 720.0) * d_dy4[i]);
+				
+			d_dy4[i] = d_dy3[i];
+			d_dy3[i] = d_dy2[i];
+			d_dy2[i] = d_dy1[i];
+			d_dy1[i] = d_dy[i];
 		});
 	Kokkos::fence();
 
 	return true;	
-}
+}
+
+bool intScattered
+(
+	real dt, 
+	realx3Field_D& y,
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	realx3PointField_D& dy4,
+	real damping = 1.0
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto d_dy2		= dy2.deviceView();
+	auto d_dy3		= dy3.deviceView();
+	auto d_dy4		= dy4.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth5::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] += damping*  dt*( 
+				  static_cast<real>(1901.0 / 720.0) * d_dy[i]
+				- static_cast<real>(2774.0 / 720.0) * d_dy1[i]
+				+ static_cast<real>(2616.0 / 720.0) * d_dy2[i]
+				- static_cast<real>(1274.0 / 720.0) * d_dy3[i]
+				+ static_cast<real>( 251.0 / 720.0) * d_dy4[i]);
+
+				d_dy4[i] = d_dy3[i];
+				d_dy3[i] = d_dy2[i];
+				d_dy2[i] = d_dy1[i];
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+	return true;
+}
+
+}
+
+pFlow::AdamsBashforth5::AdamsBashforth5
+(
+	const word &baseName, 
+	pointStructure &pStruct, 
+	const word &method, 
+	const realx3Field_D &initialValField,
+	bool  keepHistory
+)
+:
+	AdamsBashforth4(baseName, pStruct, method, initialValField, keepHistory),
+	dy4_
+	(
+		objectFile
+		(
+			groupNames(baseName,"dy4"),
+			pStruct.time().integrationFolder(),
+			objectFile::READ_IF_PRESENT,
+			keepHistory?objectFile::WRITE_ALWAYS:objectFile::WRITE_NEVER
+		),
+		pStruct,
+		zero3,
+		zero3
+	)
+{}
+
+void pFlow::AdamsBashforth5::updateBoundariesSlaveToMasterIfRequested()
+{
+	AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested();
+	dy4_.updateBoundariesSlaveToMasterIfRequested();
+}
+
+bool pFlow::AdamsBashforth5::correct
+(
+	real dt, 
+	realx3PointField_D &y, 
+	realx3PointField_D &dy, 
+	real damping
+)
+{
+	
+	bool success = false;
+	if(y.isAllActive())
+	{
+		success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3(), dy4(), damping);
+	}
+	else
+	{
+		success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3(), dy4(), damping);
+	}
+
+	success = success && boundaryList().correct(dt, y, dy);
+
+	return success;
+}
+
+bool pFlow::AdamsBashforth5::correctPStruct
+(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel
+)
+{
+bool success = false;
+	if(dy2().isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3(), dy4());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3(), dy4());
+	}
+
+	success = success && boundaryList().correctPStruct(dt, pStruct, vel);
+
+	return success;
+}
+
+

src/Integration/AdamsBashforth5/AdamsBashforth5.hpp

@@ -22,54 +22,12 @@ Licence:
 #define __AdamsBashforth5_hpp__
 
 
-#include "integration.hpp"
+#include "AdamsBashforth4.hpp"
 #include "pointFields.hpp"
 
 namespace pFlow
 {
 
-struct AB5History
-{
-	TypeInfoNV("AB5History");
-
-	realx3 dy1_={0,0,0};
-	realx3 dy2_={0,0,0};
-	realx3 dy3_={0,0,0};
-	realx3 dy4_={0,0,0};
-};
-
-
-INLINE_FUNCTION
-iIstream& operator>>(iIstream& str, AB5History& ab5)
-{
-	str.readBegin("AB5History");
-
-    str >> ab5.dy1_;
-    str >> ab5.dy2_;
-    str >> ab5.dy3_;
-    str >> ab5.dy4_;
-    
-    str.readEnd("AB5History");
-
-    str.check(FUNCTION_NAME);
-
-	return str;
-	
-}
-
-INLINE_FUNCTION
-iOstream& operator<<(iOstream& str, const AB5History& ab5)
-{
-	str << token::BEGIN_LIST << ab5.dy1_
-	    << token::SPACE << ab5.dy2_
-	    << token::SPACE << ab5.dy3_
-	    << token::SPACE << ab5.dy4_
-		<< token::END_LIST;	    
-
-	str.check(FUNCTION_NAME);
-
-	return str;
-}
 
 /**
  * Fifth order Adams-Bashforth integration method for solving ODE
@@ -78,19 +36,25 @@ iOstream& operator<<(iOstream& str, const AB5History& ab5)
  */
 class AdamsBashforth5
 :
-	public integration
+	public AdamsBashforth4
 {
+private:
+
+	friend class processorAB4BoundaryIntegration;
+
+	realx3PointField_D 		dy4_;
+
 protected:
 
-	/// Integration history
-	pointField<VectorSingle,AB5History>& history_;
+	const auto& dy4()const
+	{
+		return dy4_;
+	}
 
-	/// Range policy for integration kernel
-	using rpIntegration = Kokkos::RangePolicy<
-			DefaultExecutionSpace,
-			Kokkos::Schedule<Kokkos::Static>,
-			Kokkos::IndexType<int32>
-			>;
+	auto& dy4()
+	{
+		return dy4_;
+	}
 
 public:
 
@@ -99,20 +63,20 @@ public:
 
 	// - Constructors
 		
+		/// Construct from components
 		AdamsBashforth5(
 			const word& baseName,
-			repository& owner,
-			const pointStructure& pStruct,
-			const word& method);
+			pointStructure& pStruct,
+			const word& method,
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
+		
+		
 
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth5>(*this);
-		}
+		/// Destructor
+		~AdamsBashforth5() override =default;
 
-		virtual ~AdamsBashforth5()=default;
-
-		/// Add this to the virtual constructor table 
+		/// Add a this to the virtual constructor table 
 		add_vCtor(
 			integration,
 			AdamsBashforth5,
@@ -121,44 +85,29 @@ public:
 
 	// - Methods
 		
-		bool predict(
-			real UNUSED(dt),
-			realx3Vector_D & UNUSED(y),
-			realx3Vector_D& UNUSED(dy)) override;
+		void updateBoundariesSlaveToMasterIfRequested()override;
+
+		/// return integration method 
+		word method()const override
+		{
+			return "AdamsBashforth5";
+		}
 
 		bool correct(
 			real dt, 
-			realx3Vector_D & y,
-			realx3Vector_D& dy) override;
+			realx3PointField_D& y, 
+			realx3PointField_D& dy,
+			real damping = 1.0) override;
 
-		bool setInitialVals(
-			const int32IndexContainer& newIndices,
-			const realx3Vector& y) override;
-
-		bool needSetInitialVals()const override
-		{
-			return false;
-		}
-
-		/// Integrate on all points in the active range 
-		bool intAll(
+		bool correctPStruct(
 			real dt, 
-			realx3Vector_D& y, 
-			realx3Vector_D& dy, 
-			range activeRng);
-
-		/// Integrate on active points in the active range 
-		template<typename activeFunctor>
-		bool intRange(
-			real dt,
-			realx3Vector_D& y,
-			realx3Vector_D& dy,
-			activeFunctor activeP );
+			pointStructure& pStruct, 
+			realx3PointField_D& vel) override;
 
 };
 
 
-template<typename activeFunctor>
+/*template<typename activeFunctor>
 bool pFlow::AdamsBashforth5::intRange(
 	real dt, 
 	realx3Vector_D& y,
@@ -190,7 +139,7 @@ bool pFlow::AdamsBashforth5::intRange(
 	Kokkos::fence();
 
 	return true;
-}
+}*/
 
 } // pFlow
 

src/Integration/CMakeLists.txt

@@ -1,15 +1,23 @@
 
 list(APPEND SourceFiles 
-integration/integration.cpp 
+integration/integration.cpp
+boundaries/boundaryIntegration.cpp
+boundaries/boundaryIntegrationList.cpp 
 AdamsBashforth2/AdamsBashforth2.cpp
-#AdamsBashforth5/AdamsBashforth5.cpp
-#AdamsBashforth4/AdamsBashforth4.cpp
-#AdamsBashforth3/AdamsBashforth3.cpp
+AdamsBashforth3/AdamsBashforth3.cpp
+AdamsBashforth4/AdamsBashforth4.cpp
+AdamsBashforth5/AdamsBashforth5.cpp
+
 #AdamsMoulton3/AdamsMoulton3.cpp
 #AdamsMoulton4/AdamsMoulton4.cpp
 #AdamsMoulton5/AdamsMoulton5.cpp
 )
 
+if(pFlow_Build_MPI)
+ list(APPEND SourceFiles 
+ AdamsBashforth2/processorAB2BoundaryIntegration.cpp)
+endif()
+
 set(link_libs Kokkos::kokkos phasicFlow)
 
 pFlow_add_library_install(Integration SourceFiles link_libs)

src/Integration/boundaries/boundaryIntegration.hpp

@@ -51,9 +51,7 @@ public:
     );
 
     bool hearChanges(
-        real t,
-        real dt,
-        uint32 iter,
+        const timeInfo& ti,
         const message &msg,
         const anyList &varList) override
     {
@@ -77,6 +75,11 @@ public:
         return true;
     }
 
+    bool isActive()const override
+    {
+        return false;
+    }
+    
     static
     uniquePtr<boundaryIntegration> create(
         const boundaryBase& boundary,

src/Integration/boundaries/boundaryIntegrationList.cpp

@@ -7,11 +7,11 @@ pFlow::boundaryIntegrationList::boundaryIntegrationList(
     integration &intgrtn
 )
 :
-    ListPtr<boundaryIntegration>(6),
+    boundaryListPtr<boundaryIntegration>(),
 	boundaries_(pStruct.boundaries())
 {
-    
-    for(uint32 i=0; i<6; i++)
+	
+	ForAllBoundariesPtr(i, this)
 	{
 		this->set(
 			i,
@@ -19,23 +19,40 @@ pFlow::boundaryIntegrationList::boundaryIntegrationList(
 				boundaries_[i],
 				pStruct, 
 				method,
-                intgrtn
-            )
-        );
+				intgrtn
+			)
+		);
 	}	
-
 }
 
 bool pFlow::boundaryIntegrationList::correct(real dt, realx3PointField_D &y, realx3PointField_D &dy)
 {
-	for(auto& bndry:*this)
+	ForAllActiveBoundariesPtr(i,this)
 	{
-		if(!bndry->correct(dt, y, dy))
+		if(!boundaryPtr(i)->correct(dt, y, dy))
 		{
-			fatalErrorInFunction;
+			fatalErrorInFunction<<"Error in correcting boundary "<<
+			boundaryPtr(i)->boundaryName()<<endl;
 			return false;
 		}
-		
-    }
+	}
+    return true;
+}
+
+bool pFlow::boundaryIntegrationList::correctPStruct(
+	real dt, 
+	pointStructure &pStruct, 
+	const realx3PointField_D &vel)
+{
+	ForAllActiveBoundariesPtr(i,this)
+	{
+		if(!boundaryPtr(i)->correctPStruct(dt, vel))
+		{
+			fatalErrorInFunction<<"Error in correcting boundary "<<
+			boundaryPtr(i)->boundaryName()<<" in pointStructure."<<endl;
+			return false;
+		}
+	}
+	
     return true;
 }

src/Integration/boundaries/boundaryIntegrationList.hpp

@@ -4,7 +4,7 @@
 
 
 #include "boundaryList.hpp"
-#include "ListPtr.hpp"
+#include "boundaryListPtr.hpp"
 #include "boundaryIntegration.hpp"
 
 
@@ -15,7 +15,7 @@ class integration;
 
 class boundaryIntegrationList
 :
-    public ListPtr<boundaryIntegration>
+    public boundaryListPtr<boundaryIntegration>
 {
 private:
 
@@ -35,6 +35,12 @@ public:
 		real dt, 
 		realx3PointField_D& y, 
 		realx3PointField_D& dy);
+
+    bool correctPStruct(
+        real dt, 
+        pointStructure& pStruct, 
+        const realx3PointField_D& vel);
+
 };
 
 

src/Integration/integration/integration.cpp

@@ -22,17 +22,41 @@ Licence:
 #include "pointStructure.hpp"
 #include "repository.hpp"
 
-pFlow::integration::integration
+bool pFlow::integration::insertValues
 (
-	const word& baseName,
-	pointStructure& pStruct,
-	const word&,
-	const realx3Field_D&
+	const anyList& varList, 
+	const deviceViewType1D<realx3>& srcVals,
+	realx3PointField_D& dstFeild
 )
-:
-	owner_(*pStruct.owner()),
-	pStruct_(pStruct),
-	baseName_(baseName)	
+{
+	const word eventName = message::eventName(message::ITEMS_INSERT);
+
+	if(!varList.checkObjectType<uint32IndexContainer>(eventName))
+	{
+		fatalErrorInFunction<<"Insertion failed in integration for "<< baseName_<<
+		", variable "<< eventName <<
+		" does not exist or the type is incorrect"<<endl;
+		return false;
+	}
+
+	const auto& indices = varList.getObject<uint32IndexContainer>(
+		eventName);
+	
+	dstFeild.field().insertSetElement(indices, srcVals);
+
+    return true;
+}
+
+pFlow::integration::integration(
+    const word &baseName,
+    pointStructure &pStruct,
+    const word &,
+    const realx3Field_D &,
+	bool keepHistory)
+    : owner_(*pStruct.owner()),
+      pStruct_(pStruct),
+      baseName_(baseName),
+      keepHistory_(keepHistory)
 {}
 
 
@@ -42,12 +66,13 @@ pFlow::uniquePtr<pFlow::integration>
 		const word& baseName,
 		pointStructure& pStruct,
 		const word& method,
-		const realx3Field_D& initialValField
+		const realx3Field_D& initialValField,
+		bool  keepHistory
 )
 {
 	if( wordvCtorSelector_.search(method) )
 	{
-		return wordvCtorSelector_[method] (baseName, pStruct, method, initialValField);
+		return wordvCtorSelector_[method] (baseName, pStruct, method, initialValField, keepHistory);
 	}
 	else
 	{

src/Integration/integration/integration.hpp

@@ -24,6 +24,7 @@ Licence:
 
 #include "virtualConstructor.hpp"
 #include "pointFields.hpp"
+#include "Logical.hpp"
 
 
 namespace pFlow
@@ -63,6 +64,15 @@ private:
 		/// The base name for integration 
 		const word baseName_;
 
+		bool 	keepHistory_;
+
+protected:
+
+	bool insertValues(
+		const anyList& varList, 
+		const deviceViewType1D<realx3>& srcVals,
+		realx3PointField_D& dstFeild);
+
 public:
 
 	/// Type info
@@ -76,7 +86,8 @@ public:
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 
 		/// Copy constructor 
 		integration(const integration&) = default;
@@ -102,9 +113,10 @@ public:
 				const word& baseName,
 				pointStructure& pStruct,
 				const word& method,
-				const realx3Field_D& initialValField
+				const realx3Field_D& initialValField,
+				bool keepHistory
 			),
-			(baseName, pStruct, method, initialValField)
+			(baseName, pStruct, method, initialValField, keepHistory)
 		);
 
 
@@ -131,6 +143,13 @@ public:
 			return owner_;
 		}
 
+		bool keepHistory()const
+		{
+			return keepHistory_;
+		}
+
+		virtual 
+		void updateBoundariesSlaveToMasterIfRequested() = 0;
 		/// return integration method 
 		virtual 
 		word method()const = 0 ;
@@ -144,21 +163,10 @@ public:
 
 		/// Correction/main integration step
 		virtual 
-		bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0;
+		bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy, real damping = 1.0) = 0;
 
 		virtual 
-		bool correct(real dt, realx3Field_D& y, realx3PointField_D& dy) = 0;
-
-		/// Set the initial values for new indices 
-		virtual 
-		bool setInitialVals(
-				const int32IndexContainer& newIndices,
-				const realx3Vector& y) = 0;
-
-		/// Check if the method requires any set initial vals
-		virtual 
-		bool needSetInitialVals()const = 0;
-
+		bool correctPStruct(real dt, pointStructure& pStruct, realx3PointField_D& vel) = 0;
 
 	/// Create the polymorphic object based on inputs
 	static
@@ -166,7 +174,8 @@ public:
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 
 }; // integration
 

src/Interaction/CMakeLists.txt

@@ -21,6 +21,12 @@ interaction/interaction.cpp
 sphereInteraction/sphereInteractionsLinearModels.cpp
 sphereInteraction/sphereInteractionsNonLinearModels.cpp
 sphereInteraction/sphereInteractionsNonLinearModModels.cpp
+
+grainInteraction/grainInteractionsLinearModels.cpp
+grainInteraction/grainInteractionsNonLinearModels.cpp
+grainInteraction/grainInteractionsNonLinearModModels.cpp
+
+
 )
 
 if(pFlow_Build_MPI)

src/Interaction/Models/contactForce/cGAbsoluteLinearCF.hpp

@@ -0,0 +1,339 @@
+/*------------------------------- 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 __cGAbsoluteLinearCF_hpp__
+#define __cGAbsoluteLinearCF_hpp__
+
+#include "types.hpp"
+#include "symArrays.hpp"
+
+
+namespace pFlow::cfModels
+{
+
+template<bool limited=true>
+class cGAbsoluteLinear
+{
+public:
+
+	struct contactForceStorage
+	{
+		realx3 overlap_t_ = 0.0;
+	};
+	
+	
+	struct linearProperties
+	{		
+		real 	kn_ 	= 1000.0;
+		real 	kt_	= 800.0;
+		real 	en_ 	= 0.0;
+		real 	ethat_ 	= 0.0;
+		real    mu_ 	= 0.00001;
+
+		INLINE_FUNCTION_HD
+		linearProperties(){}
+
+		INLINE_FUNCTION_HD
+		linearProperties(real kn, real kt, real en, real etha_t, real mu ):
+			kn_(kn), kt_(kt), en_(en),ethat_(etha_t), mu_(mu) 
+		{}		
+
+		INLINE_FUNCTION_HD
+		linearProperties(const linearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		linearProperties& operator=(const linearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		~linearProperties() = default;
+	};
+	
+	
+
+protected:
+
+	using LinearArrayType = symArray<linearProperties>;
+	
+	int32 					numMaterial_ = 0;
+
+	ViewType1D<real>  		 rho_;
+
+	LinearArrayType  	linearProperties_;
+	
+	int32 		addDissipationModel_;	
+	
+
+	bool readLinearDictionary(const dictionary& dict)
+	{
+	
+		
+		auto kn = dict.getVal<realVector>("kn");
+		auto kt = dict.getVal<realVector>("kt");
+		auto en = dict.getVal<realVector>("en");
+		auto et = dict.getVal<realVector>("et");
+		auto mu = dict.getVal<realVector>("mu");
+
+		auto nElem = kn.size();
+
+		if(nElem != kt.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of kn("<<nElem<<") and kt("<<kt.size()<<") do not match.\n";
+			return false;
+		}
+
+		if(nElem != kt.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of kn("<<nElem<<") and kt("<<kt.size()<<") do not match.\n";
+			return false;
+		}
+
+		if(nElem != en.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of kn("<<nElem<<") and en("<<en.size()<<") do not match.\n";
+			return false;
+		}
+
+		if(nElem != et.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of kn("<<nElem<<") and et("<<et.size()<<") do not match.\n";
+			return false;
+		}
+
+		if(nElem != mu.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of kn("<<nElem<<") and mu("<<mu.size()<<") do not match.\n";
+			return false;
+		}
+
+		
+		// check if size of vector matchs a symetric array
+		uint32 nMat;
+		if( !LinearArrayType::getN(nElem, nMat) )
+		{
+			fatalErrorInFunction<<
+			"sizes of properties do not match a symetric array with size ("<<
+			numMaterial_<<"x"<<numMaterial_<<").\n";
+			return false;
+		}
+		else if( numMaterial_ != nMat)
+		{
+			fatalErrorInFunction<<
+			"size mismatch for porperties. \n"<<
+			"you supplied "<< numMaterial_<<" items in materials list and "<<
+			nMat << " for other properties.\n";
+			return false;
+		}
+
+
+		realVector etha_t("etha_t", nElem);
+
+		ForAll(i , kn)
+		{
+
+			etha_t[i] = -2.0*log( et[i])*sqrt(kt[i]*2/7) /
+					sqrt(log(pow(et[i],2))+ pow(Pi,2));
+		}
+
+		Vector<linearProperties> prop("prop", nElem);
+		ForAll(i,kn)
+		{
+			prop[i] = {kn[i], kt[i], en[i], etha_t[i], mu[i] };
+		}
+
+		linearProperties_.assign(prop);
+		
+		auto adm = dict.getVal<word>("additionalDissipationModel");
+
+		
+		if(adm == "none")
+		{
+			addDissipationModel_ = 1;
+		}
+		else if(adm == "LU")
+		{
+			addDissipationModel_ = 2;
+		}
+		else if (adm == "GB")
+		{
+			addDissipationModel_ = 3;
+		}
+		else
+		{
+			addDissipationModel_ = 1;
+		}
+		
+		return true;
+
+	}
+
+	static const char* modelName()
+	{
+		if constexpr (limited)
+		{
+			return "cGAbsoluteLinearLimited";
+		}
+		else
+		{
+			return "cGAbsoluteLinearNonLimited";
+		}
+		return "";
+	}
+
+public:
+
+
+	TypeInfoNV(modelName());
+
+	INLINE_FUNCTION_HD
+	cGAbsoluteLinear(){}
+
+	cGAbsoluteLinear(int32 nMaterial, const ViewType1D<real>& rho, const dictionary& dict)
+	:
+		numMaterial_(nMaterial),
+		rho_("rho",nMaterial),
+		linearProperties_("linearProperties",nMaterial)
+	{
+
+		Kokkos::deep_copy(rho_,rho);
+		if(!readLinearDictionary(dict))
+		{
+			fatalExit;
+		}
+	}
+	
+	INLINE_FUNCTION_HD
+	cGAbsoluteLinear(const cGAbsoluteLinear&) = default;
+
+	INLINE_FUNCTION_HD
+	cGAbsoluteLinear(cGAbsoluteLinear&&) = default;
+
+	INLINE_FUNCTION_HD
+	cGAbsoluteLinear& operator=(const cGAbsoluteLinear&) = default;
+
+	INLINE_FUNCTION_HD
+	cGAbsoluteLinear& operator=(cGAbsoluteLinear&&) = default;
+
+
+	INLINE_FUNCTION_HD
+	~cGAbsoluteLinear()=default;
+
+	INLINE_FUNCTION_HD
+	int32 numMaterial()const
+	{
+		return numMaterial_;
+	}
+
+	//// - Methods
+	INLINE_FUNCTION_HD
+	void contactForce
+	(
+		const real dt,
+		const uint32 i,
+		const uint32 j,
+		const uint32 propId_i,
+		const uint32 propId_j,
+		const real Ri,
+		const real Rj,
+		const real cGFi,
+		const real cGFj,
+		const real ovrlp_n,
+		const realx3& Vr,
+		const realx3& Nij,
+		contactForceStorage& history,
+		realx3& FCn,
+		realx3& FCt
+	)const
+	{
+
+		
+		auto prop = linearProperties_(propId_i,propId_j);		
+
+
+		real f_ = ( cGFi + cGFj )/2 ;
+
+
+		real vrn = dot(Vr, Nij);	
+		realx3 Vt = Vr - vrn*Nij;
+
+		history.overlap_t_ += Vt*dt;
+
+		real mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i];
+		real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j];
+
+		real sqrt_meff = sqrt((mi*mj)/(mi+mj));
+
+
+		// disipation model
+		if (addDissipationModel_==2)
+		{
+			prop.en_ = sqrt(1+((pow(prop.en_,2)-1)*f_));
+		}
+		else if (addDissipationModel_==3)
+		{
+			auto pie =3.14;
+			prop.en_ = exp((pow(f_,static_cast<real>(1.5))*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) ));
+		}
+
+		real ethan_  = -2.0*log(prop.en_)*sqrt(prop.kn_)/
+					sqrt(pow(log(prop.en_),2)+ pow(Pi,2));
+
+			//REPORT(0)<<"\n en n is  :  "<<END_REPORT;
+			//REPORT(0)<< prop.en_ <<END_REPORT;
+
+		FCn = (  -pow(f_,3)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,static_cast<real>(1.5)) * ethan_ * vrn)*Nij;
+		FCt = ( -pow(f_,3)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,static_cast<real>(1.5)) *  prop.ethat_*Vt);
+				
+		real ft = length(FCt);
+		real ft_fric = prop.mu_ * length(FCn);
+		
+		if(ft > ft_fric)
+		{
+			if( length(history.overlap_t_) >zero)
+			{
+				if constexpr (limited)
+				{
+					FCt *= (ft_fric/ft);
+                	history.overlap_t_ = - (FCt/prop.kt_);
+				}
+				else
+				{
+					FCt = (FCt/ft)*ft_fric;	
+				}
+				//cout<<"friction is applied here \n";
+				
+			}
+			else
+			{
+				FCt = 0.0;
+			}
+		}
+
+	}
+	
+};
+
+} //pFlow::cfModels
+
+#endif

src/Interaction/Models/contactForce/cGNonLinearCF.hpp

@@ -0,0 +1,307 @@
+/*------------------------------- 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 __cGNonLinearCF_hpp__
+#define __cGNonLinearCF_hpp__
+
+#include "types.hpp"
+
+namespace pFlow::cfModels
+{
+
+template<bool limited=true>
+class cGNonLinear
+{
+public:
+
+	struct contactForceStorage
+	{
+		realx3 overlap_t_ = 0.0;
+	};
+
+	struct cGNonLinearProperties
+	{		
+		real 	Yeff_	= 1000000.0;
+		real 	Geff_	= 8000000.0;
+		real 	en_ 	= 1.0;
+		real  	mu_ 	= 0.00001;
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties(){}
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties(real Yeff, real Geff, real en, real mu ):
+			Yeff_(Yeff), Geff_(Geff), en_(en), mu_(mu)
+		{}		
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties(const cGNonLinearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties& operator=(const cGNonLinearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		~cGNonLinearProperties() = default;
+	};
+
+protected:
+
+	using cGNonLinearArrayType = 	symArray<cGNonLinearProperties>;
+
+	int32 						numMaterial_ = 0;
+
+	ViewType1D<real>  		 	rho_;
+
+	int32 						addDissipationModel_;
+
+	cGNonLinearArrayType  		nonlinearProperties_;
+
+	bool readNonLinearDictionary(const dictionary& dict)
+	{
+		auto Yeff = dict.getVal<realVector>("Yeff");
+		auto Geff = dict.getVal<realVector>("Geff");
+		auto nu = dict.getVal<realVector>("nu");
+		auto en = dict.getVal<realVector>("en");
+		auto mu = dict.getVal<realVector>("mu");
+
+		auto nElem = Yeff.size();
+
+		if(nElem != nu.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and nu("<<nu.size()<<") do not match.\n";
+			return false;
+		}
+
+		if(nElem != en.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and en("<<en.size()<<") do not match.\n";
+			return false;
+		}
+
+
+		if(nElem != mu.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and mu("<<mu.size()<<") do not match.\n";
+			return false;
+		}
+
+		
+		// check if size of vector matchs a symetric array
+		uint32 nMat;
+		if( !cGNonLinearArrayType::getN(nElem, nMat) )
+		{
+			fatalErrorInFunction<<
+			"sizes of properties do not match a symetric array with size ("<<
+			numMaterial_<<"x"<<numMaterial_<<").\n";
+			return false;
+		}
+		else if( numMaterial_ != nMat)
+		{
+			fatalErrorInFunction<<
+			"size mismatch for porperties. \n";
+			return false;
+		}
+
+		Vector<cGNonLinearProperties> prop("prop",nElem);
+		ForAll(i,Yeff)
+		{
+			prop[i] = {Yeff[i], Geff[i], en[i], mu[i]};
+		}
+
+		nonlinearProperties_.assign(prop);
+
+		auto adm = dict.getVal<word>("additionalDissipationModel");
+		if(adm == "none")
+		{
+			addDissipationModel_ = 1;
+		}
+		else if(adm == "LU")
+		{
+			addDissipationModel_ = 2;
+		}
+		else if (adm == "GB")
+		{
+			addDissipationModel_ = 3;
+		}
+		else
+		{
+			addDissipationModel_ = 1;
+		}
+
+		return true;
+
+	}
+
+	static const char* modelName()
+	{
+		if constexpr (limited)
+		{
+			return "cGNonLinearLimited";
+		}
+		else
+		{
+			return "cGNonLinearNonLimited";
+		}
+		return "";
+	}
+
+public:
+
+	TypeInfoNV(modelName());
+
+	INLINE_FUNCTION_HD
+	cGNonLinear(){}
+
+	cGNonLinear(
+		int32 nMaterial,
+		const ViewType1D<real>& rho,
+		const dictionary& dict)
+	:
+		numMaterial_(nMaterial),
+		rho_("rho",nMaterial),
+		nonlinearProperties_("cGNonLinearProperties",nMaterial)
+	{
+
+		Kokkos::deep_copy(rho_,rho);
+		if(!readNonLinearDictionary(dict))
+		{
+			fatalExit;
+		}
+	}
+	
+	INLINE_FUNCTION_HD
+	cGNonLinear(const cGNonLinear&) = default;
+
+	INLINE_FUNCTION_HD
+	cGNonLinear(cGNonLinear&&) = default;
+
+	INLINE_FUNCTION_HD
+	cGNonLinear& operator=(const cGNonLinear&) = default;
+
+	INLINE_FUNCTION_HD
+	cGNonLinear& operator=(cGNonLinear&&) = default;
+
+
+	INLINE_FUNCTION_HD
+	~cGNonLinear()=default;
+
+	INLINE_FUNCTION_HD
+	int32 numMaterial()const
+	{
+		return numMaterial_;
+	}
+
+	//// - Methods
+
+	INLINE_FUNCTION_HD
+	void contactForce
+	(
+		const real dt,
+		const uint32 i,
+		const uint32 j,
+		const uint32 propId_i,
+		const uint32 propId_j,
+		const real Ri,
+		const real Rj,
+		const real cGFi,
+		const real cGFj,
+		const real ovrlp_n,
+		const realx3& Vr,
+		const realx3& Nij,
+		contactForceStorage& history,
+		realx3& FCn,
+		realx3& FCt
+	)const
+	{
+
+		const auto prop = nonlinearProperties_(propId_i,propId_j);
+
+		const real f = 2/( 1/cGFi + 1/cGFj );
+
+		real vrn = dot(Vr, Nij);	
+		realx3 Vt = Vr - vrn*Nij;
+
+		history.overlap_t_ += Vt*dt;
+
+		real mi = 3*Pi/4*pow(Ri,static_cast<real>(3))*rho_[propId_i];
+		real mj = 3*Pi/4*pow(Rj,static_cast<real>(3))*rho_[propId_j];
+		real Reff = 1/(1/Ri + 1/Rj);
+
+		real K_hertz = 4.0/3.0*prop.Yeff_*sqrt(Reff);
+		real sqrt_meff_K_hertz = sqrt((mi*mj)/(mi+mj) * K_hertz);
+
+		real en = prop.en_;
+		if (addDissipationModel_==2)
+		{
+			en = sqrt(1+((pow(prop.en_,2)-1)*f));
+		}
+		else if (addDissipationModel_==3)
+		{
+			
+			en = exp((pow(f,static_cast<real>(1.5))*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2))))/(1-(pow(f,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2)))) ) ));
+		}
+		
+		real Kn = static_cast<real>(4.0/3.0) * prop.Yeff_ * sqrt(Reff*ovrlp_n);
+
+		real ethan_  = -2.0*log(en)*sqrt(Kn)/
+			sqrt(pow(log(en),2)+ pow(Pi,2));
+
+		FCn = ( - Kn*ovrlp_n - 
+			   sqrt_meff_K_hertz*ethan_*pow(ovrlp_n,static_cast<real>(0.25))*vrn)*Nij;
+
+		FCt = (f*f)*(- static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n) ) * history.overlap_t_;
+
+		real ft = length(FCt);
+		real ft_fric = prop.mu_ * length(FCn);
+
+		// apply friction 
+		if(ft > ft_fric)
+		{
+			if( length(history.overlap_t_) >0.0)
+			{
+				if constexpr (limited)
+				{
+					real kt = static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n);
+					FCt *= (ft_fric/ft);
+                	history.overlap_t_ = - FCt/(f*f*kt);
+				}
+				else
+				{
+					FCt = (FCt/ft)*ft_fric;	
+				}
+				
+			}
+			else
+			{
+				FCt = 0.0;
+			}
+		}
+
+	}
+   
+
+};
+
+} //pFlow::CFModels
+
+#endif

src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp

@@ -0,0 +1,326 @@
+/*------------------------------- 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 __cGRelativeLinearCF_hpp__
+#define __cGRelativeLinearCF_hpp__
+
+#include "types.hpp"
+#include "symArrays.hpp"
+
+
+
+namespace pFlow::cfModels
+{
+
+template<bool limited=true>
+class cGRelativeLinear
+{
+public:
+
+	struct contactForceStorage
+	{
+		realx3 overlap_t_ = 0.0;
+	};
+	
+	
+	struct linearProperties
+	{		
+		real 	kn_ 	= 1000.0;
+		real 	kt_	= 800.0;
+		real 	en_ 	= 1.0;
+		real  mu_ 	= 0.00001;
+
+		INLINE_FUNCTION_HD
+		linearProperties(){}
+
+		INLINE_FUNCTION_HD
+		linearProperties(real kn, real kt, real en, real mu  ):
+			kn_(kn), kt_(kt), en_(en), mu_(mu)
+		{}		
+
+		INLINE_FUNCTION_HD
+		linearProperties(const linearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		linearProperties& operator=(const linearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		~linearProperties() = default;
+	};
+	
+	
+
+protected:
+
+	using LinearArrayType = symArray<linearProperties>;
+	
+	int32 					numMaterial_ = 0;
+
+	ViewType1D<real>  		 rho_;
+
+	LinearArrayType  	linearProperties_;
+	
+	int32 		addDissipationModel_;	
+	
+
+	bool readLinearDictionary(const dictionary& dict)
+	{
+	
+		
+		auto kn = dict.getVal<realVector>("kn");
+		auto kt = dict.getVal<realVector>("kt");
+		auto en = dict.getVal<realVector>("en");
+		auto mu = dict.getVal<realVector>("mu");
+		
+
+		auto nElem = kn.size();
+
+
+		if(nElem != kt.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of kn("<<nElem<<") and kt("<<kt.size()<<") do not match.\n";
+			return false;
+		}
+
+		if(nElem != en.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of kn("<<nElem<<") and en("<<en.size()<<") do not match.\n";
+			return false;
+		}
+
+
+		if(nElem != mu.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of kn("<<nElem<<") and mu("<<mu.size()<<") do not match.\n";
+			return false;
+		}
+
+		
+		// check if size of vector matchs a symetric array
+		uint32 nMat;
+		if( !LinearArrayType::getN(nElem, nMat) )
+		{
+			fatalErrorInFunction<<
+			"sizes of properties do not match a symetric array with size ("<<
+			numMaterial_<<"x"<<numMaterial_<<").\n";
+			return false;
+		}
+		else if( numMaterial_ != nMat)
+		{
+			fatalErrorInFunction<<
+			"size mismatch for porperties. \n"<<
+			"you supplied "<< numMaterial_<<" items in materials list and "<<
+			nMat << " for other properties.\n";
+			return false;
+		}
+
+
+		Vector<linearProperties> prop("prop", nElem);
+		ForAll(i,kn)
+		{
+			prop[i] = {kn[i], kt[i], en[i], mu[i]  };
+		}
+
+		linearProperties_.assign(prop);
+		
+
+		auto adm = dict.getVal<word>("additionalDissipationModel");
+		if(adm == "none")
+		{
+			addDissipationModel_ = 1;
+		}
+		else if(adm == "LU")
+		{
+			addDissipationModel_ = 2;
+		}
+		else if (adm == "GB")
+		{
+			addDissipationModel_ = 3;
+		}
+		else
+		{
+			addDissipationModel_ = 1;
+		}
+
+		return true;
+
+	}
+
+	static const char* modelName()
+	{
+		if constexpr (limited)
+		{
+			return "cGRelativeLinearLimited";
+		}
+		else
+		{
+			return "cGRelativeLinearNonLimited";
+		}
+		return "";
+	}
+
+public:
+
+
+	TypeInfoNV(modelName());
+
+	INLINE_FUNCTION_HD
+	cGRelativeLinear(){}
+
+	cGRelativeLinear(int32 nMaterial, const ViewType1D<real>& rho, const dictionary& dict)
+	:
+		numMaterial_(nMaterial),
+		rho_("rho",nMaterial),
+		linearProperties_("linearProperties",nMaterial)
+	{
+
+		Kokkos::deep_copy(rho_,rho);
+		if(!readLinearDictionary(dict))
+		{
+			fatalExit;
+		}
+	}
+	
+	INLINE_FUNCTION_HD
+	cGRelativeLinear(const cGRelativeLinear&) = default;
+
+	INLINE_FUNCTION_HD
+	cGRelativeLinear(cGRelativeLinear&&) = default;
+
+	INLINE_FUNCTION_HD
+	cGRelativeLinear& operator=(const cGRelativeLinear&) = default;
+
+	INLINE_FUNCTION_HD
+	cGRelativeLinear& operator=(cGRelativeLinear&&) = default;
+
+
+	INLINE_FUNCTION_HD
+	~cGRelativeLinear()=default;
+
+	INLINE_FUNCTION_HD
+	int32 numMaterial()const
+	{
+		return numMaterial_;
+	}
+
+	//// - Methods
+	INLINE_FUNCTION_HD
+	void contactForce
+	(
+		const real dt,
+		const uint32 i,
+		const uint32 j,
+		const uint32 propId_i,
+		const uint32 propId_j,
+		const real Ri,
+		const real Rj,
+		const real cGFi,
+		const real cGFj,
+		const real ovrlp_n,
+		const realx3& Vr,
+		const realx3& Nij,
+		contactForceStorage& history,
+		realx3& FCn,
+		realx3& FCt
+	)const
+	{
+		
+		auto prop = linearProperties_(propId_i,propId_j);		
+
+		real f_ = ( cGFi + cGFj )/2 ;
+
+
+		real vrn = dot(Vr, Nij);	
+		realx3 Vt = Vr - vrn*Nij;
+
+		history.overlap_t_ += Vt*dt;
+
+		real mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i];
+		real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j];
+
+		real sqrt_meff = sqrt((mi*mj)/(mi+mj));
+		
+		real en = prop.en_;
+		if (addDissipationModel_==2)
+		{
+			en = sqrt(1+((pow(prop.en_,2)-1)*f_));
+		}
+		else if (addDissipationModel_==3)
+		{
+			
+			en = exp(
+				(
+					pow(f_,static_cast<real>(1.5))*
+					log(prop.en_)*
+					sqrt( 
+						(1-((pow(log(prop.en_),static_cast<real>(2.0))
+						)/
+						(
+							pow(log(prop.en_),static_cast<real>(2.0))+
+							pow(Pi,static_cast<real>(2.0)))))/
+							(1-(pow(f_,3)*(pow(log(prop.en_),2))/
+							(pow(log(prop.en_),static_cast<real>(2.0))+
+							pow(Pi,static_cast<real>(2.0))))) ) ));
+		}
+
+		real ethan_  = -2.0*log(en)*sqrt(prop.kn_)/
+			sqrt(pow(log(en),2)+ pow(Pi,2));
+		
+
+		FCn = ( -f_*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,half) * ethan_ * vrn)*Nij;
+		FCt = ( -f_*prop.kt_ * history.overlap_t_);
+		
+
+		
+		real ft = length(FCt);
+		real ft_fric = prop.mu_ * length(FCn);
+		
+		if(ft > ft_fric)
+		{
+			if( length(history.overlap_t_) >static_cast<real>(0.0))
+			{
+				if constexpr (limited)
+				{
+					FCt *= (ft_fric/ft);
+                	history.overlap_t_ = - (FCt/prop.kt_);
+				}
+				else
+				{
+					FCt = (FCt/ft)*ft_fric;	
+				}
+				//cout<<"friction is applied here \n";
+				
+			}
+			else
+			{
+				FCt = 0.0;
+			}
+		}
+
+	}
+	
+};
+
+} //pFlow::cfModels
+
+#endif

src/Interaction/Models/contactForce/linearCF.hpp

@@ -139,10 +139,10 @@ protected:
 		ForAll(i , kn)
 		{
 			etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/
-					sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0));
+					sqrt(pow(log(en[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
 
 			etha_t[i] = -2.0*log( et[i]*sqrt(kt[i]) )/
-					sqrt(pow(log(et[i]),2.0)+ pow(Pi,2.0));
+					sqrt(pow(log(et[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
 		}
 
 		Vector<linearProperties> prop("prop", nElem);
@@ -243,8 +243,8 @@ public:
 
 		history.overlap_t_ += Vt*dt;
 
-		real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i];
-		real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j];
+		real mi = 3*Pi/4*pow(Ri,static_cast<real>(3.0))*rho_[propId_i];
+		real mj = 3*Pi/4*pow(Rj,static_cast<real>(3.0))*rho_[propId_j];
 
 		real sqrt_meff = sqrt((mi*mj)/(mi+mj));
 

src/Interaction/Models/contactForce/nonLinearCF.hpp

@@ -131,7 +131,7 @@ protected:
 			// we take out sqrt(meff*K_hertz) here and then consider this term 
 			// when calculating damping part. 
 			etha_n[i] = -2.2664*log(en[i])/
-					sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0));
+					sqrt(pow(log(en[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
 
 			// no damping for tangential part 
 
@@ -255,7 +255,7 @@ public:
 		// apply friction 
 		if(ft > ft_fric)
 		{
-			if( length(history.overlap_t_) >0.0)
+			if( length(history.overlap_t_) >zero)
 			{
 				if constexpr (limited)
 				{

src/Interaction/Models/grainContactForceModels.hpp

@@ -0,0 +1,50 @@
+/*------------------------------- 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 __grainContactForceModels_hpp__
+#define __grainContactForceModels_hpp__
+
+#include "cGAbsoluteLinearCF.hpp"
+#include "cGRelativeLinearCF.hpp"
+#include "cGNonLinearCF.hpp"
+
+#include "grainRolling.hpp"
+
+
+
+namespace pFlow::cfModels 
+{
+
+
+using limitedCGAbsoluteLinearGrainRolling 	= grainRolling<cGAbsoluteLinear<true>>;
+using nonLimitedCGAbsoluteLinearGrainRolling = grainRolling<cGAbsoluteLinear<false>>;
+
+using limitedCGRelativeLinearGrainRolling   = grainRolling<cGRelativeLinear<true>>;
+using nonLimitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<false>>;
+
+using limitedCGNonLinearGrainRolling   = grainRolling<cGNonLinear<true>>;
+using nonLimitedCGNonLinearGrainRolling = grainRolling<cGNonLinear<false>>;
+
+
+}
+
+
+
+#endif //__grainContactForceModels_hpp__

src/Interaction/Models/rolling/grainRolling.hpp

@@ -0,0 +1,119 @@
+/*------------------------------- 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 __grainRolling_hpp__
+#define __grainRolling_hpp__
+
+
+namespace pFlow::cfModels
+{
+
+template<typename contactForceModel>
+class grainRolling
+:
+	public contactForceModel
+{
+public:
+	
+	using contactForceStorage = 
+		typename contactForceModel::contactForceStorage;
+
+
+	realSymArray_D  	mur_;
+
+	bool readGrainDict(const dictionary& dict)
+	{
+		auto mur = dict.getVal<realVector>("mur");
+
+		uint32 nMat;
+
+		if(!realSymArray_D::getN(mur.size(),nMat) || nMat != this->numMaterial())
+		{
+			fatalErrorInFunction<<
+			"wrong number of values supplied in mur. \n";
+			return false;
+		}
+
+		mur_.assign(mur);
+
+		return true; 
+	}
+
+public:
+
+	TypeInfoNV(word("normal<"+contactForceModel::TYPENAME()+">"));
+	
+
+	grainRolling(int32 nMaterial, const ViewType1D<real>& rho, const dictionary& dict)
+	:
+		contactForceModel(nMaterial, rho, dict),
+		mur_("mur", nMaterial)
+	{
+		if(!readGrainDict(dict))
+		{
+			fatalExit;
+		}
+	}
+
+	INLINE_FUNCTION_HD
+	void rollingFriction
+	(
+		const real dt,
+		const uint32 i,
+		const uint32 j,
+		const uint32 propId_i,
+		const uint32 propId_j,
+		const real Ri,
+		const real Rj,
+		const real cGFi,
+		const real cGFj,
+		const realx3& wi,
+		const realx3& wj,
+		const realx3& Nij,
+		const realx3& FCn,
+		realx3& Mri,
+		realx3& Mrj
+	)const
+	{
+
+		realx3 w_hat = wi-wj;
+		real w_hat_mag = length(w_hat);
+
+		if( !equal(w_hat_mag,static_cast<real>(0.0)) )
+			w_hat /= w_hat_mag;
+		else
+			w_hat = static_cast<real>(0.0);
+
+		auto Reff = (Ri*Rj)/(Ri+Rj);
+
+		Mri = ( -mur_(propId_i,propId_j) *length(FCn) * Reff ) * w_hat ;
+       
+       	//removing the normal part 
+       	// Mri = Mri - ( (Mri .dot. nij)*nij )
+
+       	Mrj = -Mri;
+	}
+
+
+};
+
+}
+ 
+#endif 

src/Interaction/Models/rolling/normalRolling.hpp

@@ -94,10 +94,10 @@ public:
 		realx3 w_hat = wi-wj;
 		real w_hat_mag = length(w_hat);
 
-		if( !equal(w_hat_mag,0.0) )
+		if( !equal(w_hat_mag,static_cast<real>(0.0)) )
 			w_hat /= w_hat_mag;
 		else
-			w_hat = 0.0;
+			w_hat = static_cast<real>(0.0);
 
 		auto Reff = (Ri*Rj)/(Ri+Rj);
 

src/Interaction/contactLists/sortedContactList.hpp

@@ -80,13 +80,17 @@ public:
 
 	TypeInfoNV("sortedContactList");
 
-
-	explicit sortedContactList(uint32 initialSize =1)
+	sortedContactList(uint32 initialSize =1)
 	:
-		SortedPairs(initialSize),
-		values_("values", SortedPairs::capacity()),
-		sortedPairs0_("sortedPairs0", SortedPairs::capacity()),
-		values0_("values0", SortedPairs::capacity())
+		sortedContactList("sortedContactList", initialSize)
+	{}
+	
+	sortedContactList(const word& name, uint32 initialSize =1)
+	:
+		SortedPairs(name, initialSize),
+		values_(groupNames(name, "values"), SortedPairs::capacity()),
+		sortedPairs0_(groupNames(name, "sortedPairs0"), SortedPairs::capacity()),
+		values0_(groupNames(name, "values0"), SortedPairs::capacity())
 	{}
 
 	bool beforeBroadSearch()

src/Interaction/contactLists/sortedPairs.hpp

@@ -110,11 +110,11 @@ public:
 
 
 	// constructors
-	explicit sortedPairs(uint32 initialSize =1)
+	explicit sortedPairs(const word& name, uint32 initialSize =1)
 	:
 		UnsortedPairs(initialSize),
-		flags_("flags_",UnsortedPairs::capacity()+1),
-		sortedPairs_("sortedPairs_",UnsortedPairs::capacity())
+		flags_( groupNames(name, "flags_"), UnsortedPairs::capacity()+1),
+		sortedPairs_(groupNames(name, "sortedPairs_"), UnsortedPairs::capacity())
 	{}
 
 	
@@ -193,7 +193,7 @@ public:
 
 		if( capacity+1 > flags_.size() )
 		{
-			reallocNoInit(flags_, capacity+1);
+			reallocInit(flags_, capacity+1);
 		}
 
 		// fill the flags
@@ -219,7 +219,7 @@ public:
 		{
 			// get more space to prevent reallocations in next iterations
 			uint32 len = size_*1.1+1; 
-			reallocNoInit(sortedPairs_, len);
+			reallocInit(sortedPairs_, len);
 		}
 		
 		Kokkos::parallel_for(
@@ -231,6 +231,7 @@ public:
 		// - sort paris based on the first and second
 		sort(sortedPairs_, 0, size_ );
 		
+		
 
 	}
 

src/Interaction/contactLists/unsortedContactList.hpp

@@ -82,11 +82,16 @@ public:
 	TypeInfoNV("unsortedContactList");
 
 	explicit unsortedContactList(uint32 capacity=1)
+	:
+	unsortedContactList("unsortedContactList", capacity)
+	{}
+	
+	unsortedContactList(const word& name, uint32 capacity=1)
 	:
 		UnsortedPairs(capacity),
-		values_("values", UnsortedPairs::capacity()),
+		values_(groupNames(name, "values"), UnsortedPairs::capacity()),
 		container0_(capacity),
-		values0_("values0",container0_.capacity())
+		values0_(groupNames(name, "values0"),container0_.capacity())
 	{}
 
 

src/Interaction/contactLists/unsortedPairs.hpp

@@ -194,7 +194,9 @@ public:
 	{
 		uint newCap = container_.capacity()+len;
 		this->clear();
+		//output<<"----------------before "<<capacity()<< " " << size()<<endl;
 		container_.rehash(newCap);
+		//output<<"----------------after "<<capacity()<< " " << size()<<endl;
 	}
 
 	INLINE_FUNCTION_H

src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp

@@ -53,6 +53,47 @@ private:
 
 	boundaryContactSearchList 		csBoundaries_;
 
+	bool BroadSearch(
+		const timeInfo& ti,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force = false) override
+	{
+		
+		auto position = Particles().pointPosition().deviceViewAll();
+		auto flags = Particles().dynPointStruct().activePointsMaskDevice();
+		auto diam = Particles().boundingSphere().deviceViewAll();
+
+		return ppwContactSearch_().broadSearch(
+			ti.iter(),
+			ti.t(),
+			ti.dt(),
+			ppPairs, 
+			pwPairs,
+			position, 
+			flags, 
+			diam,
+			force
+		);
+	}
+
+	bool BoundaryBroadSearch(
+		uint32 bndryIndex,
+		const timeInfo& ti,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force = false)override
+	{		
+		return csBoundaries_[bndryIndex].broadSearch(
+			ti.iter(), 
+			ti.t(), 
+			ti.dt(), 
+			ppPairs, 
+			pwPairs, 
+			force
+		);
+	}
+
 public:
 
 	TypeInfoTemplate11("ContactSearch", SearchMethodType);
@@ -93,7 +134,7 @@ public:
 		ppwContactSearch_ = 
 			makeUnique<SearchMethodType>
 			(
-				dict(),
+				csDict,
 				this->extendedDomainBox(),
 				minD,
 				maxD,
@@ -108,95 +149,19 @@ public:
 			);
 	}
 
-
 	add_vCtor(
 		contactSearch,
 		ContactSearch,
 		dictionary);
 
+	bool enterBroadSearchBoundary(const timeInfo& ti, bool force=false)const override;
 	
-	bool broadSearch(
-		uint32 iter,
-		real t,
-		real dt,
-		csPairContainerType& ppPairs,
-		csPairContainerType& pwPairs,
-		bool force = false) override
-	{
-		
-
-		ppTimer().start();
-	
-		const auto& position = Particles().pointPosition().deviceViewAll();
-		const auto& flags = Particles().dynPointStruct().activePointsMaskDevice();
-		const auto& diam = Particles().boundingSphere().deviceViewAll();
-		
-		
-		if( !ppwContactSearch_().broadSearch(
-			iter,
-			t,
-			dt,
-			ppPairs, 
-			pwPairs,
-			position, 
-			flags, 
-			diam,
-			force) )
-		{
-			fatalErrorInFunction;
-			return false;
-		}
-		ppTimer().end();
-
-		return true;
-	}
-
-	bool boundaryBroadSearch(
-		uint32 i,
-		uint32 iter,
-		real t,
-		real dt,
-		csPairContainerType& ppPairs,
-		csPairContainerType& pwPairs,
-		bool force = false)override
-	{
-		if(i==0u)
-			Particles().boundingSphere().updateBoundaries(DataDirection::SlaveToMaster);
-		return csBoundaries_[i].broadSearch(
-			iter, 
-			t, 
-			dt, 
-			ppPairs, 
-			pwPairs, 
-			force);
-	}
-
-	bool enterBroadSearch(uint32 iter, real t, real dt)const override
-	{
-		if(ppwContactSearch_)
-		{
-			return ppwContactSearch_().performSearch(iter); 
-		}
-		return false;
-	}
-
-	bool performedBroadSearch()const override
-	{
-		return ppwContactSearch_().performedSearch();
-	}
-
-	
-	uint32 updateInterval()const override
-	{
-		return ppwContactSearch_().updateInterval();
-	}
 
 	real sizeRatio()const override
 	{
 		return ppwContactSearch_().sizeRatio();
 	}
 
-	 
 	real cellExtent()const override
 	{
 		return ppwContactSearch_().cellExtent();
@@ -204,7 +169,12 @@ public:
 	
 };
 
+template <class searchMethod>
+inline bool ContactSearch<searchMethod>::enterBroadSearchBoundary(const timeInfo &ti, bool force) const
+{
+	return performSearch(ti.iter(),force) || csBoundaries_.boundariesUpdated();
 }
 
+}
 
 #endif //__ContactSearch_hpp__

src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.hpp

@@ -73,9 +73,7 @@ public:
     }
 
     bool hearChanges(
-        real t,
-        real dt,
-        uint32 iter,
+        const timeInfo& ti,
         const message &msg,
         const anyList &varList) override
     {
@@ -83,8 +81,10 @@ public:
         if (msg.equivalentTo(message::BNDR_RESET))
         {
             // do nothing
+            return true;
         }
-        return true;
+        fatalErrorInFunction;
+        return false;
     }
 
     virtual bool broadSearch(
@@ -98,6 +98,11 @@ public:
         return true;
     }
 
+    bool isActive()const override
+    {
+        return false;
+    }
+
     static uniquePtr<boundaryContactSearch> create(
         const dictionary &dict,
         const boundaryBase &boundary,

src/Interaction/contactSearch/boundaries/boundaryContactSearchList.cpp

@@ -1,3 +1,23 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions.
+
+  phasicFlow is distributed to help others in their research in the field of
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
 #include "boundaryContactSearchList.hpp"
 #include "boundaryList.hpp"
 
@@ -5,7 +25,7 @@ void pFlow::boundaryContactSearchList::setList(
     const dictionary &dict, 
     const contactSearch &cSearch)
 {
-    for(auto i=0; i<boundaries_.size(); i++)
+    ForAllBoundariesPtr(i, this)
     {
         this->set
         (
@@ -20,10 +40,13 @@ pFlow::boundaryContactSearchList::boundaryContactSearchList(
     const boundaryList& bndrs,
     const contactSearch &cSearch)
     :
-    ListPtr(bndrs.size()), 
+    boundaryListPtr(), 
     boundaries_(bndrs)
 {
     setList(dict, cSearch);
 }
 
-
+bool pFlow::boundaryContactSearchList::boundariesUpdated() const
+{
+    return boundaries_.boundariesUpdated();
+}

src/Interaction/contactSearch/boundaries/boundaryContactSearchList.hpp

@@ -1,4 +1,24 @@
-#include "ListPtr.hpp"
+/*------------------------------- 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 "boundaryListPtr.hpp"
 #include "boundaryContactSearch.hpp"
 
 namespace pFlow
@@ -9,7 +29,7 @@ class contactSearch;
 
 class boundaryContactSearchList
 :
-    public ListPtr<boundaryContactSearch>
+    public boundaryListPtr<boundaryContactSearch>
 {
 private:
 
@@ -28,6 +48,14 @@ public:
         const contactSearch& cSearch);
     
     ~boundaryContactSearchList()=default;
+
+    inline
+    const boundaryList& boundaries()const
+    {
+        return boundaries_;
+    } 
+
+    bool boundariesUpdated()const;
     
 };
 

src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/periodicBoundaryContactSearch.hpp

@@ -74,6 +74,11 @@ public:
 		csPairContainerType &ppPairs,
 		csPairContainerType &pwPairs,
 		bool force = false) override;
+	
+	bool isActive()const override
+	{
+		return true;
+	}
 };
 
 }

src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp

@@ -11,7 +11,7 @@ pFlow::wallBoundaryContactSearch::wallBoundaryContactSearch
     const ViewType1D<realx3, memory_space> &normals
 )
 :
-    cellExtent_( max(cellExtent, 0.5 ) ),
+    cellExtent_( max(cellExtent, half ) ),
     numElements_(numElements),
     numPoints_(numPoints),
     vertices_(vertices),

src/Interaction/contactSearch/contactSearch/contactSearch.cpp

@@ -31,11 +31,11 @@ pFlow::contactSearch::contactSearch(
  	Timers& timers)
 :
 	extendedDomainBox_(extDomain),
+	updateInterval_(dict.getValMax<uint32>("updateInterval", 1u)),
 	particles_(prtcl),
 	geometry_(geom),
-	ppTimer_("particle-particle contact search", &timers),
-	pwTimer_("particle-wall contact search", &timers),
-	dict_(dict)
+	bTimer_("Boundary particles contact search", &timers),
+	ppTimer_("Internal particles contact search", &timers)
 {
 
 }
@@ -45,6 +45,78 @@ const pFlow::pointStructure &pFlow::contactSearch::pStruct() const
    return particles_.pStruct();
 }
 
+bool pFlow::contactSearch::broadSearch
+(
+	const timeInfo &ti, 
+	csPairContainerType &ppPairs, 
+	csPairContainerType &pwPairs, 
+	bool force
+)
+{
+
+	if(enterBroadSearch(ti, force))
+	{
+		ppTimer_.start();
+		if( !BroadSearch(
+			ti,
+			ppPairs,
+			pwPairs,
+			force ) )
+		{
+			fatalErrorInFunction;
+			performedSearch_ = false;
+			return false;
+		}
+		ppTimer_.end();
+		performedSearch_ = true;
+		lastUpdated_ = ti.currentIter();
+	}
+	else
+	{
+		performedSearch_ = false;
+	}
+
+	return true;
+}
+
+bool pFlow::contactSearch::boundaryBroadSearch
+(
+	uint32 bndryIndex, 
+	const timeInfo &ti, 
+	csPairContainerType &ppPairs, 
+	csPairContainerType &pwPairs, 
+	bool force
+)
+{
+	if(enterBroadSearchBoundary(ti, force))
+	{
+		bTimer_.start();
+		for(uint32 i=0u; i<6u; i++)
+		{
+			//output<<" boundarySearch "<< i <<" for iter "<< ti.iter()<<endl;
+			if(!BoundaryBroadSearch(
+				i,
+				ti,
+				ppPairs,
+				pwPairs,
+				force))
+			{
+				performedSearchBoundary_ = false;
+				return false;
+			}
+		}
+		bTimer_.end();
+		performedSearchBoundary_ = true;
+	}
+	else
+	{
+
+		performedSearchBoundary_ = false;
+	}
+
+    return true;
+}
+
 pFlow::uniquePtr<pFlow::contactSearch> pFlow::contactSearch::create(
     const dictionary &dict,
     const box &extDomain,

src/Interaction/contactSearch/contactSearch/contactSearch.hpp

@@ -26,6 +26,7 @@ Licence:
 #include "contactSearchGlobals.hpp"
 #include "dictionary.hpp"
 #include "virtualConstructor.hpp"
+#include "timeInfo.hpp"
 #include "Timer.hpp"
 
 namespace pFlow
@@ -44,15 +45,44 @@ private:
 
 	const box& 	extendedDomainBox_;
 
+	/// @brief update interval in terms of iteration numebr 
+	uint32 		updateInterval_= 1;
+
+	/// @brief  last iteration number which contact search has been performed
+	uint32 		lastUpdated_ 	= 0;
+
+	/// @brief performed search? 
+	bool 		performedSearch_ = false;
+
+	/// @brief performed search in boundaries
+	bool 		performedSearchBoundary_ = false;
+
+	/// const ref to particles
 	const particles& 		particles_;
 
+	/// const ref to geometry 
 	const geometry&			geometry_;
 
+	Timer 		bTimer_;
+
 	Timer 		ppTimer_;
 
-	Timer 		pwTimer_;
+	virtual
+	bool BroadSearch(
+		const timeInfo& ti,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force
+	)=0;
 
-	dictionary 	dict_;
+	virtual
+	bool  BoundaryBroadSearch(
+		uint32 bndryIndex,
+		const timeInfo& ti,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force = false
+	)=0;
 
 public:
 
@@ -82,66 +112,86 @@ public:
 	 	(dict, domain, prtcl, geom, timers)
 	);
 
-	
-	const auto& dict()const
+	inline
+	bool performedSearch()const
 	{
-		return dict_;
+		return performedSearch_;
 	}
 
-	const auto& extendedDomainBox()const
+	inline 
+	bool performedSearchBoundary()const 
+	{
+		return performedSearchBoundary_;
+	}
+
+	inline
+	bool performSearch(uint32 iter, bool force = false)const
+	{
+		if((iter-lastUpdated_) % updateInterval_ == 0 || iter == 0 || force )
+		{
+			return true;
+		}
+		return false;		
+	}
+
+	bool enterBroadSearch(const timeInfo& ti, bool force = false)const
+	{
+		return performSearch(ti.iter(), force);
+	}
+
+	virtual 
+	bool enterBroadSearchBoundary(const timeInfo& ti, bool force=false)const = 0;
+
+	inline
+	uint32 updateInterval()const
+	{
+		return updateInterval_;
+	}
+
+	inline
+	const box& extendedDomainBox()const
 	{
 		return extendedDomainBox_;
 	}
 
-	const auto& Particles()const
+	inline
+	const particles& Particles()const
 	{
 		return particles_;
 	}
 
 	const pointStructure& pStruct()const;
 
-	const auto& Geometry()const
+	inline
+	const geometry& Geometry()const
 	{
 		return geometry_;
 	}
 
-	auto& ppTimer()
+	inline
+	Timer& ppTimer()
 	{
 		return ppTimer_;
 	}
 
-	auto& pwTimer()
+	inline
+	Timer& bTimer()
 	{
-		return pwTimer_;
+		return bTimer_;
 	}
-
-	virtual 
+ 
 	bool broadSearch(
-		uint32 iter,
-		real t,
-		real dt,
+		const timeInfo& ti,
 		csPairContainerType& ppPairs,
 		csPairContainerType& pwPairs,
-		bool force = false) = 0;
+		bool force = false);
 	
-	virtual 
 	bool boundaryBroadSearch(
-		uint32 i,
-		uint32 iter,
-		real t,
-		real dt,
+		uint32 bndryIndex,
+		const timeInfo& ti,
 		csPairContainerType& ppPairs,
 		csPairContainerType& pwPairs,
-		bool force = false)=0;
-
-	virtual 
-	bool enterBroadSearch(uint32 iter, real t, real dt)const = 0;
-
-	virtual 
-	bool performedBroadSearch()const = 0;
-
-	virtual
-	uint32 updateInterval()const = 0;
+		bool force = false);
 
 	virtual 
 	real sizeRatio()const = 0;

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

@@ -44,9 +44,9 @@ pFlow::NBS::NBS
         position,
         flags,
         diam),
-    sizeRatio_(max(dict.getVal<real>("sizeRatio"), 1.0)),
-    cellExtent_(max(dict.getVal<real>("cellExtent"), 0.5)),
-    adjustableBox_(dict.getVal<Logical>("adjustableBox")),
+    sizeRatio_(max(dict.getVal<real>("sizeRatio"), one)),
+    cellExtent_(max(dict.getVal<real>("cellExtent"), half)),
+    adjustableBox_(false),//adjustableBox_(dict.getVal<Logical>("adjustableBox")),
     NBSLevel0_
 	(
         this->domainBox_,

src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp

@@ -31,7 +31,7 @@ pFlow::cellsWallLevel0::cellsWallLevel0
 	const ViewType1D<realx3, memory_space>& normals
 )
 :
-    cellExtent_( max(cellExtent, 0.5 ) ),
+    cellExtent_( max(cellExtent, half ) ),
     numElements_(numElements),
     numPoints_(numPoints),
     vertices_(vertices),

src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBS.cpp

@@ -23,7 +23,7 @@ Licence:
 #include "streams.hpp"
 
 pFlow::uint32 pFlow::mapperNBS::checkInterval_ = 1000;
-pFlow::real pFlow::mapperNBS::enlargementFactor_ = 1.1;
+pFlow::real pFlow::mapperNBS::enlargementFactor_ = 1.5;
 
 bool pFlow::mapperNBS::setSearchBox
 (

src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp

@@ -21,7 +21,7 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "mapperNBSKernels.hpp"
-
+#include "streams.hpp"
 
 void pFlow::mapperNBSKernels::findPointExtends
 (
@@ -153,16 +153,19 @@ bool pFlow::mapperNBSKernels::buildLists
     const pFlagTypeDevice &flags
 )
 {
+    
     auto aRange = flags.activeRange();
+    
+    auto pp = points;
     if(flags.isAllActive() )
-    {
+    { 
         Kokkos::parallel_for
         (
             "pFlow::mapperNBSKernels::buildLists",
             deviceRPolicyStatic(aRange.start(), aRange.end()),
             LAMBDA_HD(uint32 i)
-            {
-                auto ind = searchCell.pointIndex(points[i]);	
+            {               
+                auto ind = searchCell.pointIndex(pp[i]);
                 uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
                 next[i] = old;
             }

src/Interaction/contactSearch/methods/particleWallContactSearchs.cpp

@@ -31,48 +31,36 @@ pFlow::particleWallContactSearchs<method>::particleWallContactSearchs
     const ViewType1D<real, memory_space> &diam
 )
 :
-    domainBox_(domain),
-    updateInterval_
-    (
-        max(dict.getValOrSet<uint32>("updateInterval", 1),1u)
-    )
-{
-    
-}
+    domainBox_(domain)
+{}
 
 template <typename method>
 bool pFlow::particleWallContactSearchs<method>::broadSearch
-	(
-		uint32 iter,
-		real t,
-		real dt,
-		csPairContainerType& ppPairs,
-		csPairContainerType& pwPairs,
-		const deviceViewType1D<realx3>& pointPos,
-		const pFlagTypeDevice& flags,
-		const deviceViewType1D<real>& diameter,
-		bool  force
-	)
+(
+	uint32 iter,
+	real t,
+	real dt,
+	csPairContainerType& ppPairs,
+	csPairContainerType& pwPairs,
+	const deviceViewType1D<realx3>& pointPos,
+	const pFlagTypeDevice& flags,
+	const deviceViewType1D<real>& diameter,
+	bool  force
+)
+{
+	
+	if(!getMethod().impl_broadSearch(
+		ppPairs,
+		pwPairs,
+		pointPos,
+		flags,
+		diameter))
 	{
-		
-		performedSearch_ = false;
-		if( !performSearch(iter, force) ) return true;
-
-		if(!getMethod().impl_broadSearch(
-			ppPairs,
-			pwPairs,
-			pointPos,
-			flags,
-			diameter))
-		{
-			fatalErrorInFunction<<
-			"Error in performing particle-particle broadSearch in method"<<
-			getMethod().typeName()<<endl;
-			return false;
-		}
-
-		lastUpdated_ = iter;
-		
-		performedSearch_ = true;		
-		return true;
-	}
+		fatalErrorInFunction<<
+		"Error in performing particle-particle broadSearch in method"<<
+		getMethod().typeName()<<endl;
+		return false;
+	}
+	
+	return true;
+}

src/Interaction/contactSearch/methods/particleWallContactSearchs.hpp

@@ -48,18 +48,6 @@ private:
 	/// @brief  box enclosing the simulation domain (local to processor)
 	box 		domainBox_;
 
-	/*/// @brief box enclosing the area for contact search (region with points)
-	box 		searchBox_;*/
-
-	/// @brief update interval in terms of iteration numebr 
-	uint32 		updateInterval_= 1;
-
-	/// @brief  last iteration number which contact search has been performed
-	uint32 		lastUpdated_ 	= 0;
-
-	/// @brief performed search? 
-	bool 		performedSearch_ = false;
-
 protected:
 
 	inline
@@ -98,25 +86,6 @@ public:
     	const deviceViewType1D<real>& diameter,
 		bool  force = false
 	);
-		
-	bool performedSearch()const
-	{
-		return performedSearch_;
-	}
-
-	bool performSearch(uint32 iter, bool force = false)const
-	{
-		if((iter-lastUpdated_) % updateInterval_ == 0 || iter == 0 || force )
-		{
-			return true;
-		}
-		return false;		
-	}
-
-	uint32 updateInterval()const
-	{
-		return updateInterval_;
-	}
 
 	real sizeRatio()const
 	{

src/Interaction/grainInteraction/boundaries/boundaryGrainInteraction/boundaryGrainInteraction.cpp

@@ -0,0 +1,104 @@
+#include "boundarySphereInteraction.hpp"
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+template <typename cFM, typename gMM>
+void pFlow::boundaryGrainInteraction<cFM, gMM>::allocatePPPairs()
+{
+	ppPairs_.reset(nullptr);
+	ppPairs_ = makeUnique<ContactListType>(1);
+}
+
+template <typename cFM, typename gMM>
+void pFlow::boundaryGrainInteraction<cFM, gMM>::allocatePWPairs()
+{
+	pwPairs_.reset(nullptr);
+	pwPairs_ = makeUnique<ContactListType>(1);
+}
+
+
+template <typename cFM, typename gMM>
+pFlow::boundaryGrainInteraction<cFM, gMM>::boundaryGrainInteraction(
+	const boundaryBase &boundary,
+	const grainParticles &grnPrtcls,
+	const GeometryMotionModel &geomMotion)
+	: generalBoundary(boundary, grnPrtcls.pStruct(), "", ""),
+	  geometryMotion_(geomMotion),
+	  grnParticles_(grnPrtcls)
+{
+}
+
+template <typename cFM, typename gMM>
+pFlow::uniquePtr<pFlow::boundaryGrainInteraction<cFM, gMM>>
+pFlow::boundaryGrainInteraction<cFM, gMM>::create(
+	const boundaryBase &boundary,
+	const grainParticles &grnPrtcls,
+	const GeometryMotionModel &geomMotion)
+{
+	word cfTypeName = ContactForceModel::TYPENAME();
+	word gmTypeName = MotionModel::TYPENAME();
+	word bType = boundary.type();
+
+	word boundaryTypeName = angleBracketsNames3(
+		"boundaryGrainInteraction",
+		bType,
+		cfTypeName,
+		gmTypeName);
+
+	word altBTypeName = angleBracketsNames2(
+		"boundaryGrainInteraction",
+		cfTypeName,
+		gmTypeName);
+
+	if (boundaryBasevCtorSelector_.search(boundaryTypeName))
+	{
+		pOutput.space(4) << "Creating boundry type " << Green_Text(boundaryTypeName) << 
+		" for boundary " << boundary.name() << " . . ." << END_REPORT;
+		return boundaryBasevCtorSelector_[boundaryTypeName](
+			boundary,
+			grnPrtcls,
+			geomMotion);
+	}
+	else if(boundaryBasevCtorSelector_[altBTypeName])
+	{
+		// if boundary condition is not implemented, the default is used
+		
+		pOutput.space(4) << "Creating boundry type " << Green_Text(altBTypeName) << 
+		" for boundary " << boundary.name() << " . . ." << END_REPORT;
+		return boundaryBasevCtorSelector_[altBTypeName](
+			boundary,
+			grnPrtcls,
+			geomMotion);
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< boundaryTypeName<<
+			" and "<< altBTypeName << " do not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			boundaryBasevCtorSelector_
+		);
+		fatalExit;
+	}
+
+	return nullptr;
+}

src/Interaction/grainInteraction/boundaries/boundaryGrainInteraction/boundaryGrainInteraction.hpp

@@ -0,0 +1,188 @@
+/*------------------------------- 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 __boundaryGrainInteraction_hpp__
+#define __boundaryGrainInteraction_hpp__
+
+#include "virtualConstructor.hpp"
+#include "generalBoundary.hpp"
+#include "sortedContactList.hpp"
+#include "grainParticles.hpp"
+
+namespace pFlow
+{
+
+template<typename contactForceModel,typename geometryMotionModel>
+class boundaryGrainInteraction
+:
+    public generalBoundary
+{
+public:
+	
+	using BoundaryGrainInteractionType
+			= boundaryGrainInteraction<contactForceModel, geometryMotionModel>;
+
+    using GeometryMotionModel 	= geometryMotionModel;
+	
+	using ContactForceModel 	= contactForceModel;
+
+	using MotionModel 			= typename geometryMotionModel::MotionModel;
+
+	using ModelStorage 			= typename ContactForceModel::contactForceStorage;
+	
+	using IdType 			    = uint32;
+
+	using IndexType    		    = uint32;
+    
+    using ContactListType = 
+        sortedContactList<ModelStorage, DefaultExecutionSpace, IdType>;
+
+private:
+
+	const GeometryMotionModel& 			geometryMotion_;
+
+	/// const reference to sphere particles 
+	const grainParticles& 				grnParticles_;
+    
+	uniquePtr<ContactListType> 			ppPairs_ = nullptr;
+
+	uniquePtr<ContactListType>			pwPairs_ = nullptr;
+
+protected:
+
+	void allocatePPPairs();
+
+	void allocatePWPairs();
+
+public:
+
+	TypeInfoTemplate12("boundaryGrainInteraction", ContactForceModel, MotionModel);
+
+	boundaryGrainInteraction(
+		const boundaryBase& boundary,
+		const grainParticles& grnPrtcls,
+		const GeometryMotionModel& geomMotion);
+
+	create_vCtor
+	(
+		BoundaryGrainInteractionType,
+		boundaryBase,
+		(
+			const boundaryBase& boundary,
+			const grainParticles& grnPrtcls,
+			const GeometryMotionModel& geomMotion
+		),
+		(boundary, grnPrtcls, geomMotion)
+	);
+
+	add_vCtor
+	(
+		BoundaryGrainInteractionType,
+		BoundaryGrainInteractionType,
+		boundaryBase
+	);
+
+	~boundaryGrainInteraction()override=default;
+
+	const auto& grnParticles()const
+	{
+		return grnParticles_;
+	}
+
+	const auto& geometryMotion()const
+	{
+		return geometryMotion_;
+	}
+
+	ContactListType& ppPairs()
+	{
+		return ppPairs_();
+	}
+
+	const ContactListType& ppPairs()const
+	{
+		return ppPairs_();
+	}
+
+	ContactListType& pwPairs()
+	{
+		return pwPairs_();
+	}
+
+	const ContactListType& pwPairs()const
+	{
+		return pwPairs_();
+	}
+
+	bool ppPairsAllocated()const
+	{
+		if( ppPairs_)return true;
+		return false;
+	}
+
+	bool pwPairsAllocated()const
+	{
+		if( pwPairs_)return true;
+		return false;
+	}
+
+	virtual
+	bool grainGrainInteraction(
+		real dt,
+		const ContactForceModel& cfModel,
+		uint32 step)
+	{
+		// for default boundary, no thing to be done 
+		return false;
+	}
+
+	bool hearChanges
+	(
+		const timeInfo& ti,
+		const message& msg, 
+    	const anyList& varList
+	) override
+    {
+		if(msg.equivalentTo(message::BNDR_RESET))
+		{
+			return true;
+		}
+		fatalErrorInFunction<<"Event "<< msg.eventNames() << " is not handled !"<<endl;
+		return false;
+	}
+
+	bool isActive()const override
+	{
+		return false;
+	}
+
+	static
+	uniquePtr<BoundaryGrainInteractionType> create(
+		const boundaryBase& boundary,
+		const grainParticles& sphPrtcls,
+		const GeometryMotionModel& geomMotion
+	);
+
+};
+
+}
+
+#include "boundaryGrainInteraction.cpp"
+
+#endif //__boundaryGrainInteraction_hpp__

src/Interaction/grainInteraction/boundaries/boundaryGrainInteractionList.cpp

@@ -0,0 +1,23 @@
+
+
+template <typename CFModel, typename gMModel>
+pFlow::boundaryGrainInteractionList<CFModel, gMModel>::boundaryGrainInteractionList
+(
+    const grainParticles &grnPrtcls, 
+    const gMModel &geomMotion
+)
+:
+	boundaryListPtr<boundaryGrainInteraction<CFModel,gMModel>>(),
+	boundaries_(grnPrtcls.pStruct().boundaries())
+{
+	//gSettings::sleepMiliSeconds(1000*pFlowProcessors().localRank());
+	ForAllBoundariesPtr(i, this)
+	{
+		this->set(
+			i,
+			boundaryGrainInteraction<CFModel, gMModel>::create(
+				boundaries_[i],
+				grnPrtcls, 
+				geomMotion));
+	}	
+}

src/Interaction/grainInteraction/boundaries/boundaryGrainInteractionList.hpp

@@ -0,0 +1,40 @@
+#ifndef __boundaryGrainInteractionList_hpp__
+#define __boundaryGrainInteractionList_hpp__
+
+
+#include "boundaryList.hpp"
+#include "boundaryListPtr.hpp"
+#include "boundaryGrainInteraction.hpp"
+
+
+namespace pFlow
+{
+
+
+template<typename contactForceModel,typename geometryMotionModel>
+class boundaryGrainInteractionList
+:
+    public boundaryListPtr<boundaryGrainInteraction<contactForceModel,geometryMotionModel>>
+{
+private:
+
+    const boundaryList&         boundaries_;
+
+public:
+
+    boundaryGrainInteractionList(
+        const grainParticles& grnPrtcls,
+		const geometryMotionModel& geomMotion
+    );
+
+    ~boundaryGrainInteractionList()=default;
+
+};
+
+
+
+}
+
+#include "boundaryGrainInteractionList.cpp"
+
+#endif //__boundaryGrainInteractionList_hpp__

src/Interaction/grainInteraction/boundaries/createBoundaryGrainInteraction.hpp

@@ -15,51 +15,17 @@ Licence:
   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 "boundaryGrainInteraction.hpp"
+#include "periodicBoundaryGrainInteraction.hpp"
 
-#ifndef __eventObserver_hpp__
-#define __eventObserver_hpp__
+#define createBoundaryGrainInteraction(ForceModel,GeomModel) 	\
+	template class pFlow::boundaryGrainInteraction<			\
+	ForceModel,													\
+	GeomModel>;													\
+																\
+	template class pFlow::periodicBoundaryGrainInteraction<	\
+	ForceModel,													\
+	GeomModel>;
 
-#include "eventMessage.hpp"
-
-namespace pFlow
-{
-
-class eventSubscriber;
-
-
-class eventObserver
-{
-protected:
-	
-	const eventSubscriber* subscriber_ = nullptr;
-	
-	// if this object is linked to subscriber 
-	bool subscribed_ = false;
-
-public:
-
-	eventObserver();
-
-	eventObserver(const eventSubscriber& subscriber, bool subscribe = true );
-
-	virtual ~eventObserver();
-
-	inline bool subscribed()const {return subscribed_;}
-
-	bool subscribe(const eventSubscriber& subscriber);
-	
-	inline void invalidateSubscriber()
-	{
-		subscribed_ = false;
-	}
-
-	virtual bool update(const eventMessage& msg)=0;
-};
-
-} // pFlow
-
-
-#endif // __eventObserver_hpp__

src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundaryGrainInteraction.cpp

@@ -0,0 +1,70 @@
+/*------------------------------- 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 "periodicBoundarySIKernels.hpp"
+
+template <typename cFM, typename gMM>
+pFlow::periodicBoundaryGrainInteraction<cFM, gMM>::periodicBoundaryGrainInteraction(
+	const boundaryBase &boundary,
+	const grainParticles &grnPrtcls,
+	const GeometryMotionModel &geomMotion)
+	: boundaryGrainInteraction<cFM,gMM>(boundary, grnPrtcls, geomMotion),
+	transferVec_(boundary.mirrorBoundary().displacementVectroToMirror())
+{
+	if(boundary.thisBoundaryIndex()%2==1)
+    {
+    	masterInteraction_ = true;
+		this->allocatePPPairs();
+		this->allocatePWPairs();
+	   
+    }
+    else
+    {
+        masterInteraction_ = false;
+    }
+}
+
+template <typename cFM, typename gMM>
+bool pFlow::periodicBoundaryGrainInteraction<cFM, gMM>::grainGrainInteraction
+(
+	real dt,
+	const ContactForceModel &cfModel,
+	uint32 step
+)
+{
+	if(!masterInteraction_) return false;
+	
+	pFlow::periodicBoundarySIKernels::grainGrainInteraction(
+		dt,
+		this->ppPairs(),
+		cfModel,
+		transferVec_,
+		this->boundary().thisPoints(),
+		this->mirrorBoundary().thisPoints(),
+		this->grnParticles().diameter().deviceViewAll(),
+		this->grnParticles().coarseGrainFactor().deviceViewAll(),
+		this->grnParticles().propertyId().deviceViewAll(),
+		this->grnParticles().velocity().deviceViewAll(),
+		this->grnParticles().rVelocity().deviceViewAll(),
+		this->grnParticles().contactForce().deviceViewAll(),
+		this->grnParticles().contactTorque().deviceViewAll());
+
+	return false;
+}

src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundaryGrainInteraction.hpp

@@ -0,0 +1,98 @@
+/*------------------------------- 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 __periodicBoundaryGrainInteraction_hpp__
+#define __periodicBoundaryGrainInteraction_hpp__
+
+#include "boundaryGrainInteraction.hpp"
+
+namespace pFlow
+{
+
+template<typename contactForceModel,typename geometryMotionModel>
+class periodicBoundaryGrainInteraction
+:
+    public boundaryGrainInteraction<contactForceModel, geometryMotionModel>
+{
+public:
+
+    using PBSInteractionType = 
+        periodicBoundaryGrainInteraction<contactForceModel,geometryMotionModel>;
+
+    using BSInteractionType = 
+        boundaryGrainInteraction<contactForceModel, geometryMotionModel>;
+
+    using GeometryMotionModel 	= typename BSInteractionType::GeometryMotionModel;
+	
+	using ContactForceModel 	= typename BSInteractionType::ContactForceModel;
+
+	using MotionModel 			= typename geometryMotionModel::MotionModel;
+
+	using ModelStorage 			= typename ContactForceModel::contactForceStorage;
+	
+	using IdType 			    = typename BSInteractionType::IdType;
+
+	using IndexType    		    = typename BSInteractionType::IndexType;
+    
+    using ContactListType       = typename BSInteractionType::ContactListType; 
+
+private:
+
+    realx3  transferVec_;
+    
+    bool    masterInteraction_;
+public:
+
+    TypeInfoTemplate22("boundaryGrainInteraction", "periodic",ContactForceModel, MotionModel);
+
+	
+    periodicBoundaryGrainInteraction(
+        const boundaryBase& boundary,
+		const grainParticles& grnPrtcls,
+		const GeometryMotionModel& geomMotion
+    );
+
+    add_vCtor
+    (
+        BSInteractionType,
+        PBSInteractionType,
+        boundaryBase
+    );
+    
+    ~periodicBoundaryGrainInteraction()override = default;
+
+    bool grainGrainInteraction(
+        real dt,
+		const ContactForceModel& cfModel,
+        uint32 step)override;
+    
+    bool isActive()const override
+	{
+		return true;
+	}
+	
+};
+
+}
+
+#include "periodicBoundaryGrainInteraction.cpp"
+
+
+#endif //__periodicBoundaryGrainInteraction_hpp__

src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundarySIKernels.hpp

@@ -0,0 +1,121 @@
+
+namespace pFlow::periodicBoundarySIKernels
+{
+
+template<typename ContactListType, typename ContactForceModel>
+inline
+void grainGrainInteraction
+(
+    real dt,
+    const ContactListType&              cntctList,
+    const ContactForceModel&            forceModel,
+    const realx3&                       transferVec,
+    const deviceScatteredFieldAccess<realx3>& thisPoints,
+    const deviceScatteredFieldAccess<realx3>& mirrorPoints,
+    const deviceViewType1D<real>& diam,
+    const deviceViewType1D<real>& coarseGrainFactor,
+    const deviceViewType1D<uint32>& propId,
+    const deviceViewType1D<realx3>& vel,
+    const deviceViewType1D<realx3>& rVel,
+    const deviceViewType1D<realx3>& cForce,
+    const deviceViewType1D<realx3>& cTorque
+)
+{
+    
+    using ValueType = typename ContactListType::ValueType;
+    uint32 ss = cntctList.size();
+    uint32 lastItem = cntctList.loopCount();
+    if(lastItem == 0u)return;
+
+    Kokkos::parallel_for(
+        "pFlow::periodicBoundarySIKernels::grainGrainInteraction",
+        deviceRPolicyDynamic(0,lastItem),
+        LAMBDA_HD(uint32 n)
+        {
+
+        if(!cntctList.isValid(n))return;
+
+            auto [i,j] = cntctList.getPair(n);
+            uint32 ind_i = thisPoints.index(i);
+            uint32 ind_j = mirrorPoints.index(j);
+
+            real Ri = 0.5*diam[ind_i];
+            real Rj = 0.5*diam[ind_j];
+            real cGFi = coarseGrainFactor[ind_i];
+            real cGFj = coarseGrainFactor[ind_j];
+            realx3 xi = thisPoints.field()[ind_i];
+            realx3 xj = mirrorPoints.field()[ind_j]+transferVec;
+            real dist = length(xj-xi);
+            real ovrlp = (Ri+Rj) - dist;
+            
+            if( ovrlp >0.0 )
+            {
+                auto Nij = (xj-xi)/dist;
+                auto wi = rVel[ind_i];
+                auto wj = rVel[ind_j];
+                auto Vr = vel[ind_i] - vel[ind_j] + cross((Ri*wi+Rj*wj), Nij);
+                
+                auto history = cntctList.getValue(n);
+
+                int32 propId_i = propId[ind_i];
+                int32 propId_j = propId[ind_j];
+
+                realx3 FCn, FCt, Mri, Mrj, Mij, Mji;
+
+                // calculates contact force 
+                forceModel.contactForce(
+                    dt, i, j,
+                    propId_i, propId_j,
+                    Ri, Rj, cGFi , cGFj ,
+                    ovrlp,
+                    Vr, Nij,
+                    history,
+                    FCn, FCt);
+
+                forceModel.rollingFriction(
+                    dt, i, j,
+                    propId_i, propId_j,
+                    Ri, Rj, cGFi , cGFj ,
+                    wi, wj,
+                    Nij,
+                    FCn,
+                    Mri, Mrj);
+
+                auto M = cross(Nij,FCt);
+                Mij = Ri*M+Mri;
+                Mji = Rj*M+Mrj;
+                
+                auto FC = FCn + FCt;
+                
+
+                Kokkos::atomic_add(&cForce[ind_i].x_,FC.x_);
+                Kokkos::atomic_add(&cForce[ind_i].y_,FC.y_);
+                Kokkos::atomic_add(&cForce[ind_i].z_,FC.z_);
+
+                Kokkos::atomic_add(&cForce[ind_j].x_,-FC.x_);
+                Kokkos::atomic_add(&cForce[ind_j].y_,-FC.y_);
+                Kokkos::atomic_add(&cForce[ind_j].z_,-FC.z_);
+
+                Kokkos::atomic_add(&cTorque[ind_i].x_, Mij.x_);
+                Kokkos::atomic_add(&cTorque[ind_i].y_, Mij.y_);
+                Kokkos::atomic_add(&cTorque[ind_i].z_, Mij.z_);
+
+                Kokkos::atomic_add(&cTorque[ind_j].x_, Mji.x_);
+                Kokkos::atomic_add(&cTorque[ind_j].y_, Mji.y_);
+                Kokkos::atomic_add(&cTorque[ind_j].z_, Mji.z_);
+                
+
+                cntctList.setValue(n,history);
+
+            }
+            else
+            {
+                cntctList.setValue(n, ValueType());
+            }
+
+        });
+    Kokkos::fence();
+}
+
+
+} //pFlow::periodicBoundarySIKernels

src/Interaction/grainInteraction/grainInteraction/grainInteraction.cpp

@@ -0,0 +1,355 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+template<typename cFM,typename gMM,template <class, class, class> class cLT>
+bool pFlow::grainInteraction<cFM,gMM, cLT>::createGrainInteraction()
+{
+	
+	realVector_D rhoD("densities", this->densities());
+
+	auto modelDict = this->subDict("model");
+
+	forceModel_ = makeUnique<ContactForceModel>(
+		this->numMaterials(),
+		rhoD.deviceView(),
+		modelDict );
+
+
+	uint32 nPrtcl = grnParticles_.size();
+
+	contactSearch_ = contactSearch::create(
+	subDict("contactSearch"),
+	grnParticles_.extendedDomain().domainBox(),
+	grnParticles_,
+	geometryMotion_,
+	timers());
+
+	ppContactList_ = makeUnique<ContactListType>("Grain-Grain",nPrtcl+1);
+	
+	pwContactList_ = makeUnique<ContactListType>("Grain-wall",nPrtcl/5+1);
+
+	return true;
+}
+
+
+template<typename cFM,typename gMM,template <class, class, class> class cLT>
+bool pFlow::grainInteraction<cFM,gMM, cLT>::grainGrainInteraction()
+{
+	auto lastItem = ppContactList_().loopCount();
+
+	// create the kernel functor 
+	pFlow::grainInteractionKernels::ppInteractionFunctor 
+		ppInteraction(
+			this->dt(),
+			this->forceModel_(),
+			ppContactList_(), // to be read
+			grnParticles_.diameter().deviceViewAll(),
+			grnParticles_.coarseGrainFactor().deviceViewAll(),
+			grnParticles_.propertyId().deviceViewAll(),
+			grnParticles_.pointPosition().deviceViewAll(),
+			grnParticles_.velocity().deviceViewAll(),
+			grnParticles_.rVelocity().deviceViewAll(),
+			grnParticles_.contactForce().deviceViewAll(),
+			grnParticles_.contactTorque().deviceViewAll()
+			);
+	
+	Kokkos::parallel_for(
+		"ppInteraction",
+		rpPPInteraction(0,lastItem),
+		ppInteraction
+		);
+
+	Kokkos::fence();	
+
+	return true;
+}
+
+
+template<typename cFM,typename gMM,template <class, class, class> class cLT>
+bool pFlow::grainInteraction<cFM,gMM, cLT>::grainWallInteraction()
+{
+	
+	uint32 lastItem = pwContactList_().loopCount();
+	uint32 iter = this->currentIter();
+	real t = this->currentTime();
+	real dt = this->dt();
+
+	pFlow::grainInteractionKernels::pwInteractionFunctor
+		pwInteraction(
+			dt,
+			this->forceModel_(),
+			pwContactList_(),
+			geometryMotion_.getTriangleAccessor(), 
+			geometryMotion_.getModel(iter, t, dt) , 
+			grnParticles_.diameter().deviceViewAll() ,
+			grnParticles_.coarseGrainFactor().deviceViewAll() ,
+			grnParticles_.propertyId().deviceViewAll(),
+			grnParticles_.pointPosition().deviceViewAll(),
+			grnParticles_.velocity().deviceViewAll(),
+			grnParticles_.rVelocity().deviceViewAll() ,
+			grnParticles_.contactForce().deviceViewAll(),
+			grnParticles_.contactTorque().deviceViewAll() ,
+			geometryMotion_.triMotionIndex().deviceViewAll(),
+			geometryMotion_.propertyId().deviceViewAll(),
+			geometryMotion_.contactForceWall().deviceViewAll()
+		);
+
+	Kokkos::parallel_for(
+		"",
+		rpPWInteraction(0,lastItem),
+		pwInteraction
+		);
+
+	
+	Kokkos::fence();
+
+	return true;
+}
+
+
+template<typename cFM,typename gMM,template <class, class, class> class cLT>
+pFlow::grainInteraction<cFM,gMM, cLT>::grainInteraction
+(
+	systemControl& control,
+	const particles& prtcl,
+	const geometry& geom
+)
+:
+	interaction(control, prtcl, geom),
+	geometryMotion_(dynamic_cast<const GeometryMotionModel&>(geom)),
+	grnParticles_(dynamic_cast<const grainParticles&>(prtcl)),
+	boundaryInteraction_(grnParticles_,geometryMotion_),
+	ppInteractionTimer_("grain-graine interaction (internal)", &this->timers()),
+	pwInteractionTimer_("grain-wall interaction (internal)", &this->timers()),
+	boundaryInteractionTimer_("grain-grain interaction (boundary)",&this->timers()),
+	contactListMangementTimer_("list management (interal)", &this->timers()),
+	contactListMangementBoundaryTimer_("list management (boundary)", &this->timers())
+{
+	
+	if(!createGrainInteraction())
+	{
+		fatalExit;
+	}
+
+	for(uint32 i=0; i<6; i++)
+	{
+		activeBoundaries_[i] = boundaryInteraction_[i].ppPairsAllocated();
+	}
+}
+
+template<typename cFM,typename gMM,template <class, class, class> class cLT>
+bool pFlow::grainInteraction<cFM,gMM, cLT>::beforeIteration()
+{
+	return true;
+}
+
+template<typename cFM,typename gMM,template <class, class, class> class cLT>
+bool pFlow::grainInteraction<cFM,gMM, cLT>::iterate()
+{
+	
+	timeInfo ti = this->TimeInfo();
+	auto iter = ti.iter();
+	auto t = ti.t();
+	auto dt = ti.dt();
+
+	auto& contactSearchRef = contactSearch_();
+
+	bool broadSearch = contactSearchRef.enterBroadSearch(ti);
+	bool broadSearchBoundary = contactSearchRef.enterBroadSearchBoundary(ti);
+
+	/// update boundaries of the fields that are being used by inreaction 
+	grnParticles_.diameter().updateBoundaries(DataDirection::SlaveToMaster);
+	grnParticles_.velocity().updateBoundaries(DataDirection::SlaveToMaster);
+	grnParticles_.rVelocity().updateBoundaries(DataDirection::SlaveToMaster);
+	grnParticles_.propertyId().updateBoundaries(DataDirection::SlaveToMaster);
+	
+	/// lists 
+	if(broadSearch)
+	{
+		contactListMangementTimer_.start();
+		ComputationTimer().start();
+			ppContactList_().beforeBroadSearch();
+			pwContactList_().beforeBroadSearch();
+		ComputationTimer().end();
+		contactListMangementTimer_.pause();
+	} 
+
+	if(broadSearchBoundary)
+	{
+		contactListMangementBoundaryTimer_.start();
+		ComputationTimer().start();
+		for(uint32 i=0; i<6u; i++)
+		{
+			if(activeBoundaries_[i])
+			{
+				auto& BI = boundaryInteraction_[i];
+				BI.ppPairs().beforeBroadSearch();
+				BI.pwPairs().beforeBroadSearch();
+			} 
+		}
+		ComputationTimer().end();
+		contactListMangementBoundaryTimer_.pause();	
+	}
+	
+	if( grnParticles_.numActive()<=0)return true;
+	
+	ComputationTimer().start();
+	if( !contactSearchRef.broadSearch(
+			ti,		
+			ppContactList_(),
+			pwContactList_()) )
+	{
+		fatalErrorInFunction<<
+		"unable to perform broadSearch.\n";
+		fatalExit;
+	}
+	
+	for(uint32 i=0; i<6u; i++)
+	{
+		if(activeBoundaries_[i])
+		{
+			auto& BI = boundaryInteraction_[i];
+			if(!contactSearchRef.boundaryBroadSearch(
+				i,
+				ti,
+				BI.ppPairs(),
+				BI.pwPairs())
+			)
+			{
+				fatalErrorInFunction<<
+				"failed to perform broadSearch for boundary index "<<i<<endl;
+				return false;
+			}
+		}
+	}
+	ComputationTimer().end();
+	
+	if(broadSearch && contactSearchRef.performedSearch())
+	{
+		contactListMangementTimer_.resume();
+		ComputationTimer().start();
+			ppContactList_().afterBroadSearch();
+			pwContactList_().afterBroadSearch();
+		ComputationTimer().end();
+		contactListMangementTimer_.end();
+	}
+
+	if(broadSearchBoundary && contactSearchRef.performedSearchBoundary())
+	{
+		contactListMangementBoundaryTimer_.resume();
+		ComputationTimer().start();
+		for(uint32 i=0; i<6u; i++)
+		{
+			if(activeBoundaries_[i])
+			{
+				auto& BI = boundaryInteraction_[i];
+				BI.ppPairs().afterBroadSearch();
+				BI.pwPairs().afterBroadSearch();
+			}
+		}
+		ComputationTimer().end();
+		contactListMangementBoundaryTimer_.end();
+	}
+	
+	/// peform contact search on boundareis with active contact search (master boundaries)
+	auto requireStep = boundariesMask<6>(true);
+	const auto& cfModel = this->forceModel_();
+	uint32 step=1;
+	boundaryInteractionTimer_.start();
+	ComputationTimer().start();
+	while(requireStep.anyElement(true) && step <= 10)
+	{
+		for(uint32 i=0; i<6u; i++)
+		{
+			if(requireStep[i] )
+			{
+				requireStep[i] = boundaryInteraction_[i].grainGrainInteraction(
+					dt, 
+					this->forceModel_(),
+					step
+				);
+			}		
+		}
+		step++;
+	}
+	ComputationTimer().end();
+	boundaryInteractionTimer_.pause();
+	
+
+	ppInteractionTimer_.start();
+	ComputationTimer().start();
+		grainGrainInteraction();
+	ComputationTimer().end();
+	ppInteractionTimer_.end();
+
+	
+	pwInteractionTimer_.start();
+	ComputationTimer().start();
+		grainWallInteraction();
+	ComputationTimer().start();
+	pwInteractionTimer_.end();
+
+	{
+	boundaryInteractionTimer_.resume();
+	ComputationTimer().start();
+	auto requireStep = boundariesMask<6>(true);
+
+	uint32 step = 11;
+	const auto& cfModel = this->forceModel_();
+	while( requireStep.anyElement(true) && step < 20 )
+	{
+		for(uint32 i=0; i<6u; i++)
+		{
+			if(requireStep[i])
+			{
+				requireStep[i] = boundaryInteraction_[i].grainGrainInteraction(
+					dt, 
+					cfModel,
+					step
+				);
+			}		
+		}
+		step++;
+	}
+	ComputationTimer().end();
+	boundaryInteractionTimer_.end();
+	}
+	
+	return true;
+}
+
+template<typename cFM,typename gMM,template <class, class, class> class cLT>
+bool pFlow::grainInteraction<cFM,gMM, cLT>::afterIteration()
+{
+	return true;
+}
+
+template<typename cFM,typename gMM,template <class, class, class> class cLT>
+bool pFlow::grainInteraction<cFM,gMM, cLT>::hearChanges
+(
+	const timeInfo& ti,
+	const message& msg, 
+	const anyList& varList
+)
+{
+	fatalErrorInFunction<<"Event "<< msg.eventNames()<<
+	" is not handled in grainInteraction"<<endl;
+	return false;
+}

src/Interaction/grainInteraction/grainInteraction/grainInteraction.hpp

@@ -0,0 +1,167 @@
+/*------------------------------- 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 __grainInteraction_hpp__ 
+#define __grainInteraction_hpp__
+
+#include "interaction.hpp"
+#include "grainParticles.hpp"
+#include "boundaryGrainInteractionList.hpp"
+#include "grainInteractionKernels.hpp"
+#include "boundariesMask.hpp"
+#include "MPITimer.hpp"
+//#include "unsortedContactList.hpp"
+
+
+
+namespace pFlow
+{
+
+template<
+	typename contactForceModel,
+	typename geometryMotionModel,
+	template <class, class, class> class contactListType>
+class grainInteraction
+:
+	public interaction
+{
+public:
+
+	using GeometryMotionModel 	= geometryMotionModel;
+	
+	using ContactForceModel 	= contactForceModel;
+
+	using MotionModel 			= typename geometryMotionModel::MotionModel;
+
+	using ModelStorage 			= typename ContactForceModel::contactForceStorage;
+
+	using BoundaryListType 		= boundaryGrainInteractionList<ContactForceModel,GeometryMotionModel>;
+	
+	using IdType 			= uint32;
+
+	using IndexType    		= uint32;
+
+	using ContactListType 	= 
+			contactListType<ModelStorage, DefaultExecutionSpace, IdType>;
+
+	//using BoundaryContactListType = unsortedContactList<ModelStorage, DefaultExecutionSpace, IdType>;
+	
+
+
+private:
+	
+	/// const reference to geometry
+	const GeometryMotionModel& 			geometryMotion_;
+
+	/// const reference to particles 
+	const grainParticles& 				grnParticles_;
+
+	/// particle-particle and particle-wall interactions at boundaries 
+	BoundaryListType 					boundaryInteraction_;	
+
+	/// a mask for active boundaries (boundaries with intreaction)
+	boundariesMask<6> 					activeBoundaries_;
+
+	/// contact search object for pp and pw interactions 
+	uniquePtr<contactSearch> 			contactSearch_ = nullptr;
+
+	/// contact force model 
+	uniquePtr<ContactForceModel>  		forceModel_ 	= nullptr;
+
+	/// contact list for particle-particle interactoins (keeps the history)
+	uniquePtr<ContactListType> 			ppContactList_ = nullptr;
+
+	/// contact list for particle-wall interactions (keeps the history)
+	uniquePtr<ContactListType> 			pwContactList_ = nullptr;
+
+
+	/// timer for particle-particle interaction computations
+	Timer 		ppInteractionTimer_;
+
+	/// timer for particle-wall interaction computations
+	Timer       pwInteractionTimer_;
+
+	/// timer for boundary interaction time 
+	Timer 		boundaryInteractionTimer_;
+
+	/// timer for managing contact lists (only inernal points)
+	Timer 		contactListMangementTimer_;	
+
+	Timer 		contactListMangementBoundaryTimer_;
+
+	
+
+	bool createGrainInteraction();
+
+	bool grainGrainInteraction();
+
+	bool grainWallInteraction();
+
+	/// range policy for p-p interaction execution 
+	using rpPPInteraction = 
+		Kokkos::RangePolicy<Kokkos::IndexType<uint32>, Kokkos::Schedule<Kokkos::Dynamic>>;
+
+	/// range policy for p-w interaction execution 
+	using rpPWInteraction = rpPPInteraction;
+
+public:
+
+	TypeInfoTemplate13("grainInteraction", ContactForceModel, MotionModel, ContactListType);
+
+	/// Constructor from components
+	grainInteraction(
+		systemControl& control,
+		const particles& prtcl,
+		const geometry& geom);
+	
+
+	/// Add virtual constructor 
+	add_vCtor
+	(
+		interaction,
+		grainInteraction, 
+		systemControl
+	);
+
+	/// This is called in time loop, before iterate. (overriden from demComponent)
+	bool beforeIteration() override;
+	
+	/// This is called in time loop. Perform the main calculations 
+	/// when the component should evolve along time. (overriden from demComponent)
+	bool iterate() override; 
+	
+	/// This is called in time loop, after iterate. (overriden from demComponent)
+	bool afterIteration() override;
+
+	/// Check for changes in the point structures. (overriden from observer)
+	bool hearChanges(
+		const timeInfo& ti,
+		const message& msg, 
+		const anyList& varList)override;
+	
+	
+};
+
+
+}
+
+#include "grainInteraction.cpp"
+
+#endif //__grainInteraction_hpp__

src/Interaction/grainInteraction/grainInteraction/grainInteractionKernels.hpp

@@ -0,0 +1,337 @@
+/*------------------------------- 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 __grainInteractionKernels_hpp__
+#define __grainInteractionKernels_hpp__
+
+
+#include "grainTriSurfaceContact.hpp"
+
+namespace pFlow::grainInteractionKernels
+{
+
+template<
+	typename ContactForceModel, 
+	typename ContactListType>
+struct ppInteractionFunctor
+{
+
+	using PairType = typename ContactListType::PairType;
+	using ValueType = typename ContactListType::ValueType;
+
+	real dt_;
+
+	ContactForceModel forceModel_;
+	ContactListType   tobeFilled_;
+
+	deviceViewType1D<real>  		diam_;
+	deviceViewType1D<real>  		coarseGrainFactor_;
+	deviceViewType1D<uint32> 		propId_;
+	deviceViewType1D<realx3>  		pos_;
+	deviceViewType1D<realx3>  		lVel_;
+	deviceViewType1D<realx3>  		rVel_;
+	deviceViewType1D<realx3>	cForce_;
+	deviceViewType1D<realx3>	cTorque_;
+
+
+	ppInteractionFunctor(
+		real dt,
+		ContactForceModel 		forceModel,
+		ContactListType 		tobeFilled,
+		deviceViewType1D<real>  		diam,
+		deviceViewType1D<real>  coarseGrainFactor,
+		deviceViewType1D<uint32> 		propId,
+		deviceViewType1D<realx3>  		pos,
+		deviceViewType1D<realx3>  		lVel,
+		deviceViewType1D<realx3>  		rVel,
+		deviceViewType1D<realx3>		cForce,
+		deviceViewType1D<realx3> 		cTorque )
+	:
+		dt_(dt),
+		forceModel_(forceModel),
+		tobeFilled_(tobeFilled),
+		diam_(diam),
+		coarseGrainFactor_(coarseGrainFactor),
+		propId_(propId),
+		pos_(pos),
+		lVel_(lVel),
+		rVel_(rVel),
+		cForce_(cForce), // this is converted to an atomic vector 
+		cTorque_(cTorque) // this is converted to an atomic vector 
+	{}
+
+	INLINE_FUNCTION_HD
+	void operator()(const uint32 n)const
+	{
+		
+		if(!tobeFilled_.isValid(n))return;
+
+		auto [i,j] = tobeFilled_.getPair(n);
+		
+		real Ri = 0.5*diam_[i];
+		real Rj = 0.5*diam_[j];
+
+		real cGFi = coarseGrainFactor_[j];
+		real cGFj = coarseGrainFactor_[j];
+
+		realx3 xi = pos_[i];
+		realx3 xj = pos_[j];
+		real dist = length(xj-xi);
+		real ovrlp = (Ri+Rj) - dist;
+		
+		if( ovrlp >0.0 )
+		{
+			
+			auto Vi = lVel_[i];
+			auto Vj = lVel_[j];
+			auto wi = rVel_[i];
+			auto wj = rVel_[j];
+			auto Nij = (xj-xi)/dist;
+			auto Vr = Vi - Vj + cross((Ri*wi+Rj*wj), Nij);
+			
+			auto history = tobeFilled_.getValue(n);
+
+			int32 propId_i = propId_[i];
+			int32 propId_j = propId_[j];
+
+			realx3 FCn, FCt, Mri, Mrj, Mij, Mji;
+
+			// calculates contact force 
+			forceModel_.contactForce(
+				dt_, i, j,
+				propId_i, propId_j,
+				Ri, Rj, cGFi , cGFj ,
+				ovrlp,
+				Vr, Nij,
+				history,
+				FCn, FCt
+				);
+
+			forceModel_.rollingFriction(
+				dt_, i, j,
+				propId_i, propId_j,
+				Ri, Rj, cGFi , cGFj , 
+				wi, wj,
+				Nij,
+				FCn,
+				Mri, Mrj
+				);
+
+			auto M = cross(Nij,FCt);
+			Mij = Ri*M+Mri;
+			Mji = Rj*M+Mrj;
+			
+			auto FC = FCn + FCt;
+			
+			Kokkos::atomic_add(&cForce_[i].x_,FC.x_);
+			Kokkos::atomic_add(&cForce_[i].y_,FC.y_);
+			Kokkos::atomic_add(&cForce_[i].z_,FC.z_);
+
+			Kokkos::atomic_add(&cForce_[j].x_,-FC.x_);
+			Kokkos::atomic_add(&cForce_[j].y_,-FC.y_);
+			Kokkos::atomic_add(&cForce_[j].z_,-FC.z_);
+
+			Kokkos::atomic_add(&cTorque_[i].x_, Mij.x_);
+			Kokkos::atomic_add(&cTorque_[i].y_, Mij.y_);
+			Kokkos::atomic_add(&cTorque_[i].z_, Mij.z_);
+
+			Kokkos::atomic_add(&cTorque_[j].x_, Mji.x_);
+			Kokkos::atomic_add(&cTorque_[j].y_, Mji.y_);
+			Kokkos::atomic_add(&cTorque_[j].z_, Mji.z_);
+			
+
+			tobeFilled_.setValue(n,history);
+
+		}
+		else
+		{
+			tobeFilled_.setValue(n, ValueType());
+		}
+
+	}
+};
+
+
+template<
+	typename ContactForceModel, 
+	typename ContactListType,
+	typename TraingleAccessor,
+	typename MotionModel>
+struct pwInteractionFunctor
+{
+	using PairType = typename ContactListType::PairType;
+	using ValueType = typename ContactListType::ValueType;
+
+	real dt_;
+
+	ContactForceModel 	forceModel_;
+	ContactListType   	tobeFilled_;
+
+	TraingleAccessor 		triangles_;
+	MotionModel  			motionModel_;
+
+	deviceViewType1D<real>  	diam_;
+	deviceViewType1D<real>  	coarseGrainFactor_;
+	deviceViewType1D<uint32>	propId_;
+	deviceViewType1D<realx3>  	pos_;
+	deviceViewType1D<realx3>  	lVel_;
+	deviceViewType1D<realx3>  	rVel_;
+	deviceViewType1D<realx3>	cForce_;
+	deviceViewType1D<realx3>	cTorque_;
+	deviceViewType1D<uint32> 	wTriMotionIndex_;
+	deviceViewType1D<uint32>    wPropId_;
+	deviceViewType1D<realx3>    wCForce_;
+	
+
+	pwInteractionFunctor(
+		real dt,
+		ContactForceModel 		forceModel,
+		ContactListType   		tobeFilled,
+		TraingleAccessor		triangles,
+		MotionModel 			motionModel ,
+		deviceViewType1D<real>  	diam ,
+		deviceViewType1D<real>  	coarseGrainFactor,
+		deviceViewType1D<uint32>  	propId,
+		deviceViewType1D<realx3>   	pos ,
+		deviceViewType1D<realx3>   	lVel,
+		deviceViewType1D<realx3>   	rVel,
+		deviceViewType1D<realx3> 	cForce,
+		deviceViewType1D<realx3> 	cTorque ,
+		deviceViewType1D<uint32> 	wTriMotionIndex,
+		deviceViewType1D<uint32>    wPropId,
+		deviceViewType1D<realx3>    wCForce)
+	:
+		dt_(dt),
+		forceModel_(forceModel),
+		tobeFilled_(tobeFilled),
+		triangles_(triangles) ,
+		motionModel_(motionModel) ,
+		diam_(diam) ,
+		coarseGrainFactor_(coarseGrainFactor),
+		propId_(propId),
+		pos_(pos) ,
+		lVel_(lVel),
+		rVel_(rVel) ,
+		cForce_(cForce),
+		cTorque_(cTorque) ,
+		wTriMotionIndex_(wTriMotionIndex) ,
+		wPropId_(wPropId),
+		wCForce_(wCForce)
+	{}
+
+	INLINE_FUNCTION_HD
+	void operator()(const int32 n)const
+	{
+
+		if(!tobeFilled_.isValid(n))return;
+
+		auto [i,tj] = tobeFilled_.getPair(n);
+		
+		real Ri = 0.5*diam_[i];
+		real Rj = 10000.0;
+		real cGFi = coarseGrainFactor_[i];
+		real cGFj = coarseGrainFactor_[i];
+		realx3 xi = pos_[i];
+		
+		realx3x3 tri = triangles_(tj); 
+		real ovrlp;
+		realx3 Nij, cp;
+
+		if(	pFlow::grnTriInteraction::isGrainInContactBothSides(
+			tri, xi, Ri, ovrlp, Nij, cp) )
+		{
+			auto  Vi = lVel_[i];
+			auto wi = rVel_[i];
+
+			int32 mInd	= wTriMotionIndex_[tj];
+			
+			auto Vw = motionModel_(mInd, cp);
+			
+			//output<< "par-wall index "<< i<<" - "<< tj<<endl;
+			
+			auto Vr = Vi - Vw + cross(Ri*wi, Nij);
+			
+			auto history = tobeFilled_.getValue(n);
+
+			int32 propId_i = propId_[i];
+			int32 wPropId_j = wPropId_[tj];
+
+			realx3 FCn, FCt, Mri, Mrj, Mij;
+			//output<< "before "<<history.overlap_t_<<endl;
+			// calculates contact force 
+			forceModel_.contactForce(
+				dt_, i, tj,
+				propId_i, wPropId_j,
+				Ri, Rj, cGFi , cGFj ,
+				ovrlp,
+				Vr, Nij,
+				history,
+				FCn, FCt
+				);
+
+			//output<< "after "<<history.overlap_t_<<endl;
+
+			forceModel_.rollingFriction(
+				dt_, i, tj,
+				propId_i, wPropId_j,
+				Ri, Rj, cGFi , cGFj ,
+				wi, 0.0,
+				Nij,
+				FCn,
+				Mri, Mrj
+				);
+
+			auto M = cross(Nij,FCt);
+			Mij = Ri*M+Mri;
+			//output<< "FCt = "<<FCt <<endl;
+			//output<< "FCn = "<<FCn <<endl;
+			//output<<"propId i, tj"<< propId_i << " "<<wPropId_j<<endl;
+			//output<<"par i , wj"<< i<<" "<< tj<<endl;
+
+			auto FC = FCn + FCt;
+
+			Kokkos::atomic_add(&cForce_[i].x_,FC.x_);
+			Kokkos::atomic_add(&cForce_[i].y_,FC.y_);
+			Kokkos::atomic_add(&cForce_[i].z_,FC.z_);
+
+			Kokkos::atomic_add(&wCForce_[tj].x_,-FC.x_);
+			Kokkos::atomic_add(&wCForce_[tj].y_,-FC.y_);
+			Kokkos::atomic_add(&wCForce_[tj].z_,-FC.z_);
+
+			
+			Kokkos::atomic_add(&cTorque_[i].x_, Mij.x_);
+			Kokkos::atomic_add(&cTorque_[i].y_, Mij.y_);
+			Kokkos::atomic_add(&cTorque_[i].z_, Mij.z_);
+
+			tobeFilled_.setValue(n,history);
+
+		}
+		else
+		{
+			tobeFilled_.setValue(n,ValueType());
+		}
+		
+	}
+
+};
+
+} 
+
+#endif //__grainInteractionKernels_hpp__

src/Interaction/grainInteraction/grainInteraction/grainTriSurfaceContact.hpp

@@ -0,0 +1,231 @@
+/*------------------------------- 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 __grainTriSurfaceContact_hpp__
+#define __grainTriSurfaceContact_hpp__
+
+#include "triWall.hpp"
+#include "pLine.hpp"
+
+namespace pFlow::grnTriInteraction
+{
+
+INLINE_FUNCTION_HD
+bool pointInPlane(
+	const realx3& p1,
+	const realx3& p2,
+	const realx3& p3,
+	const realx3& p )
+{
+
+	realx3 p1p = p1 - p;
+	realx3 p2p = p2 - p;
+	realx3 p3p = p3 - p;
+
+	real p1p2 = dot(p1p, p2p);
+	real p2p3 = dot(p2p, p3p);
+	real p2p2 = dot(p2p, p2p);
+	real p1p3 = dot(p1p, p3p);
+
+	// first condition u.v < 0
+	// u.v = [(p1-p)x(p2-p)].[(p2-p)x(p3-p)] = (p1p.p2p)(p2p.p3p) - (p2p.p2p)(p1p.p3p) 
+	if (p1p2*p2p3 - p2p2*p1p3 < 0.0) return false;
+
+
+	// second condition u.w < 0
+	// u.w = [(p1-p)x(p2-p)].[(p3-p)x(p1-p)] = (p1p.p3p)(p2p.p1p) - (p2p.p3p)(p1p.p1p)
+	real p1p1 = dot(p1p, p1p);
+	if (p1p3*p1p2 - p2p3*p1p1 < (0.0)) return false;
+
+	return true;
+}
+
+INLINE_FUNCTION_HD
+void cramerRule2(real A[2][2], real B[2], real & x1, real &x2)
+{
+	real det = (A[0][0] * A[1][1] - A[1][0]*A[0][1]);
+	x1 = (B[0]*A[1][1] - B[1]*A[0][1]) / det;
+	x2 = (A[0][0] * B[1] - A[1][0] * B[0])/ det;
+}
+
+INLINE_FUNCTION_HD
+bool pointInPlane(
+	const realx3& p1,
+	const realx3& p2,
+	const realx3& p3,
+	const realx3 &p,
+	int32 &Ln)
+{
+	realx3 v0 = p2 - p1;
+	realx3 v1 = p3 - p1;
+	realx3 v2 = p - p1;
+
+	real A[2][2] = { dot(v0, v0), dot(v0, v1), dot(v1, v0), dot(v1, v1) };
+	real B[2] = { dot(v0, v2), dot(v1, v2) };
+	real nu, w;
+
+	cramerRule2(A, B, nu, w);
+	real nuW = nu + w;
+		
+	
+	if (nuW > 1) 
+	{
+		Ln = 2;	return true;
+	}
+	else if (nuW >= 0)
+	{
+		if (nu >= 0 && w >= 0)
+		{
+			Ln = 0;	return true;
+		}
+		if (nu > 0 && w < 0)
+		{
+			Ln = 1;	return true;
+		}
+		if (nu < 0 && w > 0)
+		{
+			Ln = 3;	return true;
+		}
+	}
+	else
+	{
+		Ln = 1;	return true;
+	}
+	
+	return false;
+}
+
+INLINE_FUNCTION_HD
+bool isGrainInContactActiveSide(
+	const realx3& p1,
+	const realx3& p2,
+	const realx3& p3,
+	const realx3& cntr,
+	real 	rad,
+	real& 	ovrlp,
+	realx3& norm,
+	realx3& cp)
+{
+	
+	triWall wall(true, p1,p2,p3);
+
+	real dist = wall.normalDistFromWall(cntr);
+	
+	if(dist < 0.0 )return false;
+
+	ovrlp = rad - dist;
+
+	if (ovrlp > 0)
+	{
+		realx3 ptOnPlane = wall.nearestPointOnWall(cntr);
+		
+		if (pointInPlane(p1, p2, p3, ptOnPlane))
+		{
+			cp = ptOnPlane;
+			norm = -wall.n_;
+			return true;
+		}
+
+		realx3 lnv;
+		
+		if (pLine(p1,p2).lineGrainCheck(cntr, rad, lnv, cp, ovrlp))
+		{
+			norm = -lnv;
+			return true;
+		}
+
+		if ( pLine(p2,p3).lineGrainCheck(cntr, rad, lnv, cp, ovrlp))
+		{
+			norm = -lnv;
+			return true;
+		}
+
+		if ( pLine(p3,p1).lineGrainCheck(cntr, rad, lnv, cp, ovrlp))
+		{
+			norm = -lnv;
+			return true;
+		}
+	}
+			
+	return false;
+}
+
+INLINE_FUNCTION_HD
+bool isGrainInContactBothSides(
+	const realx3x3& tri,
+	const realx3& cntr,
+	real 	Rad,
+	real& 	ovrlp,
+	realx3& norm,
+	realx3& cp)
+{
+	
+	triWall wall(true, tri.x_,tri.y_,tri.z_);
+
+	real dist = wall.normalDistFromWall(cntr);
+	
+	
+	ovrlp = Rad - abs(dist);
+
+	if (ovrlp > 0)
+	{
+		realx3 ptOnPlane = wall.nearestPointOnWall(cntr);
+		
+		if (pointInPlane(tri.x_,tri.y_,tri.z_,ptOnPlane))
+		{
+			cp = ptOnPlane;
+			
+			if(dist >= 0.0)
+				norm = -wall.n_;
+			else
+				norm = wall.n_;
+			return true;
+		}
+
+		realx3 lnv;
+		
+		if (pLine(tri.x_, tri.y_).lineGrainCheck(cntr, Rad, lnv, cp, ovrlp))
+		{
+			norm = -lnv;
+			return true;
+		}
+
+		if ( pLine(tri.y_, tri.z_).lineGrainCheck(cntr, Rad, lnv, cp, ovrlp))
+		{
+			norm = -lnv;
+			return true;
+		}
+
+		if ( pLine(tri.z_, tri.x_).lineGrainCheck(cntr, Rad, lnv, cp, ovrlp))
+		{
+			norm = -lnv;
+			return true;
+		}
+	}
+			
+	return false;
+}
+
+
+} // pFlow::grnTriInteraction
+
+
+#endif //__grainTriSurfaceContact_hpp__

src/Interaction/grainInteraction/grainInteraction/pLine.hpp

@@ -0,0 +1,107 @@
+/*------------------------------- 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 __pLine_hpp__
+#define __pLine_hpp__
+
+#include "types.hpp"
+
+namespace pFlow::grnTriInteraction
+{
+
+struct pLine
+{
+
+	realx3 p1_; // point 1
+	realx3 p2_; // piont 2
+	realx3 v_;  // direction vector 
+	real  L_;  // line lenght
+
+	INLINE_FUNCTION_HD
+	pLine(){}
+
+	INLINE_FUNCTION_HD 
+	pLine(const realx3 &p1, const realx3 &p2)
+	:
+		p1_(p1),
+		p2_(p2),
+		v_(p2-p1),
+		L_(length(v_))
+	{}
+	 	
+	// get a point on the line based on input  0<= t <= 1
+	INLINE_FUNCTION_HD 
+	realx3 point(real t)const {
+		return v_ * t + p1_;
+		}
+
+	// return the projected point of point p on line
+	INLINE_FUNCTION_HD
+	realx3 projectPoint(const realx3 &p) const
+	{
+		return point(projectNormLength(p));
+	}
+
+	// calculates the normalized distance between projected p and p1 
+	INLINE_FUNCTION_HD
+	real projectNormLength(realx3 p) const 
+	{
+		realx3 w = p - p1_;
+		return dot(w,v_) / (L_*L_);
+	}
+
+	INLINE_FUNCTION_HD
+	bool lineGrainCheck(
+		const realx3 pos,
+		real Rad,
+		realx3 &nv,
+		realx3 &cp,
+		real &ovrlp)const
+	{
+
+		
+		real t = projectNormLength(pos);
+
+		if(t >= 0.0 && t <= 1.0) 			 	 cp = point(t);
+		else if(t >= (-Rad / L_) && t < 0.0) 	 cp = point(0.0);
+		else if(t>1.0 && t >= (1.0 + Rad / L_)) cp = point(1.0);
+		else return false;
+		
+		realx3 vec = pos - cp; // from cp to pos
+
+		real dist = length(vec);
+		ovrlp = Rad - dist;
+
+		if (ovrlp >= 0.0)
+		{
+			if (dist > 0)
+				nv = vec / dist;
+			else
+				nv = v_;
+			return true;
+		}
+
+		return false;
+	}
+};
+
+} //pFlow::grnTriInteractio
+
+#endif 

src/Interaction/grainInteraction/grainInteraction/triWall.hpp

@@ -0,0 +1,108 @@
+/*------------------------------- 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 __triWall_hpp__
+#define __triWall_hpp__
+
+#include "types.hpp"
+
+namespace pFlow::grnTriInteraction
+{
+
+struct triWall
+{
+	
+	realx3 n_;
+	real   offset_;
+
+	INLINE_FUNCTION_H
+	triWall(const realx3& p1, const realx3& p2, const realx3& p3)
+	{
+		if(!makeWall(p1,p2,p3, n_, offset_))
+		{
+			fatalErrorInFunction<<
+			"bad input for the wall.\n";
+			fatalExit;
+		}
+	}
+
+	INLINE_FUNCTION_HD
+	triWall(bool, const realx3& p1, const realx3& p2, const realx3& p3)
+	{
+		makeWall(p1,p2,p3,n_,offset_);
+	}
+
+	INLINE_FUNCTION_HD
+	triWall(const realx3x3& tri)
+	{
+		makeWall(tri.x_, tri.y_, tri.z_, n_, offset_);
+	}
+
+	INLINE_FUNCTION_HD
+	triWall(const triWall&) = default;
+
+	INLINE_FUNCTION_HD
+	triWall& operator=(const triWall&) = default;
+
+	INLINE_FUNCTION_HD
+	triWall(triWall&&) = default;
+
+	INLINE_FUNCTION_HD
+	triWall& operator=(triWall&&) = default;
+
+	INLINE_FUNCTION_HD
+	~triWall()=default;
+
+	
+	INLINE_FUNCTION_HD 
+	real normalDistFromWall(const realx3 &p) const
+	{
+		return dot(n_, p) + offset_;
+	}
+
+	INLINE_FUNCTION_HD
+	realx3 nearestPointOnWall(const realx3 &p) const
+	{
+		real t = -(dot(n_, p) + offset_);
+		return realx3(n_.x_*t + p.x_, n_.y_*t + p.y_, n_.z_*t + p.z_);
+	}
+
+	INLINE_FUNCTION_HD static 
+	bool makeWall(
+		const realx3& p1,
+		const realx3& p2,
+		const realx3& p3,
+		realx3& n, real& offset)
+	{
+		n = cross(p2 - p1, p3 - p1);
+		real len = length(n);
+	
+		if (len < 0.00000000001) return false;
+		n /= len;
+		offset = -dot(n, p1);
+		return true;
+	}
+
+
+}; 
+
+}
+
+#endif

src/Interaction/grainInteraction/grainInteractionsLinearModels.cpp

@@ -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.
+
+-----------------------------------------------------------------------------*/
+
+#include "grainInteraction.hpp"
+#include "geometryMotions.hpp"
+#include "grainContactForceModels.hpp"
+#include "unsortedContactList.hpp"
+#include "sortedContactList.hpp"
+#include "createBoundaryGrainInteraction.hpp"
+
+
+
+#define createInteraction(ForceModel,GeomModel) 	\
+													\
+	template class pFlow::grainInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::unsortedContactList>;				\
+													\
+	template class pFlow::grainInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::sortedContactList>;					\
+	createBoundaryGrainInteraction(ForceModel, GeomModel)
+
+
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGAbsoluteLinear<true>>, pFlow::rotationAxisMotionGeometry);
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinear<true>>, pFlow::rotationAxisMotionGeometry);
+
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGAbsoluteLinear<false>>, pFlow::rotationAxisMotionGeometry);
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinear<false>>, pFlow::rotationAxisMotionGeometry);
+
+
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGAbsoluteLinear<true>>, pFlow::stationaryGeometry);
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinear<true>>, pFlow::stationaryGeometry);
+
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGAbsoluteLinear<false>>, pFlow::stationaryGeometry);
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinear<false>>, pFlow::stationaryGeometry);
+
+
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGAbsoluteLinear<true>>, pFlow::vibratingMotionGeometry);
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinear<true>>, pFlow::vibratingMotionGeometry);
+
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGAbsoluteLinear<false>>, pFlow::vibratingMotionGeometry);
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinear<false>>, pFlow::vibratingMotionGeometry);
+
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGAbsoluteLinear<true>>, pFlow::conveyorBeltMotionGeometry);
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinear<true>>, pFlow::conveyorBeltMotionGeometry);
+
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGAbsoluteLinear<false>>, pFlow::conveyorBeltMotionGeometry);
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinear<false>>, pFlow::conveyorBeltMotionGeometry);

src/Interaction/grainInteraction/grainInteractionsNonLinearModModels.cpp

@@ -18,39 +18,26 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-#include "eventObserver.hpp"
-#include "eventSubscriber.hpp" 
+#include "grainInteraction.hpp"
+#include "geometryMotions.hpp"
+#include "grainContactForceModels.hpp"
+#include "unsortedContactList.hpp"
+#include "sortedContactList.hpp"
+#include "createBoundaryGrainInteraction.hpp"
 
-pFlow::eventObserver::eventObserver():
-	subscriber_(nullptr),
-	subscribed_(false)
-{}
 
-pFlow::eventObserver::eventObserver
-(
-	const eventSubscriber& subscriber,
-	bool subscribe
-)
-:
-	subscriber_(&subscriber),
-	subscribed_(false)
-{
-	if(subscribe && subscriber_)
-	{
-		subscribed_ = subscriber_->subscribe(this);
-	}
-}
 
-pFlow::eventObserver::~eventObserver()
-{
-	if(subscribed_ && subscriber_)
-		subscriber_->unsubscribe(this);
-}
+#define createInteraction(ForceModel,GeomModel) 	\
+													\
+	template class pFlow::grainInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::unsortedContactList>;				\
+													\
+	template class pFlow::grainInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::sortedContactList>;					\
+	createBoundaryGrainInteraction(ForceModel, GeomModel)
 
-bool pFlow::eventObserver::subscribe(const eventSubscriber& subscriber)
-{
-	subscriber_ = &subscriber;
-	subscribed_ = subscriber_->subscribe(this);
-	return subscribed_;
-}
 

src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp

@@ -0,0 +1,49 @@
+/*------------------------------- 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 "grainInteraction.hpp"
+#include "geometryMotions.hpp"
+#include "grainContactForceModels.hpp"
+#include "unsortedContactList.hpp"
+#include "sortedContactList.hpp"
+#include "createBoundaryGrainInteraction.hpp"
+
+
+
+#define createInteraction(ForceModel,GeomModel) 	\
+													\
+	template class pFlow::grainInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::unsortedContactList>;				\
+													\
+	template class pFlow::grainInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::sortedContactList>;					\
+	createBoundaryGrainInteraction(ForceModel, GeomModel)
+
+
+createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry);
+
+createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry);
+
+createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::stationaryGeometry);
+createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::stationaryGeometry);

src/Interaction/interaction/interaction.hpp

@@ -51,7 +51,7 @@ private:
 	const geometry& 	geometry_;
 
 	static inline
-	const message msg_ = message::ITEM_REARRANGE;
+	const message msg_ = message::ITEMS_REARRANGE;
 
 public:
 

src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.hpp

@@ -110,7 +110,7 @@ public:
 		return geometryMotion_;
 	}
 
-	ContactListType& ppPairs()
+	ContactListType& ppPairs() 
 	{
 		return ppPairs_();
 	}
@@ -148,7 +148,7 @@ public:
 		const ContactForceModel& cfModel,
 		uint32 step)
 	{
-		// for default boundary, no thing to be done 
+		// for default boundary, nothing to be done 
 		return false;
 	}
 
@@ -156,20 +156,26 @@ public:
 
 	bool hearChanges
 	(
-		real t,
-		real dt,
-		uint32 iter,
+		const timeInfo& ti,
 		const message& msg, 
     	const anyList& varList
 	) override
     {
-		
+		if(msg.equivalentTo(message::BNDR_RESET))
+		{
+			// do nothing 
+			return true;
+		}
 		pOutput<<"Function (hearChanges in boundarySphereInteractions)is not implmented Message "<<
-		 msg <<endl<<" name "<< this->name() <<" type "<< this->type()<<endl;;
-		//notImplementedFunction;
+		 msg <<endl<<" name "<< this->boundaryName() <<" type "<< this->type()<<endl;
 		return true;
 	}
 
+	bool isActive()const override
+	{
+		return false;
+	}
+
 	static
 	uniquePtr<BoundarySphereInteractionType> create(
 		const boundaryBase& boundary,

src/Interaction/sphereInteraction/boundaries/boundarySphereInteractionList.cpp

@@ -7,11 +7,11 @@ pFlow::boundarySphereInteractionList<CFModel, gMModel>::boundarySphereInteractio
     const gMModel &geomMotion
 )
 :
-	ListPtr<boundarySphereInteraction<CFModel,gMModel>>(6),
+	boundaryListPtr<boundarySphereInteraction<CFModel,gMModel>>(),
 	boundaries_(sphPrtcls.pStruct().boundaries())
 {
-	//gSettings::sleepMiliSeconds(1000*pFlowProcessors().localRank());
-	for(uint32 i=0; i<6; i++)
+	output<<boundaries_.size()<<endl;
+	ForAllBoundariesPtr(i, this)
 	{
 		this->set(
 			i,

src/Interaction/sphereInteraction/boundaries/boundarySphereInteractionList.hpp

@@ -3,7 +3,7 @@
 
 
 #include "boundaryList.hpp"
-#include "ListPtr.hpp"
+#include "boundaryListPtr.hpp"
 #include "boundarySphereInteraction.hpp"
 
 
@@ -14,7 +14,7 @@ namespace pFlow
 template<typename contactForceModel,typename geometryMotionModel>
 class boundarySphereInteractionList
 :
-    public ListPtr<boundarySphereInteraction<contactForceModel,geometryMotionModel>>
+    public boundaryListPtr<boundarySphereInteraction<contactForceModel,geometryMotionModel>>
 {
 private:
 

src/Interaction/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.hpp

@@ -76,15 +76,17 @@ public:
         boundaryBase
     );
     
-    ~periodicBoundarySphereInteraction()override = default;
-
-
-    
+    ~periodicBoundarySphereInteraction()override = default;  
 
     bool sphereSphereInteraction(
         real dt,
 		const ContactForceModel& cfModel,
         uint32 step)override;
+    
+    bool isActive()const override
+	{
+		return true;
+	}
 	
 };
 

src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp

@@ -41,9 +41,9 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::createSphereInteraction()
 	geometryMotion_,
 	timers());
 
-	ppContactList_ = makeUnique<ContactListType>(nPrtcl+1);
+	ppContactList_ = makeUnique<ContactListType>("sphere-sphere",nPrtcl+1);
 	
-	pwContactList_ = makeUnique<ContactListType>(nPrtcl/5+1);
+	pwContactList_ = makeUnique<ContactListType>("sphere-wall",nPrtcl/5+1);
 
 	return true;
 }
@@ -134,18 +134,22 @@ pFlow::sphereInteraction<cFM,gMM, cLT>::sphereInteraction
 	geometryMotion_(dynamic_cast<const GeometryMotionModel&>(geom)),
 	sphParticles_(dynamic_cast<const sphereParticles&>(prtcl)),
 	boundaryInteraction_(sphParticles_,geometryMotion_),
-	ppInteractionTimer_("sphere-sphere interaction", &this->timers()),
-	pwInteractionTimer_("sphere-wall interaction", &this->timers()),
-	contactListMangementTimer_("contact-list management", &this->timers()),
-	boundaryContactSearchTimer_("contact search for boundary", &this->timers()),
-	boundaryInteractionTimer_("interaction for boundary", &this->timers()),
-	contactListBoundaryTimer_("contact-list management for boundary", &this->timers())
+	ppInteractionTimer_("sphere-sphere interaction (internal)", &this->timers()),
+	pwInteractionTimer_("sphere-wall interaction (internal)", &this->timers()),
+	boundaryInteractionTimer_("sphere-sphere interaction (boundary)",&this->timers()),
+	contactListMangementTimer_("list management (interal)", &this->timers()),
+	contactListMangementBoundaryTimer_("list management (boundary)", &this->timers())
 {
 	
 	if(!createSphereInteraction())
 	{
 		fatalExit;
 	}
+
+	for(uint32 i=0; i<6; i++)
+	{
+		activeBoundaries_[i] = boundaryInteraction_[i].ppPairsAllocated();
+	}
 }
 
 template<typename cFM,typename gMM,template <class, class, class> class cLT>
@@ -158,45 +162,58 @@ template<typename cFM,typename gMM,template <class, class, class> class cLT>
 bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 {
 	
-	auto iter = this->currentIter();
-	auto t = this->currentTime();
-	auto dt = this->dt();
+	timeInfo ti = this->TimeInfo();
+	auto iter = ti.iter();
+	auto t = ti.t();
+	auto dt = ti.dt();
 
-	
-	bool broadSearch = contactSearch_().enterBroadSearch(iter, t, dt);
+	auto& contactSearchRef = contactSearch_();
 
+	bool broadSearch = contactSearchRef.enterBroadSearch(ti);
+	bool broadSearchBoundary = contactSearchRef.enterBroadSearchBoundary(ti);
+
+	/// update boundaries of the fields that are being used by inreaction 
 	sphParticles_.diameter().updateBoundaries(DataDirection::SlaveToMaster);
 	sphParticles_.velocity().updateBoundaries(DataDirection::SlaveToMaster);
 	sphParticles_.rVelocity().updateBoundaries(DataDirection::SlaveToMaster);
-	sphParticles_.mass().updateBoundaries(DataDirection::SlaveToMaster);
-	sphParticles_.I().updateBoundaries(DataDirection::SlaveToMaster);
 	sphParticles_.propertyId().updateBoundaries(DataDirection::SlaveToMaster);
 	
-	
+	/// lists 
 	if(broadSearch)
 	{
 		contactListMangementTimer_.start();
-		ppContactList_().beforeBroadSearch();
-		pwContactList_().beforeBroadSearch();
+		ComputationTimer().start();
+			ppContactList_().beforeBroadSearch();
+			pwContactList_().beforeBroadSearch();
+		ComputationTimer().end();
 		contactListMangementTimer_.pause();
 	} 
 
-	contactListBoundaryTimer_.start();
-	for(uint32 i=0; i<6u; i++)
+	if(broadSearchBoundary)
 	{
-		auto& BI = boundaryInteraction_[i];
-		if(BI.ppPairsAllocated()) BI.ppPairs().beforeBroadSearch();
-		if(BI.pwPairsAllocated()) BI.pwPairs().beforeBroadSearch();
+		contactListMangementBoundaryTimer_.start();
+		ComputationTimer().start();
+
+		ForAllActiveBoundaries(i, boundaryInteraction_)
+		//for(size_t i=0; i<6; i++)
+		{
+			if(activeBoundaries_[i])
+			{
+				auto& BI = boundaryInteraction_[i];
+				BI.ppPairs().beforeBroadSearch();
+				BI.pwPairs().beforeBroadSearch();
+			}
+		}
+
+		ComputationTimer().end();
+		contactListMangementBoundaryTimer_.pause();	
 	}
-	contactListBoundaryTimer_.pause();
 	
 	if( sphParticles_.numActive()<=0)return true;
 	
-	
-	if( !contactSearch_().broadSearch(
-			iter,
-			t,
-			dt,		
+	ComputationTimer().start();
+	if( !contactSearchRef.broadSearch(
+			ti,		
 			ppContactList_(),
 			pwContactList_()) )
 	{
@@ -205,62 +222,69 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 		fatalExit;
 	}
 	
-	boundaryContactSearchTimer_.start();
-	for(uint32 i=0; i<6u; i++)
+	ForAllActiveBoundaries(i, boundaryInteraction_)
+	//for(size_t i=0; i<6; i++)
 	{
-		auto& BI =boundaryInteraction_[i];
-		if(BI.ppPairsAllocated())
+		if(activeBoundaries_[i])
 		{
-			if( !contactSearch_().boundaryBroadSearch(
+			auto& BI = boundaryInteraction_[i];
+			if(!contactSearchRef.boundaryBroadSearch(
 				i,
-				iter,
-				t,
-				dt,
+				ti,
 				BI.ppPairs(),
-				BI.pwPairs()))
+				BI.pwPairs())
+			)
 			{
 				fatalErrorInFunction<<
 				"failed to perform broadSearch for boundary index "<<i<<endl;
 				return false;
 			}
-		} 
+		}
 	}
-	boundaryContactSearchTimer_.end();
+	ComputationTimer().end();
 	
-	if(broadSearch && contactSearch_().performedBroadSearch())
+	if(broadSearch && contactSearchRef.performedSearch())
 	{
 		contactListMangementTimer_.resume();
-		ppContactList_().afterBroadSearch();
-		pwContactList_().afterBroadSearch();
+		ComputationTimer().start();
+			ppContactList_().afterBroadSearch();
+			pwContactList_().afterBroadSearch();
+		ComputationTimer().end();
 		contactListMangementTimer_.end();
 	}
 
-	contactListBoundaryTimer_.resume();
-	for(uint32 i=0; i<6u; i++)
+	if(broadSearchBoundary && contactSearchRef.performedSearchBoundary())
 	{
-		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++)
+		contactListMangementBoundaryTimer_.resume();
+		ComputationTimer().start();
+		
+		ForAllActiveBoundaries(i, boundaryInteraction_)
+		//for(size_t i=0; i<6; i++)
 		{
-			if(step==1u || requireStep[i] )
+			if(activeBoundaries_[i])
+			{
+				auto& BI = boundaryInteraction_[i];
+				BI.ppPairs().afterBroadSearch();
+				BI.pwPairs().afterBroadSearch();
+			}
+		}
+
+		ComputationTimer().end();
+		contactListMangementBoundaryTimer_.end();
+	}
+	
+	/// peform contact search on boundareis with active contact search (master boundaries)
+	auto requireStep = boundariesMask<6>(true);
+	const auto& cfModel = this->forceModel_();
+	uint32 step=1;
+	boundaryInteractionTimer_.start();
+	ComputationTimer().start();
+	while(requireStep.anyElement(true) && step <= 10)
+	{
+		ForAllBoundaries(i, boundaryInteraction_)
+		//for(size_t i=0; i<6; i++)
+		{
+			if(requireStep[i] )
 			{
 				requireStep[i] = boundaryInteraction_[i].sphereSphereInteraction(
 					dt, 
@@ -271,44 +295,47 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 		}
 		step++;
 	}
+	ComputationTimer().end();
 	boundaryInteractionTimer_.pause();
-	}
+	
+
 	ppInteractionTimer_.start();
+	ComputationTimer().start();
 		sphereSphereInteraction();
+	ComputationTimer().end();
 	ppInteractionTimer_.end();
 
 	
 	pwInteractionTimer_.start();
+	ComputationTimer().start();
 		sphereWallInteraction();
+	ComputationTimer().start();
 	pwInteractionTimer_.end();
 
 	{
 	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()};
+	ComputationTimer().start();
+	auto requireStep = boundariesMask<6>(true);
 
-	int step = 2;
+	uint32 step = 11;
 	const auto& cfModel = this->forceModel_();
-	while(std::any_of(requireStep.begin(), requireStep.end(), [](bool v) { return v==true; }))
+	while( requireStep.anyElement(true) && step < 20 )
 	{
-		for(uint32 i=0; i<6u; i++)
+		ForAllBoundaries(i, boundaryInteraction_)
+		//for(size_t i=0; i<6; i++)
 		{
 			if(requireStep[i])
 			{
 				requireStep[i] = boundaryInteraction_[i].sphereSphereInteraction(
 					dt, 
-					this->forceModel_(),
+					cfModel,
 					step
 				);
 			}		
 		}
 		step++;
 	}
+	ComputationTimer().end();
 	boundaryInteractionTimer_.end();
 	}
 	
@@ -324,17 +351,18 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::afterIteration()
 template<typename cFM,typename gMM,template <class, class, class> class cLT>
 bool pFlow::sphereInteraction<cFM,gMM, cLT>::hearChanges
 (
-	real t,
-	real dt,
-	uint32 iter,
+	const timeInfo& ti,
 	const message& msg, 
 	const anyList& varList
 )
 {
-	if(msg.equivalentTo(message::ITEM_REARRANGE))
+	if(msg.equivalentTo(message::ITEMS_REARRANGE))
 	{
 		notImplementedFunction;
+		return false;
 	}
-	return true;
-}
-
+	
+	fatalErrorInFunction<<"Event "<< msg.eventNames()<<
+	" is not handled in sphereInteraction"<<endl;
+	return false;
+}

src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.hpp

@@ -25,6 +25,8 @@ Licence:
 #include "sphereParticles.hpp"
 #include "boundarySphereInteractionList.hpp"
 #include "sphereInteractionKernels.hpp"
+#include "boundariesMask.hpp"
+#include "MPITimer.hpp"
 //#include "unsortedContactList.hpp"
 
 
@@ -74,6 +76,9 @@ private:
 	/// particle-particle and particle-wall interactions at boundaries 
 	BoundaryListType 					boundaryInteraction_;	
 
+	/// a mask for active boundaries (boundaries with intreaction)
+	boundariesMask<6> 					activeBoundaries_;
+
 	/// contact search object for pp and pw interactions 
 	uniquePtr<contactSearch> 			contactSearch_ = nullptr;
 
@@ -93,14 +98,13 @@ private:
 	/// timer for particle-wall interaction computations
 	Timer       pwInteractionTimer_;
 
-	/// timer for managing contact lists (only inernal points)
-	Timer 		contactListMangementTimer_;
-
-	Timer 		boundaryContactSearchTimer_;
 	/// timer for boundary interaction time 
 	Timer 		boundaryInteractionTimer_;
 
-	Timer 		contactListBoundaryTimer_;
+	/// timer for managing contact lists (only inernal points)
+	Timer 		contactListMangementTimer_;	
+
+	Timer 		contactListMangementBoundaryTimer_;
 
 	
 
@@ -110,10 +114,6 @@ private:
 
 	bool sphereWallInteraction();
 
-	//bool managePPContactLists();
-
-	//bool managePWContactLists();
-
 	/// range policy for p-p interaction execution 
 	using rpPPInteraction = 
 		Kokkos::RangePolicy<Kokkos::IndexType<uint32>, Kokkos::Schedule<Kokkos::Dynamic>>;
@@ -152,9 +152,7 @@ public:
 
 	/// Check for changes in the point structures. (overriden from observer)
 	bool hearChanges(
-		real t,
-		real dt,
-		uint32 iter,
+		const timeInfo& ti,
 		const message& msg, 
 		const anyList& varList)override;
 	

src/Interaction/sphereInteraction/sphereInteractionsLinearModels.cpp

@@ -53,6 +53,10 @@ createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::rotation
 createInteraction(pFlow::cfModels::limitedLinearNormalRolling, pFlow::vibratingMotionGeometry);
 createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::vibratingMotionGeometry);
 
+// conveyorBeltGeometry
+createInteraction(pFlow::cfModels::limitedLinearNormalRolling, pFlow::conveyorBeltMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::conveyorBeltMotionGeometry);
+
 // multiRotationAxisMotionGeometry
 //createInteraction(pFlow::cfModels::limitedLinearNormalRolling, pFlow::multiRotationAxisMotionGeometry);
 //createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::multiRotationAxisMotionGeometry);

src/Interaction/sphereInteraction/sphereInteractionsNonLinearModModels.cpp

@@ -53,6 +53,10 @@ createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::ro
 createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::vibratingMotionGeometry);
 createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::vibratingMotionGeometry);
 
+// conveyorBeltMotionGeometry
+createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::conveyorBeltMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::conveyorBeltMotionGeometry);
+
 // multiRotationAxisMotionGeometry
-//createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::multiRotationAxisMotionGeometry);
-//createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::multiRotationAxisMotionGeometry);
+createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::multiRotationAxisMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::multiRotationAxisMotionGeometry);

src/Interaction/sphereInteraction/sphereInteractionsNonLinearModels.cpp

@@ -53,6 +53,10 @@ createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::rotat
 createInteraction(pFlow::cfModels::limitedNonLinearNormalRolling, pFlow::vibratingMotionGeometry);
 createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::vibratingMotionGeometry);
 
+// conveyorBeltMotionGeometry
+createInteraction(pFlow::cfModels::limitedNonLinearNormalRolling, pFlow::conveyorBeltMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::conveyorBeltMotionGeometry);
+
 // multiRotationAxisMotionGeometry
 //createInteraction(pFlow::cfModels::limitedNonLinearNormalRolling, pFlow::multiRotationAxisMotionGeometry);
 //createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::multiRotationAxisMotionGeometry);

src/MotionModel/CMakeLists.txt

@@ -11,8 +11,11 @@ vibratingMotion/vibratingMotion.cpp
 stationaryWall/stationaryWall.cpp
 entities/stationary/stationary.cpp
 
-#entities/multiRotatingAxis/multiRotatingAxis.cpp
-#multiRotatingAxisMotion/multiRotatingAxisMotion.cpp
+conveyorBeltMotion/conveyorBeltMotion.cpp
+entities/conveyorBelt/conveyorBelt.cpp
+
+entities/multiRotatingAxis/multiRotatingAxis.cpp
+multiRotatingAxisMotion/multiRotatingAxisMotion.cpp
 
 )
 

src/MotionModel/MotionModel/MotionModel.cpp

@@ -32,7 +32,6 @@ bool pFlow::MotionModel<Model, Component>::impl_nameToIndex(const word& name, ui
 		indx = static_cast<uint32>(i);
 		return true;
 	}
-	
 }
 
 template<typename Model, typename Component>

src/MotionModel/conveyorBeltMotion/conveyorBeltMotion.cpp

@@ -0,0 +1,75 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#include "conveyorBeltMotion.hpp"
+
+
+pFlow::conveyorBeltMotion::conveyorBeltMotion
+(
+	const objectFile &objf, 
+	repository *owner
+)
+:
+	fileDictionary(objf, owner)
+{
+
+	if(!impl_readDictionary(*this) )
+	{
+		fatalErrorInFunction;
+		fatalExit;
+	}
+}
+
+pFlow::conveyorBeltMotion::conveyorBeltMotion
+(
+	const objectFile &objf, 
+	const dictionary &dict, 
+	repository *owner
+)
+:
+	fileDictionary(objf, dict, owner)
+{
+	if(!impl_readDictionary(*this) )
+	{
+		fatalErrorInFunction;
+		fatalExit;
+	}
+}
+
+bool pFlow::conveyorBeltMotion::write
+(
+	iOstream &os, 
+	const IOPattern &iop
+) const
+{
+	// a global dictionary
+	dictionary newDict(fileDictionary::dictionary::name(), true);
+	if( iop.thisProcWriteData() )
+	{
+		if( !this->impl_writeDictionary(newDict) || 
+			!newDict.write(os))
+		{
+			fatalErrorInFunction<<
+			" error in writing to dictionary "<< newDict.globalName()<<endl;
+			return false;
+		}
+	}
+    return true;
+}

src/MotionModel/conveyorBeltMotion/conveyorBeltMotion.hpp

@@ -0,0 +1,104 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __conveyorBeltMotion_hpp__
+#define __conveyorBeltMotion_hpp__
+
+
+#include "MotionModel.hpp"
+#include "conveyorBelt.hpp"
+#include "fileDictionary.hpp"
+
+namespace pFlow
+{
+
+
+/**
+ * conveyor belt model for walls
+ * 
+ * This class is used for simulaiton that wall components are conveyor belt.
+ * 
+\verbatim
+
+// In geometryDict file, this will defines stationary walls
+...
+motionModel conveyorBelt;
+
+// this dictionary is optional
+conveyorBeltInfo
+{
+  conveyorBelt1
+  {
+    // the definition based on class conveyorBelt
+  }
+}
+...
+\endverbatim
+ *
+ */
+class conveyorBeltMotion
+:
+    public fileDictionary,
+    public MotionModel<conveyorBeltMotion, conveyorBelt>
+{
+protected:
+
+	friend MotionModel<conveyorBeltMotion, conveyorBelt>;
+  
+	bool impl_isMoving()const
+	{
+		return false;
+	}
+
+	bool impl_move(uint32, real, real)const
+	{
+		return true;
+	}
+	
+	void impl_setTime(uint32 ,real ,real )const
+	{}
+
+public:
+
+	TypeInfo("conveyorBeltMotion");
+
+	conveyorBeltMotion(const objectFile& objf, repository* owner);
+
+	conveyorBeltMotion(
+		const objectFile& objf, 
+		const dictionary& dict, 
+		repository* owner);
+
+
+    using fileDictionary::write;
+
+	bool write(iOstream& os, const IOPattern& iop)const override;
+
+	static
+    auto noneComponent()
+    {
+      return conveyorBelt();
+    }
+};
+
+} // pFlow
+
+
+#endif // __conveyorBeltMotion_hpp__

src/MotionModel/entities/conveyorBelt/conveyorBelt.cpp

@@ -18,75 +18,53 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+#include "conveyorBelt.hpp"
+#include "dictionary.hpp"
 
-#include "eventSubscriber.hpp"
-#include "Set.hpp"
 
-pFlow::eventSubscriber::~eventSubscriber()
+FUNCTION_H
+pFlow::conveyorBelt::conveyorBelt(const dictionary& dict)
 {
-	for(auto& observer:observerList_)
+  if(!read(dict))
 	{
-		observer->invalidateSubscriber();
+		fatalErrorInFunction<<
+		"  error in reading conveyorBelt from dictionary "<< dict.globalName()<<endl;
+		fatalExit;
 	}
 }
 
-bool pFlow::eventSubscriber::subscribe
-(
-	eventObserver* observer
-)const
+FUNCTION_H
+bool pFlow::conveyorBelt::read(const dictionary& dict)
+{	
+
+	tangentVelocity_ 		= dict.getVal<realx3>("tangentVelocity");
+	
+	return true;
+}
+
+FUNCTION_H
+bool pFlow::conveyorBelt::write(dictionary& dict) const
 {
-	if(observer)
-	{
-		observerList_.push_back(observer);
-		return true;
-	}
-	else
+	if( !dict.add("tangentVelocity", tangentVelocity_) )
 	{
+		fatalErrorInFunction<<
+		"  error in writing tangentVelocity to dictionary "<< dict.globalName()<<endl;
 		return false;
 	}
-}
-
-bool pFlow::eventSubscriber::unsubscribe
-(
-	eventObserver* observer
-)const
-{
-	if(observer)
-	{
-		observerList_.remove(observer);
-	}
 	return true;
 }
 
-bool pFlow::eventSubscriber::notify
-(
-	const eventMessage &msg
-)
+FUNCTION_H
+bool pFlow::conveyorBelt::read(iIstream& is)
 {
-	for ( auto& observer:observerList_ )
-	{
-		if(observer)
-			if( !observer->update(msg) ) return false;
-	}
-
+	
+	notImplementedFunction;
 	return true;
 }
 
-bool pFlow::eventSubscriber::notify
-(
-	const eventMessage& msg,
-	const List<eventObserver*>& exclutionList
-)
+FUNCTION_H
+bool pFlow::conveyorBelt::write(iOstream& os)const
 {
-	Set<eventObserver*> sortedExcList(exclutionList.begin(),exclutionList.end());
-
-	for(auto& observer:observerList_)
-	{
-		if( observer && sortedExcList.count(observer) == 0 )
-		{
-			if(!observer->update(msg)) return false;
-		}
-	}
-
-	return true;
+  os.writeWordEntry("tangentVelocity", tangentVelocity_);
+  return true;
 }

src/MotionModel/entities/conveyorBelt/conveyorBelt.hpp

@@ -0,0 +1,108 @@
+/*------------------------------- 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 __conveyorBelt_hpp__
+#define __conveyorBelt_hpp__
+
+
+#include "types.hpp"
+#include "typeInfo.hpp"
+#include "streams.hpp"
+
+namespace pFlow
+{
+
+// forward
+class dictionary;
+
+
+/**
+ * conveyor belt model for a wall
+ * 
+ */
+class conveyorBelt
+{
+private:
+
+	realx3  tangentVelocity_{0, 0, 0};
+
+public:
+
+	TypeInfoNV("conveyorBelt");
+
+	FUNCTION_HD
+	conveyorBelt()=default;
+
+	FUNCTION_H
+	explicit conveyorBelt(const dictionary& dict);
+	
+
+	FUNCTION_HD
+	conveyorBelt(const conveyorBelt&) = default;
+
+	conveyorBelt& operator=(const conveyorBelt&) = default;
+
+	INLINE_FUNCTION_HD
+	void setTime(real t)
+	{}
+
+	INLINE_FUNCTION_HD
+	realx3 linVelocityPoint(const realx3 &)const
+	{
+		return tangentVelocity_;
+	}
+
+	INLINE_FUNCTION_HD
+	realx3 transferPoint(const realx3& p, real)const
+	{
+		return p;
+	}
+
+	// - IO operation
+	FUNCTION_H 
+	bool read(const dictionary& dict);
+
+	FUNCTION_H
+	bool write(dictionary& dict) const;
+
+	FUNCTION_H
+	bool read(iIstream& is);
+
+	FUNCTION_H
+	bool write(iOstream& os)const;
+
+};
+
+inline iOstream& operator <<(iOstream& os, const conveyorBelt& obj)
+{
+	
+	return os;
+}
+
+inline iIstream& operator >>(iIstream& is, conveyorBelt& obj)
+{
+	
+	return is;
+}
+
+}
+
+
+#endif

src/MotionModel/entities/multiRotatingAxis/multiRotatingAxis.cpp

@@ -22,31 +22,32 @@ Licence:
 #include "multiRotatingAxisMotion.hpp"
 #include "dictionary.hpp"
 
+/// Construct from dictionary
 FUNCTION_H
-pFlow::multiRotatingAxis::multiRotatingAxis
-(
-	multiRotatingAxisMotion* axisMotion
-)
-{
-	//axisList_ = axisMotion->getAxisListPtr();	
-}
+pFlow::multiRotatingAxis::multiRotatingAxis(const dictionary& dict)
+:
+	rotatingAxis(dict)
+{}
+
 
 FUNCTION_H
 pFlow::multiRotatingAxis::multiRotatingAxis
 (
-	multiRotatingAxisMotion* axisMotion,
+	multiRotatingAxis* axisListPtr,
+	const wordList& componentsNames,
 	const dictionary& dict
 )
+:
+	rotatingAxis(dict),
+	axisList_(axisListPtr)
 {
 
-	if(!read(axisMotion, dict))
+	if(!read(dict, componentsNames))
 	{
 		fatalErrorInFunction<<
 		"  error in reading rotatingAxis from dictionary "<< dict.globalName()<<endl;
 		fatalExit;
 	}
-
-	//axisList_ = axisMotion->getAxisListPtr();
 }
 
 
@@ -54,22 +55,29 @@ pFlow::multiRotatingAxis::multiRotatingAxis
 FUNCTION_H
 bool pFlow::multiRotatingAxis::read
 (
-	multiRotatingAxisMotion* axisMotion,
-	const dictionary& dict
+	const dictionary& dict,
+	const wordList& componentNames
 )
 {
-
-	if(!rotatingAxis::read(dict))return false;
 			
 	word rotAxis = dict.getValOrSet<word>("rotationAxis", "none");
 
 	if(rotAxis == "none")
 	{
-		parentAxisIndex_ = -1;
+		parentAxisIndex_ = static_cast<uint32>(-1);
 	}
 	else
 	{
-		parentAxisIndex_ = axisMotion-> nameToIndex(rotAxis);
+		if( auto i = componentNames.findi(rotAxis); i != -1 )
+		{
+			parentAxisIndex_ = i;
+		}
+		else
+		{
+			fatalErrorInFunction<<"crotationAxis "<< rotAxis<<" in dictionary "<<
+			dict.globalName()<<" is not found in list of axis names "<< componentNames<<endl;
+			return false;
+		}
 	}
 
 	return true;
@@ -78,8 +86,8 @@ bool pFlow::multiRotatingAxis::read
 FUNCTION_H
 bool pFlow::multiRotatingAxis::write
 (
-	const multiRotatingAxisMotion* axisMotion,
-	dictionary& dict
+	dictionary& dict, 
+	const wordList& componentNames
 ) const
 {
 	if( !rotatingAxis::write(dict) ) return false;
@@ -90,10 +98,8 @@ bool pFlow::multiRotatingAxis::write
 	}
 	else
 	{
-		auto rotAxis  = axisMotion->indexToName(parentAxisIndex_);
-		dict.add("rotationAxis", rotAxis);	
+		dict.add("rotationAxis", componentNames[parentAxisIndex_]);	
 	}
 	
 	return true;
 }
-

src/MotionModel/entities/multiRotatingAxis/multiRotatingAxis.hpp

@@ -24,7 +24,7 @@ Licence:
 
 #include "rotatingAxis.hpp"
 #include "KokkosTypes.hpp"
-
+#include "List.hpp"
 
 namespace pFlow
 {
@@ -79,26 +79,31 @@ class multiRotatingAxis
 protected:
 
 	/// This is device pointer to all axes
-	multiRotatingAxis*		axisList_;
+	multiRotatingAxis*		axisList_ = nullptr;
 
 	/// Index of parent axis
-	int32 					parentAxisIndex_ = -1;
+	uint32 					parentAxisIndex_ = static_cast<uint32>(-1);
 
 public:
 
+	TypeInfoNV("multiRotatingAxis");
+
 	// - Constructors
 
 		/// Empty Constructor
-		INLINE_FUNCTION_HD
-		multiRotatingAxis(){}
+		FUNCTION_HD
+		multiRotatingAxis() = default;
 
-		/// Empty with list of axes
+		/// Construct from dictionary
 		FUNCTION_H
-		multiRotatingAxis(multiRotatingAxisMotion* axisMotion);
+		explicit multiRotatingAxis(const dictionary& dict);
 
 		/// Construct from dictionary and list of axes
 		FUNCTION_H
-		multiRotatingAxis(multiRotatingAxisMotion* axisMotion, const dictionary& dict);
+		multiRotatingAxis(
+			multiRotatingAxis* axisListPtr, 
+			const wordList& componentsNames,
+			const dictionary& dict);
 		
 		/// Copy constructor 
 		FUNCTION_HD
@@ -123,11 +128,11 @@ public:
 			while(parIndex != -1)
 			{
 				auto& ax = axisList_[parIndex];
-				parentVel += ax.linTangentialVelocityPoint(p);
+				parentVel += ax.linVelocityPoint(p);
 				parIndex = ax.parentAxisIndex();
 			}
 
-			return parentVel + rotatingAxis::linTangentialVelocityPoint(p);
+			return parentVel + rotatingAxis::linVelocityPoint(p);
 		}
 
 		/// Translate point p for dt seconds based on the axis information
@@ -143,7 +148,7 @@ public:
 			}
 
 			auto parIndex = parentAxisIndex_;
-			while(parIndex != -1)
+			while(parIndex != static_cast<uint32>(-1))
 			{
 				auto& ax = axisList_[parIndex];
 				newP = pFlow::rotate(newP, ax, dt);
@@ -157,12 +162,12 @@ public:
 		INLINE_FUNCTION_HD
 		bool hasParent()const
 		{
-			return parentAxisIndex_ > -1;
+			return parentAxisIndex_ != static_cast<uint32>(-1);
 		}
 
 		/// Return the index of parent axis
 		INLINE_FUNCTION_HD
-		int32 parentAxisIndex()const
+		uint32 parentAxisIndex()const
 		{
 			return parentAxisIndex_;
 		}
@@ -182,6 +187,7 @@ public:
 		 * It is assumed that the axis with deepest level (with more parrents) is
 		 * moved first and then the axis with lower levels.
 		 */
+		INLINE_FUNCTION_HD
 		void move(real dt)
 		{
 			
@@ -201,11 +207,12 @@ public:
 
 		/// Read from dictionary
 		FUNCTION_H 
-		bool read(multiRotatingAxisMotion* axisMotion, const dictionary& dict);
+		bool read(const dictionary& dict, const wordList& componentNames);
 
 		/// Write to dictionary
 		FUNCTION_H
-		bool write(const multiRotatingAxisMotion* axisMotion, dictionary& dict) const;
+		bool write(dictionary& dict, const wordList& componentNames) const;
+
 
 };
 

src/MotionModel/entities/rotatingAxis/rotatingAxis.cpp

@@ -47,7 +47,7 @@ pFlow::rotatingAxis::rotatingAxis
 	timeInterval(),
 	line(p1, p2),
 	omega_(omega),
-	rotating_(!equal(omega,0.0))
+	rotating_(!equal(omega,static_cast<real>(0.0)))
 {
 
 }
@@ -58,7 +58,7 @@ pFlow::real pFlow::rotatingAxis::setOmega(real omega)
 {
 	auto tmp = omega_;
 	omega_ = omega;
-	rotating_ = !equal(omega, 0.0);
+	rotating_ = !equal(omega, static_cast<real>(0.0));
 	return tmp;
 }
 
@@ -72,7 +72,7 @@ bool pFlow::rotatingAxis::read
 	if(!timeInterval::read(dict))return false;
 	if(!line::read(dict)) return false;
 		
-	real omega = dict.getValOrSet("omega", 0.0);
+	real omega = dict.getValOrSet("omega", static_cast<real>(0.0));
 	
 	setOmega(omega);
 	return true;

src/MotionModel/multiRotatingAxisMotion/multiRotatingAxisMotion.cpp

@@ -19,40 +19,63 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "multiRotatingAxisMotion.hpp"
-#include "dictionary.hpp"
-#include "vocabs.hpp"
 
-
-bool pFlow::multiRotatingAxisMotion::readDictionary
+void pFlow::multiRotatingAxisMotion::impl_setTime
 (
-	const dictionary& dict
-)
+	uint32 iter, 
+	real t, 
+	real dt
+)const
 {
+	auto motion = motionComponents_.deviceViewAll();
+	Kokkos::parallel_for(
+		"multiRotatingAxisMotion::impl_setTime",
+		deviceRPolicyStatic(0, numComponents_),
+		LAMBDA_D(uint32 i){
+			motion[i].setTime(t);
+		});
+	Kokkos::fence();
+}
 
-	auto motionModel = dict.getVal<word>("motionModel");
+bool pFlow::multiRotatingAxisMotion::impl_move(uint32 iter, real t , real dt ) const
+{
+	auto motion = motionComponents_.deviceViewAll();
+	Kokkos::parallel_for(
+		"multiRotatingAxisMotion::impl_move",
+		deviceRPolicyStatic(0, numComponents_),
+		LAMBDA_D(uint32 i){
+			motion[i].move(dt);
+		});
+	Kokkos::fence();
+	return true;
+}
 
-	if(motionModel != "multiRotatingAxisMotion")
+bool pFlow::multiRotatingAxisMotion::impl_readDictionary(const dictionary &dict)
+{
+	auto modelName = dict.getVal<word>("motionModel");
+
+	if(modelName != getTypeName<ModelComponent>())
 	{
 		fatalErrorInFunction<<
-		"  motionModel should be multiRotatingAxisMotion, but found "
-		 << motionModel <<endl;
+		"  motionModel should be "<< Yellow_Text(getTypeName<ModelComponent>())<<
+		", but found "<< Yellow_Text(modelName)<<endl;
 		return false;
 	}
 
-	auto& motionInfo = dict.subDict("multiRotatingAxisMotionInfo");
-	auto axisNames = motionInfo.dictionaryKeywords();
-	wordList rotationAxis;
-
-	// first check if 
+	auto& motionInfo = dict.subDict(modelName+"Info");
+	auto compNames = motionInfo.dictionaryKeywords();
 	
-			
-	for(auto& aName: axisNames)
+	wordList rotationAxisNames;
+
+	
+	// in the first round read all dictionaries 			
+	for(auto& compName: compNames)
 	{
-		auto& axDict = motionInfo.subDict(aName);
+		auto& axDict = motionInfo.subDict(compName);
 		
 		if(auto axPtr = makeUnique<rotatingAxis>(axDict); axPtr)
 		{
-			rotationAxis.push_back(
+			rotationAxisNames.push_back(
 				axDict.getValOrSet<word>("rotationAxis", "none"));
 		}
 		else
@@ -63,26 +86,26 @@ bool pFlow::multiRotatingAxisMotion::readDictionary
 		}
 	}
 
-	if( !axisNames.search("none") )
+	if( !compNames.search("none") )
 	{
-		axisNames.push_back("none");
-		rotationAxis.push_back("none");
+		compNames.push_back("none");
+		rotationAxisNames.push_back("none");
 	}
 
 	using intPair = std::pair<int32, int32>;
 
 	std::vector<intPair> numRotAxis;
 
-	for(size_t i=0; i< axisNames.size(); i++)
+	for(size_t i=0; i< compNames.size(); i++)
 	{
-		word rotAxis = rotationAxis[i];
+		word rotAxis = rotationAxisNames[i];
 		int32 n=0;
 		while(rotAxis != "none")
 		{
 			n++;
-			if(int32 iAxis = axisNames.findi(rotAxis) ; iAxis != -1)
+			if(int32 iAxis = compNames.findi(rotAxis) ; iAxis != -1)
 			{
-				rotAxis = rotationAxis[iAxis];	
+				rotAxis = rotationAxisNames[iAxis];	
 			}else
 			{
 				fatalErrorInFunction<<
@@ -98,60 +121,73 @@ bool pFlow::multiRotatingAxisMotion::readDictionary
 	auto compareFunc = [](const intPair& a, const intPair& b) 
 		{ return a.first > b.first; };
 
-	algorithms::STD::sort(numRotAxis.data(), numRotAxis.size(), compareFunc);
+	std::sort(numRotAxis.begin(), numRotAxis.end(), compareFunc);
+	Vector<int> sortedIndex; 
+	componentNames_.clear();
 
-	sortedIndex_.clear();
-	axisName_.clear();
+	output<<compNames<<endl;
 	
-
 	for(auto ax:numRotAxis)
 	{
-		axisName_.push_back(axisNames[ax.second]);
-		sortedIndex_.push_back(ax.second);
+		componentNames_.push_back(compNames[ax.second]);
+		sortedIndex.push_back(ax.second);
 	}
 
-	numAxis_ = axisName_.size();
-	axis_.clear();
-	axis_.reserve(numAxis_);
+	numComponents_ = componentNames_.size();
+	motionComponents_.reserve(numComponents_);
+	sortedIndex_.assign(sortedIndex);
 
-
-	// create the actual axis vector 
-	for(auto& aName: axisName_)
+	Vector<ModelComponent> components("Read::modelComponent", 
+		compNames.size()+1, 
+		0, 
+		RESERVE());
+	
+		
+	for(auto& compName: componentNames_)
 	{
-		if(aName != "none")
+		
+		if(compName != "none")
 		{
-			auto& axDict = motionInfo.subDict(aName);
-			axis_.push_back(
-				multiRotatingAxis(this, axDict));	
+			auto& compDict = motionInfo.subDict(compName);
+			components.emplace_back(
+				motionComponents_.data(), 
+				componentNames_, 
+				compDict);
 		}
 		else
 		{
-			axis_.push_back(
-				multiRotatingAxis(this));
+			components.emplace_back(impl_noneComponent());
 		}
 		
 	}
 
-	return true;
+	motionComponents_.assign(components);
+    return true;
 }
 
-bool pFlow::multiRotatingAxisMotion::writeDictionary
+
+bool pFlow::multiRotatingAxisMotion::impl_writeDictionary
 (
 	dictionary& dict
 )const
 {
-	dict.add("motionModel", "multiRotatingAxisMotion");
+	word modelName = "multiRotatingAxis";
 
-	auto& motionInfo = dict.subDictOrCreate("multiRotatingAxisMotionInfo");
-	
-	ForAll(i, axis_)
+	dict.add("motionModel", modelName );
+
+	auto modelDictName = modelName+"Info";
+
+	auto& motionInfo = dict.subDictOrCreate(modelDictName);
+	auto hostComponents = motionComponents_.hostView();
+
+	ForAll(i, motionComponents_)
 	{
 		
-		auto& axDict = motionInfo.subDictOrCreate(axisName_[i]);
-		if( !axis_.hostVectorAll()[i].write(this,axDict))
+		auto& axDict = motionInfo.subDictOrCreate(componentNames_[i]);
+		if( !hostComponents[i].write(axDict, componentNames_))
 		{
 			fatalErrorInFunction<<
-			"  error in writing axis "<< axisName_[i] << " to dicrionary "
+			"  error in writing axis "<< componentNames_[i] << " to dicrionary "
 			<< motionInfo.globalName()<<endl;
 			return false;
 		}
@@ -160,79 +196,52 @@ bool pFlow::multiRotatingAxisMotion::writeDictionary
 	return true;
 }
 
-pFlow::multiRotatingAxisMotion::multiRotatingAxisMotion()
-{}
-
-pFlow::multiRotatingAxisMotion::multiRotatingAxisMotion
-(
-	const dictionary& dict
-)
+pFlow::multiRotatingAxisMotion::multiRotatingAxisMotion(
+    const objectFile &objf,
+    repository *owner)
+    : fileDictionary(objf, owner)
 {
-	if(! readDictionary(dict) )
+
+	if(! getModel().impl_readDictionary(*this) )
 	{
+		fatalErrorInFunction;
 		fatalExit;
 	}
 }
 
-FUNCTION_H
-bool pFlow::multiRotatingAxisMotion::move(real t, real dt)
-{
-
-	// every thing is done on host 
-	for(int32 i=0; i<numAxis_; i++)
-	{	
-		auto& ax = axis_[sortedIndex_[i]];
-		ax.setTime(t);
-		ax.setAxisList(getAxisListPtrHost());
-		ax.move(dt);
-	}
-
-	// transfer to device
-	axis_.modifyOnHost();
-	axis_.syncViews();
-
-	return true;
-}
-
-bool pFlow::multiRotatingAxisMotion::read
+pFlow::multiRotatingAxisMotion::multiRotatingAxisMotion
 (
-	iIstream& is
+	const objectFile &objf, 
+	const dictionary &dict, 
+	repository *owner
 )
+:
+fileDictionary(objf, dict, owner)
 {
-	// create an empty file dictionary
-	dictionary motionInfo(motionModelFile__, true);
-
-	// read dictionary from stream
-	if( !motionInfo.read(is) )
+	if(!getModel().impl_readDictionary(*this) )
 	{
-		ioErrorInFile(is.name(), is.lineNumber()) <<
-		"  error in reading dictionray " << motionModelFile__ <<" from file. \n";
-		return false;
+		fatalErrorInFunction;
+		fatalExit;
 	}
-
-	if( !readDictionary(motionInfo) ) return false;
-	
-	return true;
 }
 
 bool pFlow::multiRotatingAxisMotion::write
 (
-	iOstream& os
-)const
+	iOstream &os, 
+	const IOPattern &iop
+) const
 {
-	// create an empty file dictionary
-	dictionary motionInfo(motionModelFile__, true);
-
-	if( !writeDictionary(motionInfo))
+	// a global dictionary
+	dictionary newDict(fileDictionary::dictionary::name(), true);
+	if( iop.thisProcWriteData() )
 	{
-		return false;
+		if( !getModel().impl_writeDictionary(newDict) || 
+			!newDict.write(os))
+		{
+			fatalErrorInFunction<<
+			" error in writing to dictionary "<< newDict.globalName()<<endl;
+			return false;
+		}
 	}
-
-	if( !motionInfo.write(os) )
-	{
-		ioErrorInFile( os.name(), os.lineNumber() )<<
-		"  error in writing dictionray to file. \n";
-		return false;
-	}
-	return true;
+    return true;
 }

src/MotionModel/multiRotatingAxisMotion/multiRotatingAxisMotion.hpp

@@ -22,18 +22,16 @@ Licence:
 #define __multiRotatingAxisMotion_hpp__
 
 
-#include "types.hpp"
-#include "typeInfo.hpp"
-#include "VectorDual.hpp"
-#include "List.hpp"
+#include "MotionModel.hpp"
 #include "multiRotatingAxis.hpp"
+#include "fileDictionary.hpp"
 
 
 namespace pFlow
 {
 
 // forward
-class dictionary;
+// class dictionary;
 
 /**
  * Rotating axis motion model for walls
@@ -63,200 +61,55 @@ multiRotatingAxisMotionInfo
  *
  */
 class multiRotatingAxisMotion
+:
+    public fileDictionary,
+    public MotionModel<multiRotatingAxisMotion, multiRotatingAxis>
 {
-public:
-	
-	/** Motion model class to be passed to computational units/kernels for
-	 *  transfing points and returning velocities at various positions
-	 */
-	class Model
-	{
-	protected:
-		
-		deviceViewType1D<multiRotatingAxis> 	axis_;
-		int32 									numAxis_=0;
-
-	public:
-
-		INLINE_FUNCTION_HD
-		Model(deviceViewType1D<multiRotatingAxis> axis, int32 numAxis):
-			axis_(axis),
-			numAxis_(numAxis)
-		{}
-
-		INLINE_FUNCTION_HD
-		Model(const Model&) = default;
-		
-		
-		INLINE_FUNCTION_HD
-		Model& operator=(const Model&) = default;
-
-
-		INLINE_FUNCTION_HD
-		realx3 pointVelocity(int32 n, const realx3& p)const
-		{
-			return axis_[n].pointTangentialVel(p);
-		}
-
-		INLINE_FUNCTION_HD
-		realx3 operator()(int32 n, const realx3& p)const
-		{
-			return pointVelocity(n,p);
-		}
-
-		INLINE_FUNCTION_HD
-		realx3 transferPoint(int32 n, const realx3 p, real dt)const
-		{
-			return axis_[n].transferPoint(p, dt);	
-		}
-
-		INLINE_FUNCTION_HD int32 numComponents()const
-		{
-			return numAxis_;
-		}
-	};
-
 protected:
 
-	using axisVector_HD = VectorDual<multiRotatingAxis>;
+	VectorSingle<int32> 	sortedIndex_;
 
-	/// Vector of multiRotaingAxis axes
-	axisVector_HD 		axis_;
-	
-	/// Sorted index based on number of parrents each axis ha
-	VectorDual<int32> 	sortedIndex_;
+	friend MotionModel<multiRotatingAxisMotion, multiRotatingAxis>;
 
-	/// List of axes names
-	wordList  			axisName_;
+	/// is the geometry attached to this component moving
+	bool impl_isMoving()const
+	{
+		return true;
+	}
 
-	/// Number of axes
-	label  				numAxis_= 0;
 
-	/// Read from a dictionary
-	bool readDictionary(const dictionary& dict);
 
-	/// Write to a dictionary 
-	bool writeDictionary(dictionary& dict)const;
+	/// Read from dictionary 
+	bool impl_readDictionary(const dictionary& dict);
+
+	bool impl_writeDictionary(dictionary& dict)const;
 
 public:
 	
-	/// Type info
-	TypeInfoNV("multiRotatingAxisMotion");
+    TypeInfo("multiRotatingAxisMotion");
 
-	// - Constructor 
+	multiRotatingAxisMotion(const objectFile& objf, repository* owner);
 
-		/// Empty constructor 
-		FUNCTION_H
-		multiRotatingAxisMotion();
+	multiRotatingAxisMotion(
+		const objectFile& objf, 
+		const dictionary& dict, 
+		repository* owner);
 
-		/// Construct with dictionary 
-		FUNCTION_H
-		multiRotatingAxisMotion(const dictionary& dict);
+    using fileDictionary::write;
 
-		/// Copy constructor 
-		FUNCTION_H
-		multiRotatingAxisMotion(const multiRotatingAxisMotion&) = default;
+	bool write(iOstream& os, const IOPattern& iop)const override;
 
-		/// No Move 
-		multiRotatingAxisMotion(multiRotatingAxisMotion&&) = delete;
+	static
+    multiRotatingAxis noneComponent()
+    {
+        return multiRotatingAxis();
+    }
 
-		/// Copy assignment 
-		FUNCTION_H
-		multiRotatingAxisMotion& operator=(const multiRotatingAxisMotion&) = default;
+    // TODO: make this method protected
+    void impl_setTime(uint32 iter, real t, real dt)const;
 
-		/// No move assignment 
-		multiRotatingAxisMotion& operator=(multiRotatingAxisMotion&&) = delete;
-
-		/// Destructor 
-		FUNCTION_H
-		~multiRotatingAxisMotion() = default;
-
-	// - Methods 
-
-		/// Retrun motion model at time t 
-		Model getModel(real t)
-		{
-			for(int32 i= 0; i<numAxis_; i++ )
-			{
-				axis_[i].setTime(t);
-				axis_[i].setAxisList(getAxisListPtrDevice());
-			}
-			axis_.modifyOnHost();
-			axis_.syncViews();
-
-			return Model(axis_.deviceVector(), numAxis_);
-		}
-
-		/// Pointer to axis list on host side 
-		INLINE_FUNCTION_H
-		multiRotatingAxis* getAxisListPtrHost()
-		{
-			return axis_.hostVectorAll().data();
-		}
-
-		/// Pointer to axis list on device 
-		INLINE_FUNCTION_H
-		multiRotatingAxis* getAxisListPtrDevice()
-		{
-			return axis_.deviceVectorAll().data();
-		}
-
-		/// Name of motion component to index 
-		INLINE_FUNCTION_H
-		int32 nameToIndex(const word& name)const
-		{
-			if( auto i = axisName_.findi(name); i == -1)
-			{
-				fatalErrorInFunction<<
-				"axis name " << name << " does not exist. \n";
-				fatalExit;
-				return i;
-			}
-			else
-			{
-				return i;
-			}
-			
-		}
-
-		/// Index of motion component to component name 
-		INLINE_FUNCTION_H
-		word indexToName(label i)const
-		{
-			if(i < numAxis_ )
-				return axisName_[i];
-			else
-			{
-				fatalErrorInFunction<<
-				"out of range access to the list of axes " << i <<endl<<
-				" size of axes_ is "<<numAxis_<<endl;
-				fatalExit;
-				return "";
-			}
-		}
-
-		/// Is moving 
-		INLINE_FUNCTION_HD
-		bool isMoving()const
-		{
-			return true;
-		}
-
-		/// Move points 
-		FUNCTION_H
-		bool move(real t, real dt);
-		
-	// - IO operation 
-
-		/// Read from input stream is 
-		FUNCTION_H
-		bool read(iIstream& is);
-
-		/// Write to output stream os 
-		FUNCTION_H
-		bool write(iOstream& os)const;
-
-	
+	/// move the component itself
+	bool impl_move(uint32 iter, real t, real dt)const;
 };
 
 } // pFlow

src/MotionModel/stationaryWall/stationaryWall.cpp

@@ -29,7 +29,7 @@ pFlow::stationaryWall::stationaryWall
 :
 	fileDictionary(objf, owner)
 {
-
+	const auto& dummy = this->subDictOrCreate("stationaryInfo");
 	if(!impl_readDictionary(*this) )
 	{
 		fatalErrorInFunction;
@@ -46,6 +46,8 @@ pFlow::stationaryWall::stationaryWall
 :
 	fileDictionary(objf, dict, owner)
 {
+	const auto& dummy = this->subDictOrCreate("stationaryInfo");
+
 	if(!impl_readDictionary(*this) )
 	{
 		fatalErrorInFunction;

src/Particles/CMakeLists.txt

@@ -6,9 +6,23 @@ particles/shape/shape.cpp
 particles/particles.cpp
 particles/particleIdHandler/particleIdHandler.cpp
 particles/regularParticleIdHandler/regularParticleIdHandler.cpp
+
+globalDamping/globalDamping.cpp
+
 SphereParticles/sphereShape/sphereShape.cpp
 SphereParticles/sphereParticles/sphereParticles.cpp
 SphereParticles/sphereParticles/sphereParticlesKernels.cpp
+
+GrainParticles/grainShape/grainShape.cpp
+GrainParticles/grainParticles/grainParticles.cpp
+GrainParticles/grainParticles/grainParticlesKernels.cpp
+
+SphereParticles/boundarySphereParticles.cpp
+SphereParticles/boundarySphereParticlesList.cpp
+
+GrainParticles/boundaryGrainParticles.cpp
+GrainParticles/boundaryGrainParticlesList.cpp
+
 Insertion/collisionCheck/collisionCheck.cpp
 Insertion/insertionRegion/insertionRegion.cpp
 Insertion/insertion/insertion.cpp
@@ -18,7 +32,8 @@ Insertion/Insertion/Insertions.cpp
 
 if(pFlow_Build_MPI)
     list(APPEND SourceFiles 
-    particles/MPIParticleIdHandler/MPIParticleIdHandler.cpp)
+    particles/MPIParticleIdHandler/MPIParticleIdHandler.cpp
+    SphereParticles/processorBoundarySphereParticles.cpp)
 endif()
 
 set(link_libs Kokkos::kokkos phasicFlow Integration Property)

src/Particles/GrainParticles/boundaryGrainParticles.cpp

@@ -0,0 +1,64 @@
+#include "boundaryGrainParticles.hpp"
+#include "boundaryBase.hpp"
+#include "grainParticles.hpp"
+
+
+pFlow::boundaryGrainParticles::boundaryGrainParticles(
+    const boundaryBase &boundary, 
+    grainParticles &prtcls
+)
+:
+    generalBoundary(boundary, prtcls.pStruct(), "", ""),
+    particles_(prtcls)
+{
+
+}
+
+pFlow::grainParticles &pFlow::boundaryGrainParticles::Particles()
+{
+    return particles_;
+}
+
+const pFlow::grainParticles &pFlow::boundaryGrainParticles::Particles() const
+{
+    return particles_;
+}
+
+pFlow::uniquePtr<pFlow::boundaryGrainParticles> pFlow::boundaryGrainParticles::create(
+    const boundaryBase &boundary, 
+    grainParticles &prtcls
+)
+{
+	
+	word bType = angleBracketsNames2( 
+		"boundaryGrainParticles",
+		pFlowProcessors().localRunTypeName(), 
+		boundary.type());
+
+	word altBType{"boundaryGrainParticles<none>"};
+
+	if( boundaryBasevCtorSelector_.search(bType) )
+	{
+		pOutput.space(4)<<"Creating boundary "<< Green_Text(bType)<<
+		" for "<<boundary.name()<<endl;
+		return boundaryBasevCtorSelector_[bType](boundary, prtcls);
+	}
+	else if(boundaryBasevCtorSelector_.search(altBType))
+	{
+		pOutput.space(4)<<"Creating boundary "<< Green_Text(altBType)<<
+		" for "<<boundary.name()<<endl;
+		return boundaryBasevCtorSelector_[altBType](boundary, prtcls);
+	}
+	else
+	{
+		printKeys( 
+			fatalError << "Ctor Selector "<< bType<<
+			" and "<< altBType << " do not exist. \n"
+			<<"Avaiable ones are: \n",
+			boundaryBasevCtorSelector_
+		);
+		fatalExit;
+	}
+	
+	return nullptr;
+}

src/Particles/GrainParticles/boundaryGrainParticles.hpp

@@ -0,0 +1,83 @@
+
+
+#ifndef __boundaryGrainParticles_hpp__
+#define __boundaryGrainParticles_hpp__
+
+#include "generalBoundary.hpp"
+#include "virtualConstructor.hpp"
+#include "timeInfo.hpp"
+
+namespace pFlow
+{
+
+class grainParticles;
+
+class boundaryGrainParticles
+    : public generalBoundary
+{
+private:
+
+    grainParticles&        particles_;
+
+public:
+
+    /// type info
+    TypeInfo("boundaryGrainParticles<none>");
+
+    boundaryGrainParticles(
+        const boundaryBase &boundary,
+        grainParticles&  prtcls
+    );
+
+    create_vCtor(
+        boundaryGrainParticles,
+        boundaryBase,
+        (
+            const boundaryBase &boundary,
+            grainParticles&  prtcls
+        ),
+        (boundary, prtcls)
+    );
+
+    add_vCtor(
+        boundaryGrainParticles,
+        boundaryGrainParticles,
+        boundaryBase
+    );
+
+    grainParticles& Particles();
+
+    const grainParticles& Particles()const;
+
+    bool hearChanges(
+        const timeInfo& ti,
+        const message &msg,
+        const anyList &varList) override
+    {
+        return true;
+    }
+
+    virtual
+    bool acceleration(const timeInfo& ti, const realx3& g)
+    {
+        return true;
+    }
+
+    bool isActive()const override
+    {
+        return false;
+    }
+    
+    static
+    uniquePtr<boundaryGrainParticles> create(
+        const boundaryBase &boundary,
+        grainParticles&  prtcls);
+    
+};
+
+
+}
+
+
+
+#endif

src/Particles/GrainParticles/boundaryGrainParticlesList.cpp

@@ -0,0 +1,20 @@
+#include "boundaryGrainParticlesList.hpp"
+
+pFlow::boundaryGrainParticlesList::boundaryGrainParticlesList(
+    const boundaryList &bndrs, 
+    grainParticles &prtcls
+)
+:
+    boundaryListPtr(), 
+    boundaries_(bndrs)
+{
+    ForAllBoundariesPtr(i, this)
+    {
+        this->set
+        (
+            i, 
+            boundaryGrainParticles::create(boundaries_[i], prtcls)
+        );
+    }
+    
+}

src/Particles/GrainParticles/boundaryGrainParticlesList.hpp

@@ -0,0 +1,36 @@
+
+
+#ifndef __boundaryGrainParticlesList_hpp__
+#define __boundaryGrainParticlesList_hpp__
+
+#include "boundaryListPtr.hpp"
+#include "boundaryList.hpp"
+#include "boundaryGrainParticles.hpp"
+
+namespace pFlow
+{
+
+class boundaryGrainParticlesList
+:
+    public boundaryListPtr<boundaryGrainParticles>
+{
+private:
+
+    const boundaryList&         boundaries_;
+
+public:
+
+    boundaryGrainParticlesList(
+        const boundaryList& bndrs,
+        grainParticles& prtcls
+    );
+
+    ~boundaryGrainParticlesList()=default;
+
+};
+
+}
+
+
+
+#endif 

src/Particles/GrainParticles/grainParticles/grainParticles.cpp

@@ -0,0 +1,421 @@
+/*------------------------------- 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 "grainParticles.hpp"
+#include "systemControl.hpp"
+#include "vocabs.hpp"
+#include "grainParticlesKernels.hpp"
+
+bool pFlow::grainParticles::initializeParticles()
+{
+	
+	using exeSpace = typename realPointField_D::execution_space;
+	using policy = Kokkos::RangePolicy<
+			exeSpace,
+			Kokkos::IndexType<uint32>>;
+	
+	auto [minIndex, maxIndex] = minMax(shapeIndex().internal());		
+
+	if( !grains_.indexValid(maxIndex) )
+	{
+		fatalErrorInFunction<< 
+		"the maximum value of shapeIndex is "<< maxIndex << 
+		" which is not valid."<<endl;
+		return false;
+	}
+
+	auto aPointsMask 	= dynPointStruct().activePointsMaskDevice();
+	auto aRange 		= aPointsMask.activeRange();
+	
+	auto field_shapeIndex 	= shapeIndex().deviceView();
+	auto field_diameter 	= grainDiameter_.deviceView();
+	auto field_coarseGrainFactor 	= coarseGrainFactor_.deviceView();
+	auto field_mass 		= mass_.deviceView();
+	auto field_propId 		= propertyId_.deviceView();
+	auto field_I 			= I_.deviceView();
+
+	// get info from grains shape 
+	realVector_D d("diameter", grains_.boundingDiameter());
+	realVector_D coarseGrainFactor("coarse Grain Factor", grains_.coarseGrainFactor());
+	realVector_D mass("mass",grains_.mass());
+	uint32Vector_D propId("propId", grains_.shapePropertyIds());
+	realVector_D  I("I", grains_.Inertia());
+
+	auto d_d 		= d.deviceView();
+	auto d_coarseGrainFactor 		= coarseGrainFactor.deviceView();
+	auto d_mass 	= mass.deviceView();
+	auto d_propId 	= propId.deviceView();
+	auto d_I 		= I.deviceView();
+		
+	Kokkos::parallel_for(
+		"particles::initInertia",
+		policy(aRange.start(), aRange.end()),
+		LAMBDA_HD(uint32 i)
+		{
+			if(aPointsMask(i))
+			{
+				uint32 index = field_shapeIndex[i];
+				field_I[i] = d_I[index];
+				field_diameter[i] = d_d[index];
+				field_coarseGrainFactor[i] = d_coarseGrainFactor[index];
+				field_mass[i] = d_mass[index];
+				field_propId[i] = d_propId[index];
+			}
+		});
+	Kokkos::fence();
+	
+	return true;
+}
+
+bool
+pFlow::grainParticles::getParticlesInfoFromShape(
+  const wordVector& shapeNames,
+  uint32Vector&     propIds,
+  realVector&       diams,
+  realVector&       coarseGrainFactors,
+
+  realVector&       m,
+  realVector&       Is,
+  uint32Vector&     shIndex
+)
+{
+	auto numNew = static_cast<uint32>(shapeNames.size());
+
+	propIds.clear();
+	propIds.reserve(numNew);
+
+	diams.clear();
+	diams.reserve(numNew);
+
+	coarseGrainFactors.clear();
+	coarseGrainFactors.reserve(numNew);
+
+	m.clear();
+	m.reserve(numNew);
+
+	Is.clear();
+	Is.reserve(numNew);
+
+	shIndex.clear();
+	shIndex.reserve(numNew);
+	
+
+	for(const auto& name:shapeNames)
+	{
+		uint32 indx;
+		if(grains_.shapeNameToIndex(name,indx))
+		{
+			shIndex.push_back(indx);
+			Is.push_back( grains_.Inertia(indx));
+			m.push_back(grains_.mass(indx));
+			diams.push_back(grains_.boundingDiameter(indx));
+			coarseGrainFactors.push_back(grains_.coarseGrainFactor(indx));
+			propIds.push_back( grains_.propertyId(indx));
+		}
+		else
+		{	
+			fatalErrorInFunction<<"Shape name "<< name << 
+			"does not exist. The list is "<<grains_.shapeNameList()<<endl;
+			return false;
+		}
+	}
+	
+	return true;
+}
+
+pFlow::grainParticles::grainParticles(
+	systemControl &control,
+	const grainShape& gShape
+)
+:
+	particles(control, gShape),
+	grains_(gShape),
+	propertyId_
+	(
+		objectFile
+		(
+			"propertyId",
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_NEVER
+		),
+		dynPointStruct(),
+		0u
+	),
+	grainDiameter_
+	(
+		objectFile(
+			"diameter",
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_NEVER),
+		dynPointStruct(),
+		0.00000000001
+	),
+	coarseGrainFactor_
+	(
+		objectFile(
+			"coarseGrainFactor",
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_NEVER),
+		dynPointStruct(),
+		0.00000000001
+	),
+	mass_
+	(
+		objectFile(
+			"mass",
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_NEVER),
+		dynPointStruct(),
+		0.0000000001
+	),
+	I_
+	(
+		objectFile
+		(
+			"I",
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_NEVER
+		),
+		dynPointStruct(),
+		static_cast<real>(0.0000000001)
+	),
+	rVelocity_
+	(
+		objectFile
+		(
+			"rVelocity",
+			"",
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_ALWAYS
+		),
+		dynPointStruct(),
+		zero3
+	),
+	rAcceleration_
+	(
+		objectFile(
+			"rAcceleration",
+			"",
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_ALWAYS
+			),
+		dynPointStruct(),
+		zero3
+	),
+	boundaryGrainParticles_
+	(
+		dynPointStruct().boundaries(),
+		*this
+	),
+	accelerationTimer_(
+		"Acceleration", &this->timers() ),
+	intPredictTimer_(
+		"Integration-predict", &this->timers() ),
+	intCorrectTimer_(
+		"Integration-correct", &this->timers() ),
+	fieldUpdateTimer_(
+		"fieldUpdate", &this->timers() )
+{
+	
+	auto intMethod = control.settingsDict().getVal<word>("integrationMethod");
+	REPORT(1)<<"Creating integration method "<<Green_Text(intMethod)
+		  << " for rotational velocity."<<END_REPORT;
+		  
+	rVelIntegration_ = integration::create
+		(
+			"rVelocity",
+			dynPointStruct(),
+			intMethod,
+			rAcceleration_.field(),
+			control.keepIntegrationHistory()
+		);
+
+	if( !rVelIntegration_ )
+	{
+		fatalErrorInFunction<<
+		"  error in creating integration object for rVelocity. \n";
+		fatalExit;
+	}
+
+	WARNING<<"setFields for rVelIntegration_"<<END_WARNING;
+	
+	if(!initializeParticles())
+	{
+		fatalErrorInFunction;
+		fatalExit;
+	}
+	
+} 
+
+
+
+bool pFlow::grainParticles::beforeIteration()
+{
+	particles::beforeIteration();
+	intPredictTimer_.start();
+		auto dt = this->dt();
+		dynPointStruct().predict(dt);
+		rVelIntegration_().predict(dt,rVelocity_, rAcceleration_);
+	intPredictTimer_.end();
+
+	fieldUpdateTimer_.start();
+	propertyId_.updateBoundariesSlaveToMasterIfRequested();
+	grainDiameter_.updateBoundariesSlaveToMasterIfRequested();
+	coarseGrainFactor_.updateBoundariesSlaveToMasterIfRequested();
+	mass_.updateBoundariesSlaveToMasterIfRequested();
+	I_.updateBoundariesSlaveToMasterIfRequested();
+	rVelocity_.updateBoundariesSlaveToMasterIfRequested();
+	rAcceleration_.updateBoundariesSlaveToMasterIfRequested();
+	rVelIntegration_().updateBoundariesSlaveToMasterIfRequested();
+	fieldUpdateTimer_.end();
+	
+	return true;
+}
+
+bool pFlow::grainParticles::iterate() 
+{
+
+	const timeInfo ti = TimeInfo();
+	const realx3 g = control().g();
+
+	particles::iterate();
+	accelerationTimer_.start();
+		pFlow::grainParticlesKernels::acceleration(
+			g,
+			mass().deviceViewAll(),
+			contactForce().deviceViewAll(),
+			I().deviceViewAll(),
+			contactTorque().deviceViewAll(),
+			dynPointStruct().activePointsMaskDevice(),
+			acceleration().deviceViewAll(),
+			rAcceleration().deviceViewAll()
+			);
+		
+		ForAllActiveBoundaries(i,boundaryGrainParticles_)
+		{
+			boundaryGrainParticles_[i].acceleration(ti, g);
+		}
+	accelerationTimer_.end();
+	
+	intCorrectTimer_.start();
+	
+		if(!dynPointStruct().correct(ti.dt()))
+		{
+			return false;
+		}
+		real damping = dynPointStruct().dampingFactor(ti);
+		if(!rVelIntegration_().correct(
+			ti.dt(), 
+			rVelocity_, 
+			rAcceleration_,
+			damping))
+		{
+			return false;
+		}
+	
+	intCorrectTimer_.end();
+	
+	return true;
+}
+
+bool pFlow::grainParticles::insertParticles
+(
+	const realx3Vector &position, 
+	const wordVector &shapesNames, 
+	const anyList &setVarList
+)
+{
+	anyList newVarList(setVarList);
+
+	realVector mass("mass");
+	realVector I("I");
+	realVector diameter("diameter");
+	realVector coarseGrainFactor("coarseGrainFactor");
+	uint32Vector propId("propId");
+	uint32Vector shapeIdx("shapeIdx");
+
+	if(!getParticlesInfoFromShape(
+		shapesNames,
+		propId,
+		diameter,
+		coarseGrainFactor,
+		mass,
+		I,
+		shapeIdx))
+	{
+		return false;
+	}
+
+	newVarList.emplaceBack(
+		mass_.name()+"Vector",
+		std::move(mass));
+
+	newVarList.emplaceBack(
+		I_.name()+"Vector",
+		std::move(I));
+	
+	newVarList.emplaceBack(
+		grainDiameter_.name()+"Vector",
+		std::move(diameter));
+
+	newVarList.emplaceBack(
+		coarseGrainFactor_.name()+"Vector",
+		std::move(coarseGrainFactor));
+
+	newVarList.emplaceBack(
+		propertyId_.name()+"Vector",
+		std::move(propId));
+
+	newVarList.emplaceBack(
+		shapeIndex().name()+"Vector",
+		std::move(shapeIdx));
+	
+	if(!dynPointStruct().insertPoints(position, newVarList))
+	{
+		return false;
+	}
+
+    return true;
+}
+
+pFlow::word pFlow::grainParticles::shapeTypeName()const
+{
+	return "grain";
+}
+
+const pFlow::shape &pFlow::grainParticles::getShapes() const
+{
+    return grains_;
+}
+
+void pFlow::grainParticles::boundingSphereMinMax
+(
+	real & minDiam, 
+	real& maxDiam
+)const
+{
+	minDiam = grains_.minBoundingSphere();
+	maxDiam = grains_.maxBoundingSphere();
+}

src/Particles/GrainParticles/grainParticles/grainParticles.hpp

@@ -0,0 +1,247 @@
+/*------------------------------- 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.
+
+-----------------------------------------------------------------------------*/
+/**
+ * @class pFlow::sphereParticles
+ *
+ * @brief Class for managing spherical particles
+ *
+ * This is a top-level class that contains the essential components for
+ * defining spherical prticles in a DEM simulation.
+ */
+
+#ifndef __grainParticles_hpp__
+#define __grainParticles_hpp__
+
+#include "indexContainer.hpp"
+#include "particles.hpp"
+#include "property.hpp"
+#include "grainShape.hpp"
+#include "boundaryGrainParticlesList.hpp"
+#include "systemControl.hpp"
+
+namespace pFlow
+{
+
+class grainParticles : public particles
+{
+public:
+
+	using ShapeType = grainShape;
+
+private:
+
+	/// reference to shapes
+	const ShapeType&       grains_;
+
+	/// property id on device
+	uint32PointField_D     propertyId_;
+
+	/// diameter / boundig sphere size of particles on device
+	realPointField_D       grainDiameter_;
+
+	realPointField_D       coarseGrainFactor_;
+
+
+	/// mass of particles field
+	realPointField_D       mass_;
+
+	/// pointField of inertial of particles
+	realPointField_D       I_;
+
+	/// pointField of rotational Velocity of particles on device
+	realx3PointField_D     rVelocity_;
+
+	/// pointField of rotational acceleration of particles on device
+	realx3PointField_D     rAcceleration_;
+
+	/// boundaries 
+	boundaryGrainParticlesList boundaryGrainParticles_;
+
+	/// rotational velocity integrator
+	uniquePtr<integration> rVelIntegration_ = nullptr;
+
+	/// timer for acceleration computations
+	Timer                  accelerationTimer_;
+
+	/// timer for integration computations (prediction step)
+	Timer                  intPredictTimer_;
+
+	/// timer for integration computations (correction step)
+	Timer                  intCorrectTimer_;
+
+	Timer  				   fieldUpdateTimer_;
+
+private:
+	
+
+	bool getParticlesInfoFromShape(
+	  const wordVector& shapeNames,
+	  uint32Vector&     propIds,
+	  realVector&       diams,
+	  realVector&       coarseGrainFactors,
+	  realVector&       m,
+	  realVector&       Is,
+	  uint32Vector&     shIndex
+	);
+	
+protected:
+
+	Timer& accelerationTimer()
+	{
+		return accelerationTimer_;
+	}
+
+	Timer& intCorrectTimer()
+	{
+		return intCorrectTimer_;
+	}
+
+	integration& rVelIntegration()
+	{
+		return rVelIntegration_();
+	}
+
+public:
+
+	/// construct from systemControl and property
+	grainParticles(systemControl& control, const grainShape& gShape);
+
+	~grainParticles() override = default;
+
+	/**
+	 * Insert new particles in position with specified shapes
+	 *
+	 * This function is involked by inserted object to insert new set of
+	 * particles into the simulation. \param position position of new particles
+	 * \param shape shape of new particles
+	 * \param setField initial value of the selected fields for new particles
+	 */
+	/*bool insertParticles
+	(
+    const realx3Vector& position,
+	    const wordVector&  shapes,
+	    const setFieldList& setField
+	) override ;*/
+
+    // TODO: make this method private later 
+    bool initializeParticles();
+
+	/// const reference to shapes object
+	const auto& grains() const
+	{
+		return grains_;
+	}
+
+	/// const reference to inertia pointField
+	const auto& I() const
+	{
+		return I_;
+	}
+
+	/// reference to inertia pointField
+	auto& I()
+	{
+		return I_;
+	}
+
+	const auto& rVelocity() const
+	{
+		return rVelocity_;
+	}
+
+	auto& rVelocity()
+	{
+		return rVelocity_;
+	}
+
+	bool hearChanges(
+	  const timeInfo& ti,
+	  const message& msg,
+	  const anyList& varList
+	) override
+	{
+		notImplementedFunction;
+		return false;
+	}
+
+	const uint32PointField_D& propertyId() const override
+	{
+		return propertyId_;
+	}
+
+	const realPointField_D& diameter() const override
+	{
+		return grainDiameter_;
+	}
+
+	const realPointField_D& grainDiameter()const
+	{
+		return grainDiameter_;
+	}
+
+	const realPointField_D& coarseGrainFactor() const 
+	{
+		return coarseGrainFactor_;
+	}
+
+
+	const realPointField_D& mass() const override
+	{
+		return mass_;
+	}
+
+	/// before iteration step
+	bool beforeIteration() override;
+
+	/// iterate particles
+	bool iterate() override;
+
+	bool insertParticles(
+	  const realx3Vector& position,
+	  const wordVector&   shapesNames,
+	  const anyList&      setVarList
+	) override;
+
+	realx3PointField_D& rAcceleration() override
+	{
+		return rAcceleration_;
+	}
+
+	const realx3PointField_D& rAcceleration() const override
+	{
+		return rAcceleration_;
+	}
+
+	const realPointField_D& boundingSphere() const override
+	{
+		return diameter();
+	}
+
+	word         shapeTypeName() const override;
+
+	const shape& getShapes() const override;
+
+	void boundingSphereMinMax(real& minDiam, real& maxDiam) const override;
+
+};// grainParticles
+
+} // pFlow
+
+#endif //__sphereParticles_hpp__

src/Particles/GrainParticles/grainParticles/grainParticlesKernels.cpp

@@ -0,0 +1,101 @@
+/*------------------------------- 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 "grainParticlesKernels.hpp"
+
+using policy = Kokkos::RangePolicy<
+			pFlow::DefaultExecutionSpace,
+			Kokkos::Schedule<Kokkos::Static>,
+			Kokkos::IndexType<pFlow::uint32>>;
+
+void pFlow::grainParticlesKernels::addMassDiamInertiaProp
+(
+    deviceViewType1D<uint32>    shapeIndex,
+    deviceViewType1D<real>  	mass,
+	deviceViewType1D<real>  	diameter,
+	deviceViewType1D<real>  	coarseGrainFactor,
+	deviceViewType1D<real>  	I,
+    deviceViewType1D<uint32>  	propertyId,
+    pFlagTypeDevice 		    incld,
+    deviceViewType1D<real>  	src_mass,
+	deviceViewType1D<real>  	src_diameter,
+	deviceViewType1D<real>  	src_I,
+    deviceViewType1D<uint32>  	src_propertyId
+)
+{
+    auto aRange = incld.activeRange();
+
+    Kokkos::parallel_for(
+		"particles::initInertia",
+		policy(aRange.start(), aRange.end()),
+		LAMBDA_HD(uint32 i)
+		{
+			if(incld(i))
+			{
+				uint32 index = shapeIndex[i];
+				I[i] = src_I[index];
+				diameter[i] = src_diameter[index];
+				mass[i] = src_mass[index];
+				propertyId[i] = src_propertyId[index];
+			}
+		});
+
+}
+
+void pFlow::grainParticlesKernels::acceleration
+( 
+	const realx3&					g,
+	const deviceViewType1D<real>&  	mass,
+	const deviceViewType1D<realx3>& force,
+	const deviceViewType1D<real>&  	I,
+	const deviceViewType1D<realx3>& torque,
+	const pFlagTypeDevice& 		    incld,
+	deviceViewType1D<realx3>  		lAcc,
+	deviceViewType1D<realx3>  		rAcc
+)
+{
+
+	auto activeRange = incld.activeRange();
+	if(incld.isAllActive())
+	{
+		Kokkos::parallel_for(
+		"pFlow::grainParticlesKernels::acceleration",
+		policy(activeRange.start(), activeRange.end()),
+		LAMBDA_HD(uint32 i){
+				lAcc[i] = force[i]/mass[i] + g;
+				rAcc[i] = torque[i]/I[i];
+		});
+	}
+	else
+	{
+		Kokkos::parallel_for(
+		"pFlow::grainParticlesKernels::acceleration",
+		policy(activeRange.start(), activeRange.end()),
+		LAMBDA_HD(uint32 i){
+			if(incld(i))
+			{
+				lAcc[i] = force[i]/mass[i] + g;
+				rAcc[i] = torque[i]/I[i];
+			}
+		});
+
+	}
+	
+	Kokkos::fence();
+}

src/Particles/GrainParticles/grainParticles/grainParticlesKernels.hpp

@@ -0,0 +1,59 @@
+/*------------------------------- 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 __grainParticlesKernels_hpp__
+#define __grainParticlesKernels_hpp__
+
+#include "types.hpp"
+#include "pointFlag.hpp"
+
+namespace pFlow::grainParticlesKernels
+{
+
+void addMassDiamInertiaProp(
+    deviceViewType1D<uint32>    shapeIndex,
+    deviceViewType1D<real>  	mass,
+    deviceViewType1D<real>  	diameter,
+        deviceViewType1D<real>    coarseGrainFactor,
+
+	deviceViewType1D<real>  	I,
+    deviceViewType1D<uint32>  	propertyId,
+    pFlagTypeDevice 		    incld,
+    deviceViewType1D<real>  	src_mass,
+	deviceViewType1D<real>  	src_grainDiameter,
+	deviceViewType1D<real>  	src_I,
+    deviceViewType1D<uint32>  	src_propertyId
+);
+
+void acceleration( 
+	const realx3&					g,
+	const deviceViewType1D<real>&  	mass,
+	const deviceViewType1D<realx3>& force,
+	const deviceViewType1D<real>&  	I,
+	const deviceViewType1D<realx3>& torque,
+	const pFlagTypeDevice& 		    incld,
+	deviceViewType1D<realx3>  		lAcc,
+	deviceViewType1D<realx3>  		rAcc
+);
+
+
+}
+
+#endif 

src/Particles/GrainParticles/grainShape/grainShape.cpp

@@ -0,0 +1,261 @@
+/*------------------------------- 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 "grainShape.hpp"
+
+
+bool pFlow::grainShape::readFromDictionary3()
+{
+
+	grainDiameters_ = getVal<realVector>("grainDiameters");
+
+	sphereDiameters_ = getVal<realVector>("sphereDiameters");
+
+	coarseGrainFactor_ = grainDiameters_ / sphereDiameters_ ;
+
+
+	if(grainDiameters_.size() != numShapes() )
+	{
+		fatalErrorInFunction<< 
+		" number of elements in grain diameters in "<< globalName()<<" is not consistent"<<endl;
+		return false;
+	}
+
+	if(sphereDiameters_.size() != numShapes() )
+	{
+		fatalErrorInFunction<< 
+		" number of elements in sphere diameters in "<< globalName()<<" is not consistent"<<endl;
+		return false;
+	}
+
+	return true;
+}
+
+bool pFlow::grainShape::writeToDict(dictionary& dict)const 
+{
+	
+	if(!shape::writeToDict(dict))return false;
+
+	return true;
+}
+
+pFlow::grainShape::grainShape
+(
+	const word& fileName,
+	repository* owner,
+	const property& prop
+)
+:
+	shape(fileName, owner, prop)
+{
+
+	if(!readFromDictionary3())
+	{
+		fatalExit;
+		fatalErrorInFunction;
+	}
+}
+
+pFlow::grainShape::grainShape
+(
+	const word &shapeType, 
+	const word &fileName, 
+	repository *owner, 
+	const property &prop
+)
+:
+	grainShape(fileName, owner, prop)
+{
+}
+
+pFlow::real pFlow::grainShape::maxBoundingSphere() const
+{
+    return max(grainDiameters_);
+}
+
+pFlow::real pFlow::grainShape::minBoundingSphere() const
+{
+    return min(grainDiameters_);
+}
+
+bool pFlow::grainShape::boundingDiameter(uint32 index, real &bDiam) const
+{
+	if( indexValid(index))
+	{
+		bDiam = grainDiameters_[index];
+		return true;
+	}
+    return false;
+}
+
+pFlow::real pFlow::grainShape::boundingDiameter(uint32 index) const
+{
+    if(indexValid(index))
+	{
+		return grainDiameters_[index];
+	}
+
+	fatalErrorInFunction
+		<<"Invalid index for diameter "
+		<<index<<endl;
+	fatalExit;
+	
+	return 0.0;
+}
+
+pFlow::realVector pFlow::grainShape::boundingDiameter() const
+{
+    return grainDiameters_;
+}
+
+pFlow::real pFlow::grainShape::coarseGrainFactor(uint32 index) const
+{
+    if(indexValid(index))
+	{
+		return coarseGrainFactor_[index];
+	}
+	fatalErrorInFunction<<"Invalid index for coarse Grain Factor "<<
+	index<<endl;
+	fatalExit;
+	return 0.0;
+}
+
+pFlow::realVector pFlow::grainShape::volume() const
+{
+    return realVector("volume", Pi/6*pow(grainDiameters_,(real)3.0));
+}
+
+pFlow::realVector pFlow::grainShape::coarseGrainFactor() const
+{
+    return coarseGrainFactor_;
+}
+
+pFlow::real pFlow::grainShape::originalDiameter(uint32 index) const
+{
+	if(indexValid(index))
+	{
+    return sphereDiameters_[index];
+	}
+    fatalErrorInFunction<<"Invalid index for sphere diameter "<<
+	index<<endl;
+	fatalExit;
+	return 0.0;
+}
+
+
+
+pFlow::realVector pFlow::grainShape::originalDiameter() const
+{
+    return sphereDiameters_;
+}
+
+bool pFlow::grainShape::mass(uint32 index, real &m) const
+{
+    if( indexValid(index) )
+	{
+		real d = grainDiameters_[index];
+		real rho = indexToDensity(index);
+		m = Pi/6.0*pow(d,3)*rho;
+		return true;
+	}
+	return false;
+}
+
+pFlow::real pFlow::grainShape::mass(uint32 index) const
+{
+    if(real m; mass(index, m))
+	{
+		return m;
+	}
+	fatalErrorInFunction<<"bad index for mass "<< index<<endl;
+	fatalExit;
+	return 0;
+}
+
+pFlow::realVector pFlow::grainShape::mass() const
+{
+    return realVector ("mass", Pi/6*pow(grainDiameters_,(real)3.0)*density());
+}
+
+pFlow::realVector pFlow::grainShape::density()const 
+{
+	auto pids = shapePropertyIds();
+	realVector rho("rho", numShapes());
+	ForAll(i, pids)
+	{
+		rho[i] = properties().density(pids[i]);
+	}
+	return rho;
+}
+
+bool pFlow::grainShape::Inertia(uint32 index, real &I) const
+{
+	if( indexValid(index) )
+	{
+		I = 0.4 * mass(index) * pow(grainDiameters_[index]/2.0,2.0);
+		return true;
+	}
+    return false;
+}
+
+pFlow::real pFlow::grainShape::Inertia(uint32 index) const
+{
+	if(real I; Inertia(index, I))
+	{
+		return I;
+	}
+    fatalExit;
+	return 0;
+}
+
+pFlow::realVector pFlow::grainShape::Inertia() const
+{
+    return realVector("I", (real)0.4*mass()*pow((real)0.5*grainDiameters_,(real)2.0));
+}
+
+bool pFlow::grainShape::Inertia_xx(uint32 index, real &Ixx) const
+{
+    return Inertia(index,Ixx);
+}
+
+pFlow::real pFlow::grainShape::Inertial_xx(uint32 index) const
+{
+    return Inertia(index);
+}
+
+bool pFlow::grainShape::Inertia_yy(uint32 index, real &Iyy) const
+{
+    return Inertia(index,Iyy);
+}
+
+pFlow::real pFlow::grainShape::Inertial_yy(uint32 index) const
+{
+    return Inertia(index);
+}
+
+bool pFlow::grainShape::Inertia_zz(uint32 index, real &Izz) const
+{
+    return Inertia(index,Izz);
+}
+
+pFlow::real pFlow::grainShape::Inertial_zz(uint32 index) const
+{
+    return Inertia(index);
+}

src/Particles/GrainParticles/grainShape/grainShape.hpp

@@ -0,0 +1,128 @@
+/*------------------------------- 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 __grainShape_hpp__
+#define __grainShape_hpp__
+
+#include "shape.hpp"
+
+namespace pFlow
+{
+
+class grainShape
+:
+	public shape
+{
+private:
+	
+	/// diameter of grains
+	realVector 			grainDiameters_;
+
+	/// diameter of spheres
+	realVector 			sphereDiameters_;
+	
+	/// course-grain factor 
+	realVector 			coarseGrainFactor_;
+
+	
+	bool readFromDictionary3();
+
+protected:
+
+	bool writeToDict(dictionary& dict)const override;
+
+public:
+
+	// - type info
+	TypeInfo("shape<grain>");
+
+	grainShape(
+		const word& fileName,
+		repository* owner,
+		const property& prop);
+
+	grainShape(
+		const word& shapeType,
+		const word& fileName,
+		repository* owner,
+		const property& prop);
+
+	
+	~grainShape() override = default;
+
+	add_vCtor
+	(
+		shape, 
+		grainShape,
+		word
+	);
+	//// - Methods
+
+	real maxBoundingSphere()const override;
+
+    real minBoundingSphere()const override;
+
+	bool boundingDiameter(uint32 index, real& bDiam)const override;
+
+	real boundingDiameter(uint32 index)const override;
+
+	realVector boundingDiameter()const override;
+	
+	realVector volume()const override;
+
+	real coarseGrainFactor(uint32 index)const ;
+
+	realVector coarseGrainFactor()const ;
+
+	real originalDiameter(uint32 index)const ;
+
+	realVector originalDiameter()const ;
+
+	bool mass(uint32 index, real& m)const override;
+
+	real mass(uint32 index) const override;
+
+	realVector mass()const override;
+
+    realVector density() const override;
+
+	bool Inertia(uint32 index, real& I)const override;
+
+	real Inertia(uint32 index)const override;
+
+	realVector Inertia()const override;
+	 
+    bool Inertia_xx(uint32 index, real& Ixx)const override;
+
+    real Inertial_xx(uint32 index)const override;
+	
+    bool Inertia_yy(uint32 index, real& Iyy)const override;
+
+    real Inertial_yy(uint32 index)const override;
+
+    bool Inertia_zz(uint32 index, real& Izz)const override;
+
+    real Inertial_zz(uint32 index)const override;
+
+};
+
+} // pFlow
+
+#endif //__grainShape_hpp__

src/Particles/Insertion/Insertion/Insertions.cpp

@@ -21,4 +21,5 @@ Licence:
 
 #include "Insertions.hpp"
 
-template class pFlow::Insertion<pFlow::sphereShape>;
+template class pFlow::Insertion<pFlow::sphereShape>;
+template class pFlow::Insertion<pFlow::grainShape>;

src/Particles/Insertion/Insertion/Insertions.hpp

@@ -24,11 +24,15 @@ Licence:
 
 #include "Insertion.hpp"
 #include "sphereShape.hpp"
+#include "grainShape.hpp"
+
 
 namespace pFlow
 {
 
 using sphereInsertion = Insertion<sphereShape> ;
+using grainInsertion = Insertion<grainShape> ;
+
 
 }
 

src/Particles/Insertion/insertion/insertion.cpp

@@ -81,7 +81,7 @@ pFlow::insertion::readInsertionDict()
 
 	if (active_)
 	{
-    	checkForCollision_ = getVal<Logical>("checkForCollision");
+    	//checkForCollision_ = getVal<Logical>("checkForCollision");
 
 		REPORT(1) << "Particle insertion mechanism is " << Yellow_Text("active")
 		          << " in the simulation." << END_REPORT;

src/Particles/Insertion/insertion/insertion.hpp

@@ -42,7 +42,7 @@ private:
 	Logical    active_ = "No";
 
 	/// Check for collision? It is not active now
-	Logical    checkForCollision_ = "No";
+	Logical    checkForCollision_ = "Yes";
 
 	/// if increase velocity in case particles are failed
 	/// to be inserted due to collision

src/Particles/Insertion/insertionRegion/insertionRegion.cpp

@@ -24,6 +24,7 @@ Licence:
 #include "particles.hpp"
 #include "twoPartEntry.hpp"
 #include "types.hpp"
+#include "shape.hpp"
 
 namespace pFlow
 {