course graining added

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,186 @@
+/*------------------------------- 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
+	(
+		real t,
+		real dt,
+		uint32 iter,
+		const message& msg, 
+    	const anyList& varList
+	) override
+    {
+		
+		pOutput<<"Function (hearChanges in boundarySphereInteractions)is not implmented Message "<<
+		 msg <<endl<<" name "<< this->boundaryName() <<" type "<< this->type()<<endl;;
+		//notImplementedFunction;
+		return true;
+	}
+
+	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
+)
+:
+	ListPtr<boundaryGrainInteraction<CFModel,gMModel>>(6),
+	boundaries_(grnPrtcls.pStruct().boundaries())
+{
+	//gSettings::sleepMiliSeconds(1000*pFlowProcessors().localRank());
+	for(uint32 i=0; i<6; i++)
+	{
+		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 "ListPtr.hpp"
+#include "boundaryGrainInteraction.hpp"
+
+
+namespace pFlow
+{
+
+
+template<typename contactForceModel,typename geometryMotionModel>
+class boundaryGrainInteractionList
+:
+    public ListPtr<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

@@ -0,0 +1,31 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+-----------------------------------------------------------------------------*/
+
+#include "boundaryGrainInteraction.hpp"
+#include "periodicBoundaryGrainInteraction.hpp"
+
+#define createBoundaryGrainInteraction(ForceModel,GeomModel) 	\
+	template class pFlow::boundaryGrainInteraction<			\
+	ForceModel,													\
+	GeomModel>;													\
+																\
+	template class pFlow::periodicBoundaryGrainInteraction<	\
+	ForceModel,													\
+	GeomModel>;
+

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,96 @@
+/*------------------------------- 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;
+	
+};
+
+}
+
+#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,359 @@
+/*------------------------------- 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>(nPrtcl+1);
+	
+	pwContactList_ = makeUnique<ContactListType>(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
+(
+	real t,
+	real dt,
+	uint32 iter,
+	const message& msg, 
+	const anyList& varList
+)
+{
+	if(msg.equivalentTo(message::ITEM_REARRANGE))
+	{
+		notImplementedFunction;
+	}
+	return true;
+}

src/Interaction/grainInteraction/grainInteraction/grainInteraction.hpp

@@ -0,0 +1,169 @@
+/*------------------------------- 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(
+		real t,
+		real dt,
+		uint32 iter,
+		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,62 @@
+/*------------------------------- 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);

src/Interaction/grainInteraction/grainInteractionsNonLinearModModels.cpp

@@ -0,0 +1,43 @@
+/*------------------------------- 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)
+
+

src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp

@@ -0,0 +1,43 @@
+/*------------------------------- 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)
+
+