Boundary conditions are created and tested.

- exit and periodic, follows the previous commit.

src/Interaction/CMakeLists.txt

@@ -5,11 +5,21 @@ contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp
 contactSearch/methods/cellBased/NBS/NBSLevel0.cpp
 contactSearch/methods/cellBased/NBS/NBS.cpp
 contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
+
+contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.cpp
+contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp
+contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearch.cpp
+contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp
+contactSearch/boundaries/periodicBoundaryContactSearch/periodicBoundaryContactSearch.cpp
+contactSearch/boundaries/boundaryContactSearchList.cpp
+
 contactSearch/contactSearch/contactSearch.cpp
-contactSearch/boundaries/searchBoundary.cpp
 contactSearch/ContactSearch/ContactSearchs.cpp
+
 interaction/interaction.cpp
-sphereInteraction/sphereInteractions.cpp
+sphereInteraction/sphereInteractionsLinearModels.cpp
+sphereInteraction/sphereInteractionsNonLinearModels.cpp
+sphereInteraction/sphereInteractionsNonLinearModModels.cpp
 )
 
 set(link_libs Kokkos::kokkos phasicFlow Property Particles Geometry)

src/Interaction/contactLists/unsortedContactList.hpp

@@ -17,10 +17,11 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
 #ifndef __unsortedContactList_hpp__
 #define __unsortedContactList_hpp__ 
 
+#include "unsortedPairs.hpp"
+
 namespace pFlow
 {
 

src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp

@@ -27,6 +27,7 @@ Licence:
 #include "box.hpp"
 #include "particles.hpp"
 #include "geometry.hpp"
+#include "boundaryContactSearchList.hpp"
 
 namespace pFlow
 {
@@ -50,6 +51,8 @@ private:
 	
 	uniquePtr<SearchMethodType> 	ppwContactSearch_ = nullptr;
 
+	boundaryContactSearchList 		csBoundaries_;
+
 public:
 
 	TypeInfoTemplate11("ContactSearch", SearchMethodType);
@@ -61,14 +64,24 @@ public:
 	 	const geometry& geom,
 	 	Timers& timers)
 	:
-		contactSearch(csDict, extDomain, prtcl, geom, timers)
+		contactSearch(
+			csDict, 
+			extDomain, 
+			prtcl, 
+			geom, 
+			timers),
+		csBoundaries_(
+			csDict, 
+			Particles().pStruct().boundaries(),
+			*this)
 	{
 
-		auto method = dict().getVal<word>("method");
+		/*auto method = dict().getVal<word>("method");
 				
-		auto nbDict = dict().subDict(method+"Info");
+		auto nbDict = dict().subDict(method+"Info");*/
 
-		real minD, maxD;
+		real minD;
+		real maxD;
 		this->Particles().boundingSphereMinMax(minD, maxD);
 		
 		const auto& position = this->Particles().pointPosition().deviceViewAll();
@@ -84,7 +97,7 @@ public:
 		ppwContactSearch_ = 
 			makeUnique<SearchMethodType>
 			(
-				nbDict,
+				dict(),
 				this->extendedDomainBox(),
 				minD,
 				maxD,
@@ -141,11 +154,27 @@ public:
 		}
 		ppTimer().end();
 
-				
-	
 		return true;
 	}
 
+	bool boundaryBroadSearch(
+		uint32 i,
+		uint32 iter,
+		real t,
+		real dt,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force = false)override
+	{
+		return csBoundaries_[i].broadSearch(
+			iter, 
+			t, 
+			dt, 
+			ppPairs, 
+			pwPairs, 
+			force);
+	}
+
 
 	bool enterBroadSearch(uint32 iter, real t, real dt)const override
 	{

src/Interaction/contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.cpp

@@ -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.
+
+-----------------------------------------------------------------------------*/
+
+#include "boundaryContactSearch.hpp"
+#include "contactSearch.hpp"
+
+pFlow::boundaryContactSearch::boundaryContactSearch(
+    const dictionary &dict,
+    const boundaryBase &boundary,
+    const contactSearch &cSearch)
+    : generalBoundary(
+          boundary,
+          cSearch.pStruct(),
+          "",
+          ""),
+      contactSearch_(cSearch),
+      updateInterval_(dict.getVal<uint32>("updateInterval"))
+{
+}
+
+pFlow::uniquePtr<pFlow::boundaryContactSearch>
+pFlow::boundaryContactSearch::create(
+	const dictionary &dict,
+	const boundaryBase &boundary,
+	const contactSearch &cSearch)
+{
+	
+	word bType = angleBracketsNames2( 
+		"boundaryContactSearch",
+		pFlowProcessors().localRunTypeName(), 
+		boundary.type());
+	word altBType{"boundaryContactSearch<none>"};
+
+	if( boundaryBasevCtorSelector_.search(bType) )
+	{
+		REPORT(2)<<"Creating contact search boundary "<< Green_Text(bType)<<
+		" for "<<boundary.name()<<endl;
+		return boundaryBasevCtorSelector_[bType](dict, boundary, cSearch);
+	}
+	else if(boundaryBasevCtorSelector_.search(altBType))
+	{
+		REPORT(2)<<"Creating contact search boundary "<< Green_Text(altBType)<<
+		" for "<<boundary.name()<<endl;
+		return boundaryBasevCtorSelector_[altBType](dict, boundary, cSearch);
+	}
+  else
+  {
+    printKeys
+		( 
+			fatalError << "Ctor Selector "<< bType<<
+      " and "<< altBType << " do not exist. \n"
+			<<"Avaiable ones are: \n"
+			,
+			boundaryBasevCtorSelector_
+		);
+		fatalExit;
+  }
+	
+	return nullptr;
+}

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

@@ -0,0 +1,114 @@
+/*------------------------------- 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 __boundaryContactSearch_hpp__
+#define __boundaryContactSearch_hpp__
+
+#include "generalBoundary.hpp"
+#include "contactSearchGlobals.hpp"
+#include "virtualConstructor.hpp"
+
+namespace pFlow
+{
+
+class contactSearch;
+
+class boundaryContactSearch
+    : public generalBoundary
+{
+private:
+    const contactSearch &contactSearch_;
+
+    /// @brief update interval in terms of iteration numebr
+    uint32 updateInterval_;
+
+    /// @brief  last iteration number which contact search has been performed
+    uint32 lastUpdated_ = 0;
+
+    /// @brief performed search?
+    bool performedSearch_ = false;
+
+public:
+    // type info
+    TypeInfo("boundaryContactSearch<none>");
+
+    boundaryContactSearch(
+        const dictionary &dict,
+        const boundaryBase &boundary,
+        const contactSearch &cSearch);
+
+    create_vCtor(
+        boundaryContactSearch,
+        boundaryBase,
+        (
+            const dictionary &dict,
+            const boundaryBase &boundary,
+            const contactSearch &cSearch),
+        (dict, boundary, cSearch));
+
+    add_vCtor(
+        boundaryContactSearch,
+        boundaryContactSearch,
+        boundaryBase);
+
+    const contactSearch &cSearch() const
+    {
+        return contactSearch_;
+    }
+
+    void fill(const std::any &val) override
+    {
+        return;
+    }
+
+    bool hearChanges(
+        real t,
+        real dt,
+        uint32 iter,
+        const message &msg,
+        const anyList &varList) override
+    {
+
+        if (msg.equivalentTo(message::BNDR_RESET))
+        {
+            // do nothing
+        }
+        return true;
+    }
+
+    virtual bool broadSearch(
+        uint32 iter,
+        real t,
+        real dt,
+        csPairContainerType &ppPairs,
+        csPairContainerType &pwPairs,
+        bool force = false)
+    {
+        return true;
+    }
+
+    static uniquePtr<boundaryContactSearch> create(
+        const dictionary &dict,
+        const boundaryBase &boundary,
+        const contactSearch &cSearch);
+};
+
+}
+
+#endif //__boundaryContactSearch_hpp__

src/Interaction/contactSearch/boundaries/boundaryContactSearchList.cpp

@@ -0,0 +1,29 @@
+#include "boundaryContactSearchList.hpp"
+#include "boundaryList.hpp"
+
+void pFlow::boundaryContactSearchList::setList(
+    const dictionary &dict, 
+    const contactSearch &cSearch)
+{
+    for(auto i=0; i<boundaries_.size(); i++)
+    {
+        this->set
+        (
+            i, 
+            boundaryContactSearch::create(dict, boundaries_[i], cSearch)
+        );
+    }
+}
+
+pFlow::boundaryContactSearchList::boundaryContactSearchList(
+    const dictionary &dict,
+    const boundaryList& bndrs,
+    const contactSearch &cSearch)
+    :
+    ListPtr(bndrs.size()), 
+    boundaries_(bndrs)
+{
+    setList(dict, cSearch);
+}
+
+

src/Interaction/contactSearch/boundaries/boundaryContactSearchList.hpp

@@ -0,0 +1,35 @@
+#include "ListPtr.hpp"
+#include "boundaryContactSearch.hpp"
+
+namespace pFlow
+{
+
+class boundaryList;
+class contactSearch;
+
+class boundaryContactSearchList
+:
+    public ListPtr<boundaryContactSearch>
+{
+private:
+
+    const boundaryList&         boundaries_;
+
+    void setList(
+        const dictionary& dict,
+        const contactSearch& cSearch);
+public:
+
+    TypeInfoNV("boundaryContactSearchList");
+
+    boundaryContactSearchList(
+        const dictionary& dict,
+        const boundaryList& bndrs,
+        const contactSearch& cSearch);
+    
+    ~boundaryContactSearchList()=default;
+    
+};
+
+
+}

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

@@ -0,0 +1,131 @@
+/*------------------------------- 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 "periodicBoundaryContactSearch.hpp"
+#include "contactSearch.hpp"
+#include "particles.hpp"
+#include "pointStructure.hpp"
+#include "geometry.hpp"
+
+
+void pFlow::periodicBoundaryContactSearch::setSearchBox()
+{
+    
+    auto db = pStruct().thisDomain().domainBox();
+    auto n1 = boundary().mirrorBoundary().boundaryPlane().normal();
+    auto l1 = boundary().mirrorBoundary().neighborLength();
+    auto n2 = boundary().boundaryPlane().normal();
+    auto l2 = boundary().neighborLength();
+
+    realx3 minP = db.minPoint() + (db.maxPoint()-db.minPoint())* n1+(n2*l2);
+    realx3 maxP = db.maxPoint() + (n1*l1);
+
+    searchBox_={minP, maxP};
+}
+
+pFlow::periodicBoundaryContactSearch::periodicBoundaryContactSearch(
+    const dictionary &dict,
+    const boundaryBase &boundary,
+    const contactSearch &cSearch)
+: 
+    boundaryContactSearch(dict, boundary, cSearch),
+    transferVec_(boundary.mirrorBoundary().displacementVectroToMirror()),
+    thisIndex_(thisBoundaryIndex()),
+    mirrorIndex_(mirrorBoundaryindex()),
+    diameter_(cSearch.Particles().boundingSphere())
+{
+    
+    if(thisIndex_%2==1)
+    {
+        masterSearch_ = true;
+        
+        setSearchBox();
+        
+        real minD;
+        real maxD;
+        cSearch.Particles().boundingSphereMinMax(minD, maxD);
+        
+        ppContactSearch_ = makeUnique<ppwBndryContactSearch>(
+            searchBox_,
+            maxD);
+
+        const auto& geom = cSearch.Geometry();
+
+        pwContactSearch_ = makeUnique<wallBoundaryContactSearch>(
+            0.5,
+            geom.numPoints(),
+            geom.size(),
+            geom.points().deviceViewAll(),
+            geom.vertices().deviceViewAll(),
+            geom.normals().deviceViewAll());
+    }
+    else
+    {
+        masterSearch_ = false;
+        searchBox_={{0,0,0},{0,0,0}};
+    }
+}
+
+bool pFlow::periodicBoundaryContactSearch::broadSearch
+(
+    uint32 iter, 
+    real t, 
+    real dt, 
+    csPairContainerType &ppPairs, 
+    csPairContainerType &pwPairs, 
+    bool force
+)
+{
+    if(masterSearch_)
+    {
+        
+        auto thisP = boundary().thisPoints();
+        auto thisDiams = diameter_.BoundaryField(thisIndex_).thisField();
+        auto mirrorP = mirrorBoundary().thisPoints();
+        auto mirrorDiams = diameter_.BoundaryField(mirrorIndex_).thisField();
+
+        ppContactSearch_().broadSearchPP(
+            ppPairs, 
+            thisP, 
+            thisDiams,
+            mirrorP,
+            mirrorDiams,
+            transferVec_);
+        
+        /*pwContactSearch_().broadSearch(
+            pwPairs, 
+            ppContactSearch_().searchCells(),
+            thisP,
+            thisDiams,
+            mirrorP,
+            mirrorDiams,
+            transferVec_,
+            ppContactSearch_().sizeRatio());*/
+
+            //output<<t<<"  boundary pp size "<< ppPairs.size()<<endl;
+            //output<<t<<"  boundary pw size "<< pwPairs.size()<<endl;
+        
+        return true;
+
+    }else
+    {
+        return true;
+    }
+}

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

@@ -0,0 +1,81 @@
+/*------------------------------- 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 __periodicBoundaryContactSearch_hpp__
+#define __periodicBoundaryContactSearch_hpp__
+
+#include "boundaryContactSearch.hpp"
+#include "box.hpp"
+#include "ppwBndryContactSearch.hpp"
+#include "pointFields.hpp"
+#include "wallBoundaryContactSearch.hpp"
+
+namespace pFlow
+{
+
+
+class periodicBoundaryContactSearch
+	: public boundaryContactSearch
+{
+private:
+	box searchBox_;
+
+	realx3 transferVec_;
+
+	uint32 thisIndex_;
+
+	uint32 mirrorIndex_;
+
+	uniquePtr<ppwBndryContactSearch> 	ppContactSearch_ = nullptr;
+
+	uniquePtr<wallBoundaryContactSearch> pwContactSearch_ = nullptr;
+
+	const realPointField_D 				&diameter_;
+
+	bool masterSearch_ = false;
+
+	void setSearchBox();
+
+public:
+	TypeInfo("boundaryContactSearch<regular,periodic>")
+
+		periodicBoundaryContactSearch(
+			const dictionary &dict,
+			const boundaryBase &boundary,
+			const contactSearch &cSearch);
+
+	~periodicBoundaryContactSearch() override = default;
+
+	add_vCtor(
+		boundaryContactSearch,
+		periodicBoundaryContactSearch,
+		boundaryBase);
+
+	bool broadSearch(
+		uint32 iter,
+		real t,
+		real dt,
+		csPairContainerType &ppPairs,
+		csPairContainerType &pwPairs,
+		bool force = false) override;
+};
+
+}
+
+#endif //__periodicBoundaryContactSearch_hpp__

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

@@ -0,0 +1,109 @@
+
+#include "ppwBndryContactSearch.hpp"
+#include "ppwBndryContactSearchKernels.hpp"
+#include "phasicFlowKokkos.hpp"
+#include "streams.hpp"
+
+void pFlow::ppwBndryContactSearch::checkAllocateNext(uint32 n)
+{
+	if( nextCapacity_ < n)
+	{
+		nextCapacity_ = n;	
+		reallocNoInit(next_, n);
+	}
+}
+
+void pFlow::ppwBndryContactSearch::nullifyHead()
+{
+	fill(head_, static_cast<uint32>(-1));
+}
+
+void pFlow::ppwBndryContactSearch::nullifyNext(uint32 n)
+{
+	fill(next_, 0u, n, static_cast<uint32>(-1));
+}
+
+void pFlow::ppwBndryContactSearch::buildList(
+	const deviceScatteredFieldAccess<realx3> &points)
+{
+	if(points.empty())return;
+	uint32 n = points.size();
+	checkAllocateNext(n);
+	nullifyNext(n);
+	nullifyHead();
+	
+	pFlow::pweBndryContactSearchKernels::buildNextHead(
+		points,
+		searchCells_,
+		head_,
+		next_
+	);
+}
+
+pFlow::ppwBndryContactSearch::ppwBndryContactSearch
+(
+    const box &domain,
+    real cellSize,
+	real sizeRatio
+)
+:
+	searchCells_(domain, cellSize),
+	head_("periodic:head",searchCells_.nx(), searchCells_.ny(), searchCells_.nz()),
+	sizeRatio_(sizeRatio)
+{
+
+}
+
+bool pFlow::ppwBndryContactSearch::broadSearchPP
+(
+	csPairContainerType &ppPairs, 
+	const deviceScatteredFieldAccess<realx3> &points, 
+	const deviceScatteredFieldAccess<real>& diams, 
+	const deviceScatteredFieldAccess<realx3> &mirrorPoints, 
+	const deviceScatteredFieldAccess<real>& mirrorDiams, 
+	const realx3& transferVec
+)
+{
+	
+	buildList(points);
+
+	uint32 nNotInserted = 1;
+
+	// loop until the container size fits the numebr of contact pairs
+	while (nNotInserted > 0)
+	{
+
+		nNotInserted = pFlow::pweBndryContactSearchKernels::broadSearchPP
+		(
+			ppPairs,
+			points,
+			diams,
+			mirrorPoints,
+			mirrorDiams,
+			transferVec,
+			head_,
+			next_,
+			searchCells_,
+			sizeRatio_
+		);
+		
+	
+		if(nNotInserted)
+		{
+			// - resize the container
+			//   note that getFull now shows the number of failed insertions.
+			uint32 len = max(nNotInserted,100u) ;
+			
+			auto oldCap = ppPairs.capacity();
+			
+			ppPairs.increaseCapacityBy(len);
+
+			INFORMATION<< "Particle-particle contact pair container capacity increased from "<<
+			oldCap << " to "<<ppPairs.capacity()<<" in peiodicBoundaryContactSearch."<<END_INFO;
+			
+		}
+
+	}	
+	
+    return true;
+}

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

@@ -0,0 +1,85 @@
+/*------------------------------- 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 __ppwBndryContactSearch_hpp__
+#define __ppwBndryContactSearch_hpp__
+
+#include "contactSearchGlobals.hpp"
+#include "scatteredFieldAccess.hpp"
+#include "cells.hpp"
+
+namespace pFlow
+{
+
+class ppwBndryContactSearch
+{
+public:
+	using HeadType = deviceViewType3D<uint32>;
+
+	using NextType = deviceViewType1D<uint32>;
+
+private:
+	cells searchCells_;
+
+	HeadType head_{"periodic::head", 1, 1, 1};
+
+	NextType next_{"periodic::next", 1};
+
+	real sizeRatio_ = 1.0;
+
+	uint32 nextCapacity_ = 0;
+
+	void checkAllocateNext(uint32 n);
+
+	void nullifyHead();
+
+	void nullifyNext(uint32 n);
+
+	void buildList(
+		const deviceScatteredFieldAccess<realx3> &points);
+
+public:
+	ppwBndryContactSearch(
+		const box &domain,
+		real cellSize,
+		real sizeRatio = 1.0);
+
+	bool broadSearchPP(
+		csPairContainerType &ppPairs,
+		const deviceScatteredFieldAccess<realx3> &points,
+		const deviceScatteredFieldAccess<real> &diams,
+		const deviceScatteredFieldAccess<realx3> &mirrorPoints,
+		const deviceScatteredFieldAccess<real> &mirrorDiams,
+		const realx3 &transferVec);
+
+	const auto& searchCells()const
+	{
+		return searchCells_;
+	}
+
+	real sizeRatio()const
+	{
+		return sizeRatio_;
+	}
+};
+
+}
+
+#endif //__ppwBndryContactSearch_hpp__

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

@@ -0,0 +1,105 @@
+#include "ppwBndryContactSearchKernels.hpp"
+
+INLINE_FUNCTION_HD
+bool sphereSphereCheckB(const pFlow::realx3& p1, const pFlow::realx3 p2, pFlow::real d1, pFlow::real d2)
+{
+	return pFlow::length(p2-p1) < 0.5*(d2+d1);
+}
+
+void pFlow::pweBndryContactSearchKernels::buildNextHead
+(
+    const deviceScatteredFieldAccess<realx3> &points,
+    const cells              &searchCells,
+    deviceViewType3D<uint32> &head, 
+    deviceViewType1D<uint32> &next 
+)
+{
+    if(points.empty())return;
+    
+    uint32 n= points.size();
+
+    Kokkos::parallel_for(
+		"pFlow::ppwBndryContactSearch::buildList",
+		deviceRPolicyStatic(0,n),
+		LAMBDA_HD(uint32 i)
+		{
+			int32x3 ind;
+			if( searchCells.pointIndexInDomain(points[i], ind) )
+			{
+				// discards points out of searchCell
+				uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
+				next[i] = old;
+			}
+		}
+	);
+	Kokkos::fence();
+}
+
+pFlow::uint32 pFlow::pweBndryContactSearchKernels::broadSearchPP
+(
+    csPairContainerType &ppPairs, 
+    const deviceScatteredFieldAccess<realx3> &points, 
+    const deviceScatteredFieldAccess<real> &diams, 
+    const deviceScatteredFieldAccess<realx3> &mirrorPoints, 
+    const deviceScatteredFieldAccess<real> &mirrorDiams, 
+    const realx3 &transferVec, 
+    const deviceViewType3D<uint32> &head, 
+    const deviceViewType1D<uint32> &next, 
+    const cells &searchCells,
+    const real sizeRatio
+)
+{
+    
+    if(points.empty()) return 0;
+	if(mirrorPoints.empty())return 0;
+
+	auto nMirror = mirrorPoints.size();
+    
+    uint32 getFull = 0;
+
+    Kokkos::parallel_reduce(
+		"pFlow::pweBndryContactSearchKernels::broadSearchPP",
+		deviceRPolicyStatic(0, nMirror),
+		LAMBDA_HD(const uint32 mrrI, uint32 &getFullUpdate)
+		{
+			realx3 p_m = mirrorPoints(mrrI) + transferVec;
+			
+			int32x3 ind_m;
+			if( !searchCells.pointIndexInDomain(p_m, ind_m))return;
+			
+			real   d_m = sizeRatio*mirrorDiams[mrrI];
+			
+			for(int ii=-1; ii<2; ii++)
+			{
+			for(int jj=-1; jj<2; jj++)
+			{
+			for(int kk =-1; kk<2; kk++)
+			{
+				auto ind = ind_m + int32x3{ii,jj,kk};
+				
+				if(!searchCells.inCellRange(ind))continue;
+
+				uint32 thisI = head(ind.x(),ind.y(),ind.z());
+				while (thisI!=-1)
+				{				
+					
+					auto d_n = sizeRatio*diams[thisI];
+					
+                    // first item is for this boundary and second itme, for mirror 
+					if(sphereSphereCheckB(p_m, points[thisI], d_m, d_n)&&
+                       ppPairs.insert(thisI,mrrI) == -1)
+					{
+						getFullUpdate++;	
+					}
+
+					thisI = next(thisI);
+				}
+			}
+			}
+			}
+		},
+        getFull
+	);
+
+    return getFull;
+}

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

@@ -0,0 +1,31 @@
+
+#include "contactSearchGlobals.hpp"
+#include "cells.hpp"
+#include "contactSearchFunctions.hpp"
+#include "scatteredFieldAccess.hpp"
+
+namespace pFlow::pweBndryContactSearchKernels
+{
+
+void buildNextHead(
+    const deviceScatteredFieldAccess<realx3> &points,
+    const cells              &searchCells,
+    deviceViewType3D<uint32> &head, 
+    deviceViewType1D<uint32> &next );
+
+
+uint32 broadSearchPP
+(
+	csPairContainerType 						&ppPairs,
+	const deviceScatteredFieldAccess<realx3> 	&points, 
+	const deviceScatteredFieldAccess<real> 		&diams, 
+	const deviceScatteredFieldAccess<realx3>	&mirrorPoints, 
+	const deviceScatteredFieldAccess<real> 		&mirrorDiams, 
+	const realx3 								&transferVec,
+    const deviceViewType3D<uint32>				&head,
+    const deviceViewType1D<uint32>				&next,
+    const cells									&searchCells,
+	real sizeRatio
+);
+
+}

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

@@ -0,0 +1,141 @@
+#include "wallBoundaryContactSearch.hpp"
+#include "streams.hpp"
+
+pFlow::wallBoundaryContactSearch::wallBoundaryContactSearch
+(
+    real cellExtent, 
+    uint32 numPoints, 
+    uint32 numElements, 
+    const ViewType1D<realx3, memory_space> &points, 
+    const ViewType1D<uint32x3, memory_space> &vertices, 
+    const ViewType1D<realx3, memory_space> &normals
+)
+:
+    cellExtent_( max(cellExtent, 0.5 ) ),
+    numElements_(numElements),
+    numPoints_(numPoints),
+    vertices_(vertices),
+    points_(points),
+	normals_(normals)
+{
+    allocateArrays();
+}
+
+bool pFlow::wallBoundaryContactSearch::build(const cells &searchBox, const realx3& transferVec)
+{
+    Kokkos::parallel_for(
+		"pFlow::cellsWallLevel0::build",
+		deviceRPolicyStatic(0,numElements_),
+		CLASS_LAMBDA_HD(uint32 i)
+		{
+			auto v = vertices_[i];
+			auto p1 = points_[v.x()]+transferVec;
+			auto p2 = points_[v.y()]+transferVec;
+			auto p3 = points_[v.z()]+transferVec;
+			
+			realx3 minP;
+			realx3 maxP;
+
+			searchBox.extendBox(p1, p2, p3, cellExtent_, minP, maxP);
+			elementBox_[i] = iBoxType(searchBox.pointIndex(minP), searchBox.pointIndex(maxP));
+            auto d = elementBox_[i].maxPoint()-elementBox_[i].minPoint();
+            validBox_[i] = (d.x()*d.y()*d.z())==0? 0:1; 
+		});
+	Kokkos::fence();
+
+	return true;
+}
+
+bool pFlow::wallBoundaryContactSearch::broadSearch
+(
+    csPairContainerType &pairs, 
+    const cells &searchCells, 
+    const deviceScatteredFieldAccess<realx3> &thisPoints, 
+    const deviceScatteredFieldAccess<real> &thisDiams, 
+    const deviceScatteredFieldAccess<realx3> &mirrorPoints, 
+    const deviceScatteredFieldAccess<real> &mirroDiams, 
+    const realx3 &transferVec, 
+    real sizeRatio
+)
+{
+    uint32 nNotInserted = 1;
+
+    while (nNotInserted>0u)
+    {
+        build(searchCells,{0,0,0});
+        nNotInserted = findPairsElementRangeCount(
+            pairs,
+            searchCells,
+            thisPoints,
+            thisDiams,
+            {0,0,0},
+            0
+        );
+
+        build(searchCells, transferVec);
+        nNotInserted += findPairsElementRangeCount(
+            pairs,
+            searchCells,
+            mirrorPoints,
+            mirroDiams,
+            transferVec,
+            BASE_MIRROR_WALL_INDEX
+        );
+
+        if(nNotInserted>0u)
+        {
+            //   note that getFull now shows the number of failed insertions.
+			uint32 incCap = max(nNotInserted,50u) ;
+			
+			auto oldCap = pairs.capacity();
+			
+			pairs.increaseCapacityBy(incCap);
+
+			INFORMATION<< "The contact pair container capacity increased from "<<
+			oldCap << " to "<<pairs.capacity()<<" in wallBoundaryContactSearch."<<END_INFO;
+        }
+    }
+    return true;
+}
+
+pFlow::uint32 pFlow::wallBoundaryContactSearch::findPairsElementRangeCount
+(
+    csPairContainerType &pairs, 
+    const cells &searchCells, 
+    const deviceScatteredFieldAccess<realx3> &pPoints, 
+    const deviceScatteredFieldAccess<real> &pDiams, 
+    const realx3 &transferVec, 
+    uint  baseTriIndex
+)
+{
+    
+    if(pPoints.empty())return 0u;
+
+    uint32 nNotInserted = 0;
+    uint32 nThis = pPoints.size();
+    
+    Kokkos::parallel_reduce(
+        "pFlow::wallBoundaryContactSearch::findPairsElementRangeCount",
+        deviceRPolicyDynamic(0,nThis),
+        LAMBDA_HD(uint32 i, uint32 &notInsertedUpdate)
+        {
+            auto p = pPoints[i]+transferVec;
+            int32x3 ind;
+            if( searchCells.pointIndexInDomain(p, ind) )
+            {
+                for(uint32 nTri=0; nTri<numElements_; nTri++)
+                {
+                    if( validBox_[nTri]== 0)continue;
+                    if( elementBox_[nTri].isInside(ind)&& 
+                        pairs.insert(i,nTri+baseTriIndex) == -1)
+                    {
+                        notInsertedUpdate++;    
+                    }
+                }
+            }
+        },
+        nNotInserted
+    );
+
+    return nNotInserted;
+}

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

@@ -0,0 +1,131 @@
+/*------------------------------- 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 __wallBoundaryContactSearch_hpp__
+#define __wallBoundaryContactSearch_hpp__
+
+
+#include "contactSearchGlobals.hpp"
+#include "contactSearchFunctions.hpp"
+#include "scatteredFieldAccess.hpp"
+#include "iBox.hpp"
+#include "cells.hpp"
+
+
+
+namespace pFlow
+{
+
+class wallBoundaryContactSearch
+{
+public:
+
+	using execution_space = csExecutionSpace;
+
+	using memory_space 	= typename execution_space::memory_space;
+
+	using iBoxType 		= iBox<int32>;	
+
+private:
+
+	// - box extent
+	real 			cellExtent_ = 0.5;
+
+	// - number of triangle elements 
+	uint32			numElements_ = 0;
+
+	// - number of points
+	uint32 			numPoints_ = 0;
+
+	// - ref to vectices (borrowed)
+	ViewType1D<uint32x3, memory_space> 	vertices_;
+
+	// - ref to points in the trisurface (borrowed)
+	ViewType1D<realx3, memory_space> 	points_;
+
+	// - ref to normal vectors of triangles (borrowed)
+	ViewType1D<realx3, memory_space> 	normals_;
+
+	// cell range of element/triangle bounding box
+	ViewType1D<iBoxType, memory_space> 	elementBox_;
+
+    ViewType1D<uint8, memory_space>     validBox_;
+	
+	
+	FUNCTION_H
+	void allocateArrays()
+	{
+		reallocNoInit( elementBox_, numElements_);
+        reallocNoInit( validBox_, numElements_);
+	}
+
+public:
+
+	TypeInfoNV("wallBoundaryContactSearch");
+
+	INLINE_FUNCTION_HD
+	wallBoundaryContactSearch()=default;
+
+	FUNCTION_H
+	wallBoundaryContactSearch(
+		real 					cellExtent,
+		uint32 					numPoints,
+		uint32 					numElements,
+		const ViewType1D<realx3,memory_space>& 		points, 
+		const ViewType1D<uint32x3,memory_space>& 	vertices,
+		const ViewType1D<realx3, memory_space>& normals);
+	
+
+
+	INLINE_FUNCTION_HD
+	uint32 numElements()const
+	{
+		return numElements_;
+	}
+
+	
+	bool build(const cells& searchBox, const realx3& transferVec);
+	
+	bool broadSearch(
+		csPairContainerType& pairs, 
+		const cells &searchCells, 
+		const deviceScatteredFieldAccess<realx3>& thisPoints,
+		const deviceScatteredFieldAccess<real>&  thisDiams,
+        const deviceScatteredFieldAccess<realx3>& mirrorPoints,
+		const deviceScatteredFieldAccess<real>&  mirroDiams,
+		const realx3& transferVec,
+		real sizeRatio);
+
+	uint32 findPairsElementRangeCount(
+		csPairContainerType &pairs, 
+        const cells &searchCells, 
+        const deviceScatteredFieldAccess<realx3> &pPoints, 
+        const deviceScatteredFieldAccess<real> &pDiams, 
+        const realx3 &transferVec, 
+        uint  baseTriIndex);
+	
+
+
+}; // wallBoundaryContactSearch
+
+} // pFlow
+
+
+#endif // __wallBoundaryContactSearch_hpp__

src/Interaction/contactSearch/boundaries/searchBoundary.cpp

@@ -1,36 +0,0 @@
-
-#include "searchBoundary.hpp"
-#include "contactSearch.hpp"
-
-
-pFlow::searchBoundary::searchBoundary
-(
-    const dictionary& dict,
-	const boundaryBase& boundary,
-	const contactSearch& cSearch
-)
-:
-    generalBoundary
-    (
-        boundary,
-        cSearch.pStruct(),
-        "",
-        ""
-    ),
-    contactSearch_(cSearch),
-    updateInterval_(dict.getVal<uint32>("updateInterval"))
-{
-
-}
-
-pFlow::uniquePtr<pFlow::searchBoundary> 
-    pFlow::searchBoundary::create
-    (
-        const dictionary &dict, 
-        const boundaryBase &boundary, 
-        const contactSearch &cSearch
-    )
-{
-    return nullptr;    
-}
-

src/Interaction/contactSearch/boundaries/searchBoundary.hpp

@@ -1,113 +0,0 @@
-/*------------------------------- 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 __searchBoundary_hpp__
-#define __searchBoundary_hpp__
-
-
-#include "generalBoundary.hpp"
-#include "virtualConstructor.hpp"
-
-namespace pFlow
-{
-
-class contactSearch;
-
-class searchBoundary
-:
-    public generalBoundary
-{
-private:
-
-    const contactSearch&  contactSearch_;
-
-    /// @brief update interval in terms of iteration numebr 
-    uint32      updateInterval_;
-
-    /// @brief  last iteration number which contact search has been performed
-	uint32 		lastUpdated_ 	= 0;
-
-	/// @brief performed search? 
-	bool 		performedSearch_ = false;
-
-public:
-
-    // type info
-    TypeInfo("searchBoundary<regular,none>");
-
-
-    searchBoundary(
-		const dictionary& dict,
-		const boundaryBase& boundary,
-		const contactSearch& cSearch);
-
-    create_vCtor
-    (
-        searchBoundary,
-        boundaryBase,
-        (
-            const dictionary& dict,
-			const boundaryBase& boundary,
-			const contactSearch& cSearch
-        ),
-        (dict, boundary, cSearch)
-    );
-
-    add_vCtor
-    (
-        searchBoundary,
-        searchBoundary,
-        boundaryBase
-    );
-
-    
-	void fill(const std::any& val)override
-    {
-        return;
-    }
-
-    bool hearChanges
-	(
-		real t,
-		real dt,
-		uint32 iter,
-		const message& msg, 
-    	const anyList& varList
-	) override
-    {
-		
-		if(msg.equivalentTo(message::BNDR_RESET))
-		{
-			//do nothing
-		}
-		return true;
-	}
-
-	static
-	uniquePtr<searchBoundary> create(
-		const dictionary& dict,
-		const boundaryBase& boundary,
-		const contactSearch& cSearch);
-
-};
-
-}
-
-#endif //__searchBoundary_hpp__
-

src/Interaction/contactSearch/contactSearch/contactSearch.cpp

@@ -68,7 +68,7 @@ pFlow::uniquePtr<pFlow::contactSearch> pFlow::contactSearch::create(
 	{
 		printKeys
 		( 
-			fatalError << "Ctor Selector "<< model << " dose not exist. \n"
+			fatalError << "Ctor Selector "<< model << " does not exist. \n"
 			<<"Avaiable ones are: \n\n"
 			,
 			dictionaryvCtorSelector_

src/Interaction/contactSearch/contactSearch/contactSearch.hpp

@@ -123,6 +123,16 @@ public:
 		csPairContainerType& ppPairs,
 		csPairContainerType& pwPairs,
 		bool force = false) = 0;
+	
+	virtual 
+	bool boundaryBroadSearch(
+		uint32 i,
+		uint32 iter,
+		real t,
+		real dt,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force = false)=0;
 
 	virtual 
 	bool enterBroadSearch(uint32 iter, real t, real dt)const = 0;

src/Interaction/contactSearch/contactSearchGlobals.hpp

@@ -33,6 +33,9 @@ using csIdType = uint32;
 
 using csPairContainerType   = unsortedPairs<DefaultExecutionSpace, uint32>;
 
+inline 
+const uint32 BASE_MIRROR_WALL_INDEX = 1000000;
+
 }
 
 #endif

src/Interaction/interaction/interaction.cpp

@@ -71,14 +71,15 @@ pFlow::uniquePtr<pFlow::interaction> pFlow::interaction::create
 			clType);
 
 
-	REPORT(1)<< "Selecting interaction model..."<<END_REPORT;
+	
+	REPORT(1)<<"Creating interaction "<<Green_Text(interactionModel)<<" . . ."<<END_REPORT;
 	if( systemControlvCtorSelector_.search(interactionModel) )
 	{
 		auto objPtr = 
 			systemControlvCtorSelector_[interactionModel]
 			(control, prtcl, geom);
 
-		REPORT(2)<<"Model "<<Green_Text(interactionModel)<<" is created."<<END_REPORT;
+		
 		return objPtr;
 	}
 	else

src/Interaction/interaction/interaction.hpp

@@ -79,6 +79,17 @@ public:
 			(control, prtcl, geom)
 		);
 	
+	inline
+	const auto& Particles()const
+	{
+		return particles_;
+	}
+
+	inline
+	const auto& Geometry()const
+	{
+		return geometry_;
+	}
 
 	static 
 	uniquePtr<interaction> create(

src/Interaction/sphereInteraction/boundaries/boundarySphereInteraction/boundarySphereInteraction.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.
+
+-----------------------------------------------------------------------------*/
+
+
+template <typename cFM, typename gMM>
+pFlow::boundarySphereInteraction<cFM, gMM>::boundarySphereInteraction(
+	const boundaryBase &boundary,
+	const sphereParticles &sphPrtcls,
+	const GeometryMotionModel &geomMotion)
+	: generalBoundary(boundary, sphPrtcls.pStruct(), "", ""),
+	  geometryMotion_(geomMotion),
+	  sphParticles_(sphPrtcls)
+{
+	ppPairs_ = makeUnique<ContactListType>(1);
+	pwPairs_ = makeUnique<ContactListType>(1);
+}
+
+template <typename cFM, typename gMM>
+pFlow::uniquePtr<pFlow::boundarySphereInteraction<cFM, gMM>>
+pFlow::boundarySphereInteraction<cFM, gMM>::create(
+	const boundaryBase &boundary,
+	const sphereParticles &sphPrtcls,
+	const GeometryMotionModel &geomMotion)
+{
+	word cfTypeName = ContactForceModel::TYPENAME();
+	word gmTypeName = MotionModel::TYPENAME();
+	word bType = boundary.type();
+
+	word boundaryTypeName = angleBracketsNames3(
+		"boundarySphereInteraction",
+		bType,
+		cfTypeName,
+		gmTypeName);
+
+	word altBTypeName = angleBracketsNames2(
+		"boundarySphereInteraction",
+		cfTypeName,
+		gmTypeName);
+
+	if (boundaryBasevCtorSelector_.search(boundaryTypeName))
+	{
+		REPORT(2) << "Creating boundry type " << Green_Text(boundaryTypeName) << 
+		" for boundary " << boundary.name() << " . . ." << END_REPORT;
+		return boundaryBasevCtorSelector_[boundaryTypeName](
+			boundary,
+			sphPrtcls,
+			geomMotion);
+	}
+	else if(boundaryBasevCtorSelector_[altBTypeName])
+	{
+		// if boundary condition is not implemented, the default is used
+		
+		REPORT(2) << "Creating boundry type " << Green_Text(altBTypeName) << 
+		" for boundary " << boundary.name() << " . . ." << END_REPORT;
+		return boundaryBasevCtorSelector_[altBTypeName](
+			boundary,
+			sphPrtcls,
+			geomMotion);
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< boundaryTypeName<<
+			" and "<< altBTypeName << " do not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			boundaryBasevCtorSelector_
+		);
+		fatalExit;
+	}
+
+	return nullptr;
+}

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

@@ -0,0 +1,168 @@
+/*------------------------------- 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 __boundarySphereInteraction_hpp__
+#define __boundarySphereInteraction_hpp__
+
+#include "virtualConstructor.hpp"
+#include "generalBoundary.hpp"
+#include "unsortedContactList.hpp"
+#include "sphereParticles.hpp"
+
+namespace pFlow
+{
+
+template<typename contactForceModel,typename geometryMotionModel>
+class boundarySphereInteraction
+:
+    public generalBoundary
+{
+public:
+	
+	using BoundarySphereInteractionType
+			= boundarySphereInteraction<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 = 
+        unsortedContactList<ModelStorage, DefaultExecutionSpace, IdType>;
+
+private:
+
+	const GeometryMotionModel& 			geometryMotion_;
+
+	/// const reference to sphere particles 
+	const sphereParticles& 				sphParticles_;
+    
+	uniquePtr<ContactListType> 			ppPairs_;
+
+	uniquePtr<ContactListType>			pwPairs_;
+
+public:
+
+	TypeInfoTemplate12("boundarySphereInteraction", ContactForceModel, MotionModel);
+
+	boundarySphereInteraction(
+		const boundaryBase& boundary,
+		const sphereParticles& sphPrtcls,
+		const GeometryMotionModel& geomMotion);
+
+	create_vCtor
+	(
+		BoundarySphereInteractionType,
+		boundaryBase,
+		(
+			const boundaryBase& boundary,
+			const sphereParticles& sphPrtcls,
+			const GeometryMotionModel& geomMotion
+		),
+		(boundary, sphPrtcls, geomMotion)
+	);
+
+	add_vCtor
+	(
+		BoundarySphereInteractionType,
+		BoundarySphereInteractionType,
+		boundaryBase
+	);
+
+	~boundarySphereInteraction()override=default;
+
+	const auto& sphParticles()const
+	{
+		return sphParticles_;
+	}
+
+	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_();
+	}
+
+	virtual
+	bool sphereSphereInteraction(
+		real dt,
+		const ContactForceModel& cfModel)
+	{
+		// for default boundary, no thing to be done 
+		return true;
+	}
+
+	bool hearChanges
+	(
+		real t,
+		real dt,
+		uint32 iter,
+		const message& msg, 
+    	const anyList& varList
+	) override
+    {
+		
+		notImplementedFunction;
+		return true;
+	}
+
+	void fill(const std::any& val)override
+	{
+		notImplementedFunction;
+	}
+
+	static
+	uniquePtr<BoundarySphereInteractionType> create(
+		const boundaryBase& boundary,
+		const sphereParticles& sphPrtcls,
+		const GeometryMotionModel& geomMotion
+	);
+
+};
+
+}
+
+#include "boundarySphereInteraction.cpp"
+
+#endif //__boundarySphereInteraction_hpp__

src/Interaction/sphereInteraction/boundaries/boundarySphereInteractionList.cpp

@@ -0,0 +1,23 @@
+
+
+template <typename CFModel, typename gMModel>
+pFlow::boundarySphereInteractionList<CFModel, gMModel>::boundarySphereInteractionList
+(
+    const sphereParticles &sphPrtcls, 
+    const gMModel &geomMotion
+)
+:
+	ListPtr<boundarySphereInteraction<CFModel,gMModel>>(6),
+	boundaries_(sphPrtcls.pStruct().boundaries())
+{
+
+	for(uint32 i=0; i<6; i++)
+	{
+		this->set(
+			i,
+			boundarySphereInteraction<CFModel, gMModel>::create(
+				boundaries_[i],
+				sphPrtcls, 
+				geomMotion));
+	}	
+}

src/Interaction/sphereInteraction/boundaries/boundarySphereInteractionList.hpp

@@ -0,0 +1,40 @@
+#ifndef __boundarySphereInteractionList_hpp__
+#define __boundarySphereInteractionList_hpp__
+
+
+#include "boundaryList.hpp"
+#include "ListPtr.hpp"
+#include "boundarySphereInteraction.hpp"
+
+
+namespace pFlow
+{
+
+
+template<typename contactForceModel,typename geometryMotionModel>
+class boundarySphereInteractionList
+:
+    public ListPtr<boundarySphereInteraction<contactForceModel,geometryMotionModel>>
+{
+private:
+
+    const boundaryList&         boundaries_;
+
+public:
+
+    boundarySphereInteractionList(
+        const sphereParticles& sphPrtcls,
+		const geometryMotionModel& geomMotion
+    );
+
+    ~boundarySphereInteractionList()=default;
+
+};
+
+
+
+}
+
+#include "boundarySphereInteractionList.cpp"
+
+#endif //__boundarySphereInteractionList_hpp__

src/Interaction/sphereInteraction/boundaries/createBoundarySphereInteraction.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 "boundarySphereInteraction.hpp"
+#include "periodicBoundarySphereInteraction.hpp"
+
+#define createBoundarySphereInteraction(ForceModel,GeomModel) 	\
+	template class pFlow::boundarySphereInteraction<			\
+	ForceModel,													\
+	GeomModel>;													\
+																\
+	template class pFlow::periodicBoundarySphereInteraction<	\
+	ForceModel,													\
+	GeomModel>;
+

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

@@ -0,0 +1,125 @@
+
+namespace pFlow::periodicBoundarySIKernels
+{
+
+template<typename ContactListType, typename ContactForceModel>
+void sphereSphereInteraction
+(
+    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<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::sphereSphereInteraction",
+        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];
+            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 Vi = thisVel[i];
+                auto Vj = mirrorVel[j];
+                auto wi = thisRVel[i];
+                auto wj = mirrorRVel[j];
+                auto Nij = (xj-xi)/dist;
+                auto Vr = Vi - Vj + cross((Ri*wi+Rj*wj), Nij);*/
+
+                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,
+                    ovrlp,
+                    Vr, Nij,
+                    history,
+                    FCn, FCt);
+
+                forceModel.rollingFriction(
+                    dt, i, j,
+                    propId_i, propId_j,
+                    Ri, Rj,
+                    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/sphereInteraction/boundaries/periodicBoundarySphereInteraction/periodicBoundarySphereInteraction.cpp

@@ -0,0 +1,65 @@
+/*------------------------------- 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::periodicBoundarySphereInteraction<cFM, gMM>::periodicBoundarySphereInteraction(
+	const boundaryBase &boundary,
+	const sphereParticles &sphPrtcls,
+	const GeometryMotionModel &geomMotion)
+	: boundarySphereInteraction<cFM,gMM>(boundary, sphPrtcls, geomMotion),
+	transferVec_(boundary.mirrorBoundary().displacementVectroToMirror())
+{
+	if(boundary.thisBoundaryIndex()%2==1)
+    {
+       masterInteraction_ = true;
+    }
+    else
+    {
+        masterInteraction_ = false;
+    }
+}
+
+template <typename cFM, typename gMM>
+bool pFlow::periodicBoundarySphereInteraction<cFM, gMM>::sphereSphereInteraction
+(
+	real dt,
+	const ContactForceModel &cfModel
+)
+{
+	if(!masterInteraction_) return true;
+	
+	pFlow::periodicBoundarySIKernels::sphereSphereInteraction(
+		dt,
+		this->ppPairs(),
+		cfModel,
+		transferVec_,
+		this->boundary().thisPoints(),
+		this->mirrorBoundary().thisPoints(),
+		this->sphParticles().diameter().deviceViewAll(),
+		this->sphParticles().propertyId().deviceViewAll(),
+		this->sphParticles().velocity().deviceViewAll(),
+		this->sphParticles().rVelocity().deviceViewAll(),
+		this->sphParticles().contactForce().deviceViewAll(),
+		this->sphParticles().contactTorque().deviceViewAll());
+
+	return true;
+}

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

@@ -0,0 +1,95 @@
+/*------------------------------- 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 __periodicBoundarySphereInteraction_hpp__
+#define __periodicBoundarySphereInteraction_hpp__
+
+#include "boundarySphereInteraction.hpp"
+
+namespace pFlow
+{
+
+template<typename contactForceModel,typename geometryMotionModel>
+class periodicBoundarySphereInteraction
+:
+    public boundarySphereInteraction<contactForceModel, geometryMotionModel>
+{
+public:
+
+    using PBSInteractionType = 
+        periodicBoundarySphereInteraction<contactForceModel,geometryMotionModel>;
+
+    using BSInteractionType = 
+        boundarySphereInteraction<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("boundarySphereInteraction", "periodic",ContactForceModel, MotionModel);
+
+	
+    periodicBoundarySphereInteraction(
+        const boundaryBase& boundary,
+		const sphereParticles& sphPrtcls,
+		const GeometryMotionModel& geomMotion
+    );
+
+    add_vCtor
+    (
+        BSInteractionType,
+        PBSInteractionType,
+        boundaryBase
+    );
+    
+    ~periodicBoundarySphereInteraction()override = default;
+
+
+    
+
+    bool sphereSphereInteraction(
+        real dt,
+		const ContactForceModel& cfModel)override;
+	
+};
+
+}
+
+#include "periodicBoundarySphereInteraction.cpp"
+
+
+#endif //__periodicBoundarySphereInteraction_hpp__

src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp

@@ -35,9 +35,18 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::createSphereInteraction()
 
 	uint32 nPrtcl = sphParticles_.size();
 
+	contactSearch_ = contactSearch::create(
+	subDict("contactSearch"),
+	sphParticles_.extendedDomain().domainBox(),
+	sphParticles_,
+	geometryMotion_,
+	timers());
+
 	ppContactList_ = makeUnique<ContactListType>(nPrtcl+1);
 	
-	pwContactList_ = makeUnique<ContactListType>(nPrtcl/4+1);
+	pwContactList_ = makeUnique<ContactListType>(nPrtcl/5+1);
+
+	
 
 	return true;
 }
@@ -127,19 +136,13 @@ pFlow::sphereInteraction<cFM,gMM, cLT>::sphereInteraction
 	interaction(control, prtcl, geom),
 	geometryMotion_(dynamic_cast<const GeometryMotionModel&>(geom)),
 	sphParticles_(dynamic_cast<const sphereParticles&>(prtcl)),
+	boundaryInteraction_(sphParticles_,geometryMotion_),
 	ppInteractionTimer_("sphere-sphere interaction", &this->timers()),
 	pwInteractionTimer_("sphere-wall interaction", &this->timers()),
-	contactListTimer_("contact list management", &this->timers()),
+	contactListMangementTimer_("contact-list management", &this->timers()),
-	contactListTimer0_("contact list clear", &this->timers())
+	boundaryInteractionTimer_("interaction for boundary", &this->timers())
 {
-	contactSearch_ = contactSearch::create(
+	
-	subDict("contactSearch"),
-	prtcl.extendedDomain().domainBox(),
-	prtcl,
-	geom,
-	timers());
-
-
 	if(!createSphereInteraction())
 	{
 		fatalExit;
@@ -163,14 +166,20 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 	//output<<"iter, t, dt "<< iter<<" "<< t << " "<<dt<<endl;
 	bool broadSearch = contactSearch_().enterBroadSearch(iter, t, dt);
 
-
+	
 	if(broadSearch)
 	{
-		contactListTimer0_.start();
+		contactListMangementTimer_.start();
 		ppContactList_().beforeBroadSearch();
 		pwContactList_().beforeBroadSearch();
-		contactListTimer0_.end();
+		contactListMangementTimer_.pause();
 	} 
+
+	for(uint32 i=0; i<6u; i++)
+	{
+		boundaryInteraction_[i].ppPairs().beforeBroadSearch();
+		boundaryInteraction_[i].pwPairs().beforeBroadSearch();
+	}
 	
 	if( sphParticles_.numActive()<=0)return true;
 	
@@ -187,14 +196,34 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 		fatalExit;
 	}
 	
-
+	for(uint32 i=0; i<6u; i++)
+	{
+		if( !contactSearch_().boundaryBroadSearch(
+			i,
+			iter,
+			t,
+			dt,
+			boundaryInteraction_[i].ppPairs(),
+			boundaryInteraction_[i].pwPairs()))
+		{
+			fatalErrorInFunction<<
+			"failed to perform broadSearch for boundary index "<<i<<endl;
+			return false;
+		}
+	}
 	
 	if(broadSearch && contactSearch_().performedBroadSearch())
 	{
-		contactListTimer_.start();
+		contactListMangementTimer_.resume();
 		ppContactList_().afterBroadSearch();
 		pwContactList_().afterBroadSearch();
-		contactListTimer_.end();
+		contactListMangementTimer_.end();
+	}
+
+	for(uint32 i=0; i<6u; i++)
+	{
+		boundaryInteraction_[i].ppPairs().afterBroadSearch();
+		boundaryInteraction_[i].pwPairs().afterBroadSearch();
 	}
 	
 	ppInteractionTimer_.start();
@@ -206,7 +235,14 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 		sphereWallInteraction();
 	pwInteractionTimer_.end();
 
-	
+	boundaryInteractionTimer_.start();
+	for(uint32 i=0; i<6u; i++)
+	{
+		boundaryInteraction_[i].sphereSphereInteraction(
+			dt, 
+			this->forceModel_());
+	}
+	boundaryInteractionTimer_.end();
 	return true;
 }
 

src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.hpp

@@ -23,7 +23,11 @@ Licence:
 
 #include "interaction.hpp"
 #include "sphereParticles.hpp"
+#include "boundarySphereInteractionList.hpp"
 #include "sphereInteractionKernels.hpp"
+//#include "unsortedContactList.hpp"
+
+
 
 namespace pFlow
 {
@@ -45,6 +49,8 @@ public:
 	using MotionModel 			= typename geometryMotionModel::MotionModel;
 
 	using ModelStorage 			= typename ContactForceModel::contactForceStorage;
+
+	using BoundaryListType 		= boundarySphereInteractionList<ContactForceModel,GeometryMotionModel>;
 	
 	using IdType 			= uint32;
 
@@ -53,6 +59,9 @@ public:
 	using ContactListType 	= 
 			contactListType<ModelStorage, DefaultExecutionSpace, IdType>;
 
+	//using BoundaryContactListType = unsortedContactList<ModelStorage, DefaultExecutionSpace, IdType>;
+	
+
 
 private:
 	
@@ -62,6 +71,9 @@ private:
 	/// const reference to particles 
 	const sphereParticles& 				sphParticles_;
 
+	/// particle-particle and particle-wall interactions at boundaries 
+	BoundaryListType 					boundaryInteraction_;	
+
 	/// contact search object for pp and pw interactions 
 	uniquePtr<contactSearch> 			contactSearch_ = nullptr;
 
@@ -72,7 +84,8 @@ private:
 	uniquePtr<ContactListType> 			ppContactList_ = nullptr;
 
 	/// contact list for particle-wall interactions (keeps the history)
-	uniquePtr<ContactListType> 			pwContactList_ = nullptr;	
+	uniquePtr<ContactListType> 			pwContactList_ = nullptr;
+
 
 	/// timer for particle-particle interaction computations
 	Timer 		ppInteractionTimer_;
@@ -80,9 +93,13 @@ private:
 	/// timer for particle-wall interaction computations
 	Timer       pwInteractionTimer_;
 
-	Timer 		contactListTimer_;
+	/// timer for managing contact lists (only inernal points)
+	Timer 		contactListMangementTimer_;
 
-	Timer 		contactListTimer0_;
+	/// timer for boundary interaction time 
+	Timer 		boundaryInteractionTimer_;
+
+	
 
 	bool createSphereInteraction();
 

src/Interaction/sphereInteraction/sphereInteractions.cpp

@@ -1,277 +0,0 @@
-/*------------------------------- 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 "sphereInteraction.hpp"
-
-#include "geometryMotions.hpp"
-
-#include "contactForceModels.hpp"
-#include "unsortedContactList.hpp"
-#include "sortedContactList.hpp"
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedLinearNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedLinearNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedLinearNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedLinearNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedLinearNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedLinearNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedLinearNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedLinearNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::sortedContactList>;
-
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedLinearNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedLinearNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedLinearNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedLinearNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::sortedContactList>;
-
-
-/*template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedLinearNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedLinearNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedLinearNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedLinearNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::sortedContactList>;*/
-
-/// non-linear models 
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::sortedContactList>;
-
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::sortedContactList>;
-
-/*template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::sortedContactList>;*/
-
-
-// - nonLinearMod models
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearModNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearModNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
-	pFlow::stationaryGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearModNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearModNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
-	pFlow::rotationAxisMotionGeometry,
-	pFlow::sortedContactList>;
-
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearModNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearModNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
-	pFlow::vibratingMotionGeometry,
-	pFlow::sortedContactList>;
-
-
-/*template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearModNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::limitedNonLinearModNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::sortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::unsortedContactList>;
-
-template class pFlow::sphereInteraction<
-	pFlow::cfModels::nonLimitedNonLinearModNormalRolling,
-	pFlow::multiRotationAxisMotionGeometry,
-	pFlow::sortedContactList>;
-*/

src/Interaction/sphereInteraction/sphereInteractionsLinearModels.cpp

@@ -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.
+
+-----------------------------------------------------------------------------*/
+
+#include "sphereInteraction.hpp"
+#include "geometryMotions.hpp"
+#include "contactForceModels.hpp"
+#include "unsortedContactList.hpp"
+#include "sortedContactList.hpp"
+#include "createBoundarySphereInteraction.hpp"
+
+
+
+#define createInteraction(ForceModel,GeomModel) 	\
+													\
+	template class pFlow::sphereInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::unsortedContactList>;				\
+													\
+	template class pFlow::sphereInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::sortedContactList>;					\
+	createBoundarySphereInteraction(ForceModel, GeomModel)
+
+
+// stationaryGeometry
+createInteraction(pFlow::cfModels::limitedLinearNormalRolling, pFlow::stationaryGeometry);
+createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::stationaryGeometry);
+
+// rotationAxisMotionGeometry
+createInteraction(pFlow::cfModels::limitedLinearNormalRolling, pFlow::rotationAxisMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::rotationAxisMotionGeometry);
+
+// vibratingMotionGeometry
+createInteraction(pFlow::cfModels::limitedLinearNormalRolling, pFlow::vibratingMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::vibratingMotionGeometry);
+
+// multiRotationAxisMotionGeometry
+//createInteraction(pFlow::cfModels::limitedLinearNormalRolling, pFlow::multiRotationAxisMotionGeometry);
+//createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::multiRotationAxisMotionGeometry);
+

src/Interaction/sphereInteraction/sphereInteractionsNonLinearModModels.cpp

@@ -0,0 +1,58 @@
+/*------------------------------- 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 "sphereInteraction.hpp"
+#include "geometryMotions.hpp"
+#include "contactForceModels.hpp"
+#include "unsortedContactList.hpp"
+#include "sortedContactList.hpp"
+#include "createBoundarySphereInteraction.hpp"
+
+
+
+#define createInteraction(ForceModel,GeomModel) 	\
+													\
+	template class pFlow::sphereInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::unsortedContactList>;				\
+													\
+	template class pFlow::sphereInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::sortedContactList>;					\
+	createBoundarySphereInteraction(ForceModel, GeomModel)
+
+
+// stationaryGeometry
+createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::stationaryGeometry);
+createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::stationaryGeometry);
+
+// rotationAxisMotionGeometry
+createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::rotationAxisMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::rotationAxisMotionGeometry);
+
+// vibratingMotionGeometry
+createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::vibratingMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::vibratingMotionGeometry);
+
+// multiRotationAxisMotionGeometry
+//createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::multiRotationAxisMotionGeometry);
+//createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::multiRotationAxisMotionGeometry);

src/Interaction/sphereInteraction/sphereInteractionsNonLinearModels.cpp

@@ -0,0 +1,58 @@
+/*------------------------------- 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 "sphereInteraction.hpp"
+#include "geometryMotions.hpp"
+#include "contactForceModels.hpp"
+#include "unsortedContactList.hpp"
+#include "sortedContactList.hpp"
+#include "createBoundarySphereInteraction.hpp"
+
+
+
+#define createInteraction(ForceModel,GeomModel) 	\
+													\
+	template class pFlow::sphereInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::unsortedContactList>;				\
+													\
+	template class pFlow::sphereInteraction< 		\
+		ForceModel,									\
+		GeomModel,									\
+		pFlow::sortedContactList>;					\
+	createBoundarySphereInteraction(ForceModel, GeomModel)
+
+
+// stationaryGeometry
+createInteraction(pFlow::cfModels::limitedNonLinearNormalRolling, pFlow::stationaryGeometry);
+createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::stationaryGeometry);
+
+// rotationAxisMotionGeometry
+createInteraction(pFlow::cfModels::limitedNonLinearNormalRolling, pFlow::rotationAxisMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::rotationAxisMotionGeometry);
+
+// vibratingMotionGeometry
+createInteraction(pFlow::cfModels::limitedNonLinearNormalRolling, pFlow::vibratingMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::vibratingMotionGeometry);
+
+// multiRotationAxisMotionGeometry
+//createInteraction(pFlow::cfModels::limitedNonLinearNormalRolling, pFlow::multiRotationAxisMotionGeometry);
+//createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::multiRotationAxisMotionGeometry);