merge from main

.gitignore

@@ -1,6 +1,12 @@
+# files 
+.clang-format
+.vscode
+.dependencygraph
+
 # Prerequisites
 *.d
 
+
 # Compiled Object files
 *.slo
 *.lo

CMakeLists.txt

@@ -1,75 +1,73 @@
 cmake_minimum_required(VERSION 3.16 FATAL_ERROR)
 
 # set the project name and version
-project(phasicFlow VERSION 0.1 )
+project(phasicFlow VERSION 1.0 )
 
 set(CMAKE_CXX_STANDARD 17 CACHE STRING "" FORCE)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 set(CMAKE_INSTALL_PREFIX ${phasicFlow_SOURCE_DIR} CACHE PATH "Install path of phasicFlow" FORCE)
-set(CMAKE_BUILD_TYPE Release CACHE STRING "build type" FORCE)
+set(CMAKE_BUILD_TYPE Debug CACHE STRING "build type" FORCE)
+set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
+mark_as_advanced(FORCE var BUILD_SHARED_LIBS)
 
 message(STATUS ${CMAKE_INSTALL_PREFIX})
 
-mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
-mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMP)
-mark_as_advanced(FORCE var Kokkos_ENABLE_SERIAL)
-mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
-mark_as_advanced(FORCE var BUILD_SHARED_LIBS)
+include(cmake/globals.cmake)
 
-option(USE_STD_PARALLEL_ALG "Use TTB std parallel algorithms" ON)
+#Kokkos directory to be included 
+set(Kokkos_Source_DIR)
+
+if(DEFINED ENV{Kokkos_DIR})
+   set(Kokkos_Source_DIR $ENV{Kokkos_DIR})
+else()
+  set(Kokkos_Source_DIR $ENV{HOME}/Kokkos/kokkos)
+endif()
+message(STATUS "Kokkos source directory is ${Kokkos_Source_DIR}")
+add_subdirectory(${Kokkos_Source_DIR} ./kokkos)
+Kokkos_cmake_settings()
+
+
+option(pFlow_STD_Parallel_Alg "Use TTB std parallel algorithms" ON)
 option(pFlow_Build_Serial "Build phasicFlow and backends for serial execution" OFF)
 option(pFlow_Build_OpenMP "Build phasicFlow and backends for OpenMP execution" OFF)
 option(pFlow_Build_Cuda   "Build phasicFlow and backends for Cuda execution" OFF)
-option(pFlow_Build_Double "Build phasicFlow with double precision variables" ON)
+option(pFlow_Build_Double "Build phasicFlow with double precision floating-oint variables" ON)
+option(pFlow_Build_MPI    "Build for MPI parallelization. This will enable multi-gpu run, CPU run on clusters (distributed memory machine). Use this combination Cuda+MPI, OpenMP + MPI or Serial+MPI " OFF)
 
-set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
 
 if(pFlow_Build_Serial)
 	set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE) 
 	set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "OpenMP execution" FORCE)
 	set(Kokkos_ENABLE_CUDA OFF CACHE BOOL "Cuda execution" FORCE)
 	set(Kokkos_ENABLE_CUDA_LAMBDA OFF CACHE BOOL "Cuda execution" FORCE)
+	set(Kokkos_ENABLE_CUDA_CONSTEXPR OFF CACHE BOOL "Enable constexpr on cuda code" FORCE)
 elseif(pFlow_Build_OpenMP )
 	set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE) 
 	set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "OpenMP execution" FORCE)
 	set(Kokkos_ENABLE_CUDA OFF CACHE BOOL "Cuda execution" FORCE)
 	set(Kokkos_ENABLE_CUDA_LAMBDA OFF CACHE BOOL "Cuda execution" FORCE)
 	set(Kokkos_DEFAULT_HOST_PARALLEL_EXECUTION_SPACE SERIAL CACHE STRING "" FORCE)
+	set(Kokkos_ENABLE_CUDA_CONSTEXPR OFF CACHE BOOL "Enable constexpr on cuda code" FORCE)
 elseif(pFlow_Build_Cuda)
 	set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "Serial execution" FORCE) 
 	set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "OpenMP execution" FORCE)
 	set(Kokkos_ENABLE_CUDA ON CACHE BOOL "Cuda execution" FORCE)
 	set(Kokkos_ENABLE_CUDA_LAMBDA ON CACHE BOOL "Cuda execution" FORCE)
+	set(Kokkos_ENABLE_CUDA_CONSTEXPR ON CACHE BOOL "Enable constexpr on cuda code" FORCE)
 endif()
 
-
-
-include(cmake/globals.cmake)
-message(STATUS "Valid file extensions are ${validFiles}")
+if(pFlow_Build_MPI)
+	find_package(MPI REQUIRED)
+endif()
 
 include(cmake/makeLibraryGlobals.cmake)
 include(cmake/makeExecutableGlobals.cmake)
 
 configure_file(phasicFlowConfig.H.in phasicFlowConfig.H)
-
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 #add a global include directory 
 include_directories(src/setHelpers src/demComponent "${PROJECT_BINARY_DIR}")
 
-#main subdirectories of the code 
-set(Kokkos_Source_DIR)
-
-if(DEFINED ENV{Kokkos_DIR})
-   set(Kokkos_Source_DIR $ENV{Kokkos_DIR})
-#  add_subdirectory($ENV{Kokkos_DIR} ${phasicFlow_BINARY_DIR}/kokkos)
-#  message(STATUS "Kokkos directory is $ENV{Kokkos_DIR}")
-else()
-#  add_subdirectory($ENV{HOME}/Kokkos/kokkos ${phasicFlow_BINARY_DIR}/kokkos)
-  set(Kokkos_Source_DIR $ENV{HOME}/Kokkos/kokkos)
-endif()
-
-message(STATUS "Kokkos source directory is ${Kokkos_Source_DIR}")
-add_subdirectory(${Kokkos_Source_DIR} ${phasicFlow_BINARY_DIR}/kokkos)
-
 add_subdirectory(src)
 
 add_subdirectory(solvers)
@@ -77,12 +75,9 @@ add_subdirectory(solvers)
 add_subdirectory(utilities)
 
 add_subdirectory(DEMSystems)
-#add_subdirectory(testIO)
-
 
 install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H"
-  DESTINATION include
-)
+  DESTINATION include)
 
 include(InstallRequiredSystemLibraries)
 set(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENSE")

DEMSystems/CMakeLists.txt

@@ -1,9 +1,11 @@
 
 set(SourceFiles 
-DEMSystem/DEMSystem.cpp
-sphereDEMSystem/sphereDEMSystem.cpp
-sphereDEMSystem/sphereFluidParticles.cpp
 domainDistribute/domainDistribute.cpp
+DEMSystem/DEMSystem.cpp
+sphereDEMSystem/sphereFluidParticles.cpp
+sphereDEMSystem/sphereDEMSystem.cpp
+grainDEMSystem/grainFluidParticles.cpp
+grainDEMSystem/grainDEMSystem.cpp
 )
 
 set(link_libs Kokkos::kokkos phasicFlow Particles Geometry Property Interaction Interaction Utilities)

DEMSystems/DEMSystem/DEMSystem.cpp

@@ -33,14 +33,14 @@ pFlow::DEMSystem::DEMSystem(
 	ControlDict_()
 {
   
-  REPORT(0)<<"Initializing host/device execution spaces . . . \n";
-	REPORT(1)<<"Host execution space is "<< greenText(DefaultHostExecutionSpace::name())<<endREPORT;
-	REPORT(1)<<"Device execution space is "<<greenText(DefaultExecutionSpace::name())<<endREPORT;
+	REPORT(0)<<"Initializing host/device execution spaces . . . \n";
+	REPORT(1)<<"Host execution space is "<< Green_Text(DefaultHostExecutionSpace::name())<<END_REPORT;
+	REPORT(1)<<"Device execution space is "<<Green_Text(DefaultExecutionSpace::name())<<END_REPORT;
 
  	// initialize Kokkos
 	Kokkos::initialize( argc, argv ); 
 
-	REPORT(0)<<"\nCreating Control repository . . ."<<endREPORT;
+	REPORT(0)<<"\nCreating Control repository . . ."<<END_REPORT;
 	Control_ = makeUnique<systemControl>(
 		ControlDict_.startTime(),
 		ControlDict_.endTime(),
@@ -87,4 +87,3 @@ pFlow::uniquePtr<pFlow::DEMSystem>
 
 	return nullptr;
 }
-

DEMSystems/DEMSystem/DEMSystem.hpp

@@ -25,6 +25,7 @@ Licence:
 
 #include "types.hpp"
 #include "span.hpp"
+#include "box.hpp"
 #include "virtualConstructor.hpp"
 #include "uniquePtr.hpp"
 #include "systemControl.hpp"
@@ -60,6 +61,7 @@ public:
 
 	DEMSystem(const DEMSystem&)=delete;
 
+	/// @brief no assignment 
 	DEMSystem& operator = (const DEMSystem&)=delete;
 
 	create_vCtor(
@@ -111,19 +113,34 @@ public:
 	int32 numParInDomain(int32 di)const = 0;
 	
 	virtual
-	std::vector<int32> numParInDomain()const = 0;
+	std::vector<int32> numParInDomains()const = 0;
 
 	virtual 
-	span<const int32> parIndexInDomain(int32 di)const = 0;
+	span<const int32> parIndexInDomain(int32 domIndx)const = 0;
 
 	virtual
-	span<real> parDiameter() = 0; 
+	span<real> diameter() = 0; 
 
 	virtual 
-	span<realx3> parVelocity() = 0;
+	span<real> courseGrainFactor() = 0;
 
 	virtual
-	span<realx3> parPosition() = 0;
+	span<realx3> acceleration()=0;
+
+	virtual
+	span<realx3> velocity() = 0;
+
+	virtual
+	span<realx3> position() = 0;
+
+	virtual
+	span<realx3> rAcceleration()=0;
+
+	virtual
+	span<realx3> rVelocity() = 0;
+
+	virtual
+	span<realx3> rPosition() = 0;
 
 	virtual
 	span<realx3> parFluidForce() = 0;
@@ -153,7 +170,6 @@ public:
 	bool iterate(real upToTime) = 0;
 	
 	
-
 	static 
 	uniquePtr<DEMSystem>
 		create(
@@ -162,8 +178,6 @@ public:
 			int argc, 
 			char* argv[]);
 
-
-
 };
 
 

DEMSystems/domainDistribute/domainDistribute.cpp

@@ -33,7 +33,7 @@ void pFlow::domainDistribute::clcDomains(const std::vector<box>& domains)
 
 
 pFlow::domainDistribute::domainDistribute(
-	const Vector<box>& domains,
+	const std::vector<box>& domains,
 	real maxBoundingBox)
 :
 numDomains_(domains.size()),
@@ -47,10 +47,10 @@ maxBoundingBoxSize_(maxBoundingBox)
 }
 
 bool pFlow::domainDistribute::locateParticles(
-		ViewType1D<realx3,HostSpace> points, includeMask mask)
+		ViewType1D<realx3,HostSpace> points, const pFlagTypeHost& mask)
 {
 	
-	range activeRange = mask.activeRange();
+	const rangeU32 activeRange = mask.activeRange();
 
 
 	for(int32 di=0; di<numDomains_; di++)
@@ -59,7 +59,7 @@ bool pFlow::domainDistribute::locateParticles(
 	}
 
 	
-	for(int32 i=activeRange.first; i<activeRange.second; i++)
+	for(int32 i=activeRange.start(); i<activeRange.end(); i++)
 	{
 		if(mask(i))
 		{

DEMSystems/domainDistribute/domainDistribute.hpp

@@ -43,19 +43,16 @@ protected:
 	
 	int32Vector_H 							numParInDomain_;
 
-
 	real maxBoundingBoxSize_;
 
 	real domainExtension_ = 1.0;
 
-	using includeMask = typename pointStructure::activePointsHost;
-
 	void clcDomains(const std::vector<box>& domains);
 
 public:
 
 	domainDistribute(
-		const Vector<box>& domains,
+		const std::vector<box>& domains,
 		real maxBoundingBox);
 
 	~domainDistribute()=default;
@@ -78,7 +75,7 @@ public:
 	{
 		return
 		span<const int32>(
-			particlesInDomains_[di].hostVectorAll().data(),
+			particlesInDomains_[di].hostViewAll().data(),
 			numParInDomain_[di]
 			);
 	}
@@ -91,7 +88,7 @@ public:
 
 	//template<typename includeMask>
 	bool locateParticles(
-		ViewType1D<realx3,HostSpace> points, includeMask mask);
+		ViewType1D<realx3,HostSpace> points, const pFlagTypeHost& mask);
 	
 };
 

DEMSystems/grainDEMSystem/grainDEMSystem.cpp

@@ -0,0 +1,268 @@
+/*------------------------------- 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 "grainDEMSystem.hpp"
+#include "vocabs.hpp"
+
+bool pFlow::grainDEMSystem::loop()
+{
+
+	do 
+	{
+		//
+		if(! insertion_().insertParticles( 
+			Control().time().currentIter(),
+			Control().time().currentTime(),
+			Control().time().dt()	) )
+		{
+			fatalError<<
+			"particle insertion failed in grainDEMSystem.\n";
+			return false;
+		}	
+
+		geometry_->beforeIteration();
+		
+		interaction_->beforeIteration();
+
+		particles_->beforeIteration();
+		
+		
+		interaction_->iterate();
+
+		particles_->iterate();
+
+		geometry_->iterate();
+
+		particles_->afterIteration();
+		
+		geometry_->afterIteration();
+
+	}while(Control()++);
+
+	return true;
+}
+
+pFlow::grainDEMSystem::grainDEMSystem(
+		word  demSystemName,
+		const std::vector<box>& domains,
+		int argc, 
+		char* argv[],
+		bool requireRVel)
+:
+	DEMSystem(demSystemName, domains, argc, argv),
+	requireRVel_(requireRVel)
+{ 
+	
+	REPORT(0)<<"\nReading proprties . . . "<<END_REPORT;
+	property_ = makeUnique<property>(
+		propertyFile__,
+    	Control().caseSetup().path());
+
+	REPORT(0)<< "\nCreating surface geometry for grainDEMSystem . . . "<<END_REPORT;
+	geometry_ = geometry::create(Control(), Property());
+
+
+	REPORT(0)<<"\nReading grain particles . . ."<<END_REPORT;
+	particles_ = makeUnique<grainFluidParticles>(Control(), Property());
+
+
+	insertion_ = makeUnique<grainInsertion>( 
+		particles_(), 
+		particles_().grains());
+	
+	REPORT(0)<<"\nCreating interaction model for grain-grain contact and grain-wall contact . . ."<<END_REPORT;
+	interaction_ = interaction::create(
+		Control(),
+		Particles(),
+		Geometry());
+
+	auto maxD = maxBounndingSphereSize();
+
+	particleDistribution_ = makeUnique<domainDistribute>(domains, maxD);
+
+}
+
+
+pFlow::grainDEMSystem::~grainDEMSystem()
+{
+	
+}
+
+bool pFlow::grainDEMSystem::updateParticleDistribution(
+	real extentFraction,
+	const std::vector<box> domains)
+{
+
+	if(!particleDistribution_->changeDomainsSize(
+		extentFraction,
+		maxBounndingSphereSize(),
+		domains))
+	{
+		fatalErrorInFunction<<
+		"cannot change the domain size"<<endl;
+		return false;
+	}
+
+	if(!particleDistribution_->locateParticles(
+			positionHost_, 
+			particles_->pStruct().activePointsMaskHost()))
+	{
+		fatalErrorInFunction<<
+		"error in locating particles among sub-domains"<<endl;
+		return false;
+	}
+
+	return true;
+}
+
+pFlow::int32 
+	pFlow::grainDEMSystem::numParInDomain(int32 di)const
+{
+	return particleDistribution_().numParInDomain(di);
+}
+
+std::vector<pFlow::int32> 
+	pFlow::grainDEMSystem::numParInDomains()const
+{
+	const auto& distribute = particleDistribution_();
+	int32 numDomains = distribute.numDomains();
+	std::vector<int32> nums(numDomains);
+	for(int32 i=0; i<numDomains; i++)
+	{
+		nums[i] = distribute.numParInDomain(i);
+	}
+
+	return nums;
+}
+
+pFlow::span<const pFlow::int32> 
+pFlow::grainDEMSystem::parIndexInDomain(int32 di)const
+{
+	return particleDistribution_->particlesInDomain(di);
+}
+
+pFlow::span<pFlow::real> pFlow::grainDEMSystem::diameter() 
+{	
+	return span<real>(diameterHost_.data(), diameterHost_.size());
+}
+
+pFlow::span<pFlow::real> pFlow::grainDEMSystem::courseGrainFactor()
+{
+    return span<real>(particles_->courseGrainFactorHost().data(), particles_->courseGrainFactorHost().size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::acceleration()
+{
+    return span<realx3>(nullptr, 0);
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::velocity()  
+{
+	return span<realx3>(velocityHost_.data(), velocityHost_.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::position()  
+{
+	return span<realx3>(positionHost_.data(), positionHost_.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::rAcceleration()
+{
+    return span<realx3>(nullptr, 0);
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::rVelocity()
+{
+    return span<realx3>(rVelocityHost_.data(), rVelocityHost_.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::rPosition()
+{
+    return span<realx3>(nullptr, 0);
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::parFluidForce() 
+{
+	auto& hVec =  particles_->fluidForceHost();
+	
+	return span<realx3>(hVec.data(), hVec.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::parFluidTorque() 
+{
+	auto& hVec =  particles_->fluidTorqueHost();
+	return span<realx3>(hVec.data(), hVec.size());
+}
+
+bool pFlow::grainDEMSystem::sendFluidForceToDEM() 
+{
+	particles_->fluidForceHostUpdatedSync();
+	return true;
+}
+
+bool pFlow::grainDEMSystem::sendFluidTorqueToDEM()
+{
+	particles_->fluidTorqueHostUpdatedSync();
+	return true;
+}
+
+bool pFlow::grainDEMSystem::beforeIteration()
+{
+	velocityHost_ = std::as_const(particles_()).velocity().hostView();
+	positionHost_ = std::as_const(particles_()).pointPosition().hostView();
+	diameterHost_ = particles_->diameter().hostView();
+
+	if(requireRVel_)
+		rVelocityHost_ = std::as_const(particles_()).rVelocity().hostView();
+	
+
+	return true;
+}
+
+
+bool pFlow::grainDEMSystem::iterate(
+	real upToTime, 
+	real timeToWrite, 
+	word timeName) 
+{
+
+	Control().time().setStopAt(upToTime);
+	Control().time().setOutputToFile(timeToWrite, timeName);
+	
+	return loop();
+
+	return true;
+}
+
+bool pFlow::grainDEMSystem::iterate(real upToTime)
+{
+	Control().time().setStopAt(upToTime);
+	return loop();
+	return true;
+}
+
+pFlow::real 
+	pFlow::grainDEMSystem::maxBounndingSphereSize()const
+{
+	real minD, maxD;
+	particles_->boundingSphereMinMax(minD, maxD);
+
+	return maxD;
+}

DEMSystems/grainDEMSystem/grainDEMSystem.hpp

@@ -0,0 +1,163 @@
+/*------------------------------- 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 __grainDEMSystem_hpp__
+#define __grainDEMSystem_hpp__
+
+#include "DEMSystem.hpp"
+#include "property.hpp"
+#include "uniquePtr.hpp"
+#include "geometry.hpp"
+#include "grainFluidParticles.hpp"
+#include "interaction.hpp"
+#include "Insertions.hpp"
+#include "domainDistribute.hpp"
+
+
+namespace pFlow
+{
+
+class grainDEMSystem
+:
+	public DEMSystem
+{
+protected:
+	
+	// protected members 
+
+	uniquePtr<property> 			property_ = nullptr;
+
+	uniquePtr<geometry> 			geometry_ = nullptr;
+
+	uniquePtr<grainFluidParticles>  particles_ = nullptr;
+
+	uniquePtr<grainInsertion> 		insertion_ = nullptr;
+
+	uniquePtr<interaction> 			interaction_ = nullptr;
+
+	uniquePtr<domainDistribute> 	particleDistribution_=nullptr;
+
+	// to be used for CPU communications 
+	ViewType1D<realx3, HostSpace> 	velocityHost_;
+
+	ViewType1D<realx3, HostSpace> 	positionHost_;
+
+	ViewType1D<real, HostSpace> 	diameterHost_;
+
+	bool 							requireRVel_ = false;
+
+	ViewType1D<realx3, HostSpace> 	rVelocityHost_;
+
+
+//  protected member functions
+	auto& Property()
+	{
+		return property_();
+	}
+
+	auto& Geometry()
+	{
+		return geometry_();
+	}
+
+	auto& Particles()
+	{
+		return particles_();
+	}
+
+	auto& Interaction()
+	{
+		return interaction_();
+	} 
+
+	bool loop();
+
+public:
+
+	TypeInfo("grainDEMSystem");
+
+	grainDEMSystem(
+		word  demSystemName,
+		const std::vector<box>& domains,
+		int argc, 
+		char* argv[],
+		bool requireRVel=false);
+
+	virtual ~grainDEMSystem();
+
+	grainDEMSystem(const grainDEMSystem&)=delete;
+
+	grainDEMSystem& operator = (const grainDEMSystem&)=delete;
+
+	add_vCtor(
+		DEMSystem,
+		grainDEMSystem,
+		word);
+
+
+	bool updateParticleDistribution(real extentFraction, const std::vector<box> domains) override;
+	 
+	int32 numParInDomain(int32 di)const override;
+	
+	std::vector<int32> numParInDomains()const override;
+
+	span<const int32> parIndexInDomain(int32 di)const override;
+
+	span<real> diameter() override;
+
+	span<real> courseGrainFactor() override;
+
+	span<realx3> acceleration() override;
+
+	span<realx3> velocity() override;
+
+	span<realx3> position() override;
+
+	span<realx3> rAcceleration() override;
+
+	span<realx3> rVelocity() override;
+
+	span<realx3> rPosition() override;
+
+	span<realx3> parFluidForce() override;
+
+	span<realx3> parFluidTorque() override;
+
+	bool sendFluidForceToDEM() override;
+
+	bool sendFluidTorqueToDEM() override;
+
+	bool beforeIteration() override;
+
+	bool iterate(
+		real upToTime, 
+		real timeToWrite, 
+		word timeName) override;
+
+	bool iterate(real upToTime) override;
+	
+	real maxBounndingSphereSize()const override;
+
+};
+
+
+} // pFlow
+
+#endif

DEMSystems/grainDEMSystem/grainFluidParticles.cpp

@@ -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.
+
+-----------------------------------------------------------------------------*/
+
+#include "grainFluidParticles.hpp"
+#include "grainFluidParticlesKernels.hpp"
+
+void pFlow::grainFluidParticles::checkHostMemory()
+{
+	if(fluidForce_.size()!=fluidForceHost_.size())
+	{
+		resizeNoInit(fluidForceHost_, fluidForce_.size());
+		resizeNoInit(fluidTorqueHost_, fluidTorque_.size());
+	}
+
+	// copy the data (if required) from device to host
+	courseGrainFactorHost_ = coarseGrainFactor().hostView();
+}
+
+pFlow::grainFluidParticles::grainFluidParticles(
+    systemControl &control,
+    const property &prop)
+    : grainParticles(control, prop),
+      fluidForce_(
+          objectFile(
+              "fluidForce",
+              "",
+              objectFile::READ_IF_PRESENT,
+              objectFile::WRITE_ALWAYS),
+          dynPointStruct(),
+          zero3),
+      fluidTorque_(
+          objectFile(
+              "fluidTorque",
+              "",
+              objectFile::READ_IF_PRESENT,
+              objectFile::WRITE_NEVER),
+          dynPointStruct(),
+          zero3)
+{
+	checkHostMemory();
+} 
+
+bool pFlow::grainFluidParticles::beforeIteration()
+{
+	grainParticles::beforeIteration();
+	checkHostMemory();
+	return true;
+}
+
+
+bool pFlow::grainFluidParticles::iterate() 
+{
+
+	accelerationTimer().start();
+		pFlow::grainFluidParticlesKernels::acceleration(
+			control().g(),
+			mass().deviceViewAll(),
+			contactForce().deviceViewAll(),
+			fluidForce_.deviceViewAll(),
+			I().deviceViewAll(),
+			contactTorque().deviceViewAll(),
+			fluidTorque_.deviceViewAll(),
+			pStruct().activePointsMaskDevice(),
+			accelertion().deviceViewAll(),
+			rAcceleration().deviceViewAll()
+			);
+	accelerationTimer().end();
+	
+	intCorrectTimer().start();
+	
+		dynPointStruct().correct(this->dt(), accelertion());
+	
+		rVelIntegration().correct(this->dt(), rVelocity(), rAcceleration());
+	
+	intCorrectTimer().end();
+	
+	return true;
+}
+
+void pFlow::grainFluidParticles::fluidForceHostUpdatedSync()
+{
+	copy(fluidForce_.deviceView(), fluidForceHost_);
+	return;
+}
+
+void pFlow::grainFluidParticles::fluidTorqueHostUpdatedSync()
+{
+	copy(fluidTorque_.deviceView(), fluidTorqueHost_);
+	return;
+}

DEMSystems/grainDEMSystem/grainFluidParticles.hpp

@@ -0,0 +1,105 @@
+/*------------------------------- 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::grainFluidParticles
+    
+   @brief Class for managing grain particles with fluid interactions 
+   
+   This is a top-level class that contains the essential components for
+   defining grain prticles with fluid interactions in a DEM simulation.
+ */
+
+#ifndef __grainFluidParticles_hpp__
+#define __grainFluidParticles_hpp__ 
+
+#include "grainParticles.hpp"
+
+
+namespace pFlow
+{
+
+class grainFluidParticles
+:
+	public grainParticles 
+{
+protected:
+	
+	/// pointField of rotational acceleration of particles on device  
+	realx3PointField_D 		fluidForce_;
+
+	hostViewType1D<realx3> 	fluidForceHost_;
+
+	realx3PointField_D 		fluidTorque_;
+
+	hostViewType1D<realx3>	fluidTorqueHost_;
+
+	hostViewType1D<real> 	courseGrainFactorHost_;
+
+
+	void checkHostMemory();
+
+public:
+
+	/// construct from systemControl and property 
+	grainFluidParticles(systemControl &control, const property& prop);
+
+	~grainFluidParticles() override = default;
+	
+	/// before iteration step 
+	bool beforeIteration() override;
+	
+	/// iterate particles 
+	bool iterate() override;	
+
+	auto& fluidForce()
+	{
+		return fluidForce_;
+	}
+
+	auto& fluidTorque()
+	{
+		return fluidTorque_;
+	}
+
+	auto& fluidForceHost()
+	{
+		return fluidForceHost_; 
+	}
+
+	auto& fluidTorqueHost()
+	{
+		return fluidTorqueHost_;
+	}
+
+	auto& courseGrainFactorHost()
+	{
+		return courseGrainFactorHost_;
+	}
+
+	void fluidForceHostUpdatedSync();
+
+	void fluidTorqueHostUpdatedSync();
+	
+}; //grainFluidParticles
+
+} // pFlow
+
+#endif //__sphereFluidParticles_hpp__

DEMSystems/grainDEMSystem/grainFluidParticlesKernels.hpp

@@ -0,0 +1,79 @@
+/*------------------------------- 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 __grainFluidParticlesKernels_hpp__
+#define __grainFluidParticlesKernels_hpp__
+
+namespace pFlow::grainFluidParticlesKernels
+{
+
+using rpAcceleration = Kokkos::RangePolicy<
+			DefaultExecutionSpace,
+			Kokkos::Schedule<Kokkos::Static>,
+			Kokkos::IndexType<int32>
+			>;
+
+template<typename IncludeFunctionType>
+void acceleration( 
+	realx3		g,
+	deviceViewType1D<real>  	mass,
+	deviceViewType1D<realx3>  	force,
+	deviceViewType1D<realx3>  	fluidForce,
+	deviceViewType1D<real>  	I,
+	deviceViewType1D<realx3>  	torque,
+	deviceViewType1D<realx3>  	fluidTorque,
+	IncludeFunctionType 		incld,
+	deviceViewType1D<realx3> 	lAcc,
+	deviceViewType1D<realx3> 	rAcc
+	)
+{
+
+	auto activeRange = incld.activeRange();
+	if(incld.isAllActive())
+	{
+		Kokkos::parallel_for(
+		"pFlow::grainParticlesKernels::acceleration",
+		rpAcceleration(activeRange.first, activeRange.second),
+		LAMBDA_HD(int32 i){
+				lAcc[i] = (force[i]+fluidForce[i])/mass[i] + g;
+				rAcc[i] = (torque[i]+fluidTorque[i])/I[i];
+		});
+	}
+	else
+	{
+		Kokkos::parallel_for(
+		"pFlow::grainParticlesKernels::acceleration",
+		rpAcceleration(activeRange.first, activeRange.second),
+		LAMBDA_HD(int32 i){
+			if(incld(i))
+			{
+				lAcc[i] = (force[i]+fluidForce[i])/mass[i] + g;
+				rAcc[i] = (torque[i]+fluidTorque[i])/I[i];
+			}
+		});
+
+	}
+	
+	Kokkos::fence();
+}
+
+}
+
+#endif 

DEMSystems/sphereDEMSystem/sphereDEMSystem.cpp

@@ -19,6 +19,7 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "sphereDEMSystem.hpp"
+#include "vocabs.hpp"
 
 bool pFlow::sphereDEMSystem::loop()
 {
@@ -26,17 +27,16 @@ bool pFlow::sphereDEMSystem::loop()
 	do 
 	{
 		//
-		
 		if(! insertion_().insertParticles( 
+			Control().time().currentIter(),
 			Control().time().currentTime(),
 			Control().time().dt()	) )
 		{
 			fatalError<<
-			"particle insertion failed in sphereDFlow solver.\n";
+			"particle insertion failed in sphereDEMSystem.\n";
 			return false;
 		}	
 
-
 		geometry_->beforeIteration();
 		
 		interaction_->beforeIteration();
@@ -63,29 +63,31 @@ pFlow::sphereDEMSystem::sphereDEMSystem(
 		word  demSystemName,
 		const std::vector<box>& domains,
 		int argc, 
-		char* argv[])
+		char* argv[],
+		bool requireRVel)
 :
-	DEMSystem(demSystemName, domains, argc, argv)
+	DEMSystem(demSystemName, domains, argc, argv),
+	requireRVel_(requireRVel)
 { 
 	
-	REPORT(0)<<"\nReading proprties . . . "<<endREPORT;
+	REPORT(0)<<"\nReading proprties . . . "<<END_REPORT;
 	property_ = makeUnique<property>(
-    	Control().caseSetup().path()+propertyFile__);
+		propertyFile__,
+    	Control().caseSetup().path());
 
-	REPORT(0)<< "\nCreating surface geometry for sphereDEMSystem . . . "<<endREPORT;
+	REPORT(0)<< "\nCreating surface geometry for sphereDEMSystem . . . "<<END_REPORT;
 	geometry_ = geometry::create(Control(), Property());
 
 
-	REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
+	REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
 	particles_ = makeUnique<sphereFluidParticles>(Control(), Property());
 
 
 	insertion_ = makeUnique<sphereInsertion>( 
-		Control().caseSetup().path()+insertionFile__, 
 		particles_(), 
-		particles_().shapes());
+		particles_().spheres());
 	
-	REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<endREPORT;
+	REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<END_REPORT;
 	interaction_ = interaction::create(
 		Control(),
 		Particles(),
@@ -98,6 +100,7 @@ pFlow::sphereDEMSystem::sphereDEMSystem(
 
 }
 
+
 pFlow::sphereDEMSystem::~sphereDEMSystem()
 {
 	
@@ -119,8 +122,8 @@ bool pFlow::sphereDEMSystem::updateParticleDistribution(
 	}
 
 	if(!particleDistribution_->locateParticles(
-			parPosition_, 
-			particles_->pStruct().activePointsMaskH()))
+			positionHost_, 
+			particles_->pStruct().activePointsMaskHost()))
 	{
 		fatalErrorInFunction<<
 		"error in locating particles among sub-domains"<<endl;
@@ -137,7 +140,7 @@ pFlow::int32
 }
 
 std::vector<pFlow::int32> 
-	pFlow::sphereDEMSystem::numParInDomain()const
+	pFlow::sphereDEMSystem::numParInDomains()const
 {
 	const auto& distribute = particleDistribution_();
 	int32 numDomains = distribute.numDomains();
@@ -156,31 +159,56 @@ pFlow::sphereDEMSystem::parIndexInDomain(int32 di)const
 	return particleDistribution_->particlesInDomain(di);
 }
 
-pFlow::span<pFlow::real> pFlow::sphereDEMSystem::parDiameter() 
+pFlow::span<pFlow::real> pFlow::sphereDEMSystem::diameter() 
 {	
-	return span<real>(parDiameter_.data(), parDiameter_.size());
+	return span<real>(diameterHost_.data(), diameterHost_.size());
 }
 
-pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parVelocity()  
+pFlow::span<pFlow::real> pFlow::sphereDEMSystem::courseGrainFactor()
 {
-	return span<realx3>(parVelocity_.data(), parVelocity_.size());
+    return span<real>(nullptr, 0);
 }
 
-pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parPosition()  
+pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::acceleration()
 {
-	return span<realx3>(parPosition_.data(), parPosition_.size());
+    return span<realx3>(nullptr, 0);
+}
+
+pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::velocity()  
+{
+	return span<realx3>(velocityHost_.data(), velocityHost_.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::position()  
+{
+	return span<realx3>(positionHost_.data(), positionHost_.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::rAcceleration()
+{
+    return span<realx3>(nullptr, 0);
+}
+
+pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::rVelocity()
+{
+    return span<realx3>(rVelocityHost_.data(), rVelocityHost_.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::rPosition()
+{
+    return span<realx3>(nullptr, 0);
 }
 
 pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parFluidForce() 
 {
-	auto& hVec =  particles_->fluidForceHostAll();
+	auto& hVec =  particles_->fluidForceHost();
 	
 	return span<realx3>(hVec.data(), hVec.size());
 }
 
 pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parFluidTorque() 
 {
-	auto& hVec =  particles_->fluidTorqueHostAll();
+	auto& hVec =  particles_->fluidTorqueHost();
 	return span<realx3>(hVec.data(), hVec.size());
 }
 
@@ -198,9 +226,14 @@ bool pFlow::sphereDEMSystem::sendFluidTorqueToDEM()
 
 bool pFlow::sphereDEMSystem::beforeIteration()
 {
-	parVelocity_ = particles_->dynPointStruct().velocityHostAll();
-	parPosition_ = particles_->dynPointStruct().pointPositionHostAll();
-	parDiameter_ = particles_->diameter().hostVectorAll();
+	velocityHost_ = std::as_const(particles_()).velocity().hostView();
+	positionHost_ = std::as_const(particles_()).pointPosition().hostView();
+	diameterHost_ = particles_->diameter().hostView();
+
+	if(requireRVel_)
+		rVelocityHost_ = std::as_const(particles_()).rVelocity().hostView();
+	
+
 	return true;
 }
 

DEMSystems/sphereDEMSystem/sphereDEMSystem.hpp

@@ -42,24 +42,29 @@ protected:
 	
 	// protected members 
 
-	uniquePtr<property> 		property_ = nullptr;
+	uniquePtr<property> 			property_ = nullptr;
 
-	uniquePtr<geometry> 		geometry_ = nullptr;
+	uniquePtr<geometry> 			geometry_ = nullptr;
 
-	uniquePtr<sphereFluidParticles> 	particles_ = nullptr;
+	uniquePtr<sphereFluidParticles> particles_ = nullptr;
 
-	uniquePtr<sphereInsertion> 	insertion_ = nullptr;
+	uniquePtr<sphereInsertion> 		insertion_ = nullptr;
 
-	uniquePtr<interaction> 		interaction_ = nullptr;
+	uniquePtr<interaction> 			interaction_ = nullptr;
 
-	uniquePtr<domainDistribute> particleDistribution_=nullptr;
+	uniquePtr<domainDistribute> 	particleDistribution_=nullptr;
 
 	// to be used for CPU communications 
-	ViewType1D<realx3, HostSpace>  parVelocity_;
+	ViewType1D<realx3, HostSpace> 	velocityHost_;
 
-	ViewType1D<realx3, HostSpace>  parPosition_;
+	ViewType1D<realx3, HostSpace> 	positionHost_;
+
+	ViewType1D<real, HostSpace> 	diameterHost_;
+
+	bool 							requireRVel_ = false;
+
+	ViewType1D<realx3, HostSpace> 	rVelocityHost_;
 
-	ViewType1D<real, HostSpace> parDiameter_;
 
 //  protected member functions
 	auto& Property()
@@ -92,7 +97,8 @@ public:
 		word  demSystemName,
 		const std::vector<box>& domains,
 		int argc, 
-		char* argv[]);
+		char* argv[],
+		bool requireRVel=false);
 
 	virtual ~sphereDEMSystem();
 
@@ -110,15 +116,25 @@ public:
 	 
 	int32 numParInDomain(int32 di)const override;
 	
-	std::vector<int32> numParInDomain()const override;
+	std::vector<int32> numParInDomains()const override;
 
 	span<const int32> parIndexInDomain(int32 di)const override;
 
-	span<real> parDiameter() override;
+	span<real> diameter() override;
 
-	span<realx3> parVelocity() override;
+	span<real> courseGrainFactor() override;
 
-	span<realx3> parPosition() override;
+	span<realx3> acceleration() override;
+
+	span<realx3> velocity() override;
+
+	span<realx3> position() override;
+
+	span<realx3> rAcceleration() override;
+
+	span<realx3> rVelocity() override;
+
+	span<realx3> rPosition() override;
 
 	span<realx3> parFluidForce() override;
 

DEMSystems/sphereDEMSystem/sphereFluidParticles.cpp

@@ -21,42 +21,44 @@ Licence:
 #include "sphereFluidParticles.hpp"
 #include "sphereFluidParticlesKernels.hpp"
 
+void pFlow::sphereFluidParticles::checkHostMemory()
+{
+	if(fluidForce_.size()!=fluidForceHost_.size())
+	{
+		resizeNoInit(fluidForceHost_, fluidForce_.size());
+		resizeNoInit(fluidTorqueHost_, fluidTorque_.size());
+	}
+}
 
 pFlow::sphereFluidParticles::sphereFluidParticles(
-	systemControl &control,
-	const property& prop
-)
-:
-	sphereParticles(control, prop),
-	fluidForce_(
-		this->time().emplaceObject<realx3PointField_HD>(
-			objectFile(
-				"fluidForce",
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS
-				),
-			pStruct(),
-			zero3
-			)
-		),
-	fluidTorque_(
-		this->time().emplaceObject<realx3PointField_HD>(
-			objectFile(
-				"fluidTorque",
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS
-				),
-			pStruct(),
-			zero3
-			)
-		)
-{} 
+    systemControl &control,
+    const property &prop)
+    : sphereParticles(control, prop),
+      fluidForce_(
+          objectFile(
+              "fluidForce",
+              "",
+              objectFile::READ_IF_PRESENT,
+              objectFile::WRITE_ALWAYS),
+          dynPointStruct(),
+          zero3),
+      fluidTorque_(
+          objectFile(
+              "fluidTorque",
+              "",
+              objectFile::READ_IF_PRESENT,
+              objectFile::WRITE_NEVER),
+          dynPointStruct(),
+          zero3)
+{
+	checkHostMemory();
+} 
 
 bool pFlow::sphereFluidParticles::beforeIteration()
 {
 	sphereParticles::beforeIteration();
+	checkHostMemory();
+
 	return true;
 }
 
@@ -64,42 +66,40 @@ bool pFlow::sphereFluidParticles::beforeIteration()
 bool pFlow::sphereFluidParticles::iterate() 
 {
 
-	accelerationTimer_.start();
+	accelerationTimer().start();
 		pFlow::sphereFluidParticlesKernels::acceleration(
 			control().g(),
-			mass().deviceVectorAll(),
-			contactForce().deviceVectorAll(),
-			fluidForce().deviceVectorAll(),
-			I().deviceVectorAll(),
-			contactTorque().deviceVectorAll(),
-			fluidTorque().deviceVectorAll(),
-			pStruct().activePointsMaskD(),
-			accelertion().deviceVectorAll(),
-			rAcceleration().deviceVectorAll()
+			mass().deviceViewAll(),
+			contactForce().deviceViewAll(),
+			fluidForce_.deviceViewAll(),
+			I().deviceViewAll(),
+			contactTorque().deviceViewAll(),
+			fluidTorque_.deviceViewAll(),
+			pStruct().activePointsMaskDevice(),
+			accelertion().deviceViewAll(),
+			rAcceleration().deviceViewAll()
 			);
-	accelerationTimer_.end();
+	accelerationTimer().end();
 	
-	intCorrectTimer_.start();
+	intCorrectTimer().start();
 	
-		dynPointStruct_.correct(this->dt(), accelertion_);
+		dynPointStruct().correct(this->dt(), accelertion());
 	
-		rVelIntegration_().correct(this->dt(), rVelocity_, rAcceleration_);
+		rVelIntegration().correct(this->dt(), rVelocity(), rAcceleration());
 	
-	intCorrectTimer_.end();
+	intCorrectTimer().end();
 	
 	return true;
 }
 
 void pFlow::sphereFluidParticles::fluidForceHostUpdatedSync()
 {
-	fluidForce_.modifyOnHost();
-	fluidForce_.syncViews();
+	copy(fluidForce_.deviceView(), fluidForceHost_);
 	return;
 }
 
 void pFlow::sphereFluidParticles::fluidTorqueHostUpdatedSync()
 {
-	fluidTorque_.modifyOnHost();
-	fluidTorque_.syncViews();
+	copy(fluidTorque_.deviceView(), fluidTorqueHost_);
 	return;
 }

DEMSystems/sphereDEMSystem/sphereFluidParticles.hpp

@@ -43,12 +43,17 @@ class sphereFluidParticles
 protected:
 	
 	/// pointField of rotational acceleration of particles on device  
-	realx3PointField_HD& 	fluidForce_;
+	realx3PointField_D 		fluidForce_;
 
-	realx3PointField_HD& 	fluidTorque_;
+	hostViewType1D<realx3> 	fluidForceHost_;
 
+	realx3PointField_D 		fluidTorque_;
 
-	void zeroFluidForce_H()
+	hostViewType1D<realx3>	fluidTorqueHost_;
+
+	void checkHostMemory();
+
+	/*void zeroFluidForce_H()
 	{
 		fluidForce_.fillHost(zero3);
 	}
@@ -56,7 +61,7 @@ protected:
 	void zeroFluidTorque_H()
 	{
 		fluidTorque_.fillHost(zero3);
-	}
+	}*/
 
 public:
 
@@ -81,17 +86,16 @@ public:
 	}
 
 
-	auto& fluidForceHostAll()
+	auto& fluidForceHost()
 	{
-		return fluidForce_.hostVectorAll();
+		return fluidForceHost_; 
 	}
 
-	auto& fluidTorqueHostAll()
+	auto& fluidTorqueHost()
 	{
-		return fluidTorque_.hostVectorAll();
+		return fluidTorqueHost_;
 	}
 
-
 	void fluidForceHostUpdatedSync();
 
 	void fluidTorqueHostUpdatedSync();

DEMSystems/sphereDEMSystem/sphereFluidParticlesKernels.hpp

@@ -46,7 +46,7 @@ void acceleration(
 {
 
 	auto activeRange = incld.activeRange();
-	if(incld.allActive())
+	if(incld.isAllActive())
 	{
 		Kokkos::parallel_for(
 		"pFlow::sphereParticlesKernels::acceleration",

README.md

@@ -3,10 +3,18 @@
 </div>
 
 
-**PhasicFlow** is a parallel C++ code for performing DEM simulations. It can run on shared-memory multi-core computational units such as multi-core CPUs or GPUs (for now it works on CUDA-enabled GPUs). The parallelization method mainly relies on loop-level parallelization on a shared-memory computational unit. You can build and run PhasicFlow in serial mode on regular PCs, in parallel mode for multi-core CPUs, or build it for a GPU device to off-load computations to a GPU. In its current statues you can simulate millions of particles (up to 32M particles tested) on a single desktop computer. You can see the [performance tests of PhasicFlow](https://github.com/PhasicFlow/phasicFlow/wiki/Performance-of-phasicFlow) in the wiki page. 
+**PhasicFlow** is a parallel C++ code for performing DEM simulations. It can run on shared-memory multi-core computational units such as multi-core CPUs or GPUs (for now it works on CUDA-enabled GPUs). The parallelization method mainly relies on loop-level parallelization on a shared-memory computational unit. You can build and run PhasicFlow in serial mode on regular PCs, in parallel mode for multi-core CPUs, or build it for a GPU device to off-load computations to a GPU. In its current statues you can simulate millions of particles (up to 80M particles tested) on a single desktop computer. You can see the [performance tests of PhasicFlow](https://github.com/PhasicFlow/phasicFlow/wiki/Performance-of-phasicFlow) in the wiki page.
 
+**MPI** parallelization with dynamic load balancing is under development. With this level of parallelization, PhasicFlow can leverage the computational power of **multi-gpu** workstations or clusters with distributed memory CPUs. 
+In summary PhasicFlow can have 6 execution modes:
+1. Serial on a single CPU core,
+2. Parallel on a multi-core computer/node (using OpenMP),
+3. Parallel on an nvidia-GPU (using Cuda),
+4. Parallel on distributed memory workstation (Using MPI)
+5. Parallel on distributed memory workstations with multi-core nodes (using MPI+OpenMP)
+6. Parallel on workstations with multiple GPUs (using MPI+Cuda).
 ## How to build?
-You can build PhasicFlow for CPU and GPU executions. [Here is a complete step-by-step procedure](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-Build-PhasicFlow).
+You can build PhasicFlow for CPU and GPU executions. The latest release of PhasicFlow is v-0.1. [Here is a complete step-by-step procedure for building phasicFlow-v-0.1.](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-Build-PhasicFlow).
 
 ## Online code documentation
 You can find a full documentation of the code, its features, and other related materials on [online documentation of the code](https://phasicflow.github.io/phasicFlow/)
@@ -22,3 +30,19 @@ PhasicFlowPlus is and extension to PhasicFlow for simulating particle-fluid syst
 * [Kokkos](https://github.com/kokkos/kokkos) from National Technology & Engineering Solutions of Sandia, LLC (NTESS)
 * [CLI11 1.8](https://github.com/CLIUtils/CLI11) from University of Cincinnati.
 
+## How to cite PhasicFlow
+If you are using PhasicFlow in your research or industrial work, cite the following [article](https://www.sciencedirect.com/science/article/pii/S0010465523001662):
+```
+@article{NOROUZI2023108821,
+title = {PhasicFlow: A parallel, multi-architecture open-source code for DEM simulations},
+journal = {Computer Physics Communications},
+volume = {291},
+pages = {108821},
+year = {2023},
+issn = {0010-4655},
+doi = {https://doi.org/10.1016/j.cpc.2023.108821},
+url = {https://www.sciencedirect.com/science/article/pii/S0010465523001662},
+author = {H.R. Norouzi},
+keywords = {Discrete element method, Parallel computing, CUDA, GPU, OpenMP, Granular flow}
+}
+```

cmake/autoComplete

@@ -0,0 +1,33 @@
+PF_cFlags="--description --help --version"
+
+_pFlowToVTK(){
+    if [ "$3" == "--time" ]; then
+      COMPREPLY=( $(ls) )
+    else
+      COMPREPLY=( $(compgen -W "$PF_cFlags --binary --no-geometry --no-particles --out-folder --time --separate-surfaces --fields" -- "$2") )
+    fi
+}
+complete -F _pFlowToVTK pFlowToVTK
+
+_postprocessPhasicFlow(){
+    if [ "$3" == "--time" ]; then
+      COMPREPLY=( $(ls) )
+    else
+      COMPREPLY=( $(compgen -W "$PF_cFlags --out-folder --time --zeroFolder" -- "$2") )
+    fi
+}
+complete -F _postprocessPhasicFlow postprocessPhasicFlow
+
+complete -W "$PF_cFlags --positionParticles-only --setFields-only --coupling" particlesPhasicFlow
+
+complete -W "$PF_cFlags --coupling" geometryPhasicFlow 
+
+complete -W "$PF_cFlags --coupling" iterateGeometry
+
+complete -W "$PF_cFlags" iterateSphereParticles
+
+complete -W "$PF_cFlags" sphereGranFlow
+
+complete -W "$PF_cFlags" grainGranFlow
+
+complete -W "$PF_cFlags" checkPhasicFlow

cmake/bashrc

@@ -1,8 +1,7 @@
 
-export pFlow_PROJECT_VERSION=v0.1
-
-export pFlow_PROJECT=phasicFlow
+export pFlow_PROJECT_VERSION=v-1.0
 
+export pFlow_PROJECT="phasicFlow-$pFlow_PROJECT_VERSION"
 
 projectDir="$HOME/PhasicFlow" 
 kokkosDir="$HOME/Kokkos/kokkos"
@@ -20,6 +19,8 @@ export pFlow_SRC_DIR="$pFlow_PROJECT_DIR/src"
 
 export Kokkos_DIR="$kokkosDir"
 
+export Zoltan_DIR="$projectDir/Zoltan"
+
 # Cleanup variables (done as final statement for a clean exit code)
 unset projectDir
 
@@ -27,5 +28,7 @@ export PATH="$pFlow_BIN_DIR:$PATH"
 
 export LD_LIBRARY_PATH="$pFlow_LIB_DIR:$LD_LIBRARY_PATH"
 
+export LD_LIBRARY_PATH="$Zoltan_DIR/lib:$LD_LIBRARY_PATH"
 
+source $pFlow_PROJECT_DIR/cmake/autoComplete 
 #------------------------------------------------------------------------------

cmake/globals.cmake

@@ -2,3 +2,52 @@
 set(validFiles)
 list(APPEND validFiles *.C *.cpp *.cxx *.c *.cu *.H *.hpp *.hxx *.h *.cuh)
 
+macro(Kokkos_cmake_settings)
+
+#mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LAMBDA)
+mark_as_advanced(FORCE var Kokkos_CXX_STANDARD)
+#mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_CONSTEXPR)
+mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMP)
+mark_as_advanced(FORCE var Kokkos_ENABLE_SERIAL)
+mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA)
+mark_as_advanced(FORCE var Kokkos_ENABLE_HIP)
+mark_as_advanced(FORCE var Kokkos_ENABLE_AGGRESSIVE_VECTORIZATION)
+mark_as_advanced(FORCE var Kokkos_ENABLE_BENCHMARKS)
+mark_as_advanced(FORCE var Kokkos_ENABLE_COMPILE_AS_CMAKE_LANGUAGE)
+mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_LDG_INTRINSIC)
+mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_RELOCATABLE_DEVICE_CODE)
+mark_as_advanced(FORCE var Kokkos_ENABLE_CUDA_UVM)
+mark_as_advanced(FORCE var Kokkos_ENABLE_DEPRECATED_CODE_3)
+mark_as_advanced(FORCE var Kokkos_ENABLE_DEPRECATED_CODE_4)
+mark_as_advanced(FORCE var Kokkos_ENABLE_DEPRECATION_WARNINGS)
+mark_as_advanced(FORCE var Kokkos_ENABLE_DESUL_ATOMICS_EXTERNAL)
+mark_as_advanced(FORCE var Kokkos_ENABLE_EXAMPLES)
+mark_as_advanced(FORCE var Kokkos_ENABLE_HEADER_SELF_CONTAINMENT_TESTS)
+mark_as_advanced(FORCE var Kokkos_ENABLE_HIP_MULTIPLE_KERNEL_INSTANTIATIONS)
+mark_as_advanced(FORCE var Kokkos_ENABLE_HIP_RELOCATABLE_DEVICE_CODE)
+mark_as_advanced(FORCE var Kokkos_ENABLE_HPX)
+mark_as_advanced(FORCE var Kokkos_ENABLE_HWLOC)
+mark_as_advanced(FORCE var Kokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC)
+mark_as_advanced(FORCE var Kokkos_ENABLE_IMPL_HPX_ASYNC_DISPATCH)
+mark_as_advanced(FORCE var Kokkos_ENABLE_LARGE_MEM_TESTS)
+mark_as_advanced(FORCE var Kokkos_ENABLE_DEBUG_DUALVIEW_MODIFY_CHECK)
+mark_as_advanced(FORCE var Kokkos_ENABLE_LIBQUADMATH)
+mark_as_advanced(FORCE var Kokkos_ENABLE_LIBRT)
+mark_as_advanced(FORCE var Kokkos_ENABLE_MEMKIND)
+mark_as_advanced(FORCE var Kokkos_ENABLE_ONEDPL)
+mark_as_advanced(FORCE var Kokkos_ENABLE_OPENACC)
+mark_as_advanced(FORCE var Kokkos_ENABLE_OPENMPTARGET)
+mark_as_advanced(FORCE var Kokkos_ENABLE_ROCM)
+mark_as_advanced(FORCE var Kokkos_ENABLE_SYCL)
+mark_as_advanced(FORCE var Kokkos_ENABLE_TESTS)
+mark_as_advanced(FORCE var Kokkos_ENABLE_THREADS)
+mark_as_advanced(FORCE var Kokkos_ENABLE_TUNING)
+mark_as_advanced(FORCE var Kokkos_ENABLE_UNSUPPORTED_ARCHS)
+mark_as_advanced(FORCE var Kokkos_HPX_DIR)
+mark_as_advanced(FORCE var Kokkos_HWLOC_DIR)
+mark_as_advanced(FORCE var Kokkos_MEMKIND_DIR)
+mark_as_advanced(FORCE var Kokkos_ROCM_DIR)
+mark_as_advanced(FORCE var Kokkos_THREADS_DIR)
+
+endmacro()
+

cmake/makeExecutableGlobals.cmake

@@ -31,8 +31,8 @@ target_include_directories(${target_name}
 
 message(STATUS "\nCreating make file for executable ${target_name}")
 message(STATUS "  ${target_name} link libraries are: ${${target_link_libs}}")
-message(STATUS "  ${target_name} source files are: ${source_files}")
-message(STATUS "  ${target_name} include dirs are: ${includeDirs}\n")
+#message(STATUS "  ${target_name} source files are: ${${source_files}}")
+#message(STATUS "  ${target_name} include dirs are: ${includeDirs}\n")
 
 
 install(TARGETS ${target_name} DESTINATION bin)

cmake/makeLibraryGlobals.cmake

@@ -41,8 +41,8 @@ target_include_directories(${target_name}
 
 message(STATUS "\nCreating make file for library ${target_name}")
 message(STATUS "  ${target_name} link libraries are: ${${target_link_libs}}")
-message(STATUS "  ${target_name} source files are: ${source_files}")
-message(STATUS "  ${target_name} include dirs are: ${includeDirs}\n")
+#message(STATUS "  ${target_name} source files are: ${source_files}")
+#message(STATUS "  ${target_name} include dirs are: ${includeDirs}\n \n")
 
 install(TARGETS ${target_name} DESTINATION lib)
 install(FILES ${includeFiles}  DESTINATION include/${target_name})

phasicFlowConfig.H.in

@@ -5,5 +5,6 @@
 #cmakedefine pFlow_Build_Serial
 #cmakedefine pFlow_Build_OpenMP
 #cmakedefine pFlow_Build_Cuda
-#cmakedefine USE_STD_PARALLEL_ALG
+#cmakedefine pFlow_STD_Parallel_Alg
 #cmakedefine pFlow_Build_Double
+#cmakedefine pFlow_Build_MPI

solvers/CMakeLists.txt

@@ -1,7 +1,9 @@
 
 
-#add_subdirectory(iterateSphereParticles)
+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();
+
+
+}	
+

solvers/iterateGeometry/iterateGeometry.cpp

@@ -29,17 +29,16 @@ Licence:
  * folder.
  */
 
+#include "vocabs.hpp"
 #include "systemControl.hpp"
-#include "geometryMotion.hpp"
+#include "geometry.hpp"
 #include "commandLine.hpp"
-#include "readControlDict.hpp"
-
-using pFlow::commandLine;
+//#include "readControlDict.hpp"
 
 int main( int argc, char* argv[] )
 {
 
-commandLine cmds(
+pFlow::commandLine cmds(
     "iterateGeometry",
     "Performs simulation without particles, only geometry is solved");
 
@@ -54,6 +53,8 @@ commandLine cmds(
    
 
 // this should be palced in each main 
+pFlow::processors::initProcessors(argc, argv);
+pFlow::initialize_pFlowProcessors();
 #include "initialize_Control.hpp"
 	
 	#include "setProperty.hpp"
@@ -68,6 +69,7 @@ commandLine cmds(
 
 // this should be palced in each main 
 #include "finalize.hpp"
+pFlow::processors::finalizeProcessors();
 
 }
 

solvers/iterateSphereParticles/CMakeLists.txt

@@ -0,0 +1,7 @@
+
+set(source_files
+iterateSphereParticles.cpp
+)
+set(link_lib  Particles)
+
+pFlow_make_executable_install(iterateSphereParticles source_files link_lib)

solvers/iterateSphereParticles/createDEMComponents.hpp

@@ -18,10 +18,8 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+//
+REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
+pFlow::sphereParticles sphParticles(Control, proprties);
 
-#include "peakableRegions.hpp"
-
- 
-#ifdef BUILD_SHARED_LIBS
-	#include "peakableRegionInstantiate.cpp"
-#endif
+WARNING<<"Particle insertion has not been set yet!"<<END_WARNING;

solvers/iterateSphereParticles/iterateSphereParticles.cpp

@@ -0,0 +1,89 @@
+
+/*------------------------------- 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 "property.hpp"
+#include "sphereParticles.hpp"
+#include "systemControl.hpp"
+#include "commandLine.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(
+		"sphereGranFlow",
+		"DEM solver for non-cohesive spherical 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 "createDEMComponents.hpp"
+	
+	REPORT(0)<<"\nStart of time loop . . .\n"<<END_REPORT;
+
+	do 
+	{
+		
+		sphParticles.beforeIteration();
+	
+		sphParticles.iterate();
+
+		sphParticles.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();
+
+
+}	
+

solvers/sphereGranFlow/createDEMComponents.hpp

@@ -19,11 +19,11 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 //
-REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
-sphereParticles sphParticles(Control, proprties);
+REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
+pFlow::sphereParticles sphParticles(Control, proprties);
 
 //
-REPORT(0)<<"\nCreating particle insertion object . . ."<<endREPORT;
+REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
 /*auto& sphInsertion = 
 	Control.caseSetup().emplaceObject<sphereInsertion>(
 		objectFile(
@@ -36,13 +36,12 @@ REPORT(0)<<"\nCreating particle insertion object . . ."<<endREPORT;
 		sphParticles.shapes()
 	);*/
 
-auto sphInsertion =  sphereInsertion(
-	Control.caseSetup().path()+insertionFile__, 
+auto sphInsertion =  pFlow::sphereInsertion(
 	sphParticles, 
-	sphParticles.shapes());
+	sphParticles.spheres());
 
-REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<endREPORT;
-auto interactionPtr = interaction::create(
+REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<END_REPORT;
+auto interactionPtr = pFlow::interaction::create(
 	Control,
 	sphParticles,
 	surfGeometry

solvers/sphereGranFlow/sphereGranFlow.cpp

@@ -30,29 +30,17 @@ Licence:
  * (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 "sphereParticles.hpp"
+#include "interaction.hpp"
 #include "Insertions.hpp"
-#include "systemControl.hpp"
-#include "contactSearch.hpp"
-#include "sphereInteraction.hpp"
-#include "commandLine.hpp"
-#include "readControlDict.hpp"
-
-using pFlow::output;
-using pFlow::endl;
-using pFlow::property;
-using pFlow::sphereParticles;
-using pFlow::objectFile;
-using pFlow::sphereInsertion;
-using pFlow::insertionFile__;
-using pFlow::interactionFile__;
-using pFlow::contactSearch;
-using pFlow::interaction;
-using pFlow::commandLine;
-
 
+//#include "readControlDict.hpp"
 
 /**
  * DEM solver for simulating granular flow of cohesion-less particles.
@@ -63,7 +51,7 @@ using pFlow::commandLine;
 int main( int argc, char* argv[])
 {
 
-commandLine cmds(
+pFlow::commandLine cmds(
 		"sphereGranFlow",
 		"DEM solver for non-cohesive spherical particles with particle insertion "
 		"mechanism and moving geometry");
@@ -73,20 +61,22 @@ 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"<<endREPORT;
+	REPORT(0)<<"\nStart of time loop . . .\n"<<END_REPORT;
 
 	do 
 	{
 		
 		if(! sphInsertion.insertParticles( 
+			Control.time().currentIter(),
 			Control.time().currentTime(),
 			Control.time().dt()	) )
 		{
@@ -95,30 +85,34 @@ if(!cmds.parse(argc, argv)) return 0;
 			return 1;
 		}
 		 
+		// set force to zero
 		surfGeometry.beforeIteration();
 
+		// set force to zero, predict, particle deletion and etc. 
+		sphParticles.beforeIteration();
+
 		sphInteraction.beforeIteration();
 		
-		sphParticles.beforeIteration();
-		
-		
 		sphInteraction.iterate();	
 		
-		sphParticles.iterate();
-
 		surfGeometry.iterate();
 
-		sphParticles.afterIteration();
+		sphParticles.iterate();
+
+		sphInteraction.afterIteration();
 
 		surfGeometry.afterIteration();
 
+		sphParticles.afterIteration();
+				
 
 	}while(Control++);
 
-	REPORT(0)<<"\nEnd of time loop.\n"<<endREPORT;
+	REPORT(0)<<"\nEnd of time loop.\n"<<END_REPORT;
 
 // this should be palced in each main 
 #include "finalize.hpp"
+pFlow::processors::finalizeProcessors();
 
 
 }	

src/CMakeLists.txt

@@ -7,9 +7,9 @@ add_subdirectory(Property)
 
 add_subdirectory(Particles)
 
+add_subdirectory(Geometry)
+
 add_subdirectory(Interaction)
 
 add_subdirectory(MotionModel)
 
-add_subdirectory(Geometry)
-

src/Geometry/geometry/geometry.cpp

@@ -19,32 +19,40 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "geometry.hpp"
+#include "systemControl.hpp"
 #include "vocabs.hpp"
 
 
-bool pFlow::geometry::findPropertyId()
+bool pFlow::geometry::createPropertyId()
 {
-
-	int8Vector propId(0, surface().capacity(),RESERVE());
-	propId.clear();
-
-	uint32 pId;
-	ForAll(matI, materialName_)
+	if(materialName_.size() != numSurfaces() )
 	{
+		fatalErrorInFunction<<
+		"number of subSurface and material names do not match"<<endl;
+		return false;
+	}
 
-		if( !wallProperty_.nameToIndex( materialName_[matI], pId ) )
+	uint32Vector propId(
+		"propId", 
+		capacity(), 
+		size(),
+		RESERVE());
+
+	ForAll(i, materialName_)
+	{
+		uint32 pIdx =0;
+		
+		if( !wallProperty_.nameToIndex(materialName_[i], pIdx) )
 		{
-		 	fatalErrorInFunction<<
-		 	"material name for the geometry is invalid: "<< materialName_[matI]<<endl;
-		 	return false;
+			fatalErrorInFunction<<
+			"Property/material name is invalid "<<materialName_[i]<<
+			". A list of valid names are \n"<< wallProperty_.materials()<<endl;
+			return false;
 		}
 		
-		int32 surfSize = surface().surfNumTriangles(matI);
+		auto triRange = subSurfaceRange(i);
+		propId.fill(triRange.start(), triRange.end(), pIdx);
 		
-		for(int32 i=0; i<surfSize; i++)
-		{
-			propId.push_back(pId);
-		}
 	}
 
 	propertyId_.assign(propId);
@@ -53,83 +61,363 @@ bool pFlow::geometry::findPropertyId()
 
 }
 
+void pFlow::geometry::zeroForce()
+{
+	contactForceWall_.fill(zero3);
+}
+
 pFlow::geometry::geometry
 (
 	systemControl& control,
 	const property& prop
 )
 :
-	demGeometry(control),
+	multiTriSurface
+	(
+		objectFile
+		(
+			triSurfaceFile__,
+			"",
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_ALWAYS
+		),
+		&control.geometry()
+	),
+	demComponent
+	(
+		"geometry",
+		control
+	),
 	wallProperty_(prop),
-	geometryRepository_(control.geometry()),
-	triSurface_(
-		control.geometry().emplaceObject<multiTriSurface>(
-			objectFile(
-				"triSurface",
-				"",
-				objectFile::READ_ALWAYS,
-				objectFile::WRITE_ALWAYS
-				)
-			)
+	propertyId_
+	(
+		objectFile
+		(
+			"propertyId",
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_NEVER
 		),
-	motionComponentName_(
-		control.geometry().emplaceObject<wordField>(
-			objectFile(
-				"motionComponentName",
-				"",
-				objectFile::READ_ALWAYS,
-				objectFile::WRITE_ALWAYS	
-				),
-			"motionNamesList"
-			)
+		*this,
+		0u
+	),
+	contactForceWall_
+	(
+		objectFile
+		(
+			"contactForcWall",
+			"",
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_NEVER
 		),
-	materialName_(
-		control.geometry().emplaceObject<wordField>(
-			objectFile(
-				"materialName",
-				"",
-				objectFile::READ_ALWAYS,
-				objectFile::WRITE_ALWAYS
-				),
-			"materialNamesList"
-			)
+		*this,
+		zero3
+	),
+	normalStressWall_
+	(
+		objectFile
+		(
+			"normalStressWall",
+			"",
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_NEVER
 		),
-	propertyId_(
-		control.geometry().emplaceObject<int8TriSurfaceField_D>(
-			objectFile(
-				"propertyId",
-				"",
-				objectFile::READ_NEVER,
-				objectFile::WRITE_NEVER),
-			surface(),
-			0 ) ),
-	contactForceWall_(
-		control.geometry().emplaceObject<realx3TriSurfaceField_D>(
-			objectFile(
-				"contactForceWall",
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			surface(),
-			zero3) ),
-	stressWall_(
-		control.geometry().emplaceObject<realx3TriSurfaceField_D>(
-			objectFile(
-				"stressWall",
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			surface(),
-			zero3) )
+		*this,
+		zero3
+	),
+	shearStressWall_
+	(
+		objectFile
+		(
+			"shearStressWall",
+			"",
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_NEVER
+		),
+		*this,
+		zero3
+	)
 {
 
-	if(!findPropertyId())
+	readWholeObject_ = false;
+	if( !IOobject::readObject() )
+    {
+        fatalErrorInFunction<<
+        "Error in reading from file "<<IOobject::path()<<endl;
+        fatalExit;
+    }
+	readWholeObject_ = true;
+
+	if( this->numSurfaces() != motionComponentName_.size() )
+	{
+		fatalErrorInFunction<<
+		"Number of surfaces is not equal to number of motion component names"<<endl;
+		fatalExit;
+	}
+	
+	if(!createPropertyId())
 	{
 		fatalExit;
 	}
 }
 
 pFlow::geometry::geometry
+(
+	systemControl &control, 
+	const property &prop, 
+	multiTriSurface &surf, 
+	const wordVector &motionCompName, 
+	const wordVector &materialName,
+	const dictionary& motionDict
+)
+:
+	multiTriSurface
+	(
+		objectFile
+		(
+			triSurfaceFile__,
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_ALWAYS
+		),
+		&control.geometry(),
+		surf
+	),
+	demComponent
+	(
+		"geometry",
+		control
+	),
+	wallProperty_
+	(
+		prop
+	),
+	propertyId_
+	(
+		objectFile
+		(
+			"propertyId",
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_NEVER
+		),
+		*this,
+		0u
+	),
+	contactForceWall_
+	(
+		objectFile
+		(
+			"contactForcWall",
+			"",
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_NEVER
+		),
+		*this,
+		zero3
+	),
+	normalStressWall_
+	(
+		objectFile
+		(
+			"normalStressWall",
+			"",
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_NEVER
+		),
+		*this,
+		zero3
+	),
+	shearStressWall_
+	(
+		objectFile
+		(
+			"shearStressWall",
+			"",
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_NEVER
+		),
+		*this,
+		zero3
+	)
+	
+{
+	motionComponentName_.assign(motionCompName);
+	materialName_.assign(materialName);
+
+
+	if( this->numSurfaces() != motionComponentName_.size() )
+	{
+		fatalErrorInFunction<<
+		"Number of surfaces ("<< this->numSurfaces() <<
+        ") is not equal to number of motion component names("<<
+        motionComponentName_.size()<<")"<<endl;
+		fatalExit;
+	}
+
+	if(!createPropertyId())
+	{
+		fatalExit;
+	}
+}
+
+bool pFlow::geometry::beforeIteration()
+{
+	zeroForce();
+    return true;
+}
+
+bool pFlow::geometry::iterate()
+{
+    return true;
+}
+
+bool pFlow::geometry::afterIteration()
+{
+	auto numTris = size();
+	auto& normalsD = normals().deviceViewAll();
+	auto& areaD = area().deviceViewAll();
+	auto& cForceD = contactForceWall_.deviceViewAll();
+	auto& sStressD = shearStressWall_.deviceViewAll();
+	auto& nStressD = normalStressWall_.deviceViewAll();
+
+	Kokkos::parallel_for(
+		"pFlow::geometry::afterIteration",
+		deviceRPolicyStatic(0, numTris),
+		LAMBDA_HD(uint32 i)
+		{	
+			realx3 n = normalsD[i];
+			real A = max(areaD[i],smallValue);
+			realx3 nF = dot(cForceD[i],n)*n;
+			realx3 sF = cForceD[i] - nF;
+			nStressD[i] = nF/A;
+			sStressD[i] = sF/A;
+		});
+	Kokkos::fence();
+    
+	return true;
+}
+
+bool pFlow::geometry::read(iIstream &is, const IOPattern &iop)
+{
+	motionComponentName_.read(is, iop);
+	
+	materialName_.read(is, iop);
+	
+	if( readWholeObject_ )
+	{
+		return multiTriSurface::read(is, iop);
+	}
+
+    return true;
+}
+
+bool pFlow::geometry::write(iOstream &os, const IOPattern &iop) const
+{
+    
+	if( !motionComponentName_.write(os, iop) )
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+	if( !materialName_.write(os,iop))
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+	
+	return multiTriSurface::write(os,iop);
+}
+
+
+pFlow::uniquePtr<pFlow::geometry> 
+	pFlow::geometry::create
+(
+	systemControl& control,
+	const property& prop
+)
+{
+	
+	//
+	fileDictionary dict( motionModelFile__, control.time().geometry().path());
+
+	word model = dict.getVal<word>("motionModel");
+	
+	auto geomModel = angleBracketsNames("geometry", model);
+
+	REPORT(1)<< "Selecting geometry model . . ."<<END_REPORT;
+	if( systemControlvCtorSelector_.search(geomModel) )
+	{
+		auto objPtr = systemControlvCtorSelector_[geomModel] (control, prop);
+		REPORT(2)<<"Model "<< Green_Text(geomModel)<<" is created.\n"<<END_REPORT;
+		return objPtr;
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< Yellow_Text(geomModel) << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			systemControlvCtorSelector_
+		);
+		fatalExit;
+	}
+
+	return nullptr;
+}
+
+pFlow::uniquePtr<pFlow::geometry> 
+	pFlow::geometry::create
+	(
+		systemControl& control, 
+		const property& prop,
+		multiTriSurface& surf,
+		const wordVector& motionCompName,
+		const wordVector& materialName,
+		const dictionary& motionDic
+	)
+{
+
+	word model = motionDic.getVal<word>("motionModel");
+
+	auto geomModel = angleBracketsNames("geometry", model);
+
+	REPORT(1)<< "Selecting geometry model . . ."<<END_REPORT;
+
+	if( dictionaryvCtorSelector_.search(geomModel) )
+	{
+		auto objPtr = dictionaryvCtorSelector_[geomModel] 
+			(
+				control,
+				prop, 
+				surf,
+				motionCompName,
+				materialName,
+				motionDic
+			);
+		REPORT(2)<<"Model "<< Green_Text(geomModel)<<" is created.\n"<<END_REPORT;
+		return objPtr;
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< Yellow_Text(geomModel) << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			dictionaryvCtorSelector_
+		);
+		fatalExit;
+	}
+	return nullptr;
+}
+
+
+
+
+/*pFlow::geometry::geometry
 (
 	systemControl& control,
 	const property& prop,
@@ -224,52 +512,7 @@ pFlow::geometry::geometry
 {}
 
 
-pFlow::uniquePtr<pFlow::geometry> 
-	pFlow::geometry::create
-(
-	systemControl& control,
-	const property& prop
-)
-{
-	//motionModelFile__
-	auto motionDictPtr = IOobject::make<dictionary>
-	(
-		objectFile
-		(
-			motionModelFile__,
-			control.geometry().path(),
-			objectFile::READ_ALWAYS,
-			objectFile::WRITE_NEVER
-		),
-		motionModelFile__,
-		true
-	);
 
-	word model = motionDictPtr().getObject<dictionary>().getVal<word>("motionModel");
-	
-	auto geomModel = angleBracketsNames("geometry", model);
-
-	REPORT(1)<< "Selecting geometry model . . ."<<endREPORT;
-	if( systemControlvCtorSelector_.search(geomModel) )
-	{
-		auto objPtr = systemControlvCtorSelector_[geomModel] (control, prop);
-		REPORT(2)<<"Model "<< greenText(geomModel)<<" is created.\n"<<endREPORT;
-		return objPtr;
-	}
-	else
-	{
-		printKeys
-		( 
-			fatalError << "Ctor Selector "<< yellowText(geomModel) << " dose not exist. \n"
-			<<"Avaiable ones are: \n\n"
-			,
-			systemControlvCtorSelector_
-		);
-		fatalExit;
-	}
-
-	return nullptr;
-}
 
 bool pFlow::geometry::beforeIteration()
 { 
@@ -296,48 +539,4 @@ bool pFlow::geometry::afterIteration()
 	Kokkos::fence();
 	return true;
 	
-}
-
-pFlow::uniquePtr<pFlow::geometry> 
-	pFlow::geometry::create(
-			systemControl& control,
-			const property& prop,
-			const dictionary& dict,
-			const multiTriSurface& triSurface,
-			const wordVector& motionCompName,
-			const wordVector& propName)
-{
-
-	word model = dict.getVal<word>("motionModel");
-
-	auto geomModel = angleBracketsNames("geometry", model);
-
-	REPORT(1)<< "Selecting geometry model . . ."<<endREPORT;
-
-	if( dictionaryvCtorSelector_.search(geomModel) )
-	{
-		auto objPtr = dictionaryvCtorSelector_[geomModel] 
-			(
-				control,
-				prop, 
-				dict,
-				triSurface,
-				motionCompName,
-				propName
-			);
-		REPORT(2)<<"Model "<< greenText(geomModel)<<" is created.\n"<<endREPORT;
-		return objPtr;
-	}
-	else
-	{
-		printKeys
-		( 
-			fatalError << "Ctor Selector "<< yellowText(geomModel) << " dose not exist. \n"
-			<<"Avaiable ones are: \n\n"
-			,
-			dictionaryvCtorSelector_
-		);
-		fatalExit;
-	}
-	return nullptr;
-}
+}*/

src/Geometry/geometry/geometry.hpp

@@ -23,13 +23,14 @@ Licence:
 
 
 #include "virtualConstructor.hpp"
-#include "demGeometry.hpp"
+#include "demComponent.hpp"
 #include "property.hpp"
-#include "Fields.hpp"
-#include "Vectors.hpp"
 #include "multiTriSurface.hpp"
 #include "triSurfaceFields.hpp"
-#include "dictionary.hpp"
+//#include "Fields.hpp"
+//#include "Vectors.hpp"
+
+
 
 namespace pFlow
 {
@@ -42,48 +43,49 @@ namespace pFlow
  */
 class geometry
 :
-	public demGeometry
+	public multiTriSurface,
+	public demComponent
 {
-protected:
+private:
 
 	// - Protected members 
 
 		/// Const reference to physical property of materials 
 		const property& 	wallProperty_;
 
-		/// Repository to store geometry data at each simulation moment
-		repository&  		geometryRepository_;
-
-		/// All triangles in the set of wall surfaces 
-		multiTriSurface&  	triSurface_;
-
 		/// The name of motion component of each wall surface 
-		wordField& 			motionComponentName_;
+		wordField_H 		motionComponentName_{
+								"motionComponentName",	
+								"motionComponentName"};
 
 		/// Material name of each wall surface  
-		wordField&   		materialName_;
+		wordField_H   		materialName_{
+								"materialName",
+								"materialName"};
 		
 		/// Property id of each triangle in the set of wall surfaces
-		int8TriSurfaceField_D& 		propertyId_;
+		uint32TriSurfaceField_D 		propertyId_;
 
 		/// Contact force on each triangle in the set of wall surfaces
-		realx3TriSurfaceField_D& 	contactForceWall_;
+		realx3TriSurfaceField_D 		contactForceWall_;
 
-		/// Stress on ech triangle in the set of wall surfaces 
-		realx3TriSurfaceField_D&    stressWall_;
+		/// Stress on each triangle in the set of wall surfaces 
+		realx3TriSurfaceField_D 		normalStressWall_;
+
+		/// Stress on each triangle in the set of wall surfaces 
+		realx3TriSurfaceField_D 		shearStressWall_;
+
+
+		bool readWholeObject_ = true; 			
 
 	// - Protected member functions
 	
 		/// Find property id of each triangle based on the supplied material name
 		/// and the surface wall that the triangle belongs to. 
-		bool findPropertyId();
+		bool createPropertyId();
 
 		/// Initialize contact force to zero 
-		void zeroForce()
-		{
-			contactForceWall_.fill(zero3);
-		}
-
+		void zeroForce();
 		
 public:
 
@@ -95,8 +97,15 @@ public:
 		/// Construct from controlSystem and property, for reading from file
 		geometry(systemControl& control, const property& prop);
 
-		/// Construct from components
 		geometry(systemControl& control, 
+			const property& prop,
+			multiTriSurface& surf,
+			const wordVector& motionCompName,
+			const wordVector& materialName,
+			const dictionary& motionDict);
+
+		/// Construct from components
+		/*geometry(systemControl& control,
 				 const property& prop,
 				 const multiTriSurface& triSurface,
 				 const wordVector& motionCompName,
@@ -110,10 +119,10 @@ public:
 				 const dictionary& dict,
 				 const multiTriSurface& triSurface,
 				 const wordVector& motionCompName,
-				 const wordVector& propName);
+				 const wordVector& propName);*/
 
 		/// Destructor
-		virtual ~geometry() = default;
+		~geometry()override = default;
 
 		/// Virtual constructor 
 		create_vCtor
@@ -132,156 +141,96 @@ public:
 		(
 			geometry,
 			dictionary,
-			(systemControl& control,
-			 const property& prop,
-			 const dictionary& dict,
-			 const multiTriSurface& triSurface,
-			 const wordVector& motionCompName,
-			 const wordVector& propName),
-			(control, prop, dict, triSurface, motionCompName, propName)
+			(
+				systemControl& control, 
+				const property& prop,
+				multiTriSurface& surf,
+				const wordVector& motionCompName,
+				const wordVector& materialName,
+				const dictionary& motionDic),
+			(control, prop, surf, motionCompName, materialName, motionDic)
 		);
 
 	//- Methods 	
-
-		/// Size of tri-surface 
 		inline 
-		auto size()const
+		const auto& motionComponentName()const
 		{
-			return triSurface_.size();
-		}
-
-		/// Number of points in the set of surface walls 
-		inline
-		auto numPoints()const
-		{
-			return triSurface_.numPoints();
-		}
-
-		/// Number of triangles in the set of surface walls 
-		inline
-		auto numTriangles()const
-		{
-			return size();
-		}
-
-		/// Access to the points
-		inline
-		const auto& points()const
-		{
-			return triSurface_.points();
-		}
-
-		/// Access to the vertices 
-		inline
-		const auto& vertices()const
-		{
-			return triSurface_.vertices();
-		}
-
-		/// Obtain an object for accessing triangles 
-		inline auto getTriangleAccessor()const
-		{
-			return triSurface_.getTriangleAccessor();
-		}
-
-		/// Surface 
-		inline auto& surface()
-		{
-			return triSurface_;
-		}
-
-		/// Surface 
-		inline const auto& surface()const
-		{
-			return triSurface_;
+			return motionComponentName_;
 		}
 		
 		/// Access to contact force
 		inline
-		realx3TriSurfaceField_D& contactForceWall()
+		auto& contactForceWall()
 		{
 			return contactForceWall_;
 		}
 
 		/// Access to contact force
 		inline
-		const realx3TriSurfaceField_D& contactForceWall() const
+		const auto& contactForceWall() const
 		{
 			return contactForceWall_;
 		}
 
+		/// Property ide of triangles 
+		inline
+		const auto& propertyId()const
+		{
+			return propertyId_;
+		}
+
 		/// Access to property 
 		inline const auto& wallProperty()const
 		{
 			return wallProperty_;
 		}
 
-		/// Owner repository
-		inline
-		const repository& owner()const
-		{
-			return geometryRepository_;
-		}
-
-		/// Owner repository
-		inline
-		repository& owner()
-		{
-			return geometryRepository_;
-		}
-
-		/// Path to the repository folder 
-		inline auto path()
-		{
-			return owner().path();
-		}
-
 		/// The name of motion model 
 		virtual 
 		word motionModelTypeName()const = 0;
 
 		/// Motion model index of triangles 
 		virtual
-		const int8Vector_HD& 			triMotionIndex() const =0;
+		const uint32Field_D& triMotionIndex() const =0;
 
 		/// Motion model index of points 
 		virtual 
-		const int8Vector_HD& 			pointMotionIndex()const = 0;
+		const uint32Field_D& pointMotionIndex()const = 0;
 
-		/// Property ide of triangles 
-		const int8TriSurfaceField_D& 	propertyId() const
-		{
-			return propertyId_;
-		}
 		
-		/// Operations before each iteration 
 		bool beforeIteration() override;
 		
-		/// Operations after each iteration
+		/// This is called in time loop. Perform the main calculations 
+		/// when the component should evolve along time.
+		bool iterate() override;
+		
+		/// This is called in time loop, after iterate.
 		bool afterIteration() override;
 		
 	//- IO
 
+		bool read(iIstream& is, const IOPattern& iop) override;
+
 		/// write 
-		bool write()const
-		{
-			return owner().write();
-		}
+		bool write( iOstream& os, const IOPattern& iop )const override;
+		
 
 
 	//- Static members 
 
 		static
-		uniquePtr<geometry> create(systemControl& control, const property& prop);
+		uniquePtr<geometry> create(
+			systemControl& control, 
+			const property& prop);
 
 		static
 		uniquePtr<geometry> create(
-				systemControl& control,
-				const property& prop,
-				const dictionary& dict,
-				const multiTriSurface& triSurface,
-				const wordVector& motionCompName,
-				const wordVector& propName);
+			systemControl& control, 
+			const property& prop,
+			multiTriSurface& surf,
+			const wordVector& motionCompName,
+			const wordVector& materialName,
+			const dictionary& motionDic);
 
 };
 

src/Geometry/geometryMotion/geometryMotion.cpp

@@ -21,41 +21,116 @@ Licence:
 template<typename MotionModel>
 bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
 {
-	motionIndex_.clear();
-	triMotionIndex_.reserve( this->surface().capacity() );
-	triMotionIndex_.clear();
 	
-	ForAll( surfI, motionComponentName_)
+	if(motionComponentName().size() != numSurfaces() )
 	{
-		auto mName = motionComponentName_[surfI];
-		auto mInd = motionModel_.nameToIndex(mName);
-		motionIndex_.push_back(mInd);
-		// fill motionIndex for triangles of the surface
-		int32 surfSize = this->surface().surfNumTriangles(surfI);
-		for(int32 i=0; i<surfSize; i++)
+		fatalErrorInFunction<<
+		"size of motion component names in the triSurface is not"<<
+		" equal to size of number of sub-surfaces"<<endl;
+		 return false;
+	}
+
+	uint32Vector surfMotionIndex("surfMotionIndex");
+	uint32Vector triMotionIndex("triMotionIndex");
+	uint32Vector pointMotionIndex("pointMotionIndex");
+
+	ForAll( surfI, motionComponentName())
+	{
+		auto mName = motionComponentName()[surfI];
+		uint32 mInd=0;
+		
+		if( !motionModel_.nameToIndex(mName, mInd) )
 		{
-			triMotionIndex_.push_back(mInd);
+			fatalErrorInFunction<<
+			mName<< " does not exist in the list of motion names -> "<<
+			motionModel_.componentNames();
+            return false;
+		}
+		surfMotionIndex.push_back(mInd);
+		
+		auto surfRange = subSurfaceRange(surfI);
+		
+		for(uint32 i=0; i<surfRange.numElements(); i++)
+		{
+			triMotionIndex.push_back(mInd);
+		}
+
+		auto pointRange = subSurfacePointRange(surfI);
+		for(uint32 n=0; n<pointRange.numElements(); n++)
+		{
+			pointMotionIndex.push_back(mInd);
 		}
 	}
-	motionIndex_.syncViews();
-	triMotionIndex_.syncViews();
 	
-	pointMotionIndex_.reserve(triSurface_.numPoints());
-	pointMotionIndex_.clear();
-
-	ForAll(surfI, motionIndex_)
-	{
-		auto nP = triSurface_.surfNumPoints(surfI);
-		for(int32 i=0; i<nP; i++)
-		{
-			pointMotionIndex_.push_back(motionIndex_[surfI]);
-		}
-	}
-	pointMotionIndex_.syncViews();
+	surfMotionIndex_.assign(surfMotionIndex);
+	triMotionIndex_.assign(triMotionIndex);
+	pointMotionIndex_.assign(pointMotionIndex);
 
 	return true;
 }
 
+namespace pFlow::GMotion
+{
+    template<typename ModelInterface>
+    void moveGeometry(
+        real dt,
+        uint32 nPoints,
+        const ModelInterface& mModel,
+        const deviceViewType1D<uint32>& pointMIndexD,
+        const deviceViewType1D<realx3>& pointsD
+    )
+    {
+        Kokkos::parallel_for(
+         "geometryMotion<MotionModel>::movePoints",
+         deviceRPolicyStatic(0, nPoints),
+         LAMBDA_HD(uint32 i){
+             auto newPos = mModel.transferPoint(pointMIndexD[i], pointsD[i], dt);
+             pointsD[i] = newPos;
+         });
+
+        Kokkos::fence();
+    }
+}
+
+template<typename MotionModel>
+ bool pFlow::geometryMotion<MotionModel>::moveGeometry()
+ {
+
+    uint32 iter = this->currentIter();
+    real t = this->currentTime();
+    real dt = this->dt();
+        
+    auto mModel = motionModel_.getModelInterface(iter, t, dt);
+    
+    auto& pointMIndexD= pointMotionIndex_.deviceViewAll();
+    auto& pointsD = points().deviceViewAll();
+    
+    pFlow::GMotion::moveGeometry(
+        dt, 
+        numPoints(), 
+        motionModel_.getModelInterface(iter, t, dt),
+        pointMotionIndex_.deviceViewAll(),
+        points().deviceViewAll()
+    );
+
+    /*Kokkos::parallel_for(
+         "geometryMotion<MotionModel>::movePoints",
+         deviceRPolicyStatic(0, numPoints()),
+         LAMBDA_HD(uint32 i){
+             auto newPos = mModel.transferPoint(pointMIndexD[i], pointsD[i], dt);
+             pointsD[i] = newPos;
+         });
+
+    Kokkos::fence();*/
+
+    // move the motion components
+    motionModel_.move(iter, t,dt);
+
+    // end of calculations
+     
+
+     return true;
+ }
 
 template<typename MotionModel>
 pFlow::geometryMotion<MotionModel>::geometryMotion
@@ -65,144 +140,75 @@ pFlow::geometryMotion<MotionModel>::geometryMotion
 )
 :
 	geometry(control, prop),
-	motionModel_(
-		this->owner().template emplaceObject<MotionModel>(
-			objectFile(
-				motionModelFile__,
-				"",
-				objectFile::READ_ALWAYS,
-				objectFile::WRITE_ALWAYS
-				)
-			)
+	motionModel_
+	(
+		objectFile
+		(
+			motionModelFile__,
+			"",
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_ALWAYS
 		),
+		owner()
+	),
 	moveGeomTimer_("move geometry", &this->timers())
 {
-	findMotionIndex();
-}
-
-template<typename MotionModel>
-pFlow::geometryMotion<MotionModel>::geometryMotion
-(
-	systemControl& control,
-	const property& prop,
-	const multiTriSurface& triSurface,
-	const wordVector& motionCompName,
-	const wordVector& propName,
-	const MotionModel& motionModel
-)
-:
-	geometry(
-		control,
-		prop,
-		triSurface,
-		motionCompName,
-		propName
-		),
-	motionModel_(
-		this->owner().template emplaceObject<MotionModel>(
-			objectFile(
-				motionModelFile__,
-				"",
-				objectFile::READ_NEVER,
-				objectFile::WRITE_ALWAYS
-				),
-			motionModel
-			)
-		),
-	moveGeomTimer_("move geometry", &this->timers())
-{
-	findMotionIndex();
-}
-
-template<typename MotionModel>
-pFlow::geometryMotion<MotionModel>::geometryMotion
-(
-	systemControl& control,
-	const property& prop,
-	const dictionary& dict,
-	const multiTriSurface& triSurface,
-	const wordVector& motionCompName,
-	const wordVector& propName
-)
-:
-	geometry(
-		control,
-		prop,
-		dict,
-		triSurface,
-		motionCompName,
-		propName
-		),
-	motionModel_(
-		this->owner().template emplaceObject<MotionModel>(
-			objectFile(
-				motionModelFile__,
-				"",
-				objectFile::READ_NEVER,
-				objectFile::WRITE_ALWAYS
-				),
-			dict
-			)
-		),
-	moveGeomTimer_("move geometry", &this->timers())
-{
-	findMotionIndex();
-}
-
-template<typename MotionModel>
-bool pFlow::geometryMotion<MotionModel>::beforeIteration() 
-{ 
-	geometry::beforeIteration();
-	return true;
-}
-
-template<typename MotionModel>
-bool pFlow::geometryMotion<MotionModel>::iterate()
-{
-	if( motionModel_.isMoving() )
+	if(!findMotionIndex())
 	{
-		moveGeomTimer_.start();
-		moveGeometry();
-		moveGeomTimer_.end();
+		fatalExit;
 	}
-	return true;
 }
 
-template<typename MotionModel>
-bool pFlow::geometryMotion<MotionModel>::afterIteration() 
+template <typename MotionModelType>
+pFlow::geometryMotion<MotionModelType>::geometryMotion
+(
+	systemControl &control, 
+	const property &prop, 
+	multiTriSurface &surf, 
+	const wordVector &motionCompName, 
+	const wordVector &materialName, 
+	const dictionary &motionDict
+)
+:
+	geometry
+	(
+		control, 
+		prop, 
+		surf, 
+		motionCompName, 
+		materialName,
+        motionDict
+	),
+	motionModel_
+	(
+		objectFile
+		(
+			motionModelFile__,
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_ALWAYS
+		),
+		motionDict,
+		owner()
+	),
+	moveGeomTimer_("move geometry", &this->timers())
 {
-	geometry::afterIteration();
-	return true;
+	if(!findMotionIndex())
+	{
+		fatalExit;
+	}
 } 
 
-template<typename MotionModel>
-bool pFlow::geometryMotion<MotionModel>::moveGeometry()
-{
 
-	real dt = this->dt();
-	real t = this->currentTime();
-
-	auto pointMIndex= pointMotionIndex_.deviceVector();
-	auto mModel = motionModel_.getModel(t);
-	realx3* points = triSurface_.pointsData_D();
-	auto numPoints = triSurface_.numPoints();
-
-
-	Kokkos::parallel_for(
-		"geometryMotion<MotionModel>::movePoints",
-		numPoints,
-		LAMBDA_HD(int32 i){
-			auto newPos = mModel.transferPoint(pointMIndex[i], points[i], dt);
-			points[i] = newPos;
-		});
-
-	Kokkos::fence();
-
-	// move the motion components 
-	motionModel_.move(t,dt);
-
-	// end of calculations 
-	moveGeomTimer_.end();
-
-	return true;
-}
+ template<typename MotionModel>
+ bool pFlow::geometryMotion<MotionModel>::iterate()
+ {
+     if( motionModel_.isMoving() )
+     {
+        moveGeomTimer_.start();
+          moveGeometry();
+          this->calculateNormals();
+        moveGeomTimer_.end();
+     }
+     return true;
+ }

src/Geometry/geometryMotion/geometryMotion.hpp

@@ -20,9 +20,8 @@ Licence:
 #ifndef __geometryMotion_hpp__
 #define __geometryMotion_hpp__
 
-
+#include "vocabs.hpp"
 #include "geometry.hpp"
-#include "VectorDuals.hpp"
 
 namespace pFlow
 {
@@ -39,21 +38,23 @@ class	geometryMotion
 {
 public:
 	
-	using MotionModel = MotionModelType;
+	using MotionModel 	 = MotionModelType;
 
-protected: 
+	using ModelComponent = typename MotionModelType::ModelComponent;
+
+private: 
 	
 	/// Ref to motion model 
-	MotionModel& 		motionModel_;
+	MotionModelType 	motionModel_;
 	
 	/// motion indext mapped on each surface 
-	int32Vector_HD  	motionIndex_;
+	uint32Field_D  		surfMotionIndex_{"triMotionIndex"};
 
 	/// motion index mapped on each triangle 
-	int8Vector_HD 		triMotionIndex_; 
+	uint32Field_D 		triMotionIndex_ {"surfMotionIndex"}; 
 
 	/// motion index mapped on each point 
-	int8Vector_HD		pointMotionIndex_;
+	uint32Field_D		pointMotionIndex_{"pointMotionIndex"};
 
 	/// timer for moveGeometry
 	Timer 				moveGeomTimer_;			
@@ -61,32 +62,25 @@ protected:
 	/// determine the motion index of each triangle 
 	bool findMotionIndex();
 	
+	/// Move geometry 
+	bool moveGeometry();
+
 public:
 
 	/// Type info
-	TypeInfoTemplate("geometry", MotionModel);
+	TypeInfoTemplate11("geometry", ModelComponent);
 
 	// - Constructors
 		
 		geometryMotion(systemControl& control, const property& prop);
 
-		
 		geometryMotion(
 			systemControl& control, 
 			const property& prop,
-			const multiTriSurface& triSurface,
+			multiTriSurface& surf,
 			const wordVector& motionCompName,
-			const wordVector& propName,
-			const MotionModel& motionModel);
-
-		/// Construct from components and dictionary that contains 
-		/// motionModel
-		geometryMotion(systemControl& control,
-				 const property& prop,
-				 const dictionary& dict,
-				 const multiTriSurface& triSurface,
-				 const wordVector& motionCompName,
-				 const wordVector& propName);
+			const wordVector& materialName,
+			const dictionary& motionDict);
 
 		/// Add virtual constructor 
 		add_vCtor
@@ -107,9 +101,9 @@ public:
 	// - Methods
 
 		/// Obtain motion model at time t 
-		auto getModel(real t)const
+		auto getModel(uint32 iter, real t, real dt)const
 		{
-			return motionModel_.getModel(t);
+			return motionModel_.getModelInterface(iter, t, dt);
 		}
 
 		/// TypeName / TypeInfo of motion model 
@@ -119,28 +113,21 @@ public:
 		}
 
 		/// Access to motion model index of triangles 
-		const int8Vector_HD& triMotionIndex()const override
+		const uint32Field_D& triMotionIndex()const override
 		{
 			return triMotionIndex_;
 		}
 
 		/// Access to motion model index of points 
-		const int8Vector_HD& pointMotionIndex()const override
+		const uint32Field_D& pointMotionIndex()const override
 		{
 			return pointMotionIndex_;
 		}
 
-		/// Operations before each iteration 
-		bool beforeIteration() override;
-
 		/// Iterate geometry one time step  
 		bool iterate() override ;
 
-		/// Operations after each iteration 
-		bool afterIteration() override;
 		
-		/// Move geometry 
-		bool moveGeometry();
 
 };
 
@@ -148,9 +135,6 @@ public:
 
 #include "geometryMotion.cpp"
 
-#ifndef BUILD_SHARED_LIBS
-	#include "geometryMotionsInstantiate.cpp"
-#endif
 
 
 #endif  //__geometryMotion_hpp__

src/Geometry/geometryMotion/geometryMotions.cpp

@@ -17,9 +17,15 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
 #include "geometryMotions.hpp"
 
-#ifdef BUILD_SHARED_LIBS
-#include "geometryMotionsInstantiate.cpp"
-#endif
+
+template class pFlow::geometryMotion<pFlow::vibratingMotion>;
+
+template class pFlow::geometryMotion<pFlow::rotatingAxisMotion>;
+
+template class pFlow::geometryMotion<pFlow::stationaryWall>;
+
+template class pFlow::geometryMotion<pFlow::conveyorBeltMotion>;
+
+//template class pFlow::geometryMotion<pFlow::multiRotatingAxisMotion>;

src/Geometry/geometryMotion/geometryMotions.hpp

@@ -22,22 +22,26 @@ Licence:
 #define __geometryMotions_hpp__
 
 #include "geometryMotion.hpp"
-#include "fixedWall.hpp"
+#include "stationaryWall.hpp"
 #include "rotatingAxisMotion.hpp"
-#include "multiRotatingAxisMotion.hpp"
+#include "conveyorBeltMotion.hpp"
+//#include "multiRotatingAxisMotion.hpp"
 #include "vibratingMotion.hpp"
 
 
 namespace pFlow
 {
 
-typedef geometryMotion<vibratingMotion> vibratingMotionGeometry;
+using vibratingMotionGeometry = geometryMotion<vibratingMotion>;
 
-typedef geometryMotion<rotatingAxisMotion> rotationAxisMotionGeometry;
+using rotationAxisMotionGeometry = geometryMotion<rotatingAxisMotion>;
 
-typedef geometryMotion<multiRotatingAxisMotion> multiRotationAxisMotionGeometry;
+using stationaryGeometry = geometryMotion<stationaryWall>;
+
+using conveyorBeltMotionGeometry = geometryMotion<conveyorBeltMotion>;
+
+//typedef geometryMotion<multiRotatingAxisMotion> multiRotationAxisMotionGeometry;
 
-typedef geometryMotion<fixedWall> fixedGeometry;
 
 
 

src/Integration/AdamsBashforth2/AB2Kernels.hpp

@@ -0,0 +1,63 @@
+
+
+#ifndef __AB2Kernels_hpp__
+#define __AB2Kernels_hpp__
+
+#include "KokkosTypes.hpp"
+#include "types.hpp"
+
+namespace pFlow::AB2Kernels
+{
+inline
+bool intAllActive(
+    const word& name,
+	real dt, 
+    rangeU32 activeRng,
+	const deviceViewType1D<realx3>& y, 
+	const deviceViewType1D<realx3>& dy,
+	const deviceViewType1D<realx3>& dy1
+)
+{
+	Kokkos::parallel_for(
+		name,
+		deviceRPolicyStatic (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			y[i] +=  dt*(static_cast<real>(1.5) * dy[i] - static_cast<real>(0.5) * dy1[i]);
+			dy1[i] = dy[i];
+		});
+	Kokkos::fence();
+
+	return true;	
+}
+
+inline
+bool intScattered
+(
+	const word& name,
+	real dt, 
+    const pFlagTypeDevice& activePoints,
+	const deviceViewType1D<realx3>& y, 
+	const deviceViewType1D<realx3>& dy,
+	const deviceViewType1D<realx3>& dy1
+)
+{
+
+	Kokkos::parallel_for(
+		name,
+		deviceRPolicyStatic (activePoints.activeRange().start(), activePoints.activeRange().end()),
+		LAMBDA_HD(uint32 i){
+			if( activePoints(i))
+			{
+				y[i] +=  dt*(static_cast<real>(1.5) * dy[i] - static_cast<real>(0.5) * dy1[i]);
+				dy1[i] = dy[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+}
+
+#endif

src/Integration/AdamsBashforth2/AdamsBashforth2.cpp

@@ -19,61 +19,171 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "AdamsBashforth2.hpp"
+#include "pointStructure.hpp"
+#include "Time.hpp"
+#include "vocabs.hpp"
 
-//const real AB2_coef[] = { 3.0 / 2.0, 1.0 / 2.0};
+namespace pFlow
+{
+
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
+
+bool intAllActive(
+	real dt, 
+	realx3Field_D& y, 
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1)
+{
+
+	auto d_dy = dy.deviceView();
+	auto d_y  = y.deviceView();
+	auto d_dy1= dy1.deviceView();
+	auto activeRng = dy1.activeRange();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth2::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			d_y[i] +=  dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
+			d_dy1[i] = d_dy[i];
+		});
+	Kokkos::fence();
+
+	return true;	
+}
+
+bool intScattered
+(
+	real dt, 
+	realx3Field_D& y,
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth2::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] +=  dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+}
 
 pFlow::AdamsBashforth2::AdamsBashforth2
 (
 	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
+	pointStructure& pStruct,
+	const word& method,
+	const realx3Field_D& initialValField
 )
 :
-	integration(baseName, owner, pStruct, method),
-	dy1_(
-		owner.emplaceObject<realx3PointField_D>(
-			objectFile(
-				groupNames(baseName,"dy1"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			zero3))
-{
+	integration(baseName, pStruct, method, initialValField),
+	realx3PointField_D
+	(
+		objectFile
+		(
+			groupNames(baseName,"dy1"),
+			pStruct.time().integrationFolder(),
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_ALWAYS
+		),
+		pStruct,
+		zero3,
+		zero3
+	),
+	boundaryList_(pStruct, method, *this)
+{}
 
+void pFlow::AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested()
+{
+	realx3PointField_D::updateBoundariesSlaveToMasterIfRequested();
+}
+
+bool pFlow::AdamsBashforth2::predict(
+    real UNUSED(dt),
+    realx3PointField_D &UNUSED(y),
+    realx3PointField_D &UNUSED(dy))
+{
+	return true;
 }
 
 bool pFlow::AdamsBashforth2::predict
 (
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
+	real dt, 
+	realx3Field_D &y, 
+	realx3PointField_D &dy
 )
 {
-
-	return true;
+    return true;
 }
 
 bool pFlow::AdamsBashforth2::correct
 (
 	real dt,
-	realx3Vector_D& y,
-	realx3Vector_D& dy
+	realx3PointField_D& y,
+	realx3PointField_D& dy
 )
 {
-	if(this->pStruct().allActive())
+	auto& dy1l = dy1();
+	bool success = false;
+	if(dy1l.isAllActive())
 	{
-		return intAll(dt, y, dy, this->pStruct().activeRange());
+		success = intAllActive(dt, y.field(), dy, dy1l);
 	}
 	else
 	{
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
+		success = intScattered(dt, y.field(), dy, dy1l);
 	}
 
-	return true;
+	success = success && boundaryList_.correct(dt, y, dy);
+
+	return success;
+
 }
 
+bool pFlow::AdamsBashforth2::correctPStruct(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel)
+{
+	auto& dy1l = dy1();
+	bool success = false;
+	if(dy1l.isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1l);
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1l);
+	}
+
+	success = success && boundaryList_.correctPStruct(dt, pStruct, vel);
+
+	return success;
+}
+
+
 bool pFlow::AdamsBashforth2::setInitialVals(
 	const int32IndexContainer& newIndices,
 	const realx3Vector& y)
@@ -81,25 +191,3 @@ bool pFlow::AdamsBashforth2::setInitialVals(
 	return true;
 }
 
-bool pFlow::AdamsBashforth2::intAll(
-	real dt, 
-	realx3Vector_D& y, 
-	realx3Vector_D& dy, 
-	range activeRng)
-{
-
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
-	auto d_dy1= dy1_.deviceVectorAll();
-
-	Kokkos::parallel_for(
-		"AdamsBashforth2::correct",
-		rpIntegration (activeRng.first, activeRng.second),
-		LAMBDA_HD(int32 i){
-			d_y[i] +=  dt*(static_cast<real>(3.0 / 2.0) * d_dy[i] - static_cast<real>(1.0 / 2.0) * d_dy1[i]);
-			d_dy1[i] = d_dy[i];
-		});
-	Kokkos::fence();
-
-	return true;	
-}

src/Integration/AdamsBashforth2/AdamsBashforth2.hpp

@@ -17,13 +17,13 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
 #ifndef __AdamsBashforth2_hpp__
 #define __AdamsBashforth2_hpp__
 
 
 #include "integration.hpp"
 #include "pointFields.hpp"
+#include "boundaryIntegrationList.hpp"
 
 namespace pFlow
 {
@@ -36,41 +36,41 @@ namespace pFlow
  */
 class AdamsBashforth2
 :
-	public integration
+	public integration,
+	public realx3PointField_D
 {
-protected:
+private:
 
-	/// dy at t-dt
-	realx3PointField_D& dy1_;
+	boundaryIntegrationList boundaryList_;
 
-	/// Range policy for integration kernel (alias)
-	using rpIntegration = Kokkos::RangePolicy<
-			DefaultExecutionSpace,
-			Kokkos::Schedule<Kokkos::Static>,
-			Kokkos::IndexType<int32>
-			>;
+	friend class processorAB2BoundaryIntegration;
+
+	const auto& dy1()const
+	{
+		return static_cast<const realx3PointField_D&>(*this);
+	}
+
+	auto& dy1()
+	{
+		return static_cast<realx3PointField_D&>(*this);
+	}
 	
 public:
 
-	/// Type info
-	TypeInfo("AdamsBashforth2");
+	/// Class info
+	ClassInfo("AdamsBashforth2");
 
 	// - Constructors
 
 		/// Construct from components
 		AdamsBashforth2(
 			const word& baseName,
-			repository& owner,
-			const pointStructure& pStruct,
-			const word& method);
-
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth2>(*this);
-		}
+			pointStructure& pStruct,
+			const word& method,
+			const realx3Field_D& initialValField);
 
 		/// Destructor 
-		virtual ~AdamsBashforth2()=default;
+		~AdamsBashforth2()final = default;
 
 		/// Add this to the virtual constructor table 
 		add_vCtor(
@@ -81,70 +81,55 @@ public:
 
 	// - Methods
 
+		void updateBoundariesSlaveToMasterIfRequested()override;
+
+		/// return integration method 
+		word method()const override
+		{
+			return "AdamsBashforth2";
+		}
+
 		bool predict(
 			real UNUSED(dt), 
-			realx3Vector_D& UNUSED(y), 
-			realx3Vector_D& UNUSED(dy)) override;
+			realx3PointField_D& UNUSED(y), 
+			realx3PointField_D& UNUSED(dy)) final;
+		
+		bool predict(
+			real dt, 
+			realx3Field_D& y, 
+			realx3PointField_D& dy) final;
 
 		bool correct(
 			real dt, 
-			realx3Vector_D& y, 
-			realx3Vector_D& dy) override;
+			realx3PointField_D& y, 
+			realx3PointField_D& dy) final;
+
+		bool correctPStruct(
+			real dt, 
+			pointStructure& pStruct, 
+			realx3PointField_D& vel) final;
+			
+		
+		/*bool hearChanges
+		(
+			real t,
+			real dt,
+			uint32 iter,
+			const message& msg, 
+			const anyList& varList
+		) override;*/
 
 		bool setInitialVals(
 			const int32IndexContainer& newIndices,
-			const realx3Vector& y) override;
+			const realx3Vector& y) final;
 
-		bool needSetInitialVals()const override
+		bool needSetInitialVals()const final
 		{
 			return false;
 		}
 
-		/// Integrate on all points in the active range 
-		bool intAll(
-			real dt, 
-			realx3Vector_D& y, 
-			realx3Vector_D& dy, 
-			range activeRng);
-
-		/// Integrate on active points in the active range 
-		template<typename activeFunctor>
-		bool intRange(
-			real dt, 
-			realx3Vector_D& y, 
-			realx3Vector_D& dy, 
-			activeFunctor activeP );
-
 };
 
-template<typename activeFunctor>
-bool pFlow::AdamsBashforth2::intRange(
-	real dt, 
-	realx3Vector_D& y,
-	realx3Vector_D& dy,
-	activeFunctor activeP )
-{
-
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
-	auto d_dy1= dy1_.deviceVectorAll();
-	auto activeRng = activeP.activeRange();
-
-	Kokkos::parallel_for(
-		"AdamsBashforth2::correct",
-		rpIntegration (activeRng.first, activeRng.second),
-		LAMBDA_HD(int32 i){
-			if( activeP(i))
-			{
-				d_y[i] +=  dt*(static_cast<real>(3.0 / 2.0) * d_dy[i] - static_cast<real>(1.0 / 2.0) * d_dy1[i]);
-				d_dy1[i] = d_dy[i];
-			}
-		});
-	Kokkos::fence();
-
-
-	return true;
-}
 
 } // pFlow
 

src/Integration/AdamsBashforth3/AdamsBashforth3.cpp

@@ -89,9 +89,9 @@ bool pFlow::AdamsBashforth3::intAll(
 	realx3Vector_D& dy, 
 	range activeRng)
 {
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
-	auto d_history = history_.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
+	auto d_history = history_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth3::correct",

src/Integration/AdamsBashforth3/AdamsBashforth3.hpp

@@ -160,9 +160,9 @@ bool pFlow::AdamsBashforth3::intRange(
 	realx3Vector_D& dy,
 	activeFunctor activeP )
 {
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
-	auto d_history = history_.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
+	auto d_history = history_.deviceViewAll();
 	auto activeRng = activeP.activeRange();
 
 	Kokkos::parallel_for(

src/Integration/AdamsBashforth4/AdamsBashforth4.cpp

@@ -89,9 +89,9 @@ bool pFlow::AdamsBashforth4::intAll(
 	realx3Vector_D& dy, 
 	range activeRng)
 {
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
-	auto d_history = history_.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
+	auto d_history = history_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth4::correct",

src/Integration/AdamsBashforth4/AdamsBashforth4.hpp

@@ -165,9 +165,9 @@ bool pFlow::AdamsBashforth4::intRange(
 	realx3Vector_D& dy,
 	activeFunctor activeP )
 {
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
-	auto d_history = history_.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
+	auto d_history = history_.deviceViewAll();
 	auto activeRng = activeP.activeRange();
 
 	Kokkos::parallel_for(

src/Integration/AdamsBashforth5/AdamsBashforth5.cpp

@@ -89,9 +89,9 @@ bool pFlow::AdamsBashforth5::intAll(
 	realx3Vector_D& dy, 
 	range activeRng)
 {
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
-	auto d_history = history_.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
+	auto d_history = history_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth5::correct",

src/Integration/AdamsBashforth5/AdamsBashforth5.hpp

@@ -165,9 +165,9 @@ bool pFlow::AdamsBashforth5::intRange(
 	realx3Vector_D& dy,
 	activeFunctor activeP )
 {
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
-	auto d_history = history_.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
+	auto d_history = history_.deviceViewAll();
 	auto activeRng = activeP.activeRange();
 
 	Kokkos::parallel_for(

src/Integration/AdamsMoulton3/AdamsMoulton3.cpp

@@ -124,11 +124,11 @@ bool pFlow::AdamsMoulton3::predictAll(
 	range activeRng)
 {
 
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
-	auto d_y0 = y0_.deviceVectorAll();
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_dy1= dy1_.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
+	auto d_y0 = y0_.deviceViewAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_dy1= dy1_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsMoulton3::predict",
@@ -150,12 +150,12 @@ bool pFlow::AdamsMoulton3::intAll(
 	range activeRng)
 {
 
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
 
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_y0  = y0_.deviceVectorAll();
-	auto d_dy1 = dy1_.deviceVectorAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_y0  = y0_.deviceViewAll();
+	auto d_dy1 = dy1_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsMoulton3::correct",

src/Integration/AdamsMoulton3/AdamsMoulton3.hpp

@@ -144,11 +144,11 @@ bool AdamsMoulton3::predictRange(
 	realx3Vector_D& dy,
 	activeFunctor activeP)
 {
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
-	auto d_y0 = y0_.deviceVectorAll();
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_dy1= dy1_.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
+	auto d_y0 = y0_.deviceViewAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_dy1= dy1_.deviceViewAll();
 	
 	auto activeRng = activeP.activeRange();
 
@@ -182,12 +182,12 @@ bool pFlow::AdamsMoulton3::intRange(
 	activeFunctor activeP)
 {
 
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
 
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_y0  = y0_.deviceVectorAll();
-	auto d_dy1 = dy1_.deviceVectorAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_y0  = y0_.deviceViewAll();
+	auto d_dy1 = dy1_.deviceViewAll();
 
 	auto activeRng = activeP.activeRange();
 

src/Integration/AdamsMoulton4/AdamsMoulton4.cpp

@@ -135,13 +135,13 @@ bool pFlow::AdamsMoulton4::predictAll(
 	range activeRng)
 {
 
-	auto d_dy  = dy.deviceVectorAll();
-	auto d_y   = y.deviceVectorAll();
+	auto d_dy  = dy.deviceViewAll();
+	auto d_y   = y.deviceViewAll();
 
-	auto d_y0  = y0_.deviceVectorAll();
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_dy1 = dy1_.deviceVectorAll();
-	auto d_dy2 = dy2_.deviceVectorAll();
+	auto d_y0  = y0_.deviceViewAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_dy1 = dy1_.deviceViewAll();
+	auto d_dy2 = dy2_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsMoulton4::predict",
@@ -165,13 +165,13 @@ bool pFlow::AdamsMoulton4::intAll(
 	range activeRng)
 {
 
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
 
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_y0  = y0_.deviceVectorAll();
-	auto d_dy1 = dy1_.deviceVectorAll();
-	auto d_dy2 = dy2_.deviceVectorAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_y0  = y0_.deviceViewAll();
+	auto d_dy1 = dy1_.deviceViewAll();
+	auto d_dy2 = dy2_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsMoulton4::correct",

src/Integration/AdamsMoulton4/AdamsMoulton4.hpp

@@ -147,13 +147,13 @@ bool AdamsMoulton4::predictRange(
 	realx3Vector_D& dy,
 	activeFunctor activeP )
 {
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
 
-	auto d_y0 = y0_.deviceVectorAll();
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_dy1 = dy1_.deviceVectorAll();
-	auto d_dy2 = dy2_.deviceVectorAll();
+	auto d_y0 = y0_.deviceViewAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_dy1 = dy1_.deviceViewAll();
+	auto d_dy2 = dy2_.deviceViewAll();
 	
 	auto activeRng = activeP.activeRange();
 
@@ -185,13 +185,13 @@ bool pFlow::AdamsMoulton4::intRange(
 	activeFunctor activeP )
 {
 
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
 
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_y0  = y0_.deviceVectorAll();
-	auto d_dy1 = dy1_.deviceVectorAll();
-	auto d_dy2 = dy2_.deviceVectorAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_y0  = y0_.deviceViewAll();
+	auto d_dy1 = dy1_.deviceViewAll();
+	auto d_dy2 = dy2_.deviceViewAll();
 
 	auto activeRng = activeP.activeRange();
 

src/Integration/AdamsMoulton5/AdamsMoulton5.cpp

@@ -145,14 +145,14 @@ bool pFlow::AdamsMoulton5::predictAll(
 	range activeRng)
 {
 
-	auto d_dy  = dy.deviceVectorAll();
-	auto d_y   = y.deviceVectorAll();
+	auto d_dy  = dy.deviceViewAll();
+	auto d_y   = y.deviceViewAll();
 
-	auto d_y0  = y0_.deviceVectorAll();
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_dy1 = dy1_.deviceVectorAll();
-	auto d_dy2 = dy2_.deviceVectorAll();
-	auto d_dy3 = dy3_.deviceVectorAll();
+	auto d_y0  = y0_.deviceViewAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_dy1 = dy1_.deviceViewAll();
+	auto d_dy2 = dy2_.deviceViewAll();
+	auto d_dy3 = dy3_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsMoulton5::predict",
@@ -178,14 +178,14 @@ bool pFlow::AdamsMoulton5::intAll(
 	range activeRng)
 {
 
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
 
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_y0  = y0_.deviceVectorAll();
-	auto d_dy1 = dy1_.deviceVectorAll();
-	auto d_dy2 = dy2_.deviceVectorAll();
-	auto d_dy3 = dy3_.deviceVectorAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_y0  = y0_.deviceViewAll();
+	auto d_dy1 = dy1_.deviceViewAll();
+	auto d_dy2 = dy2_.deviceViewAll();
+	auto d_dy3 = dy3_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsMoulton5::correct",

src/Integration/AdamsMoulton5/AdamsMoulton5.hpp

@@ -150,14 +150,14 @@ bool AdamsMoulton5::predictRange(
 	realx3Vector_D& dy,
 	activeFunctor activeP )
 {
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
 
-	auto d_y0 = y0_.deviceVectorAll();
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_dy1 = dy1_.deviceVectorAll();
-	auto d_dy2 = dy2_.deviceVectorAll();
-	auto d_dy3 = dy3_.deviceVectorAll();
+	auto d_y0 = y0_.deviceViewAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_dy1 = dy1_.deviceViewAll();
+	auto d_dy2 = dy2_.deviceViewAll();
+	auto d_dy3 = dy3_.deviceViewAll();
 	
 	auto activeRng = activeP.activeRange();
 
@@ -189,14 +189,14 @@ bool pFlow::AdamsMoulton5::intRange(
 	activeFunctor activeP )
 {
 
-	auto d_dy = dy.deviceVectorAll();
-	auto d_y  = y.deviceVectorAll();
+	auto d_dy = dy.deviceViewAll();
+	auto d_y  = y.deviceViewAll();
 
-	auto d_dy0 = dy0_.deviceVectorAll();
-	auto d_y0  = y0_.deviceVectorAll();
-	auto d_dy1 = dy1_.deviceVectorAll();
-	auto d_dy2 = dy2_.deviceVectorAll();
-	auto d_dy3 = dy3_.deviceVectorAll();
+	auto d_dy0 = dy0_.deviceViewAll();
+	auto d_y0  = y0_.deviceViewAll();
+	auto d_dy1 = dy1_.deviceViewAll();
+	auto d_dy2 = dy2_.deviceViewAll();
+	auto d_dy3 = dy3_.deviceViewAll();
 
 	auto activeRng = activeP.activeRange();
 

src/Integration/CMakeLists.txt

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

src/Integration/boundaries/boundaryIntegration.cpp

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

src/Integration/boundaries/boundaryIntegration.hpp

@@ -0,0 +1,91 @@
+
+#ifndef __boundaryIntegration_hpp__
+#define __boundaryIntegration_hpp__
+
+
+#include "generalBoundary.hpp"
+#include "virtualConstructor.hpp"
+#include "pointFields.hpp"
+
+namespace pFlow
+{
+
+class integration;
+
+class boundaryIntegration
+:
+    public generalBoundary
+{
+private:
+
+	const integration& 		integration_;
+
+public:
+
+    TypeInfo("boundaryIntegration<none>");
+
+    boundaryIntegration(
+        const boundaryBase& boundary,
+        const pointStructure& pStruct, 
+        const word& method,
+        integration& intgrtn);
+
+    ~boundaryIntegration()override = default;
+
+    create_vCtor(
+        boundaryIntegration,
+        boundaryBase,
+        (
+            const boundaryBase& boundary,
+            const pointStructure& pStruct, 
+            const word& method,
+			integration& intgrtn
+        ),
+        (boundary, pStruct, method, intgrtn)
+    );
+    
+    add_vCtor(
+        boundaryIntegration,
+        boundaryIntegration,
+        boundaryBase
+    );
+
+    bool hearChanges(
+        real t,
+        real dt,
+        uint32 iter,
+        const message &msg,
+        const anyList &varList) override
+    {
+        return true;
+    }
+
+	const integration& Integration()const
+	{
+		return integration_;
+	}
+
+	virtual
+	bool correct(real dt, const realx3PointField_D& y, const realx3PointField_D& dy)
+	{
+		return true;
+	}
+
+    virtual
+    bool correctPStruct(real dt, const realx3PointField_D& vel)
+    {
+        return true;
+    }
+
+    static
+    uniquePtr<boundaryIntegration> create(
+        const boundaryBase& boundary,
+        const pointStructure& pStruct, 
+        const word& method,
+		integration& intgrtn);
+
+};
+
+}
+
+#endif

src/Integration/boundaries/boundaryIntegrationList.cpp

@@ -0,0 +1,59 @@
+#include "boundaryIntegrationList.hpp"
+#include "integration.hpp"
+
+pFlow::boundaryIntegrationList::boundaryIntegrationList(
+    const pointStructure &pStruct, 
+    const word &method, 
+    integration &intgrtn
+)
+:
+    ListPtr<boundaryIntegration>(6),
+	boundaries_(pStruct.boundaries())
+{
+    
+    for(uint32 i=0; i<6; i++)
+	{
+		this->set(
+			i,
+			boundaryIntegration::create(
+				boundaries_[i],
+				pStruct, 
+				method,
+                intgrtn
+            )
+        );
+	}	
+
+}
+
+bool pFlow::boundaryIntegrationList::correct(real dt, realx3PointField_D &y, realx3PointField_D &dy)
+{
+	for(auto& bndry:*this)
+	{
+		if(!bndry->correct(dt, y, dy))
+		{
+			fatalErrorInFunction<<"Error in correcting boundary "<<
+			bndry->boundaryName()<<endl;
+			return false;
+		}
+		
+    }
+    return true;
+}
+
+bool pFlow::boundaryIntegrationList::correctPStruct(
+	real dt, 
+	pointStructure &pStruct, 
+	const realx3PointField_D &vel)
+{
+	for(auto& bndry:*this)
+	{
+		if(!bndry->correctPStruct(dt, vel))
+		{
+			fatalErrorInFunction<<"Error in correcting boundary "<<
+			bndry->boundaryName()<<" in pointStructure."<<endl;
+			return false;
+		}
+	}
+    return true;
+}

src/Integration/boundaries/boundaryIntegrationList.hpp

@@ -0,0 +1,51 @@
+
+#ifndef __boundaryIntegrationList_hpp__
+#define __boundaryIntegrationList_hpp__
+
+
+#include "boundaryList.hpp"
+#include "ListPtr.hpp"
+#include "boundaryIntegration.hpp"
+
+
+namespace pFlow
+{
+
+class integration;
+
+class boundaryIntegrationList
+:
+    public ListPtr<boundaryIntegration>
+{
+private:
+
+    const boundaryList&         boundaries_;
+
+public:
+
+    boundaryIntegrationList(
+        const pointStructure& pStruct, 
+        const word& method,
+		integration& intgrtn
+    );
+
+    ~boundaryIntegrationList()=default;
+
+    bool correct(
+		real dt, 
+		realx3PointField_D& y, 
+		realx3PointField_D& dy);
+
+    bool correctPStruct(
+        real dt, 
+        pointStructure& pStruct, 
+        const realx3PointField_D& vel);
+
+};
+
+
+
+}
+
+
+#endif //__boundaryIntegrationList_hpp__

src/Integration/integration/integration.cpp

@@ -19,33 +19,35 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "integration.hpp"
+#include "pointStructure.hpp"
+#include "repository.hpp"
 
 pFlow::integration::integration
 (
 	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
+	pointStructure& pStruct,
+	const word&,
+	const realx3Field_D&
 )
 :
-	owner_(owner),
-	baseName_(baseName),
-	pStruct_(pStruct)
-{
-	CONSUME_PARAM(method);
-}
+	owner_(*pStruct.owner()),
+	pStruct_(pStruct),
+	baseName_(baseName)	
+{}
 
 
 pFlow::uniquePtr<pFlow::integration> 
-	pFlow::integration::create(
+	pFlow::integration::create
+(
 		const word& baseName,
-		repository& owner,
-		const pointStructure& pStruct,
-		const word& method)
+		pointStructure& pStruct,
+		const word& method,
+		const realx3Field_D& initialValField
+)
 {
 	if( wordvCtorSelector_.search(method) )
 	{
-		return wordvCtorSelector_[method] (baseName, owner, pStruct, method);
+		return wordvCtorSelector_[method] (baseName, pStruct, method, initialValField);
 	}
 	else
 	{

src/Integration/integration/integration.hpp

@@ -23,14 +23,16 @@ Licence:
 
 
 #include "virtualConstructor.hpp"
-#include "Vectors.hpp"
-#include "pointStructure.hpp"
-#include "repository.hpp"
+#include "pointFields.hpp"
 
 
 namespace pFlow
 {
 
+// - Forward
+class repository;
+class pointStructure;
+
 /**
  * Base class for integrating the first order ODE (IVP) 
  * 
@@ -48,19 +50,19 @@ namespace pFlow
  */
 class integration
 {
-protected:
+private:
 
 	// - Protected data members 
 
 		/// The owner repository that all fields are storred in
 		repository& owner_;
 
-		/// The base name for integration 
-		const word baseName_;
-
 		/// A reference to pointStructure
 		const pointStructure& pStruct_;
 
+		/// The base name for integration 
+		const word baseName_;
+
 public:
 
 	/// Type info
@@ -72,9 +74,9 @@ public:
 		/// Construct from components
 		integration(
 			const word& baseName,
-			repository& owner,
-			const pointStructure& pStruct,
-			const word& method);
+			pointStructure& pStruct,
+			const word& method,
+			const realx3Field_D& initialValField);
 
 		/// Copy constructor 
 		integration(const integration&) = default;
@@ -88,22 +90,22 @@ public:
 		/// Move assignment
 		integration& operator = (integration&&) = default;
 
-		/// Polymorphic copy/cloning
-		virtual 
-		uniquePtr<integration> clone()const=0;
-
 		/// Destructor 
 		virtual ~integration()=default;
 
 		/// Add a virtual constructor 
-		create_vCtor(
+		create_vCtor
+		(
 			integration,
 			word,
-			(const word& baseName,
-			repository& owner,
-			const pointStructure& pStruct,
-			const word& method),
-			(baseName, owner, pStruct, method) );
+			(
+				const word& baseName,
+				pointStructure& pStruct,
+				const word& method,
+				const realx3Field_D& initialValField
+			),
+			(baseName, pStruct, method, initialValField)
+		);
 
 
 	// - Methods
@@ -129,13 +131,25 @@ public:
 			return owner_;
 		}
 
+		virtual 
+		void updateBoundariesSlaveToMasterIfRequested() = 0;
+		/// return integration method 
+		virtual 
+		word method()const = 0 ;
+
 		/// Prediction step in integration
 		virtual 
-		bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
+		bool predict(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0;
+
+		virtual 
+		bool predict(real dt, realx3Field_D& y, realx3PointField_D& dy) = 0;
 
 		/// Correction/main integration step
 		virtual 
-		bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
+		bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0;
+
+		virtual 
+		bool correctPStruct(real dt, pointStructure& pStruct, realx3PointField_D& vel) = 0;
 
 		/// Set the initial values for new indices 
 		virtual 
@@ -152,9 +166,9 @@ public:
 	static
 		uniquePtr<integration> create(
 			const word& baseName,
-			repository& owner,
-			const pointStructure& pStruct,
-			const word& method);
+			pointStructure& pStruct,
+			const word& method,
+			const realx3Field_D& initialValField);
 
 }; // integration
 

src/Interaction/CMakeLists.txt

@@ -1,11 +1,43 @@
 
 set(SourceFiles 
+contactSearch/methods/cellBased/NBS/mapperNBS.cpp
+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/ContactSearch/ContactSearchs.cpp
+
 interaction/interaction.cpp
-sphereInteraction/sphereInteractions.cpp
+sphereInteraction/sphereInteractionsLinearModels.cpp
+sphereInteraction/sphereInteractionsNonLinearModels.cpp
+sphereInteraction/sphereInteractionsNonLinearModModels.cpp
+
+grainInteraction/grainInteractionsLinearModels.cpp
+grainInteraction/grainInteractionsNonLinearModels.cpp
+grainInteraction/grainInteractionsNonLinearModModels.cpp
+
+
 )
 
+if(pFlow_Build_MPI)
+    list(APPEND SourceFiles 
+    contactSearch/boundaries/processorBoundaryContactSearch/processorBoundaryContactSearch.cpp
+    sphereInteraction/boundaries/processorBoundarySphereInteraction/processorBoundarySphereInteractions.cpp
+    contactSearch/boundaries/twoPartContactSearch/twoPartContactSearchKernels.cpp
+    contactSearch/boundaries/twoPartContactSearch/twoPartContactSearch.cpp
+    )
+endif()
+
 set(link_libs Kokkos::kokkos phasicFlow Property Particles Geometry)
 
 pFlow_add_library_install(Interaction SourceFiles link_libs)

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/contactForce/linearCF.hpp

@@ -133,8 +133,9 @@ protected:
 			return false;
 		}
 
-		realVector etha_n(nElem);
-		realVector etha_t(nElem);
+		realVector etha_n("etha_n", nElem);
+		realVector etha_t("etha_t", nElem);
+
 		ForAll(i , kn)
 		{
 			etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/
@@ -144,7 +145,7 @@ protected:
 					sqrt(pow(log(et[i]),2.0)+ pow(Pi,2.0));
 		}
 
-		Vector<linearProperties> prop(nElem);
+		Vector<linearProperties> prop("prop", nElem);
 		ForAll(i,kn)
 		{
 			prop[i] = {kn[i], kt[i], etha_n[i], etha_t[i], mu[i]};
@@ -219,10 +220,10 @@ public:
 	void contactForce
 	(
 		const real dt,
-		const int32 i,
-		const int32 j,
-		const int32 propId_i,
-		const int32 propId_j,
+		const uint32 i,
+		const uint32 j,
+		const uint32 propId_i,
+		const uint32 propId_j,
 		const real Ri,
 		const real Rj,
 		const real ovrlp_n,

src/Interaction/Models/contactForce/nonLinearCF.hpp

@@ -121,7 +121,7 @@ protected:
 		}
 
 
-		realVector etha_n(nElem);
+		realVector etha_n("etha_n",nElem);
 		
 		ForAll(i , en)
 		{
@@ -137,7 +137,7 @@ protected:
 
 		}
 
-		Vector<nonLinearProperties> prop(nElem);
+		Vector<nonLinearProperties> prop("prop",nElem);
 		ForAll(i,Yeff)
 		{
 			prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]};
@@ -214,10 +214,10 @@ public:
 	void contactForce
 	(
 		const real dt,
-		const int32 i,
-		const int32 j,
-		const int32 propId_i,
-		const int32 propId_j,
+		const uint32 i,
+		const uint32 j,
+		const uint32 propId_i,
+		const uint32 propId_j,
 		const real Ri,
 		const real Rj,
 		const real ovrlp_n,

src/Interaction/Models/contactForce/nonLinearMod.hpp

@@ -121,7 +121,7 @@ protected:
 		}
 
 
-		Vector<nonLinearProperties> prop(nElem);
+		Vector<nonLinearProperties> prop("prop",nElem);
 		ForAll(i,Yeff)
 		{
 			prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]};
@@ -198,10 +198,10 @@ public:
 	void contactForce
 	(
 		const real dt,
-		const int32 i,
-		const int32 j,
-		const int32 propId_i,
-		const int32 propId_j,
+		const uint32 i,
+		const uint32 j,
+		const uint32 propId_i,
+		const uint32 propId_j,
 		const real Ri,
 		const real Rj,
 		const real ovrlp_n,

src/Interaction/Models/grainContactForceModels.hpp

@@ -18,30 +18,28 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+#ifndef __grainContactForceModels_hpp__
+#define __grainContactForceModels_hpp__
 
-#ifndef __Control_hpp__
-#define __Control_hpp__
+#include "cGAbsoluteLinearCF.hpp"
+#include "cGRelativeLinearCF.hpp"
+
+#include "grainRolling.hpp"
 
 
-// top-level entity repository for the whole application 
-// Each application that is executed in pFlow, should have 
-// settings/systemControl file in it. 
-// This  repository holds two main repositories: Time and settings 
 
-#include "systemControl.hpp"
-#include "timeFolder.hpp"
-
-namespace pFlow
+namespace pFlow::cfModels 
 {
 
 
-inline systemControl& Control()
-{
-	static systemControl control_;
-	return control_;
+using limitedCGAbsoluteLinearGrainRolling 	= grainRolling<cGAbsoluteLinear<true>>;
+using nonLimitedCGAbsoluteLinearGrainRolling = grainRolling<cGAbsoluteLinear<false>>;
+
+using limitedCGRelativeLinearGrainRolling   = grainRolling<cGRelativeLinear<true>>;
+using nonLimitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<false>>;
+
 }
 
-} // pFlow
 
 
-#endif // __Control_hpp__
+#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/Models/rolling/normalRolling.hpp

@@ -76,10 +76,10 @@ public:
 	void rollingFriction
 	(
 		const real dt,
-		const int32 i,
-		const int32 j,
-		const int32 propId_i,
-		const int32 propId_j,
+		const uint32 i,
+		const uint32 j,
+		const uint32 propId_i,
+		const uint32 propId_j,
 		const real Ri,
 		const real Rj,
 		const realx3& wi,

src/Interaction/contactLists/sortedContactList.hpp

@@ -53,7 +53,7 @@ protected:
 
 	ViewType1D<ValueType,ExecutionSpace> 	values_;
 
-	int32 									size0_ = 0;
+	uint32 									size0_ = 0;
 
 	ViewType1D<PairType,ExecutionSpace> 	sortedPairs0_;
 
@@ -73,7 +73,7 @@ protected:
 	using rpReFillPairs = Kokkos::RangePolicy<
 		ExecutionSpace,
 		Kokkos::Schedule<Kokkos::Static>,
-		Kokkos::IndexType<int32>,
+		Kokkos::IndexType<uint32>,
 		TagReFillPairs>;
 
 public:
@@ -81,7 +81,7 @@ public:
 	TypeInfoNV("sortedContactList");
 
 
-	sortedContactList(int32 initialSize =1)
+	explicit sortedContactList(uint32 initialSize =1)
 	:
 		SortedPairs(initialSize),
 		values_("values", SortedPairs::capacity()),
@@ -114,31 +114,31 @@ public:
 	}
 
 	INLINE_FUNCTION_HD
-	ValueType getValue(int32 idx)const
+	ValueType getValue(uint32 idx)const
 	{
 		return values_[idx];
 	}
 
 	INLINE_FUNCTION_HD
-	void setValue(int32 idx, const ValueType& val)const
+	void setValue(uint32 idx, const ValueType& val)const
 	{
 		values_[idx] = val;
 	}
 
 	INLINE_FUNCTION_HD
-	void operator()(TagReFillPairs, int32 idx)const
+	void operator()(TagReFillPairs, uint32 idx)const
 	{
 
-		auto searchLen = max(size0_/1000,10);
-		auto start 	= max(0,idx-searchLen);
-		auto end 	= min(size0_,idx+searchLen);
+		uint32 searchLen = max(size0_/1000u,10u);
+		uint32 start 	= idx-min(searchLen,idx);
+		uint32 end 	= min(size0_,idx+searchLen);
 		auto newPair = this->sortedPairs_[idx];
 		if(	auto idx0 = binarySearch(
 				sortedPairs0_,
 				start,
 				end,
 				newPair);
-				idx0>=0)
+				idx0!=static_cast<uint32>(-1))
 		{
 			values_[idx] = values0_[idx0];
 		}
@@ -147,7 +147,7 @@ public:
 			start, 
 			end,
 			newPair);
-			idx0>=0)
+			idx0!= static_cast<uint32>(-1) )
 		{
 			values_[idx] = values0_[idx0];
 			

src/Interaction/contactLists/sortedPairs.hpp

@@ -52,24 +52,24 @@ public:
 	{
 		using PairType = typename sortedPairs::PairType;
 
-		int32 size_;
+		uint32 size_;
 
 		ViewType1D<PairType,ExecutionSpace> sortedParis_;
 
 		INLINE_FUNCTION_HD
-		int32 size()const { return size_; }
+		uint32 size()const { return size_; }
 
 		INLINE_FUNCTION_HD
-		int32 loopCount()const { return size_; }
+		uint32 loopCount()const { return size_; }
 
 		INLINE_FUNCTION_HD
-		bool isValid(int32 i)const { return i<size_; }
+		bool isValid(uint32 i)const { return i<size_; }
 
 		INLINE_FUNCTION_HD
-		PairType getPair(int i)const {	return sortedParis_[i];	}
+		PairType getPair(uint32 i)const {	return sortedParis_[i];	}
 
 		INLINE_FUNCTION_HD
-		bool getPair(int32 i, PairType& pair)const {
+		bool getPair(uint32 i, PairType& pair)const {
 			if(i<size_)	{
 				pair = sortedParis_[i];
 				return true;
@@ -85,22 +85,22 @@ public:
 protected:
 
 	/// size of pair list
-	int32 			size_ = 0;
+	uint32 			size_ = 0;
 
-	ViewType1D<int32,ExecutionSpace>		flags_;
+	ViewType1D<uint32,ExecutionSpace>		flags_;
 
 	ViewType1D<PairType,ExecutionSpace> 	sortedPairs_;
 
 	using rpFillFlag = Kokkos::RangePolicy<
 		ExecutionSpace,
 		Kokkos::Schedule<Kokkos::Static>,
-		Kokkos::IndexType<int32>,
+		Kokkos::IndexType<uint32>,
 		TagFillFlag	>;
 
 	using rpFillPairs = Kokkos::RangePolicy<
 		ExecutionSpace,
 		Kokkos::Schedule<Kokkos::Static>,
-		Kokkos::IndexType<int32>,
+		Kokkos::IndexType<uint32>,
 		TagFillPairs>;
 	
 public:
@@ -110,7 +110,7 @@ public:
 
 
 	// constructors
-	sortedPairs(int32 initialSize =1)
+	explicit sortedPairs(uint32 initialSize =1)
 	:
 		UnsortedPairs(initialSize),
 		flags_("flags_",UnsortedPairs::capacity()+1),
@@ -134,7 +134,7 @@ public:
 	//   return the pair at index idx
 	//   perform no check for size and existance
 	INLINE_FUNCTION_HD
-	PairType getPair(int32 idx)const
+	PairType getPair(uint32 idx)const
 	{
 		return sortedPairs_[idx];
 	}
@@ -142,7 +142,7 @@ public:
 	// - Device/host call
 	//   return the pair at index idx
 	INLINE_FUNCTION_HD
-	bool getPair(int32 idx, PairType& p)const
+	bool getPair(uint32 idx, PairType& p)const
 	{
 		if(isValid(idx))
 		{
@@ -156,7 +156,7 @@ public:
 	}
 
 	INLINE_FUNCTION_HD
-	bool isValid(int32 idx)const
+	bool isValid(uint32 idx)const
 	{
 		return idx < size_;	
 	}
@@ -164,12 +164,12 @@ public:
 
 	//use this when the value of size_ is updated
 	INLINE_FUNCTION_H
-	int32 size()const
+	uint32 size()const
 	{
 		return size_;
 	}
 
-	int32 loopCount()const
+	uint32 loopCount()const
 	{
 		return size_;
 	}
@@ -189,7 +189,7 @@ public:
 	void prepareSorted()
 	{
 		// first check the size of flags_ 
-		int32 capacity = UnsortedPairs::capacity();
+		uint32 capacity = UnsortedPairs::capacity();
 
 		if( capacity+1 > flags_.size() )
 		{
@@ -218,7 +218,7 @@ public:
 		if( size_>sortedPairs_.size() )
 		{
 			// get more space to prevent reallocations in next iterations
-			int32 len = size_*1.1+1; 
+			uint32 len = size_*1.1+1; 
 			reallocNoInit(sortedPairs_, len);
 		}
 		
@@ -235,7 +235,7 @@ public:
 	}
 
 	INLINE_FUNCTION_HD
-	void operator()(TagFillFlag, int32 i)const
+	void operator()(TagFillFlag, uint32 i)const
 	{
 		if(this->container_.valid_at(i) )
 			flags_[i] = 1;
@@ -244,7 +244,7 @@ public:
 	}
 
 	INLINE_FUNCTION_HD
-	void operator()(TagFillPairs, int32 i)const
+	void operator()(TagFillPairs, uint32 i)const
 	{
 		auto fi = flags_[i];
 		if(fi!=flags_[i+1])

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
 {
 
@@ -72,7 +73,7 @@ protected:
 	using rpFillPairs = Kokkos::RangePolicy<
 		ExecutionSpace,
 		Kokkos::Schedule<Kokkos::Static>,
-		Kokkos::IndexType<int32>,
+		Kokkos::IndexType<uint32>,
 		TagReFillPairs>;
 
 
@@ -80,7 +81,7 @@ public:
 
 	TypeInfoNV("unsortedContactList");
 
-	unsortedContactList(int32 capacity=1)
+	explicit unsortedContactList(uint32 capacity=1)
 	:
 		UnsortedPairs(capacity),
 		values_("values", UnsortedPairs::capacity()),
@@ -94,7 +95,8 @@ public:
 		// swap conainer and values 
 		swapViews(values0_, values_);
 		swapViews(container0_, this->container_);
-		return UnsortedPairs::beforeBroadSearch();
+		UnsortedPairs::beforeBroadSearch();
+		return true;
 	}
 
 	bool afterBroadSearch()
@@ -122,7 +124,7 @@ public:
 	INLINE_FUNCTION_HD
 	bool getValue(const PairType& p, ValueType& val)const
 	{
-		if(auto idx = this->find(p); idx>=0)
+		if(auto idx = this->find(p); idx!=static_cast<uint32>(-1))
 		{
 			val = getValue(idx); 
 			return true;
@@ -131,7 +133,7 @@ public:
 	}
 
 	INLINE_FUNCTION_HD
-	void setValue(int32 idx, const ValueType& val)const
+	void setValue(uint32 idx, const ValueType& val)const
 	{
 		values_[idx] = val;
 	}
@@ -139,7 +141,7 @@ public:
 	INLINE_FUNCTION_HD
 	bool setValue(const PairType& p, const ValueType& val)const
 	{
-		if(auto idx = this->find(p); idx>=0)
+		if(uint32 idx = this->find(p); idx!=static_cast<uint32>(-1))
 		{
 			setValue(idx, val);
 			return true;;
@@ -148,13 +150,13 @@ public:
 	}
 
 	INLINE_FUNCTION_HD
-	void operator()(TagReFillPairs, int32 idx)const
+	void operator()(TagReFillPairs, uint32 idx)const
 	{
 		if( this->isValid(idx) )
 		{
-			if( int32 idx0 = 
+			if( uint32 idx0 = 
 					container0_.find(this->getPair(idx));
-					idx0>=0 )
+					idx0!= static_cast<uint32>(-1) )
 			{
 				values_[idx] = values0_[idx0];
 			}

src/Interaction/contactLists/unsortedPairs.hpp

@@ -41,7 +41,7 @@ public:
 
 	using memory_space 		= typename ExecutionSpace::memory_space;
 
-	using PairType 			= kPair<idType,idType>;
+	using PairType 			= Pair<idType,idType>;
 
 	using ContainerType 	= unorderedSet<PairType, ExecutionSpace>;
 
@@ -52,19 +52,19 @@ public:
 		ContainerType 		Container_;
 
 		INLINE_FUNCTION_HD
-		int32 size()const { return Container_.size(); }
+		uint32 size()const { return Container_.size(); }
 
 		INLINE_FUNCTION_HD
-		int32 loopCount()const { return Container_.capacity(); }
+		uint32 loopCount()const { return Container_.capacity(); }
 
 		INLINE_FUNCTION_HD
-		bool isValid(int32 idx)const { return Container_.valid_at(idx); }
+		bool isValid(uint32 idx)const { return Container_.valid_at(idx); }
 
 		INLINE_FUNCTION_HD
-		PairType getPair(int idx)const {	return Container_.key_at(idx); }
+		PairType getPair(uint32 idx)const {	return Container_.key_at(idx); }
 
 		INLINE_FUNCTION_HD
-		bool getPair(int32 idx, PairType& pair)const {
+		bool getPair(uint32 idx, PairType& pair)const {
 			if(Container_.valid_at(idx))	{
 				pair = Container_.key_at(idx);
 				return true;
@@ -84,7 +84,7 @@ public:
 	TypeInfoNV("unsorderedPairs");
 
 	// constructor
-	unsortedPairs(int32 capacity=1)
+	explicit unsortedPairs(uint32 capacity=1)
 	:
 		container_(capacity) // the minimum capacity would be 128 
 	{}
@@ -102,20 +102,20 @@ public:
 	
 	// - Device call
 	INLINE_FUNCTION_HD
-	int32 insert(idType i, idType j)const
+	uint32 insert(idType i, idType j)const
 	{
 		if(auto insertResult = container_.insert(PairType(i,j)); insertResult.failed())
-			return -1;
+			return static_cast<uint32>(-1);
 		else
 			return insertResult.index();
 			
 	}
 
 	INLINE_FUNCTION_HD
-	int32 insert(const PairType& p)const
+	uint32 insert(const PairType& p)const
 	{
 		if(auto insertResult = container_.insert(p); insertResult.failed())
-			return -1;
+			return static_cast<uint32>(-1);
 		else
 			return insertResult.index();
 
@@ -125,7 +125,7 @@ public:
 	//   return the pair at index idx
 	//   perform no check for size and existance
 	INLINE_FUNCTION_HD
-	PairType getPair(int32 idx)const
+	PairType getPair(uint32 idx)const
 	{
 		return container_.key_at(idx);
 	}
@@ -133,7 +133,7 @@ public:
 	// - Device call
 	//   return the pair at index idx
 	INLINE_FUNCTION_HD
-	bool getPair(int32 idx, PairType& p)const
+	bool getPair(uint32 idx, PairType& p)const
 	{
 		if(container_.valid_at(idx))
 		{
@@ -148,36 +148,36 @@ public:
 	}
 
 	INLINE_FUNCTION_HD
-	int32 find(const PairType & p)const
+	uint32 find(const PairType & p)const
 	{
 		if( auto idx = container_.find(p); 
 			idx != Kokkos::UnorderedMapInvalidIndex )
 			return idx;
 		else
-			return -1;
+			return static_cast<uint32>(-1);
 	}
 
 	INLINE_FUNCTION_HD
-	bool isValid(int32 idx)const
+	bool isValid(uint32 idx)const
 	{
 		return container_.valid_at(idx);	
 	}
 
 
 	INLINE_FUNCTION_HD
-	int32 capacity() const
+	uint32 capacity() const
 	{
 		return container_.capacity();
 	}
 
-	int32 loopCount()const
+	uint32 loopCount()const
 	{
 		return container_.capacity();
 	}
 	
 	//use this when the value of size_ is updated
 	INLINE_FUNCTION_H
-	int32 size()const
+	uint32 size()const
 	{
 		return container_.size();
 	}
@@ -190,7 +190,7 @@ public:
 	/// increase the capacity of the container by at-least len 
 	/// the content will be erased. 
 	INLINE_FUNCTION_H
-	void increaseCapacityBy(int32 len)
+	void increaseCapacityBy(uint32 len)
 	{
 		uint newCap = container_.capacity()+len;
 		this->clear();

src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp

@@ -22,16 +22,18 @@ Licence:
 #ifndef __ContactSearch_hpp__
 #define __ContactSearch_hpp__
 
-
+#include "contactSearchGlobals.hpp"
 #include "contactSearch.hpp"
 #include "box.hpp"
+#include "particles.hpp"
+#include "geometry.hpp"
+#include "boundaryContactSearchList.hpp"
 
 namespace pFlow
 {
 
 template<
-	template<class> class BaseMethod,
-	template<class> class WallMapping	
+	class searchMethod
 >
 class ContactSearch
 :
@@ -39,207 +41,134 @@ class ContactSearch
 {
 public:
 	
-	using IdType 			= typename contactSearch::IdType;
+	using IdType 			= uint32;
 
-	using IndexType 		= typename contactSearch::IndexType;
+	using ExecutionSpace 	= DefaultExecutionSpace;
 
-	using ExecutionSpace 	= typename contactSearch::ExecutionSpace;
+	using SearchMethodType = searchMethod;
 
-	using PairContainerType = typename contactSearch::PairContainerType;
+private:
 	
-	using ParticleContactSearchType = 
-		BaseMethod<
-			ExecutionSpace>;
+	uniquePtr<SearchMethodType> 	ppwContactSearch_ = nullptr;
 
-	using WallMappingType =
-		WallMapping<
-			ExecutionSpace>;
+	boundaryContactSearchList 		csBoundaries_;
 
-protected:
+	bool BroadSearch(
+		const timeInfo& ti,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force = false) override
+	{
 		
+		const auto& position = Particles().pointPosition().deviceViewAll();
+		const auto& flags = Particles().dynPointStruct().activePointsMaskDevice();
+		const auto& diam = Particles().boundingSphere().deviceViewAll();
 
-	uniquePtr<ParticleContactSearchType> particleContactSearch_ = nullptr;
+		return ppwContactSearch_().broadSearch(
+			ti.iter(),
+			ti.t(),
+			ti.dt(),
+			ppPairs, 
+			pwPairs,
+			position, 
+			flags, 
+			diam,
+			force
+		);
+	}
 
-	uniquePtr<WallMappingType> 			 wallMapping_ = nullptr;
+	bool BoundaryBroadSearch(
+		uint32 bndryIndex,
+		const timeInfo& ti,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force = false)override
+	{		
+		return csBoundaries_[bndryIndex].broadSearch(
+			ti.iter(), 
+			ti.t(), 
+			ti.dt(), 
+			ppPairs, 
+			pwPairs, 
+			force
+		);
+	}
 
 public:
 
-	TypeInfoTemplate2("ContactSearch", ParticleContactSearchType, WallMappingType);
+	TypeInfoTemplate11("ContactSearch", SearchMethodType);
 
 	ContactSearch(
 		const dictionary& csDict,
-		const box& domain,
+		const box& extDomain,
 	 	const particles& prtcl,
 	 	const geometry& geom,
 	 	Timers& timers)
 	:
-		contactSearch(csDict, domain, prtcl, geom, timers)
-		
+		contactSearch(
+			csDict, 
+			extDomain, 
+			prtcl, 
+			geom, 
+			timers),
+		csBoundaries_(
+			csDict, 
+			Particles().pStruct().boundaries(),
+			*this)
 	{
 
-		auto method = dict().getVal<word>("method");
-		auto wmMethod = dict().getVal<word>("wallMapping");
-		
-		auto nbDict = dict().subDict(method+"Info");
-
-		real minD, maxD;
+		real minD;
+		real maxD;
 		this->Particles().boundingSphereMinMax(minD, maxD);
 		
-		const auto& position = this->Particles().pointPosition().deviceVectorAll();
-		const auto& diam = this->Particles().boundingSphere().deviceVectorAll();
+		const auto& position = this->Particles().pointPosition().deviceViewAll();
+		const auto& flags = this->Particles().dynPointStruct().activePointsMaskDevice();
+		const auto& diam = this->Particles().boundingSphere().deviceViewAll();
 
-		particleContactSearch_ = 
-			makeUnique<ParticleContactSearchType>
+		uint32 wnPoints = this->Geometry().numPoints();
+		uint32 wnTri    = this->Geometry().size();
+		const auto& wPoints = this->Geometry().points().deviceViewAll();
+		const auto& wVertices = this->Geometry().vertices().deviceViewAll();
+		const auto& wNormals = this->Geometry().normals().deviceViewAll();
+
+		ppwContactSearch_ = 
+			makeUnique<SearchMethodType>
 			(
-				nbDict,
-				this->domain(),
+				dict(),
+				this->extendedDomainBox(),
 				minD,
 				maxD,
 				position,
-				diam
-			);
-		REPORT(2)<<"Contact search algorithm for particle-particle is "<<
-				 greenText(particleContactSearch_().typeName())<<endREPORT;
-
-
-		auto wmDict = dict().subDict(wmMethod+"Info");
-
-		int32 wnPoints = this->Geometry().numPoints();
-		int32 wnTri    = this->Geometry().size();
-
-		const auto& wPoints = this->Geometry().points().deviceVectorAll();
-		const auto& wVertices = this->Geometry().vertices().deviceVectorAll();
-
-		wallMapping_ = 
-			makeUnique<WallMappingType>(
-				wmDict,
-				particleContactSearch_().numLevels(),
-				particleContactSearch_().getCellsLevels(),
+				flags,
+				diam,
 				wnPoints,
 				wnTri,
 				wPoints,
-				wVertices
-				);
-		REPORT(2)<<"Wall mapping algorithm for particle-wall is "<<
-				 greenText(wallMapping_().typeName())<< endREPORT;
-
+				wVertices,
+				wNormals
+			);
 	}
 
-
 	add_vCtor(
 		contactSearch,
 		ContactSearch,
 		dictionary);
 
-	bool broadSearch(
-		PairContainerType& ppPairs,
-		PairContainerType& pwPairs,
-		bool force = false) override
+	bool enterBroadSearchBoundary(const timeInfo& ti, bool force=false)const override
 	{
-
-		
-		if(particleContactSearch_)
-		{
-			auto activeRange = this->Particles().activeRange();
-
-			sphereSphereTimer_.start();
-
-			if(this->Particles().allActive())
-			{		
-				particleContactSearch_().broadSearch(ppPairs, activeRange, force);
-			}
-			else
-			{		
-				particleContactSearch_().broadSearch(ppPairs, activeRange, this->Particles().activePointsMaskD(), force);
-			}
-
-			sphereSphereTimer_.end();
-
-		}
-		else
-			return false;
-		
-		if(wallMapping_)
-		{
-			sphereWallTimer_.start();
-			wallMapping_().broadSearch(pwPairs, particleContactSearch_(), force);
-			sphereWallTimer_.end();
-		}
-		else
-			return false;
-
-		
-		
-		return true;
+		return enterBroadSearch(ti, force) || csBoundaries_.boundariesUpdated();
 	}
 
-
-	bool ppEnterBroadSearch()const override
+	real sizeRatio()const override
 	{
-		if(particleContactSearch_)
-		{
-			return particleContactSearch_().enterBoadSearch(); 
-		}
-		return false;
+		return ppwContactSearch_().sizeRatio();
 	}
 
-	bool pwEnterBroadSearch()const override
+	real cellExtent()const override
 	{
-		if(wallMapping_)
-		{
-			return wallMapping_().enterBoadSearch();
-		}
-		return false;
+		return ppwContactSearch_().cellExtent();
 	}
 	
-
-	bool ppPerformedBroadSearch()const override
-	{
-		if(particleContactSearch_)
-		{
-			return particleContactSearch_().performedSearch();
-		}
-		return false;
-	}
-
-	
-	bool pwPerformedBroadSearch()const override
-	{
-		if(wallMapping_)
-		{
-			return wallMapping_().performedSearch();
-		}
-		return false;
-	}
-
-	/*bool update(const eventMessage& msg) 
-	{
-		if(msg.isSizeChanged() ) 
-		{
-			auto newSize = this->prtcl().size();
-			if(!particleContactSearch_().objectSizeChanged(newSize))
-			{
-				fatalErrorInFunction<<
-				"erro in changing the size for particleContactSearch_ \n";
-				return false;
-			}
-		}
-
-		if(msg.isCapacityChanged() ) 
-		{
-			auto newSize = this->prtcl().capacity();
-			if(!particleContactSearch_().objectSizeChanged(newSize))
-			{
-				fatalErrorInFunction<<
-				"erro in changing the capacity for particleContactSearch_ \n";
-				return false;
-			}
-		}	
-		
-		return true;
-	}*/
-	
-	
 };
 
 }

src/Interaction/contactSearch/ContactSearch/ContactSearchs.cpp

@@ -20,11 +20,11 @@ Licence:
 
 #include "ContactSearch.hpp"
 
-#include "cellMapping.hpp"
+//#include "cellMapping.hpp"
 #include "NBS.hpp"
-#include "multiGridNBS.hpp"
-#include "multiGridMapping.hpp"
+//#include "multiGridNBS.hpp"
+//#include "multiGridMapping.hpp"
 
 
-template class pFlow::ContactSearch<pFlow::NBS, pFlow::cellMapping>;
-template class pFlow::ContactSearch<pFlow::multiGridNBS, pFlow::multiGridMapping>;
+template class pFlow::ContactSearch<pFlow::NBS>;
+//template class pFlow::ContactSearch<pFlow::multiGridNBS, pFlow::multiGridMapping>;

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) )
+	{
+		pOutput.space(4)<<"Creating contact search boundary "<< Green_Text(bType)<<
+		" for "<<boundary.name()<<endl;
+		return boundaryBasevCtorSelector_[bType](dict, boundary, cSearch);
+	}
+	else if(boundaryBasevCtorSelector_.search(altBType))
+	{
+		pOutput.space(4)<<"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,109 @@
+/*------------------------------- 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_;
+    }
+
+    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,32 @@
+#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);
+}
+
+bool pFlow::boundaryContactSearchList::boundariesUpdated() const
+{
+    return boundaries_.boundariesUpdated();
+}

src/Interaction/contactSearch/boundaries/boundaryContactSearchList.hpp

@@ -0,0 +1,43 @@
+#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;
+
+    inline
+    const boundaryList& boundaries()const
+    {
+        return boundaries_;
+    } 
+
+    bool boundariesUpdated()const;
+    
+};
+
+
+}

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,86 @@
+/*------------------------------- 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!=static_cast<uint32>(-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) == static_cast<uint32>(-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,144 @@
+#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();
+    const auto& numElements = numElements_;
+    const auto& elementBox = elementBox_;
+    const auto& validBox = validBox_;
+
+    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) == static_cast<uint32>(-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/contactSearch/contactSearch.cpp

@@ -19,54 +19,124 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "contactSearch.hpp"
-
+#include "streams.hpp"
+#include "particles.hpp"
 
 
 pFlow::contactSearch::contactSearch(
 	const dictionary& dict,
-	const box& domain,
+	const box& extDomain,
  	const particles& prtcl,
  	const geometry& geom,
  	Timers& timers)
 :
-	interactionBase(prtcl, geom),
-	domain_(domain),
-	dict_(dict),
-	sphereSphereTimer_("particle-particle contact search", &timers),
-	sphereWallTimer_("particle-wall contact search", &timers)
+	extendedDomainBox_(extDomain),
+	particles_(prtcl),
+	geometry_(geom),
+	bTimer_("Boundary particles contact search", &timers),
+	ppTimer_("Internal particles contact search", &timers),
+	dict_(dict)
 {
 
 }
 
+const pFlow::pointStructure &pFlow::contactSearch::pStruct() const
+{
+   return particles_.pStruct();
+}
 
-pFlow::uniquePtr<pFlow::contactSearch> pFlow::contactSearch::create(
-	const dictionary& dict,
-	const box& domain,
- 	const particles& prtcl,
- 	const geometry& geom,
- 	Timers& timers)
+bool pFlow::contactSearch::broadSearch
+(
+	const timeInfo &ti, 
+	csPairContainerType &ppPairs, 
+	csPairContainerType &pwPairs, 
+	bool force
+)
 {
 
+	if(enterBroadSearch(ti, force))
+	{
+		ppTimer_.start();
+		if( !BroadSearch(
+			ti,
+			ppPairs,
+			pwPairs,
+			force ) )
+		{
+			fatalErrorInFunction;
+			performedSearch_ = false;
+			return false;
+		}
+		ppTimer_.end();
+		performedSearch_ = true;
+	}
+	else
+	{
+		performedSearch_ = false;
+	}
+
+	return true;
+}
+
+bool pFlow::contactSearch::boundaryBroadSearch
+(
+	uint32 bndryIndex, 
+	const timeInfo &ti, 
+	csPairContainerType &ppPairs, 
+	csPairContainerType &pwPairs, 
+	bool force
+)
+{
+	if(enterBroadSearchBoundary(ti, force))
+	{
+		bTimer_.start();
+		for(uint32 i=0u; i<6u; i++)
+		{
+			if(!BoundaryBroadSearch(
+				i,
+				ti,
+				ppPairs,
+				pwPairs,
+				force))
+			{
+				performedSearchBoundary_ = false;
+				return false;
+			}
+		}
+		bTimer_.end();
+		performedSearchBoundary_ = true;
+	}
+	else
+	{
+
+		performedSearchBoundary_ = false;
+	}
+
+    return true;
+}
+
+pFlow::uniquePtr<pFlow::contactSearch> pFlow::contactSearch::create(
+    const dictionary &dict,
+    const box &extDomain,
+    const particles &prtcl,
+    const geometry &geom,
+    Timers &timers)
+{
 	word baseMethName	= dict.getVal<word>("method");	
-	word wallMethod 	= dict.getVal<word>("wallMapping");
-
-	auto model = angleBracketsNames2("ContactSearch", baseMethName, wallMethod);
-
-
-	REPORT(1)<<"Selecting contact search model . . ."<<endREPORT;
 
+	auto model = angleBracketsNames("ContactSearch", baseMethName);
 	
+	pOutput.space(2)<<"Selecting contact search model "<<Green_Text(model)<<endl;	
 	if( dictionaryvCtorSelector_.search(model))
 	{
-		auto objPtr = dictionaryvCtorSelector_[model] (dict, domain, prtcl, geom, timers);
-		REPORT(2)<<"Model "<< greenText(model)<<" is created."<<endREPORT;
+		auto objPtr = dictionaryvCtorSelector_[model] (dict, extDomain, prtcl, geom, timers);
 		return objPtr;
 	}
 	else
 	{
 		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

@@ -23,42 +23,68 @@ Licence:
 #define __contactSearch_hpp__
 
 
-#include "interactionBase.hpp"
-#include "unsortedPairs.hpp"
-#include "box.hpp"
+#include "contactSearchGlobals.hpp"
 #include "dictionary.hpp"
+#include "virtualConstructor.hpp"
+#include "timeInfo.hpp"
+#include "Timer.hpp"
 
 namespace pFlow
 {
 
+// - forward 
+class box;
+class particles;
+class geometry;
+class pointStructure;
+
 
 class contactSearch
-:
-	public interactionBase
 {
-public:
-	using IdType 			= typename interactionBase::IdType;
+private:
 
-	using IndexType 		= typename interactionBase::IndexType;
+	const box& 	extendedDomainBox_;
 
-	using ExecutionSpace 	= typename interactionBase::ExecutionSpace;
+	/// @brief update interval in terms of iteration numebr 
+	uint32 		updateInterval_= 1;
 
-	using PairContainerType   = unsortedPairs<ExecutionSpace, IdType>;
+	/// @brief  last iteration number which contact search has been performed
+	uint32 		lastUpdated_ 	= 0;
 
-protected:
+	/// @brief performed search? 
+	bool 		performedSearch_ = false;
 
-	const box& 			domain_;
+	/// @brief performed search in boundaries
+	bool 		performedSearchBoundary_ = false;
+
+	/// const ref to particles
+	const particles& 		particles_;
+
+	/// const ref to geometry 
+	const geometry&			geometry_;
+
+	Timer 		bTimer_;
+
+	Timer 		ppTimer_;
 
 	dictionary 	dict_;
 
-	Timer 		sphereSphereTimer_;
+	virtual
+	bool BroadSearch(
+		const timeInfo& ti,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force
+	)=0;
 
-	Timer 		sphereWallTimer_;
-
-	auto& dict()
-	{
-		return dict_;
-	}
+	virtual
+	bool  BoundaryBroadSearch(
+		uint32 bndryIndex,
+		const timeInfo& ti,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force = false
+	)=0;
 
 public:
 
@@ -66,7 +92,7 @@ public:
 
 	contactSearch(
 		const dictionary& dict,
-		const box& domain,
+		const box& extDomain,
 	 	const particles& prtcl,
 	 	const geometry& geom,
 	 	Timers& timers);
@@ -88,40 +114,103 @@ public:
 	 	(dict, domain, prtcl, geom, timers)
 	);
 
-	const auto& domain()const
+	inline
+	bool performedSearch()const
 	{
-		return domain_;
+		return performedSearch_;
 	}
 
-	const auto& dict()const
+	inline 
+	bool performedSearchBoundary()const 
+	{
+		return performedSearchBoundary_;
+	}
+
+	inline
+	bool performSearch(uint32 iter, bool force = false)const
+	{
+		if((iter-lastUpdated_) % updateInterval_ == 0 || iter == 0 || force )
+		{
+			return true;
+		}
+		return false;		
+	}
+
+	bool enterBroadSearch(const timeInfo& ti, bool force = false)const
+	{
+		return performSearch(ti.iter(), force);
+	}
+
+	virtual 
+	bool enterBroadSearchBoundary(const timeInfo& ti, bool force=false)const = 0;
+
+	inline
+	uint32 updateInterval()const
+	{
+		return updateInterval_;
+	}
+
+	inline
+	const dictionary& dict()const
 	{
 		return dict_;
 	}
 
+	inline
+	const box& extendedDomainBox()const
+	{
+		return extendedDomainBox_;
+	}
+
+	inline
+	const particles& Particles()const
+	{
+		return particles_;
+	}
+
+	const pointStructure& pStruct()const;
+
+	inline
+	const geometry& Geometry()const
+	{
+		return geometry_;
+	}
+
+	inline
+	Timer& ppTimer()
+	{
+		return ppTimer_;
+	}
+
+	inline
+	Timer& bTimer()
+	{
+		return bTimer_;
+	}
  
-	virtual 
 	bool broadSearch(
-		PairContainerType& ppPairs,
-		PairContainerType& pwPairs,
-		bool force = false) = 0;
+		const timeInfo& ti,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force = false);
+	
+	bool boundaryBroadSearch(
+		uint32 bndryIndex,
+		const timeInfo& ti,
+		csPairContainerType& ppPairs,
+		csPairContainerType& pwPairs,
+		bool force = false);
 
 	virtual 
-	bool ppEnterBroadSearch()const = 0;
+	real sizeRatio()const = 0;
 
 	virtual 
-	bool pwEnterBroadSearch()const = 0;
-
-	virtual 
-	bool ppPerformedBroadSearch()const = 0;
-
-	virtual 
-	bool pwPerformedBroadSearch()const = 0;
-	
+	real cellExtent()const = 0;
 	
 	static 
 	uniquePtr<contactSearch> create(
 		const dictionary& dict,
-		const box& domain,
+		const box& extDomain,
 	 	const particles& prtcl,
 	 	const geometry& geom,
 	 	Timers& timers);

src/Interaction/contactSearch/contactSearchFunctions.hpp

@@ -96,17 +96,6 @@ void indexToCell(const indexType idx, const iBox<cellIndexType>& box, triple<cel
 	cell+= box.minPoint();
 }
 
-
-
-
-
-
-INLINE_FUNCTION_HD
-bool sphereSphereCheck(const realx3& p1, const realx3 p2, real d1, real d2)
-{
-	return length(p2-p1) < 0.5*(d2+d1);
-}
-
 }
 
 #endif //__broadSearchFunctions_hpp__

src/Interaction/contactSearch/contactSearchGlobals.hpp

@@ -18,28 +18,24 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-#ifndef __interactionTypes_hpp__
-#define __interactionTypes_hpp__
-
-
-
 #include "types.hpp"
+#include "unsortedPairs.hpp"
 
+#ifndef __contactSearchGlobals_hpp__
+#define __contactSearchGlobals_hpp__
 
 namespace pFlow
 {
 
+using csExecutionSpace = DefaultExecutionSpace;
 
-// always use signed types
-using CELL_INDEX_TYPE 	= int32;
+using csIdType = uint32;
 
-using ID_TYPE 			= int32;
+using csPairContainerType   = unsortedPairs<DefaultExecutionSpace, uint32>;
 
-//constexpr int32 minCellIndex = largestNegative<CELL_INDEX_TYPE>();
-
-//constexpr int32 maxCellIndex = largestPositive<CELL_INDEX_TYPE>();
+inline 
+const uint32 BASE_MIRROR_WALL_INDEX = 1000000;
 
 }
 
-
-#endif //__interactionTypes_hpp__
+#endif

src/Interaction/contactSearch/methods/NBS.hpp

@@ -1,209 +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 __NBS_hpp__
-#define __NBS_hpp__
-
-#include "NBSLevel0.hpp"
-
-namespace pFlow
-{
-
-
-template<typename executionSpace>
-class NBS
-{
-public:
-
-	using NBSLevel0Type	= NBSLevel0<executionSpace>;
-
-	using cellIterator 	= typename NBSLevel0Type::cellIterator;
-
-	using IdType 		= typename NBSLevel0Type::IdType;
-	
-	using IndexType 	= typename NBSLevel0Type::IndexType;
-
-	using Cells 		= typename NBSLevel0Type::Cells;
-
-	using CellType 		= typename Cells::CellType;
-
-	using execution_space 	= typename NBSLevel0Type::execution_space;
-
-	using memory_space 		= typename NBSLevel0Type::memory_space;
-
-
-protected:
-
-	real 		sizeRatio_ = 1.0;
-
-	int32 		updateFrequency_= 1;
-
-	int32 		currentIter_ 	= 0;
-
-	bool 		performedSearch_ = false;
-
-	NBSLevel0Type NBSLevel0_;
-
-private:
-
-	bool performSearch()
-	{
-		if(currentIter_ % updateFrequency_ == 0)
-		{
-			currentIter_++;
-			return true;
-
-		}else
-		{
-			currentIter_++;
-			return false;
-		}
-	}
-
-public:
-
-	TypeInfoNV("NBS");
-
-		NBS(
-			const dictionary& dict,
-			const 	box& domain,
-			real  minSize,
-			real 	maxSize,
-			const 	ViewType1D<realx3, memory_space>& position,
-			const 	ViewType1D<real, memory_space>& diam)
-		:
-			sizeRatio_( 
-				max(
-					dict.getValOrSet<real>("sizeRatio", 1.0),
-					1.0
-					)),
-			updateFrequency_(
-				max( 
-					dict.getValOrSet<int32>("updateFrequency", 1),
-					1
-					)),
-			NBSLevel0_(
-				domain, 
-				maxSize*sizeRatio_,
-				sizeRatio_,
-				position,
-				diam)
-		{}
-			
-		INLINE_FUNCTION_HD
-		NBS(const NBS&) = default;
-
-		INLINE_FUNCTION_HD
-		NBS& operator = (const NBS&) = default;
-
-		INLINE_FUNCTION_HD
-		~NBS()=default;
-
-	//// - Methods 
-
-	bool enterBoadSearch()const
-	{
-		return currentIter_%updateFrequency_==0;
-	}
-
-	bool performedSearch()const
-	{
-		return performedSearch_;
-	}
-
-	Vector<cellIterator> getCellIteratorLevels()
-	{
-		return Vector<cellIterator>("cellIterator", 1, NBSLevel0_.getCellIterator());
-	}
-
-	auto getCellIterator(int32 lvl)const
-	{
-		return NBSLevel0_.getCellIterator();
-	}
-
-	int32 numLevels()const
-	{
-		return 1;
-	}
-
-	Vector<Cells> getCellsLevels()const
-	{
-		return Vector<Cells>("Cells", 1, NBSLevel0_.getCells());
-	}
-
-	auto getCells()const
-	{
-		return NBSLevel0_.getCells();
-	}
-
-	bool objectSizeChanged(int32 newSize)
-	{
-		NBSLevel0_.checkAllocateNext(newSize);
-		return true;
-	}
-
-	// - Perform the broad search to find pairs
-	//	 with force = true, perform broad search regardless of 
-	//   updateFrequency_ value
-	//   on all the points in the range of [0,numPoints_)
-	template<typename PairsContainer>
-	bool broadSearch(PairsContainer& pairs, range activeRange, bool force=false)
-	{
-		
-		if(force) currentIter_ = 0;
-		performedSearch_ = false;
-
-		if( !performSearch() ) return true;
-		
-		
-		NBSLevel0_.build(activeRange);
-		
-		NBSLevel0_.findPairs(pairs);
-		
-		performedSearch_ = true;		
-		return true;
-	}
-
-	// - Perform the broad search to find pairs,
-	//   ignore particles with incld(i) = true,
-	//	 with force = true, perform broad search regardless of 
-	//   updateFrequency_ value
-	template<typename PairsContainer, typename IncludeFunction>
-	bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld, bool force = false)
-	{
-		if(force) currentIter_ = 0;
-		performedSearch_ = false;
-
-		if( !performSearch() ) return true;
-
-		NBSLevel0_.build(activeRange, incld);
-	
-		NBSLevel0_.findPairs(pairs);
-		
-		performedSearch_ = true;
-		return true;
-	}
-
-};
-
-}
-
-#endif

src/Interaction/contactSearch/methods/NBSCrossLoop.hpp

@@ -1,82 +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.
-
------------------------------------------------------------------------------*/
-
-
-int32 m = this->head_(i,j,k);
-CellType currentCell(i,j,k);
-int32 n = -1;
-
-while( m > -1 )
-{
-	
-	auto p_m = this->pointPosition_[m];	
-	auto d_m = this->sizeRatio_* this->diameter_[m];
-	
-	int32x3 crossIndex = mapIndexLevels(
-		int32x3(i,j,k),
-		level_, 
-		upperLevel.level());
-
-	
-	for(int32 ii = -1; ii<2; ii++)
-	{
-		for(int32 jj=-1; jj<2; jj++)
-		{
-			int32 kk=-1;
-			while( kk < 2)
-			{
-				int32x3 nghbrCI = crossIndex + int32x3(ii,jj,kk);
-				
-				if( upperLevel.isInRange(nghbrCI) )
-				{
-					n = upperLevel.head_(
-						nghbrCI.x(),
-						nghbrCI.y(),
-						nghbrCI.z());
-
-					while( n >-1)
-					{
-						auto p_n = this->pointPosition_[n];
-						auto d_n = this->sizeRatio_*this->diameter_[n];
-
-						if( sphereSphereCheck(p_m, p_n, d_m, d_n) )
-						{
-							auto ln = n;
-							auto lm = m;
-							if(lm>ln) this->Swap(lm,ln);
-							if( auto res = pairs.insert(lm,ln); res <0)
-							{
-								getFullUpdate++;
-							}
-						}
-
-						n = this->next_[n];
-					}
-
-				}
-
-				kk++;
-			}
-		}
-	}
-
-	m = this->next_[m];
-}
-

src/Interaction/contactSearch/methods/NBSLevel.hpp

@@ -1,127 +0,0 @@
-#ifndef __NBSLevel_hpp__
-#define __NBSLevel_hpp__
-
-
-#include "NBSLevel0.hpp"
-
-
-namespace pFlow
-{
-
-INLINE_FUNCTION_HD
-int32x3 mapIndexLevels(const int32x3& ind, int32 lowerLevel, int32 upperLevel);
-
-template<typename executionSpace>
-class
-NBSLevel
-:
-	public NBSLevel0<executionSpace>
-{
-public:
-
-	using NBSLevel0Type	= NBSLevel0<executionSpace>;
-
-	using cellIterator 	= typename NBSLevel0Type::cellIterator;
-
-	using IdType 		= typename NBSLevel0Type::IdType;
-	
-	using IndexType 	= typename NBSLevel0Type::IndexType;
-
-	using Cells 		= typename NBSLevel0Type::Cells;
-
-	using CellType 		= typename Cells::CellType;
-
-	using execution_space 	= typename NBSLevel0Type::execution_space;
-
-	using memory_space 		= typename NBSLevel0Type::memory_space;
-
-	using mdrPolicyFindPairs = typename NBSLevel0Type::mdrPolicyFindPairs;
-
-	using HeadType 			= typename NBSLevel0Type::HeadType;
-
-	using NextType 			= typename NBSLevel0Type::NextType;
-
-	template<typename exeSpace>
-	friend class NBSLevels;
-
-protected:
-
-	int32 									level_ = 0;
-
-	
-public:
-	
-	TypeInfoNV("NBSLevel0");
-
-	INLINE_FUNCTION_HD
-	NBSLevel(){}
-
-	NBSLevel(
-		int32 	lvl,
-		const 	box& domain,
-		real 	cellSize,
-		real 	sizeRatio,
-		const 	ViewType1D<realx3, memory_space>& position,
-		const 	ViewType1D<real, memory_space>& diam)
-	:
-		NBSLevel0Type(
-			domain,
-			cellSize,
-			sizeRatio,
-			position,
-			diam,
-			lvl==0),
-		level_(lvl)
-	{}
-
-	INLINE_FUNCTION_HD
-	NBSLevel(const NBSLevel&) = default;
-	
-	INLINE_FUNCTION_HD
-	NBSLevel& operator = (const NBSLevel&) = default;
-
-	INLINE_FUNCTION_HD
-	~NBSLevel() = default;
-
-	INLINE_FUNCTION_HD
-	auto level()const
-	{
-		return level_;
-	}
-
-	template<typename PairsContainer>	
-	INLINE_FUNCTION_H
-	int32 findPairsCountCross(PairsContainer& pairs, NBSLevel& upperLevel)
-	{
-
-		mdrPolicyFindPairs
-			mdrPolicy(
-				{0,0,0},
-				{this->nx(),this->ny(),this->nz()} );
-
-		int32 notInsertedPairs;
-
-		Kokkos::parallel_reduce	(
-			"NBSLevel::findPairsCountCross",
-			mdrPolicy,
-			CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){
-				#include "NBSCrossLoop.hpp"
-			}, notInsertedPairs);
-			
-		return notInsertedPairs;	
-	}
-	 
-	
-};
-
-INLINE_FUNCTION_HD
-int32x3 mapIndexLevels( const int32x3& ind, int32 lowerLevel, int32 upperLevel)
-{
-	int32 a = pow(2, static_cast<int32>(upperLevel-lowerLevel));
-	return ind/a;
-}
-
-
-}
-
-#endif

src/Interaction/contactSearch/methods/NBSLevel0.hpp

@@ -1,240 +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 __NBSLevel0_hpp__
-#define __NBSLevel0_hpp__
-
-#include "mapperNBS.hpp"
-
-namespace pFlow
-{
-
-
-template<typename executionSpace>
-class NBSLevel0
-:
-	public mapperNBS<executionSpace>
-{
-public:
-
-	using MapperType 		= mapperNBS<executionSpace>;
-
-	using cellIterator  	= typename MapperType::cellIterator;
-
-	using IdType 			= typename MapperType::IdType;
-	
-	using IndexType 		= typename MapperType::IndexType;
-
-	using Cells 			= typename MapperType::Cells;
-
-	using CellType 			= typename Cells::CellType;
-
-	using execution_space 	= typename MapperType::execution_space;
-
-	using memory_space 		= typename MapperType::memory_space;
-
-	using HeadType 			= typename MapperType::HeadType;
-
-	using NextType 			= typename MapperType::NextType;	
-
-	struct TagFindPairs{};
-
-
-protected:
-
-	real 			sizeRatio_ 		= 1.0;
-
-	// borrowed ownership
-	ViewType1D<real, memory_space>			diameter_;
-
-
-	using mdrPolicyFindPairs = 
-		Kokkos::MDRangePolicy<
-			Kokkos::Rank<3>,
-			Kokkos::Schedule<Kokkos::Dynamic>,
-			execution_space>;
-
-	static INLINE_FUNCTION_HD
-	void Swap(int32& x, int32& y)
-	{
-		int32 tmp = x;
-		x = y;
-		y = tmp;
-	}
-
-public:
-
-	TypeInfoNV("NBSLevel0");
-
-		INLINE_FUNCTION_HD
-		NBSLevel0(){}
-
-		NBSLevel0(
-			const 	box& domain,
-			real 		cellSize,
-			const 	ViewType1D<realx3, memory_space>& position,
-			const 	ViewType1D<real, memory_space>& diam)
-		:
-			MapperType(domain, cellSize, position),
-			diameter_(diam)
-		{}
-
-		NBSLevel0(
-			const 	box& domain,
-			int32 	nx,
-			int32 	ny,
-			int32 	nz,
-			const 	ViewType1D<realx3, memory_space>& position,
-			const 	ViewType1D<real, memory_space>& diam)
-		:
-			MapperType(domain, nx, ny, nz, position),
-			diameter_(diam)
-		{	}
-
-		NBSLevel0(
-			const 	box& domain,
-			real 		cellSize,
-			real 		sizeRatio,
-			const 	ViewType1D<realx3, memory_space>& position,
-			const 	ViewType1D<real, memory_space>& diam,
-			bool 		nextOwner = true)
-		:
-			MapperType(domain, cellSize, position, nextOwner),
-			sizeRatio_(sizeRatio),
-			diameter_(diam)
-		{}
-
-		INLINE_FUNCTION_HD
-		NBSLevel0(const NBSLevel0&) = default;
-
-		INLINE_FUNCTION_HD
-		NBSLevel0& operator = (const NBSLevel0&) = default;
-
-		INLINE_FUNCTION_HD
-		~NBSLevel0()=default;
-
-
-	INLINE_FUNCTION_HD
-	auto sizeRatio()const
-	{
-		return sizeRatio_;
-	}
-
-	INLINE_FUNCTION_HD
-	auto& diameter()
-	{
-		return diameter_;
-	}
-
-	// - Perform the broad search to find pairs
-	//	 with force = true, perform broad search regardless of 
-	//   updateFrequency_ value
-	//   on all the points in the range of [0,numPoints_)
-	template<typename PairsContainer>
-	bool broadSearch(PairsContainer& pairs, range activeRange)
-	{
-		
-		
-		this->build(activeRange);
-		
-		findPairs(pairs);
-			
-		return true;
-	}
-
-	// - Perform the broad search to find pairs,
-	//   ignore particles with incld(i) = true,
-	//	 with force = true, perform broad search regardless of 
-	//   updateFrequency_ value
-	template<typename PairsContainer, typename IncludeFunction>
-	bool broadSearch(PairsContainer& pairs, range activeRange, IncludeFunction incld)
-	{
-		
-		this->build(activeRange, incld);
-	
-		findPairs(pairs);
-		
-		return true;
-	}
-
-	template<typename PairsContainer>	
-	INLINE_FUNCTION_H
-	bool findPairs(PairsContainer& pairs)
-	{
-		
-
-		int32 getFull = 1;
-		
-		
-		// loop until the container size fits the numebr of contact pairs
-		while (getFull > 0)
-		{
-
-			getFull = findPairsCount(pairs);		
-		
-			if(getFull)
-			{
-				// - resize the container
-				//   note that getFull now shows the number of failed insertions.
-				uint32 len = max(getFull,500) ;
-				
-				auto oldCap = pairs.capacity();
-				
-				pairs.increaseCapacityBy(len);
-
-				INFORMATION<< "The contact pair container capacity increased from "<<
-				oldCap << " to "<<pairs.capacity()<<" in NBSLevel0."<<endINFO;
-				
-			}
-
-			Kokkos::fence();
-		}
-		
-		return true;
-	}
-
-	template<typename PairsContainer>	
-	INLINE_FUNCTION_H
-	int32 findPairsCount(PairsContainer& pairs)
-	{
-		mdrPolicyFindPairs
-			mdrPolicy(
-				{0,0,0},
-				{this->nx(),this->ny(),this->nz()} );
-
-		int32 notInsertedPairs;
-			
-		Kokkos::parallel_reduce	(
-			"NBSLevel0::findPairs",
-			mdrPolicy,
-			CLASS_LAMBDA_HD(int32 i, int32 j, int32 k, int32& getFullUpdate){
-				#include "NBSLoop.hpp"
-			}, notInsertedPairs);
-			
-		return notInsertedPairs;
-			
-	}
-
-};
-
-} // pFlow
-
-#endif // __NBSLevel0_hpp__

src/Interaction/contactSearch/methods/NBSLevels.hpp

@@ -1,438 +0,0 @@
-#ifndef __NBSLevels_hpp__
-#define __NBSLevels_hpp__
-
-#include "NBSLevel.hpp"
-#include "NBSLevel0.hpp"
-#include "KokkosTypes.hpp"
-
-namespace pFlow
-{
-
-template<typename executionSpace>
-class NBSLevels
-{
-
-public:
-
-	using NBSLevelType 	= NBSLevel<executionSpace>;
-
-	using cellIterator 	= typename NBSLevelType::cellIterator;
-
-	using IdType 		= typename NBSLevelType::IdType;
-	
-	using IndexType 	= typename NBSLevelType::IndexType;
-
-	using Cells 		= typename NBSLevelType::Cells;
-
-	using CellType 		= typename Cells::CellType;
-
-	using execution_space 	= typename NBSLevelType::execution_space;
-
-	using memory_space 		= typename NBSLevelType::memory_space;
-
-	using realRange 	= kPair<real,real>;
-
-protected:
-	
-	real 		minSize_;
-
-	real 		maxSize_;
-	
-	int32 		numLevels_=1;
-
-	Vector<NBSLevelType> 					nbsLevels_;
-
-	ViewType1D<realRange, memory_space> 	sizeRangeLevels_;
-
-	ViewType1D<realRange, HostSpace> 		sizeRangeLevelsHost_;
-
-	ViewType1D<real, memory_space> 			maxSizeLevels_;
-
-	ViewType1D<real, HostSpace> 			maxSizeLevelsHost_;
-
-	ViewType1D<int8, memory_space> 			particleLevel_;
-
-	range 									activeRange_{0,0};
-
-	using rangePolicyType	=
-		Kokkos::RangePolicy<
-			Kokkos::IndexType<int32>,
-			Kokkos::Schedule<Kokkos::Static>,
-			execution_space>;
-
-	int32 setNumLevels()
-	{
-
-		int32 maxOvermin = static_cast<int32>(maxSize_/minSize_);
-
-		if (maxOvermin <=1)
-			return 1;
-		else if(maxOvermin<=3)
-			return 2;
-		else if(maxOvermin<=7)
-			return 3;
-		else if(maxOvermin<15)
-			return 4;
-		else if(maxOvermin<31)
-			return 5;
-		else if(maxOvermin<63)
-			return 6;
-		else if(maxOvermin <127)
-			return 7;
-		else
-		{
-			fatalErrorInFunction<<
-			"size ratio is not valid for multi-grid NBS "<< maxOvermin<<endl;
-			fatalExit;
-		}
-
-		return -1;
-	}
-
-	bool setDiameterRange(real sizeRatio)
-	{
-		real lvl_maxD  = sizeRatio* maxSize_;
-		real lvl_minD  = lvl_maxD/2.0;		
-
-		for(int32 lvl=numLevels_-1; lvl>=0; lvl--)
-		{
-
-			if(lvl == 0 ) lvl_minD = 0.01*minSize_;
-			
-			sizeRangeLevelsHost_[lvl] = {lvl_minD, lvl_maxD};
-			maxSizeLevelsHost_[lvl] = lvl_maxD;
-			lvl_maxD = lvl_minD;
-			lvl_minD /= 2.0;
-		}
-		
-		copy(sizeRangeLevels_, sizeRangeLevelsHost_);
-		copy(maxSizeLevels_, maxSizeLevelsHost_);
-
-		REPORT(2)<<"Grids with "<< yellowText(numLevels_)<< " levels have been created."<<endREPORT;
-		for(int32 lvl=0; lvl<numLevels_; lvl++)
-		{
-			REPORT(3)<<"Cell gird No "<< yellowText(lvl)<<" with size range ("
-				<<sizeRangeLevelsHost_[lvl].first<<","<<sizeRangeLevelsHost_[lvl].second<<"]."<<endREPORT;
-		}
-		
-		return true;
-	}
-
-	bool initLevels(
-		const 	box& domain,
-		real 	sizeRatio,
-		const 	ViewType1D<realx3, memory_space>& position,
-		const 	ViewType1D<real, memory_space>& diam)
-	{
-		
-
-		for(int32 lvl = 0; lvl<numLevels_; lvl++)
-		{
-			nbsLevels_[lvl] = NBSLevelType(
-				lvl, 
-				domain, 
-				maxSizeLevelsHost_[lvl]*sizeRatio,
-				sizeRatio,
-				position,
-				diam );
-		}
-
-		auto next0 = nbsLevels_[0].next();
-		for(int32 lvl=1; lvl<numLevels_; lvl++)
-		{
-			nbsLevels_[lvl].setNext(next0);
-		}
-
-		return true;
-	}
-
-
-	void manageAllocateNext(range active)
-	{
-		activeRange_ = active;
-
-		if(activeRange_.second > nbsLevels_[0].capacity())
-		{
-			nbsLevels_[0].checkAllocateNext(activeRange_.second);
-
-			auto next0 = nbsLevels_[0].next();
-
-			for(int32 lvl=1; lvl<numLevels_; lvl++)
-			{
-				nbsLevels_[lvl].setNext(next0);
-			}
-		}
-	}
-
-	void nullify( range active)
-	{
-		for(int32 lvl=0; lvl<numLevels_; lvl++)
-		{
-			nbsLevels_[lvl].nullify(active);
-		}
-	}
-
-public:
-
-	NBSLevels(
-		const 	box& domain,
-		real 	minSize,
-		real 	maxSize,
-		real 	sizeRatio,
-		const 	ViewType1D<realx3, memory_space>& position,
-		const 	ViewType1D<real, memory_space>& diam)
-	:
-	minSize_(minSize),
-	maxSize_(maxSize),
-	numLevels_(setNumLevels()),
-	nbsLevels_("nbsLevels", numLevels_, numLevels_, RESERVE()),
-	sizeRangeLevels_("sizeRangeLevels", numLevels_),
-	sizeRangeLevelsHost_("sizeRangeLevelsHost", numLevels_),
-	maxSizeLevels_("maxSizeLevels", numLevels_),
-	maxSizeLevelsHost_("maxSizeLevelsHost", numLevels_)
-	{
-
-		setDiameterRange(sizeRatio);
-
-		initLevels(domain, sizeRatio, position, diam);
-	}
-
-	auto getCellIterator(int32 lvl)const
-	{
-		return nbsLevels_[lvl].getCellIterator();
-	}
-
-	int32 numLevels()const
-	{
-		return numLevels_;
-	}
-
-	Cells getCells(int32 lvl)const
-	{
-		return nbsLevels_[lvl].getCells();
-	}
-
-	template<typename PairsContainer>	
-	INLINE_FUNCTION_H
-	bool findPairs(PairsContainer& pairs)
-	{
-		
-		int32 getFull = 1;
-		
-		// loop until the container size fits the numebr of contact pairs
-		while (getFull > 0)
-		{
-
-			getFull = findPairsCount(pairs);		
-		
-			if(getFull)
-			{
-				// - resize the container
-				//   note that getFull now shows the number of failed insertions.
-				uint32 len = max(getFull,100) ;
-				
-				auto oldCap = pairs.capacity();
-				
-				pairs.increaseCapacityBy(len);
-
-				INFORMATION<< "The contact pair container capacity increased from "<<
-				oldCap << " to "<<pairs.capacity()<<" in NBSLevels."<<endINFO;
-				
-			}
-
-			Kokkos::fence();
-		}
-		
-		return true;
-	}
-
-
-	template<typename PairsContainer>	
-	INLINE_FUNCTION_H
-	int32 findPairsCount(PairsContainer& pairs)
-	{
-
-		int32 notInsertedCount = 0;
-
-		for(int32 lvl=0; lvl<numLevels_; lvl++)
-		{
-			
-			// the same level
-			notInsertedCount+= nbsLevels_[lvl].findPairsCount(pairs);
-			
-			for(int32 crsLvl = lvl+1; crsLvl<numLevels_; crsLvl++)
-			{
-				// cross levels 
-				notInsertedCount+= 
-					nbsLevels_[lvl].findPairsCountCross(pairs, nbsLevels_[crsLvl]);
-				
-			}
-		}
-		
-		return notInsertedCount;
-	}
-
-	INLINE_FUNCTION_H
-	void build(range activeRange)
-	{
-				
-		// nullify next and heads
-		findParticleLevel(activeRange.first, activeRange.second);
-		manageAllocateNext(activeRange);
-		nullify(activeRange);
-		
-		//
-
-		rangePolicyType rPolicy(activeRange.first, activeRange.second);
-		
-		auto nbsLevel0 = nbsLevels_[0];
-		auto pointPosition = nbsLevel0.pointPosition();
-		auto particleLevel = particleLevel_;
-		auto next = nbsLevel0.next();
-		auto head0 = nbsLevel0.head();
-		
-		typename NBSLevelType::HeadType head1, head2, head3, head4, head5, head6;
-
-		if(numLevels_>1) head1 = nbsLevels_[1].head();
-		if(numLevels_>2) head2 = nbsLevels_[2].head();
-		if(numLevels_>3) head3 = nbsLevels_[3].head();
-		if(numLevels_>4) head4 = nbsLevels_[4].head();
-		if(numLevels_>5) head5 = nbsLevels_[5].head();
-		if(numLevels_>6) head6 = nbsLevels_[6].head();
-
-
-		
-		Kokkos::parallel_for(
-			"NBSLevels::build",
-			rPolicy,
-			LAMBDA_HD(int32 i){
-				
-				int8 lvl = particleLevel[i];
-				auto ind = nbsLevel0.pointIndex(pointPosition[i]);
-				ind = mapIndexLevels(ind, 0, lvl);
-				int32 old;
-				if(lvl==0)
-					old =Kokkos::atomic_exchange(&head0(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==1)
-					old =Kokkos::atomic_exchange(&head1(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==2)
-					old =Kokkos::atomic_exchange(&head2(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==3)
-					old =Kokkos::atomic_exchange(&head3(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==4)
-					old =Kokkos::atomic_exchange(&head4(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==5)
-					old =Kokkos::atomic_exchange(&head5(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==6)
-					old =Kokkos::atomic_exchange(&head6(ind.x(), ind.y(), ind.z()),i);
-
-				next(i) = old;
-			});
-
-		Kokkos::fence();
-	}
-
-	template<typename IncludeFunction>
-	INLINE_FUNCTION_H
-	void build(range activeRange, IncludeFunction incld)
-	{
-		// nullify next and heads
-		findParticleLevel(activeRange.first, activeRange.second);
-		manageAllocateNext(activeRange);
-		nullify(activeRange);
-		
-
-		rangePolicyType rPolicy(activeRange.first, activeRange.second);
-
-		auto nbsLevel0 = nbsLevels_[0];
-		auto pointPosition = nbsLevel0.pointPosition();
-		auto particleLevel = particleLevel_;
-		auto next = nbsLevel0.next();
-		auto head0 = nbsLevel0.head();
-
-		typename NBSLevelType::HeadType head1, head2, head3, head4, head5, head6;
-
-		if(numLevels_>1) head1 = nbsLevels_[1].head();
-		if(numLevels_>2) head2 = nbsLevels_[2].head();
-		if(numLevels_>3) head3 = nbsLevels_[3].head();
-		if(numLevels_>4) head4 = nbsLevels_[4].head();
-		if(numLevels_>5) head5 = nbsLevels_[5].head();
-		if(numLevels_>6) head6 = nbsLevels_[6].head();
-
-		Kokkos::parallel_for(
-			"NBSLevels::build",
-			rPolicy,
-			LAMBDA_HD(int32 i){
-				if(!incld(i)) return;
-				
-				int8 lvl = particleLevel[i];
-				auto ind = nbsLevel0.pointIndex(pointPosition[i]);
-				
-				ind = mapIndexLevels(ind, 0, lvl);
-				int32 old;
-				if(lvl==0)
-					old =Kokkos::atomic_exchange(&head0(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==1)
-					old =Kokkos::atomic_exchange(&head1(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==2)
-					old =Kokkos::atomic_exchange(&head2(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==3)
-					old =Kokkos::atomic_exchange(&head3(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==4)
-					old =Kokkos::atomic_exchange(&head4(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==5)
-					old =Kokkos::atomic_exchange(&head5(ind.x(), ind.y(), ind.z()),i);
-				else if(lvl==6)
-					old =Kokkos::atomic_exchange(&head6(ind.x(), ind.y(), ind.z()),i);
-				
-				next(i) = old;
-				
-			});
-
-		Kokkos::fence();
-
-	}
-
-	bool findParticleLevel(int32 first, int32 last)
-	{
-
-		if(last > particleLevel_.size())
-		{
-			reallocNoInit(particleLevel_,last);
-		}
-
-		auto diameter  = nbsLevels_[0].diameter();
-		auto const maxSizes = maxSizeLevels_;
-		auto particleLevel = particleLevel_;
-		auto const sizeRatio = 0.999*nbsLevels_[0].sizeRatio();
-
-		int8 maxLvl = sizeRangeLevels_.size();
-
-		rangePolicyType rPolicy(first, last);
-
-		Kokkos::parallel_for(
-			"NBSLevels::findParticleLevel",
-			rPolicy,
-			LAMBDA_HD(int32 i)
-			{
-				for(int8 lvl = 0; lvl<maxLvl; lvl++)
-				{
-					if( sizeRatio*diameter[i]<= maxSizes[lvl] )
-					{
-						particleLevel[i] = lvl;
-						return;
-					}
-				}
-				particleLevel[i] = static_cast<int8>(-1);
-			});
-		Kokkos::fence();
-
-		return true;
-	}
-
-
-}; //NBSLevels
-
-}
-
-#endif

src/Interaction/contactSearch/methods/NBSLoop.hpp

@@ -1,101 +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.
-
------------------------------------------------------------------------------*/
-
-
-int32 m = this->head_(i,j,k);
-CellType currentCell(i,j,k);
-int32 n = -1;
-
-while( m > -1 )
-{
-	
-	auto p_m = this->pointPosition_[m];	
-	auto d_m = sizeRatio_* diameter_[m];
-	
-	// the same cell
-	n = this->next_(m);
-	
-	while(n >-1)
-	{
-		
-		auto p_n = this->pointPosition_[n];
-		auto d_n = sizeRatio_*diameter_[n];
-		
-		if( sphereSphereCheck(p_m, p_n, d_m, d_n) )
-		{
-			auto ln = n;
-			auto lm = m;
-
-			if(lm>ln) Swap(lm,ln);
-			if( auto res = pairs.insert(lm,ln); res <0)
-			{
-				getFullUpdate++;
-			}
-		}
-		
-		n = this->next_(n);
-	}
-
-	// neighbor cells
-	CellType neighborCell;
-	for(int32 ni=0; ni<13; ni++)
-	{
-		if(ni==0)		neighborCell = currentCell + CellType( 0, 0,-1);
-		else if(ni==1)	neighborCell = currentCell + CellType(-1, 0,-1);
-		else if(ni==2)	neighborCell = currentCell + CellType(-1, 0, 0); 
-		else if(ni==3)	neighborCell = currentCell + CellType(-1, 0, 1); 
-		else if(ni==4)	neighborCell = currentCell + CellType( 0,-1,-1); 
-		else if(ni==5)	neighborCell = currentCell + CellType( 0,-1, 0); 
-		else if(ni==6)	neighborCell = currentCell + CellType( 0,-1, 1); 
-		else if(ni==7)	neighborCell = currentCell + CellType(-1,-1,-1); 
-		else if(ni==8)	neighborCell = currentCell + CellType(-1,-1, 0); 
-		else if(ni==9)	neighborCell = currentCell + CellType(-1,-1, 1);
-		else if(ni==10)	neighborCell = currentCell + CellType( 1,-1,-1);  
-		else if(ni==11)	neighborCell = currentCell + CellType( 1,-1, 0);  
-		else if(ni==12)	neighborCell = currentCell + CellType( 1,-1, 1);  
-		 
-		if( this->isInRange(neighborCell) )
-		{
-			
-			n = this->head_(neighborCell.x(), neighborCell.y(), neighborCell.z());
-			while( n>-1)
-			{
-				
-				auto p_n = this->pointPosition_[n];
-				auto d_n = sizeRatio_*diameter_[n];
-				
-				if(sphereSphereCheck(p_m, p_n, d_m, d_n))
-				{
-					auto ln = n;
-					auto lm = m;
-					if(lm>ln) Swap(lm,ln);
-					if( auto res = pairs.insert(lm,ln); res <0)
-					{
-						getFullUpdate++;
-					}
-				}
-				n = this->next_[n];
-			}
-		}
-
-	}
-	m = this->next_[m];
-}
-