Merge branch 'main' into develop

CMakeLists.txt

@@ -76,7 +76,6 @@ add_subdirectory(solvers)
 add_subdirectory(utilities)
 
 #add_subdirectory(DEMSystems)
-#add_subdirectory(testIO)
 
 install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H"
   DESTINATION include)

README.md

@@ -5,8 +5,16 @@
 
 **PhasicFlow** is a parallel C++ code for performing DEM simulations. It can run on shared-memory multi-core computational units such as multi-core CPUs or GPUs (for now it works on CUDA-enabled GPUs). The parallelization method mainly relies on loop-level parallelization on a shared-memory computational unit. You can build and run PhasicFlow in serial mode on regular PCs, in parallel mode for multi-core CPUs, or build it for a GPU device to off-load computations to a GPU. In its current statues you can simulate millions of particles (up to 80M particles tested) on a single desktop computer. You can see the [performance tests of PhasicFlow](https://github.com/PhasicFlow/phasicFlow/wiki/Performance-of-phasicFlow) in the wiki page.
 
+**MPI** parallelization with dynamic load balancing is under development. With this level of parallelization, PhasicFlow can leverage the computational power of **multi-gpu** workstations or clusters with distributed memory CPUs. 
+In summary PhasicFlow can have 6 execution modes:
+1. Serial on a single CPU core,
+2. Parallel on a multi-core computer/node (using OpenMP),
+3. Parallel on an nvidia-GPU (using Cuda),
+4. Parallel on distributed memory workstation (Using MPI)
+5. Parallel on distributed memory workstations with multi-core nodes (using MPI+OpenMP)
+6. Parallel on workstations with multiple GPUs (using MPI+Cuda).
 ## How to build?
-You can build PhasicFlow for CPU and GPU executions. [Here is a complete step-by-step procedure](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-Build-PhasicFlow).
+You can build PhasicFlow for CPU and GPU executions. The latest release of PhasicFlow is v-0.1. [Here is a complete step-by-step procedure for building phasicFlow-v-0.1.](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-Build-PhasicFlow).
 
 ## Online code documentation
 You can find a full documentation of the code, its features, and other related materials on [online documentation of the code](https://phasicflow.github.io/phasicFlow/)

cmake/bashrc

@@ -1,7 +1,7 @@
 
 export pFlow_PROJECT_VERSION=v-1.0
 
-export pFlow_PROJECT=phasicFlow
+export pFlow_PROJECT="phasicFlow-$pFlow_PROJECT_VERSION"
 
 
 projectDir="$HOME/PhasicFlow" 

solvers/CMakeLists.txt

@@ -5,3 +5,5 @@ add_subdirectory(iterateSphereParticles)
 add_subdirectory(iterateGeometry)
 
 add_subdirectory(sphereGranFlow)
+
+add_subdirectory(grainGranFlow)

solvers/grainGranFlow/CMakeLists.txt

@@ -0,0 +1,7 @@
+
+set(source_files
+grainGranFlow.cpp
+)
+set(link_lib Kokkos::kokkos phasicFlow Particles Geometry Property Interaction Interaction Utilities)
+
+pFlow_make_executable_install(grainGranFlow source_files link_lib)

solvers/grainGranFlow/createDEMComponents.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.
+
+-----------------------------------------------------------------------------*/
+
+//
+REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
+pFlow::grainParticles grnParticles(Control, proprties);
+
+//
+REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
+/*auto& sphInsertion = 
+	Control.caseSetup().emplaceObject<sphereInsertion>(
+		objectFile(
+			insertionFile__,
+			"",
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_ALWAYS
+			),
+		sphParticles,
+		sphParticles.shapes()
+	);*/
+
+auto grnInsertion =  pFlow::grainInsertion(
+	grnParticles, 
+	grnParticles.grains());
+
+REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<END_REPORT;
+auto interactionPtr = pFlow::interaction::create(
+	Control,
+	grnParticles,
+	surfGeometry
+	);
+
+auto& grnInteraction = interactionPtr();

solvers/grainGranFlow/grainGranFlow.cpp

@@ -0,0 +1,123 @@
+
+/*------------------------------- 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.
+
+-----------------------------------------------------------------------------*/
+/** 
+ * \file sphereGranFlow.cpp
+ * \brief sphereGranFlow solver 
+ *
+ * This solver simulate granular flow of cohesion-less, spherical particles. 
+ * Particle insertion can be activated in the solver to gradually insert 
+ * particles into the simulation. Geometry can be defined generally with 
+ * built-in features of the code or through ASCII stl files or a combination
+ * of both. For more information refer to [tutorials/sphereGranFlow/]
+ * (https://github.com/PhasicFlow/phasicFlow/tree/main/tutorials/sphereGranFlow) folder.
+ */
+
+#include "vocabs.hpp"
+#include "phasicFlowKokkos.hpp"
+#include "systemControl.hpp"
+#include "commandLine.hpp"
+#include "property.hpp"
+#include "geometry.hpp"
+#include "grainParticles.hpp"
+#include "interaction.hpp"
+#include "Insertions.hpp"
+
+
+
+//#include "readControlDict.hpp"
+
+
+/**
+ * DEM solver for simulating granular flow of cohesion-less particles.
+ *
+ * In the root case directory just simply enter the following command to 
+ * run the simulation. For command line options use flag -h. 
+ */
+int main( int argc, char* argv[])
+{
+
+pFlow::commandLine cmds(
+		"grainGranFlow",
+		"DEM solver for non-cohesive spherical grain particles with particle insertion "
+		"mechanism and moving geometry");
+
+bool isCoupling = false;
+
+if(!cmds.parse(argc, argv)) return 0;
+
+// this should be palced in each main 
+pFlow::processors::initProcessors(argc, argv);
+pFlow::initialize_pFlowProcessors();
+#include "initialize_Control.hpp"
+	 
+	#include "setProperty.hpp"
+	#include "setSurfaceGeometry.hpp"
+	
+
+	#include "createDEMComponents.hpp"
+	
+	REPORT(0)<<"\nStart of time loop . . .\n"<<END_REPORT;
+
+	do 
+	{
+		
+		if(! grnInsertion.insertParticles( 
+			Control.time().currentIter(),
+			Control.time().currentTime(),
+			Control.time().dt()	) )
+		{
+			fatalError<<
+			"particle insertion failed in grain DFlow solver.\n";
+			return 1;
+		}
+		 
+		// set force to zero
+		surfGeometry.beforeIteration();
+
+		// set force to zero, predict, particle deletion and etc. 
+		grnParticles.beforeIteration();
+
+		grnInteraction.beforeIteration();
+		
+		grnInteraction.iterate();	
+		
+		surfGeometry.iterate();
+
+		grnParticles.iterate();
+
+		grnInteraction.afterIteration();
+
+		surfGeometry.afterIteration();
+
+		grnParticles.afterIteration();
+				
+
+	}while(Control++);
+
+	REPORT(0)<<"\nEnd of time loop.\n"<<END_REPORT;
+
+// this should be palced in each main 
+#include "finalize.hpp"
+pFlow::processors::finalizeProcessors();
+
+
+}	
+

src/Integration/CMakeLists.txt

@@ -4,7 +4,6 @@ integration/integration.cpp
 boundaries/boundaryIntegration.cpp
 boundaries/boundaryIntegrationList.cpp 
 AdamsBashforth2/AdamsBashforth2.cpp
-AdamsBashforth2/processorAB2BoundaryIntegration.cpp
 #AdamsBashforth5/AdamsBashforth5.cpp
 #AdamsBashforth4/AdamsBashforth4.cpp
 #AdamsBashforth3/AdamsBashforth3.cpp
@@ -13,6 +12,11 @@ AdamsBashforth2/processorAB2BoundaryIntegration.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/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,350 @@
+/*------------------------------- 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.0))+ pow(Pi,2.0));
+		}
+
+		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.0)*rho_[propId_i];
+		real mj = 3*Pi/4*pow(Rj,3.0)*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_,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.0)+ pow(Pi,2.0));
+
+			//REPORT(0)<<"\n en n is  :  "<<END_REPORT;
+			//REPORT(0)<< prop.en_ <<END_REPORT;
+
+		FCn = (  -pow(f_,3.0)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,1.5) * ethan_ * vrn)*Nij;
+		FCt = ( -pow(f_,3.0)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,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_) >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/cGRelativeLinearCF.hpp

@@ -0,0 +1,340 @@
+/*------------------------------- 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_ 	= 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 != 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.0))+ pow(Pi,2.0));
+					
+					
+		}
+
+		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 "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.0)*rho_[propId_i];
+		real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j];
+
+		real sqrt_meff = sqrt((mi*mj)/(mi+mj));
+		
+		
+				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_,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.0)+ pow(Pi,2.0));
+		
+
+		FCn = (  -pow(f_,1.0)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,0.5) * ethan_ * vrn)*Nij;
+		FCt = ( -pow(f_,1.0)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,0.5) *  prop.ethat_*Vt);
+		
+
+		
+		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/grainContactForceModels.hpp

@@ -0,0 +1,45 @@
+/*------------------------------- 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 "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>>;
+
+}
+
+
+
+#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,0.0) )
+			w_hat /= w_hat_mag;
+		else
+			w_hat = 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/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)
+
+

src/Particles/CMakeLists.txt

@@ -9,8 +9,17 @@ particles/regularParticleIdHandler/regularParticleIdHandler.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

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,80 @@
+
+
+#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(
+        real t,
+        real dt,
+        uint32 iter,
+        const message &msg,
+        const anyList &varList) override
+    {
+        return true;
+    }
+
+    virtual
+    bool acceleration(const timeInfo& ti, const realx3& g)
+    {
+        return true;
+    }
+
+    static
+    uniquePtr<boundaryGrainParticles> create(
+        const boundaryBase &boundary,
+        grainParticles&  prtcls);
+    
+};
+
+
+}
+
+
+
+#endif

src/Particles/GrainParticles/boundaryGrainParticlesList.cpp

@@ -0,0 +1,19 @@
+#include "boundaryGrainParticlesList.hpp"
+
+pFlow::boundaryGrainParticlesList::boundaryGrainParticlesList(
+    const boundaryList &bndrs, 
+    grainParticles &prtcls
+)
+:
+    ListPtr(bndrs.size()), 
+    boundaries_(bndrs)
+{
+    for(auto i=0; i<boundaries_.size(); i++)
+    {
+        this->set
+        (
+            i, 
+            boundaryGrainParticles::create(boundaries_[i], prtcls)
+        );
+    }
+}

src/Particles/GrainParticles/boundaryGrainParticlesList.hpp

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

src/Particles/GrainParticles/grainParticles/grainParticles.cpp

@@ -0,0 +1,422 @@
+/*------------------------------- 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 property& prop
+)
+:
+	particles(control),
+	grains_
+	(
+		shapeFile__,
+		&control.caseSetup(),
+		prop
+	),
+	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,
+			rVelocity_.field()
+		);
+
+	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, accelertion());
+		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() 
+{
+
+	timeInfo ti = TimeInfo();
+	realx3 g = control().g();
+
+	particles::iterate();
+	accelerationTimer_.start();
+		pFlow::grainParticlesKernels::acceleration(
+			g,
+			mass().deviceViewAll(),
+			contactForce().deviceViewAll(),
+			I().deviceViewAll(),
+			contactTorque().deviceViewAll(),
+			dynPointStruct().activePointsMaskDevice(),
+			accelertion().deviceViewAll(),
+			rAcceleration().deviceViewAll()
+			);
+		for(auto& bndry:boundaryGrainParticles_)
+		{
+			bndry->acceleration(ti, g);
+		}
+	accelerationTimer_.end();
+	
+	intCorrectTimer_.start();
+	
+		if(!dynPointStruct().correct(dt(), accelertion()))
+		{
+			return false;
+		}
+		if(!rVelIntegration_().correct(
+			dt(), 
+			rVelocity_, 
+			rAcceleration_))
+		{
+			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,236 @@
+/*------------------------------- 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
+	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
+	);
+	/*bool initializeParticles();
+
+	bool insertSphereParticles(
+	    const wordVector& names,
+	    const int32IndexContainer& indices,
+	    bool setId = true);
+
+	virtual uniquePtr<List<eventObserver*>> getFieldObjectList()const override;
+	*/
+
+public:
+
+	/// construct from systemControl and property
+	grainParticles(systemControl& control, const property& prop);
+
+	~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(
+	  real           t,
+	  real           dt,
+	  uint32         iter,
+	  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& 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,242 @@
+/*------------------------------- 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::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::coarseGrainFactor() const
+{
+    return coarseGrainFactor_;
+}
+
+pFlow::real pFlow::grainShape::orginalDiameter(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::orginalDiameter() 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,110 @@
+/*------------------------------- 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 spheres
+	realVector 			grainDiameters_;
+	realVector 			sphereDiameters_;
+	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() override = default;
+
+	//// - 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;
+
+		real coarseGrainFactor(uint32 index)const ;
+
+	realVector coarseGrainFactor()const ;
+
+	real orginalDiameter(uint32 index)const ;
+
+	realVector orginalDiameter()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

@@ -22,3 +22,4 @@ Licence:
 #include "Insertions.hpp"
 
 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/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp

@@ -164,10 +164,12 @@ protected:
 		return true;
 	}
 	
-	uint32 markInNegativeSide(const word& name, uint32Vector_D& markedIndices )const;
 
 public:
 
+	uint32 markInNegativeSide(const word& name, uint32Vector_D& markedIndices )const;
+
+
 	TypeInfo("boundaryBase");
 	
 	boundaryBase(

tutorials/sphereGranFlow/binarySystemOfParticles/README.md

@@ -1,6 +1,6 @@
 # Problem definition
-A rotating drum with two particle sizes is randomly filled and let it rotate to see the segregation of particles.
-The focus of this tutorial is to show how to use the pre-processing tool, `particlesPhasicFlow`, to create the initial mixture of small and large particles. 
+
+A rotating drum is randomly filled with two particle sizes and rotated to observe particle segregation. The focus of this tutorial is to show how to use the preprocessing tool `particlesPhasicFlow` to create the initial mixture of small and large particles. 
 
 **Note:** It is supposed that you have reviewed [simulating a rotating drum](https://github.com/PhasicFlow/phasicFlow/wiki/Simulating-a-rotating-drum) tutorial before starting this tutorial.
 
@@ -14,11 +14,13 @@ a view of the rotating drum with small and large particles after 7 seconds of ro
 ***
 
 # Case setup
-PhasicFlow simulation case setup is based on the text-based files that we provide in two folders located in the simulation case folder: `settings` and `caseSetup`. Here we will have a look at some important files and the rest can be found in the tutorial folder of this case setup. 
+
+In the file `caseSetup/sphereShape` two particle types with the names `smallSphere` and `largeSphere` and the diameters 3 and 5 mm are defined. 
 
 [Simulation case setup files can be found in tutorials/sphereGranFlow folder.](https://github.com/PhasicFlow/phasicFlow/tree/main/tutorials/sphereGranFlow/binarySystemOfParticles)
 ### Shape definition 
-In file `caseSetup/sphereShape`, two particle types with names `smallSphere` and `largeSphere` and diameters 3 and 5 mm are defined. 
+
+In the file `caseSetup/sphereShape` two particle types with the names `smallSphere` and `largeSphere` and the diameters 3 and 5 mm are defined. 
 
 <div align="center"> 
 in <b>caseSetup/sphereShape</b> file
@@ -31,7 +33,7 @@ materials    (prop1  prop1);       // material names for shapes
 ```
 ### Positioning and initial mixture 
 
-In dictionary `positionParticles` located in file `settings/particlesDict`, 30000 particles are located in a cylindrical region. These particles are positioned in order along `z`, `x` and then `y` axis with 0.005 m distance between their centers. 
+In the dictionary `positionParticles` located in file `settings/particlesDict`, 30000 particles are located in a cylindrical region. These particles are positioned in order along `z`, `x` and then `y` axis with 0.005 m distance between their centers. 
 
 <div align="center"> 
 in <b>settings/particlesDict</b> file
@@ -42,10 +44,16 @@ in <b>settings/particlesDict</b> file
 // positions particles 
 positionParticles
 {
-    method positionOrdered;     // ordered positioning
+    method ordered;     // other options: random or empty
 
-    maxNumberOfParticles 30001; // maximum number of particles in the simulation
-    mortonSorting Yes;          // perform initial sorting based on morton code?   
+    orderedInfo
+    {        
+        diameter 0.005;     // minimum space between centers of particles
+        numPoints 30000;    // number of particles in the simulation 
+        axisOrder (z x y);  // axis order for filling the space with particles
+    }
+
+    regionType cylinder;                                   // other options: box and sphere 
 
     cylinder  // cylinder region for positioning particles 
     {
@@ -53,17 +61,9 @@ positionParticles
         p2 (0.0 0.0 0.097);     // end point of cylinder axis (m m m)
         radius       0.117;     // radius of cylinder (m)
     }
-        
-    positionOrderedInfo
-    {        
-        diameter 0.005;     // minimum space between centers of particles
-        numPoints 30000;    // number of particles in the simulation 
-        axisOrder (z x y);  // axis order for filling the space with particles
-    }
 }
 ```
-
-In dictionary `setFields` located in file `settings/particlesDict`, you define the initial `velocity`, `acceleration`, `rotVelocity`, and  `shapeName` fields for all 30000 particles in the simulation. In `selectors` dictionary, you can select  subsets of particles and set the field value for these subsets. In `shapeAssigne` sub-dictionary, the `selectRange` selector is defined. It defines a range with `begin` (begin index), `end` (end index) and `stride` to select particles. And in `fieldValue` sub-dictionary, the fields values for selected particles are set (any number of field values can be set here).
+In the `setFields` dictionary, located in the `settings/particlesDict` file, you define the initial `velocity`, `acceleration`, `rotVelocity` and `shapeName` fields for all 30000 particles in the simulation. In the `selectors' dictionary, you can select subsets of particles and set the field value for those subsets. The `selectRange` selector is defined in the `shapeAssigne` subdictionary. It defines a range with `begin`, `end` and `stride` to select particles. And in the `fieldValue` subdictionary the field values for selected particles are set (any number of field values can be set here).
 
 **Note:** Other selectors are: `selectBox` that selects particles inside a box and `randomSelect` that selects particles randomly from a given index range. 
 
@@ -77,7 +77,8 @@ setFields
     /*
         Default value for fields defined for particles
     These fields should always be defined for simulations with 
-    spherical particles.*/
+    spherical particles.
+*/
 
     defaultValue 
     {
@@ -91,8 +92,9 @@ setFields
     {
         shapeAssigne
         {
-            selector     selectRange; // type of point selector
-            selectRangeInfo
+            selector     stridedRange; 				       // other options: box, cylinder, sphere, randomPoints
+
+            stridedRangeInfo
             {
                 begin   0;            // begin index of points
                 end     30000;        // end index of points 
@@ -101,10 +103,10 @@ setFields
             fieldValue  // fields that the selector is applied to 
             {
                 /*
-                sets shapeName of the selected points to largeSphere*/
+                    sets shapeName of the selected points to largeSphere
+                */
                 shapeName     word     largeSphere;
             }
-        
     }
 }
 

tutorials/sphereGranFlow/binarySystemOfParticles/caseSetup/interaction

@@ -6,49 +6,46 @@ objectName  interaction;
 objectType  dicrionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
+materials     (prop1);          // properties of material
 
+densities     (1000.0);         // density of materials [kg/m3]
 
-materials      (prop1);    // a list of materials names
-densities      (1000.0);   // density of materials [kg/m3]
-
-contactListType   sortedContactList; 
-
-model
-{
-   contactForceModel nonLinearNonLimited;
-   rollingFrictionModel normal;
-
-   Yeff  (1.0e6);       // Young modulus [Pa]
-
-   Geff  (0.8e6);       // Shear modulus [Pa]
-
-   nu    (0.25);        // Poisson's ratio [-]
-
-   en    (0.7);         // coefficient of normal restitution
-
-   et    (1.0);         // coefficient of tangential restitution 
-
-   mu    (0.3);         // dynamic friction 
-
-   mur   (0.1);         // rolling friction 
-        
-}
+contactListType     sortedContactList;
 
 contactSearch
 {
-   method         multiGridNBS;       // method for broad search particle-particle
-   wallMapping    multiGridMapping;   // method for broad search particle-wall 
   
-   multiGridNBSInfo
-   {
-      updateFrequency 10;        // each 10 timesteps, update neighbor list 
-      sizeRatio      1.1;        // bounding box size to particle diameter (max)
-   }
+   method          NBS;          // method for broad search 
     
-   multiGridMappingInfo
-   {
-      updateFrequency 10;        // each 10 timesteps, update neighbor list  
-      cellExtent     0.6;        // bounding box for particle-wall search (> 0.5)
-   }
+   updateInterval   10;
 
+   sizeRatio       1.1;
+
+   cellExtent     0.55;
+
+   adjustableBox   Yes;
 } 
+
+model
+{
+   contactForceModel     nonLinearNonLimited;
+
+   rollingFrictionModel               normal;
+
+   // Property (solid-solid Properties)
+
+   Yeff     (1.0e6);          // Young modulus [Pa]
+
+   Geff     (0.8e6);          // Shear modulus [Pa]
+
+   nu       (0.25);           // Poisson's ratio [-]
+
+   en       (0.7);            // coefficient of normal restitution
+
+   et       (1.0);            // coefficient of tangential restitution 
+
+   mu       (0.3);            // dynamic friction 
+
+   mur      (0.1);            // rolling friction     
+}
+

tutorials/sphereGranFlow/binarySystemOfParticles/caseSetup/particleInsertion

@@ -6,10 +6,6 @@ objectName 	particleInsertion;
 objectType 	dicrionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
+collisionCheck     No;          // is checked          -> yes or no
 
-
-active  no;          // is insertion active?
-
-collisionCheck No;   // not implemented for yes
-
-
+active             No;          // is insertion active -> yes or no

tutorials/sphereGranFlow/binarySystemOfParticles/caseSetup/shapes

@@ -6,7 +6,8 @@ objectName 	sphereDict;
 objectType 	sphereShape;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
+names 	  (smallSphere largeSphere); 	         // names of shapes 
 
-names 		(lightSphere heavySphere); 	// names of shapes 
-diameters 	(0.007 0.007);				// diameter of shapes 
-materials	(lightMat heavyMat);		// material names for shapes 
+diameters (0.003 0.005);	                     // diameter of shapes 
+
+materials (prop1 prop1);	                     // material names for shapes 

tutorials/sphereGranFlow/binarySystemOfParticles/settings/domainDict

@@ -0,0 +1,63 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName 	domainDict;
+objectType 	dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox          // Simulation domain: every particles that goes outside this domain will be deleted
+{
+	min     (-0.12 -0.12 0);
+
+	max     (0.12 0.12 0.1);
+}
+
+decomposition
+{
+	direction z;
+}
+
+boundaries
+{
+	// Determines how often (how many iterations) do you want to 
+	// rebuild the list of particles in the neighbor list 
+	// of all boundaries in the simulation domain 
+	neighborListUpdateInterval     50; 
+
+	// Determines how often do you want to update the new changes in the boundary
+	updateInterval                 10;
+
+	// The distance from the boundary plane within which particles are marked to be in the boundary list 
+	neighborLength              0.004;
+
+	left
+	{
+		type     exit;	        // other options: periodict, reflective 
+	}
+
+	right
+	{
+		type     exit;          // other options: periodict, reflective 
+	}
+
+	bottom
+	{
+		type     exit;          // other options: periodict, reflective 
+	}
+
+	top
+	{
+		type     exit;          // other options: periodict, reflective 
+	}
+
+	rear
+	{
+		type     exit;          // other options: periodict, reflective 
+	}
+
+	front
+	{
+		type     exit;          // other options: periodict, reflective 
+	}
+}

tutorials/sphereGranFlow/binarySystemOfParticles/settings/geometryDict

@@ -6,66 +6,80 @@ objectName 	geometryDict;
 objectType 	dictionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
+motionModel     rotatingAxis;          // motion model can be rotatingAxis or stationary or vibrating
 
-// motion model: rotating object around an axis
-motionModel rotatingAxisMotion; 
-
-surfaces
-{
-	/*
-	    A cylinder with begin and end radii 0.12 m and axis points at (0 0 0) 
-	and (0 0 0.1)
-	*/
-	cylinder
-	{
-		type cylinderWall;  // type of the wall
-		p1 (0.0 0.0 0.0);   // begin point of cylinder axis
-		p2 (0.0 0.0 0.1);   // end point of cylinder axis
-		radius1 0.12;		// radius at p1
-		radius2 0.12;		// radius at p2
-		resolution 24;      // number of divisions
-		material prop1;     // material name of this wall
-		motion rotAxis;		// motion component name 
-	}
-
-	/*
-	This is a plane wall at the rear end of cylinder
-	*/
-	wall1
-	{
-		type planeWall;			// type of the wall
-		p1 (-0.12 -0.12 0.0);	// first point of the wall
-		p2 ( 0.12 -0.12 0.0);   // second point
-		p3 ( 0.12  0.12 0.0);   // third point
-		p4 (-0.12  0.12 0.0);   // fourth point 
-		material prop1;         // material name of the wall  
-		motion rotAxis;			// motion component name 
-	}
-
-	/*
-	This is a plane wall at the front end of cylinder
-	*/
-	wall2
-	{
-		type planeWall;
-		p1 (-0.12 -0.12 0.1);
-		p2 ( 0.12 -0.12 0.1);
-		p3 ( 0.12  0.12 0.1);
-		p4 (-0.12  0.12 0.1);
-		material prop1;
-		motion rotAxis;
-	}
-		
-	
-}
-
-// information for rotatingAxisMotion motion model 
-rotatingAxisMotionInfo
+rotatingAxisInfo          // information for rotatingAxis motion model 
 {
 	rotAxis 
 	{
-		p1 (0.0 0.0 0.0);	// first point for the axis of rotation 
-		p2 (0.0 0.0 1.0);	// second point for the axis of rotation
-		omega 1.214; 		// rotation speed (rad/s)
+		p1 (0.0 0.0 0.0);	       // first point for the axis of rotation 
+
+		p2 (0.0 0.0 1.0);	       // second point for the axis of rotation
+		
+		omega      1.214; 		   // rotation speed (rad/s)
+	}
+}
+
+surfaces
+{
+	cylinder
+	{
+		type     cylinderWall;          // other options: cuboidWall and planeWall
+
+		p1 (0.0 0.0 0.0);          // begin point of cylinder axis
+
+		p2 (0.0 0.0 0.1);          // end point of cylinder axis
+
+		radius1     0.12;		   // radius at p1
+
+		radius2     0.12;		   // radius at p2
+
+		resolution    24;          // number of divisions
+
+		material   prop1;          // material name of this wall
+
+		motion   rotAxis;		   // motion component name 
+	}
+
+	/*
+		This is a plane wall at the rear end of cylinder
+	*/
+
+	wall1
+	{
+		type     planeWall;			 // other options: cuboidWall and  cylinderWall
+
+		p1 (-0.12 -0.12 0.0);	       // first point of the wall
+
+		p2  (0.12 -0.12 0.0);          // second point of the wall
+
+		p3   (0.12 0.12 0.0);          // third point of the wall
+
+		p4  (-0.12 0.12 0.0);          // fourth point of the wall
+
+		material       prop1;          // material name of the wall 
+
+		motion       rotAxis;		   // motion component name 
+	}
+
+	/*
+		This is a plane wall at the front end of cylinder
+	*/
+
+	wall2
+	{
+		type     planeWall;          // other options: cuboidWall and  cylinderWall 
+
+		p1 (-0.12 -0.12 0.1);          // first point of the wall
+
+		p2  (0.12 -0.12 0.1);          // second point of the wall
+
+		p3   (0.12 0.12 0.1);          // third point of the wall
+
+		p4  (-0.12 0.12 0.1);          // fourth point of the wall
+
+		material       prop1;          // material name of the wall 
+
+		motion       rotAxis;          // motion component name 
 	}
 }

tutorials/sphereGranFlow/binarySystemOfParticles/settings/particlesDict

@@ -6,64 +6,68 @@ objectName 	particlesDict;
 objectType 	dictionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
-
-// positions particles 
-positionParticles
-{
-	method positionOrdered;     // ordered positioning
-
-	maxNumberOfParticles 30001; // maximum number of particles in the simulation
-	mortonSorting Yes;          // perform initial sorting based on morton code?   
-
-	cylinder  // cylinder region for positioning particles 
-	{
-		p1 (0.0 0.0 0.003); 	// begin point of cylinder axis
-		p2 (0.0 0.0 0.097);		// end point of cylinder axis
-		radius 		0.117;		// radius of cylinder 
-	}
-		
-	positionOrderedInfo
-	{		
-		diameter 0.005; 	// minimum space between centers of particles
-		numPoints 30000; 	// number of particles in the simulation 
-		axisOrder (z x y);  // axis order for filling the space with particles
-	}
-}
-
 setFields
 {
 	/*
-	    Default value for fields defined for particles
-	These fields should always be defined for simulations with 
-	spherical particles.*/
+	    	Default value for fields defined for particles:
+		These fields should always be defined for simulations with spherical particles
+	*/
 
 	defaultValue 
 	{
-		velocity 		realx3 	(0 0 0); // linear velocity (m/s)
-		acceleration 	realx3 	(0 0 0); // linear acceleration (m/s2)
-		rVelocity 		realx3 	(0 0 0); // rotational velocity (rad/s)
-		shapeName 		word	smallSphere; // name of the particle shape 
+		velocity 		realx3 	(0 0 0); 			       // linear velocity (m/s)
+
+		acceleration 	realx3 	(0 0 0); 			       // linear acceleration (m/s2)
+
+		rVelocity 		realx3 	(0 0 0); 			       // rotational velocity (rad/s)
+
+		shapeName 		word	smallSphere; 		       // name of the particle shape 
 	}
 
 	selectors
 	{
 		shapeAssigne
 		{
-			selector 	selectRange; // type of point selector
-			selectRangeInfo
+			selector 	stridedRange; 				       // other options: box, cylinder, sphere, randomPoints 
+ 
+			stridedRangeInfo
 			{
-				begin 	0;			// begin index of points
-				end 	30000;		// end index of points 
-				stride  3;			// stride for selector 
+				begin 	0;			 				       // begin index of points
+
+				end 	30000;		 				       // end index of points 
+
+				stride  3;			 				       // stride for selector 
 			}
-			fieldValue  // fields that the selector is applied to 
+
+			fieldValue  							       // fields that the selector is applied to 
 			{
-				/*
-				sets shapeName of the selected points to largeSphere*/
-				shapeName 	word 	largeSphere;
+				shapeName 	word 	largeSphere;           // sets shapeName of the selected points to largeSphere
 			}
 		}
 	}
 }
 
+positionParticles 				                           // positions particles 
+{
+	method ordered;                                        // other options: random and empty
 
+	orderedInfo
+	{		
+		diameter  0.005; 	                               // diameter of particles
+
+		numPoints 30000; 	                               // number of particles in the simulation
+
+		axisOrder (z x y);                                 // axis order for filling the space with particles
+	}
+	
+	regionType cylinder;                                   // other options: box and sphere  
+
+	cylinderInfo  				                           // cylinder information for positioning particles 
+	{
+		p1 (0.0 0.0 0.003); 	                           // begin point of cylinder axis
+
+		p2 (0.0 0.0 0.097);		                           // end point of cylinder axis
+
+		radius 		0.117;		                           // radius of cylinder 
+	}
+}

tutorials/sphereGranFlow/binarySystemOfParticles/settings/settingsDict

@@ -6,36 +6,30 @@ objectName 	settingsDict;
 objectType 	dictionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
+run 	binarySystemofParticles;
 
+dt       0.00001; 	       // time step for integration (seconds)
 
-run 			binarySystemofParticles;
+startTime    0.0; 		   // start time for simulation 
 
-dt 				0.00001; 	// time step for integration (s)
+endTime     10.0; 		   // end time for simulation 
 
-startTime 		0; 			// start time for simulation 
+saveInterval 0.1; 		   // time interval for saving the simulation
 
-endTime 		10; 		// end time for simulation 
+timePrecision  6;		   // maximum number of digits for time folder 
 
-saveInterval 	0.1; 		// time interval for saving the simulation
+g 	  (0 -9.8 0);          // gravity vector (m/s2) 
 
-timePrecision   6;			// maximum number of digits for time folder 
+// save necessary (i.e., required) data on disk
+includeObjects (diameter);                                                 
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
 
-g 				(0 -9.8 0); // gravity vector (m/s2) 
+integrationMethod 		AdamsBashforth2; 	      // integration method 
 
-/*
-	Simulation domain 
-	every particles that goes outside this domain is deleted.
-*/
-domain 
-{
-	min (-0.12 -0.12 0);
-	max (0.12   0.12 0.1);
-}
+writeFormat 			          ascii;          // data writting format (ascii or binary)
 
-integrationMethod 		AdamsBashforth2; 	// integration method 
+timersReport 			            Yes;  	      // report timers
 
-writeFormat 			ascii;
+timersReportInterval   	           0.01;	      // time interval for reporting timers
  
-timersReport 			Yes;  	// report timers?
-
-timersReportInterval   	0.01;	// time interval for reporting timers

tutorials/sphereGranFlow/layeredSiloFilling/caseSetup/interaction

@@ -6,15 +6,29 @@ objectName interaction;
 objectType dicrionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
+materials   (lightMat heavyMat wallMat);          // a list of materials names
 
-materials   (lightMat heavyMat wallMat); // a list of materials names
-densities   (1000 1500.0 2500);      // density of materials [kg/m3]
+densities            (1000 1500.0 2500);          // density of materials [kg/m3]
 
 contactListType   sortedContactList;
 
+contactSearch
+{
+   method         NBS;                           // method for broad search particle-particle
+   
+   updateInterval 10;
+
+   sizeRatio      1.1;
+
+   cellExtent     0.55;
+
+   adjustableBox  Yes;
+} 
+
 model
 {
    contactForceModel nonLinearLimited;
+
    rollingFrictionModel normal;
 
    /*
@@ -23,51 +37,34 @@ model
                                                       wallMat-wallMat );
    */
 
-   Yeff (1.0e6 1.0e6 1.0e6  // Young modulus [Pa]
+   Yeff (1.0e6 1.0e6 1.0e6            // Young modulus [Pa]
                1.0e6 1.0e6
                      1.0e6);
    
-   Geff (0.8e6 0.8e6 0.8e6  // Shear modulus [Pa]
+   Geff (0.8e6 0.8e6 0.8e6            // Shear modulus [Pa]
                0.8e6 0.8e6
                      0.8e6);
 
-   nu    (0.25 0.25  0.25  // Poisson's ratio [-]
+   nu    (0.25 0.25  0.25             // Poisson's ratio [-]
                0.25  0.25
                      0.25);
 
-   en (0.97  0.97    0.85   // coefficient of normal restitution
+   en (0.97  0.97    0.85             // coefficient of normal restitution
              0.97    0.85
                      1.00);
 
-   et (1.0   1.0  1.0   // coefficient of tangential restitution 
+   et (1.0   1.0  1.0                 // coefficient of tangential restitution 
              1.0  1.0
                   1.0);
 
-   mu (0.65   0.65 0.35   // dynamic friction 
+   mu (0.65   0.65 0.35               // dynamic friction 
               0.65 0.35
                    0.35);
 
-   mur (0.1  0.1 0.1    // rolling friction 
+   mur (0.1  0.1 0.1                  // rolling friction 
              0.1 0.1
-                  0.1);
-        
+                 0.1);       
 }
 
-contactSearch
-{
-   method         NBS;           // method for broad search particle-particle
-   wallMapping    cellMapping;   // method for broad search particle-wall 
 
-   NBSInfo
-   {
-      updateFrequency 10;        // each 20 timesteps, update neighbor list 
-      sizeRatio      1.1;        // bounding box size to particle diameter (max)
-   }
 
-   cellMappingInfo
-   {
-      updateFrequency 10;        // each 20 timesteps, update neighbor list  
-      cellExtent     0.6;        // bounding box for particle-wall search (> 0.5)
-   }
-
-}

tutorials/sphereGranFlow/layeredSiloFilling/caseSetup/particleInsertion

@@ -6,147 +6,49 @@ objectName particleInsertion;
 objectType dicrionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
+active           Yes;          // is insertion active -> yes or no
 
-active  yes;          // is insertion active?
-
-collisionCheck No;   // not implemented for yes
+checkForCollision No;          // is checked          -> yes or no
 
 /*
-five layers of particles are packed one-by-one using 5 insertion steps. 
+   one layers of particles are packed 
 */
 
 layer0
 {  
-   type		      cylinderRegion; // type of insertion region
-   rate 	      15000;  // insertion rate (particles/s)
-   startTime 	  0;      // (s)
-   endTime   	  0.5;    // (s)
-   interval       0.025;  //s  
+   timeControl          simulationTime;
 
-   cylinderRegionInfo 
+   regionType                 cylinder;          // type of insertion region
+
+   rate                          15000;          // insertion rate (particles/s)
+
+   startTime                         0;          // (s)
+
+   endTime                         0.5;          // (s)
+
+   insertionInterval             0.025;          // s 
+
+   cylinderInfo 
    {
-   	radius 	0.09;       // radius of cylinder (m)
-   	p1 ( 0.0  0.0 0.1 ); // (m,m,m)
-   	p2 ( 0.0  0.0 0.11); // (m,m,m)
+   	radius 	                   0.09;          // radius of cylinder (m)
+
+   	p1              ( 0.0  0.0 0.1 );          // (m,m,m)
+
+   	p2              ( 0.0  0.0 0.11);          // (m,m,m)
    }
    
    setFields
    {
-     velocity    realx3 (0.0 0.0 -0.6); // initial velocity of inserted particles 
+     velocity     realx3 (0.0 0.0 -0.6);         // initial velocity of inserted particles 
    }
    
    mixture
    {
-      lightSphere 1; // mixture composition of inserted particles 
-      
+      lightSphere                    1;          // mixture composition of inserted particles 
    }
 }
 
-layer1
-{
-   type         cylinderRegion;
-   rate         15000; // (particles/s)
-   startTime    0.7; // (s)
-   endTime      1.2; // (s)
-   interval     0.025; //s  
 
-   cylinderRegionInfo 
-   {
-      radius   0.09;
-      p1 ( 0.0  0.0 0.16 ); // (m,m,m)
-      p2 ( 0.0  0.0 0.17); // (m,m,m)
-   }
 
-   setFields
-   {
-     velocity    realx3 (0.0 0.0 -0.6);
-   }
 
-   mixture
-   {
-      heavySphere 1; // only heavySphere
-      
-   }
-}
-
-layer2
-{
-   type         cylinderRegion;
-   rate         15000; // (particles/s)
-   startTime    1.4; // (s)
-   endTime      1.9; // (s)
-   interval     0.025; //s  
-
-   cylinderRegionInfo 
-   {
-      radius   0.09;
-      p1 ( 0.0  0.0 0.2 ); // (m,m,m)
-      p2 ( 0.0  0.0 0.21); // (m,m,m)
-   }
-   
-   setFields
-   {
-     velocity    realx3 (0.0 0.0 -0.6);
-   }
-   
-   mixture
-   {
-      lightSphere 1; // only lightSphere
-      
-   }
-}
-
-layer3
-{
-   type         cylinderRegion;
-   rate         15000; // (particles/s)
-   startTime    2.1; // (s)
-   endTime      2.6; // (s)
-   interval     0.025; //s  
-
-   cylinderRegionInfo 
-   {
-      radius   0.09;
-      p1 ( 0.0  0.0 0.28 ); // (m,m,m)
-      p2 ( 0.0  0.0 0.29); // (m,m,m)
-   }
-   
-   setFields
-   {
-     velocity   realx3 (0.0 0.0 -0.6);
-   }
-   
-   mixture
-   {
-      heavySphere 1;
-      
-   }
-}
-
-layer4
-{
-   type         cylinderRegion;
-   rate         15000; // (particles/s)
-   startTime    2.8; // (s)
-   endTime      3.3; // (s)
-   interval     0.025; //s  
-
-   cylinderRegionInfo 
-   {
-      radius   0.09;
-      p1 ( 0.0  0.0 0.37 ); // (m,m,m)
-      p2 ( 0.0  0.0 0.38); // (m,m,m)
-   }
-   
-   setFields
-   {
-     velocity    realx3 (0.0 0.0 -0.6);
-   }
-   
-   mixture
-   {
-      lightSphere 1;
-      
-   }
-}
 

tutorials/sphereGranFlow/layeredSiloFilling/caseSetup/shapes

@@ -6,7 +6,10 @@ objectName 	sphereDict;
 objectType 	sphereShape;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
+names 		 (lightSphere heavySphere); 	     // names of shapes 
+
+diameters 	             (0.007 0.007);		     // diameter of shapes 
+
+materials	       (lightMat heavyMat);		     // material names for shapes 
+
 
-names 		(smallSphere largeSphere); 	// names of shapes 
-diameters 	(0.003 	0.005);	// diameter of shapes 
-materials	(prop1  prop1);	// material names for shapes 

tutorials/sphereGranFlow/layeredSiloFilling/settings/domainDict

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName 	domainDict;
+objectType 	dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                    // Simulation domain: every particles that goes outside this domain will be deleted
+{
+	min (-0.11 -0.11 -0.11);
+	
+	max ( 0.11  0.11  0.41);
+}
+
+boundaries
+{
+	// Determines how often (how many iterations) do you want to 
+	// rebuild the list of particles in the neighbor list 
+	// of all boundaries in the simulation domain 
+	neighborListUpdateInterval 30;
+	
+	// Determines how often do you want to update the new changes in the boundary
+	updateInterval             10;
+
+	// The distance from the boundary plane within which particles are marked to be in the boundary list 
+	neighborLength          0.004;
+
+	left
+	{
+		type exit;	         // other options: periodict, reflective 
+	}
+
+	right
+	{
+		type exit;           // other options: periodict, reflective 
+	}
+
+	bottom
+	{
+		type exit;           // other options: periodict, reflective 
+	}
+
+	top
+	{
+		type exit;           // other options: periodict, reflective 
+	}
+
+	rear
+	{
+		type exit;           // other options: periodict, reflective 
+	}
+
+	front
+	{
+		type exit;           // other options: periodict, reflective 
+	}
+}    
+
+

tutorials/sphereGranFlow/layeredSiloFilling/settings/geometryDict

@@ -6,47 +6,71 @@ objectName 	geometryDict;
 objectType 	dictionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
+motionModel stationary;          // motion model can be rotatingAxis or stationary or vibrating
 
-// motion model: all surfaces are fixed
-motionModel fixedWall; 
+stationaryInfo
+{
+	
+}
 
 surfaces
 {
-	
 	cylinderShell
 	{
-		type 		cylinderWall;  	// type of the wall
-		p1 			(0.0 0.0 0.0);  // begin point of cylinder axis
-		p2 			(0.0 0.0 0.4);  // end point of cylinder axis
-		radius1 	0.1;		// radius at p1
-		radius2 	0.1;		// radius at p2
-		resolution 	36;      	// number of divisions
-		material 	wallMat;   	// material name of this wall
+		type 		 cylinderWall;  	    // other options: cuboidWall and planeWall
+
+		p1 			(0.0 0.0 0.0);          // begin point of cylinder axis
+
+		p2 			(0.0 0.0 0.4);          // end point of cylinder axis
+
+		radius1 	          0.1;		    // radius at p1
+ 
+		radius2 	          0.1;		    // radius at p2
+
+		resolution 	           36;      	// number of divisions
+
+		material 	      wallMat;   	    // material name of this wall
 	}
 
 	coneShell
 	{
-		type 		cylinderWall;  	// type of the wall
-		p1 			(0.0 0.0 -0.1);  // begin point of cylinder axis
-		p2 			(0.0 0.0 0.0);  // end point of cylinder axis
-		radius1 	0.02;		// radius at p1
-		radius2 	0.1;		// radius at p2
-		resolution 	36;      	// number of divisions
-		material 	wallMat;   	// material name of this wall
+		type 		 cylinderWall;  	    // other options: cuboidWall and planeWall
+
+		p1 		   (0.0 0.0 -0.1);          // begin point of cylinder axis
+
+		p2 			(0.0 0.0 0.0);          // end point of cylinder axis
+
+		radius1 	         0.02;		    // radius at p1
+
+		radius2 	          0.1;		    // radius at p2
+
+		resolution 	           36;      	// number of divisions
+
+		material 	      wallMat;   	    // material name of this wall
 	}
 
 	/*
-	This is a plane wall at the exit of silo
+		This is a plane wall at the exit of silo
 	*/
+
 	exitGate
 	{
-		type planeWall;
-		p1 (-0.02 -0.02 -0.1);
-		p2 ( 0.02 -0.02 -0.1);
-		p3 ( 0.02  0.02 -0.1);
-		p4 (-0.02  0.02 -0.1);
-		material wallMat;
-	}
+		type planeWall;                     // other options: cuboidWall and cylinderWall
 
+		p1 (-0.02 -0.02 -0.1);              // first point of the wall
+
+		p2 ( 0.02 -0.02 -0.1);              // second point of the wall
+
+		p3 ( 0.02  0.02 -0.1);              // third point of the wall
+
+		p4 (-0.02  0.02 -0.1);              // fourth point of the wall
+		 
+		material wallMat;                   // material name of the wall 
+	}
 }
 
+
+
+
+
+

tutorials/sphereGranFlow/layeredSiloFilling/settings/particlesDict

@@ -6,27 +6,40 @@ objectName 	particlesDict;
 objectType 	dictionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
-
 setFields
 {
+	/*
+	    	Default value for fields defined for particles
+		These fields should always be defined for simulations with 
+		spherical particles.
+	*/
+
 	defaultValue 
 	{
-		velocity 		realx3 	(0 0 0); 		// linear velocity (m/s)
-		acceleration 	realx3 	(0 0 0); 		// linear acceleration (m/s2)
-		rVelocity 		realx3 	(0 0 0);  		// rotational velocity (rad/s)
-		shapeName 		word	lightSphere; 	// name of the particle shape 
+		velocity 		realx3 	    (0 0 0); 		  // linear velocity (m/s)
+
+		acceleration 	realx3 	    (0 0 0); 		  // linear acceleration (m/s2)
+
+		rVelocity 		realx3 	    (0 0 0);  		  // rotational velocity (rad/s)
+
+		shapeName 		word	lightSphere; 	      // name of the particle shape 
 	}
 
 	selectors
 	{}
 }
 
-// positions particles 
-positionParticles
+positionParticles                                          // positions particles 
 {
-	method empty;     		// creates the required fields with zero particles (empty). 
+	method                empty;          // other options: ordered and random 
 
-	maxNumberOfParticles 	50000; // maximum number of particles in the simulation
-	mortonSorting 			Yes;    // perform initial sorting based on morton code?   
+	regionType              box;          // other options: cylinder and sphere  
 
+	boxInfo                               // box region for positioning particles 
+	{
+		min (-0.08 -0.08 0.015);          // lower corner point of the box 
+
+		max ( 0.08  0.08 0.098);	      // upper corner point of the box 
+	} 
 }
+

tutorials/sphereGranFlow/layeredSiloFilling/settings/settingsDict

@@ -6,35 +6,33 @@ objectName 	settingsDict;
 objectType 	dictionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
+run   layerdSiloFilling;
 
-run 			layerdSiloFilling;
+dt 				0.00001; 	      // time step for integration (s)
 
-dt 				0.00001; 	// time step for integration (s)
+startTime 		    0.0; 		  // start time for simulation 
 
-startTime 		0; 			// start time for simulation 
+endTime 		    5.0;	      // end time for simulation 
 
-endTime 		5;	 		// end time for simulation 
+saveInterval 	   0.05; 		  // time interval for saving the simulation
 
-saveInterval 	0.05; 		// time interval for saving the simulation
+timePrecision         6;	      // maximum number of digits for time folder 
 
-timePrecision   6;			// maximum number of digits for time folder 
+g 			 (0 0 -9.8);          // gravity vector (m/s2) 
 
-g 				(0 0 -9.8); // gravity vector (m/s2) 
+// save data objects that are not automatically saved on disk.
+// overrides the default behavior 
+includeObjects (diameter);                                                 
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
 
-/*
-	Simulation domain 
-	every particles that goes outside this domain is deleted.
-*/
-domain 
-{
-	min (-0.1 -0.1 -0.1);
-	max ( 0.1  0.1  0.40);
-}
+integrationMethod 		AdamsBashforth2; 	      // integration method 
 
-integrationMethod 		AdamsBashforth3; 	// integration method 
+writeFormat				          ascii;          // data writting format (ascii or binary)
+
+timersReport 			            Yes;  	      // report timers
+
+timersReportInterval   	           0.01;	      // time interval for reporting timers
 
-writeFormat				ascii;
 
-timersReport 			Yes;  	// report timers?
 
-timersReportInterval   	0.01;	// time interval for reporting timers