Merge pull request #125 from PhasicFlow/develop

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)

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

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

src/Particles/Insertion/Insertion/Insertions.hpp

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

src/phasicFlow/repository/repository/repository.cpp

@@ -144,7 +144,8 @@ bool pFlow::repository::addToRepository(IOobject* io)
     if( !objects_.insertIf(io->name(), io ) )
 	{
 		warningInFunction<<
-    	"Failed to add object " << io->name() <<
+    	"Failed to add object " << io->name() << 
+		" with type "<< lookupObjectTypeName(io->name())<<
     	" to repository " << this->name() << 
     	". It is already in this repository. \n";
     	return false;

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/grainGranFlow/rotatingDrum/caseSetup/interaction

@@ -0,0 +1,56 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName interaction;
+objectType dicrionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+materials      (prop1);    // a list of materials names
+densities      (2500.0);   // density of materials [kg/m3]
+
+contactListType   sortedContactList; 
+
+model
+{
+   contactForceModel cGRelativeLinearLimited; // cGAbsoluteLinearLimited   or   cGRelativeLinearLimited
+   
+   rollingFrictionModel normal;
+   
+   grainDissipationModel  none; // simplified or KTGF or none 
+
+   Yeff  (1.0e6);       // Young modulus [Pa]
+
+   Geff  (0.8e6);       // Shear modulus [Pa]
+
+   nu    (0.25);        // Poisson's ratio [-]
+
+   en    (0.8);         // coefficient of normal restitution
+
+   et    (1.0);         // coefficient of tangential restitution 
+
+   mu    (0.3);         // dynamic friction 
+
+   mur   (0.1);         // rolling friction 
+   
+   kn    (1050);
+   
+   kt    (800);
+   
+   additionalDissipationModel GB;
+   
+   
+}
+
+contactSearch
+{
+   method         NBS;           // method for broad search particle-particle
+   
+   
+   updateInterval 10;
+   sizeRatio      1.1;
+   cellExtent     0.55;
+   adjustableBox  Yes;
+
+}

tutorials/grainGranFlow/rotatingDrum/caseSetup/particleInsertion

@@ -0,0 +1,14 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName particleInsertion;
+objectType dicrionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+active  no;          // is insertion active?
+
+collisionCheck No;   // not implemented for yes
+
+

tutorials/grainGranFlow/rotatingDrum/caseSetup/shapes

@@ -0,0 +1,13 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName 	grainDict;
+objectType 	grainShape;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+names 		(sphere1); 	// names of shapes 
+grainDiameters	(0.0008);	// diameter of shapes 
+sphereDiameters	(0.0004);	// diameter of shapes 
+materials	(prop1);	// material names for shapes 

tutorials/grainGranFlow/rotatingDrum/runThisCase

@@ -0,0 +1,21 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+echo "\n<--------------------------------------------------------------------->"
+echo "1) Creating particles"
+echo "<--------------------------------------------------------------------->\n"
+particlesPhasicFlow
+
+echo "\n<--------------------------------------------------------------------->"
+echo "2) Creating geometry"
+echo "<--------------------------------------------------------------------->\n"
+geometryPhasicFlow
+
+echo "\n<--------------------------------------------------------------------->"
+echo "3) Running the case"
+echo "<--------------------------------------------------------------------->\n"
+grainGranFlow
+
+
+
+
+#------------------------------------------------------------------------------

tutorials/grainGranFlow/rotatingDrum/settings/domainDict

@@ -0,0 +1,55 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName 	domainDict;
+objectType 	dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+
+globalBox 
+{
+	min (-0.037 -0.037 -0.001);
+	max (0.037   0.037 0.026);
+}
+
+boundaries
+{
+	timeControl runTime;
+	
+	neighborListUpdateInterval 0.0001;
+
+	neighborLength 0.004;
+
+
+	left
+	{
+		type exit;	
+	}
+
+	right
+	{
+		type exit;
+	}
+
+	bottom
+	{
+		type exit;
+	}
+
+	top
+	{
+		type exit;
+	}
+
+	rear
+	{
+		type exit;
+	}
+
+	front
+	{
+		type exit;
+	}
+}

tutorials/grainGranFlow/rotatingDrum/settings/geometryDict

@@ -0,0 +1,71 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName 	geometryDict;
+objectType 	dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+// motion model: rotating object around an axis
+motionModel rotatingAxis; 
+
+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.025);   // end point of cylinder axis
+		radius1 0.036;		// radius at p1
+		radius2 0.036;		// 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.036 -0.036 0.0);	// first point of the wall
+		p2 ( 0.036 -0.036 0.0);   // second point
+		p3 ( 0.036  0.036 0.0);   // third point
+		p4 (-0.036  0.036 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.036 -0.036 0.025);	// first point of the wall
+		p2 ( 0.036 -0.036 0.025);   // second point
+		p3 ( 0.036  0.036 0.025);   // third point
+		p4 (-0.036  0.036 0.025); 
+		material prop1;
+		motion rotAxis;
+	}
+		
+	
+}
+
+// information for rotatingAxisMotion motion model 
+rotatingAxisInfo
+{
+	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 0.8; 		// rotation speed (rad/s)
+	}
+}

tutorials/grainGranFlow/rotatingDrum/settings/particlesDict

@@ -0,0 +1,59 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+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	sphere1; // name of the particle shape 
+	}
+
+	selectors
+	{}
+}
+
+// positions particles 
+positionParticles
+{
+	method ordered;     		// ordered positioning
+
+	maxNumberOfParticles 	200000; // maximum number of particles in the simulation
+	mortonSorting 			Yes;    // perform initial sorting based on morton code?   
+	regionType box;
+	
+	boxInfo  // box for positioning particles 
+	{
+		min (-0.025 -0.025 0.001);    // lower corner point of the box 
+		max ( 0.025  0.025 0.024);	// upper corner point of the box 
+	}
+	
+	orderedInfo
+	{
+		// minimum space between centers of particles
+		diameter 0.001;
+		
+		// number of particles in the simulation 	 	
+		numPoints 50000;
+	
+		// axis order for filling the space with particles		 	
+		axisOrder (x y z);  
+	}
+	
+
+
+}

tutorials/grainGranFlow/rotatingDrum/settings/settingsDict

@@ -0,0 +1,37 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName 	settingsDict;
+objectType 	dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+run 			rotatingDrumMedium;
+
+dt 				0.000005; 	// time step for integration (s)
+
+startTime 		0.0; 			// start time for simulation 
+
+endTime 		3.0; 		// end time for simulation 
+
+saveInterval 	0.05; 		// time interval for saving the simulation
+
+timePrecision   6;			// maximum number of digits for time folder 
+
+g 				(0 -9.8 0); // gravity vector (m/s2) 
+
+/*
+	Simulation domain 
+	every particles that goes outside this domain is deleted.
+*/
+
+includeObjects (diameter);
+
+integrationMethod 		AdamsBashforth2; 	// integration method 
+
+writeFormat				ascii;
+
+timersReport 			Yes;  	// report timers?
+
+timersReportInterval   	0.01;	// time interval for reporting timers

utilities/CMakeLists.txt

@@ -9,5 +9,5 @@ add_subdirectory(pFlowToVTK)
 
 add_subdirectory(Utilities)
 
-#add_subdirectory(postprocessPhasicFlow)
+add_subdirectory(postprocessPhasicFlow)
 

utilities/Utilities/readControlDict.hpp

@@ -46,8 +46,8 @@ protected:
 
 	int32 precision_;
 
-	inline static fileSystem defaultRootPath = CWD();
-	inline static fileSystem defaultCDPath = CWD()/"system"+"controlDict";
+	inline static fileSystem defaultRootPath = pFlow::CWD();
+	inline static fileSystem defaultCDPath = pFlow::CWD()/word("system")+word("controlDict");
 	
 	word convertTimeToName(const real t, const int32 precision)const;
 

utilities/Utilities/readFromTimeFolder.hpp

@@ -113,7 +113,7 @@ public:
 	}
 
 	template<typename T>
-	uniquePtr<pointField_H<T>>  createUniformPointField_H(word const& name, T value)
+	uniquePtr<pointField_H<T>>  createUniformPointField_H(word const& name, T value, bool regRep = true)
 	{
 		if( !checkForPointStructure() )
 		{
@@ -124,21 +124,43 @@ public:
 
 		auto& pStruct = repository_.lookupObject<pointStructure>(pointStructureFile__);
 
-
-		auto Ptr = makeUnique<pointField_H<T>>
-		(
-			objectFile
+		if(regRep)
+		{
+			auto Ptr = makeUnique<pointField_H<T>>
 			(
-				name,
-				"",
-				objectFile::READ_NEVER,
-				objectFile::WRITE_NEVER
-			),
-			pStruct,
-			T{},
-			value
-		);
-		return Ptr;
+				objectFile
+				(
+					name,
+					"",
+					objectFile::READ_NEVER,
+					objectFile::WRITE_NEVER
+				),
+				pStruct,
+				T{},
+				value
+			);
+			return Ptr;
+		}
+		else
+		{
+			auto Ptr = makeUnique<pointField_H<T>>
+			(
+				objectFile
+				(
+					name,
+					"",
+					objectFile::READ_NEVER,
+					objectFile::WRITE_NEVER
+				),
+				nullptr,
+				pStruct,
+				value
+			);
+			Ptr().fill(value);
+			return Ptr;
+		}
+		
+		return nullptr;
 	}
 
 	template<typename T>

utilities/postprocessPhasicFlow/CMakeLists.txt

@@ -1,12 +1,15 @@
 
 set(source_files 
-postprocessPhasicFlow.cpp 
+cellMapper.cpp
+rectangleMesh.cpp
+countField.cpp
+countFields.cpp
 postprocess.cpp
 processField.cpp
 ProcessFields.cpp
 includeMask.cpp
 IncludeMasks.cpp
-
+postprocessPhasicFlow.cpp
 )
 set(link_lib phasicFlow Interaction Kokkos::kokkos Utilities)
 

utilities/postprocessPhasicFlow/IncludeMask.hpp

@@ -186,11 +186,12 @@ protected:
 		
 	Operator 				operator_;
 	
-	pointField_H<T> 		field_;
+	uniquePtr<pointField_H<T>> 		fieldPtr_;
 
+	hostViewType1D<T> 			field_;
 public:
 
-	TypeInfoTemplate2("IncludeMask", T, Operator);
+	TypeInfoTemplate12("IncludeMask", T, Operator);
 
 	IncludeMask(
 		const dictionary& dict,
@@ -199,8 +200,10 @@ public:
 	:
 		includeMask(dict, opType, timeFolder),
 		operator_(dict),
-		field_(timeFolder.readPointField_H<T>(this->fieldName()))
-	{}
+		fieldPtr_(timeFolder.readPointField_H<T>(this->fieldName())),
+		field_(fieldPtr_().hostView())
+	{
+	}
 
 	add_vCtor(
 		includeMask,
@@ -212,6 +215,11 @@ public:
 		return operator_(field_[n]);
 	}
 
+	uint32 size()const override
+	{
+		return field_.size();
+	}
+
 };
 
 
@@ -221,7 +229,7 @@ class IncludeMask<T,allOp<T>>
 	public includeMask
 {
 public:
-	TypeInfoTemplate2("IncludeMask", T, allOp<int8>);
+	TypeInfoTemplate12("IncludeMask", T, allOp<int8>);
 
 	IncludeMask(
 		const dictionary& dict,
@@ -240,6 +248,11 @@ public:
 	{
 		return true;
 	}
+
+	uint32 size()const override
+	{
+		return 0;
+	}
 };
 
 

utilities/postprocessPhasicFlow/ProcessField.hpp

@@ -40,15 +40,14 @@ class ProcessField
 
 protected:
 
-	pointField_H<T>& 	field_;
+	uniquePtr<pointField_H<T>> 	field_;
 
 
-	rectMeshField_H<T>& 	processedField_;
+	rectMeshField_H<T> 	processedField_;
 
 public:
 
-	TypeInfoTemplate("ProcessField", T);
-
+	TypeInfoTemplate11("ProcessField", T);
 
 	ProcessField(
 		const dictionary& dict,
@@ -58,24 +57,14 @@ public:
 		processField(dict, pToCell, rep),
 		field_(
 			this->isUniform()?
-			timeFolder().createUniformPointField_H(this->fieldName(), getUniformValue() ):
+			timeFolder().createUniformPointField_H(this->fieldName(), getUniformValue(), false ):
 			timeFolder().readPointField_H<T>(this->fieldName()) 
 			),
 		processedField_
 		(
-			processedRepository().emplaceObject<rectMeshField_H<T>>
-			(
-				objectFile
-				(
-					processedFieldName(),
-					"",
-					objectFile::READ_NEVER,
-					objectFile::WRITE_ALWAYS
-				),
-				mesh(),
-				processedFieldName(),
-				T{}
-			)
+			mesh(),
+			processedFieldName(),
+			T{}
 		)
 	{
 		
@@ -99,33 +88,59 @@ public:
 		
 		const includeMask& incMask = includeMask_();
 		
-		auto numerator = sumMaksOp( field_ , this->pointToCell(), incMask);
+		auto numeratorPtr =  sumMaksOp( field_() , this->pointToCell(), incMask);
+		uniquePtr<rectMeshField_H<real>> denomeratorPtr;
 		
-		rectMeshField_H<real> denomerator( this->mesh(), real{} );
-
 		if(operation() == "sum")
 		{
-			denomerator = rectMeshField_H<real>(this->mesh(), static_cast<real>(1.0));
+			denomeratorPtr = makeUnique<rectMeshField_H<real>>(this->mesh(), static_cast<real>(1.0));
 
 		}else if(operation() == "average")
 		{
 
-			pointField_H<real> oneFld(field_.pStruct(), static_cast<real>(1.0), static_cast<real>(1.0));
+			pointField_H<real> oneFld(
+				objectFile
+				(
+					"oneField",
+					"", 
+					objectFile::READ_NEVER, 
+					objectFile::WRITE_NEVER
+				),
+				const_cast<pointStructure&>(field_().pStruct()), 
+				static_cast<real>(1.0), 
+				static_cast<real>(1.0)
+			);
 
-			denomerator = sumOp(oneFld, this->pointToCell());
+			denomeratorPtr = sumOp(oneFld, this->pointToCell());
 
 		}else if(operation() == "averageMask")
 		{
-			pointField_H<real> oneFld(field_.pStruct(), static_cast<real>(1.0), static_cast<real>(1.0));
+			//pointField_H<real> oneFld(field_().pStruct(), static_cast<real>(1.0), static_cast<real>(1.0));
+			pointField_H<real> oneFld(
+				objectFile
+				(
+					"oneField",
+					"", 
+					objectFile::READ_NEVER, 
+					objectFile::WRITE_NEVER
+				),
+				const_cast<pointStructure&>(field_().pStruct()), 
+				static_cast<real>(1.0), 
+				static_cast<real>(1.0)
+			);
+			
+			denomeratorPtr = sumMaksOp(oneFld, this->pointToCell(), incMask);
 
-			denomerator = sumMaksOp(oneFld, this->pointToCell(), incMask);
-		}else
+		}
+		else
 		{
 			fatalErrorInFunction<<"operation is not known: "<< operation()<<endl;
 			fatalExit;
 		}
 
-		
+		auto& numerator = numeratorPtr();
+		auto& denomerator = denomeratorPtr();
+
 		for(int32 i=0; i<this->mesh().nx(); i++ )
 		{
 			for(int32 j=0; j<this->mesh().ny(); j++ )

utilities/postprocessPhasicFlow/cellMapper.cpp

@@ -0,0 +1,91 @@
+/*------------------------------- 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 "cellMapper.hpp"
+//#include "cellMapperKernels.hpp"
+#include "streams.hpp"
+
+
+void pFlow::cellMapper::allocateArrays(rangeU32 nextRng)
+{
+    checkAllocateNext(nextRng);
+    nullifyNext(nextRng);
+    reallocFill(head_, domainCells_.nx(), domainCells_.ny(), domainCells_.nz(), NoPos);
+}
+
+void pFlow::cellMapper::checkAllocateNext(rangeU32 nextRng)
+{
+
+    auto newCap = nextRng.end();
+    
+    if( nextCapacity_ < newCap)
+    {
+        nextCapacity_ = newCap;	
+        reallocNoInit(next_, nextCapacity_);
+    }
+}
+
+void pFlow::cellMapper::nullifyHead()
+{
+    fill(head_, NoPos);	
+}
+
+void pFlow::cellMapper::nullifyNext(rangeU32 nextRng)
+{
+	fill(next_, nextRng, NoPos);
+}
+
+pFlow::cellMapper::cellMapper(
+    const rectangleMesh& rectMesh,
+	const hostViewType1D<realx3>& pointPos,
+	const pFlagTypeHost& flags
+)
+: 
+   	domainCells_(rectMesh)
+{
+
+    allocateArrays(flags.activeRange());
+}
+
+bool pFlow::cellMapper::build
+(
+    const hostViewType1D<realx3>& pointPos, 
+    const pFlagTypeHost & flags
+)
+{
+    auto aRange = flags.activeRange();
+    
+    checkAllocateNext(aRange);
+    nullifyHead();
+    nullifyNext(aRange);
+	int32x3 ind;
+	for(uint32 i=aRange.start(); i<aRange.end(); i++)
+	{
+		if( domainCells_.isInsideIndex(pointPos[i], ind) )
+		{
+			next_[i] = head_(ind.x(), ind.y(), ind.z());
+			head_(ind.x(), ind.y(), ind.z()) = i;
+		}
+	}     
+    
+    return true;
+
+}
+

utilities/postprocessPhasicFlow/cellMapper.hpp

@@ -0,0 +1,137 @@
+/*------------------------------- 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 __cellMapper_hpp__
+#define __cellMapper_hpp__
+
+#include "phasicFlowKokkos.hpp"
+#include "pointFlag.hpp"
+#include "rectangleMesh.hpp"
+
+
+namespace pFlow
+{
+
+class cellMapper
+{
+public:
+
+	using HeadType 			= hostViewType3D<uint32>;
+
+	using NextType 			= hostViewType1D<uint32>;	
+
+	
+	static constexpr uint32 NoPos 	= 0xFFFFFFFF;
+
+	class CellIterator
+	{
+	private:
+		HeadType 		head_;
+
+		NextType 		next_;
+
+	public:
+
+		CellIterator(const HeadType& head, const NextType& next)
+		:
+			head_(head),
+			next_(next)
+		{}
+		
+		static constexpr uint32 NoPos = 0xFFFFFFFF;
+
+		INLINE_FUNCTION_H
+		int32x3 numCells()const {
+			return int32x3(head_.extent(0), head_.extent(1), head_.extent(2));}
+
+		INLINE_FUNCTION_H
+		uint32 start(int32 i, int32 j, int32 k)const {
+			return head_(i,j,k); }
+
+		INLINE_FUNCTION_H
+		uint32 getNext(uint32 n)const {
+			if(n == NoPos ) return NoPos;
+			return next_(n); }
+
+		INLINE_FUNCTION_H
+		uint32 next(uint32 n)const{
+			return next_(n);}
+	};
+
+private:
+	
+	const rectangleMesh&  	domainCells_;
+
+	HeadType 			head_{"NBS::head",1,1,1};
+
+	NextType 			next_{"NBS::next", 1};
+
+	uint32 				nextCapacity_ = 0;
+
+	void allocateArrays(rangeU32 nextRng);
+
+	void checkAllocateNext(rangeU32 nextRng);
+	
+	void nullifyHead();
+	
+	void nullifyNext(rangeU32 nextRng);
+	
+public:
+
+	TypeInfoNV("cellMapper");
+
+		
+		cellMapper(
+			const rectangleMesh& rectMesh,
+			const hostViewType1D<realx3>& pointPos,
+			const pFlagTypeHost& flags);
+					
+		INLINE_FUNCTION_HD
+		cellMapper(const cellMapper&) = default;
+
+		INLINE_FUNCTION_HD
+		cellMapper(cellMapper&&) = default;
+		
+		INLINE_FUNCTION_HD
+		cellMapper& operator = (const cellMapper&) = default;
+
+		INLINE_FUNCTION_HD
+		cellMapper& operator = (cellMapper&&) = default;
+
+		INLINE_FUNCTION_HD
+		~cellMapper()=default;
+
+	//// - Methods
+
+		auto getCellIterator()const
+		{
+			return CellIterator(head_, next_);
+		}
+
+		
+		bool build(
+			const hostViewType1D<realx3>& pointPos,
+			const pFlagTypeHost& flags);
+};
+
+} // pFlow
+
+#endif // __cellMapper_hpp__

utilities/postprocessPhasicFlow/countField.cpp

@@ -0,0 +1,94 @@
+/*------------------------------- 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 "countField.hpp"
+#include "repository.hpp"
+#include "twoPartEntry.hpp"
+
+bool pFlow::countField::getFieldType(
+	const dictionary& dict,
+	readFromTimeFolder& timeFolder,
+	word& fieldName,
+	word& fieldType)
+{
+	if(dict.containsDataEntry("field"))
+	{
+		const dataEntry& entry = dict.dataEntryRef("field");
+
+		if( isTwoPartEntry(entry))
+		{
+			twoPartEntry tpEntry(entry);
+			fieldName = "uniformField";
+			fieldType = tpEntry.firstPart();
+		}
+		else
+		{
+			fieldName = dict.getVal<word>("field");
+			if( !timeFolder.pointFieldFileGetType(fieldName, fieldType) )
+			{
+				fatalErrorInFunction<<"error in reading field type from file "<< fieldName<<
+				"in folder "<< timeFolder.path()<<endl;
+				return false;
+			}
+		}
+	}
+	else
+	{
+		fatalErrorInFunction<< "dictionary "<< dict.globalName()<<
+		"does not contain field keyword"<<endl;
+		return false;
+	}
+
+	return true;
+}
+
+
+pFlow::countField::countField(const dictionary& dict, repository& rep)
+:
+	dict_(dict),
+	timeFolder_(rep)
+{
+
+	word includeMaskType = dict_.getVal<word>("includeMask");
+
+	auto& incDict = dict_.subDictOrCreate(includeMaskType+"Info");
+	
+	includeMask_ = includeMask::create(incDict, includeMaskType, timeFolder_);
+
+}
+
+
+bool pFlow::countField::process(uint32& countedValue)
+{
+	auto& incMask = includeMask_();
+
+	countedValue = 0;
+	uint32 n = incMask.size();
+
+	for(uint32 i=0; i<n; i++)
+	{
+		if( incMask(i) )
+		{
+			countedValue++;
+		}
+	}
+
+	return true;
+}

utilities/postprocessPhasicFlow/countField.hpp

@@ -0,0 +1,102 @@
+/*------------------------------- 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 __countField_hpp__
+#define __countField_hpp__
+
+
+#include "virtualConstructor.hpp"
+#include "dictionary.hpp"
+#include "readFromTimeFolder.hpp"
+#include "includeMask.hpp"
+
+
+namespace pFlow
+{
+
+
+class repository;
+
+class countField
+{
+protected:
+
+	dictionary 	dict_;
+	
+	mutable readFromTimeFolder 	timeFolder_;
+
+	uniquePtr<includeMask>  	includeMask_ = nullptr;
+	
+	bool static getFieldType(
+			const dictionary& dict,
+			readFromTimeFolder& timeFolder,
+			word& fieldName,
+			word& fieldType);
+
+public:
+
+	TypeInfo("countField");
+
+	countField(const dictionary& dict, repository& rep);
+
+	
+	auto& dict()
+	{
+		return dict_;
+	}
+
+	const auto& dict()const
+	{
+		return dict_;
+	}
+
+	auto& timeFolderRepository()
+	{
+		return timeFolder_.db();
+	}
+	
+	
+	auto& timeFolder()
+	{
+		return timeFolder_;
+	}
+
+	word variableName()const
+	{
+		return dict_.name();
+	}
+
+	// requires a class to read pointField from timefolder 
+	bool process(uint32& countedValue);
+
+
+
+	//virtual bool writeToFile(iOstream& is) const = 0;
+
+	/*static 
+	uniquePtr<countField> create(
+		const dictionary& dict,
+		repository& rep);*/
+};
+
+
+}
+
+
+#endif //__countField_hpp__

utilities/postprocessPhasicFlow/countFields.cpp

@@ -0,0 +1,53 @@
+/*------------------------------- 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 "countFields.hpp"
+#include "countField.hpp"
+#include "repository.hpp"
+#include "includeMask.hpp"
+
+
+pFlow::countFields::countFields(const dictionary& dict)
+:
+	dict_(dict),
+	variableNames_(dict.dictionaryKeywords()),
+	countedValues_(variableNames_.size(), 0u)
+{
+
+}
+
+
+bool pFlow::countFields::process(repository& rep)
+{
+	uint32List counted;
+
+	for(const auto& name: variableNames_)
+	{
+		const dictionary& varDict = dict_.subDict(name);
+		uint32 count;
+		countField cField(varDict, rep);
+		cField.process(count);
+		counted.push_back(count);
+	}
+
+	countedValues_ = counted;
+
+	return true;
+}

utilities/postprocessPhasicFlow/countFields.hpp

@@ -0,0 +1,77 @@
+/*------------------------------- 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 __countFields_hpp__
+#define __countFields_hpp__
+
+
+#include "dictionary.hpp"
+
+
+namespace pFlow
+{
+
+
+class repository;
+
+class countFields
+{
+protected:
+
+	dictionary 	dict_;
+	
+	wordList 	variableNames_;
+
+	uint32List 	countedValues_;
+
+public:
+
+	TypeInfo("countFields");
+
+	countFields(const dictionary& dict);
+
+	auto& dict()
+	{
+		return dict_;
+	}
+
+	const auto& dict()const
+	{
+		return dict_;
+	}
+
+	const wordList& variableNames()const
+	{
+		return variableNames_;
+	}
+	const uint32List& countedValues()const
+	{
+		return countedValues_;
+	}
+	
+	// requires a class to read pointField from timefolder 
+	bool process(repository& rep);
+
+};
+
+
+}
+
+
+#endif //__countFields_hpp__

utilities/postprocessPhasicFlow/fieldOperations.hpp

@@ -31,25 +31,26 @@ namespace pFlow
 
 
 template<typename T>
-rectMeshField_H<T> sumOp( const pointField_H<T> field, const pointRectCell& pointToCell)
+uniquePtr<rectMeshField_H<T>> sumOp( pointField_H<T>& field, pointRectCell& pointToCell)
 {
 	// create field
-	const auto& mesh = pointToCell.mesh();
+	auto& mesh = pointToCell.mesh();
 	auto iterator = pointToCell.getCellIterator();
+	auto f = field.deviceView();
 
-	rectMeshField_H<T> results(mesh, T(0));
-
+	auto resultsPtr = makeUnique<rectMeshField_H<T>>(mesh, T(0));
+	auto& results = resultsPtr();
 	for(int32 i=0; i<mesh.nx(); i++)
 	{
 		for(int32 j=0; j<mesh.ny(); j++)
 		{
 			for(int32 k=0; k<mesh.nz(); k++)
 			{
-				auto n = iterator.start(i,j,k);	
+				uint32 n = iterator.start(i,j,k);	
 				T res (0);
-				while(n>-1)
+				while(n != cellMapper::NoPos)
 				{
-					res += field[n];
+					res += f[n];
 					n = iterator.getNext(n);
 				}
 
@@ -58,17 +59,19 @@ rectMeshField_H<T> sumOp( const pointField_H<T> field, const pointRectCell& poin
 		}
 	}
 
-	return results;
+	return resultsPtr;
 }
 
 template<typename T, typename incMask>
-rectMeshField_H<T> sumMaksOp( const pointField_H<T> field, const pointRectCell& pointToCell, const incMask& mask)
+uniquePtr<rectMeshField_H<T>> sumMaksOp( pointField_H<T>& field, pointRectCell& pointToCell, const incMask& mask)
 {
 	// create field
-	const auto& mesh = pointToCell.mesh();
+	auto& mesh = pointToCell.mesh();
 	auto iterator = pointToCell.getCellIterator();
+	auto f = field.deviceView();
 
-	rectMeshField_H<T> results(mesh, T(0));
+	auto resultsPtr = makeUnique<rectMeshField_H<T>>(mesh, T(0));
+	auto& results = resultsPtr();
 
 	for(int32 i=0; i<mesh.nx(); i++)
 	{
@@ -77,15 +80,15 @@ rectMeshField_H<T> sumMaksOp( const pointField_H<T> field, const pointRectCell&
 			for(int32 k=0; k<mesh.nz(); k++)
 			{
 				//auto [loop, n] = pointToCell.startLoop(i,j,k);
-				auto n = iterator.start(i,j,k);	
+				uint32 n = iterator.start(i,j,k);	
 				T res (0);
 				
-				while(n>-1)
+				while(n!= cellMapper::NoPos)
 				{
 					
 					if(mask(n))
 					{
-						res += field[n];
+						res += f[n];
 					}
 
 					n = iterator.getNext(n);
@@ -96,7 +99,7 @@ rectMeshField_H<T> sumMaksOp( const pointField_H<T> field, const pointRectCell&
 		}
 	}
 
-	return results;
+	return resultsPtr;
 }
 
 

utilities/postprocessPhasicFlow/includeMask.cpp

@@ -85,7 +85,7 @@ pFlow::uniquePtr<pFlow::includeMask> pFlow::includeMask::create(
 		auto objPtr = 
 			dictionaryvCtorSelector_[method]
 			(dict, opType, timeFolder);
-		REPORT(2)<< dict.name()<< " with model "<<greenText(method)<<" is processed."<<endREPORT;
+		REPORT(2)<< dict.name()<< " with model "<<Green_Text(method)<<" is processed."<<END_REPORT;
 		return objPtr;
 	}
 	else

utilities/postprocessPhasicFlow/includeMask.hpp

@@ -85,6 +85,8 @@ public:
 
 	virtual bool isIncluded(int32 n) const = 0;
 
+	virtual uint32 size()const = 0;
+
 	bool operator()(int32 n) const 
 	{
 		return isIncluded(n);

utilities/postprocessPhasicFlow/pointRectCell.hpp

@@ -21,7 +21,7 @@ Licence:
 #ifndef __pointRectCell_hpp__
 #define __pointRectCell_hpp__
 
-#include "mapperNBS.hpp"
+#include "cellMapper.hpp"
 #include "rectMeshFields.hpp"
 #include "pointStructure.hpp"
 
@@ -33,20 +33,19 @@ class pointRectCell
 {
 public:
 
-	using mapType = mapperNBS<DefaultHostExecutionSpace>;
+	using mapType = cellMapper;
 
-	using memory_space = typename mapType::memory_space;
+	//using memory_space = typename mapType::memory_space;
 
 protected:
 
 	repository&				processedRepository_;
 
-	rectangleMesh& 			mesh_;
-
+	rectangleMesh 			mesh_;
 
 	const pointStructure& 	pStruct_;
 
-	ViewType1D<realx3, memory_space> 		pointPosition_;
+	hostViewType1D<realx3> 	pointPosition_;
 
 	mapType 				map_;
 
@@ -61,33 +60,29 @@ public:
 	:
 		processedRepository_(rep),
 		mesh_
-		( 
-			processedRepository_.emplaceObject<rectangleMesh>
-			(
-				objectFile
-				(
-					"rectMesh",
-					"",
-					objectFile::READ_NEVER,
-					objectFile::WRITE_NEVER
-				),
-				dictMesh
-			)
+		(
+			dictMesh,
+			&rep
 		),
 		pStruct_(pStruct),
-		pointPosition_(pStruct_.pointPosition().hostVectorAll()),
-		map_( 
-			mesh_.domain(), 
-			mesh_.nx(),
-			mesh_.ny(),
-			mesh_.nz(),
-			pointPosition_),
-		nPointInCell_(mesh_, 0)
+		pointPosition_(pStruct_.pointPosition().hostViewAll()),
+		map_
+		( 
+			mesh_,
+			pointPosition_,
+			pStruct_.activePointsMaskHost()
+		),
+		nPointInCell_(mesh_,"nPointInCell", 0)
 	{
 
 		mapPOints();
 	}
 
+	auto& mesh()
+	{
+		return mesh_;
+	}
+
 	const auto& mesh()const
 	{
 		return mesh_;
@@ -100,25 +95,23 @@ public:
 
 	void mapPOints()
 	{
-		range activeRange = pStruct_.activeRange();
-		auto activeMask  = pStruct_.activePointsMaskH();
+		auto points = pStruct_.pointPositionHost();
+		auto activeMask  = pStruct_.activePointsMaskHost();
 
-		
-		map_.buildCheckInDomain(activeRange, activeMask);
+		map_.build(points, activeMask);
 		
 		auto iterator = map_.getCellIterator();
-		
 	
-		for(int32 i=0; i<map_.nx(); i++)
+		for(int32 i=0; i<mesh_.nx(); i++)
 		{
-			for(int32 j=0; j<map_.ny(); j++)
+			for(int32 j=0; j<mesh_.ny(); j++)
 			{
-				for(int32 k=0; k<map_.nz(); k++)
+				for(int32 k=0; k<mesh_.nz(); k++)
 				{
 
 					int32 res = 0;
-					int32 n = iterator.start(i,j,k);
-					while( n>-1)
+					uint32 n = iterator.start(i,j,k);
+					while( n!= cellMapper::NoPos)
 					{
 						res+=1;
 						n = iterator.getNext(n);

utilities/postprocessPhasicFlow/postprocess.cpp

@@ -21,6 +21,7 @@ Licence:
 #include "postprocess.hpp" 
 #include "timeFolder.hpp"
 #include "pointStructure.hpp"
+#include "countFields.hpp"
 #include "vocabs.hpp"
 #include "vtkFile.hpp"
 
@@ -29,22 +30,26 @@ pFlow::postprocess::postprocess(systemControl& control)
 	control_(control),
 	dict_(postprocessFile__, control_.settings().path()+postprocessFile__)
 {
-	REPORT(1)<<"Reading numberBased dictionary ..."<<endREPORT;
+	REPORT(1)<<"Reading numberBased dictionary ..."<<END_REPORT;
 	auto nbDict = dict_.subDictOrCreate("numberBased");
 
 	numberBasedDictNames_ = dict_.subDictOrCreate("numberBased").dictionaryKeywords();
 	if(!numberBasedDictNames_.empty())
 	{
-		REPORT(1)<< "numberBased dictionary contains " << yellowText(numberBasedDictNames_)<<endREPORT<<endl;	
+		REPORT(1)<< "numberBased dictionary contains " << Yellow_Text(numberBasedDictNames_)<<END_REPORT<<endl;	
 	}
 	
-
 	weightBasedDictNames_ = dict_.subDictOrCreate("weightBased").dictionaryKeywords();
 	if(!weightBasedDictNames_.empty())
 	{
-		REPORT(1)<< "numberBased dictionary contains " << yellowText(weightBasedDictNames_)<<endREPORT<<endl;	
+		REPORT(1)<< "numberBased dictionary contains " << Yellow_Text(weightBasedDictNames_)<<END_REPORT<<endl;	
+	}
+
+	countDictNames_ = dict_.subDictOrCreate("counting").dictionaryKeywords();
+	if(!countDictNames_.empty())
+	{
+		REPORT(1)<< "counting dictionary contains " << Yellow_Text(countDictNames_)<<END_REPORT<<endl;	
 	}
-	
 
 }
 
@@ -55,43 +60,32 @@ bool pFlow::postprocess::processTimeFolder(real time, const word& tName, const f
 	
 	time_ = time;
 
-	REPORT(0)<<"Working on time folder "<< cyanText(time)<<endREPORT; 
-	timeFolderReposiory_ = 
-		makeUnique<repository>
-		(
-			"timeFolder-"+tName,
-			localFolder,
-			&control_
-		);
+	REPORT(0)<<"Working on time folder "<< Cyan_Text(time)<<END_REPORT; 
 	
-	REPORT(1)<<"Reading pointStructure"<<endREPORT;
-	timeFolderReposiory().emplaceObject<pointStructure>(
-		objectFile
-		(
-			pFlow::pointStructureFile__,
-			"",
-			objectFile::READ_ALWAYS,
-			objectFile::WRITE_NEVER
-		));	
+	control_.time().setTime(time);
 
+	REPORT(1)<<"Reading pointStructure"<<END_REPORT;
+	pointStructure pStruct(control_);
+
+	// first delete the object to remove fields from repository 
+	pointToCell_.reset(nullptr);
 	
-	REPORT(1)<<"Creating mesh and point to cell mapper"<<endREPORT;
+	REPORT(1)<<"Creating mesh and point to cell mapper"<<END_REPORT;
 	pointToCell_ = makeUnique<pointRectCell>(
 		dict_.subDict("rectMesh"),
-		timeFolderReposiory().lookupObject<pointStructure>(pointStructureFile__),
-		timeFolderReposiory());
+		pStruct,
+		control_.time());
 
 	// first numberbased dict 
+	
 	processedFields_.clear();
 	for(word& dictName:numberBasedDictNames_)
 	{
-		
-		
 		auto fieldDict  = dict_.subDict("numberBased").subDict(dictName);
 		auto ppFieldPtr = processField::create(
 			fieldDict,
 			pointToCell_(),
-			timeFolderReposiory());
+			control_.time());
 
 		if(!ppFieldPtr->process())
 		{
@@ -99,11 +93,28 @@ bool pFlow::postprocess::processTimeFolder(real time, const word& tName, const f
 		}
 
 		processedFields_.push_back( ppFieldPtr.release() );
-
 		output<<endl;
-
 	}
 
+	countedVariableNamesList_.clear();
+	countedVairablesLists_.clear();
+
+	for(const auto& countDictName:countDictNames_)
+	{
+		REPORT(1)<<"Processing "<< Yellow_Text("counting."<<countDictName)<<END_REPORT;
+		const dictionary& countDict = dict_.subDict("counting").subDict(countDictName);
+
+		countFields cFields(countDict);
+
+		cFields.process(timeFolderReposiory());
+
+		countedVariableNamesList_.push_back(
+			cFields.variableNames());
+
+		countedVairablesLists_.push_back(cFields.countedValues());
+
+	}
+	output<<"\n";
 
 
 	return true;
@@ -125,7 +136,7 @@ bool pFlow::postprocess::writeToVTK(fileSystem destPath, word bName)const
 
 	if(!vtk) return false;
 
-	REPORT(1)<<"Writing processed fields to vtk file..."<<endREPORT;
+	REPORT(1)<<"Writing processed fields to vtk file..."<<END_REPORT;
 	// mesh
 	pointToCell_->mesh().writeToVtk(vtk());
 

utilities/postprocessPhasicFlow/postprocess.hpp

@@ -46,6 +46,12 @@ protected:
 
 	wordList 				weightBasedDictNames_;
 
+	wordList 				countDictNames_;
+
+	List<wordList>  		countedVariableNamesList_;
+
+	List<uint32List> 		countedVairablesLists_;
+
 	uniquePtr<repository> 	timeFolderReposiory_ {nullptr};
 
 	uniquePtr<pointRectCell>	pointToCell_ {nullptr};

utilities/postprocessPhasicFlow/postprocessPhasicFlow.cpp

@@ -74,7 +74,7 @@ int main(int argc, char** argv )
 	#include "initialize_Control.hpp"
 
 	
-	pFlow::postprocess post(Control);
+	
 		
 	// time folders in case 
 	timeFolder folders(Control);
@@ -93,12 +93,16 @@ int main(int argc, char** argv )
 
 	pFlow::fileSystem destFolder = pFlow::fileSystem(outFolder);
 
+	pFlow::postprocess post(Control);
+
 	do 
 	{
 
-
-		if( !validRange.isMember( folders.time() ) )continue;
-
+		if( !validRange.isMember( folders.time() ) )
+		{
+			continue;
+		}
+				
 		if( !withZeroFolder && pFlow::equal(folders.time() , 0.0))continue;
 
 		post.processTimeFolder(folders);
@@ -106,9 +110,10 @@ int main(int argc, char** argv )
 		if(!post.writeToVTK(destFolder, "processed"))
 		{
 			fatalExit;
-		}	
-		
+		}
 
+		
+			
 	}while (folders++);
 
 	#include "finalize.hpp"

utilities/postprocessPhasicFlow/processField.cpp

@@ -35,7 +35,7 @@ pFlow::processField::processField(
 	processedFieldName_(dict.name()),
 	operation_(dict.getVal<word>("operation")),
 	includeMaskType_(dict.getVal<word>("includeMask")),
-	threshold_(dict.getValOrSet<int32>("threshold", 1))
+	threshold_(dict.getValOrSetMax<int32>("threshold", 1))
 {
 	
 	if(!processField::getFieldType(
@@ -50,8 +50,7 @@ pFlow::processField::processField(
 	auto& incDict = dict_.subDictOrCreate(includeMaskType_+"Info");
 	
 	includeMask_ = includeMask::create(incDict, includeMaskType_, timeFolder_);
-	
-	
+
 }
 
 bool pFlow::processField::getFieldType(
@@ -107,7 +106,6 @@ pFlow::processField::create(
 		return nullptr;
 	}
 	
-
 	auto method = angleBracketsNames("ProcessField", fType);
 
 	if( dictionaryvCtorSelector_.search(method) )
@@ -115,7 +113,7 @@ pFlow::processField::create(
 		auto objPtr = 
 			dictionaryvCtorSelector_[method]
 			(dict, pToCell, rep);
-		REPORT(2)<<"Processing/creating " << yellowText(dict.name())<< " with model "<<greenText(method)<<"."<<endREPORT;
+		REPORT(2)<<"Processing/creating " << Yellow_Text(dict.name())<< " with model "<<Green_Text(method)<<"."<<END_REPORT;
 		return objPtr;
 	}
 	else

utilities/postprocessPhasicFlow/processField.hpp

@@ -70,6 +70,7 @@ public:
 
 	processField(const dictionary& dict, pointRectCell& pToCell, repository& rep);
 
+	virtual ~processField() = default;
 
 	create_vCtor(
 			processField,
@@ -85,11 +86,21 @@ public:
 		return pointToCell_.mesh();
 	}
 
+	auto& mesh()
+	{
+		return pointToCell_.mesh();
+	}
+
 	const auto& pointToCell()const
 	{
 		return pointToCell_;
 	}
 
+	auto& pointToCell()
+	{
+		return pointToCell_;
+	}
+
 	auto& dict()
 	{
 		return dict_;

utilities/postprocessPhasicFlow/rectMeshField.hpp

@@ -22,17 +22,19 @@ Licence:
 #define __rectMeshField_hpp__
 
 #include "rectangleMesh.hpp"
-#include "baseAlgorithms.hpp"
+#include "ViewAlgorithms.hpp"
 
 namespace pFlow
 {
 
-template<typename T, typename MemorySpace=void>
+template<typename T>
 class rectMeshField
+:
+	public IOobject
 {
 public:
 
-	using viewType = ViewType3D<T,MemorySpace>;
+	using viewType = ViewType3D<T,HostSpace>;
 
 	using memory_space = typename viewType::memory_space;
 
@@ -54,22 +56,49 @@ protected:
 public:
 
 
-	TypeInfoTemplateNV2("rectMeshField", T, memoerySpaceName());
+	TypeInfoTemplateNV111("rectMeshField", T, memoerySpaceName());
 	
-	rectMeshField(const rectangleMesh& mesh, const word& name ,const T& defVal)
+	rectMeshField(rectangleMesh& mesh, const word& name ,const T& defVal)
 	:
+		IOobject(
+			objectFile
+			(
+				name,
+				"",
+				objectFile::READ_NEVER,
+				objectFile::WRITE_NEVER
+			),
+			IOPattern::MasterProcessorOnly,
+			nullptr
+		),
 		mesh_(&mesh),
 		name_(name),
-		field_("pFlow::reactMeshField", mesh_->nx(), mesh_->ny(), mesh_->nz()),
+		field_("reactMeshField."+name, mesh_->nx(), mesh_->ny(), mesh_->nz()),
 		defaultValue_(defVal)
 	{
 		this->fill(defaultValue_);
 	}
 
-	rectMeshField(const rectangleMesh& mesh, const T& defVal)
+	rectMeshField(rectangleMesh& mesh, const T& defVal)
 	:
-		rectMeshField(mesh, "noName", defVal)
-	{}
+		IOobject(
+			objectFile
+			(
+				"noName",
+				"",
+				objectFile::READ_NEVER,
+				objectFile::WRITE_NEVER
+			),
+			IOPattern::MasterProcessorOnly,
+			nullptr
+		),
+		mesh_(&mesh),
+		name_("noName"),
+		field_("reactMeshField.noName", mesh_->nx(), mesh_->ny(), mesh_->nz()),
+		defaultValue_(defVal)
+	{
+		this->fill(defaultValue_);
+	}
 
 	rectMeshField(const rectMeshField&) = default;
 
@@ -144,13 +173,23 @@ public:
 	{
 		pFlow::fill(
 			field_,
-			range(0,mesh_->nx()),
-			range(0,mesh_->ny()),
-			range(0,mesh_->nz()),
 			val
 			);
 	}
 
+	 
+	bool write(iOstream& is, const IOPattern& iop)const override
+	{
+		notImplementedFunction;
+		return true;
+	}
+
+	bool read(iIstream& is, const IOPattern& iop) override
+	{
+		notImplementedFunction;
+		return true;
+	}
+
 	bool read(iIstream& is)
 	{
 		notImplementedFunction;

utilities/postprocessPhasicFlow/rectMeshFieldToVTK.hpp

@@ -26,7 +26,7 @@ namespace pFlow
 {
 
 template<typename T>
-bool convertRectMeshField(iOstream& os, rectMeshField_H<T>& field)
+bool convertRectMeshField(iOstream& os, const rectMeshField_H<T>& field)
 {
 	fatalErrorInFunction<< "this type is not supported "<<
 	field.typeName()<<endl;
@@ -36,7 +36,7 @@ bool convertRectMeshField(iOstream& os, rectMeshField_H<T>& field)
 
 
 template<>
-bool convertRectMeshField(iOstream& os, rectMeshField_H<real>& field)
+bool convertRectMeshField(iOstream& os, const rectMeshField_H<real>& field)
 {
 	
 	os<<"FIELD FieldData 1 " << field.name() << " 1 "<< field.size() << " float\n";
@@ -55,7 +55,7 @@ bool convertRectMeshField(iOstream& os, rectMeshField_H<real>& field)
 }
 
 template<>
-bool convertRectMeshField(iOstream& os, rectMeshField_H<realx3>& field)
+bool convertRectMeshField(iOstream& os, const rectMeshField_H<realx3>& field)
 {
 	
 	os<<"FIELD FieldData 1 " << field.name() << " 3 "<< field.size() << " float\n";
@@ -74,7 +74,7 @@ bool convertRectMeshField(iOstream& os, rectMeshField_H<realx3>& field)
 }
 
 template<>
-bool convertRectMeshField(iOstream& os, rectMeshField_H<int32>& field)
+bool convertRectMeshField(iOstream& os, const rectMeshField_H<int32>& field)
 {
 	
 	os<<"FIELD FieldData 1 " << field.name() << " 1 "<< field.size() << " int\n";

utilities/postprocessPhasicFlow/rectMeshFields.hpp

@@ -27,17 +27,17 @@ namespace pFlow
 {
 
 template<typename T>
-using rectMeshField_H 		= rectMeshField<T,HostSpace>;
+using rectMeshField_H 		= rectMeshField<T>;
 
-using int8RectMeshField_H  = rectMeshField<int8, HostSpace>;
+using int8RectMeshField_H  = rectMeshField<int8>;
 
-using int32RectMeshField_H 	= rectMeshField<int32, HostSpace>;
+using int32RectMeshField_H 	= rectMeshField<int32>;
 
-using int64RectMeshField_H  = rectMeshField<int64, HostSpace>;
+using int64RectMeshField_H  = rectMeshField<int64>;
 
-using realRectMeshField_H	= rectMeshField<real, HostSpace>;
+using realRectMeshField_H	= rectMeshField<real>;
 
-using realx3RectMeshField_H	= rectMeshField<realx3, HostSpace>;
+using realx3RectMeshField_H	= rectMeshField<realx3>;
 
 }
 

utilities/postprocessPhasicFlow/rectangleMesh.cpp

@@ -0,0 +1,52 @@
+#include "rectangleMesh.hpp"
+
+
+pFlow::rectangleMesh::rectangleMesh
+(
+    const box& mshBox, 
+    int32 nx, 
+    int32 ny, 
+    int32 nz,
+    repository* rep
+)
+:
+    IOobject(
+        objectFile
+        (
+            "rectMesh",
+            "",
+            objectFile::READ_NEVER,
+            objectFile::WRITE_NEVER
+        ),
+        IOPattern::MasterProcessorOnly,
+        rep
+    ),
+    meshBox_(mshBox),
+    numCells_(nx, ny, nz)
+{
+    if(mshBox.minPoint()>= mshBox.maxPoint())
+    {
+        fatalErrorInFunction<<"Lower corner point of mesh "<<mshBox.minPoint()<<
+        " confilicts with upper corner point of mesh "<<mshBox.maxPoint()<<endl;
+        fatalExit;
+    }
+
+    numCells_ = max( numCells_ , int32x3(1) );
+
+    dx_ = (mshBox.maxPoint() - mshBox.minPoint())/
+        realx3(numCells_.x_, numCells_.y_, numCells_.z_);
+
+}
+
+pFlow::rectangleMesh::rectangleMesh(const dictionary &dict, repository* rep)
+:
+    rectangleMesh(
+        box(dict),
+        dict.getVal<int32>("nx"),
+		dict.getVal<int32>("ny"),
+		dict.getVal<int32>("nz"),
+        rep
+    )
+{
+
+}

utilities/postprocessPhasicFlow/rectangleMesh.hpp

@@ -21,100 +21,112 @@ Licence:
 #ifndef __rectangleMesh_hpp__
 #define __rectangleMesh_hpp__
 
-#include "cells.hpp"
+#include "types.hpp"
+#include "error.hpp"
+#include "box.hpp"
+#include "IOobject.hpp"
+
+class repository;
 
 namespace pFlow
 {
 
-
-
 class rectangleMesh
 :
-	public cells<int32>
+	public IOobject
 {
-	
+private:
+
+	box 		meshBox_;
+
+	int32x3 	numCells_;
+
+	realx3 		dx_;
+
 public:
 
-	TypeInfoNV("rectangleMesh");
-	
+	TypeInfo("rectangleMesh");
 
-	INLINE_FUNCTION_HD
-	rectangleMesh(){};
-
-	INLINE_FUNCTION_HD
 	rectangleMesh(
-		const realx3& minP,
-		const realx3& maxP,
+		const box& mshBox,
 		int32 nx,
 		int32 ny,
-		int32 nz)
-	:
-	cells(
-		box(minP, maxP),
-		nx, ny, nz)
-	{}
+		int32 nz,
+		repository* rep);
+	
+	rectangleMesh(const dictionary & dict, repository* rep);
 
-	INLINE_FUNCTION_H
-	rectangleMesh(const dictionary & dict)
-	:
-	cells(
-		box(
-			dict.getVal<realx3>("min"),
-			dict.getVal<realx3>("max")),
-		dict.getVal<int32>("nx"),
-		dict.getVal<int32>("ny"),
-		dict.getVal<int32>("nz")
-		)
-	{}
-
-	INLINE_FUNCTION_HD
+	
 	rectangleMesh(const rectangleMesh&) = default;
 
-	INLINE_FUNCTION_HD
+	
 	rectangleMesh& operator = (const rectangleMesh&) = default;
 
-	INLINE_FUNCTION_HD
 	rectangleMesh(rectangleMesh&&) = default;
 
-	INLINE_FUNCTION_HD
+	
 	rectangleMesh& operator = (rectangleMesh&&) = default;
 
-	INLINE_FUNCTION_HD
-	~rectangleMesh() = default;
+	~rectangleMesh() override = default;
 
-
-	INLINE_FUNCTION_HD
 	int64 size()const
 	{
-		return this->totalCells();
+		return numCells_.x_*numCells_.y_*numCells_.z_;
+	}
+
+	int32 nx()const
+	{
+		return numCells_.x();
+	}
+
+	int32 ny()const
+	{
+		return numCells_.y();
+	}
+
+	int32 nz()const
+	{
+		return numCells_.z();
 	}
 
-	INLINE_FUNCTION_HD
 	real cellVol()const
 	{
-		auto [dx,dy,dz] = this->cellSize();
-		return dx*dy*dz;
+		return dx_.x_*dx_.y_*dx_.z_;
+	}
+	
+	realx3 minPoint()const
+	{
+		return meshBox_.minPoint();
+	}
+	
+	realx3 maxPoint()const
+	{
+		return meshBox_.maxPoint();
 	}
 
-	INLINE_FUNCTION_H
-	auto minPoint()const
+	inline
+	bool isInsideIndex(const realx3 p, int32x3 & ind )const
 	{
-		return domain().minPoint();
+		if(meshBox_.isInside(p))
+		{
+			ind = (p - meshBox_.minPoint())/dx_ ;
+			return true;
+		}
+		else
+		{
+			return false;
+		}
 	}
 
-	INLINE_FUNCTION_H
-	auto maxPoint()const
-	{
-		return domain().maxPoint();
-	}
-
-	bool read(iIstream& is)
+	bool write(iOstream& is, const IOPattern& iop)const override
 	{
+		notImplementedFunction;
 		return true;
 	}
 
-	bool write(iOstream& os)const
+	bool read(iIstream& is, const IOPattern& iop) override
 	{
+		notImplementedFunction;
 		return true;
 	}
 
@@ -124,7 +136,7 @@ public:
 		os<<"DIMENSIONS "<<nx()+1<<" "<< ny()+1 << " "<< nz()+1 <<endl;
 		
 		auto [x, y , z] = this->minPoint();
-		auto [dx, dy, dz] = this->cellSize();
+		auto [dx, dy, dz] = dx_;
 
 		os<<"X_COORDINATES "<< nx()+1 <<" float\n";
 		for(int32 i=0; i<nx()+1; i++)